1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Math/Constants.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/DISKinematics.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief H1 energy flow and charged particle spectra
 11  ///
 12  /// @author Peter Richardson
 13  /// Based on the equivalent HZTool analysis
 14  class H1_1994_I372350 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(H1_1994_I372350);
 19
 20
 21    /// @name Analysis methods
 22    /// @{
 23
 24    /// Initialise projections and histograms
 25    void init() {
 26      // Projections
 27      const DISLepton dl;
 28      declare(dl, "Lepton");
 29      declare(DISKinematics(), "Kinematics");
 30      declare(dl.remainingFinalState(), "FS");
 31
 32      // Histos
 33      book(_histEnergyFlowLowX, 1, 1, 1);
 34      book(_histEnergyFlowHighX, 1, 1, 2);
 35
 36      book(_histEECLowX, 2, 1, 1);
 37      book(_histEECHighX, 2, 1, 2);
 38
 39      book(_histSpectraW77, 3, 1, 1);
 40      book(_histSpectraW122, 3, 1, 2);
 41      book(_histSpectraW169, 3, 1, 3);
 42      book(_histSpectraW117, 3, 1, 4);
 43
 44      book(_histPT2, 4, 1, 1);
 45
 46      book(_w77 .first, "TMP/w77_1");
 47      book(_w122.first, "TMP/w122_1");
 48      book(_w169.first, "TMP/w169_1");
 49      book(_w117.first, "TMP/w117_1");
 50      book(_wEnergy.first, "TMP/wEnergy_1");
 51
 52      book(_w77 .second, "TMP/w77_2");
 53      book(_w122.second, "TMP/w122_2");
 54      book(_w169.second, "TMP/w169_2");
 55      book(_w117.second, "TMP/w117_2");
 56      book(_wEnergy.second, "TMP/wEnergy_2");
 57    }
 58
 59
 60    /// Analyse each event
 61    void analyze(const Event& event) {
 62
 63      // Get the DIS kinematics
 64      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 65      if ( dk.failed() ) vetoEvent;
 66      const double x  = dk.x();
 67      const double w2 = dk.W2();
 68      const double w = sqrt(w2);
 69
 70      // Momentum of the scattered lepton
 71      const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 72      if ( dl.failed() )  vetoEvent;
 73      const FourMomentum leptonMom = dl.out();
 74      const double ptel = leptonMom.pT();
 75      const double enel = leptonMom.E();
 76      const double thel = leptonMom.angle(dk.beamHadron().mom())/degree;
 77
 78      // Extract the particles other than the lepton
 79      Particles particles = apply<FinalState>(event, "FS").particles();
 80
 81      // Cut on the forward energy
 82      double efwd = 0.0;
 83      for (const Particle& p : particles) {
 84        const double th = p.angle(dk.beamHadron())/degree;
 85        if (inRange(th, 4.4, 15))  efwd += p.E();
 86      }
 87
 88      // Apply the cuts
 89      // Lepton energy and angle, w2 and forward energy
 90      MSG_DEBUG("enel/GeV = " << enel/GeV << ", thel = " << thel
 91                << ", w2 = " << w2 << ", efwd/GeV = " << efwd/GeV);
 92      bool cut = enel/GeV > 14. && thel > 157. && thel < 172.5 && w2 >= 3000. && efwd/GeV > 0.5;
 93      if (!cut) vetoEvent;
 94
 95      // Weight of the event
 96      if (x < 1e-3) _wEnergy.first->fill();
 97      else          _wEnergy.second->fill();
 98
 99      // Boost to hadronic CM
100      const LorentzTransform hcmboost = dk.boostHCM();
101      // Loop over the particles
102      long ncharged(0);
103      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
104        const Particle& p = particles[ip1];
105
106        const double th = p.angle(dk.beamHadron().momentum()) / degree;
107        // Boost momentum to lab
108        const FourMomentum hcmMom = hcmboost.transform(p.momentum());
109        // Angular cut
110        if (th <= 4.4) continue;
111
112        // Energy flow histogram
113        const double et = fabs(hcmMom.Et());
114        const double eta = hcmMom.eta();
115        (x < 1e-3 ? _histEnergyFlowLowX : _histEnergyFlowHighX)->fill(eta, et);
116        if (PID::charge3(p.pid()) != 0) {
117          /// @todo Use units in w comparisons... what are the units?
118          if (w > 50. && w <= 200.) {
119            double xf= 2 * hcmMom.z() / w;
120            double pt2 = hcmMom.pT2();
121            if (w > 50. && w <= 100.) {
122              _histSpectraW77 ->fill(xf);
123            } else if (w > 100. && w <= 150.) {
124              _histSpectraW122->fill(xf);
125            } else if (w > 150. && w <= 200.) {
126              _histSpectraW169->fill(xf);
127            }
128            _histSpectraW117->fill(xf);
129            /// @todo Is this profile meant to be filled with 2 weight factors?
130            _histPT2->fill(xf, pt2/GeV2);
131            ++ncharged;
132          }
133        }
134
135
136        // Energy-energy correlation
137        if (th <= 8.) continue;
138        double phi1 = p.phi(ZERO_2PI);
139        double eta1 = p.eta();
140        double et1 = fabs(p.momentum().Et());
141        for (size_t ip2 = ip1+1; ip2 < particles.size(); ++ip2) {
142          const Particle& p2 = particles[ip2];
143
144          //double th2 = beamAngle(p2.momentum(), order);
145          double th2 = p2.angle(dk.beamHadron().momentum()) / degree;
146          if (th2 <= 8.) continue;
147          double phi2 = p2.phi(ZERO_2PI);
148
149          /// @todo Use angle function
150          double deltaphi = phi1 - phi2;
151          if (fabs(deltaphi) > PI) deltaphi = fabs(fabs(deltaphi) - TWOPI);
152          double eta2 = p2.eta();
153          double omega = sqrt(sqr(eta1-eta2) + sqr(deltaphi));
154          double et2 = fabs(p2.momentum().Et());
155          double wt = et1*et2 / sqr(ptel);
156          (x < 1e-3 ? _histEECLowX : _histEECHighX)->fill(omega, wt);
157        }
158      }
159
160      // Factors for normalization
161      if (w > 50. && w <= 200.) {
162        if (w <= 100.) {
163          _w77.first ->fill(ncharged);
164          _w77.second->fill();
165        } else if (w <= 150.) {
166          _w122.first ->fill(ncharged);
167          _w122.second->fill();
168        } else {
169          _w169.first ->fill(ncharged);
170          _w169.second->fill();
171        }
172        _w117.first ->fill(ncharged);
173        _w117.second->fill();
174      }
175    }
176
177
178    // Normalize inclusive single particle distributions to the average number of charged particles per event.
179    void finalize() {
180      normalize(_histSpectraW77,  *_w77.first/ *_w77.second);
181      normalize(_histSpectraW122,  *_w122.first/ *_w122.second);
182      normalize(_histSpectraW169,  *_w169.first/ *_w169.second);
183      normalize(_histSpectraW117,  *_w117.first/ *_w117.second);
184
185      scale(_histEnergyFlowLowX , 1./ *_wEnergy.first );
186      scale(_histEnergyFlowHighX, 1./ *_wEnergy.second);
187
188      scale(_histEECLowX , 1./ *_wEnergy.first );
189      scale(_histEECHighX, 1./ *_wEnergy.second);
190    }
191
192    /// @}
193
194
195  private:
196
197    /// Polar angle with right direction of the beam
198    inline double beamAngle(const FourVector& v, bool order) {
199      double thel = v.polarAngle()/degree;
200      if (thel < 0.)  thel += 180.;
201      if (!order)     thel = 180 - thel;
202      return thel;
203    }
204
205    /// @name Histograms
206    /// @{
207    Histo1DPtr _histEnergyFlowLowX, _histEnergyFlowHighX;
208    Histo1DPtr _histEECLowX, _histEECHighX;
209    Histo1DPtr _histSpectraW77, _histSpectraW122, _histSpectraW169, _histSpectraW117;
210    Profile1DPtr _histPT2;
211    /// @}
212
213    /// @name Storage of weights to calculate averages for normalisation
214    /// @{
215    pair<CounterPtr,CounterPtr> _w77, _w122, _w169, _w117, _wEnergy;
216    /// @}
217
218  };
219
220
221
222  RIVET_DECLARE_ALIASED_PLUGIN(H1_1994_I372350, H1_1994_S2919893);
223
224}