rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

H1_1999_I481112

Measurement of D* meson cross-sections at HERA and determination of the gluon density in the proton using NLO QCD
Experiment: H1 (HERA)
Inspire ID: 481112
Status: VALIDATED
Authors:
  • Luca Marsili
  • Hannes Jung
References: Beams: e+ p+, p+ e+
Beam energies: (27.5, 820.0); (820.0, 27.5) GeV
    No run details listed

With the H1 detector at the $ep$ collider HERA, $D^*$ meson production cross sections have been measured in deep inelastic scattering with four-momentum transfers $Q^2>2$ GeV$^2$ and in photoproduction at energies around $W_{\gamma p} \sim 88$ GeV and 194 GeV. Next-to-Leading Order QCD calculations are found to describe the differential cross sections within theoretical and experimental uncertainties. Using these calculations, the NLO gluon momentum distribution in the proton, $x_g g(x_g)$, has been extracted in the momentum fraction range $7.5 \times 10^{-4}< x_g <4 \times 10^{-2}$ at average scales $mu^2 =25$ to 50 GeV$^2$. The gluon momentum fraction $x_g$ has been obtained from the measured kinematics of the scattered electron and the $D^*$ meson in the final state.

Source code: H1_1999_I481112.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/DISKinematics.hh"
  6#include "Rivet/Projections/UnstableParticles.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Measurement of D* meson cross-sections at HERA
 12  class H1_1999_I481112 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(H1_1999_I481112);
 17
 18
 19    /// @name Analysis methods
 20    ///@{
 21    
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25     
 26      // Initialise and register projections
 27
 28      // The basic final-state projection:
 29      // all final-state particles within
 30      // the given eta acceptance
 31      const FinalState fs(Cuts::abseta < 1.5);
 32      declare(fs, "fs");
 33      // The final-state particles declared above are clustered using FastJet with
 34    
 35      //Initialize quantities needed for cuts
 36      declare(DISKinematics(), "Kinematics");
 37      declare(UnstableParticles(), "DStars");
 38      
 39
 40      Histo1DPtr dummy;
 41     
 42      book(_h["211"],2,1,1);
 43      book(_h["311"],3,1,1);
 44      book(_h["411"],4,1,1);
 45      book(_h["511"],5,1,1);
 46      book(_h["611"],6,1,1);
 47      book(_h["rap194"], 7,1,1);
 48      book(_h["pt194"], 8,1,1);
 49      book(_h["rap88"], 9,1,1);
 50      book(_h["pt88"], 10,1,1);
 51      _hpt.add(2.5, 3.5, book(dummy, 11,1,1));
 52      _hpt.add(3.5, 5.0, book(dummy, 11,1,2));
 53      _hpt.add(5.0, 10.0, book(dummy, 11,1,3));
 54      book(_h["1211"],12,1,1);
 55      book(_h["1212"],12,1,2);
 56      book(_h["1311"],13,1,1);
 57    }
 58    
 59    
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62     
 63      
 64      const DISKinematics& dk = applyProjection<DISKinematics>(event, "Kinematics");
 65      
 66      
 67      bool isDIS = false;
 68      bool ETAG44 = false;
 69      bool ETAG33 = false;
 70
 71      const double y = dk.y();
 72      const double Q2 = dk.Q2();
 73
 74      if(Q2 > 2 && Q2 <100 &&  y > 0.05 && y < 0.7) isDIS = true;
 75      if(Q2 < 0.009 && y > 0.02 && y <0.32  ) ETAG44 = true;
 76      if(Q2 < 0.01 &&  y > 0.29 && y <0.62 ) ETAG33 = true;
 77      if(isDIS == false && ETAG44 == false && ETAG33 == false) vetoEvent;
 78      
 79      
 80      //Creating array of D*
 81      Particles unstables ;
 82      if (isDIS ) { unstables = apply<ParticleFinder>(event, "DStars").particles(Cuts::pT > 1.5*GeV && Cuts::absetaIn(0,1.5));}
 83      else { unstables = apply<ParticleFinder>(event, "DStars").particles(Cuts::pT > 2.0*GeV && Cuts::absrapIn(0,1.5));}
 84      const Particles dstars = filter_select(unstables, [](const Particle& p){return p.abspid() == PID::DSPLUS;});
 85      
 86      if(dstars.empty() ) vetoEvent;
 87      MSG_DEBUG("D*" << dstars.size());
 88      
 89      const Particle& dstar = dstars.front();
 90      // boosting the system
 91      const LorentzTransform hcmboost = dk.boostHCM();
 92      const FourMomentum hcmMom = hcmboost.transform(dstar.momentum());
 93
 94      //discriminate between dis and photoprod, and between ETA33 and ETA 44
 95
 96      //kinematics quantities
 97      const double E = dstar.E();
 98      const double p_z = dstar.pz();
 99      //std::cout<<"y: "<<y<<endl;
