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      
 41      if      (isCompatibleWithSqrtS(45)) offset = 1;
 42      else if (isCompatibleWithSqrtS(66)) offset = 2;
 43      else if (isCompatibleWithSqrtS(76)) offset = 3;
 44      else if (isCompatibleWithSqrtS(183)) {
 45	offset2= 0;			   
 46	offset = 1;			   
 47      }					   
 48      else if (isCompatibleWithSqrtS(189)) {
 49	offset2= 0;			   
 50	offset = 2;			   
 51      }					   
 52      else if (isCompatibleWithSqrtS(192)) {
 53	offset2= 0;			   
 54	offset = 3;			   
 55      }					   
 56      else if (isCompatibleWithSqrtS(196)) {
 57	offset2= 0;			   
 58	offset = 4;			   
 59      }					   
 60      else if (isCompatibleWithSqrtS(200)) {
 61	offset2= 1;			   
 62	offset = 1;			   
 63      }					   
 64      else if (isCompatibleWithSqrtS(202)) {
 65	offset2= 1;			   
 66	offset = 2;			   
 67      }					   
 68      else if (isCompatibleWithSqrtS(205)) {
 69	offset2= 1;			   
 70	offset = 3;			   
 71      }					   
 72      else if (isCompatibleWithSqrtS(207)) {
 73	offset2= 1;
 74	offset = 4;
 75      }
 76      else    MSG_ERROR("Beam energy not supported!");
 77      // Book the histograms
 78      if(offset2 < 0) {
 79	book(_h_thrust, 1, 1, offset);
 80	book(_h_major, 2, 1, offset);
 81	book(_h_minor, 3, 1, offset);
 82	book(_h_sphericity, 4, 1, offset);
 83	book(_h_planarity, 5, 1, offset);
 84	book(_h_oblateness, 6, 1, offset);
 85	book(_h_heavy_jet_mass, 7, 1, offset);
 86	book(_h_light_jet_mass, 9, 1, offset);
 87	book(_h_diff_jet_mass, 10, 1, offset);
 88	book(_h_total_jet_mass, 11, 1, offset);
 89	book(_h_heavy_jet_mass_E,  8, 1, offset);
 90	book(_h_total_jet_mass_E, 12, 1, offset);
 91	book(_h_wide_broading, 13, 1, offset);
 92	book(_h_narrow_broading, 14, 1, offset);
 93	book(_h_total_broading, 15, 1, offset);
 94	book(_h_diff_broading, 16, 1, offset);
 95	book(_h_CParam, 17, 1, offset);
 96      }
 97      else {
 98	book(_h_rap, 30+offset2, 1, offset);
 99	book(_h_xi, 32+offset2, 1, offset);
100	book(_h_pTIn, 34+offset2, 1, offset);
101	book(_h_pTOut, 36+offset2, 1, offset);
102	book(_h_thrust, 38+offset2, 1, offset);
103	book(_h_major, 40+offset2, 1, offset);
104	book(_h_minor, 42+offset2, 1, offset);
105	book(_h_oblateness, 44+offset2, 1, offset);
106	book(_h_wide_broading, 46+offset2, 1, offset);
107	book(_h_total_broading, 48+offset2, 1, offset);
108	book(_h_diff_broading, 50+offset2, 1, offset);
109	book(_h_CParam, 52+offset2, 1, offset);
110	book(_h_DParam, 54+offset2, 1, offset);
111	book(_h_heavy_jet_mass, 56+offset2, 1, offset);
112	book(_h_heavy_jet_mass_P, 58+offset2, 1, offset);
113	book(_h_heavy_jet_mass_E, 60+offset2, 1, offset);
114	book(_h_light_jet_mass, 62+offset2, 1, offset);
115	book(_h_diff_jet_mass, 64+offset2, 1, offset);
116	book(_h_sphericity, 66+offset2, 1, offset);
117	book(_h_planarity, 68+offset2, 1, offset);
118	book(_h_aplanarity, 70+offset2, 1, offset);
119      }
120    }
121
122
123    /// Perform the per-event analysis
124    void analyze(const Event& event) {
125
126      // Get beams and average beam momentum
127      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
128      const double meanBeamMom = ( beams.first.p3().mod() +
129                                   beams.second.p3().mod() ) / 2.0;
130      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
131
132      const Thrust& thrust = apply<Thrust>(event, "Thrust");
133      // thrust related observables
134      _h_thrust    ->fill(1.-thrust.thrust()  );
135      _h_major     ->fill(thrust.thrustMajor());
136      _h_minor     ->fill(thrust.thrustMinor());
137      _h_oblateness->fill(thrust.oblateness() );
138
139      // sphericity related
140      const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
141      _h_sphericity->fill(sphericity.sphericity());
142      _h_planarity ->fill(sphericity.planarity() );
143      if(_h_aplanarity) _h_aplanarity->fill(sphericity.aplanarity());
144      // hemisphere related
145      const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
146      // standard jet masses
147      _h_heavy_jet_mass->fill(hemi.scaledM2high());
148      _h_light_jet_mass->fill(hemi.scaledM2low() );
149      _h_diff_jet_mass ->fill(hemi.scaledM2diff());
150      if(_h_total_jet_mass) _h_total_jet_mass->fill(hemi.scaledM2low()+hemi.scaledM2high());
151      // jet broadening
152      _h_wide_broading  ->fill(hemi.Bmax() );
153      if(_h_narrow_broading) _h_narrow_broading->fill(hemi.Bmin() );
154      _h_total_broading ->fill(hemi.Bsum() );
155      _h_diff_broading  ->fill(hemi.Bdiff());
156      // E and p scheme jet masses
157      Vector3 axis = thrust.thrustAxis();
158      FourMomentum p4WithE, p4AgainstE;
159      FourMomentum p4WithP, p4AgainstP;
160      double Evis(0);
161      for (const Particle& p : apply<FinalState>(event, "FS").particles()) {
162	Vector3 p3 = p.momentum().vector3().unitVec();
163	const double   E = p.momentum().E();
164	Evis += E;
165	p3 = E*p3;
166	const double p3Para = dot(p3, axis);
167	FourMomentum p4E(E,p3.x(),p3.y(),p3.z());
168	FourMomentum p4P(p.p3().mod(),p.p3().x(),p.p3().y(),p.p3().z());
169	if (p3Para > 0)      {
170	  p4WithE    += p4E;
171	  p4WithP    += p4P;
172	}
173	else if (p3Para < 0) {
174	  p4AgainstE += p4E;
175	  p4AgainstP += p4P;
176	}
177	else {
178	  MSG_WARNING("Particle split between hemispheres");
179	  p4WithE    += 0.5 * p4E;
180	  p4AgainstE += 0.5 * p4E;
181	  p4WithP    += 0.5 * p4P;
182	  p4AgainstP += 0.5 * p4P;
183	}
184      }
185      // E scheme
186      const double mass2With_E    = p4WithE.mass2()/sqr(Evis);
187      const double mass2Against_E = p4AgainstE.mass2()/sqr(Evis);
188      // fill the histograms
189      _h_heavy_jet_mass_E->fill(max(mass2With_E,mass2Against_E));
190      if(_h_total_jet_mass_E) _h_total_jet_mass_E->fill(mass2With_E+mass2Against_E);
191      // pscheme
192      const double mass2With_P    = p4WithP.mass2()/sqr(Evis);
193      const double mass2Against_P = p4AgainstP.mass2()/sqr(Evis);
194      // fill the histograms
195      if(_h_heavy_jet_mass_P) _h_heavy_jet_mass_P->fill(max(mass2With_P,mass2Against_P));
196      
197      MSG_DEBUG("Calculating Parisi params");
198      const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
199      _h_CParam->fill(parisi.C());
200      if(_h_DParam) _h_DParam->fill(parisi.D());
201
202      // single particle distributions
203      const FinalState& fs = apply<FinalState>(event, "FS");
204      if(_h_xi) {
205	for (const Particle& p : fs.particles()) {
206	  if( ! PID::isCharged(p.pid())) continue;
207	  // Get momentum and energy of each particle.
208	  const Vector3 mom3 = p.p3();
209	  const double energy = p.E();
210	  
211	  // Scaled momenta.
212	  const double mom = mom3.mod();
213	  const double scaledMom = mom/meanBeamMom;
214	  const double logInvScaledMom = -std::log(scaledMom);
215	  _h_xi->fill(logInvScaledMom);
216	  
217	  // Get momenta components w.r.t. thrust and sphericity.
218	  const double momT = dot(thrust.thrustAxis(), mom3);
219	  const double pTinT = dot(mom3, thrust.thrustMajorAxis());
220	  const double pToutT = dot(mom3, thrust.thrustMinorAxis());
221	  _h_pTIn ->fill(fabs(pTinT/GeV));
222	  _h_pTOut->fill(fabs(pToutT/GeV));
223	  
224	  // Calculate rapidities w.r.t. thrust and sphericity.
225	  const double rapidityT = 0.5 * std::log((energy + momT) / (energy - momT));
226	  _h_rap->fill(fabs(rapidityT));
227	  MSG_TRACE(fabs(rapidityT) << " " << scaledMom/GeV);
228	}
229      }
230    }
231
232
233    /// Normalise histograms etc., after the run
234    void finalize() {
235
236      normalize(_h_thrust          );
237      normalize(_h_major           );
238      normalize(_h_minor           );
239      normalize(_h_sphericity      );
240      normalize(_h_planarity       );
241      if(_h_aplanarity) normalize(_h_aplanarity       );
242      normalize(_h_oblateness      );
243      normalize(_h_heavy_jet_mass  );
244      normalize(_h_light_jet_mass  );
245      normalize(_h_diff_jet_mass   );
246      if(_h_total_jet_mass) normalize(_h_total_jet_mass  );
247      normalize(_h_heavy_jet_mass_E);
248      if(_h_total_jet_mass_E) normalize(_h_total_jet_mass_E);
249      if(_h_heavy_jet_mass_P) normalize(_h_heavy_jet_mass_P);
250      normalize(_h_wide_broading   );
251      if(_h_narrow_broading) normalize(_h_narrow_broading );
252      normalize(_h_total_broading  );
253      normalize(_h_diff_broading   );
254      normalize(_h_CParam   );
255      if(_h_DParam) normalize(_h_DParam   );
256      if(_h_xi) {
257	scale(_h_xi   ,1./sumOfWeights());
258	scale(_h_pTIn ,1./sumOfWeights());
259	scale(_h_pTOut,1./sumOfWeights());
260	scale(_h_rap  ,1./sumOfWeights());
261      }
262    }
263
264    //@}
265
266
267    /// @name Histograms
268    //@{
269    Histo1DPtr _h_thrust,_h_major,_h_minor;
270    Histo1DPtr _h_sphericity,_h_planarity,_h_aplanarity,_h_oblateness;
271    Histo1DPtr _h_heavy_jet_mass,_h_light_jet_mass,_h_diff_jet_mass,_h_total_jet_mass;
272    Histo1DPtr _h_heavy_jet_mass_E,_h_total_jet_mass_E;
273    Histo1DPtr _h_heavy_jet_mass_P;
274    Histo1DPtr _h_wide_broading,_h_narrow_broading,_h_total_broading,_h_diff_broading;
275    Histo1DPtr _h_CParam,_h_DParam;
276    Histo1DPtr _h_xi, _h_pTIn, _h_pTOut,_h_rap;
277    //@}
278  };
279
280
281  // The hook for the plugin system
282  RIVET_DECLARE_PLUGIN(DELPHI_2003_I620250);
283
284
285}