Rivet analyses referenceH1_2006_I699835Measurement of Event Shape Variables in Deep-Inelastic Scattering at HERAExperiment: H1 (HERA) Inspire ID: 699835 Status: VALIDATED Authors:
Beam energies: (27.6, 820.0); (27.6, 920.0) GeV
Deep-inelastic ep scattering data taken with the H1 detector at HERA and corresponding to an integrated luminosity of $106 pb^{-1}$ are used to study the differential distributions of event shape variables. These include thrust, jet broadening, jet mass and the C-parameter. The four-momentum transfer Q is taken to be the relevant energy scale and ranges between 14 GeV and 200 GeV. Source code: H1_2006_I699835.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/PromptFinalState.hh"
5#include "Rivet/Projections/DISKinematics.hh"
6#include "Rivet/Projections/DISFinalState.hh"
7#include "Rivet/Projections/Thrust.hh"
8
9namespace Rivet {
10
11
12
13 /// @brief Measurement of Event Shape Variables in Deep-Inelastic Scattering at HERA (H1)
14 class H1_2006_I699835 : public Analysis {
15 public:
16
17 /// Constructor
18 RIVET_DEFAULT_ANALYSIS_CTOR(H1_2006_I699835);
19
20
21 /// @name Analysis methods
22 ///@{
23
24 /// Book histograms and initialise projections before the run
25 void init() {
26
27 //cerr << endl << endl << "Initializing --------------------------------" << endl << endl;
28
29 // Initialise and register projections
30
31 // The basic final-state projection:
32 // all final-state particles within
33 // the given eta acceptance
34 const FinalState fs(Cuts::abseta < 4.9);
35 declare(fs, "FS");
36
37 const DISFinalState DISfs(DISFinalState::BoostFrame::BREIT);
38
39 const FinalState DISfsCut(DISfs, Cuts::eta < 0);
40
41 declare(Thrust(DISfsCut), "ThrustCut");
42
43 declare(DISKinematics(), "Kinematics");
44
45 //Book histograms for different Q ranges:
46 book(_Nevt_after_cuts, "TMP/Nevt_after_cuts");
47
48 Histo1DPtr dummy;
49 for(int iQ = 0; iQ < 7; ++iQ) {
50 // cout << " iQ " << iQ << " " << QEdges[iQ] << " " << QEdges[iQ+1] << endl;
51 _h_tauc.add(QEdges[iQ], QEdges[iQ+1], book(dummy,1+iQ,1,1));
52 _h_tau.add(QEdges[iQ], QEdges[iQ+1], book(dummy,8+iQ,1,1));
53 _h_B.add(QEdges[iQ], QEdges[iQ+1], book(dummy,15+iQ,1,1));
54 _h_rho.add(QEdges[iQ], QEdges[iQ+1], book(dummy,22+iQ,1,1));
55 book(_Nevt_after_cuts_Q[iQ], "TMP/Nevt_after_cuts_Q"+ to_string(iQ));
56 }
57
58 }
59
60
61 /// Perform the per-event analysis
62 void analyze(const Event& event) {
63
64
65 const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
66
67 // The kinematic region covered by the analysis is defined by ranges of Q² and y
68 if (dk.Q2() < 196 or dk.Q2() > 40000 or dk.y() < 0.1 or dk.y() > 0.7) {
69 // cerr << "Event out of the kinematic region covered by the analysis" << endl;
70 vetoEvent;
71 }
72
73 double Q=sqrt(dk.Q2());
74 _Nevt_after_cuts -> fill();
75 for(int iQ = 0; iQ < 7; ++iQ) {
76 if (inRange(Q, QEdges[iQ],QEdges[iQ+1])) {
77 // cout << " Q " << Q << " " << iQ << endl;
78 _Nevt_after_cuts_Q[iQ] -> fill(); }
79 }
80
81
82 // Boost to hadronic Breit frame
83 const LorentzTransform breitboost = dk.boostBreit();
84
85 const FinalState& fs = apply<FinalState>(event, "FS");
86
87 /*Calculate event shape variables:
88 thrust_num is \sum |pz_h|
89 thrust_den is \sum |p_h|
90 b_num is \sum |pt_h|
91 sumE is the sum of energies
92 (All sums run through the particles in the current hemisphere)
93 */
94 double thrust_num, thrust_den, b_num, sumE;
95 thrust_num = thrust_den = b_num = sumE = 0;
96 Vector3 sumMom;
97
98 for (const Particle& p : fs.particles()) {
99 // Boost to Breit frame
100 const FourMomentum BreitMom = breitboost.transform(p.momentum());
101 if (BreitMom.eta() < 0) {
102 thrust_num += abs(BreitMom.pz());
103 thrust_den += BreitMom.p();
104 b_num += abs(BreitMom.pt());
105 sumMom.operator+=(BreitMom.p3());
106 sumE += BreitMom.E();
107 }
108 }
109
110 //The energy in the current hemisphere must exceed a certain value.
111 if (sumE <= Q/10.0) {
112 // cerr << "Energy in the current hemisphere too low" << endl;
113 vetoEvent;
114 }
115
116 /* Comment from A. Galvan:
117 Thrust here is with respect to the z axis (tau in the paper), while the
118 one from Rivet projection is with respect to the maximum thrust
119 axis (1 - tau_c in the paper)
120 */
121
122 double thrust_mine = 1.0 - ((double)thrust_num)/((double)thrust_den);
123 double b_mine = ((double)b_num)/(2.0*(double)thrust_den);
124 double rho_num = thrust_den*thrust_den - sumMom.dot(sumMom);
125 double rho_mine = ((double)rho_num)/(4.0*(double)thrust_den*(double)thrust_den);
126
127 const Thrust& thrCut = applyProjection<Thrust>(event, "ThrustCut");
128
129 //Fill histograms:
130 _h_tauc.fill(Q, 1.0 - thrCut.thrust());
131 _h_tau.fill(Q, thrust_mine);
132 _h_B.fill(Q, b_mine);
133 _h_rho.fill(Q, rho_mine);
134
135 /* Comment from A. Galvan:
136 As for the C-parameter, my results did not fit the reference data at
137 all. The formula given in the paper is a bit ambiguous, because there is
138 a sum that runs through all pairs of particles h,h' and I was not sure
139 whether each pair should be counted twice (h,h' and h',h) or not, or if
140 the pair of a particle with itself (hh) should be considered.
141 That is why I tried all of the possibilities, but none of them worked.
142 */
143 }
144
145
146 /// Normalise histograms etc., after the run
147 void finalize() {
148
149 //cerr << endl << endl << "Finalizing --------------------------------" << endl << endl;
150
151 // cerr << "Cross section: " << crossSection() << endl;
152
153 // double Nev = dbl(*_Nevt_after_cuts) ;
154 // cout << " N total " << Nev << endl;
155
156
157 int iQ=0;
158 for (Histo1DPtr histo : _h_tauc.histos()) {
159 double Nev = dbl(*_Nevt_after_cuts_Q[iQ]) ;
160 // cout << " Nev " << Nev << " " << iQ << endl;
161 if (Nev != 0) scale(histo, 1./Nev);
162 vector<YODA::HistoBin1D>& bins = histo -> bins();
163 for (auto & b : bins) b.scaleW(b.xWidth());
164 ++iQ;
165 }
166
167 iQ=0;
168 for (Histo1DPtr histo : _h_tau.histos()) {
169 double Nev = dbl(*_Nevt_after_cuts_Q[iQ]) ;
170 //cout << " Nev " << Nev << " " << iQ << endl;
171 if (Nev != 0) scale(histo, 1./Nev);
172 vector<YODA::HistoBin1D>& bins = histo -> bins();
173 for (auto & b : bins) b.scaleW(b.xWidth());
174 ++iQ;
175 }
176
177 iQ = 0;
178 for (Histo1DPtr histo : _h_B.histos()) {
179 double Nev = dbl(*_Nevt_after_cuts_Q[iQ]) ;
180 // cout << " Nev " << Nev << " " << iQ << endl;
181 if (Nev != 0) scale(histo, 1./Nev);
182 vector<YODA::HistoBin1D>& bins = histo -> bins();
183 for (auto & b : bins) b.scaleW(b.xWidth());
184 ++iQ;
185 }
186
187 iQ = 0;
188 for (Histo1DPtr histo : _h_rho.histos()) {
189 double Nev = dbl(*_Nevt_after_cuts_Q[iQ]) ;
190 // cout << " Nev " << Nev << " " << iQ << endl;
191 if (Nev != 0) scale(histo, 1./Nev);
192 vector<YODA::HistoBin1D>& bins = histo -> bins();
193 for (auto & b : bins) b.scaleW(b.xWidth());
194 ++iQ;
195 }
196
197 }
198
199 ///@}
200
201
202 /// @name Histograms
203 ///@{
204 map<string, Histo1DPtr> _h;
205 map<string, Profile1DPtr> _p;
206 map<string, CounterPtr> _c;
207 BinnedHistogram _h_tauc, _h_tau, _h_B, _h_rho;
208 CounterPtr _Nevt_after_cuts;
209 CounterPtr _Nevt_after_cuts_Q[7];
210 ///@}
211 const vector<double> QEdges {14., 16., 20., 30., 50., 70., 100., 200.};
212
213
214 };
215
216
217 RIVET_DECLARE_PLUGIN(H1_2006_I699835);
218
219}
|