rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

L3_2004_I652683

Jet rates and event shapes at LEP I+II
Experiment: L3 (LEP I+II)
Inspire ID: 652683
Status: VALIDATED
Authors:
  • Adil Jueid
References:
  • Phys.Rept. 399 (2004) 71-174
  • 10.1016/j.physrep.2004.07.002
  • arXiv: hep-ex/0406049
Beams: e+ e-
Beam energies: (20.7, 20.7); (27.6, 27.6); (32.7, 32.7); (37.9, 37.9); (41.1, 41.1); (42.5, 42.5); (45.6, 45.6); (65.0, 65.0); (68.0, 68.0); (80.7, 80.7); (86.2, 86.2); (91.4, 91.4); (94.3, 94.3); (97.2, 97.2); (100.1, 100.1); (103.1, 103.1) GeV
Run details:
  • $e^+ e^- \to q \bar{q}$ at several energies. Beam energy must be specified as analysis option "ENERGY" when rivet-merge'ing samples.

Event shapes, charged particle multiplicity and scaled momentum distributions, flavour separated at the Z-boson peak. The analysis also includes non-flavour separated event shapes, charged particle multiplicity and scaled momentum distributions at events both above and below the Z pole. There are also jet distributions in the paper which are not currently implemented. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples.

Source code: L3_2004_I652683.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/FinalState.hh"
  6#include "Rivet/Projections/ChargedFinalState.hh"
  7#include "Rivet/Projections/Thrust.hh"
  8#include "Rivet/Projections/ParisiTensor.hh"
  9#include "Rivet/Projections/Hemispheres.hh"
 10
 11#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
 12#include "Rivet/Projections/InitialQuarks.hh"
 13#include "fastjet/EECambridgePlugin.hh"
 14
 15namespace Rivet {
 16
 17
 18  /// Jet rates and event shapes at LEP I+II
 19  class L3_2004_I652683 : public Analysis {
 20  public:
 21
 22    /// Constructor
 23    RIVET_DEFAULT_ANALYSIS_CTOR(L3_2004_I652683);
 24
 25    /// Book histograms and initialise projections before the run
 26    void init() {
 27
 28      // Projections to use
 29      const FinalState FS;
 30      declare(FS, "FS");
 31      declare(Beam(), "beams");
 32      const ChargedFinalState CFS;
 33      declare(CFS, "CFS");
 34      const Thrust thrust(FS);
 35      declare(thrust, "thrust");
 36      declare(ParisiTensor(FS), "Parisi");
 37      declare(Hemispheres(thrust), "Hemispheres");
 38      declare(InitialQuarks(), "initialquarks");
 39      FastJets jadeJets = FastJets(FS, JetAlg::JADE, 0.7, JetMuons::ALL, JetInvisibles::DECAY);
 40      FastJets durhamJets = FastJets(FS, JetAlg::DURHAM, 0.7, JetMuons::ALL, JetInvisibles::DECAY);
 41      FastJets cambridgeJets = FastJets(FS, JetAlg::CAM, 0.7, JetMuons::ALL, JetInvisibles::DECAY);
 42      declare(jadeJets, "JadeJets");
 43      declare(durhamJets, "DurhamJets");
 44
 45      // Book the histograms
 46      if (isCompatibleWithSqrtS(91.2*GeV)) {
 47        // z pole
 48        book(_h_Thrust_udsc             , 47, 1, 1);
 49        book(_h_Thrust_bottom           , 47, 1, 2);
 50        book(_h_heavyJetmass_udsc       , 48, 1, 1);
 51        book(_h_heavyJetmass_bottom     , 48, 1, 2);
 52        book(_h_totalJetbroad_udsc      , 49, 1, 1);
 53        book(_h_totalJetbroad_bottom    , 49, 1, 2);
 54        book(_h_wideJetbroad_udsc       , 50, 1, 1);
 55        book(_h_wideJetbroad_bottom     , 50, 1, 2);
 56        book(_h_Cparameter_udsc         , 51, 1, 1);
 57        book(_h_Cparameter_bottom       , 51, 1, 2);
 58        book(_h_Dparameter_udsc         , 52, 1, 1);
 59        book(_h_Dparameter_bottom       , 52, 1, 2);
 60        book(_h_Ncharged                , "/TMP/NCHARGED"     , 28, 1, 57);
 61        book(_h_Ncharged_udsc           , "/TMP/NCHARGED_UDSC", 28, 1, 57);
 62        book(_h_Ncharged_bottom         , "/TMP/NCHARGED_B"   , 27, 3, 57);
 63        book(_h_scaledMomentum          , 65, 1, 1);
 64        book(_h_scaledMomentum_udsc     , 65, 1, 2);
 65        book(_h_scaledMomentum_bottom   , 65, 1, 3);
 66      }
 67      else if(sqrtS()/GeV < 90) {
 68        int i1(-1),i2(-1);
 69        if(isCompatibleWithSqrtS(41.4*GeV)) {
 70          i1=0; i2=1;
 71        }
 72        else if(isCompatibleWithSqrtS(55.3*GeV)) {
 73          i1=0; i2=2;
 74        }
 75        else if(isCompatibleWithSqrtS(65.4*GeV)) {
 76          i1=0; i2=3;
 77        }
 78        else if(isCompatibleWithSqrtS(75.7*GeV)) {
 79          i1=1; i2=1;
 80        }
 81        else if(isCompatibleWithSqrtS(82.3*GeV)) {
 82          i1=1; i2=2;
 83        }
 84        else if(isCompatibleWithSqrtS(85.1*GeV)) {
 85          i1=1; i2=3;
 86        }
 87        else {
 88          MSG_ERROR("Beam energy not supported!");
 89        }
 90
 91        book(_h_thrust , 21+i1,1,i2);
 92        book(_h_rho    , 26+i1,1,i2);
 93        book(_h_B_T    , 31+i1,1,i2);
 94        book(_h_B_W    , 36+i1,1,i2);
 95      }
 96      else if(sqrtS()/GeV > 120) {
 97        int i1(-1),i2(-1);
 98        if(isCompatibleWithSqrtS(130.1*GeV)) {
 99          i1=0; i2=1;
100        }
101        else if(isCompatibleWithSqrtS(136.1*GeV)) {
102          i1=0; i2=2;
103        }
104        else if(isCompatibleWithSqrtS(161.3*GeV)) {
105          i1=0; i2=3;
106        }
107        else if(isCompatibleWithSqrtS(172.3*GeV)) {
108          i1=1; i2=1;
109        }
110        else if(isCompatibleWithSqrtS(182.8*GeV)) {
111          i1=1; i2=2;
112        }
113        else if(isCompatibleWithSqrtS(188.6*GeV)) {
114          i1=1; i2=3;
115        }
116        else if(isCompatibleWithSqrtS(194.4*GeV)) {
117          i1=2; i2=1;
118        }
119        else if(isCompatibleWithSqrtS(200.2*GeV)) {
120          i1=2; i2=2;
121        }
122        else if(isCompatibleWithSqrtS(206.2*GeV)) {
123          i1=2; i2=3;
124        }
125        else {
126          MSG_ERROR("Beam energy not supported!");
127        }
128
129        book(_h_thrust , 23+i1,1,i2);
130        book(_h_rho    , 28+i1,1,i2);
131        book(_h_B_T    , 33+i1,1,i2);
132        book(_h_B_W    , 38+i1,1,i2);
133        book(_h_C      , 41+i1,1,i2);
134        book(_h_D      , 44+i1,1,i2);
135        book(_h_N      , "/TMP/NCHARGED", 22, 9, 53);
136        book(_h_xi     , 66+i1,1,i2);
137        // and the jets
138        int i3 = 3*i1+i2;
139        book(_s["y_2_JADE"],     i3,1,1);
140        book(_s["y_3_JADE"],     i3,1,2);
141        book(_s["y_4_JADE"],     i3,1,3);
142        book(_s["y_5_JADE"],     i3,1,4);
143        book(_s["y_2_Durham"], 9+i3,1,1);
144        book(_s["y_3_Durham"], 9+i3,1,2);
145        book(_s["y_4_Durham"], 9+i3,1,3);
146        book(_s["y_5_Durham"], 9+i3,1,4);
147        if (i3==8 || i3==9) {
148          book(_s["y_2_Cambridge"], 11+i3,1,1);
149          book(_s["y_3_Cambridge"], 11+i3,1,2);
150          book(_s["y_4_Cambridge"], 11+i3,1,3);
151          book(_s["y_5_Cambridge"], 11+i3,1,4);
152        }
153      }
154
155      book(_sumW_udsc, "_sumW_udsc");
156      book(_sumW_b, "_sumW_b");
157      book(_sumW_ch, "_sumW_ch");
158      book(_sumW_ch_udsc, "_sumW_ch_udsc");
159      book(_sumW_ch_b, "_sumW_ch_b");
160
161    }
162
163
164    /// Perform the per-event analysis
165    void analyze(const Event& event) {
166      // Get beam average momentum
167      const ParticlePair& beams = apply<Beam>(event, "beams").beams();
168      const double beamMomentum = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0;
169
170      // InitialQuarks projection to have udsc events separated from b events
171      /// @todo Yuck!!! Eliminate when possible...
172      int iflav = 0;
173      // only need the flavour at Z pole
174      if (_h_Thrust_udsc) {
175        int flavour = 0;
176        const InitialQuarks& iqf = apply<InitialQuarks>(event, "initialquarks");
177        Particles quarks;
178        if ( iqf.particles().size() == 2 ) {
179          flavour = iqf.particles().front().abspid();
180          quarks  = iqf.particles();
181        } else {
182          map<int, Particle> quarkmap;
183          for (const Particle& p : iqf.particles()) {
184            if (quarkmap.find(p.pid()) == quarkmap.end()) quarkmap[p.pid()] = p;
185            else if (quarkmap[p.pid()].E() < p.E()) quarkmap[p.pid()] = p;
186          }
187          double max_energy = 0.;
188          for (int i = 1; i <= 5; ++i) {
189            double energy = 0.;
190            if (quarkmap.find(i) != quarkmap.end())
191              energy += quarkmap[ i].E();
192            if (quarkmap.find(-i) != quarkmap.end())
193              energy += quarkmap[-i].E();
194            if (energy > max_energy)
195              flavour = i;
196          }
197          if (quarkmap.find(flavour) != quarkmap.end())
198            quarks.push_back(quarkmap[flavour]);
199          if (quarkmap.find(-flavour) != quarkmap.end())
200            quarks.push_back(quarkmap[-flavour]);
201        }
202        // Flavour label
203        /// @todo Change to a bool?
204        iflav = (flavour == PID::DQUARK || flavour == PID::UQUARK || flavour == PID::SQUARK || flavour == PID::CQUARK) ? 1 : (flavour == PID::BQUARK) ? 5 : 0;
205      }
206      // Update weight sums
207      if (iflav == 1) {
208        _sumW_udsc->fill();
209      } else if (iflav == 5) {
210        _sumW_b->fill();
211      }
212      _sumW_ch->fill();
213
214      // Charged multiplicity
215      const FinalState& cfs = apply<FinalState>(event, "CFS");
216      if(_h_Ncharged) _h_Ncharged->fill(cfs.size());
217      if (iflav == 1) {
218        _sumW_ch_udsc->fill();
219        _h_Ncharged_udsc->fill(cfs.size());
220      } else if (iflav == 5) {
221        _sumW_ch_b->fill();
222        _h_Ncharged_bottom->fill(cfs.size());
223      }
224      else if(_h_N) {
225        _h_N->fill(cfs.size());
226      }
227
228      // Scaled momentum
229      const Particles& chparticles = cfs.particlesByPt();
230      for (const Particle& p : chparticles) {
231        const Vector3 momentum3 = p.p3();
232        const double mom = momentum3.mod();
233        const double scaledMom = mom/beamMomentum;
234        const double logScaledMom = std::log(scaledMom);
235        if(_h_scaledMomentum) _h_scaledMomentum->fill(-logScaledMom);
236        if (iflav == 1) {
237          _h_scaledMomentum_udsc->fill(-logScaledMom);
238        } else if (iflav == 5) {
239          _h_scaledMomentum_bottom->fill(-logScaledMom);
240        }
241        else if(_h_xi) {
242          _h_xi->fill(-logScaledMom);
243        }
244      }
245
246      // Thrust
247      const Thrust& thrust = apply<Thrust>(event, "thrust");
248      if (iflav == 1) {
249        _h_Thrust_udsc->fill(thrust.thrust());
250      } else if (iflav == 5) {
251        _h_Thrust_bottom->fill(thrust.thrust());
252      }
253      else if(_h_thrust) {
254        _h_thrust->fill(1.-thrust.thrust());
255      }
256
257      // C and D Parisi parameters
258      const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
259      if (iflav == 1) {
260        _h_Cparameter_udsc->fill(parisi.C());
261        _h_Dparameter_udsc->fill(parisi.D());
262      } else if (iflav == 5) {
263        _h_Cparameter_bottom->fill(parisi.C());
264        _h_Dparameter_bottom->fill(parisi.D());
265      }
266      else if(_h_C) {
267        _h_C->fill(parisi.C());
268        _h_D->fill(parisi.D());
269      }
270
271      // The hemisphere variables
272      const Hemispheres& hemisphere = apply<Hemispheres>(event, "Hemispheres");
273      if (iflav == 1) {
274        _h_heavyJetmass_udsc->fill(hemisphere.scaledM2high());
275        _h_totalJetbroad_udsc->fill(hemisphere.Bsum());
276        _h_wideJetbroad_udsc->fill(hemisphere.Bmax());
277      } else if (iflav == 5) {
278        _h_heavyJetmass_bottom->fill(hemisphere.scaledM2high());
279        _h_totalJetbroad_bottom->fill(hemisphere.Bsum());
280        _h_wideJetbroad_bottom->fill(hemisphere.Bmax());
281      }
282      else if (_h_rho) {
283        _h_rho->fill(hemisphere.scaledM2high());
284        _h_B_T->fill(hemisphere.Bsum());
285        _h_B_W->fill(hemisphere.Bmax());
286      }
287      // jade jet rates
288      if (_s.count("y_2_JADE")) {
289        const FastJets& jadejet = apply<FastJets>(event, "JadeJets");
290        if (jadejet.clusterSeq()) {
291          const double y_23 = jadejet.clusterSeq()->exclusive_ymerge_max(2);
292          const double y_34 = jadejet.clusterSeq()->exclusive_ymerge_max(3);
293          const double y_45 = jadejet.clusterSeq()->exclusive_ymerge_max(4);
294          const double y_56 = jadejet.clusterSeq()->exclusive_ymerge_max(5);
295          for (auto& b : _s["y_2_JADE"]->bins()) {
296            const double ycut = std::stod(b.xEdge());
297            if (y_23 < ycut)  _s["y_2_JADE"]->fill(b.xEdge());
298          }
299          for (auto& b : _s["y_3_JADE"]->bins()) {
300            const double ycut = std::stod(b.xEdge());
301            if (y_34 < ycut && y_23 > ycut) {
302              _s["y_3_JADE"]->fill(b.xEdge());
303            }
304          }
305          for (auto& b : _s["y_4_JADE"]->bins()) {
306            const double ycut = std::stod(b.xEdge());
307            if (y_45 < ycut && y_34 > ycut) {
308              _s["y_4_JADE"]->fill(b.xEdge());
309            }
310          }
311          for (auto& b : _s["y_5_JADE"]->bins()) {
312            const double ycut = std::stod(b.xEdge());
313            if (y_56 < ycut && y_45 > ycut) {
314              _s["y_5_JADE"]->fill(b.xEdge());
315            }
316          }
317        }
318      }
319      // Durham jet rates
320      if (_s.count("y_2_Durham")) {
321        const FastJets& durhamjet = apply<FastJets>(event, "DurhamJets");
322        if (durhamjet.clusterSeq()) {
323          const double y_23 = durhamjet.clusterSeq()->exclusive_ymerge_max(2);
324          const double y_34 = durhamjet.clusterSeq()->exclusive_ymerge_max(3);
325          const double y_45 = durhamjet.clusterSeq()->exclusive_ymerge_max(4);
326          const double y_56 = durhamjet.clusterSeq()->exclusive_ymerge_max(5);
327          for (auto& b : _s["y_2_Durham"]->bins()) {
328            const double ycut = std::stod(b.xEdge());
329            if (y_23 < ycut)  _s["y_2_Durham"]->fill(b.xEdge());
330          }
331          for (auto& b : _s["y_3_Durham"]->bins()) {
332            const double ycut = std::stod(b.xEdge());
333            if (y_34 < ycut && y_23 > ycut) {
334              _s["y_3_Durham"]->fill(b.xEdge());
335            }
336          }
337          for (auto& b : _s["y_4_Durham"]->bins()) {
338            const double ycut = std::stod(b.xEdge());
339            if (y_45 < ycut && y_34 > ycut) {
340              _s["y_4_Durham"]->fill(b.xEdge());
341            }
342          }
343          for (auto& b : _s["y_5_Durham"]->bins()) {
344            const double ycut = std::stod(b.xEdge());
345            if (y_56 < ycut && y_45 > ycut) {
346              _s["y_5_Durham"]->fill(b.xEdge());
347            }
348          }
349        }
350      }
351      // Cambridge
352      if (_s.count("y_2_Cambridge")) {
353        PseudoJets pjs;
354        const FinalState& fs = apply<FinalState>(event, "FS");
355        for (size_t i = 1; i < fs.particles().size()+1; ++i) {
356          fastjet::PseudoJet pj = fs.particles()[i];
357          pjs.push_back(pj);
358        }
359        for (size_t i = 1; i < _s["y_2_Cambridge"]->numBins()+1; ++i) {
360          const auto& b = _s["y_2_Cambridge"]->bin(i);
361          const double ycut = std::stod(b.xEdge());
362          fastjet::EECambridgePlugin plugin(ycut);
363          fastjet::JetDefinition jdef(&plugin);
364          fastjet::ClusterSequence cseq(pjs, jdef);
365          size_t njet = cseq.inclusive_jets().size();
366          if (njet==2) {
367            _s["y_2_Cambridge"]->fill(b.xEdge());
368          }
369          else if (njet==3) {
370            if (i < _s["y_3_Cambridge"]->numBins()+1) {
371              _s["y_3_Cambridge"]->fill(b.xEdge());
372            }
373          }
374          else if(njet==4) {
375            if (i < _s["y_4_Cambridge"]->numBins()+1) {
376              _s["y_4_Cambridge"]->fill(b.xEdge());
377            }
378          }
379          else if(njet==5) {
380            if (i < _s["y_5_Cambridge"]->numBins()+1) {
381              _s["y_5_Cambridge"]->fill(b.xEdge());
382            }
383          }
384        }
385      }
386    }
387
388    Estimate1DPtr convertHisto(unsigned int ix,unsigned int iy, unsigned int iz, Histo1DPtr histo) {
389      Estimate1DPtr mult;
390      book(mult, ix, iy, iz);
391      barchart(histo, mult);
392      return mult;
393    }
394
395    /// Normalise histograms etc., after the run
396    void finalize() {
397      // Z pole plots
398      if(_h_Thrust_udsc) {
399        scale(_h_Thrust_udsc,  1/_sumW_udsc->sumW());
400        scale(_h_heavyJetmass_udsc,  1/_sumW_udsc->sumW());
401        scale(_h_totalJetbroad_udsc, 1/_sumW_udsc->sumW());
402        scale(_h_wideJetbroad_udsc, 1/_sumW_udsc->sumW());
403        scale(_h_Cparameter_udsc, 1/_sumW_udsc->sumW());
404        scale(_h_Dparameter_udsc, 1/_sumW_udsc->sumW());
405        scale(_h_Thrust_bottom, 1./_sumW_b->sumW());
406        scale(_h_heavyJetmass_bottom, 1./_sumW_b->sumW());
407        scale(_h_totalJetbroad_bottom, 1./_sumW_b->sumW());
408        scale(_h_wideJetbroad_bottom, 1./_sumW_b->sumW());
409        scale(_h_Cparameter_bottom, 1./_sumW_b->sumW());
410        scale(_h_Dparameter_bottom, 1./_sumW_b->sumW());
411        scale(_h_Ncharged, 1./_sumW_ch->sumW());
412        scale(_h_Ncharged_udsc, 1./_sumW_ch_udsc->sumW());
413        scale(_h_Ncharged_bottom, 1./_sumW_ch_b->sumW());
414        convertHisto(59,1,1,_h_Ncharged       );
415        convertHisto(59,1,2,_h_Ncharged_udsc  );
416        convertHisto(59,1,3,_h_Ncharged_bottom);
417
418        scale(_h_scaledMomentum, 1/_sumW_ch->sumW());
419        scale(_h_scaledMomentum_udsc, 1/_sumW_ch_udsc->sumW());
420        scale(_h_scaledMomentum_bottom, 1/_sumW_ch_b->sumW());
421      }
422      else {
423        if(_h_thrust) normalize(_h_thrust);
424        if(_h_rho) normalize(_h_rho);
425        if(_h_B_T) normalize(_h_B_T);
426        if(_h_B_W) normalize(_h_B_W);
427        if(_h_C) normalize(_h_C);
428        if(_h_D) normalize(_h_D);
429        if(_h_N) normalize(_h_N);
430        if(_h_xi) scale(_h_xi,1./sumOfWeights());
431
432
433        Estimate1DPtr mult;
434        if(_h_N) {
435          if(isCompatibleWithSqrtS(130.1*GeV)) {
436            convertHisto(60, 1, 1, _h_N);
437          }
438          else if(isCompatibleWithSqrtS(136.1*GeV)) {
439            convertHisto(60, 1, 2, _h_N);
440          }
441          else if(isCompatibleWithSqrtS(161.3*GeV)) {
442            convertHisto(60, 1, 3, _h_N);
443          }
444          else if(isCompatibleWithSqrtS(172.3*GeV)) {
445            convertHisto(61, 1, 1, _h_N);
446          }
447          else if(isCompatibleWithSqrtS(182.8*GeV)) {
448            convertHisto(61, 1, 2, _h_N);
449          }
450          else if(isCompatibleWithSqrtS(188.6*GeV)) {
451            convertHisto(61, 1, 3, _h_N);
452          }
453          else if(isCompatibleWithSqrtS(194.4*GeV)) {
454            convertHisto(62, 1, 1, _h_N);
455          }
456          else if(isCompatibleWithSqrtS(200.2*GeV)) {
457            convertHisto(62, 1, 2, _h_N);
458          }
459          else if(isCompatibleWithSqrtS(206.2*GeV)) {
460            convertHisto(62, 1, 3, _h_N);
461          }
462        }
463        // // the jets
464        scale(_s, 1./sumOfWeights());
465      }
466    }
467
468
469    /// Weight counters
470    CounterPtr _sumW_udsc, _sumW_b, _sumW_ch, _sumW_ch_udsc, _sumW_ch_b;
471
472    /// @name Histograms
473    /// @{
474    // at the Z pole
475    Histo1DPtr _h_Thrust_udsc, _h_Thrust_bottom;
476    Histo1DPtr _h_heavyJetmass_udsc, _h_heavyJetmass_bottom;
477    Histo1DPtr _h_totalJetbroad_udsc, _h_totalJetbroad_bottom;
478    Histo1DPtr _h_wideJetbroad_udsc, _h_wideJetbroad_bottom;
479    Histo1DPtr _h_Cparameter_udsc, _h_Cparameter_bottom;
480    Histo1DPtr _h_Dparameter_udsc, _h_Dparameter_bottom;
481    Histo1DPtr _h_Ncharged, _h_Ncharged_udsc, _h_Ncharged_bottom;
482    Histo1DPtr _h_scaledMomentum, _h_scaledMomentum_udsc, _h_scaledMomentum_bottom;
483    // at other enegies
484    Histo1DPtr _h_thrust, _h_rho, _h_B_T, _h_B_W, _h_C, _h_D, _h_N, _h_xi;
485    Histo1DPtr _h_y_2_JADE,_h_y_3_JADE,_h_y_4_JADE,_h_y_5_JADE;
486    map<string, BinnedHistoPtr<string>> _s;
487    /// @}
488
489  };
490
491
492  RIVET_DECLARE_PLUGIN(L3_2004_I652683);
493
494}