Rivet analyses referenceSLD_2004_I630327Production of $\pi^+$, $\pi^-$, $K^+$, $K^-$, $p$ and $\bar p$ in Light ($uds$), $c$ and $b$ Jets from Z DecaysExperiment: SLD (SLC) Inspire ID: 630327 Status: VALIDATED Authors:
Beam energies: (45.6, 45.6) GeV Run details:
Measurements of the differential production rates of stable charged particles in hadronic $Z^0$ decays, and of charged pions, kaons and protons identified over a wide momentum range. In addition to flavour-inclusive $Z^0$ decays, measurements are made for $Z^0$ decays into light ($u$, $d$, $s$), $c$ and $b$ primary flavors. Source code: SLD_2004_I630327.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/ChargedFinalState.hh"
6#include "Rivet/Projections/Thrust.hh"
7
8#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
9#include "Rivet/Projections/InitialQuarks.hh"
10
11namespace Rivet {
12
13
14 /// @brief SLD flavour-dependent fragmentation paper
15 ///
16 /// @author Peter Richardson
17 class SLD_2004_I630327 : public Analysis {
18 public:
19
20 /// Constructor
21 RIVET_DEFAULT_ANALYSIS_CTOR(SLD_2004_I630327);
22
23
24 /// @name Analysis methods
25 /// @{
26
27 void analyze(const Event& e) {
28 // First, veto on leptonic events by requiring at least 2 charged FS particles
29 const FinalState& fs = apply<FinalState>(e, "FS");
30 const size_t numParticles = fs.particles().size();
31
32 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
33 if (numParticles < 2) {
34 MSG_DEBUG("Failed ncharged cut");
35 vetoEvent;
36 }
37 MSG_DEBUG("Passed ncharged cut");
38
39 // Get beams and average beam momentum
40 const ParticlePair& beams = apply<Beam>(e, "Beams").beams();
41 const double meanBeamMom = ( beams.first.p3().mod() +
42 beams.second.p3().mod() ) / 2.0;
43 MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
44 int flavour = 0;
45 const InitialQuarks& iqf = apply<InitialQuarks>(e, "IQF");
46
47 // If we only have two quarks (qqbar), just take the flavour.
48 // If we have more than two quarks, look for the highest energetic q-qbar pair.
49 Particles quarks;
50 if (iqf.particles().size() == 2) {
51 flavour = iqf.particles().front().abspid();
52 quarks = iqf.particles();
53 }
54 else {
55 map<int, Particle > quarkmap;
56 for (const Particle& p : iqf.particles()) {
57 if (quarkmap.find(p.pid())==quarkmap.end())
58 quarkmap[p.pid()] = p;
59 else if (quarkmap[p.pid()].E() < p.E())
60 quarkmap[p.pid()] = p;
61 }
62 double maxenergy = 0.;
63 for (int i = 1; i <= 5; ++i) {
64 double energy(0.);
65 if(quarkmap.find( i)!=quarkmap.end())
66 energy += quarkmap[ i].E();
67 if(quarkmap.find(-i)!=quarkmap.end())
68 energy += quarkmap[-i].E();
69 if (energy > maxenergy) flavour = i;
70 }
71 if(quarkmap.find( flavour)!=quarkmap.end())
72 quarks.push_back(quarkmap[ flavour]);
73 if(quarkmap.find(-flavour)!=quarkmap.end())
74 quarks.push_back(quarkmap[-flavour]);
75 }
76
77 // total multiplicities
78 switch (flavour) {
79 case PID::DQUARK:
80 case PID::UQUARK:
81 case PID::SQUARK:
82 _weightLight ->fill();
83 _weightedTotalChargedPartNumLight ->fill(numParticles);
84 break;
85 case PID::CQUARK:
86 _weightCharm ->fill();
87 _weightedTotalChargedPartNumCharm ->fill(numParticles);
88 break;
89 case PID::BQUARK:
90 _weightBottom ->fill();
91 _weightedTotalChargedPartNumBottom ->fill(numParticles);
92 break;
93 }
94 // thrust axis for projections
95 Vector3 axis = apply<Thrust>(e, "Thrust").thrustAxis();
96 double dot(0.);
97 if(!quarks.empty()) {
98 dot = quarks[0].p3().dot(axis);
99 if(quarks[0].pid()<0) dot *= -1.;
100 }
101 // spectra and individual multiplicities
102 for (const Particle& p : fs.particles()) {
103 double pcm = p.p3().mod();
104 const double xp = pcm/meanBeamMom;
105
106 // if in quark or antiquark hemisphere
107 bool quark = p.p3().dot(axis)*dot>0.;
108
109 _h_PCharged ->fill(pcm );
110 // all charged
111 switch (flavour) {
112 case PID::DQUARK:
113 case PID::UQUARK:
114 case PID::SQUARK:
115 _h_XpChargedL->fill(xp);
116 break;
117 case PID::CQUARK:
118 _h_XpChargedC->fill(xp);
119 break;
120 case PID::BQUARK:
121 _h_XpChargedB->fill(xp);
122 break;
123 }
124
125 int id = p.abspid();
126 // charged pions
127 if (id == PID::PIPLUS) {
128 _h_XpPiPlus->fill(xp);
129 _h_XpPiPlusTotal->fill(xp);
130 switch (flavour) {
131 case PID::DQUARK:
132 case PID::UQUARK:
133 case PID::SQUARK:
134 _h_XpPiPlusL->fill(xp);
135 _h_NPiPlusL->fill(sqrtS());
136 if( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 ))
137 _h_RPiPlus->fill(xp);
138 else
139 _h_RPiMinus->fill(xp);
140 break;
141 case PID::CQUARK:
142 _h_XpPiPlusC->fill(xp);
143 _h_NPiPlusC->fill(sqrtS());
144 break;
145 case PID::BQUARK:
146 _h_XpPiPlusB->fill(xp);
147 _h_NPiPlusB->fill(sqrtS());
148 break;
149 }
150 }
151 else if (id == PID::KPLUS) {
152 _h_XpKPlus->fill(xp);
153 _h_XpKPlusTotal->fill(xp);
154 switch (flavour) {
155 case PID::DQUARK:
156 case PID::UQUARK:
157 case PID::SQUARK:
158 _h_XpKPlusL->fill(xp);
159 _h_NKPlusL->fill(sqrtS());
160 if( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 ))
161 _h_RKPlus->fill(xp);
162 else
163 _h_RKMinus->fill(xp);
164 break;
165 case PID::CQUARK:
166 _h_XpKPlusC->fill(xp);
167 _h_NKPlusC->fill(sqrtS());
168 break;
169 case PID::BQUARK:
170 _h_XpKPlusB->fill(xp);
171 _h_NKPlusB->fill(sqrtS());
172 break;
173 }
174 }
175 else if (id == PID::PROTON) {
176 _h_XpProton->fill(xp);
177 _h_XpProtonTotal->fill(xp);
178 switch (flavour) {
179 case PID::DQUARK:
180 case PID::UQUARK:
181 case PID::SQUARK:
182 _h_XpProtonL->fill(xp);
183 _h_NProtonL->fill(sqrtS());
184 if( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 ))
185 _h_RProton->fill(xp);
186 else
187 _h_RPBar ->fill(xp);
188 break;
189 case PID::CQUARK:
190 _h_XpProtonC->fill(xp);
191 _h_NProtonC->fill(sqrtS());
192 break;
193 case PID::BQUARK:
194 _h_XpProtonB->fill(xp);
195 _h_NProtonB->fill(sqrtS());
196 break;
197 }
198 }
199 }
200 }
201
202
203 void init() {
204 // Projections
205 declare(Beam(), "Beams");
206 declare(ChargedFinalState(), "FS");
207 declare(InitialQuarks(), "IQF");
208 declare(Thrust(FinalState()), "Thrust");
209
210 // Book histograms
211 book(_h_PCharged , 1, 1, 1);
212 book(_h_XpPiPlus , 2, 1, 2);
213 book(_h_XpKPlus , 3, 1, 2);
214 book(_h_XpProton , 4, 1, 2);
215 book(_h_XpPiPlusTotal , 2, 2, 2);
216 book(_h_XpKPlusTotal , 3, 2, 2);
217 book(_h_XpProtonTotal , 4, 2, 2);
218 book(_h_XpPiPlusL , 5, 1, 1);
219 book(_h_XpPiPlusC , 5, 1, 2);
220 book(_h_XpPiPlusB , 5, 1, 3);
221 book(_h_XpKPlusL , 6, 1, 1);
222 book(_h_XpKPlusC , 6, 1, 2);
223 book(_h_XpKPlusB , 6, 1, 3);
224 book(_h_XpProtonL , 7, 1, 1);
225 book(_h_XpProtonC , 7, 1, 2);
226 book(_h_XpProtonB , 7, 1, 3);
227 book(_h_XpChargedL , 8, 1, 1);
228 book(_h_XpChargedC , 8, 1, 2);
229 book(_h_XpChargedB , 8, 1, 3);
230
231 book(_h_NPiPlusL , 5, 2, 1);
232 book(_h_NPiPlusC , 5, 2, 2);
233 book(_h_NPiPlusB , 5, 2, 3);
234 book(_h_NKPlusL , 6, 2, 1);
235 book(_h_NKPlusC , 6, 2, 2);
236 book(_h_NKPlusB , 6, 2, 3);
237 book(_h_NProtonL , 7, 2, 1);
238 book(_h_NProtonC , 7, 2, 2);
239 book(_h_NProtonB , 7, 2, 3);
240
241 book(_h_RPiPlus , 9, 1, 1);
242 book(_h_RPiMinus , 9, 1, 2);
243 book(_h_RKPlus ,10, 1, 1);
244 book(_h_RKMinus ,10, 1, 2);
245 book(_h_RProton ,11, 1, 1);
246 book(_h_RPBar ,11, 1, 2);
247
248 // Ratios: used as target of divide() later
249 book(_s_PiM_PiP, 9, 1, 3);
250 book(_s_KM_KP , 10, 1, 3);
251 book(_s_Pr_PBar, 11, 1, 3);
252
253 book(_weightedTotalChargedPartNumLight, "_weightedTotalChargedPartNumLight");
254 book(_weightedTotalChargedPartNumCharm, "_weightedTotalChargedPartNumCharm");
255 book(_weightedTotalChargedPartNumBottom, "_weightedTotalChargedPartNumBottom");
256 book(_weightLight, "_weightLight");
257 book(_weightCharm, "_weightCharm");
258 book(_weightBottom, "_weightBottom");
259
260 book(tmp1, 8, 2, 1);
261 book(tmp2, 8, 2, 2);
262 book(tmp3, 8, 2, 3);
263 book(tmp4, 8, 3, 2);
264 book(tmp5, 8, 3, 3);
265
266
267 }
268
269
270 /// Finalize
271 void finalize() {
272
273 // Multiplicities
274 /// @todo Include errors
275 const double avgNumPartsLight = _weightedTotalChargedPartNumLight->val() / _weightLight->val();
276 const double avgNumPartsCharm = _weightedTotalChargedPartNumCharm->val() / _weightCharm->val();
277 const double avgNumPartsBottom = _weightedTotalChargedPartNumBottom->val() / _weightBottom->val();
278 tmp1->bin(1).set(avgNumPartsLight, 0.);
279 tmp2->bin(1).set(avgNumPartsCharm, 0.);
280 tmp3->bin(1).set(avgNumPartsBottom, 0.);
281 tmp4->bin(1).set(avgNumPartsCharm - avgNumPartsLight, 0.);
282 tmp5->bin(1).set(avgNumPartsBottom - avgNumPartsLight, 0.);
283
284 // Do divisions
285 divide(*_h_RPiMinus - *_h_RPiPlus, *_h_RPiMinus + *_h_RPiPlus, _s_PiM_PiP);
286 divide(*_h_RKMinus - *_h_RKPlus, *_h_RKMinus + *_h_RKPlus, _s_KM_KP);
287 divide(*_h_RProton - *_h_RPBar, *_h_RProton + *_h_RPBar, _s_Pr_PBar);
288
289 // Scale histograms
290 scale(_h_PCharged, 1./sumOfWeights());
291 scale(_h_XpPiPlus, 1./sumOfWeights());
292 scale(_h_XpKPlus, 1./sumOfWeights());
293 scale(_h_XpProton, 1./sumOfWeights());
294 scale(_h_XpPiPlusTotal, 1./sumOfWeights());
295 scale(_h_XpKPlusTotal, 1./sumOfWeights());
296 scale(_h_XpProtonTotal, 1./sumOfWeights());
297 scale(_h_XpPiPlusL, 1. / *_weightLight);
298 scale(_h_XpPiPlusC, 1. / *_weightCharm);
299 scale(_h_XpPiPlusB, 1. / *_weightBottom);
300 scale(_h_XpKPlusL, 1. / *_weightLight);
301 scale(_h_XpKPlusC, 1. / *_weightCharm);
302 scale(_h_XpKPlusB, 1. / *_weightBottom);
303 scale(_h_XpProtonL, 1. / *_weightLight);
304 scale(_h_XpProtonC, 1. / *_weightCharm);
305 scale(_h_XpProtonB, 1. / *_weightBottom);
306
307 scale(_h_XpChargedL, 1. / *_weightLight);
308 scale(_h_XpChargedC, 1. / *_weightCharm);
309 scale(_h_XpChargedB, 1. / *_weightBottom);
310
311 scale(_h_NPiPlusL, 1. / *_weightLight);
312 scale(_h_NPiPlusC, 1. / *_weightCharm);
313 scale(_h_NPiPlusB, 1. / *_weightBottom);
314 scale(_h_NKPlusL, 1. / *_weightLight);
315 scale(_h_NKPlusC, 1. / *_weightCharm);
316 scale(_h_NKPlusB, 1. / *_weightBottom);
317 scale(_h_NProtonL, 1. / *_weightLight);
318 scale(_h_NProtonC, 1. / *_weightCharm);
319 scale(_h_NProtonB, 1. / *_weightBottom);
320
321 // Paper suggests this should be 0.5/weight but it has to be 1.0 to get normalisations right...
322 scale(_h_RPiPlus, 1. / *_weightLight);
323 scale(_h_RPiMinus, 1. / *_weightLight);
324 scale(_h_RKPlus, 1. / *_weightLight);
325 scale(_h_RKMinus, 1. / *_weightLight);
326 scale(_h_RProton, 1. / *_weightLight);
327 scale(_h_RPBar, 1. / *_weightLight);
328
329 // convert ratio to %
330 _s_PiM_PiP->scale(100.);
331 _s_KM_KP ->scale(100.);
332 _s_Pr_PBar->scale(100.);
333 }
334
335 /// @}
336
337
338 private:
339
340 Estimate1DPtr tmp1, tmp2, tmp3, tmp4, tmp5;
341
342 /// Multiplicities
343 CounterPtr _weightedTotalChargedPartNumLight, _weightedTotalChargedPartNumCharm, _weightedTotalChargedPartNumBottom;
344
345 /// Weights
346 CounterPtr _weightLight, _weightCharm, _weightBottom;
347
348 // Histograms
349 /// @{
350 Histo1DPtr _h_PCharged;
351 Histo1DPtr _h_XpPiPlus, _h_XpKPlus, _h_XpProton;
352 Histo1DPtr _h_XpPiPlusTotal, _h_XpKPlusTotal, _h_XpProtonTotal;
353 Histo1DPtr _h_XpPiPlusL, _h_XpPiPlusC, _h_XpPiPlusB;
354 Histo1DPtr _h_XpKPlusL, _h_XpKPlusC, _h_XpKPlusB;
355 Histo1DPtr _h_XpProtonL, _h_XpProtonC, _h_XpProtonB;
356 Histo1DPtr _h_XpChargedL, _h_XpChargedC, _h_XpChargedB;
357 Histo1DPtr _h_NPiPlusL, _h_NPiPlusC, _h_NPiPlusB;
358 Histo1DPtr _h_NKPlusL, _h_NKPlusC, _h_NKPlusB;
359 Histo1DPtr _h_NProtonL, _h_NProtonC, _h_NProtonB;
360 Histo1DPtr _h_RPiPlus, _h_RPiMinus, _h_RKPlus;
361 Histo1DPtr _h_RKMinus, _h_RProton, _h_RPBar;
362 Estimate1DPtr _s_PiM_PiP, _s_KM_KP, _s_Pr_PBar;
363 /// @}
364
365 };
366
367
368
369 RIVET_DECLARE_ALIASED_PLUGIN(SLD_2004_I630327, SLD_2004_S5693039);
370
371}
|