Rivet analyses referenceDELPHI_2003_I620250Measurements of event shapes by DELPHI, above and below $m_Z$Experiment: DELPHI (LEP) Inspire ID: 620250 Status: VALIDATED Authors:
Beam energies: (22.5, 22.5); (33.0, 33.0); (38.0, 38.0); (91.5, 91.5); (94.5, 94.5); (96.0, 96.0); (98.0, 98.0); (100.0, 100.0); (101.0, 101.0); (102.5, 102.5); (103.5, 103.5) GeV Run details:
Measurement of a wide range of event shapes by DELPHI at energies below the Z pole using radiative events and above $m_Z$ from LEP2. This analyses allows the energy dependence of simulations to be studied. Only the distributions and not the means are implemented. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: DELPHI_2003_I620250.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/Thrust.hh"
6#include "Rivet/Projections/Sphericity.hh"
7#include "Rivet/Projections/Hemispheres.hh"
8#include "Rivet/Projections/ParisiTensor.hh"
9
10namespace Rivet {
11
12
13 /// @brief DELPHI event shapes below the Z pole
14 class DELPHI_2003_I620250 : public Analysis {
15 public:
16
17 /// Constructor
18 RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_2003_I620250);
19
20
21 /// @name Analysis methods
22 /// @{
23
24 /// Book histograms and initialise projections before the run
25 void init() {
26
27 // Initialise and register projections.
28 declare(Beam(), "Beams");
29 const FinalState fs;
30 declare(fs, "FS");
31 const Thrust thrust(fs);
32 declare(thrust, "Thrust");
33 declare(Sphericity(fs), "Sphericity");
34 declare(ParisiTensor(fs), "Parisi");
35 declare(Hemispheres(thrust), "Hemispheres");
36
37 // Histogram booking offset numbers.
38 unsigned int offset = 0;
39 int offset2 = -1;
40 isDisc = false;
41 skipBin = false;
42 if (isCompatibleWithSqrtS(45*GeV)) offset = 1;
43 else if (isCompatibleWithSqrtS(66*GeV)) offset = 2;
44 else if (isCompatibleWithSqrtS(76*GeV)) offset = 3;
45 else if (isCompatibleWithSqrtS(183*GeV)) {
46 offset2= 0;
47 offset = 1;
48 }
49 else if (isCompatibleWithSqrtS(189*GeV)) {
50 offset2= 0;
51 offset = 2;
52 }
53 else if (isCompatibleWithSqrtS(192*GeV)) {
54 offset2= 0;
55 offset = 3;
56 }
57 else if (isCompatibleWithSqrtS(196*GeV)) {
58 offset2= 0;
59 offset = 4;
60 }
61 else if (isCompatibleWithSqrtS(200*GeV)) {
62 offset2= 1;
63 offset = 1;
64 skipBin = true;
65 }
66 else if (isCompatibleWithSqrtS(202*GeV)) {
67 offset2= 1;
68 offset = 2;
69 skipBin = true;
70 }
71 else if (isCompatibleWithSqrtS(205*GeV)) {
72 offset2= 1;
73 offset = 3;
74 skipBin = true;
75 }
76 else if (isCompatibleWithSqrtS(207*GeV)) {
77 offset2= 1;
78 offset = 4;
79 skipBin = true;
80 }
81 else MSG_ERROR("Beam energy not supported!");
82 // Book the histograms
83 if (offset2 < 0) {
84 book(_h["thrust"], 1, 1, offset);
85 book(_h["major"], 2, 1, offset);
86 book(_h["minor"], 3, 1, offset);
87 book(_h["sphericity"], 4, 1, offset);
88 book(_h["planarity"], 5, 1, offset);
89 book(_h["oblateness"], 6, 1, offset);
90 book(_h["heavy_jet_mass"], 7, 1, offset);
91 book(_h["light_jet_mass"], 9, 1, offset);
92 book(_h["diff_jet_mass"], 10, 1, offset);
93 book(_h["total_jet_mass"], 11, 1, offset);
94 book(_h["heavy_jet_mass_E"], 8, 1, offset);
95 book(_h["total_jet_mass_E"], 12, 1, offset);
96 book(_h["wide_broading"], 13, 1, offset);
97 book(_h["narrow_broading"], 14, 1, offset);
98 book(_h["total_broading"], 15, 1, offset);
99 book(_h["diff_broading"], 16, 1, offset);
100 book(_h["CParam"], 17, 1, offset);
101 }
102 else {
103 isDisc = true;
104 book(_d["rap"], 30+offset2, 1, offset);
105 book(_d["xi"], 32+offset2, 1, offset);
106 book(_d["pTIn"], 34+offset2, 1, offset);
107 book(_d["pTOut"], 36+offset2, 1, offset);
108 book(_d["thrust"], 38+offset2, 1, offset);
109 book(_d["major"], 40+offset2, 1, offset);
110 book(_d["minor"], 42+offset2, 1, offset);
111 book(_d["oblateness"], 44+offset2, 1, offset);
112 book(_d["wide_broading"], 46+offset2, 1, offset);
113 book(_d["total_broading"], 48+offset2, 1, offset);
114 book(_d["diff_broading"], 50+offset2, 1, offset);
115 book(_d["CParam"], 52+offset2, 1, offset);
116 book(_d["DParam"], 54+offset2, 1, offset);
117 book(_d["heavy_jet_mass"], 56+offset2, 1, offset);
118 book(_d["heavy_jet_mass_P"], 58+offset2, 1, offset);
119 book(_d["heavy_jet_mass_E"], 60+offset2, 1, offset);
120 book(_d["light_jet_mass"], 62+offset2, 1, offset);
121 book(_d["diff_jet_mass"], 64+offset2, 1, offset);
122 book(_d["sphericity"], 66+offset2, 1, offset);
123 book(_d["planarity"], 68+offset2, 1, offset);
124 book(_d["aplanarity"], 70+offset2, 1, offset);
125
126 _axis["rap"] = YODA::Axis<double>({0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.5});
127 _axis["xi"] = YODA::Axis<double>({0.0, 0.4, 0.8, 1.2, 1.6, 2.0, 2.4, 2.8, 3.2, 3.6, 4.0,
128 4.4, 4.8, 5.2, 5.6, 6.0, 6.4});
129 _axis["pTIn"] = YODA::Axis<double>({0.0, 0.1, 0.4, 0.65, 0.9, 1.1, 1.4, 2.0, 3.0, 4.0, 6.0, 8.0, 12.0});
130 _axis["pTOut"] = YODA::Axis<double>({0.0, 0.2, 0.4, 0.6, 0.85, 1.2, 1.6, 2.0, 3.0});
131 _axis["thrust"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09,
132 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.24, 0.28, 0.32, 0.36});
133 _axis["major"] = YODA::Axis<double>({0.0, 0.02, 0.04, 0.05, 0.06, 0.07, 0.08, 0.1, 0.12, 0.14, 0.16,
134 0.2, 0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48, 0.52, 0.56, 0.6});
135 _axis["minor"] = YODA::Axis<double>({0.0, 0.02, 0.04, 0.05, 0.06, 0.07, 0.08, 0.1, 0.12, 0.14, 0.16, 0.2, 0.24, 0.28, 0.32});
136 _axis["oblateness"] = YODA::Axis<double>({0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14, 0.16,
137 0.18, 0.20, 0.24, 0.28, 0.32, 0.36, 0.4, 0.44});
138 _axis["wide_broading"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07,
139 0.08, 0.1, 0.12, 0.14, 0.17, 0.20, 0.24, 0.28});
140 _axis["total_broading"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1,
141 0.11, 0.13, 0.15, 0.17, 0.19, 0.21, 0.24, 0.27, 0.3, 0.33, 0.36});
142 _axis["diff_broading"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08,
143 0.09, 0.1, 0.12, 0.14, 0.16, 0.18, 0.2, 0.24, 0.28});
144 _axis["CParam"] = YODA::Axis<double>({0.0, 0.04, 0.08, 0.12, 0.16, 0.2, 0.24, 0.28, 0.32, 0.36, 0.4, 0.44,
145 0.48, 0.52, 0.56, 0.6, 0.64, 0.68, 0.72, 0.76, 0.8, 0.84, 0.88});
146 _axis["DParam"] = YODA::Axis<double>({0.00, 0.02,0.04,0.06,0.08,0.10,0.12,0.14,0.16,0.20,
147 0.24,0.28,0.32,0.36,0.40,0.44,0.48,0.54});
148 _axis["heavy_jet_mass"] = YODA::Axis<double>({0.00, 0.01,0.02,0.03,0.04,0.05,0.06,
149 0.08,0.10,0.12,0.14,0.16,0.20,0.24,0.28,0.32});
150 _axis["heavy_jet_mass_P"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08,
151 0.1, 0.12, 0.14, 0.16, 0.2, 0.24, 0.28, 0.32});
152 _axis["heavy_jet_mass_E"] = YODA::Axis<double>({0.00, 0.01,0.02,0.03,0.04,0.05,0.06,
153 0.08,0.10,0.12,0.14,0.16,0.20,0.24,0.28,0.32});
154 _axis["light_jet_mass"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05});
155 _axis["diff_jet_mass"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.06, 0.08, 0.12, 0.16, 0.2, 0.25, 0.3});
156 _axis["sphericity"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1,
157 0.12, 0.16, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6});
158 _axis["planarity"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08,
159 0.1, 0.12, 0.16, 0.2, 0.25, 0.3, 0.35, 0.4});
160 _axis["aplanarity"] = YODA::Axis<double>({0.0, 0.004, 0.01, 0.016, 0.02, 0.03, 0.04, 0.06, 0.08, 0.1, 0.12, 0.16});
161 }
162 }
163
164
165 /// Perform the per-event analysis
166 void analyze(const Event& event) {
167
168 if (_edges.empty()) {
169 for (const auto& item : _axis) {
170 _edges[item.first] = _d[item.first]->xEdges();
171 }
172 }
173
174 // Get beams and average beam momentum
175 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
176 const double meanBeamMom = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0;
177 MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
178
179 const Thrust& thrust = apply<Thrust>(event, "Thrust");
180 // sphericity related
181 const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
182 // hemisphere related
183 const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
184 smartfill("thrust", 1.-thrust.thrust());
185 smartfill("major", thrust.thrustMajor());
186 smartfill("minor", thrust.thrustMinor());
187 smartfill("oblateness", thrust.oblateness() );
188 smartfill("sphericity", sphericity.sphericity());
189 smartfill("planarity", sphericity.planarity() );
190 if(isDisc) smartfill("aplanarity", sphericity.aplanarity());
191 smartfill("heavy_jet_mass", hemi.scaledM2high());
192 smartfill("light_jet_mass", hemi.scaledM2low() );
193 smartfill("diff_jet_mass", hemi.scaledM2diff());
194 smartfill("wide_broading", hemi.Bmax());
195 if(!isDisc) smartfill("narrow_broading", hemi.Bmin());
196 smartfill("total_broading", hemi.Bsum());
197 smartfill("diff_broading", hemi.Bdiff());
198 if(!isDisc) smartfill("total_jet_mass", hemi.scaledM2low()+hemi.scaledM2high());
199 // E and p scheme jet masses
200 Vector3 axis = thrust.thrustAxis();
201 FourMomentum p4WithE, p4AgainstE;
202 FourMomentum p4WithP, p4AgainstP;
203 double Evis(0);
204 for (const Particle& p : apply<FinalState>(event, "FS").particles()) {
205 Vector3 p3 = p.momentum().vector3().unitVec();
206 const double E = p.momentum().E();
207 Evis += E;
208 p3 = E*p3;
209 const double p3Para = dot(p3, axis);
210 FourMomentum p4E(E,p3.x(),p3.y(),p3.z());
211 FourMomentum p4P(p.p3().mod(),p.p3().x(),p.p3().y(),p.p3().z());
212 if (p3Para > 0) {
213 p4WithE += p4E;
214 p4WithP += p4P;
215 }
216 else if (p3Para < 0) {
217 p4AgainstE += p4E;
218 p4AgainstP += p4P;
219 }
220 else {
221 MSG_WARNING("Particle split between hemispheres");
222 p4WithE += 0.5 * p4E;
223 p4AgainstE += 0.5 * p4E;
224 p4WithP += 0.5 * p4P;
225 p4AgainstP += 0.5 * p4P;
226 }
227 }
228 // E scheme
229 const double mass2With_E = p4WithE.mass2()/sqr(Evis);
230 const double mass2Against_E = p4AgainstE.mass2()/sqr(Evis);
231 // fill the histograms
232 smartfill("heavy_jet_mass_E", max(mass2With_E,mass2Against_E));
233 if(!isDisc) smartfill("total_jet_mass_E", mass2With_E+mass2Against_E);
234 // pscheme
235 const double mass2With_P = p4WithP.mass2()/sqr(Evis);
236 const double mass2Against_P = p4AgainstP.mass2()/sqr(Evis);
237 // fill the histograms
238 if(isDisc) smartfill("heavy_jet_mass_P", max(mass2With_P, mass2Against_P));
239
240 MSG_DEBUG("Calculating Parisi params");
241 const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
242 smartfill("CParam", parisi.C());
243 if(isDisc) smartfill("DParam", parisi.D());
244
245 // single particle distributions
246 const FinalState& fs = apply<FinalState>(event, "FS");
247 if (isDisc) {
248 for (const Particle& p : fs.particles()) {
249 if ( ! PID::isCharged(p.pid())) continue;
250 // Get momentum and energy of each particle.
251 const Vector3 mom3 = p.p3();
252 const double energy = p.E();
253
254 // Scaled momenta.
255 const double mom = mom3.mod();
256 const double scaledMom = mom/meanBeamMom;
257 const double logInvScaledMom = -std::log(scaledMom);
258 smartfill("xi", logInvScaledMom);
259
260 // Get momenta components w.r.t. thrust and sphericity.
261 const double momT = dot(thrust.thrustAxis(), mom3);
262 const double pTinT = dot(mom3, thrust.thrustMajorAxis());
263 const double pToutT = dot(mom3, thrust.thrustMinorAxis());
264 smartfill("pTIn", fabs(pTinT/GeV));
265 smartfill("pTOut", fabs(pToutT/GeV));
266
267 // Calculate rapidities w.r.t. thrust and sphericity.
268 const double rapidityT = 0.5 * std::log((energy + momT) / (energy - momT));
269 smartfill("rap", fabs(rapidityT));
270 MSG_TRACE(fabs(rapidityT) << " " << scaledMom/GeV);
271 }
272 }
273 }
274
275 void smartfill(const string& tag, const double value) {
276 if (isDisc) {
277 size_t idx = _axis[tag].index(value);
278 // skip masked bin in wide broadening
279 if(tag=="wide_broading" && skipBin) {
280 if(idx==8) idx=0;
281 else if(idx>8) idx-=1;
282 }
283 if (idx && idx <= _edges[tag].size()) {
284 _d[tag]->fill(_edges[tag][idx-1]);
285 }
286 else {
287 _d[tag]->fill(string("OTHER"));
288 }
289 }
290 else {
291 _h[tag]->fill(value);
292 }
293 }
294
295 /// Normalise histograms etc., after the run
296 void finalize() {
297 normalize(_h);
298 for (auto& item : _d) {
299 if (item.first == "rap" ||
300 item.first == "xi" ||
301 item.first == "pTIn" ||
302 item.first == "pTOut") scale(item.second, 1./sumOfWeights());
303 else normalize(item.second);
304 for(auto & b: item.second->bins()) {
305 size_t idx = b.index();
306 // skip masked bin in wide broadening
307 if(item.first=="wide_broading" && skipBin) {
308 if(idx>=8) ++idx;
309 }
310 b.scaleW(1./_axis[item.first].width(idx));
311 }
312 }
313 }
314
315 /// @}
316
317
318 /// @name Histograms
319 /// @{
320 map<string, Histo1DPtr> _h;
321 map<string, BinnedHistoPtr<string>> _d;
322 map<string, YODA::Axis<double>> _axis;
323 map<string, vector<string>> _edges;
324 bool isDisc;
325 bool skipBin;
326
327 /// @}
328 };
329
330
331 RIVET_DECLARE_PLUGIN(DELPHI_2003_I620250);
332
333
334}
|