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#include "Rivet/Projections/DISKinematics.hh"
6
7namespace Rivet {
8
9
10 /// @brief Forward pi0 meson production at HERA (H1)
11 class H1_1999_I504022 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(H1_1999_I504022);
16
17
18 /// @name Analysis methods
19 ///@{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23
24 // Initialise and register projections
25
26 // The basic final-state projection:
27
28 declare(FinalState(Cuts::abseta < 7 ), "FS");
29 declare(DISKinematics(), "Kinematics");
30 declare(UnstableParticles(), "UFS");
31
32
33 // Book histograms
34 // take binning from reference data using HEPData ID (digits in "d01-x01-y01" etc.)
35
36 book(_h["p_T>2.5&x"], 1, 1, 1);
37 book(_h["Q:2.0-4.5-eta"], 2, 1, 1);
38 book(_h["Q:2.0-4.5-p_T"], 3, 1, 1);
39 book(_h["Q:4.5-15.0-x"], 4, 1, 1);
40 book(_h["Q:4.5-15.0-eta"], 5, 1, 1);
41 book(_h["Q:4.5-15.0-p_T"], 6, 1, 1);
42 book(_h["Q:15.0-70.0-x"], 7, 1, 1);
43 book(_h["Q:15.0-70.0-eta"], 8, 1, 1);
44 book(_h["Q:15.0-70.0-p_T"], 9, 1, 1);
45 book(_h["p_T>2.5&Q"], 10, 1, 1);
46 book(_h["p_T>3.5&x"], 11, 1, 1);
47 book(_h["p_T>3.5&Q"], 12, 1, 1);
48 }
49
50
51 /// Perform the per-event analysis
52 void analyze(const Event& event) {
53
54 /// @todo Do the event by event analysis here
55
56 const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
57
58 // Get the DIS kinematics
59 double xbj = dk.x();
60 double ybj = dk.y();
61 double Q2 = dk.Q2()/GeV2;
62
63 // Q2 and inelasticity cuts
64 //cout << " after xbj " << xbj << endl;
65 if (!inRange(ybj, 0.1, 0.6)) vetoEvent;
66 //cout << " after ybj " << ybj << endl;
67 if (!inRange(Q2, 2.0*GeV2, 70.0*GeV2)) vetoEvent;
68 //cout << " after Q2 " << Q2 << endl;
69
70 const FinalState& fs = apply<FinalState>(event, "FS");
71 const size_t numParticles = fs.particles().size();
72 //cout << " Num all final state particles " << numParticles << endl;
73 //proton beam energy: 820 GeV
74 double e_proton = dk.beamHadron().E()/GeV;
75 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
76 if (numParticles < 2) {
77 MSG_DEBUG("Failed leptonic event cut");
78 vetoEvent;
79 }
80
81 // Extracting the pi0
82 const UnstableParticles& ufs = apply<UnstableFinalState>(event, "UFS");
83 //Get the hadronic CMS kinematics
84 const LorentzTransform hcmboost = dk.boostHCM();
85
86 for (const Particle& p : ufs.particles(Cuts::pid==PID::PI0)) {
87 //cout << " Pid = " << p.pid() << endl;
88 //Get the LAB kinematics
89 double theta = p.theta();
90 double eta = p.pseudorapidity();
91 // double pT = p.pT()/GeV;
92 //Boost hcm
93 const FourMomentum hcmMom = hcmboost.transform(p.mom());
94 double pThcm = hcmMom.pT();
95
96 double e_pi0 = p.E()/GeV;
97 double x_pi0_proton = e_pi0/e_proton;
98
99 // epi0/e_proton, theta and pThcm cuts
100 if(x_pi0_proton < 0.01) continue;
101 //cout << " theta " << theta << " in deg " << theta/degree << endl;
102 if (!inRange(theta/degree, 5, 25)) continue;
103 if(pThcm < 2.5*GeV) continue;
104
105 //Three cuts for Q2:
106 if (Q2 > 2.0*GeV2 && Q2 < 4.5*GeV2){
107 _h["Q:2.0-4.5-eta"]->fill(eta);
108 _h["Q:2.0-4.5-p_T"]->fill(pThcm);
109 }
110 if (Q2 > 4.5*GeV2 && Q2 < 15.0*GeV2){
111 _h["Q:4.5-15.0-x"]->fill(xbj);
112 _h["Q:4.5-15.0-eta"]->fill(eta);
113 _h["Q:4.5-15.0-p_T"]->fill(pThcm);
114 }
115 if (Q2 > 15.0*GeV2 && Q2 < 70.0*GeV2){
116 _h["Q:15.0-70.0-x"]->fill(xbj);
117 _h["Q:15.0-70.0-eta"]->fill(eta);
118 _h["Q:15.0-70.0-p_T"]->fill(pThcm);
119 }
120
121 //Two cuts for p_T:
122 if (pThcm > 2.5*GeV){
123 _h["p_T>2.5&x"]->fill(xbj);
124 _h["p_T>2.5&Q"]->fill(Q2);
125 }
126 if (pThcm > 3.5*GeV){
127 _h["p_T>3.5&x"]->fill(xbj);
128 _h["p_T>3.5&Q"]->fill(Q2);
129 }
130
131 }
132
133 }
134
135
136 /// Normalise histograms etc., after the run
137 void finalize() {
138 const double sn = crossSection()/nanobarn/sumW();
139 const double sp = crossSection()/picobarn/sumW();
140 scale(_h["p_T>2.5&x"], sn);
141 scale(_h["Q:2.0-4.5-eta"], sp);
142 scale(_h["Q:2.0-4.5-p_T"], sp);
143 scale(_h["Q:4.5-15.0-x"], sn);
144 scale(_h["Q:4.5-15.0-eta"], sp);
145 scale(_h["Q:4.5-15.0-p_T"], sp);
146 scale(_h["Q:15.0-70.0-x"], sn);
147 scale(_h["Q:15.0-70.0-eta"], sp);
148 scale(_h["Q:15.0-70.0-p_T"], sp);
149 scale(_h["p_T>2.5&Q"], sp);
150 scale(_h["p_T>3.5&x"], sn);
151 scale(_h["p_T>3.5&Q"], sp);
152 }
153
154 ///@}
155
156
157 /// @name Histograms
158 ///@{
159 map<string, Histo1DPtr> _h;
160 ///@}
161
162
163 };
164
165
166 RIVET_DECLARE_PLUGIN(H1_1999_I504022);
167
168}
|