rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ALEPH_1996_S3486095

Studies of QCD with the ALEPH detector.
Experiment: ALEPH (LEP 1)
Inspire ID: 428072
Status: VALIDATED
Authors:
  • Holger Schulz
References:
  • Phys. Rept., 294, 1--165 (1998)
Beams: e+ e-
Beam energies: (45.6, 45.6) GeV
Run details:
  • Hadronic Z decay events generated on the Z pole ($\sqrt{s} = 91.2$ GeV)

Summary paper of QCD results as measured by ALEPH at LEP 1. The publication includes various event shape variables, multiplicities (identified particles and inclusive), and particle spectra.

Source code: ALEPH_1996_S3486095.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/Sphericity.hh"
  5#include "Rivet/Projections/Thrust.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7#include "Rivet/Projections/ParisiTensor.hh"
  8#include "Rivet/Projections/Hemispheres.hh"
  9#include "Rivet/Projections/FinalState.hh"
 10#include "Rivet/Projections/ChargedFinalState.hh"
 11#include "Rivet/Projections/UnstableParticles.hh"
 12
 13namespace Rivet {
 14
 15
 16  /// @brief ALEPH QCD study with event shapes and identified particles
 17  ///
 18  /// @author Holger Schulz
 19  class ALEPH_1996_S3486095 : public Analysis {
 20  public:
 21
 22    RIVET_DEFAULT_ANALYSIS_CTOR(ALEPH_1996_S3486095);
 23
 24
 25    /// @name Analysis methods
 26    /// @{
 27
 28    void init() {
 29      // Set up projections
 30      declare(Beam(), "Beams");
 31      const ChargedFinalState cfs;
 32      declare(cfs, "FS");
 33      declare(UnstableParticles(), "UFS");
 34      declare(FastJets(cfs, FastJets::DURHAM, 0.7), "DurhamJets");
 35      declare(Sphericity(cfs), "Sphericity");
 36      declare(ParisiTensor(cfs), "Parisi");
 37      const Thrust thrust(cfs);
 38      declare(thrust, "Thrust");
 39      declare(Hemispheres(thrust), "Hemispheres");
 40
 41      // Book histograms
 42      book(_histSphericity   ,1, 1, 1);
 43      book(_histAplanarity   ,2, 1, 1);
 44
 45      book(_hist1MinusT      ,3, 1, 1);
 46      book(_histTMinor       ,4, 1, 1);
 47
 48      book(_histY3           ,5, 1, 1);
 49      book(_histHeavyJetMass ,6, 1, 1);
 50      book(_histCParam       ,7, 1, 1);
 51      book(_histOblateness   ,8, 1, 1);
 52
 53      book(_histScaledMom    ,9, 1, 1);
 54      book(_histRapidityT    ,10, 1, 1);
 55
 56      book(_histPtSIn        ,11, 1, 1);
 57      book(_histPtSOut       ,12, 1, 1);
 58
 59      book(_histLogScaledMom ,17, 1, 1);
 60
 61      book(_histChMult       ,18, 1, 1);
 62      book(_histMeanChMult   ,19, 1, 1);
 63
 64      book(_histMeanChMultRapt05,20, 1, 1);
 65      book(_histMeanChMultRapt10,21, 1, 1);
 66      book(_histMeanChMultRapt15,22, 1, 1);
 67      book(_histMeanChMultRapt20,23, 1, 1);
 68
 69
 70      // Particle spectra
 71      book(_histMultiPiPlus        ,25, 1, 1);
 72      book(_histMultiKPlus         ,26, 1, 1);
 73      book(_histMultiP             ,27, 1, 1);
 74      book(_histMultiPhoton        ,28, 1, 1);
 75      book(_histMultiPi0           ,29, 1, 1);
 76      book(_histMultiEta           ,30, 1, 1);
 77      book(_histMultiEtaPrime      ,31, 1, 1);
 78      book(_histMultiK0            ,32, 1, 1);
 79      book(_histMultiLambda0       ,33, 1, 1);
 80      book(_histMultiXiMinus       ,34, 1, 1);
 81      book(_histMultiSigma1385Plus ,35, 1, 1);
 82      book(_histMultiXi1530_0      ,36, 1, 1);
 83      book(_histMultiRho           ,37, 1, 1);
 84      book(_histMultiOmega782      ,38, 1, 1);
 85      book(_histMultiKStar892_0    ,39, 1, 1);
 86      book(_histMultiPhi           ,40, 1, 1);
 87
 88      book(_histMultiKStar892Plus  ,43, 1, 1);
 89
 90      // Mean multiplicities
 91      book(_histMeanMultiPi0           ,44, 1,  2);
 92      book(_histMeanMultiEta           ,44, 1,  3);
 93      book(_histMeanMultiEtaPrime      ,44, 1,  4);
 94      book(_histMeanMultiK0            ,44, 1,  5);
 95      book(_histMeanMultiRho           ,44, 1,  6);
 96      book(_histMeanMultiOmega782      ,44, 1,  7);
 97      book(_histMeanMultiPhi           ,44, 1,  8);
 98      book(_histMeanMultiKStar892Plus  ,44, 1,  9);
 99      book(_histMeanMultiKStar892_0    ,44, 1, 10);
100      book(_histMeanMultiLambda0       ,44, 1, 11);
101      book(_histMeanMultiSigma0        ,44, 1, 12);
102      book(_histMeanMultiXiMinus       ,44, 1, 13);
103      book(_histMeanMultiSigma1385Plus ,44, 1, 14);
104      book(_histMeanMultiXi1530_0      ,44, 1, 15);
105      book(_histMeanMultiOmegaOmegaBar ,44, 1, 16);
106      book(_weightedTotalPartNum, "/TMP/TotalPartNum");
107
108      book(_weightedTotalPartNum, "/TMP/weightedTotalPartNum");
109    }
110
111
112    void analyze(const Event& e) {
113      // First, veto on leptonic events by requiring at least 4 charged FS particles
114      const FinalState& fs = apply<FinalState>(e, "FS");
115      const size_t numParticles = fs.particles().size();
116
117      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
118      if (numParticles < 2) {
119        MSG_DEBUG("Failed leptonic event cut");
120        vetoEvent;
121      }
122      MSG_DEBUG("Passed leptonic event cut");
123
124      _weightedTotalPartNum->fill(numParticles);
125
126      // Get beams and average beam momentum
127      const ParticlePair& beams = apply<Beam>(e, "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      // Thrusts
133      MSG_DEBUG("Calculating thrust");
134      const Thrust& thrust = apply<Thrust>(e, "Thrust");
135      _hist1MinusT->fill(1 - thrust.thrust());
136      _histTMinor->fill(thrust.thrustMinor());
137      _histOblateness->fill(thrust.oblateness());
138
139      // Jets
140      MSG_DEBUG("Calculating differential jet rate plots:");
141      const FastJets& durjet = apply<FastJets>(e, "DurhamJets");
142      if (durjet.clusterSeq()) {
143        double y3 = durjet.clusterSeq()->exclusive_ymerge_max(2);
144        if (y3>0.0) _histY3->fill(-1. * std::log(y3));
145      }
146
147      // Sphericities
148      MSG_DEBUG("Calculating sphericity");
149      const Sphericity& sphericity = apply<Sphericity>(e, "Sphericity");
150      _histSphericity->fill(sphericity.sphericity());
151      _histAplanarity->fill(sphericity.aplanarity());
152
153      // C param
154      MSG_DEBUG("Calculating Parisi params");
155      const ParisiTensor& parisi = apply<ParisiTensor>(e, "Parisi");
156      _histCParam->fill(parisi.C());
157
158      // Hemispheres
159      MSG_DEBUG("Calculating hemisphere variables");
160      const Hemispheres& hemi = apply<Hemispheres>(e, "Hemispheres");
161      _histHeavyJetMass->fill(hemi.scaledM2high());
162
163      // Iterate over all the charged final state particles.
164      double rapt05 = 0.;
165      double rapt10 = 0.;
166      double rapt15 = 0.;
167      double rapt20 = 0.;
168      MSG_DEBUG("About to iterate over charged FS particles");
169      for (const Particle& p : fs.particles()) {
170        // Get momentum and energy of each particle.
171        const Vector3 mom3 = p.p3();
172        const double energy = p.E();
173
174        // Scaled momenta.
175        const double mom = mom3.mod();
176        const double scaledMom = mom/meanBeamMom;
177        const double logInvScaledMom = -std::log(scaledMom);
178        _histLogScaledMom->fill(logInvScaledMom);
179        _histScaledMom->fill(scaledMom);
180
181        // Get momenta components w.r.t. thrust and sphericity.
182        const double momT = dot(thrust.thrustAxis(), mom3);
183        const double pTinS = dot(mom3, sphericity.sphericityMajorAxis());
184        const double pToutS = dot(mom3, sphericity.sphericityMinorAxis());
185        _histPtSIn->fill(fabs(pTinS/GeV));
186        _histPtSOut->fill(fabs(pToutS/GeV));
187
188        // Calculate rapidities w.r.t. thrust.
189        const double rapidityT = 0.5 * std::log((energy + momT) / (energy - momT));
190        _histRapidityT->fill(fabs(rapidityT));
191        if (std::fabs(rapidityT) <= 0.5)  {
192            rapt05 += 1.0;
193        }
194        if (std::fabs(rapidityT) <= 1.0)  {
195            rapt10 += 1.0;
196        }
197        if (std::fabs(rapidityT) <= 1.5) {
198            rapt15 += 1.0;
199        }
200        if (std::fabs(rapidityT) <= 2.0)  {
201            rapt20 += 1.0;
202        }
203
204      }
205
206      _histChMult->fill(numParticles);
207
208      _histMeanChMultRapt05->fill(_histMeanChMultRapt05->bin(0).xMid(), rapt05);
209      _histMeanChMultRapt10->fill(_histMeanChMultRapt10->bin(0).xMid(), rapt10);
210      _histMeanChMultRapt15->fill(_histMeanChMultRapt15->bin(0).xMid(), rapt15);
211      _histMeanChMultRapt20->fill(_histMeanChMultRapt20->bin(0).xMid(), rapt20);
212      _histMeanChMult->fill(_histMeanChMult->bin(0).xMid(), numParticles);
213
214
215      //// Final state of unstable particles to get particle spectra
216      const UnstableParticles& ufs = apply<UnstableParticles>(e, "UFS");
217      for (Particles::const_iterator p = ufs.particles().begin(); p != ufs.particles().end(); ++p) {
218        const Vector3 mom3 = p->momentum().p3();
219        int id = abs(p->pid());
220        const double mom = mom3.mod();
221        const double energy = p->momentum().E();
222        const double scaledMom = mom/meanBeamMom;
223        const double scaledEnergy = energy/meanBeamMom;  // meanBeamMom is approximately beam energy
224        switch (id) {
225           case 22:
226              _histMultiPhoton->fill(-1.*std::log(scaledMom));
227              break;
228           case -321:
229           case 321:
230              _histMultiKPlus->fill(scaledMom);
231              break;
232           case 211:
233           case -211:
234              _histMultiPiPlus->fill(scaledMom);
235              break;
236           case 2212:
237           case -2212:
238              _histMultiP->fill(scaledMom);
239              break;
240           case 111:
241              _histMultiPi0->fill(scaledMom);
242              _histMeanMultiPi0->fill(_histMeanMultiPi0->bin(0).xMid());
243              break;
244           case 221:
245              if (scaledMom >= 0.1) {
246                _histMultiEta->fill(scaledEnergy);
247                _histMeanMultiEta->fill(_histMeanMultiEta->bin(0).xMid());
248              }
249              break;
250           case 331:
251              if (scaledMom >= 0.1) {
252                _histMultiEtaPrime->fill(scaledEnergy);
253                _histMeanMultiEtaPrime->fill(_histMeanMultiEtaPrime->bin(0).xMid());
254              }
255              break;
256           case 130: //klong
257           case 310: //kshort
258              _histMultiK0->fill(scaledMom);
259              _histMeanMultiK0->fill(_histMeanMultiK0->bin(0).xMid());
260              break;
261           case 113:
262              _histMultiRho->fill(scaledMom);
263              _histMeanMultiRho->fill(_histMeanMultiRho->bin(0).xMid());
264              break;
265           case 223:
266              _histMultiOmega782->fill(scaledMom);
267              _histMeanMultiOmega782->fill(_histMeanMultiOmega782->bin(0).xMid());
268              break;
269           case 333:
270              _histMultiPhi->fill(scaledMom);
271              _histMeanMultiPhi->fill(_histMeanMultiPhi->bin(0).xMid());
272              break;
273           case 313:
274           case -313:
275              _histMultiKStar892_0->fill(scaledMom);
276              _histMeanMultiKStar892_0->fill(_histMeanMultiKStar892_0->bin(0).xMid());
277              break;
278           case 323:
279           case -323:
280              _histMultiKStar892Plus->fill(scaledEnergy);
281              _histMeanMultiKStar892Plus->fill(_histMeanMultiKStar892Plus->bin(0).xMid());
282              break;
283           case 3122:
284           case -3122:
285              _histMultiLambda0->fill(scaledMom);
286              _histMeanMultiLambda0->fill(_histMeanMultiLambda0->bin(0).xMid());
287              break;
288           case 3212:
289           case -3212:
290              _histMeanMultiSigma0->fill(_histMeanMultiSigma0->bin(0).xMid());
291              break;
292           case 3312:
293           case -3312:
294              _histMultiXiMinus->fill(scaledEnergy);
295              _histMeanMultiXiMinus->fill(_histMeanMultiXiMinus->bin(0).xMid());
296              break;
297           case 3114:
298           case -3114:
299           case 3224:
300           case -3224:
301              _histMultiSigma1385Plus->fill(scaledEnergy);
302              _histMeanMultiSigma1385Plus->fill(_histMeanMultiSigma1385Plus->bin(0).xMid());
303              break;
304           case 3324:
305           case -3324:
306              _histMultiXi1530_0->fill(scaledEnergy);
307              _histMeanMultiXi1530_0->fill(_histMeanMultiXi1530_0->bin(0).xMid());
308              break;
309           case 3334:
310              _histMeanMultiOmegaOmegaBar->fill(_histMeanMultiOmegaOmegaBar->bin(0).xMid());
311              break;
312        }
313      }
314
315    }
316
317
318
319    /// Finalize
320    void finalize() {
321      // Normalize inclusive single particle distributions to the average number
322      // of charged particles per event.
323      const double avgNumParts = _weightedTotalPartNum->sumW() / sumOfWeights();
324
325      normalize(_histPtSIn, avgNumParts);
326      normalize(_histPtSOut, avgNumParts);
327
328      normalize(_histRapidityT, avgNumParts);
329      normalize(_histY3);
330
331      normalize(_histLogScaledMom, avgNumParts);
332      normalize(_histScaledMom, avgNumParts);
333
334      // particle spectra
335      scale(_histMultiPiPlus        ,1./sumOfWeights());
336      scale(_histMultiKPlus         ,1./sumOfWeights());
337      scale(_histMultiP             ,1./sumOfWeights());
338      scale(_histMultiPhoton        ,1./sumOfWeights());
339      scale(_histMultiPi0           ,1./sumOfWeights());
340      scale(_histMultiEta           ,1./sumOfWeights());
341      scale(_histMultiEtaPrime      ,1./sumOfWeights());
342      scale(_histMultiK0            ,1./sumOfWeights());
343      scale(_histMultiLambda0       ,1./sumOfWeights());
344      scale(_histMultiXiMinus       ,1./sumOfWeights());
345      scale(_histMultiSigma1385Plus ,1./sumOfWeights());
346      scale(_histMultiXi1530_0      ,1./sumOfWeights());
347      scale(_histMultiRho           ,1./sumOfWeights());
348      scale(_histMultiOmega782      ,1./sumOfWeights());
349      scale(_histMultiKStar892_0    ,1./sumOfWeights());
350      scale(_histMultiPhi           ,1./sumOfWeights());
351
352      scale(_histMultiKStar892Plus  ,1./sumOfWeights());
353
354      // event shape
355      normalize(_hist1MinusT);
356      normalize(_histTMinor);
357      normalize(_histOblateness);
358
359      normalize(_histSphericity);
360      normalize(_histAplanarity);
361      normalize(_histHeavyJetMass);
362      normalize(_histCParam);
363
364
365      // mean multiplicities
366      scale(_histChMult              , 2.0/sumOfWeights()); // taking into account the binwidth of 2
367      scale(_histMeanChMult          , 1.0/sumOfWeights());
368      scale(_histMeanChMultRapt05    , 1.0/sumOfWeights());
369      scale(_histMeanChMultRapt10    , 1.0/sumOfWeights());
370      scale(_histMeanChMultRapt15    , 1.0/sumOfWeights());
371      scale(_histMeanChMultRapt20    , 1.0/sumOfWeights());
372
373
374      scale(_histMeanMultiPi0          , 1.0/sumOfWeights());
375      scale(_histMeanMultiEta          , 1.0/sumOfWeights());
376      scale(_histMeanMultiEtaPrime     , 1.0/sumOfWeights());
377      scale(_histMeanMultiK0           , 1.0/sumOfWeights());
378      scale(_histMeanMultiRho          , 1.0/sumOfWeights());
379      scale(_histMeanMultiOmega782     , 1.0/sumOfWeights());
380      scale(_histMeanMultiPhi          , 1.0/sumOfWeights());
381      scale(_histMeanMultiKStar892Plus , 1.0/sumOfWeights());
382      scale(_histMeanMultiKStar892_0   , 1.0/sumOfWeights());
383      scale(_histMeanMultiLambda0      , 1.0/sumOfWeights());
384      scale(_histMeanMultiSigma0       , 1.0/sumOfWeights());
385      scale(_histMeanMultiXiMinus      , 1.0/sumOfWeights());
386      scale(_histMeanMultiSigma1385Plus, 1.0/sumOfWeights());
387      scale(_histMeanMultiXi1530_0     , 1.0/sumOfWeights());
388      scale(_histMeanMultiOmegaOmegaBar, 1.0/sumOfWeights());
389    }
390
391    /// @}
392
393
394  private:
395
396    /// Store the weighted sums of numbers of charged / charged+neutral
397    /// particles - used to calculate average number of particles for the
398    /// inclusive single particle distributions' normalisations.
399    CounterPtr _weightedTotalPartNum;
400
401    /// @name Histograms
402    /// @{
403    Histo1DPtr _histSphericity;
404    Histo1DPtr _histAplanarity;
405
406    Histo1DPtr _hist1MinusT;
407    Histo1DPtr _histTMinor;
408
409    Histo1DPtr _histY3;
410    Histo1DPtr _histHeavyJetMass;
411    Histo1DPtr _histCParam;
412    Histo1DPtr _histOblateness;
413
414    Histo1DPtr _histScaledMom;
415    Histo1DPtr _histRapidityT;
416
417    Histo1DPtr _histPtSIn;
418    Histo1DPtr _histPtSOut;
419
420    Histo1DPtr _histJetRate2Durham;
421    Histo1DPtr _histJetRate3Durham;
422    Histo1DPtr _histJetRate4Durham;
423    Histo1DPtr _histJetRate5Durham;
424
425    Histo1DPtr _histLogScaledMom;
426
427    Histo1DPtr _histChMult;
428
429    Histo1DPtr _histMultiPiPlus;
430    Histo1DPtr _histMultiKPlus;
431    Histo1DPtr _histMultiP;
432    Histo1DPtr _histMultiPhoton;
433    Histo1DPtr _histMultiPi0;
434    Histo1DPtr _histMultiEta;
435    Histo1DPtr _histMultiEtaPrime;
436    Histo1DPtr _histMultiK0;
437    Histo1DPtr _histMultiLambda0;
438    Histo1DPtr _histMultiXiMinus;
439    Histo1DPtr _histMultiSigma1385Plus;
440    Histo1DPtr _histMultiXi1530_0;
441    Histo1DPtr _histMultiRho;
442    Histo1DPtr _histMultiOmega782;
443    Histo1DPtr _histMultiKStar892_0;
444    Histo1DPtr _histMultiPhi;
445    Histo1DPtr _histMultiKStar892Plus;
446
447    // mean multiplicities
448    Histo1DPtr _histMeanChMult;
449    Histo1DPtr _histMeanChMultRapt05;
450    Histo1DPtr _histMeanChMultRapt10;
451    Histo1DPtr _histMeanChMultRapt15;
452    Histo1DPtr _histMeanChMultRapt20;
453
454    Histo1DPtr _histMeanMultiPi0;
455    Histo1DPtr _histMeanMultiEta;
456    Histo1DPtr _histMeanMultiEtaPrime;
457    Histo1DPtr _histMeanMultiK0;
458    Histo1DPtr _histMeanMultiRho;
459    Histo1DPtr _histMeanMultiOmega782;
460    Histo1DPtr _histMeanMultiPhi;
461    Histo1DPtr _histMeanMultiKStar892Plus;
462    Histo1DPtr _histMeanMultiKStar892_0;
463    Histo1DPtr _histMeanMultiLambda0;
464    Histo1DPtr _histMeanMultiSigma0;
465    Histo1DPtr _histMeanMultiXiMinus;
466    Histo1DPtr _histMeanMultiSigma1385Plus;
467    Histo1DPtr _histMeanMultiXi1530_0;
468    Histo1DPtr _histMeanMultiOmegaOmegaBar;
469    /// @}
470
471  };
472
473
474
475  RIVET_DECLARE_ALIASED_PLUGIN(ALEPH_1996_S3486095, ALEPH_1996_I428072);
476
477}