rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

H1_1995_I394793

A Study of the Fragmentation of Quarks in ep Collisions at HERA
Experiment: H1 (HERA)
Inspire ID: 394793
Status: VALIDATED HEPDATA
Authors:
  • Narmin Rahimova
  • Hannes Jung
References:
  • Nucl.Phys.B 445 (1995) 3-21
  • DOI:10.1016/0550-3213(95)91599-H
  • arXiv: hep-ex/9505003
Beams: e- p+, p+ e-
Beam energies: (26.7, 820.0); (820.0, 26.7) GeV
    No run details listed

Deep inelastic scattering (DIS) events, selected from 1993 data taken by the H1 experiment at HERA, are studied in the Breit frame of reference. It is shown that certain aspects of the quarks emerging from within the proton in $ep$ interactions are essentially the same as those of quarks pair-created from the vacuum in $e^+e^-$ annihilation.

Source code: H1_1995_I394793.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/PromptFinalState.hh"
  4#include "Rivet/Projections/DISKinematics.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/DISLepton.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief A Study of the Fragmentation of Quarks in ep Collisions at HERA (H1)
 12  class H1_1995_I394793 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(H1_1995_I394793);
 17
 18
 19    /// @name Analysis methods
 20    ///@{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      // Initialise and register projections
 26      declare(DISLepton(), "Lepton");
 27      declare(DISKinematics(), "Kinematics");
 28
 29      // The basic final-state projection:
 30      // all final-state particles within
 31      // the given eta acceptance
 32      const FinalState fs(Cuts::abseta < 4.9);
 33      declare(fs, "FS");
 34      const ChargedFinalState cfs;
 35      declare(cfs, "CFS");
 36      
 37      // book a counter
 38      book(_Nevt_after_cuts, "TMP/Nevt_after_cuts");
 39      book(_Nevt_afterfwd_cuts, "TMP/Nevt_afterfwd_cuts");
 40      book(_Nevt_afterh_cuts, "TMP/Nevt_afterh_cuts");
 41      book(_Nevt_afterhfwd_cuts, "TMP/Nevt_afterhfwd_cuts");
 42
 43
 44      // Book histograms
 45      // specify custom binning
 46      // take binning from reference data using HEPData ID (digits in "d01-x01-y01" etc.)
 47      book(_h["costh_lowQ"], 1, 1, 1);
 48      book(_h["costh_highQ"], 1, 1, 2);     
 49      book(_h["costh_lowQ_noEfwd"], 2, 1, 1);
 50      book(_h["costh_highQ_noEfwd"], 2, 1, 2);
 51      book(_h["xp_posCharge_lowQ"], 3, 1, 1);
 52      book(_h["xp_negCharge_lowQ"], 3, 1, 2);
 53      book(_h["xp_posCharge_highQ"], 3, 1, 3);
 54      book(_h["xp_negCharge_highQ"], 3, 1, 4);
 55      book(_h["ksi_lowQ"], 4, 1, 1);
 56      book(_h["ksi_highQ"], 4, 1, 2);
 57      book(_s["Mult_vrs_Q2"], 5, 1, 1);
 58      book(_s["Mult_vrs_Q2_noEfwd"], 6, 1, 1);
 59      
 60      book(_h["Mult_vrs_Q2_nchrg"],"TMP/Mult_vrs_Q2_nchrg", refData(5,1,1));
 61      book(_h["Mult_vrs_Q2_noEfwd_nchrg"],"TMP/Mult_vrs_Q2_noEfwd_nchrg", refData(6,1,1));
 62      book(_h["Mult_vrs_Q2_count"],"TMP/Mult_vrs_Q2_count", refData(5,1,1));
 63      book(_h["Mult_vrs_Q2_noEfwd_count"],"TMP/Mult_vrs_Q2_noEfwd_count", refData(6,1,1));
 64
 65    }
 66 
 67
 68    /// Perform the per-event analysis
 69    void analyze(const Event& event) {
 70
 71     const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 72    
 73   
 74    //DIS kinematics
 75      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
 76      if ( dk.failed() ) vetoEvent;
 77      double y   = dk.y();
 78      double w2  = dk.W2();
 79      double Q2  = dk.Q2();
 80      
 81      
 82      
 83      bool cut ;
 84      
 85      cut = Q2 > 12 && y < 0.6 && w2 > 3000 ;
 86      
 87      
 88      if ( !cut ) vetoEvent ; 
 89      
 90      const DISLepton& dl = apply<DISLepton>(event,"Lepton");
 91      if ( dl.failed() ) vetoEvent;
 92      /*
 93      cout << "  scattered lepton angle " << 180.- dl.out().momentum().angle(dl.in().momentum())/degree << endl;
 94      cout << " in  lepton " << dl.in().momentum() << endl;
 95      cout << " out lepton " << dl.out().momentum() << endl;
 96      */
 97      
 98      const FinalState& fs = apply<FinalState>(event, "FS");
 99      Particles particles; particles.reserve(fs.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
108      double efwd = 0.;
109
110      for (const Particle& p : particles) {
111        const double th = 180. - p.momentum().angle(dl.in().momentum())/degree;
112
113        if (inRange(th, 4.4, 15.0)) { 
114           efwd += p.E();
115           //cout << " angle " << th << " pid " << p.pid() << " Efwd = " << efwd << endl;
116           } 
117      }
118
119      bool evcut[2];
120      evcut[0] =   efwd > 0.5;
121      
122       // fill the counter
123      _Nevt_after_cuts -> fill();
124      if (Q2 > 100 )  _Nevt_afterh_cuts -> fill();
125       
126      if ( evcut[0] && (Q2 < 80 ) ) _Nevt_afterfwd_cuts -> fill();
127      if ( evcut[0] && (Q2 > 100 ) ) _Nevt_afterhfwd_cuts -> fill();
128
129      double n_charg = 0; 
130
131      // Boost to Breit
132      const LorentzTransform breitboost = dk.boostBreit();
133       
134      for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
135        const Particle& p = particles[ip1];
136        const FourMomentum BreMom = breitboost.transform(p.momentum());
137       // cout << BreMom.pz() << endl;
138        double x = cos(BreMom.theta());
139        
140        if (Q2 < 80 ) _h["costh_lowQ"] ->fill(x);
141        if (Q2 > 100 ) _h["costh_highQ"]  ->fill(x);
142        
143        if (Q2 < 80  && evcut[0]) _h["costh_lowQ_noEfwd"] ->fill(x);
144        if (Q2 > 100 && evcut[0]) _h["costh_highQ_noEfwd"]  ->fill(x);
145
146
147        if ( BreMom.pz() > 0. ) continue;
148        double pcal= sqrt(BreMom.px2() + BreMom.py2()+ BreMom.pz2()) ;
149        double xp = 2*pcal/(sqrt(Q2));
150        double xi = log(1/xp);
151        
152        double charge = p.charge() ;
153       // cout << " charge " << charge << endl;
154       
155        if (charge > 0 ) {
156          if (Q2 < 80 ) _h["xp_posCharge_lowQ"] -> fill(xp);
157          if (Q2 > 100 ) _h["xp_posCharge_highQ"] -> fill(xp);
158        } else {
159         if (Q2 < 80 ) _h["xp_negCharge_lowQ"] -> fill(xp);
160         if (Q2 > 100 ) _h["xp_negCharge_highQ"] -> fill(xp);
161        }
162    
163        if (Q2 < 80 ) _h["ksi_lowQ"] -> fill(xi);
164        if (Q2 > 100 ) _h["ksi_highQ"] -> fill (xi);
165 
166        n_charg = n_charg + 1; 
167
168     }
169     _h["Mult_vrs_Q2_nchrg"] -> fill(Q2,n_charg) ;
170     _h["Mult_vrs_Q2_count"] -> fill(Q2) ;
171     if ( evcut[0]) _h["Mult_vrs_Q2_noEfwd_nchrg"] -> fill(Q2,n_charg) ;
172     if ( evcut[0]) _h["Mult_vrs_Q2_noEfwd_count"] -> fill(Q2) ;
173}
174
175    /// Normalise histograms etc., after the run
176    void finalize() {
177     
178   
179      normalize(_h["xp_posCharge_lowQ"]);
180      if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["xp_posCharge_highQ"], 1.0/ *_Nevt_afterh_cuts);
181
182      normalize(_h["xp_negCharge_lowQ"]);
183      if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["xp_negCharge_highQ"], 1.0/ *_Nevt_afterh_cuts);
184
185      scale(_h["costh_lowQ"], 1.0/ *_Nevt_after_cuts);
186      if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["costh_highQ"], 1.0/ *_Nevt_afterh_cuts);
187      //cout << " after fwd cuts " << dbl(*_Nevt_afterfwd_cuts) << endl;
188      
189      if(dbl(*_Nevt_afterfwd_cuts)>0) scale(_h["costh_lowQ_noEfwd"], 1.0/ *_Nevt_afterfwd_cuts);
190      if(dbl(*_Nevt_afterhfwd_cuts)>0) scale(_h["costh_highQ_noEfwd"], 1.0/ *_Nevt_afterhfwd_cuts);
191      
192      scale(_h["ksi_lowQ"], 1.0/ *_Nevt_after_cuts);
193      if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["ksi_highQ"], 1.0/ *_Nevt_afterh_cuts);
194      //cout << " Nevt " << dbl(*_Nevt_after_cuts) << endl;
195
196      divide(_h["Mult_vrs_Q2_nchrg"], _h["Mult_vrs_Q2_count"], _s["Mult_vrs_Q2"]);
197      divide(_h["Mult_vrs_Q2_noEfwd_nchrg"], _h["Mult_vrs_Q2_noEfwd_count"], _s["Mult_vrs_Q2_noEfwd"]);
198    }
199
200    ///@}
201
202
203    /// @name Histograms
204    ///@{
205    map<string, Histo1DPtr> _h;
206    map<string, Profile1DPtr> _p;
207    map<string, CounterPtr> _c;
208    map<string, Estimate1DPtr> _s;
209    CounterPtr _Nevt_after_cuts;
210    CounterPtr _Nevt_afterfwd_cuts;
211    CounterPtr _Nevt_afterh_cuts;
212    CounterPtr _Nevt_afterhfwd_cuts;
213    ///@}
214
215
216  };
217
218
219  RIVET_DECLARE_PLUGIN(H1_1995_I394793);
220
221}