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 = apply<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      Cut cuts = isDIS? (Cuts::pT > 1.5*GeV && Cuts::abseta < 1.5) : (Cuts::pT > 2*GeV && Cuts::absrap < 1.5); 
 82      Particles unstables = apply<ParticleFinder>(event, "DStars").particles(cuts);
 83      const Particles dstars = filter_select(unstables, [](const Particle& p){
 84        return p.abspid() == PID::DSPLUS;
 85      });
 86
 87      if(dstars.empty() ) vetoEvent;
 88      MSG_DEBUG("D*" << dstars.size());
 89
 90      const Particle& dstar = dstars.front();
 91      // boosting the system
 92      const LorentzTransform hcmboost = dk.boostHCM();
 93      const FourMomentum hcmMom = hcmboost.transform(dstar.momentum());
 94
 95      //discriminate between dis and photoprod, and between ETA33 and ETA 44
 96
 97      //kinematics quantities
 98      const double E = dstar.E();
 99      const double p_z = dstar.pz();
100      //std::cout<<"y: "<<y<<endl;
101      if (y<0.02)  vetoEvent;
102      const double m2 = 2.25; // charm mass^2
103      const double E_e = dk.beamLepton().E();
104      const double z = (E - p_z)/(2*y*E_e);
105      if (z > 1) {
106        MSG_DEBUG("Momentum fraction greater than unity! This should not happen. Vetoing event.");
107        vetoEvent;
108      }
109
110      const double M2 = (1.44*hcmMom.pT2() + m2)/(z*(1-z));
111      const double x_g = (M2 + Q2)/(y*dk.s());
112
113      const double y_capp = dstar.rapidity();
114      const double W = sqrt(dk.W2());
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        if ( W > 130 &&  W<230) {
135          _h["1212"]->fill(log10(x_g));
136        }
137      }
138      if (ETAG44 == true && abs(y_capp) < 1.5 && dstar.pT()> 2) { 
139
140        _h["pt88"]->fill(dstar.pT());
141        _h["rap88"]->fill(y_capp);
142        _h["1311"]->fill(log10(x_g));
143
144      }
145    }
146
147    /// Normalise histograms etc., after the run
148    void finalize() {
149      // conversion factors from ep to gamma p xsections (as given in publication)
150      const double F_etag33 = 0.0128;
151      const double F_etag44 = 0.0838;
152
153      double norm = crossSection()/nanobarn/sumW() ;
154      scale( _h["211"], norm);
155      
156      scale(_h["311"], norm);  
157      scale(_h["411"], norm);
158      scale(_h["511"], norm);
159      scale(_h["611"], norm);
160      double norm_mub = crossSection()/microbarn/sumW();
161
162      scale(_h["rap194"],norm_mub/F_etag33 );
163      scale(_h["pt194"], norm_mub/F_etag33 );
164      scale(_h["rap88"], norm_mub/F_etag44 );
165      scale(_h["pt88"],  norm_mub/F_etag44 );
166         
167      _hpt.scale( norm/F_etag33, this);
168      scale(_h["1211"], norm_mub/F_etag33);
169      scale(_h["1212"], norm_mub/F_etag33);
170      scale(_h["1311"], norm_mub/F_etag44);
171    }
172
173    ///@}
174
175
176    /// @name Histograms
177    ///@{
178    map<string, Histo1DPtr> _h;
179    BinnedHistogram  _hpt; 
180    ///@}
181
182
183  };
184
185
186  RIVET_DECLARE_PLUGIN(H1_1999_I481112);
187
188}