rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ALEPH_2004_S5765862

Jet rates and event shapes at LEP I and II
Experiment: ALEPH (LEP Run 1 and 2)
Inspire ID: 636645
Status: VALIDATED
Authors:
  • Frank Siegert
References: Beams: e+ e-
Beam energies: (45.6, 45.6); (66.5, 66.5); (80.5, 80.5); (86.0, 86.0); (91.5, 91.5); (94.5, 94.5); (98.5, 98.5); (100.0, 100.0); (103.0, 103.0) GeV
Run details:
  • $e^+ e^- \to$ jet jet (+ jets). Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples.

Jet rates, event-shape variables and inclusive charged particle spectra are measured in $e^+ e^-$ collisions at CMS energies between 91 and 209 GeV. The previously published data at 91.2 GeV and 133 GeV have been re-processed and the higher energy data are presented here for the first time. Note that the data have been corrected to include neutrinos. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples.

Source code: ALEPH_2004_S5765862.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/Thrust.hh"
  7#include "Rivet/Projections/Sphericity.hh"
  8#include "Rivet/Projections/ParisiTensor.hh"
  9#include "Rivet/Projections/Hemispheres.hh"
 10#include "Rivet/Projections/Beam.hh"
 11
 12namespace Rivet {
 13
 14
 15  /// ALEPH jet rates and event shapes at LEP 1 and 2
 16  class ALEPH_2004_S5765862 : public Analysis {
 17  public:
 18
 19    RIVET_DEFAULT_ANALYSIS_CTOR(ALEPH_2004_S5765862);
 20
 21
 22    void init() {
 23      _initialisedJets    = true;
 24      _initialisedSpectra = true;
 25      // TODO: According to the paper they seem to discard neutral particles
 26      //       between 1 and 2 GeV. That correction is included in the systematic
 27      //       uncertainties and overly complicated to program, so we ignore it.
 28      const FinalState fs;
 29      declare(fs, "FS");
 30      FastJets durhamjets(fs, FastJets::DURHAM, 0.7, JetAlg::Muons::ALL, JetAlg::Invisibles::ALL);
 31      declare(durhamjets, "DurhamJets");
 32
 33      const Thrust thrust(fs);
 34      declare(thrust, "Thrust");
 35      declare(Sphericity(fs), "Sphericity");
 36      declare(ParisiTensor(fs), "Parisi");
 37      declare(Hemispheres(thrust), "Hemispheres");
 38
 39      const ChargedFinalState cfs;
 40      declare(Beam(), "Beams");
 41      declare(cfs, "CFS");
 42
 43      // Histos
 44      // offset for the event shapes and jets
 45      int offset = 0;
 46      switch (int(sqrtS()/GeV + 0.5)) {
 47        case 91: offset = 0; break;
 48        case 133: offset = 1; break;
 49        case 161: offset = 2; break;
 50        case 172: offset = 3; break;
 51        case 183: offset = 4; break;
 52        case 189: offset = 5; break;
 53        case 200: offset = 6; break;
 54        case 206: offset = 7; break;
 55        default:
 56          _initialisedJets = false;
 57      }
 58      // event shapes
 59      if(_initialisedJets) {
 60        book(_h_thrust ,offset+54, 1, 1);
 61        book(_h_heavyjetmass ,offset+62, 1, 1);
 62        book(_h_totaljetbroadening ,offset+70, 1, 1);
 63        book(_h_widejetbroadening ,offset+78, 1, 1);
 64        book(_h_cparameter ,offset+86, 1, 1);
 65        book(_h_thrustmajor ,offset+94, 1, 1);
 66        book(_h_thrustminor ,offset+102, 1, 1);
 67        book(_h_jetmassdifference ,offset+110, 1, 1);
 68        book(_h_aplanarity ,offset+118, 1, 1);
 69        if ( offset != 0 )
 70          book(_h_planarity, offset+125, 1, 1);
 71        book(_h_oblateness ,offset+133, 1, 1);
 72        book(_h_sphericity ,offset+141, 1, 1);
 73
 74        // Durham n->m jet resolutions
 75        book(_h_y_Durham[0] ,offset+149, 1, 1);   // y12 d149 ... d156
 76        book(_h_y_Durham[1] ,offset+157, 1, 1);   // y23 d157 ... d164
 77        if (offset < 6) { // there is no y34, y45 and y56 for 200 gev
 78          book(_h_y_Durham[2] ,offset+165, 1, 1); // y34 d165 ... d172, but not 171
 79          book(_h_y_Durham[3] ,offset+173, 1, 1); // y45 d173 ... d179
 80          book(_h_y_Durham[4] ,offset+180, 1, 1); // y56 d180 ... d186
 81        }
 82        else if (offset == 6) {
 83          _h_y_Durham[2] = Histo1DPtr();
 84          _h_y_Durham[3] = Histo1DPtr();
 85          _h_y_Durham[4] = Histo1DPtr();
 86        }
 87        else if (offset == 7) {
 88          book(_h_y_Durham[2] ,172, 1, 1);
 89          book(_h_y_Durham[3] ,179, 1, 1);
 90          book(_h_y_Durham[4] ,186, 1, 1);
 91        }
 92
 93        // Durham n-jet fractions
 94        book(_h_R_Durham[0] ,offset+187, 1, 1); // R1 d187 ... d194
 95        book(_h_R_Durham[1] ,offset+195, 1, 1); // R2 d195 ... d202
 96        book(_h_R_Durham[2] ,offset+203, 1, 1); // R3 d203 ... d210
 97        book(_h_R_Durham[3] ,offset+211, 1, 1); // R4 d211 ... d218
 98        book(_h_R_Durham[4] ,offset+219, 1, 1); // R5 d219 ... d226
 99        book(_h_R_Durham[5] ,offset+227, 1, 1); // R>=6 d227 ... d234
100      }
101      // offset for the charged particle distributions
102      offset = 0;
103      switch (int(sqrtS() + 0.5)) {
104        case 133: offset = 0; break;
105        case 161: offset = 1; break;
106        case 172: offset = 2; break;
107        case 183: offset = 3; break;
108        case 189: offset = 4; break;
109        case 196: offset = 5; break;
110        case 200: offset = 6; break;
111        case 206: offset = 7; break;
112        default:
113          _initialisedSpectra = false;
114      }
115      if (_initialisedSpectra) {
116        book(_h_xp , 2+offset, 1, 1);
117        book(_h_xi ,11+offset, 1, 1);
118        book(_h_xe ,19+offset, 1, 1);
119        book(_h_pTin  ,27+offset, 1, 1);
120        if (offset == 7)
121          book(_h_pTout, 35, 1, 1);
122        book(_h_rapidityT ,36+offset, 1, 1);
123        book(_h_rapidityS ,44+offset, 1, 1);
124      }
125      book(_weightedTotalChargedPartNum, "_weightedTotalChargedPartNum");
126      if (!_initialisedSpectra && !_initialisedJets) {
127        MSG_WARNING("CoM energy of events sqrt(s) = " << sqrtS()/GeV
128                    << " doesn't match any available analysis energy .");
129      }
130
131      book(mult, 1, 1, 1);
132    }
133
134
135    void analyze(const Event& e) {
136
137      const Thrust& thrust = apply<Thrust>(e, "Thrust");
138      const Sphericity& sphericity = apply<Sphericity>(e, "Sphericity");
139
140      if(_initialisedJets) {
141        bool LEP1 = isCompatibleWithSqrtS(91.2*GeV,0.01);
142        // event shapes
143        double thr = LEP1 ? thrust.thrust() : 1.0 - thrust.thrust();
144        _h_thrust->fill(thr);
145        _h_thrustmajor->fill(thrust.thrustMajor());
146        if(LEP1)
147          _h_thrustminor->fill(log(thrust.thrustMinor()));
148        else
149          _h_thrustminor->fill(thrust.thrustMinor());
150        _h_oblateness->fill(thrust.oblateness());
151
152        const Hemispheres& hemi = apply<Hemispheres>(e, "Hemispheres");
153        _h_heavyjetmass->fill(hemi.scaledM2high());
154        _h_jetmassdifference->fill(hemi.scaledM2diff());
155        _h_totaljetbroadening->fill(hemi.Bsum());
156        _h_widejetbroadening->fill(hemi.Bmax());
157
158        const ParisiTensor& parisi = apply<ParisiTensor>(e, "Parisi");
159        _h_cparameter->fill(parisi.C());
160
161        _h_aplanarity->fill(sphericity.aplanarity());
162        if(_h_planarity)
163          _h_planarity->fill(sphericity.planarity());
164        _h_sphericity->fill(sphericity.sphericity());
165
166        // Jet rates
167        const FastJets& durjet = apply<FastJets>(e, "DurhamJets");
168        double log10e = log10(exp(1.));
169        if (durjet.clusterSeq()) {
170          double logynm1=0.;
171          double logyn;
172          for (size_t i=0; i<5; ++i) {
173            double yn = durjet.clusterSeq()->exclusive_ymerge_max(i+1);
174            if (yn<=0.0) continue;
175            logyn = -log(yn);
176            if (_h_y_Durham[i]) {
177              _h_y_Durham[i]->fill(logyn);
178            }
179            if(!LEP1) logyn *= log10e;
180            for (size_t j = 0; j < _h_R_Durham[i]->numBins(); ++j) {
181              double val   = _h_R_Durham[i]->bin(j).xMin();
182              double width = _h_R_Durham[i]->bin(j).xWidth();
183              if(-val<=logynm1) break;
184              if(-val<logyn) {
185                _h_R_Durham[i]->fill(val+0.5*width, width);
186              }
187            }
188            logynm1 = logyn;
189          }
190          for (size_t j = 0; j < _h_R_Durham[5]->numBins(); ++j) {
191            double val   = _h_R_Durham[5]->bin(j).xMin();
192            double width = _h_R_Durham[5]->bin(j).xWidth();
193            if(-val<=logynm1) break;
194            _h_R_Durham[5]->fill(val+0.5*width, width);
195          }
196        }
197        if( !_initialisedSpectra) {
198          const ChargedFinalState& cfs = apply<ChargedFinalState>(e, "CFS");
199          const size_t numParticles = cfs.particles().size();
200          _weightedTotalChargedPartNum->fill(numParticles);
201        }
202      }
203
204      // charged particle distributions
205      if(_initialisedSpectra) {
206        const ChargedFinalState& cfs = apply<ChargedFinalState>(e, "CFS");
207        const size_t numParticles = cfs.particles().size();
208        _weightedTotalChargedPartNum->fill(numParticles);
209        const ParticlePair& beams = apply<Beam>(e, "Beams").beams();
210        const double meanBeamMom = ( beams.first.p3().mod() +
211                                     beams.second.p3().mod() ) / 2.0;
212        for (const Particle& p : cfs.particles()) {
213          const double xp = p.p3().mod()/meanBeamMom;
214          _h_xp->fill(xp   );
215          const double logxp = -std::log(xp);
216          _h_xi->fill(logxp);
217          const double xe = p.E()/meanBeamMom;
218          _h_xe->fill(xe   );
219          const double pTinT  = dot(p.p3(), thrust.thrustMajorAxis());
220          const double pToutT = dot(p.p3(), thrust.thrustMinorAxis());
221          _h_pTin->fill(fabs(pTinT/GeV));
222          if(_h_pTout) _h_pTout->fill(fabs(pToutT/GeV));
223          const double momT = dot(thrust.thrustAxis()        ,p.p3());
224          const double rapidityT = 0.5 * std::log((p.E() + momT) /
225                                                  (p.E() - momT));
226          _h_rapidityT->fill(fabs(rapidityT));
227          const double momS = dot(sphericity.sphericityAxis(),p.p3());
228          const double rapidityS = 0.5 * std::log((p.E() + momS) /
229                                                  (p.E() - momS));
230          _h_rapidityS->fill(fabs(rapidityS));
231        }
232      }
233    }
234
235    void finalize() {
236      if(!_initialisedJets && !_initialisedSpectra) return;
237
238      if (_initialisedJets) {
239        normalize(_h_thrust);
240        normalize(_h_heavyjetmass);
241        normalize(_h_totaljetbroadening);
242        normalize(_h_widejetbroadening);
243        normalize(_h_cparameter);
244        normalize(_h_thrustmajor);
245        normalize(_h_thrustminor);
246        normalize(_h_jetmassdifference);
247        normalize(_h_aplanarity);
248        if(_h_planarity) normalize(_h_planarity);
249        normalize(_h_oblateness);
250        normalize(_h_sphericity);
251
252        for (size_t n=0; n<6; ++n) {
253          scale(_h_R_Durham[n], 1./sumOfWeights());
254        }
255
256        for (size_t n = 0; n < 5; ++n) {
257          if (_h_y_Durham[n]) {
258            scale(_h_y_Durham[n], 1.0/sumOfWeights());
259          }
260        }
261      }
262
263      Histo1D temphisto(refData(1, 1, 1));
264      const double avgNumParts = dbl(*_weightedTotalChargedPartNum) / sumOfWeights();
265
266
267      for (size_t b = 0; b < temphisto.numBins(); b++) {
268        const double x  = temphisto.bin(b).xMid();
269        const double ex = temphisto.bin(b).xWidth()/2.;
270        if (inRange(sqrtS()/GeV, x-ex, x+ex)) {
271          mult->addPoint(x, avgNumParts, ex, 0.);
272        }
273      }
274
275      if (_initialisedSpectra) {
276        normalize(_h_xp, avgNumParts);
277        normalize(_h_xi, avgNumParts);
278        normalize(_h_xe, avgNumParts);
279        normalize(_h_pTin , avgNumParts);
280        if (_h_pTout) normalize(_h_pTout, avgNumParts);
281        normalize(_h_rapidityT, avgNumParts);
282        normalize(_h_rapidityS, avgNumParts);
283      }
284    }
285
286
287  private:
288
289    bool _initialisedJets = false;
290    bool _initialisedSpectra = false;
291
292    Scatter2DPtr mult;
293    Histo1DPtr _h_xp;
294    Histo1DPtr _h_xi;
295    Histo1DPtr _h_xe;
296    Histo1DPtr _h_pTin;
297    Histo1DPtr _h_pTout;
298    Histo1DPtr _h_rapidityT;
299    Histo1DPtr _h_rapidityS;
300    Histo1DPtr _h_thrust;
301    Histo1DPtr _h_heavyjetmass;
302    Histo1DPtr _h_totaljetbroadening;
303    Histo1DPtr _h_widejetbroadening;
304    Histo1DPtr _h_cparameter;
305    Histo1DPtr _h_thrustmajor;
306    Histo1DPtr _h_thrustminor;
307    Histo1DPtr _h_jetmassdifference;
308    Histo1DPtr _h_aplanarity;
309    Histo1DPtr _h_planarity;
310    Histo1DPtr _h_oblateness;
311    Histo1DPtr _h_sphericity;
312
313    Histo1DPtr _h_R_Durham[6];
314    Histo1DPtr _h_y_Durham[5];
315
316    CounterPtr _weightedTotalChargedPartNum;
317
318  };
319
320
321
322  RIVET_DECLARE_ALIASED_PLUGIN(ALEPH_2004_S5765862, ALEPH_2004_I636645);
323
324}