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 (H1 1996)
 13  class H1_1996_I424463 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    DEFAULT_RIVET_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 histograms
 27      declare(DISLepton(), "Lepton");
 28      declare(DISKinematics(), "Kinematics");
 29      declare(ChargedFinalState(), "CFS");
 30      declare(FinalState(), "FS");
 31
 32
 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 = applyProjection<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      const 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++) {
 81      ibin[i] = 0; }
 82
 83      if(5.<Q2&&Q2<50.&&x>0.0001&&x<0.001)    ibin[0]=1;
 84      if(5.<Q2&&Q2<10.&&x>0.0001&&x<0.0002)   ibin[1]=1;
 85      if(6.<Q2&&Q2<10.&&x>0.0002&&x<0.0005)   ibin[2]=1;
 86      if(10.<Q2&&Q2<20.&&x>0.0002&&x<0.0005)  ibin[3]=1;
 87      if(10.<Q2&&Q2<20.&&x>0.0005&&x<0.0008)  ibin[4]=1;
 88      if(10.<Q2&&Q2<20.&&x>0.0008&&x<0.0015)  ibin[5]=1;
 89      if(10.<Q2&&Q2<20.&&x>0.0015&&x<0.004)   ibin[6]=1;
 90      if(20.<Q2&&Q2<50.&&x>0.0005&&x<0.0014)  ibin[7]=1;
 91      if(20.<Q2&&Q2<50.&&x>0.0014&&x<0.003)   ibin[8]=1;
 92      if(20.<Q2&&Q2<50.&&x>0.003&&x<0.01)     ibin[9]=1;
 93
 94      for ( int i=0; i< 10; i++) { if(ibin[i]==1) _Nevt_after_cuts[i] ->fill(); }
 95
 96      // Extract the particles other than the lepton
 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      int mult = 0 ;
110      // Loop over the particles
111      // long ncharged(0);
112      double ptmax_high[10], ptmax_low[10] ;
113      for ( int i=0; i< 10; i++) {
114       ptmax_high[i] = 0.; ptmax_low[i] = 0.; }
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            // cout << " charge " << PID::charge(p.pid()) << endl;
129            if (etahcm > 0. && etahcm < 2.0) {
130              EtSum = EtSum + hcmMom.Et();
131            }
132            if (PID::charge(p.pid()) != 0) {
133              if (etahcm > 0.5 && etahcm < 1.5) {
134               for ( int i=0; i< 10; i++) {
135                 if (ibin[i]==1) {
136                   _h_dndpt_low_eta_bin[i] ->fill(pThcm);
137                   if (pThcm > ptmax_low[i] ) ptmax_low[i] = pThcm;
138                 }
139               }
140
141              }
142              if (etahcm > 1.5 && etahcm < 2.5){
143               for ( int i=0; i< 10; i++) {
144                 if (ibin[i]==1) {
145                   _h_dndpt_high_eta_bin[i] ->fill(pThcm);
146                   if (pThcm > ptmax_high[i] ) ptmax_high[i] = pThcm;
147                 }
148               }
149              }
150              for ( int i=0; i< 10; ++i) {
151                if(ibin[i]==1) _hdndeta_bin[i] ->fill(etahcm);
152                if(ibin[i]==1 && pThcm > 1.) _hdndeta_pt1_bin[i] ->fill(etahcm);
153              }
154           }
155         }  // end of loop over the particles
156      }
157      int ii=0;
158      for ( int i=0; i< 10; ++i) {
159         if (i != 6 && i != 9 ) {
160             if( ibin[i]==1 && EtSum > 6. ) {
161              _hdndptmax_low_eta_bin[ii] ->fill(ptmax_low[i]);
162             }
163            ii=ii+1;
164         }
165      }
166
167    }
168    /// Normalise histograms etc., after the run
169    void finalize() {
170
171
172      int ii =0;
173      for (int i=0; i < 10; ++i) {
174        if (_Nevt_after_cuts[i]->val() != 0) {
175          scale(_h_dndpt_high_eta_bin[i], 1./ *_Nevt_after_cuts[i]);
176          scale(_h_dndpt_low_eta_bin[i],  1./ *_Nevt_after_cuts[i]);
177          scale(_hdndeta_bin[i],          1./ *_Nevt_after_cuts[i]);
178          scale(_hdndeta_pt1_bin[i],      1./ *_Nevt_after_cuts[i]); 
179        }
180        if ( i !=6 && i!=9) {
181          if (_Nevt_after_cuts[i]->val()  != 0) normalize(_hdndptmax_low_eta_bin[ii]); ii=ii+1; 
182        }
183      }
184    }
185
186    ///@}
187
188
189    /// @name Histograms
190    ///@{
191    map<string, Histo1DPtr> _h;
192    map<string, Profile1DPtr> _p;
193    map<string, CounterPtr> _c;
194    array<CounterPtr,10> _Nevt_after_cuts;
195    vector<Histo1DPtr> _h_dndpt_low_eta_bin, _h_dndpt_high_eta_bin;
196    vector<Histo1DPtr> _hdndeta_bin, _hdndeta_pt1_bin, _hdndptmax_low_eta_bin;
197    CounterPtr _NevAll;
198    ///@}
199
200
201  };
202
203
204  RIVET_DECLARE_ALIASED_PLUGIN(H1_1996_I424463, H1_1997_I424463);
205
206}