Rivet analyses referenceH1_1999_I481112Measurement of D* meson cross-sections at HERA and determination of the gluon density in the proton using NLO QCDExperiment: H1 (HERA) Inspire ID: 481112 Status: VALIDATED Authors:
Beam energies: (27.5, 820.0); (820.0, 27.5) GeV
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 Q2>2 GeV2 and in photoproduction at energies around Wγp∼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, xgg(xg), has been extracted in the momentum fraction range 7.5×10−4<xg<4×10−2 at average scales mu2=25 to 50 GeV2. The gluon momentum fraction xg 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 book(_hpt, {2.5, 3.5, 5.0, 10.0}, {"d11-x01-y01", "d11-x01-y02", "d11-x01-y03"});
52 book(_h["1211"], 12,1,1);
53 book(_h["1212"], 12,1,2);
54 book(_h["1311"], 13,1,1);
55 }
56
57
58 /// Perform the per-event analysis
59 void analyze(const Event& event) {
60
61
62 const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
63
64
65 bool isDIS = false;
66 bool ETAG44 = false;
67 bool ETAG33 = false;
68
69 const double y = dk.y();
70 const double Q2 = dk.Q2();
71
72 if (Q2 > 2 && Q2 <100 && y > 0.05 && y < 0.7) isDIS = true;
73 if (Q2 < 0.009 && y > 0.02 && y <0.32 ) ETAG44 = true;
74 if (Q2 < 0.01 && y > 0.29 && y <0.62 ) ETAG33 = true;
75 if (isDIS == false && ETAG44 == false && ETAG33 == false) vetoEvent;
76
77
78 //Creating array of D*
79 Cut cuts = isDIS? (Cuts::pT > 1.5*GeV && Cuts::abseta < 1.5) : (Cuts::pT > 2*GeV && Cuts::absrap < 1.5);
80 Particles unstables = apply<ParticleFinder>(event, "DStars").particles(cuts);
81 const Particles dstars = select(unstables, [](const Particle& p){
82 return p.abspid() == PID::DSPLUS;
83 });
84
85 if(dstars.empty() ) vetoEvent;
86 MSG_DEBUG("D*" << dstars.size());
87
88 const Particle& dstar = dstars.front();
89 // boosting the system
90 const LorentzTransform hcmboost = dk.boostHCM();
91 const FourMomentum hcmMom = hcmboost.transform(dstar.momentum());
92
93 //discriminate between dis and photoprod, and between ETA33 and ETA 44
94
95 //kinematics quantities
96 const double E = dstar.E();
97 const double p_z = dstar.pz();
98 //std::cout<<"y: "<<y<<endl;
99 if (y<0.02) vetoEvent;
100 const double m2 = 2.25; // charm mass^2
101 const double E_e = dk.beamLepton().E();
102 const double z = (E - p_z)/(2*y*E_e);
103 if (z > 1) {
104 MSG_DEBUG("Momentum fraction greater than unity! This should not happen. Vetoing event.");
105 vetoEvent;
106 }
107
108 const double M2 = (1.44*hcmMom.pT2() + m2)/(z*(1-z));
109 const double x_g = (M2 + Q2)/(y*dk.s());
110
111 const double y_capp = dstar.rapidity();
112 const double W = sqrt(dk.W2());
113
114 //perform the cuts
115 if (isDIS == true){
116 _h["211"]->fill(dstar.pT());
117 _h["411"]->fill(dstar.eta());
118 _h["511"]->fill(Q2);
119 _h["611"]->fill(log10(x_g));
120 //boosting to the hcm frame
121 _h["311"]->fill(hcmMom.pT());
122 }
123 if (ETAG33 == true && abs(y_capp) < 1.5 && dstar.pT() > 2.5*GeV ) {
124
125 _h["rap194"]->fill(y_capp);
126 _h["pt194"] ->fill(dstar.pT());
127 _hpt->fill(dstar.pT(), y_capp);
128
129 if ( W > 173 && W<273) {
130 _h["1211"]->fill(log10(x_g));
131 }
132 if ( W > 130 && W<230) {
133 _h["1212"]->fill(log10(x_g));
134 }
135 }
136 if (ETAG44 == true && abs(y_capp) < 1.5 && dstar.pT()> 2) {
137
138 _h["pt88"]->fill(dstar.pT());
139 _h["rap88"]->fill(y_capp);
140 _h["1311"]->fill(log10(x_g));
141
142 }
143 }
144
145 /// Normalise histograms etc., after the run
146 void finalize() {
147 // conversion factors from ep to gamma p xsections (as given in publication)
148 const double F_etag33 = 0.0128;
149 const double F_etag44 = 0.0838;
150
151 double norm = crossSection()/nanobarn/sumW() ;
152 scale( _h["211"], norm);
153
154 scale(_h["311"], norm);
155 scale(_h["411"], norm);
156 scale(_h["511"], norm);
157 scale(_h["611"], norm);
158 double norm_mub = crossSection()/microbarn/sumW();
159
160 scale(_h["rap194"],norm_mub/F_etag33 );
161 scale(_h["pt194"], norm_mub/F_etag33 );
162 scale(_h["rap88"], norm_mub/F_etag44 );
163 scale(_h["pt88"], norm_mub/F_etag44 );
164
165 scale(_hpt, norm/F_etag33);
166 scale(_h["1211"], norm_mub/F_etag33);
167 scale(_h["1212"], norm_mub/F_etag33);
168 scale(_h["1311"], norm_mub/F_etag44);
169 }
170
171 ///@}
172
173
174 /// @name Histograms
175 ///@{
176 map<string, Histo1DPtr> _h;
177 Histo1DGroupPtr _hpt;
178 ///@}
179
180
181 };
182
183
184 RIVET_DECLARE_PLUGIN(H1_1999_I481112);
185
186}
|
