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
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 (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      // Initialise and register projections
 26      //declare(FinalState(Cuts::abseta < 5 && Cuts::pT > 100*MeV), "FS");
 27
 28      // Book histograms
 29      declare(DISLepton(), "Lepton");
 30      declare(DISKinematics(), "Kinematics");
 31      declare(ChargedFinalState(), "CFS");
 32      declare(FinalState(), "FS");
 33
 34 
 35      // Book histograms
 36      // specify custom binning
 37      //book(_h["XXXX"], "myh1", 20, 0.0, 100.0);
 38      //book(_h["YYYY"], "myh2", logspace(20, 1e-2, 1e3));
 39      //book(_h["ZZZZ"], "myh3", {0.0, 1.0, 2.0, 4.0, 8.0, 16.0});
 40      // take binning from reference data using HEPData ID (digits in "d01-x01-y01" etc.)
 41      
 42      book(_NevAll, "TMP/Nev_all");
 43      _h_dndpt_high_eta_bin.resize(10);
 44      _h_dndpt_low_eta_bin.resize(10);
 45      _hdndeta_pt1_bin.resize(10);
 46      _hdndeta_bin.resize(10);
 47      _hdndptmax_low_eta_bin.resize(8);
 48      int ixx = 0 ;
 49      for (size_t ix = 0; ix < 10; ++ix) {
 50        book(_Nevt_after_cuts[ix], "TMP/Nevt_after_cuts" + to_string(ix));
 51        book(_h_dndpt_high_eta_bin[ix], ix+1, 1, 1);
 52        book(_h_dndpt_low_eta_bin[ix], ix+11, 1, 1);
 53        if (ix != 6 && ix != 9) {
 54          book(_hdndptmax_low_eta_bin[ixx], ixx+21, 1, 1);
 55          ixx=ixx+1;
 56        }
 57        book(_hdndeta_pt1_bin[ix], ix+29, 1, 1);
 58        book(_hdndeta_bin[ix],  ix+39, 1, 1);
 59      }
 60
 61
 62    }
 63
 64
 65    /// Perform the per-event analysis
 66    void analyze(const Event& event) {
 67      // const ChargedFinalState& cfs = applyProjection<ChargedFinalState>(event, "CFS");      
 68      const FinalState& fs = apply<FinalState>(event, "FS");
 69      const DISKinematics& dk = applyProjection<DISKinematics>(event, "Kinematics");
 70      const DISLepton& dl = applyProjection<DISLepton>(event,"Lepton");
 71
 72      // Get the DIS kinematics
 73      double x  = dk.x();
 74      double y = dk.y();
 75      double Q2 = dk.Q2()/GeV;
 76      double W2 = dk.W2()/GeV;
 77
 78      // Momentum of the scattered lepton
 79      FourMomentum leptonMom = dl.out().momentum();
 80      double enel = leptonMom.E();
 81      double thel = 180.-leptonMom.angle(dl.in().momentum())/degree;
 82
 83      _NevAll -> fill() ;
 84      
 85       // cout <<"enel/GeV = "<<enel/GeV<<", thel = "<<thel<<", y = "<<y<<", x = "<<x<<std::endl;
 86       bool cut = y > 0.05 && Q2 > 5. && Q2 < 100.&& enel > 12. && W2 > 4400. && thel > 157. && thel < 173.;
 87       if (!cut) vetoEvent;
 88      
 89       int ibin[10] ;
 90       for ( int i=0; i< 10; i++) {
 91       ibin[i] = 0; } 
 92      
 93       if(5.<Q2&&Q2<50.&&x>0.0001&&x<0.001)    ibin[0]=1; 
 94       if(5.<Q2&&Q2<10.&&x>0.0001&&x<0.0002)   ibin[1]=1;
 95       if(6.<Q2&&Q2<10.&&x>0.0002&&x<0.0005)   ibin[2]=1;
 96       if(10.<Q2&&Q2<20.&&x>0.0002&&x<0.0005)  ibin[3]=1;
 97       if(10.<Q2&&Q2<20.&&x>0.0005&&x<0.0008)  ibin[4]=1;
 98       if(10.<Q2&&Q2<20.&&x>0.0008&&x<0.0015)  ibin[5]=1;
 99       if(10.<Q2&&Q2<20.&&x>0.0015&&x<0.004)   ibin[6]=1;
100       if(20.<Q2&&Q2<50.&&x>0.0005&&x<0.0014)  ibin[7]=1;
101       if(20.<Q2&&Q2<50.&&x>0.0014&&x<0.003)   ibin[8]=1;
102       if(20.<Q2&&Q2<50.&&x>0.003&&x<0.01)     ibin[9]=1;
103
104       for ( int i=0; i< 10; i++) { if(ibin[i]==1) _Nevt_after_cuts[i] ->fill(); } 
105
106      // Extract the particles other than the lepton
107      Particles particles;
108      particles.reserve(fs.particles().size());
109      ConstGenParticlePtr dislepGP = dl.out().genParticle();
110      for (const Particle& p : fs.particles()) {
111        ConstGenParticlePtr loopGP = p.genParticle();
112        if (loopGP == dislepGP) continue;
113        particles.push_back(p);
114      }
115
116      // Boost to hadronic CM
117      const LorentzTransform hcmboost = dk.boostHCM();
118      
119      int mult = 0 ;
120      // Loop over the particles
121      // long ncharged(0);
122      double ptmax_high[10], ptmax_low[10] ;
123      for ( int i=0; i< 10; i++) {
124       ptmax_high[i] = 0.; ptmax_low[i] = 0.; } 
125      double EtSum = 0;
126      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
127         const Particle& p = particles[ip1];
128
129         double eta = p.momentum().pseudorapidity();
130         // Boost to hcm
131         const FourMomentum hcmMom = hcmboost.transform(p.momentum());
132   
133         // apply safety cuts
134         if (eta > -5  && eta < 10.){
135            mult = mult + 1;              
136	      double pThcm =hcmMom.pT() ; 
137            double etahcm = hcmMom.pseudorapidity();
138            // cout << " charge " << PID::charge(p.pid()) << endl;
139            if (etahcm > 0. && etahcm < 2.0){ EtSum = EtSum + hcmMom.Et();}
140            if (PID::charge(p.pid()) != 0) {
141              if (etahcm > 0.5 && etahcm < 1.5){ 
142               for ( int i=0; i< 10; i++) { 
143                  if(ibin[i]==1) {
144                    // cout << " fill histo low "<< i << " pt = " << pThcm << endl; 
145                       _h_dndpt_low_eta_bin[i] ->fill(pThcm);
146                       if(pThcm > ptmax_low[i] ) ptmax_low[i] = pThcm;
147                    } 
148               } 
149
150              }   
151              if (etahcm > 1.5 && etahcm < 2.5){
152               for ( int i=0; i< 10; i++) { 
153                  if(ibin[i]==1) 
154                    _h_dndpt_high_eta_bin[i] ->fill(pThcm); 
155                    if(pThcm > ptmax_high[i] ) ptmax_high[i] = pThcm;
156               } 
157              }                    
158              for ( int i=0; i< 10; i++) { 
159              if(ibin[i]==1) _hdndeta_bin[i] ->fill(etahcm); 
160              if(ibin[i]==1 && pThcm > 1.) _hdndeta_pt1_bin[i] ->fill(etahcm); 
161              }
162           }  
163         }  // end of loop over the particles
164      }
165      int ii=0;
166      for ( int i=0; i< 10; i++) { 
167         if (i != 6 && i != 9 ) { 
168             if( ibin[i]==1 && EtSum > 6. ) {
169              _hdndptmax_low_eta_bin[ii] ->fill(ptmax_low[i]); 
170               // cout << " filling ptmax " << ii << "  " <<  i << endl; 
171             }
172            ii=ii+1;
173         }
174      }
175
176    }
177    /// Normalise histograms etc., after the run
178    void finalize() {
179    cout << " All events: " << _NevAll->val() << " after cuts: "<<  _Nevt_after_cuts[0]->val() << endl;
180    cout << " cut1 events: " << _NevAll->val() << " after cuts: "<<  _Nevt_after_cuts[1]->val() << endl;
181     
182
183     int ii =0;
184     for ( int i=0; i< 10; i++) { 
185        if (_Nevt_after_cuts[i]->val()  != 0) {
186            scale(_h_dndpt_high_eta_bin[i], 1./ *_Nevt_after_cuts[i]);  
187            scale(_h_dndpt_low_eta_bin[i], 1./ *_Nevt_after_cuts[i]); 
188            scale(_hdndeta_bin[i], 1./ *_Nevt_after_cuts[i]); 
189            scale(_hdndeta_pt1_bin[i], 1./ *_Nevt_after_cuts[i]); }
190        if ( i !=6 && i!=9) { 
191           // if (_Nevt_after_cuts[i]->val()  != 0) scale(_hdndptmax_low_eta_bin[ii], 1./ *_Nevt_after_cuts[i]); ii=ii+1; }; 
192            if (_Nevt_after_cuts[i]->val()  != 0) normalize(_hdndptmax_low_eta_bin[ii]); ii=ii+1; }; 
193     } 
194    }
195
196    ///@}
197
198
199    /// @name Histograms
200    ///@{
201    map<string, Histo1DPtr> _h;
202    map<string, Profile1DPtr> _p;
203    map<string, CounterPtr> _c;
204    array<CounterPtr,10> _Nevt_after_cuts;
205    vector<Histo1DPtr> _h_dndpt_low_eta_bin, _h_dndpt_high_eta_bin;
206    vector<Histo1DPtr> _hdndeta_bin, _hdndeta_pt1_bin, _hdndptmax_low_eta_bin;
207    CounterPtr _NevAll ;
208   ///@}
209
210
211  };
212
213
214  DECLARE_RIVET_PLUGIN(H1_1996_I424463);
215
216}