1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Tools/ParticleIdUtils.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/DISKinematics.hh"
  6#include "Rivet/Projections/ChargedFinalState.hh"
  7#include "Rivet/Projections/DISKinematics.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Transverse momentum spectra of charged particles in DIS
 13  class H1_1996_I424463 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(H1_1996_I424463);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      // Book projections
 27      declare(DISLepton(), "Lepton");
 28      declare(DISKinematics(), "Kinematics");
 29      declare(ChargedFinalState(), "CFS");
 30      declare(FinalState(), "FS");
 31
 32      // Book histograms
 33      book(_NevAll, "TMP/Nev_all");
 34      _h_dndpt_high_eta_bin.resize(10);
 35      _h_dndpt_low_eta_bin.resize(10);
 36      _hdndeta_pt1_bin.resize(10);
 37      _hdndeta_bin.resize(10);
 38      _hdndptmax_low_eta_bin.resize(8);
 39      int ixx = 0 ;
 40      for (size_t ix = 0; ix < 10; ++ix) {
 41        book(_Nevt_after_cuts[ix], "TMP/Nevt_after_cuts" + to_string(ix));
 42        book(_h_dndpt_high_eta_bin[ix], ix+1, 1, 1);
 43        book(_h_dndpt_low_eta_bin[ix], ix+11, 1, 1);
 44        if (ix != 6 && ix != 9) {
 45          book(_hdndptmax_low_eta_bin[ixx], ixx+21, 1, 1);
 46          ixx=ixx+1;
 47        }
 48        book(_hdndeta_pt1_bin[ix], ix+29, 1, 1);
 49        book(_hdndeta_bin[ix],  ix+39, 1, 1);
 50      }
 51
 52
 53    }
 54
 55
 56    /// Perform the per-event analysis
 57    void analyze(const Event& event) {
 58      // const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 59      const FinalState& fs = apply<FinalState>(event, "FS");
 60      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 61      const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 62
 63      // Get the DIS kinematics
 64      double x  = dk.x();
 65      double y = dk.y();
 66      double Q2 = dk.Q2()/GeV;
 67      double W2 = dk.W2()/GeV;
 68
 69      // Momentum of the scattered lepton
 70      FourMomentum leptonMom = dl.out().momentum();
 71      double enel = leptonMom.E()/GeV;
 72      double thel = 180.-leptonMom.angle(dl.in().momentum())/degree;
 73
 74      _NevAll -> fill() ;
 75
 76      bool cut = y > 0.05 && Q2 > 5. && Q2 < 100.&& enel > 12. && W2 > 4400. && thel > 157. && thel < 173.;
 77      if (!cut) vetoEvent;
 78
 79      int ibin[10] ;
 80      for (int i = 0; i < 10; i++) ibin[i] = 0;
 81      if ( 5. < Q2 && Q2 < 50. && x > 0.0001 && x < 0.0010) ibin[0] = 1;
 82      if ( 5. < Q2 && Q2 < 10. && x > 0.0001 && x < 0.0002) ibin[1] = 1;
 83      if ( 6. < Q2 && Q2 < 10. && x > 0.0002 && x < 0.0005) ibin[2] = 1;
 84      if (10. < Q2 && Q2 < 20. && x > 0.0002 && x < 0.0005) ibin[3] = 1;
 85      if (10. < Q2 && Q2 < 20. && x > 0.0005 && x < 0.0008) ibin[4] = 1;
 86      if (10. < Q2 && Q2 < 20. && x > 0.0008 && x < 0.0015) ibin[5] = 1;
 87      if (10. < Q2 && Q2 < 20. && x > 0.0015 && x < 0.0040) ibin[6] = 1;
 88      if (20. < Q2 && Q2 < 50. && x > 0.0005 && x < 0.0014) ibin[7] = 1;
 89      if (20. < Q2 && Q2 < 50. && x > 0.0014 && x < 0.0030) ibin[8] = 1;
 90      if (20. < Q2 && Q2 < 50. && x > 0.0030 && x < 0.0100) ibin[9] = 1;
 91      for (int i = 0; i < 10; i++) {
 92        if (ibin[i]==1) _Nevt_after_cuts[i]->fill();
 93      }
 94
 95      // Extract the particles other than the lepton
 96      /// @todo Improve to avoid HepMC digging
 97      Particles particles;
 98      particles.reserve(fs.particles().size());
 99      ConstGenParticlePtr dislepGP = dl.out().genParticle();
100      for (const Particle& p : fs.particles()) {
101        ConstGenParticlePtr loopGP = p.genParticle();
102        if (loopGP == dislepGP) continue;
103        particles.push_back(p);
104      }
105
106      // Boost to hadronic CM
107      const LorentzTransform hcmboost = dk.boostHCM();
108
109      // Loop over the particles
110      int mult = 0 ;
111      double ptmax_high[10], ptmax_low[10];
112      for (int i = 0; i < 10; i++) {
113        ptmax_high[i] = 0.; ptmax_low[i] = 0.;
114      }
115      double EtSum = 0;
116      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
117        const Particle& p = particles[ip1];
118
119        double eta = p.momentum().pseudorapidity();
120        // Boost to hcm
121        const FourMomentum hcmMom = hcmboost.transform(p.momentum());
122
123        // Apply safety cuts
124        if (eta > -5  && eta < 10.) {
125          mult = mult + 1;
126          double pThcm =hcmMom.pT() ;
127          double etahcm = hcmMom.pseudorapidity();
128          if (etahcm > 0. && etahcm < 2.0) {
129            EtSum = EtSum + hcmMom.Et();
130          }
131          if (PID::charge(p.pid()) != 0) {
132            if (etahcm > 0.5 && etahcm < 1.5) {
133              for ( int i=0; i< 10; i++) {
134                if (ibin[i]==1) {
135                  _h_dndpt_low_eta_bin[i] ->fill(pThcm);
136                  if (pThcm > ptmax_low[i] ) ptmax_low[i] = pThcm;
137                }
138              }
139
140            }
141            if (etahcm > 1.5 && etahcm < 2.5){
142              for ( int i=0; i< 10; i++) {
143                if (ibin[i]==1) {
144                  _h_dndpt_high_eta_bin[i] ->fill(pThcm);
145                  if (pThcm > ptmax_high[i] ) ptmax_high[i] = pThcm;
146                }
147              }
148            }
149            for ( int i=0; i< 10; i++) {
150              if (ibin[i]==1) _hdndeta_bin[i] ->fill(etahcm);
151              if (ibin[i]==1 && pThcm > 1.) _hdndeta_pt1_bin[i] ->fill(etahcm);
152            }
153          }
154        }  // end of loop over the particles
155      }
156      int ii=0;
157      for ( int i=0; i< 10; i++) {
158        if (i != 6 && i != 9 ) {
159          if ( ibin[i]==1 && EtSum > 6. ) {
160            _hdndptmax_low_eta_bin[ii] ->fill(ptmax_low[i]);
161            // cout << " filling ptmax " << ii << "  " <<  i << endl;
162          }
163          ii=ii+1;
164        }
165      }
166
167    }
168
169
170    /// Normalise histograms etc., after the run
171    void finalize() {
172      MSG_DEBUG("All events: " << _NevAll->val() << " after cuts: "<<  _Nevt_after_cuts[0]->val());
173      MSG_DEBUG("Cut1 events: " << _NevAll->val() << " after cuts: "<<  _Nevt_after_cuts[1]->val());
174      int ii = 0;
175      for ( int i=0; i< 10; i++) {
176        if (_Nevt_after_cuts[i]->val()  != 0) {
177          scale(_h_dndpt_high_eta_bin[i], 1./ *_Nevt_after_cuts[i]);
178          scale(_h_dndpt_low_eta_bin[i], 1./ *_Nevt_after_cuts[i]);
179          scale(_hdndeta_bin[i], 1./ *_Nevt_after_cuts[i]);
180          scale(_hdndeta_pt1_bin[i], 1./ *_Nevt_after_cuts[i]);
181        }
182        if ( i !=6 && i!=9) {
183          // if (_Nevt_after_cuts[i]->val()  != 0) scale(_hdndptmax_low_eta_bin[ii], 1./ *_Nevt_after_cuts[i]); ii=ii+1; };
184          if (_Nevt_after_cuts[i]->val()  != 0) normalize(_hdndptmax_low_eta_bin[ii]); ii=ii+1;
185        }
186      }
187    }
188
189    /// @}
190
191
192    /// @name Histograms
193    /// @{
194    map<string, Histo1DPtr> _h;
195    map<string, Profile1DPtr> _p;
196    map<string, CounterPtr> _c;
197    array<CounterPtr,10> _Nevt_after_cuts;
198    vector<Histo1DPtr> _h_dndpt_low_eta_bin, _h_dndpt_high_eta_bin;
199    vector<Histo1DPtr> _hdndeta_bin, _hdndeta_pt1_bin, _hdndptmax_low_eta_bin;
200    CounterPtr _NevAll ;
201    /// @}
202
203  };
204
205
206  RIVET_DECLARE_ALIASED_PLUGIN(H1_1996_I424463, H1_1997_I424463);
207
208}