rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

H1_2013_I1217865

Charged particle production in deep-inelastic ep scattering at H1
Experiment: H1 (HERA)
Inspire ID: 1217865
Status: VALIDATED
Authors:
  • Anastasia Grebenyuk
  • Hannes Jung
References:
  • Eur.Phys.J. C73 (2013) 2406,
  • arXiv: 1302.1321
Beams: e+ p+, p+ e+, e- p+, p+ e-
Beam energies: ANY
Run details:
  • Inclusive DIS

Charged particle production in deep-inelastic ep scattering is measured with the H1 detector at HERA. The kinematic range of the analysis covers low photon virtualities, 5 < Q2 < 100 GeV2 and small values of Bjorken-x, 10**-2 < x < 10**-2. The analysis is performed in the hadronic centre-of-mass system. The charged particle densities are measured as a function of pseudorapidity eta* and transverse momentum pT* in the range 0 <eta* <5 and 0 < pT* < 10 GeV in bins of x and Q2.

Source code: H1_2013_I1217865.cc
  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 Charged particle production in deep-inelastic ep scattering at H1
 13  class H1_2013_I1217865 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17
 18    RIVET_DEFAULT_ANALYSIS_CTOR(H1_2013_I1217865);
 19
 20
 21    /// @name Analysis methods
 22    /// @{
 23
 24    /// Book histograms and initialise projections before the run
 25    void init() {
 26
 27      // Initialise and register projections
 28      //declare(FinalState(Cuts::abseta < 5 && Cuts::pT > 100*MeV), "FS");
 29
 30      // Book histograms
 31      declare(DISLepton(), "Lepton");
 32      declare(DISKinematics(), "Kinematics");
 33      declare(ChargedFinalState(), "CFS");
 34      declare(FinalState(), "FS");
 35      _h_dn_dpT_cen.resize(9);
 36      _h_dn_dpT_curr.resize(9);
 37      _h_dn_deta_soft.resize(9);
 38      _h_dn_deta_hard.resize(9);
 39
 40      book(_h_dn_dpT_cen[0],19,1,1);
 41      book(_h_dn_dpT_curr[0] ,20,1,1);
 42      book(_h_dn_deta_soft[0],1,1,1);
 43      book(_h_dn_deta_hard[0],2,1,1);
 44      for (size_t ix = 0; ix < 9; ++ix) {
 45        book(_Nevt_after_cuts[ix], "TMP/Nevt_after_cuts" + to_string(ix));
 46        if (ix > 0 ) {
 47          book(_h_dn_dpT_cen[ix], ix+20, 1, 1);
 48	    book(_h_dn_dpT_curr[ix],ix+28,1,1);
 49	    book(_h_dn_deta_soft[ix],ix+2,1,1);
 50	    book(_h_dn_deta_hard[ix],ix+10,1,1);
 51        }
 52      }
 53
 54
 55    }
 56
 57
 58    /// Perform the per-event analysis
 59    void analyze(const Event& event) {
 60      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 61      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 62      const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 63
 64      // Get the DIS kinematics
 65      double x  = dk.x();
 66      double y = dk.y();
 67      double Q2 = dk.Q2()/GeV;
 68
 69      // Momentum of the scattered lepton
 70      FourMomentum leptonMom = dl.out().momentum();
 71      double enel = leptonMom.E();
 72      double thel = 180.-leptonMom.angle(dl.in().momentum())/degree;
 73
 74
 75      getLog()<<Log::DEBUG<<"enel/GeV = "<<enel/GeV<<", thel = "<<thel<<", y = "<<y<<", x = "<<x<<std::endl;
 76       bool cut = y > 0.05 && y < 0.6 && Q2 > 5. && Q2 < 100.;
 77       if (!cut) vetoEvent;
 78
 79
 80       int ibin[10] ;
 81       for ( int i=0; i< 9; i++) {
 82       ibin[i] = 0; }
 83
 84       ibin[0] = 1 ;
 85       if(5.<Q2&&Q2<10.&&x>0.0001&&x<0.00024)   ibin[1]=1;
 86       if(5.<Q2&&Q2<10.&&x>0.00024&&x<0.0005)   ibin[2]=1;
 87       if(5.<Q2&&Q2<10.&&x>0.0005&&x<0.002)     ibin[3]=1;
 88       if(10.<Q2&&Q2<20.&&x>0.0002&&x<0.00052)  ibin[4]=1;
 89       if(10.<Q2&&Q2<20.&&x>0.00052&&x<0.0011)  ibin[5]=1;
 90       if(10.<Q2&&Q2<20.&&x>0.0011&&x<0.0037)   ibin[6]=1;
 91       if(20.<Q2&&Q2<100.&&x>0.0004&&x<0.0017)  ibin[7]=1;
 92       if(20.<Q2&&Q2<100.&&x>0.0017&&x<0.01)    ibin[8]=1;
 93
 94       for ( int i=0; i< 9; i++) { if(ibin[i]==1) _Nevt_after_cuts[i] ->fill(); }
 95
 96
 97      // Extract the particles other than the lepton
 98      Particles particles;
 99      particles.reserve(cfs.particles().size());
100      ConstGenParticlePtr dislepGP = dl.out().genParticle();
101      for (const Particle& p : cfs.particles()) {
102        ConstGenParticlePtr loopGP = p.genParticle();
103        if (loopGP == dislepGP) continue;
104        particles.push_back(p);
105      }
106
107      // Boost to hadronic CM
108      const LorentzTransform hcmboost = dk.boostHCM();
109
110      int mult = 0 ;
111      // Loop over the particles
112      // long ncharged(0);
113      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
114      const Particle& p = particles[ip1];
115
116      double eta = p.momentum().pseudorapidity();
117      double pT = p.momentum().pT()/GeV;
118
119      // Boost to hcm
120      const FourMomentum hcmMom = hcmboost.transform(p.momentum());
121
122      if (pT > 0.15 && eta > -2. && eta < 2.5){
123
124        mult = mult + 1;
125
126	  double pThcm =hcmMom.pT() ;
127        double etahcm = hcmMom.pseudorapidity();
128
129
130        if (etahcm > 0. && etahcm < 1.5){
131
132            _h_dn_dpT_cen[0]->fill(pThcm );
133	      for ( int i=1; i< 9; i++) {
134                  if(ibin[i]==1) { _h_dn_dpT_cen[i]->fill(pThcm );}}
135       }
136
137       if (etahcm > 1.5 && etahcm < 5.){
138           _h_dn_dpT_curr[0]->fill(pThcm );
139	      for ( int i=1; i< 9; i++) {
140                  if(ibin[i]==1) { _h_dn_dpT_curr[i]->fill(pThcm );}}
141       }
142
143       if(pThcm < 1.){
144           _h_dn_deta_soft[0]->fill(etahcm );
145	      for ( int i=1; i< 9; i++) {
146                  if(ibin[i]==1) {_h_dn_deta_soft[i]->fill(etahcm );}}
147      }
148
149      if(pThcm > 1. && pThcm < 10.){
150            _h_dn_deta_hard[0]->fill(etahcm );
151	      for ( int i=1; i< 9; i++) {
152                  if(ibin[i]==1) {_h_dn_deta_hard[i]->fill(etahcm );}}
153
154       }
155    }  // if (etahcm > 0. && etahcm < 1.5){
156  }  // end of loop over the particles
157
158
159
160    }
161
162
163    /// Normalise histograms etc., after the run
164    void finalize() {
165    //  normalize(_h_dn_dpT_cen);
166
167     if (_Nevt_after_cuts[0]->val()  != 0)  scale(_h_dn_dpT_cen[0],  1./ *_Nevt_after_cuts[0]);
168     if (_Nevt_after_cuts[0]->val()  != 0)  scale(_h_dn_dpT_curr[0],  1./ *_Nevt_after_cuts[0]);
169     if (_Nevt_after_cuts[0]->val()  != 0)  scale(_h_dn_deta_soft[0],  1./ *_Nevt_after_cuts[0]);
170     if (_Nevt_after_cuts[0]->val()  != 0)  scale(_h_dn_deta_hard[0],   1./ *_Nevt_after_cuts[0]);
171
172
173     for ( int i=1; i< 9; i++) {
174        if (_Nevt_after_cuts[i]->val()  != 0) {
175           scale(_h_dn_dpT_cen[i],  1./ *_Nevt_after_cuts[i]);
176           scale(_h_dn_dpT_curr[i], 1./ *_Nevt_after_cuts[i]);
177           scale(_h_dn_deta_soft[i],1./ *_Nevt_after_cuts[i]);
178           scale(_h_dn_deta_hard[i],1./ *_Nevt_after_cuts[i]);
179        }
180     }
181
182
183
184    }
185  private:
186
187    /**
188     *  Polar angle with right direction of the beam
189     */
190    inline double beamAngle(const FourVector& v, const bool & order) {
191      double thel = v.polarAngle()/degree;
192      if(thel<0.) thel+=180.;
193      if(!order) thel = 180.-thel;
194      return thel;
195    }
196
197    /// @}
198
199
200    /// @name Histograms
201    /// @{
202      Histo1DPtr _h_dn_dpT_2r;
203      Histo1DPtr _h_dn_dpT_2l;
204
205      vector<Histo1DPtr> _h_dn_dpT_cen;
206      vector<Histo1DPtr> _h_dn_dpT_curr;
207      vector<Histo1DPtr> _h_dn_deta_soft;
208      vector<Histo1DPtr> _h_dn_deta_hard;
209      array<CounterPtr,9> _Nevt_after_cuts;
210
211
212    /// @}
213
214
215  };
216
217
218  RIVET_DECLARE_PLUGIN(H1_2013_I1217865);
219
220
221}