100      if(y<0.02) vetoEvent;
101      const double m2 = 2.25; // charm mass^2
102      const double E_e = dk.beamLepton().E();
103      const double z = (E - p_z)/(2*y*E_e);
104      
105     
106      const double M2 = (1.44*hcmMom.pT2() + m2)/(z*(1-z));
107      const double x_g = (M2 + Q2)/(y*dk.s());
108      // std::cout<<"s:  "<<dk.s()<<endl;
109      // std::cout<<"M2: "<<M2<<endl;
110      //  std::cout<<"z: "<<z<<endl;
111      
112      const double y_capp = dstar.rapidity();
113      const double W = sqrt(dk.W2());
114      //std::cout<<"x_g: "<<x_g<<endl;
115      
116      //perform the cuts
117      if(isDIS == true){
118         _h["211"]->fill(dstar.pT());
119         _h["411"]->fill(dstar.eta());
120         _h["511"]->fill(Q2);
121         _h["611"]->fill(log10(x_g));
122		//boosting to the hcm frame
123         _h["311"]->fill(hcmMom.pT());
124      }
125      if(ETAG33 == true && abs(y_capp) < 1.5 && dstar.pT() > 2.5*GeV  ) {
126		
127         _h["rap194"]->fill(y_capp);
128         _h["pt194"] ->fill(dstar.pT());
129         _hpt.fill(dstar.pT(), y_capp);	
130	
131	if( W > 173 &&  W<273){
132	   _h["1211"]->fill(log10(x_g));
133	  
134	}
135	if( W > 130 &&  W<230){
136	   _h["1212"]->fill(log10(x_g));
137	  
138	}
139	
140	
141      }
142      if(ETAG44 == true && abs(y_capp) < 1.5 && dstar.pT()> 2) { 
143
144	   _h["pt88"]->fill(dstar.pT());
145	   _h["rap88"]->fill(y_capp);
146	   _h["1311"]->fill(log10(x_g));
147	
148      }
149    }
150
151    /// Normalise histograms etc., after the run
152    void finalize() {
153      // conversion factors from ep to gamma p xsections (as given in publication)
154      const double F_etag33 = 0.0128;
155      const double F_etag44 = 0.0838;
156
157      double norm = crossSection()/nanobarn/sumW() ;
158      scale( _h["211"], norm);
159      
160      scale(_h["311"], norm);  
161      scale(_h["411"], norm);
162      scale(_h["511"], norm);
163      scale(_h["611"], norm);
164      double norm_mub = crossSection()/microbarn/sumW();
165
166      scale(_h["rap194"],norm_mub/F_etag33 );
167      scale(_h["pt194"], norm_mub/F_etag33 );
168      scale(_h["rap88"], norm_mub/F_etag44 );
169      scale(_h["pt88"],  norm_mub/F_etag44 );
170         
171      _hpt.scale( norm/F_etag33, this);
172      scale(_h["1211"], norm_mub/F_etag33);
173      scale(_h["1212"], norm_mub/F_etag33);
174      scale(_h["1311"], norm_mub/F_etag44);
175    }
176
177    ///@}
178
179
180    /// @name Histograms
181    ///@{
182    map<string, Histo1DPtr> _h;
183    map<string, Profile1DPtr> _p;
184    map<string, CounterPtr> _c;
185    BinnedHistogram  _hpt; 
186    ///@}
187
188
189  };
190
191
192  RIVET_DECLARE_PLUGIN(H1_1999_I481112);
193
194}