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
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 }
65 else if (isCompatibleWithSqrtS(202*GeV)) {
66 offset2= 1;
67 offset = 2;
68 }
69 else if (isCompatibleWithSqrtS(205*GeV)) {
70 offset2= 1;
71 offset = 3;
72 }
73 else if (isCompatibleWithSqrtS(207*GeV)) {
74 offset2= 1;
75 offset = 4;
76 }
77 else MSG_ERROR("Beam energy not supported!");
78 // Book the histograms
79 if (offset2 < 0) {
80 book(_h["thrust"], 1, 1, offset);
81 book(_h["major"], 2, 1, offset);
82 book(_h["minor"], 3, 1, offset);
83 book(_h["sphericity"], 4, 1, offset);
84 book(_h["planarity"], 5, 1, offset);
85 book(_h["oblateness"], 6, 1, offset);
86 book(_h["heavy_jet_mass"], 7, 1, offset);
87 book(_h["light_jet_mass"], 9, 1, offset);
88 book(_h["diff_jet_mass"], 10, 1, offset);
89 book(_h["total_jet_mass"], 11, 1, offset);
90 book(_h["heavy_jet_mass_E"], 8, 1, offset);
91 book(_h["total_jet_mass_E"], 12, 1, offset);
92 book(_h["wide_broading"], 13, 1, offset);
93 book(_h["narrow_broading"], 14, 1, offset);
94 book(_h["total_broading"], 15, 1, offset);
95 book(_h["diff_broading"], 16, 1, offset);
96 book(_h["CParam"], 17, 1, offset);
97 }
98 else {
99 isDisc = true;
100 book(_d["rap"], 30+offset2, 1, offset);
101 book(_d["xi"], 32+offset2, 1, offset);
102 book(_d["pTIn"], 34+offset2, 1, offset);
103 book(_d["pTOut"], 36+offset2, 1, offset);
104 book(_d["thrust"], 38+offset2, 1, offset);
105 book(_d["major"], 40+offset2, 1, offset);
106 book(_d["minor"], 42+offset2, 1, offset);
107 book(_d["oblateness"], 44+offset2, 1, offset);
108 book(_d["wide_broading"], 46+offset2, 1, offset);
109 book(_d["total_broading"], 48+offset2, 1, offset);
110 book(_d["diff_broading"], 50+offset2, 1, offset);
111 book(_d["CParam"], 52+offset2, 1, offset);
112 book(_d["DParam"], 54+offset2, 1, offset);
113 book(_d["heavy_jet_mass"], 56+offset2, 1, offset);
114 book(_d["heavy_jet_mass_P"], 58+offset2, 1, offset);
115 book(_d["heavy_jet_mass_E"], 60+offset2, 1, offset);
116 book(_d["light_jet_mass"], 62+offset2, 1, offset);
117 book(_d["diff_jet_mass"], 64+offset2, 1, offset);
118 book(_d["sphericity"], 66+offset2, 1, offset);
119 book(_d["planarity"], 68+offset2, 1, offset);
120 book(_d["aplanarity"], 70+offset2, 1, offset);
121
122 _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});
123 _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,
124 4.4, 4.8, 5.2, 5.6, 6.0, 6.4});
125 _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});
126 _axis["pTOut"] = YODA::Axis<double>({0.0, 0.2, 0.4, 0.6, 0.85, 1.2, 1.6, 2.0, 3.0});
127 _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,
128 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.24, 0.28, 0.32, 0.36});
129 _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,
130 0.2, 0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48, 0.52, 0.56, 0.6});
131 _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});
132 _axis["oblateness"] = YODA::Axis<double>({0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14, 0.16,
133 0.18, 0.20, 0.24, 0.28, 0.32, 0.36, 0.4, 0.44});
134 _axis["wide_broading"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07,
135 0.08, 0.1, 0.12, 0.14, 0.17, 0.20, 0.24, 0.28});
136 _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,
137 0.11, 0.13, 0.15, 0.17, 0.19, 0.21, 0.24, 0.27, 0.3, 0.33, 0.36});
138 _axis["diff_broading"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08,
139 0.09, 0.1, 0.12, 0.14, 0.16, 0.18, 0.2, 0.24, 0.28});
140 _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,
141 0.48, 0.52, 0.56, 0.6, 0.64, 0.68, 0.72, 0.76, 0.8, 0.84, 0.88});
142 _axis["DParam"] = 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,
143 0.48, 0.52, 0.56, 0.6, 0.64, 0.68, 0.72, 0.76, 0.8, 0.84, 0.88, 0.92});
144 _axis["heavy_jet_mass"] = YODA::Axis<double>({0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14, 0.16,
145 0.2, 0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48, 0.54});
146 _axis["heavy_jet_mass_P"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08,
147 0.1, 0.12, 0.14, 0.16, 0.2, 0.24, 0.28, 0.32});
148 _axis["heavy_jet_mass_E"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05});
149 _axis["light_jet_mass"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05});
150 _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});
151 _axis["sphericity"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1,
152 0.12, 0.16, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6});
153 _axis["planarity"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08,
154 0.1, 0.12, 0.16, 0.2, 0.25, 0.3, 0.35, 0.4});
155 _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});
156 }
157 }
158
159
160 /// Perform the per-event analysis
161 void analyze(const Event& event) {
162
163 if (_edges.empty()) {
164 for (const auto& item : _axis) {
165 _edges[item.first] = _d[item.first]->xEdges();
166 }
167 }
168
169 // Get beams and average beam momentum
170 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
171 const double meanBeamMom = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0;
172 MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
173
174 const Thrust& thrust = apply<Thrust>(event, "Thrust");
175 // sphericity related
176 const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
177 // hemisphere related
178 const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
179 smartfill("thrust", 1.-thrust.thrust());
180 smartfill("major", thrust.thrustMajor());
181 smartfill("minor", thrust.thrustMinor());
182 smartfill("oblateness", thrust.oblateness() );
183 smartfill("sphericity", sphericity.sphericity());
184 smartfill("planarity", sphericity.planarity() );
185 smartfill("aplanarity", sphericity.aplanarity());
186 smartfill("heavy_jet_mass", hemi.scaledM2high());
187 smartfill("light_jet_mass", hemi.scaledM2low() );
188 smartfill("diff_jet_mass", hemi.scaledM2diff());
189 smartfill("wide_broading", hemi.Bmax());
190 smartfill("narrow_broading", hemi.Bmin());
191 smartfill("total_broading", hemi.Bsum());
192 smartfill("diff_broading", hemi.Bdiff());
193 smartfill("total_jet_mass", hemi.scaledM2low()+hemi.scaledM2high());
194 // E and p scheme jet masses
195 Vector3 axis = thrust.thrustAxis();
196 FourMomentum p4WithE, p4AgainstE;
197 FourMomentum p4WithP, p4AgainstP;
198 double Evis(0);
199 for (const Particle& p : apply<FinalState>(event, "FS").particles()) {
200 Vector3 p3 = p.momentum().vector3().unitVec();
201 const double E = p.momentum().E();
202 Evis += E;
203 p3 = E*p3;
204 const double p3Para = dot(p3, axis);
205 FourMomentum p4E(E,p3.x(),p3.y(),p3.z());
206 FourMomentum p4P(p.p3().mod(),p.p3().x(),p.p3().y(),p.p3().z());
207 if (p3Para > 0) {
208 p4WithE += p4E;
209 p4WithP += p4P;
210 }
211 else if (p3Para < 0) {
212 p4AgainstE += p4E;
213 p4AgainstP += p4P;
214 }
215 else {
216 MSG_WARNING("Particle split between hemispheres");
217 p4WithE += 0.5 * p4E;
218 p4AgainstE += 0.5 * p4E;
219 p4WithP += 0.5 * p4P;
220 p4AgainstP += 0.5 * p4P;
221 }
222 }
223 // E scheme
224 const double mass2With_E = p4WithE.mass2()/sqr(Evis);
225 const double mass2Against_E = p4AgainstE.mass2()/sqr(Evis);
226 // fill the histograms
227 smartfill("heavy_jet_mass_E", max(mass2With_E,mass2Against_E));
228 smartfill("total_jet_mass_E", mass2With_E+mass2Against_E);
229 // pscheme
230 const double mass2With_P = p4WithP.mass2()/sqr(Evis);
231 const double mass2Against_P = p4AgainstP.mass2()/sqr(Evis);
232 // fill the histograms
233 smartfill("heavy_jet_mass_P", max(mass2With_P, mass2Against_P));
234
235 MSG_DEBUG("Calculating Parisi params");
236 const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
237 smartfill("CParam", parisi.C());
238 smartfill("DParam", parisi.D());
239
240 // single particle distributions
241 const FinalState& fs = apply<FinalState>(event, "FS");
242 if (isDisc) {
243 for (const Particle& p : fs.particles()) {
244 if ( ! PID::isCharged(p.pid())) continue;
245 // Get momentum and energy of each particle.
246 const Vector3 mom3 = p.p3();
247 const double energy = p.E();
248
249 // Scaled momenta.
250 const double mom = mom3.mod();
251 const double scaledMom = mom/meanBeamMom;
252 const double logInvScaledMom = -std::log(scaledMom);
253 smartfill("xi", logInvScaledMom);
254
255 // Get momenta components w.r.t. thrust and sphericity.
256 const double momT = dot(thrust.thrustAxis(), mom3);
257 const double pTinT = dot(mom3, thrust.thrustMajorAxis());
258 const double pToutT = dot(mom3, thrust.thrustMinorAxis());
259 smartfill("pTIn", fabs(pTinT/GeV));
260 smartfill("pTOut", fabs(pToutT/GeV));
261
262 // Calculate rapidities w.r.t. thrust and sphericity.
263 const double rapidityT = 0.5 * std::log((energy + momT) / (energy - momT));
264 smartfill("rap", fabs(rapidityT));
265 MSG_TRACE(fabs(rapidityT) << " " << scaledMom/GeV);
266 }
267 }
268 }
269
270 void smartfill(const string& tag, const double value) {
271 if (isDisc) {
272 const size_t idx = _axis[tag].index(value);
273 if (idx && idx <= _edges[tag].size()) {
274 _d[tag]->fill(_edges[tag][idx-1]);
275 }
276 else {
277 _d[tag]->fill(string("OTHER"));
278 }
279 }
280 else {
281 _h[tag]->fill(value);
282 }
283 }
284
285 /// Normalise histograms etc., after the run
286 void finalize() {
287 normalize(_h);
288 for (auto& item : _d) {
289 if (item.first == "rap" ||
290 item.first == "xi" ||
291 item.first == "pTIn" ||
292 item.first == "pTOut") scale(item.second, 1./sumOfWeights());
293 else normalize(item.second);
294 }
295 }
296
297 /// @}
298
299
300 /// @name Histograms
301 /// @{
302 map<string, Histo1DPtr> _h;
303 map<string, BinnedHistoPtr<string>> _d;
304 map<string, YODA::Axis<double>> _axis;
305 map<string, vector<string>> _edges;
306 bool isDisc;
307
308 /// @}
309 };
310
311
312 RIVET_DECLARE_PLUGIN(DELPHI_2003_I620250);
313
314
315}
|