Rivet analyses referenceH1_2007_I736052Measurement of inclusive production of D* mesons both with and without dijet production in DIS collisions (H1)Experiment: H1 (HERA) Inspire ID: 736052 Status: VALIDATED Authors:
Beam energies: (27.5, 920.0); (920.0, 27.5) GeV
Inclusive $D^*$ production is measured in deep-inelastic $ep$ scattering at HERA with the H1 detector. In addition, the production of dijets in events with a $D^*$ meson is investigated. The analysis covers values of photon virtuality $2< Q^2 \leq 100$ GeV$^2$ and of inelasticity $0.05 \leq y \leq 0.7$. Differential cross sections are measured as a function of $Q^2$ and $x$ and of various $D^*$ meson and jet observables. Source code: H1_2007_I736052.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/DISFinalState.hh"
5#include "Rivet/Projections/FastJets.hh"
6#include "Rivet/Projections/DISKinematics.hh"
7#include "Rivet/Projections/DISLepton.hh"
8#include "Rivet/Projections/UnstableParticles.hh"
9
10namespace Rivet {
11
12
13 /// @brief Measurement of inclusive production of D* mesons both with and without dijet production in DIS collisions (H1)
14 class H1_2007_I736052 : public Analysis {
15 public:
16
17 /// Constructor
18 RIVET_DEFAULT_ANALYSIS_CTOR(H1_2007_I736052);
19
20
21 /// @name Analysis methods
22 ///@{
23
24 /// Book histograms and initialise projections before the run
25 void init() {
26
27 declare(DISLepton(), "Lepton");
28 declare(DISKinematics(), "Kinematics");
29
30 const FinalState fs;
31 declare(fs, "FS");
32 const UnstableParticles ufs;
33 declare(ufs, "UFS");
34
35 double jet_radius = 1.0;
36 const DISFinalState DISfs(DISFinalState::BoostFrame::BREIT);
37 declare(FastJets(DISfs, fastjet::JetAlgorithm::kt_algorithm, fastjet::RecombinationScheme::Et_scheme, jet_radius), "jets_fs");
38
39 Histo1DPtr tmp;
40 book(_h["411"], 4, 1, 1);
41 book(_h["511"], 5, 1, 1);
42 book(_h["611"], 6, 1, 1);
43 book(_h["711"], 7, 1, 1);
44 book(_h["811"], 8, 1, 1);
45 book(_h["911"], 9, 1, 1);
46 _h_binned["Q2xbj"].add( 2., 4.22, book(_h["1011"], 10, 1, 1));
47 _h_binned["Q2xbj"].add( 4.22, 10., book(_h["1111"], 11, 1, 1));
48 _h_binned["Q2xbj"].add( 10., 17.8, book(_h["1211"], 12, 1, 1));
49 _h_binned["Q2xbj"].add( 17.8, 31.6, book(_h["1311"], 13, 1, 1));
50 _h_binned["Q2xbj"].add( 31.6, 100., book(_h["1411"], 14, 1, 1));
51 book(_h["1511"], 15, 1, 1);
52 book(_h["1611"], 16, 1, 1);
53 book(_h["1711"], 17, 1, 1);
54 book(_h["1811"], 18, 1, 1);
55 book(_h["1911"], 19, 1, 1);
56 book(_h["2011"], 20, 1, 1);
57 book(_h["2111"], 21, 1, 1);
58 book(_h["2211"], 22, 1, 1);
59
60 _h_binned["Q2phi"].add( 2., 10., book(tmp, 23, 1, 1));
61 _h_binned["Q2phi"].add( 10., 100., book(tmp, 24, 1, 1));
62 book(_h["2511"], 25, 1, 1);
63 book(_h["2611"], 26, 1, 1);
64 book(_h["2711"], 27, 1, 1);
65 book(_h["2811"], 28, 1, 1);
66 _h_binned["Q2xgam"].add( 2., 5., book(_h["2911"], 29, 1, 1));
67 _h_binned["Q2xgam"].add( 5., 10., book(_h["2912"], 29, 1, 2));
68 _h_binned["Q2xgam"].add( 10., 100., book(_h["2913"], 29, 1, 3));
69 book(_h["3011"], 30, 1, 1);
70 _h_binned["Q2xglue"].add( 2., 5., book(_h["3111"], 31, 1, 1));
71 _h_binned["Q2xglue"].add( 5., 10., book(_h["3112"], 31, 1, 2));
72 _h_binned["Q2xglue"].add( 10., 100., book(_h["3113"], 31, 1, 3));
73 }
74
75
76 /// Perform the per-event analysis
77 void analyze(const Event& event) {
78
79 const FinalState& fs = apply<FinalState>(event, "FS");
80 const size_t numParticles = fs.particles().size();
81
82 Jets jets_fs = apply<JetAlg>(event, "jets_fs").jetsByPt(); // Jets with cut on eta
83 double jet_radius = 1.0;
84
85 const UnstableParticles& ufs = apply<UnstableFinalState>(event, "UFS");
86
87 if (numParticles < 2){
88 MSG_DEBUG("Failed leptonic event cut");
89 vetoEvent;
90 }
91
92 Particles Dstar;
93 for(const Particle& p : filter_select(ufs.particles(), Cuts::pT > 1.5*GeV and Cuts::pT < 15*GeV and Cuts::abseta < 1.5 and Cuts::abspid==413)) {
94 Dstar.push_back(p);
95 }
96 if(Dstar.size() == 0){ // Cut on Dstar
97 MSG_DEBUG("Failed Dstar cut");
98 vetoEvent;
99 }
100
101 const DISKinematics& dk = applyProjection<DISKinematics>(event, "Kinematics");
102 const DISLepton& dl = applyProjection<DISLepton>(event,"Lepton");
103
104
105 double Q2 = dk.Q2();
106 double y = dk.y();
107
108 if(y < 0.05 or y > 0.7 or Q2 < 2 or Q2 > 100){ // Cut on event kinematics
109 MSG_DEBUG("Failed kinematics cut");
110 vetoEvent;
111 }
112
113 // Extract the particles other than the lepton
114 Particles particles;
115 particles.reserve(fs.particles().size());
116 ConstGenParticlePtr dislepGP = dl.out().genParticle();
117 for (const Particle& p : fs.particles()) {
118 ConstGenParticlePtr loopGP = p.genParticle();
119 if (loopGP == dislepGP) continue;
120 particles.push_back(p);
121 }
122
123
124
125 const LorentzTransform hcmboost = dk.boostHCM(); // Hadron cm system
126 const LorentzTransform breitboost = dk.boostBreit(); //Breit system
127 const LorentzTransform labboost = breitboost.inverse(); //Labsystem from Breit
128
129 double xbj = dk.x();
130 double W2 = dk.W2();
131 double pT_cm(0);
132
133 _h["411"] -> fill(Q2);
134 _h["511"] -> fill(xbj);
135 _h["611"] -> fill(std::sqrt(W2));
136
137 _h_binned["Q2xbj"].fill(Q2, xbj);
138
139 for(const Particle& p : Dstar){
140 _h["711"] -> fill(p.pT());
141 _h["811"] -> fill(p.eta());
142 _h["911"] -> fill((p.E() - p.pz())/(2*y*dk.beamLepton().E()));
143
144 const FourMomentum hcmMom = hcmboost.transform(p.momentum());
145 if(pT_cm < hcmMom.pT()) pT_cm = hcmMom.pT();
146
147 if(hcmMom.pT() > 2){ //Cut for 1711 and 1811 histograms
148 _h["1711"] -> fill(p.pT());
149 _h["1811"] -> fill(p.eta());
150 }
151 }
152
153 if(pT_cm > 2){
154 _h["1511"] -> fill(Q2);
155 _h["1611"] -> fill(xbj);
156 }
157
158/*
159 FourMomentum gammaZ;
160 const FourMomentum leptonOUT = dl.out();
161 const FourMomentum leptonIN = dl.in();
162 gammaZ = leptonIN - leptonOUT ;
163 cout << " LAB: gammaZ " << gammaZ << endl;
164 cout << " HCM: gammaZ " << hcmboost.transform(gammaZ) << endl;
165 cout << " LAB: p-beam " << dk.beamHadron().momentum() << endl;
166 cout << " HCM: p-beam " << hcmboost.transform(dk.beamHadron().momentum()) << endl;
167 cout << " Breit: gammaZ " << breitboost.transform(gammaZ) << endl;
168 cout << " Breit: p-beam " << breitboost.transform(dk.beamHadron().momentum()) << endl;
169 for(ConstGenParticlePtr p: HepMCUtils::particles(event.genEvent())) {
170 const PdgId pid = p->pdg_id();
171 if (abs(pid) == 23) {
172 cout<< " HEPMC: photon/Z " << p->momentum() << endl;
173 }
174 }
175*/
176 bool two_jets_with_cut(false);
177
178 Jets jet_cut;
179
180 for(const Jet& j : jets_fs){
181 double etalab = labboost.transform(j.momentum()).eta() ;
182 if(j.momentum().Et() > 3 and etalab > -1 and etalab < 2.5){ // Cut on jet energy in Breit sys.
183 jet_cut.push_back(j); }
184 }
185
186 if(jet_cut.size() >= 2 && jet_cut[0].momentum().Et() > 4 ) two_jets_with_cut = true;
187
188 double delta_phi;
189 FourMomentum jet1;
190 FourMomentum jet2;
191
192 Jets jet_DJ, jet_OJ;
193 bool found_DJ(false), found_OJ(false);
194 //cout << " check D and other jets " << endl;
195 if(two_jets_with_cut){
196 jet1 = jet_cut[0].momentum(); // momentum of jet #1 in Breit sys.
197 jet2 = jet_cut[1].momentum(); // momentum of jet #2 in Breit sys.
198
199 delta_phi = deltaPhi(jet1,jet2)/degree ;
200 _h["1911"] -> fill(Q2);
201 _h["2011"] -> fill(xbj);
202 _h["2111"] -> fill(jet_cut[0].momentum().Et()/GeV);
203 _h["2211"] -> fill(FourMomentum(jet1+jet2).mass()/GeV); // Jets invariant mass in Breit sys.
204 //cout << " dphi " << delta_phi <<" " << deltaPhi(jet_cut[0].momentum(),jet_cut[1].momentum())/degree << endl;
205 _h_binned["Q2phi"].fill(Q2, delta_phi);
206 for (const Jet& jet : jet_cut) {
207 for(const Particle & p : Dstar) {
208 if(deltaR(breitboost.transform(p.momentum()), jet.momentum()) < jet_radius ) {
209 jet_DJ.push_back(jet);
210 found_DJ = true;
211 }
212 else {
213 jet_OJ.push_back(jet);
214 found_OJ = true;
215 }
216 }
217
218 if( found_DJ && found_OJ ) break ;
219 }
220 //cout << " DJ jet size " << jet_DJ.size() << "found Dj " << found_DJ << " OJ jet size " << jet_DJ.size() << " found OJ " << found_OJ << endl;
221 if(jet_DJ.size()>0 && jet_OJ.size()>0 )
222 {
223 double eta_DJ_breit = jet_DJ[0].momentum().eta();
224 double eta_OJ_breit = jet_OJ[0].momentum().eta();
225
226 _h["2511"] -> fill(eta_DJ_breit);
227 _h["2611"] -> fill(eta_OJ_breit);
228 _h["2711"] -> fill(abs(eta_DJ_breit - eta_OJ_breit));
229
230 double x_gamma, x_gluon;
231 double E_star_p_z_star(0); // E() - pz() sum for for x_gamma in \gamma p cm
232 double E_star_p_z_had(0);
233
234 // for jets in had cms: boost first back to lab and then to hcm
235 Jet jet_DJ_hcm, jet_OJ_hcm;
236 jet_DJ_hcm = hcmboost.transform(labboost(jet_DJ[0].momentum()));
237 jet_OJ_hcm = hcmboost.transform(labboost(jet_OJ[0].momentum()));
238
239// Need to change sign: by default hcmboost has gamma* in +z dir, and p in -z dir, but here we need: gamma* -in -z and proton in +z.
240// so - -> +: watch out for E-pz -> E+pz and exp(eta_OJ_hcm) -> exp(-eta_OJ_hcm)
241// this was already noted in H1_2007_I746380.cc
242
243
244
245 double Et_DJ_hcm = jet_DJ_hcm.Et();
246 double eta_DJ_hcm = jet_DJ_hcm.eta();
247 double Et_OJ_hcm = jet_OJ_hcm.Et();
248 double eta_OJ_hcm = jet_OJ_hcm.eta();
249
250 //observed fraction of the photon momentum carried by the parton involved in the hard subprocess
251
252 E_star_p_z_star = jet_DJ_hcm.momentum().E() + jet_DJ_hcm.momentum().pz() + jet_OJ_hcm.momentum().E() + jet_OJ_hcm.momentum().pz();
253
254 for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
255 const Particle& p = particles[ip1];
256 E_star_p_z_had += (hcmboost.transform(p.momentum())).E() + (hcmboost.transform(p.momentum())).pz();
257 }
258 x_gamma = E_star_p_z_star/E_star_p_z_had;
259 x_gluon = (Et_OJ_hcm*exp(-eta_OJ_hcm) + Et_DJ_hcm*exp(-eta_DJ_hcm))/(2.*hcmboost.transform(dk.beamHadron().momentum()).E());
260 //observed fraction of the proton momentum carried by the gluon
261
262 //cout << " xgamm " << x_gamma << " x_glu " << x_gluon << endl;
263 _h["2811"] -> fill(x_gamma);
264 _h_binned["Q2xgam"].fill(Q2, x_gamma);
265
266 _h["3011"] -> fill(log10(x_gluon));
267 _h_binned["Q2xglue"].fill(Q2, log10(x_gluon));
268 }
269 }
270 }
271
272
273 /// Normalise histograms etc., after the run
274 void finalize() {
275 double norm = crossSection()/nanobarn/sumW();
276 //double norm = crossSection()/nanobarn/sumOfWeights();
277 //cout << " SumOfWeights " << sumW() << " "<< sumOfWeights() << endl;
278
279 scale(_h["411"], norm);
280 scale(_h["511"], norm);
281 scale(_h["611"], norm);
282 scale(_h["711"], norm);
283 scale(_h["811"], norm);
284 scale(_h["911"], norm);
285 _h_binned["Q2xbj"].scale(norm, this);
286 scale(_h["1511"], norm);
287 scale(_h["1611"], norm);
288 scale(_h["1711"], norm);
289 scale(_h["1811"], norm);
290 scale(_h["1911"], norm);
291 scale(_h["2011"], norm);
292 scale(_h["2111"], norm);
293 scale(_h["2211"], norm);
294 _h_binned["Q2phi"].scale(norm*180./M_PI, this);
295 scale(_h["2511"], norm);
296 scale(_h["2611"], norm);
297 scale(_h["2711"], norm);
298 scale(_h["2811"], norm);
299 _h_binned["Q2xgam"].scale(norm, this);
300
301 scale(_h["3011"], norm);
302 // new scaling needed, since x bins are in log10(x)
303 vector<YODA::HistoBin1D>& bins = _h["3011"] -> bins();
304 for (auto & b : bins) {
305 double scale_new = b.xWidth()/pow(10,b.xWidth()) ;
306 // jacobian d log10/dx = dlog10 / dlogx dlogx/dx = 2.3026 /x
307 double factor = pow(10,b.xMid())/2.3026;
308 b.scaleW(scale_new/factor) ;
309 }
310
311 _h_binned["Q2xglue"].scale(norm, this);
312 // new scaling needed, since x bins are in log10(x)
313 for (Histo1DPtr histo : _h_binned["Q2xglue"].histos()) {
314 vector<YODA::HistoBin1D>& bins = histo -> bins();
315 for (auto & b : bins) {
316 double scale_new = b.xWidth()/pow(10,b.xWidth()) ;
317 // jacobian d log10/dx = dlog10 / dlogx dlogx/dx = 2.3026 /x
318 double factor = pow(10,b.xMid())/2.3026;
319 b.scaleW(scale_new/factor) ;
320 }
321
322 }
323
324 }
325
326 ///@}
327
328
329 /// @name Histograms
330 ///@{
331 map<string, Histo1DPtr> _h;
332 map<string, Profile1DPtr> _p;
333 map<string, CounterPtr> _c;
334 map<string, BinnedHistogram> _h_binned;
335 ///@}
336
337
338 };
339
340
341 RIVET_DECLARE_PLUGIN(H1_2007_I736052);
342
343}
|