Rivet analyses referenceALEPH_1996_I428072Studies of QCD with the ALEPH detector.Experiment: ALEPH (LEP 1) Inspire ID: 428072 Status: VALIDATED Authors:
Beam energies: (45.6, 45.6) GeV Run details:
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_I428072.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_I428072 : public Analysis {
20 public:
21
22 RIVET_DEFAULT_ANALYSIS_CTOR(ALEPH_1996_I428072);
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, JetAlg::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
63 book(_histMeanChMult ,19, 1, 1);
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(Ecms, rapt05);
209 _histMeanChMultRapt10->fill(Ecms, rapt10);
210 _histMeanChMultRapt15->fill(Ecms, rapt15);
211 _histMeanChMultRapt20->fill(Ecms, rapt20);
212 _histMeanChMult->fill(Ecms, 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(Ecms);
243 break;
244 case 221:
245 if (scaledMom >= 0.1) {
246 _histMultiEta->fill(scaledEnergy);
247 _histMeanMultiEta->fill(Ecms);
248 }
249 break;
250 case 331:
251 if (scaledMom >= 0.1) {
252 _histMultiEtaPrime->fill(scaledEnergy);
253 _histMeanMultiEtaPrime->fill(Ecms);
254 }
255 break;
256 case 130: //klong
257 case 310: //kshort
258 _histMultiK0->fill(scaledMom);
259 _histMeanMultiK0->fill(Ecms);
260 break;
261 case 113:
262 _histMultiRho->fill(scaledMom);
263 _histMeanMultiRho->fill(Ecms);
264 break;
265 case 223:
266 _histMultiOmega782->fill(scaledMom);
267 _histMeanMultiOmega782->fill(Ecms);
268 break;
269 case 333:
270 _histMultiPhi->fill(scaledMom);
271 _histMeanMultiPhi->fill(Ecms);
272 break;
273 case 313:
274 case -313:
275 _histMultiKStar892_0->fill(scaledMom);
276 _histMeanMultiKStar892_0->fill(Ecms);
277 break;
278 case 323:
279 case -323:
280 _histMultiKStar892Plus->fill(scaledEnergy);
281 _histMeanMultiKStar892Plus->fill(Ecms);
282 break;
283 case 3122:
284 case -3122:
285 _histMultiLambda0->fill(scaledMom);
286 _histMeanMultiLambda0->fill(Ecms);
287 break;
288 case 3212:
289 case -3212:
290 _histMeanMultiSigma0->fill(Ecms);
291 break;
292 case 3312:
293 case -3312:
294 _histMultiXiMinus->fill(scaledEnergy);
295 _histMeanMultiXiMinus->fill(Ecms);
296 break;
297 case 3114:
298 case -3114:
299 case 3224:
300 case -3224:
301 _histMultiSigma1385Plus->fill(scaledEnergy);
302 _histMeanMultiSigma1385Plus->fill(Ecms);
303 break;
304 case 3324:
305 case -3324:
306 _histMultiXi1530_0->fill(scaledEnergy);
307 _histMeanMultiXi1530_0->fill(Ecms);
308 break;
309 case 3334:
310 _histMeanMultiOmegaOmegaBar->fill(Ecms);
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 , 1.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 const string Ecms = "91.2";
402
403 /// @name Histograms
404 /// @{
405 Histo1DPtr _histSphericity;
406 Histo1DPtr _histAplanarity;
407
408 Histo1DPtr _hist1MinusT;
409 Histo1DPtr _histTMinor;
410
411 Histo1DPtr _histY3;
412 Histo1DPtr _histHeavyJetMass;
413 Histo1DPtr _histCParam;
414 Histo1DPtr _histOblateness;
415
416 Histo1DPtr _histScaledMom;
417 Histo1DPtr _histRapidityT;
418
419 Histo1DPtr _histPtSIn;
420 Histo1DPtr _histPtSOut;
421
422 Histo1DPtr _histJetRate2Durham;
423 Histo1DPtr _histJetRate3Durham;
424 Histo1DPtr _histJetRate4Durham;
425 Histo1DPtr _histJetRate5Durham;
426
427 Histo1DPtr _histLogScaledMom;
428
429 BinnedHistoPtr<int> _histChMult;
430
431 Histo1DPtr _histMultiPiPlus;
432 Histo1DPtr _histMultiKPlus;
433 Histo1DPtr _histMultiP;
434 Histo1DPtr _histMultiPhoton;
435 Histo1DPtr _histMultiPi0;
436 Histo1DPtr _histMultiEta;
437 Histo1DPtr _histMultiEtaPrime;
438 Histo1DPtr _histMultiK0;
439 Histo1DPtr _histMultiLambda0;
440 Histo1DPtr _histMultiXiMinus;
441 Histo1DPtr _histMultiSigma1385Plus;
442 Histo1DPtr _histMultiXi1530_0;
443 Histo1DPtr _histMultiRho;
444 Histo1DPtr _histMultiOmega782;
445 Histo1DPtr _histMultiKStar892_0;
446 Histo1DPtr _histMultiPhi;
447 Histo1DPtr _histMultiKStar892Plus;
448
449 // mean multiplicities
450 BinnedHistoPtr<string> _histMeanChMult;
451 BinnedHistoPtr<string> _histMeanChMultRapt05;
452 BinnedHistoPtr<string> _histMeanChMultRapt10;
453 BinnedHistoPtr<string> _histMeanChMultRapt15;
454 BinnedHistoPtr<string> _histMeanChMultRapt20;
455
456 BinnedHistoPtr<string> _histMeanMultiPi0;
457 BinnedHistoPtr<string> _histMeanMultiEta;
458 BinnedHistoPtr<string> _histMeanMultiEtaPrime;
459 BinnedHistoPtr<string> _histMeanMultiK0;
460 BinnedHistoPtr<string> _histMeanMultiRho;
461 BinnedHistoPtr<string> _histMeanMultiOmega782;
462 BinnedHistoPtr<string> _histMeanMultiPhi;
463 BinnedHistoPtr<string> _histMeanMultiKStar892Plus;
464 BinnedHistoPtr<string> _histMeanMultiKStar892_0;
465 BinnedHistoPtr<string> _histMeanMultiLambda0;
466 BinnedHistoPtr<string> _histMeanMultiSigma0;
467 BinnedHistoPtr<string> _histMeanMultiXiMinus;
468 BinnedHistoPtr<string> _histMeanMultiSigma1385Plus;
469 BinnedHistoPtr<string> _histMeanMultiXi1530_0;
470 BinnedHistoPtr<string> _histMeanMultiOmegaOmegaBar;
471 /// @}
472
473 };
474
475
476
477 RIVET_DECLARE_ALIASED_PLUGIN(ALEPH_1996_I428072, ALEPH_1996_S3486095);
478
479}
|