Rivet analyses referenceH1_1999_I504022Forward pi0 meson production at HERAExperiment: H1 (HERA) Inspire ID: 504022 Status: VALIDATED Authors:
Beam energies: (820.0, 27.5); (27.5, 820.0) GeV Run details:
High transverse momentum $\pi^0$ mesons have been measured with the H1 detector at HERA in deep-inelastic ep scattering events at low Bjorken-x, down to $x \approx 4 \cdot 10^{-5}$ .The measurement is performed in a region of small angles with respect to the proton remnant in the laboratory frame of reference, namely the forward region, and corresponds to central rapidity in the centre of mass system of the virtual photon and proton. This region is expected to be particularly sensitive to QCD effects in hadronic final states. Differential cross-sections for inclusive $\pi^0$ meson production are presented as a function of Bjorken-x and the four-momentum transfer $Q^2$, and as a function of transverse momentum and pseudorapidity. Source code: H1_1999_I504022.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/UnstableParticles.hh"
5
6#include "Rivet/Projections/DISKinematics.hh"
7
8namespace Rivet {
9
10
11 /// @brief Forward pi0 meson production at HERA (H1)
12 class H1_1999_I504022 : public Analysis {
13 public:
14
15 /// Constructor
16 RIVET_DEFAULT_ANALYSIS_CTOR(H1_1999_I504022);
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
27 // The basic final-state projection:
28
29 declare(FinalState(Cuts::abseta < 7 ), "FS");
30 declare(DISKinematics(), "Kinematics");
31 declare(UnstableParticles(), "UFS");
32
33
34 // Book histograms
35 // take binning from reference data using HEPData ID (digits in "d01-x01-y01" etc.)
36
37 book(_h["p_T>2.5&x"], 1, 1, 1);
38 book(_h["Q:2.0-4.5-eta"], 2, 1, 1);
39 book(_h["Q:2.0-4.5-p_T"], 3, 1, 1);
40 book(_h["Q:4.5-15.0-x"], 4, 1, 1);
41 book(_h["Q:4.5-15.0-eta"], 5, 1, 1);
42 book(_h["Q:4.5-15.0-p_T"], 6, 1, 1);
43 book(_h["Q:15.0-70.0-x"], 7, 1, 1);
44 book(_h["Q:15.0-70.0-eta"], 8, 1, 1);
45 book(_h["Q:15.0-70.0-p_T"], 9, 1, 1);
46 book(_h["p_T>2.5&Q"], 10, 1, 1);
47 book(_h["p_T>3.5&x"], 11, 1, 1);
48 book(_h["p_T>3.5&Q"], 12, 1, 1);
49 }
50
51
52 /// Perform the per-event analysis
53 void analyze(const Event& event) {
54
55 /// @todo Do the event by event analysis here
56
57 const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
58
59 // Get the DIS kinematics
60 double xbj = dk.x();
61 double ybj = dk.y();
62 double Q2 = dk.Q2()/GeV2;
63
64 // Q2 and inelasticity cuts
65 //cout << " after xbj " << xbj << endl;
66 if (!inRange(ybj, 0.1, 0.6)) vetoEvent;
67 //cout << " after ybj " << ybj << endl;
68 if (!inRange(Q2, 2.0*GeV2, 70.0*GeV2)) vetoEvent;
69 //cout << " after Q2 " << Q2 << endl;
70
71 const FinalState& fs = apply<FinalState>(event, "FS");
72 const size_t numParticles = fs.particles().size();
73 //cout << " Num all final state particles " << numParticles << endl;
74 //proton beam energy: 820 GeV
75 double e_proton = dk.beamHadron().E()/GeV;
76 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
77 if (numParticles < 2) {
78 MSG_DEBUG("Failed leptonic event cut");
79 vetoEvent;
80 }
81
82 // Extracting the pi0
83 const UnstableParticles& ufs = apply<UnstableFinalState>(event, "UFS");
84 //Get the hadronic CMS kinematics
85 const LorentzTransform hcmboost = dk.boostHCM();
86
87 for (const Particle& p : ufs.particles(Cuts::pid==PID::PI0)) {
88 //cout << " Pid = " << p.pid() << endl;
89 //Get the LAB kinematics
90 double theta = p.theta();
91 double eta = p.momentum().pseudorapidity();
92 // double pT = p.momentum().pT()/GeV;
93 //Boost hcm
94 const FourMomentum hcmMom = hcmboost.transform(p.momentum());
95 double pThcm =hcmMom.pT();
96
97 double e_pi0 = p.E()/GeV;
98 double x_pi0_proton = e_pi0/e_proton;
99
100 // epi0/e_proton, theta and pThcm cuts
101 if(x_pi0_proton < 0.01) continue;
102 //cout << " theta " << theta << " in deg " << theta/degree << endl;
103 if (!inRange(theta/degree, 5, 25)) continue;
104 if(pThcm < 2.5*GeV) continue;
105
106 //Three cuts for Q2:
107 if (Q2 > 2.0*GeV2 && Q2 < 4.5*GeV2){
108 _h["Q:2.0-4.5-eta"]->fill(eta);
109 _h["Q:2.0-4.5-p_T"]->fill(pThcm);
110 }
111 if (Q2 > 4.5*GeV2 && Q2 < 15.0*GeV2){
112 _h["Q:4.5-15.0-x"]->fill(xbj);
113 _h["Q:4.5-15.0-eta"]->fill(eta);
114 _h["Q:4.5-15.0-p_T"]->fill(pThcm);
115 }
116 if (Q2 > 15.0*GeV2 && Q2 < 70.0*GeV2){
117 _h["Q:15.0-70.0-x"]->fill(xbj);
118 _h["Q:15.0-70.0-eta"]->fill(eta);
119 _h["Q:15.0-70.0-p_T"]->fill(pThcm);
120 }
121
122 //Two cuts for p_T:
123 if (pThcm > 2.5*GeV){
124 _h["p_T>2.5&x"]->fill(xbj);
125 _h["p_T>2.5&Q"]->fill(Q2);
126 }
127 if (pThcm > 3.5*GeV){
128 _h["p_T>3.5&x"]->fill(xbj);
129 _h["p_T>3.5&Q"]->fill(Q2);
130 }
131
132 }
133
134 }
135
136
137 /// Normalise histograms etc., after the run
138 void finalize() {
139 const double sn = crossSection()/nanobarn/sumW();
140 const double sp = crossSection()/picobarn/sumW();
141 scale(_h["p_T>2.5&x"], sn);
142 scale(_h["Q:2.0-4.5-eta"], sp);
143 scale(_h["Q:2.0-4.5-p_T"], sp);
144 scale(_h["Q:4.5-15.0-x"], sn);
145 scale(_h["Q:4.5-15.0-eta"], sp);
146 scale(_h["Q:4.5-15.0-p_T"], sp);
147 scale(_h["Q:15.0-70.0-x"], sn);
148 scale(_h["Q:15.0-70.0-eta"], sp);
149 scale(_h["Q:15.0-70.0-p_T"], sp);
150 scale(_h["p_T>2.5&Q"], sp);
151 scale(_h["p_T>3.5&x"], sn);
152 scale(_h["p_T>3.5&Q"], sp);
153 }
154
155 ///@}
156
157
158 /// @name Histograms
159 ///@{
160 map<string, Histo1DPtr> _h;
161 map<string, Profile1DPtr> _p;
162 map<string, CounterPtr> _c;
163 ///@}
164
165
166 };
167
168
169 RIVET_DECLARE_PLUGIN(H1_1999_I504022);
170
171}
|