rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

DELPHI_2003_I620250

Measurements of event shapes by DELPHI, above and below $m_Z$
Experiment: DELPHI (LEP)
Inspire ID: 620250
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Eur.Phys.J. C29 (2003) 285-312
Beams: e+ e-
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:
  • Hadronic Z decay events generated below the Z pole. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples.

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}