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 Q2 from 3.2 to 2200 GeV2; the Bjorken scaling variable x from 8×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 Q2, 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 Q2 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}
|
