Rivet analyses referenceH1_2000_I503947H1 energy flow in DISExperiment: H1 (HERA) Inspire ID: 503947 Status: VALIDATED Authors:
Beam energies: (820.0, 27.5) GeV Run details:
Measurements of transverse energy flow for neutral current deep- inelastic scattering events produced in positron-proton collisions at HERA. The kinematic range covers squared momentum transfers $Q^2$ from 3.2 to 2200 GeV$^2$; the Bjorken scaling variable $x$ from $8 \times 10^{-5}$ to 0.11 and the hadronic mass $W$ from 66 to 233 GeV. The transverse energy flow is measured in the hadronic centre of mass frame and is studied as a function of $Q^2$, $x$, $W$ and pseudorapidity. The behaviour of the mean transverse energy in the central pseudorapidity region and an interval corresponding to the photon fragmentation region are analysed as a function of $Q^2$ and $W$. This analysis is useful for exploring the effect of photon PDFs and for tuning models of parton evolution and treatment of fragmentation and the proton remnant in DIS. Source code: H1_2000_I503947.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Math/Constants.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/DISKinematics.hh"
6
7namespace Rivet {
8
9
10 /// @brief H1 energy flow and charged particle spectra
11 ///
12 /// @author Peter Richardson
13 ///
14 /// Based on the HZTOOL analysis HZ99091
15 class H1_2000_I503947 : public Analysis {
16 public:
17
18 /// Constructor
19 RIVET_DEFAULT_ANALYSIS_CTOR(H1_2000_I503947);
20
21
22 /// @name Analysis methods
23 /// @{
24
25 /// Initialise projections and histograms
26 void init() {
27 // Projections
28 const DISLepton dl;
29 declare(dl, "Lepton");
30 declare(DISKinematics(), "Kinematics");
31 declare(dl.remainingFinalState(), "FS");
32
33 // Histograms and weight vectors for low Q^2 a
34 _histETLowQa.resize(17);
35 for (size_t ix = 0; ix < 17; ++ix) {
36 book(_histETLowQa[ix], ix+1, 1, 1);
37 book(_weightETLowQa[ix], "TMP/ETLowQa" + to_string(ix));
38 }
39
40 // Histograms and weight vectors for high Q^2 a
41 _histETHighQa.resize(7);
42 for (size_t ix = 0; ix < 7; ++ix) {
43 book(_histETHighQa[ix], ix+18, 1, 1);
44 book(_weightETHighQa[ix], "TMP/ETHighQa" + to_string(ix));
45 }
46
47 // Histograms and weight vectors for low Q^2 b
48 _histETLowQb.resize(5);
49 for (size_t ix = 0; ix < 5; ++ix) {
50 book(_histETLowQb[ix], ix+25, 1, 1);
51 book(_weightETLowQb[ix], "TMP/ETLowQb" + to_string(ix));
52 }
53
54 // Histograms and weight vectors for high Q^2 b
55 _histETHighQb.resize(5);
56 for (size_t ix = 0; ix < 3; ++ix) {
57 book(_histETHighQb[ix], 30+ix, 1, 1);
58 book(_weightETHighQb[ix], "TMP/ETHighQb" + to_string(ix));
59 }
60
61 // Histograms for the averages
62 book(_histAverETCentral, 33, 1, 1);
63 book(_histAverETFrag, 34, 1, 1);
64 }
65
66
67 /// Analyze each event
68 void analyze(const Event& event) {
69
70 // DIS kinematics
71 const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
72 if ( dk.failed() ) vetoEvent;
73 double q2 = dk.Q2();
74 double x = dk.x();
75 double y = dk.y();
76 double w2 = dk.W2();
77
78 // Kinematics of the scattered lepton
79 const DISLepton& dl = apply<DISLepton>(event,"Lepton");
80 if ( dl.failed() ) vetoEvent;
81 const FourMomentum leptonMom = dl.out();
82 const double enel = leptonMom.E();
83 const double thel = 180 - leptonMom.angle(dl.in().mom())/degree;
84
85 // Extract the particles other than the lepton
86 Particles particles= apply<FinalState>(event, "FS").particles();
87
88 // Cut on the forward energy
89 double efwd = 0.;
90 for (const Particle& p : particles) {
91 const double th = 180 - p.angle(dl.in())/degree;
92 if (inRange(th, 4.4, 15.0)) efwd += p.E();
93 }
94 // There are four possible selections for events
95 bool evcut[4];
96 // Low Q2 selection a
97 evcut[0] = enel/GeV > 12. && w2 >= 4400.*GeV2 && efwd/GeV > 0.5 && inRange(thel,157.,176.);
98 // Low Q2 selection b
99 evcut[1] = enel/GeV > 12. && inRange(y,0.3,0.5);
100 // High Q2 selection a
101 evcut[2] = inRange(thel,12.,150.) && inRange(y,0.05,0.6) && w2 >= 4400.*GeV2 && efwd > 0.5;
102 // High Q2 selection b
103 evcut[3] = inRange(thel,12.,150.) && inRange(y,0.05,0.6) && inRange(w2,27110.*GeV2,45182.*GeV2);
104
105 // Veto if fails all cuts
106 /// @todo Can we use all()?
107 if (! (evcut[0] || evcut[1] || evcut[2] || evcut[3]) ) vetoEvent;
108
109 // Find the bins
110 int bin[4] = {-1,-1,-1,-1};
111 // For the low Q2 selection a)
112 if (q2 > 2.5*GeV2 && q2 <= 5.*GeV2) {
113 if (x > 0.00005 && x <= 0.0001 ) bin[0] = 0;
114 if (x > 0.0001 && x <= 0.0002 ) bin[0] = 1;
115 if (x > 0.0002 && x <= 0.00035) bin[0] = 2;
116 if (x > 0.00035 && x <= 0.0010 ) bin[0] = 3;
117 }
118 else if (q2 > 5.*GeV2 && q2 <= 10.*GeV2) {
119 if (x > 0.0001 && x <= 0.0002 ) bin[0] = 4;
120 if (x > 0.0002 && x <= 0.00035) bin[0] = 5;
121 if (x > 0.00035 && x <= 0.0007 ) bin[0] = 6;
122 if (x > 0.0007 && x <= 0.0020 ) bin[0] = 7;
123 }
124 else if (q2 > 10.*GeV2 && q2 <= 20.*GeV2) {
125 if (x > 0.0002 && x <= 0.0005) bin[0] = 8;
126 if (x > 0.0005 && x <= 0.0008) bin[0] = 9;
127 if (x > 0.0008 && x <= 0.0015) bin[0] = 10;
128 if (x > 0.0015 && x <= 0.040 ) bin[0] = 11;
129 }
130 else if (q2 > 20.*GeV2 && q2 <= 50.*GeV2) {
131 if (x > 0.0005 && x <= 0.0014) bin[0] = 12;
132 if (x > 0.0014 && x <= 0.0030) bin[0] = 13;
133 if (x > 0.0030 && x <= 0.0100) bin[0] = 14;
134 }
135 else if (q2 > 50.*GeV2 && q2 <= 100.*GeV2) {
136 if (x >0.0008 && x <= 0.0030) bin[0] = 15;
137 if (x >0.0030 && x <= 0.0200) bin[0] = 16;
138 }
139 // check in one of the bins
140 evcut[0] &= bin[0] >= 0;
141 // For the low Q2 selection b)
142 if (q2 > 2.5*GeV2 && q2 <= 5. *GeV2) bin[1] = 0;
143 if (q2 > 5. *GeV2 && q2 <= 10. *GeV2) bin[1] = 1;
144 if (q2 > 10.*GeV2 && q2 <= 20. *GeV2) bin[1] = 2;
145 if (q2 > 20.*GeV2 && q2 <= 50. *GeV2) bin[1] = 3;
146 if (q2 > 50.*GeV2 && q2 <= 100.*GeV2) bin[1] = 4;
147 // check in one of the bins
148 evcut[1] &= bin[1] >= 0;
149 // for the high Q2 selection a)
150 if (q2 > 100.*GeV2 && q2 <= 400.*GeV2) {
151 if (x > 0.00251 && x <= 0.00631) bin[2] = 0;
152 if (x > 0.00631 && x <= 0.0158 ) bin[2] = 1;
153 if (x > 0.0158 && x <= 0.0398 ) bin[2] = 2;
154 }
155 else if (q2 > 400.*GeV2 && q2 <= 1100.*GeV2) {
156 if (x > 0.00631 && x <= 0.0158 ) bin[2] = 3;
157 if (x > 0.0158 && x <= 0.0398 ) bin[2] = 4;
158 if (x > 0.0398 && x <= 1. ) bin[2] = 5;
159 }
160 else if (q2 > 1100.*GeV2 && q2 <= 100000.*GeV2) {
161 if (x > 0. && x <= 1.) bin[2] = 6;
162 }
163 // check in one of the bins
164 evcut[2] &= bin[2] >= 0;
165 // for the high Q2 selection b)
166 if (q2 > 100.*GeV2 && q2 <= 220.*GeV2) bin[3] = 0;
167 else if (q2 > 220.*GeV2 && q2 <= 400.*GeV2) bin[3] = 1;
168 else if (q2 > 400. ) bin[3] = 2;
169 // check in one of*GeV the bins
170 evcut[3] &= bin[3] >= 0;
171
172 // Veto if fails all cuts after bin selection
173 /// @todo Can we use all()?
174 if (! (evcut[0] || evcut[1] || evcut[2] || evcut[3])) vetoEvent;
175
176 // Increment the count for normalisation
177 if (evcut[0]) _weightETLowQa [bin[0]]->fill();
178 if (evcut[1]) _weightETLowQb [bin[1]]->fill();
179 if (evcut[2]) _weightETHighQa[bin[2]]->fill();
180 if (evcut[3]) _weightETHighQb[bin[3]]->fill();
181
182 // Boost to hadronic CoM
183 const LorentzTransform hcmboost = dk.boostHCM();
184
185 // Loop over the particles
186 double etcent = 0;
187 double etfrag = 0;
188 for (const Particle& p : particles) {
189 // Boost momentum to CMS
190 const FourMomentum hcmMom = hcmboost.transform(p.momentum());
191 double et = fabs(hcmMom.Et());
192 double eta = hcmMom.eta();
193 // Averages in central and forward region
194 if (fabs(eta) < 0.5 ) etcent += et;
195 if (eta > 2 && eta <= 3.) etfrag += et;
196 // Histograms of Et flow
197 if (evcut[0]) _histETLowQa [bin[0]]->fill(eta, et);
198 if (evcut[1]) _histETLowQb [bin[1]]->fill(eta, et);
199 if (evcut[2]) _histETHighQa[bin[2]]->fill(eta, et);
200 if (evcut[3]) _histETHighQb[bin[3]]->fill(eta, et);
201 }
202 // Fill histograms for the average quantities
203 if (evcut[1] || evcut[3]) {
204 _histAverETCentral->fill(q2, etcent);
205 _histAverETFrag ->fill(q2, etfrag);
206 }
207 }
208
209
210 // Finalize
211 void finalize() {
212 // Normalization of the Et distributions to unit cross-section
213 // These histograms are filled once per particles though,
214 // so no unit area to be expected
215 /// @todo Simplify by using normalize() instead? Are all these being normalized to area=1?
216 for (size_t ix = 0; ix < _weightETLowQa.size(); ++ix) {
217 if (_weightETLowQa[ix]->val() ) scale(_histETLowQa[ix], 1/ *_weightETLowQa[ix]);
218 }
219 for (size_t ix = 0; ix < _weightETHighQa.size(); ++ix) {
220 if (_weightETHighQa[ix]->val()) scale(_histETHighQa[ix], 1/ *_weightETHighQa[ix]);
221 }
222 for (size_t ix = 0; ix < _weightETLowQb.size(); ++ix) {
223 if (_weightETLowQb[ix]->val() ) scale(_histETLowQb[ix], 1/ *_weightETLowQb[ix]);
224 }
225 for (size_t ix = 0; ix < _weightETHighQb.size(); ++ix) {
226 if (_weightETHighQb[ix]->val()) scale(_histETHighQb[ix], 1/ *_weightETHighQb[ix]);
227 }
228 }
229
230 /// @}
231
232
233 private:
234
235 /// @name Histograms
236 /// @{
237 vector<Histo1DPtr> _histETLowQa;
238 vector<Histo1DPtr> _histETHighQa;
239 vector<Histo1DPtr> _histETLowQb;
240 vector<Histo1DPtr> _histETHighQb;
241 Profile1DPtr _histAverETCentral;
242 Profile1DPtr _histAverETFrag;
243 /// @}
244
245 /// @name Storage of weights for normalisation
246 /// @{
247 array<CounterPtr,17> _weightETLowQa;
248 array<CounterPtr, 7> _weightETHighQa;
249 array<CounterPtr, 5> _weightETLowQb;
250 array<CounterPtr, 3> _weightETHighQb;
251 /// @}
252
253 };
254
255
256
257 RIVET_DECLARE_ALIASED_PLUGIN(H1_2000_I503947, H1_2000_S4129130);
258
259}
|