rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

H1_1996_I424463

Transverse momentum spectra of charged particles in DIS (H1 1996)
Experiment: H1 (HERA)
Inspire ID: 424463
Status: VALIDATED
Authors:
  • Suraj Kumar Singh
  • Hannes Jung
  • Andrii Verbytskyi
References: Beams: e+ p+, p+ e+
Beam energies: (27.5, 820.0); (820.0, 27.5) GeV
    No run details listed

Transverse momentum spectra of charged particles produced in deep inelastic scattering are measured as a function of the kinematic variables x and Q2 using the H1 detector at the epcollider HERA. The data are compared to different patton emission models, either with or without ordering of the emissions in transverse momentum. The data provide evidence for a relatively large amount of patton radiation between the current and the remnant systems.

Source code: H1_1996_I424463.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 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      int ixx = 0 ;
 35      for (size_t ix = 0; ix < 10; ++ix) {
 36        book(_Nevt_after_cuts[ix], "TMP/Nevt_after_cuts" + to_string(ix));
 37        for(unsigned int ih=0;ih<2;++ih) {
 38          book(_h_dndpt_eta_bin[ih][ix], ih*10 + ix+1 , 1, 1);
 39          book(_hdndeta_bin    [ih][ix], ih*10 + ix+29, 1, 1);
 40        }
 41        if (ix != 6 && ix != 9) {
 42          book(_hdndptmax_low_eta_bin[ixx], ixx+21, 1, 1);
 43          ixx=ixx+1;
 44        }
 45      }
 46    }
 47
 48
 49    /// Perform the per-event analysis
 50    void analyze(const Event& event) {
 51      if (_edgespT.empty()) {
 52        _edgespT    = _h_dndpt_eta_bin   [0][0]->xEdges();
 53        _edgesEta   = _hdndeta_bin       [0][0]->xEdges();
 54        _edgespTMax = _hdndptmax_low_eta_bin[0]->xEdges();
 55      }
 56      const FinalState& fs = apply<FinalState>(event, "FS");
 57      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 58      const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 59
 60      // Get the DIS kinematics
 61      double x  = dk.x();
 62      double y = dk.y();
 63      double Q2 = dk.Q2()/GeV;
 64      double W2 = dk.W2()/GeV;
 65
 66      // Momentum of the scattered lepton
 67      FourMomentum leptonMom = dl.out().momentum();
 68      double enel = leptonMom.E()/GeV;
 69      double thel = 180.-leptonMom.angle(dl.in().momentum())/degree;
 70
 71      _NevAll -> fill() ;
 72
 73      bool cut = y > 0.05 && Q2 > 5. && Q2 < 100.&& enel > 12. && W2 > 4400. && thel > 157. && thel < 173.;
 74      if (!cut) vetoEvent;
 75
 76      int ibin[10] ;
 77      for (int i = 0; i < 10; i++) ibin[i] = 0;
 78      if ( 5. < Q2 && Q2 < 50. && x > 0.0001 && x < 0.0010) ibin[0] = 1;
 79      if ( 5. < Q2 && Q2 < 10. && x > 0.0001 && x < 0.0002) ibin[1] = 1;
 80      if ( 6. < Q2 && Q2 < 10. && x > 0.0002 && x < 0.0005) ibin[2] = 1;
 81      if (10. < Q2 && Q2 < 20. && x > 0.0002 && x < 0.0005) ibin[3] = 1;
 82      if (10. < Q2 && Q2 < 20. && x > 0.0005 && x < 0.0008) ibin[4] = 1;
 83      if (10. < Q2 && Q2 < 20. && x > 0.0008 && x < 0.0015) ibin[5] = 1;
 84      if (10. < Q2 && Q2 < 20. && x > 0.0015 && x < 0.0040) ibin[6] = 1;
 85      if (20. < Q2 && Q2 < 50. && x > 0.0005 && x < 0.0014) ibin[7] = 1;
 86      if (20. < Q2 && Q2 < 50. && x > 0.0014 && x < 0.0030) ibin[8] = 1;
 87      if (20. < Q2 && Q2 < 50. && x > 0.0030 && x < 0.0100) ibin[9] = 1;
 88      for (int i = 0; i < 10; i++) {
 89        if (ibin[i]==1) _Nevt_after_cuts[i]->fill();
 90      }
 91
 92      // Extract the particles other than the lepton
 93      /// @todo Improve to avoid HepMC digging
 94      Particles particles;
 95      particles.reserve(fs.particles().size());
 96      ConstGenParticlePtr dislepGP = dl.out().genParticle();
 97      for (const Particle& p : fs.particles()) {
 98        ConstGenParticlePtr loopGP = p.genParticle();
 99        if (loopGP == dislepGP) continue;
100        particles.push_back(p);
101      }
102
103      // Boost to hadronic CM
104      const LorentzTransform hcmboost = dk.boostHCM();
105
106      // Loop over the particles
107      int mult = 0 ;
108      double ptmax_high[10], ptmax_low[10];
109      for (int i = 0; i < 10; i++) {
110        ptmax_high[i] = 0.; ptmax_low[i] = 0.;
111      }
112      double EtSum = 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        // Boost to hcm
118        const FourMomentum hcmMom = hcmboost.transform(p.momentum());
119
120        // Apply safety cuts
121        if (eta > -5  && eta < 10.) {
122          mult = mult + 1;
123          double pThcm =hcmMom.pT() ;
124          double etahcm = hcmMom.pseudorapidity();
125          if (etahcm > 0. && etahcm < 2.0) {
126            EtSum = EtSum + hcmMom.Et();
127          }
128          if (PID::charge(p.pid()) != 0) {
129            if (etahcm > 0.5 && etahcm < 1.5) {
130              for ( int i=0; i< 10; i++) {
131                if (ibin[i]==1) {
132                  fillpT(1,i,pThcm);
133                  if (pThcm > ptmax_low[i] ) ptmax_low[i] = pThcm;
134                }
135              }
136
137            }
138            if (etahcm > 1.5 && etahcm < 2.5){
139              for ( int i=0; i< 10; i++) {
140                if (ibin[i]==1) {
141                  fillpT(0,i,pThcm);
142                  if (pThcm > ptmax_high[i] ) ptmax_high[i] = pThcm;
143                }
144              }
145            }
146            for ( int i=0; i< 10; i++) {
147              if (ibin[i]==1)               fillEta(1,i,etahcm);
148              if (ibin[i]==1 && pThcm > 1.) fillEta(0,i,etahcm);
149            }
150          }
151        }  // end of loop over the particles
152      }
153      int ii=0;
154      for ( int i=0; i< 10; i++) {
155        if (i != 6 && i != 9 ) {
156          if ( ibin[i]==1 && EtSum > 6. ) {
157            fillpTMax(ii,ptmax_low[i]);
158          }
159          ii=ii+1;
160        }
161      }
162    }
163
164
165    /// Normalise histograms etc., after the run
166    void finalize() {
167      MSG_DEBUG("All events: " << _NevAll->val() << " after cuts: "<<  _Nevt_after_cuts[0]->val());
168      MSG_DEBUG("Cut1 events: " << _NevAll->val() << " after cuts: "<<  _Nevt_after_cuts[1]->val());
169      int ii = 0;
170      for ( int i=0; i< 10; i++) {
171        if (_Nevt_after_cuts[i]->val()  != 0) {
172          for(unsigned int ih=0;ih<2;++ih) {
173            scale(_h_dndpt_eta_bin[ih][i], 1./ *_Nevt_after_cuts[i]);
174            size_t ioff=0;
175            if( (ih==0 && (i==6 || i==9)) ||
176                (ih==1 && (i<=4 || i==7)))    ioff=1;
177            else if (ih==1 && (i==5 || i==8)) ioff=2;
178            for(auto & b : _h_dndpt_eta_bin[ih][i]->bins()) {
179              const size_t idx = b.index()+ioff;
180              b.scaleW(1./_axispT.width(idx));
181            }
182            scale(_hdndeta_bin    [ih][i], 1./ *_Nevt_after_cuts[i]);
183            ioff=0;
184            if( (ih==0 && (i==5 || i==6 || i==8 || i==9)) ||
185                (ih==1 && (i==0 || i==4 || i==5 || i==8))) ioff=1;
186            else if (ih==1 && (i==6 || i==9))       ioff=2;
187            for(auto & b : _hdndeta_bin[ih][i]->bins()) {
188              const size_t idx = b.index()+ioff;
189              b.scaleW(1./_axisEta.width(idx));
190            }
191          }
192        }
193        if ( i !=6 && i!=9) {
194          if (_Nevt_after_cuts[i]->val()  != 0) {
195            normalize(_hdndptmax_low_eta_bin[ii]);
196            for(auto & b : _hdndptmax_low_eta_bin[ii]->bins()) {
197              const size_t idx = b.index();
198              b.scaleW(1./_axispTMax.width(idx));
199            }
200          }
201          ii=ii+1;
202        }
203      }
204      
205    }
206    /// @}
207
208    void fillEta(const unsigned int ix, const unsigned int iy, const double value) {
209      string edge = "OTHER";
210      const size_t idx = _axisEta.index(value);
211      if (idx && idx <= _edgesEta.size()) {
212        if ( ( (ix==0 && iy>=1&& iy<=3) || (ix==0 && (iy==1 || iy==3)) ) && idx == _edgesEta.size())
213          ;
214        else if (idx==1 && ( (ix==0 && (iy==5 || iy==6 || iy==8 || iy==9)) ||
215                             (ix==1 && (iy==0 || iy==4 || iy==5 || iy==8))))
216          ;
217        else if (idx<=2 && (ix==1 && (iy==6 || iy==9)) )
218          ;
219        else {
220          edge = _edgesEta[idx-1];
221        }
222      }
223      _hdndeta_bin[ix][iy]->fill(edge);
224    }
225    
226    void fillpT(const unsigned int ix, const unsigned int iy, const double value) {
227      string edge = "OTHER";
228      const size_t idx = _axispT.index(value);
229      if (idx && idx <= _edgespT.size()) {
230        if(((ix==0 && iy==8) || (ix==1 && (iy==2 || iy==4 || iy==6))) && idx == _edgespT.size())
231          ;
232        else if (idx==1 && ( (ix==0 && (iy==6 || iy==9)) ||
233                             (ix==1 && (iy<=4 || iy==7))))
234          ;
235        else if (idx<=2 && (ix==1 && (iy==5 || iy==8)))
236          ;
237        else
238          edge = _edgespT[idx-1];
239      }
240      _h_dndpt_eta_bin[ix][iy]->fill(edge);
241    }
242
243    void fillpTMax(const unsigned int ix, const double value) {
244      string edge = "OTHER";
245      const size_t idx = _axispTMax.index(value);
246      if (idx && idx <= _edgespTMax.size()) {
247        if(idx==_edgespTMax.size() && (ix==1 || ix==2 || ix==4))
248          ;
249        else if(idx>=_edgespTMax.size()-1 && ix==5)
250          ;
251        else
252          edge = _edgespTMax[idx-1];
253      }
254      _hdndptmax_low_eta_bin[ix]->fill(edge);
255    }
256
257    /// @name Histograms
258    /// @{
259    array<CounterPtr,10> _Nevt_after_cuts;
260    BinnedHistoPtr<string> _h_dndpt_eta_bin[2][10], _hdndeta_bin[2][10], _hdndptmax_low_eta_bin[8];
261    CounterPtr _NevAll ;
262
263    vector<string> _edgespT,_edgesEta,_edgespTMax;
264    YODA::Axis<double> _axispT    = YODA::Axis<double>{ 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.825, 2.125, 2.525, 3.125, 4.0, 5.0};
265    YODA::Axis<double> _axisEta   = YODA::Axis<double>{ 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.00};
266    YODA::Axis<double> _axispTMax = YODA::Axis<double>{0.005, 0.255, 0.505, 0.755, 1.005, 1.255, 1.505, 1.755, 2.065, 2.5, 3.125, 4.0, 5.0};
267    /// @}
268
269  };
270
271
272  RIVET_DECLARE_ALIASED_PLUGIN(H1_1996_I424463, H1_1997_I424463);
273
274}