]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA8/AliDecayerPythia8.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update trigger efficiency values for JPsi and Upsilon standard tests with the default LUT
[u/mrichter/AliRoot.git] / PYTHIA8 / AliDecayerPythia8.cxx
CommitLineData
948d10f4 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18// Implementation of AliDecayer using Pythia8
19// Author: andreas.morsch@cern.ch
20#include <TMath.h>
21#include <TPDGCode.h>
22#include <TPythia8.h>
23#include "AliDecayerPythia8.h"
24#include "ParticleData.h"
25
26ClassImp(AliDecayerPythia8)
27
28AliDecayerPythia8::AliDecayerPythia8():
29 TPythia8Decayer(),
30 fDecay(kAll),
31 fHeavyFlavour(kTRUE)
32{
33 // Constructor
34}
35
36void AliDecayerPythia8::ForceDecay()
37{
38 //
39// Force a particle decay mode
40// Switch heavy flavour production off if requested
41 if (!fHeavyFlavour) SwitchOffHeavyFlavour();
42//
43 Decay_t decay = fDecay;
44 TPythia8::Instance()->ReadString("HadronLevel:Decay = on");
45
46 if (decay == kNoDecayHeavy) return;
47
48//
49// select mode
50 Int_t products[2];
51 Int_t mult[2];
52 Int_t products1[3];
53 Int_t mult1[3];
54
55 switch (decay)
56 {
57 case kHardMuons:
58 products1[0] = 13;
59 products1[1] = 443;
60 products1[2] = 100443;
61 mult1[0] = 1;
62 mult1[1] = 1;
63 mult1[2] = 1;
64 ForceParticleDecay( 511, products1, mult1, 3);
65 ForceParticleDecay( 521, products1, mult1, 3);
66 ForceParticleDecay( 531, products1, mult1, 3);
67 ForceParticleDecay( 5122, products1, mult1, 3);
68 ForceParticleDecay( 5132, products1, mult1, 3);
69 ForceParticleDecay( 5232, products1, mult1, 3);
70 ForceParticleDecay( 5332, products1, mult1, 3);
71 ForceParticleDecay( 100443, 443, 1); // Psi' -> J/Psi X
72 ForceParticleDecay( 443, 13, 2); // J/Psi -> mu+ mu-
73 ForceParticleDecay( 411,13,1); // D+/-
74 ForceParticleDecay( 421,13,1); // D0
75 ForceParticleDecay( 431,13,1); // D_s
76 ForceParticleDecay( 4122,13,1); // Lambda_c
77 ForceParticleDecay( 4132,13,1); // Xsi_c
78 ForceParticleDecay( 4232,13,1); // Sigma_c
79 ForceParticleDecay( 4332,13,1); // Omega_c
80 break;
81 case kChiToJpsiGammaToMuonMuon:
82 products[0] = 443;
83 products[1] = 22;
84 mult[0] = 1;
85 mult[1] = 1;
86 ForceParticleDecay( 20443, products, mult, 2); // Chi_1c -> J/Psi Gamma
87 ForceParticleDecay( 445, products, mult, 2); // Chi_2c -> J/Psi Gamma
88 ForceParticleDecay( 443, 13, 2); // J/Psi -> mu+ mu-
89 break;
90 case kChiToJpsiGammaToElectronElectron:
91 products[0] = 443;
92 products[1] = 22;
93 mult[0] = 1;
94 mult[1] = 1;
95 ForceParticleDecay( 20443, products, mult, 2); // Chi_1c -> J/Psi Gamma
96 ForceParticleDecay( 445, products, mult, 2); // Chi_2c -> J/Psi Gamma
97 ForceParticleDecay( 443, 11, 2); // J/Psi -> e+ e-
98 break;
99
100 case kBSemiMuonic:
101 ForceParticleDecay( 511,13,1); // B0
102 ForceParticleDecay( 521,13,1); // B+/-
103 ForceParticleDecay( 531,13,1); // B_s
104 ForceParticleDecay( 5122,13,1); // Lambda_b
105 ForceParticleDecay( 5132,13,1); // Xsi_b
106 ForceParticleDecay( 5232,13,1); // Sigma_b
107 ForceParticleDecay( 5332,13,1); // Omega_b
108 break;
109 case kSemiMuonic:
110 ForceParticleDecay( 411,13,1); // D+/-
111 ForceParticleDecay( 421,13,1); // D0
112 ForceParticleDecay( 431,13,1); // D_s
113 ForceParticleDecay( 4122,13,1); // Lambda_c
114 ForceParticleDecay( 4132,13,1); // Xsi_c
115 ForceParticleDecay( 4232,13,1); // Sigma_c
116 ForceParticleDecay( 4332,13,1); // Omega_c
117 ForceParticleDecay( 511,13,1); // B0
118 ForceParticleDecay( 521,13,1); // B+/-
119 ForceParticleDecay( 531,13,1); // B_s
120 ForceParticleDecay( 5122,13,1); // Lambda_b
121 ForceParticleDecay( 5132,13,1); // Xsi_b
122 ForceParticleDecay( 5232,13,1); // Sigma_b
123 ForceParticleDecay( 5332,13,1); // Omega_b
124 break;
125 case kDiMuon:
126 ForceParticleDecay( 113,13,2); // rho
127 ForceParticleDecay( 221,13,2); // eta
128 ForceParticleDecay( 223,13,2); // omega
129 ForceParticleDecay( 333,13,2); // phi
130 ForceParticleDecay( 443,13,2); // J/Psi
131 ForceParticleDecay(100443,13,2);// Psi'
132 ForceParticleDecay( 553,13,2); // Upsilon
133 ForceParticleDecay(100553,13,2);// Upsilon'
134 ForceParticleDecay(200553,13,2);// Upsilon''
135 break;
136 case kBSemiElectronic:
137 ForceParticleDecay( 511,11,1); // B0
138 ForceParticleDecay( 521,11,1); // B+/-
139 ForceParticleDecay( 531,11,1); // B_s
140 ForceParticleDecay( 5122,11,1); // Lambda_b
141 ForceParticleDecay( 5132,11,1); // Xsi_b
142 ForceParticleDecay( 5232,11,1); // Sigma_b
143 ForceParticleDecay( 5332,11,1); // Omega_b
144 break;
145 case kSemiElectronic:
146 ForceParticleDecay( 411,11,1); // D+/-
147 ForceParticleDecay( 421,11,1); // D0
148 ForceParticleDecay( 431,11,1); // D_s
149 ForceParticleDecay( 4122,11,1); // Lambda_c
150 ForceParticleDecay( 4132,11,1); // Xsi_c
151 ForceParticleDecay( 4232,11,1); // Sigma_c
152 ForceParticleDecay( 4332,11,1); // Omega_c
153 ForceParticleDecay( 511,11,1); // B0
154 ForceParticleDecay( 521,11,1); // B+/-
155 ForceParticleDecay( 531,11,1); // B_s
156 ForceParticleDecay( 5122,11,1); // Lambda_b
157 ForceParticleDecay( 5132,11,1); // Xsi_b
158 ForceParticleDecay( 5232,11,1); // Sigma_b
159 ForceParticleDecay( 5332,11,1); // Omega_b
160 break;
161 case kDiElectron:
162 ForceParticleDecay( 113,11,2); // rho
163 ForceParticleDecay( 333,11,2); // phi
164 ForceParticleDecay( 221,11,2); // eta
165 ForceParticleDecay( 223,11,2); // omega
166 ForceParticleDecay( 443,11,2); // J/Psi
167 ForceParticleDecay(100443,11,2);// Psi'
168 ForceParticleDecay( 553,11,2); // Upsilon
169 ForceParticleDecay(100553,11,2);// Upsilon'
170 ForceParticleDecay(200553,11,2);// Upsilon''
171 break;
172 case kBJpsiDiMuon:
173
174 products[0] = 443;
175 products[1] = 100443;
176 mult[0] = 1;
177 mult[1] = 1;
178
179 ForceParticleDecay( 511, products, mult, 2); // B0 -> J/Psi (Psi') X
180 ForceParticleDecay( 521, products, mult, 2); // B+/- -> J/Psi (Psi') X
181 ForceParticleDecay( 531, products, mult, 2); // B_s -> J/Psi (Psi') X
182 ForceParticleDecay( 5122, products, mult, 2); // Lambda_b -> J/Psi (Psi') X
183 ForceParticleDecay( 100443, 443, 1); // Psi' -> J/Psi X
184 ForceParticleDecay( 443,13,2); // J/Psi -> mu+ mu-
185 break;
186 case kBPsiPrimeDiMuon:
187 ForceParticleDecay( 511,100443,1); // B0
188 ForceParticleDecay( 521,100443,1); // B+/-
189 ForceParticleDecay( 531,100443,1); // B_s
190 ForceParticleDecay( 5122,100443,1); // Lambda_b
191 ForceParticleDecay(100443,13,2); // Psi'
192 break;
193 case kBJpsiDiElectron:
194 ForceParticleDecay( 511,443,1); // B0
195 ForceParticleDecay( 521,443,1); // B+/-
196 ForceParticleDecay( 531,443,1); // B_s
197 ForceParticleDecay( 5122,443,1); // Lambda_b
198 ForceParticleDecay( 443,11,2); // J/Psi
199 break;
200 case kBJpsi:
201 ForceParticleDecay( 511,443,1); // B0
202 ForceParticleDecay( 521,443,1); // B+/-
203 ForceParticleDecay( 531,443,1); // B_s
204 ForceParticleDecay( 5122,443,1); // Lambda_b
205 break;
206 case kBPsiPrimeDiElectron:
207 ForceParticleDecay( 511,100443,1); // B0
208 ForceParticleDecay( 521,100443,1); // B+/-
209 ForceParticleDecay( 531,100443,1); // B_s
210 ForceParticleDecay( 5122,100443,1); // Lambda_b
211 ForceParticleDecay(100443,11,2); // Psi'
212 break;
213 case kPiToMu:
214 ForceParticleDecay(211,13,1); // pi->mu
215 break;
216 case kKaToMu:
217 ForceParticleDecay(321,13,1); // K->mu
218 break;
219 case kAllMuonic:
220 ForceParticleDecay(211,13,1); // pi->mu
221 ForceParticleDecay(321,13,1); // K->mu
222 break;
223 case kWToMuon:
224 ForceParticleDecay( 24, 13,1); // W -> mu
225 break;
226 case kWToCharm:
227 ForceParticleDecay( 24, 4,1); // W -> c
228 break;
229 case kWToCharmToMuon:
230 ForceParticleDecay( 24, 4,1); // W -> c
231 ForceParticleDecay( 411,13,1); // D+/- -> mu
232 ForceParticleDecay( 421,13,1); // D0 -> mu
233 ForceParticleDecay( 431,13,1); // D_s -> mu
234 ForceParticleDecay( 4122,13,1); // Lambda_c
235 ForceParticleDecay( 4132,13,1); // Xsi_c
236 ForceParticleDecay( 4232,13,1); // Sigma_c
237 ForceParticleDecay( 4332,13,1); // Omega_c
238 break;
239 case kZDiMuon:
240 ForceParticleDecay( 23, 13,2); // Z -> mu+ mu-
241 break;
242 case kZDiElectron:
243 ForceParticleDecay( 23, 11,2); // Z -> e+ e-
244 break;
245 case kHadronicD:
246 ForceHadronicD(1);
247 break;
248 case kHadronicDWithout4Bodies:
249 ForceHadronicD(0);
250 break;
251 case kPhiKK:
252 ForceParticleDecay(333,321,2); // Phi->K+K-
253 break;
254 case kOmega:
255// ForceOmega();
256 case kAll:
257 break;
258 case kNoDecay:
259 TPythia8::Instance()->ReadString("HadronLevel:Decay = off");
260 break;
261 case kNoDecayHeavy:
262 case kNeutralPion:
263 break;
264 }
265}
266
267Float_t AliDecayerPythia8::GetPartialBranchingRatio(Int_t ipart)
268{
269 // Get the partial branching ration for the forced decay channels
270
271 Pythia8::Pythia* thePythia = TPythia8::Instance()->Pythia8();
272 Pythia8::ParticleDataTable table = thePythia->particleData;
273 Pythia8::ParticleDataEntry* pd = table.particleDataPtr(ipart);
274 Pythia8::DecayTable decays = pd->decay;
275
276
277 Int_t nc = decays.size();
278 Float_t br = 0.;
279//
280// Loop over decay channels
281 for (Int_t ic = 0; ic < nc; ic++) {
282 Pythia8::DecayChannel& decCh = decays[ic];
283 for (Int_t i = 1; i <= decCh.multiplicity(); i++) {
284 br += decCh.bRatio();
285 }
286 }
287 return (br);
288}
289
290
291Int_t AliDecayerPythia8::CountProducts(Pythia8::DecayChannel& channel , Int_t particle)
292{
293// Count decay products of a given type
294 Int_t np = 0;
295 for (Int_t i = 1; i <= channel.multiplicity(); i++) {
296 if (TMath::Abs(channel.product(i)) == particle) np++;
297 }
298 return np;
299}
300
301
302
303void AliDecayerPythia8::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
304{
305//
306// Force decay of particle into products with multiplicity mult
307
308 Pythia8::Pythia* thePythia = TPythia8::Instance()->Pythia8();
309 Pythia8::ParticleDataTable table = thePythia->particleData;
310 Pythia8::ParticleDataEntry* pd = table.particleDataPtr(particle);
311 pd->setMayDecay(true);
312 Pythia8::DecayTable decays = pd->decay;
313
314
315 Int_t nc = decays.size();
316//
317// Loop over decay channels
318 for (Int_t ic = 0; ic < nc; ic++) {
319 Pythia8::DecayChannel& decCh = decays[ic];
320 if (CountProducts(decCh, product) >= mult) {
321 decCh.onMode(1);
322 } else {
323 decCh.onMode(0);
324 }
325 }
326}
327
328void AliDecayerPythia8::ForceParticleDecay(Int_t particle, Int_t* products, Int_t* mult, Int_t npart)
329{
330//
331// Force decay of particle into products with multiplicity mult
332
333 Pythia8::Pythia* thePythia = TPythia8::Instance()->Pythia8();
334 Pythia8::ParticleDataTable table = thePythia->particleData;
335 Pythia8::ParticleDataEntry* pd = table.particleDataPtr(particle);
336 pd->setMayDecay(true);
337 Pythia8::DecayTable decays = pd->decay;
338
339 Int_t nc = decays.size();
340//
341// Loop over decay channels
342 for (Int_t ic = 0; ic < nc; ic++) {
343 Int_t nprod = 0;
344 Pythia8::DecayChannel& decCh = decays[ic];
345
346 for (Int_t i = 0; i < npart; i++) {
347 nprod += (CountProducts(decCh, products[i]) >= mult[i]);
348 }
349
350 if (nprod) {
351 decCh.onMode(1);
352 } else {
353 decCh.onMode(0);
354 }
355 }
356}
357
358
359Float_t AliDecayerPythia8::GetLifetime(Int_t kf)
360{
361 // Return lifetime of particle
362 Pythia8::Pythia* thePythia = TPythia8::Instance()->Pythia8();
363 Pythia8::ParticleDataTable table = thePythia->particleData;
364 Float_t tau = table.tau0(kf);
365 return ( tau);
366}
367
368void AliDecayerPythia8::SwitchOffHeavyFlavour()
369{
370 // Switch off heavy flavour production
371 //
372// Maximum number of quark flavours used in pdf
373 TPythia8::Instance()->ReadString("PDFinProcess:nQuarkIn = 3");
374// Maximum number of flavors that can be used in showers
375 TPythia8::Instance()->ReadString("SpaceShower:nQuarkIn = 3");
376 TPythia8::Instance()->ReadString("TimeShower:nGammaToQuark = 3");
377 TPythia8::Instance()->ReadString("TimeShower:nGluonToQuark = 3");
378}
379
380
381void AliDecayerPythia8::ForceHadronicD(Int_t optUse4Bodies)
382{
383//
384// Force golden D decay modes
385//
386 const Int_t kNHadrons = 5;
387 Int_t hadron[kNHadrons] = {411, 421, 431, 4112, 4122};
388
389 // for D+ -> K0* (-> K- pi+) pi+
390 Int_t iKstar0 = 313;
391 Int_t iKstarbar0 = -313;
392 Int_t products[2] = {kKPlus, kPiMinus}, mult[2] = {1, 1};
393 ForceParticleDecay(iKstar0, products, mult, 2);
394 // for Ds -> Phi pi+
395 Int_t iPhi = 333;
396 ForceParticleDecay(iPhi, kKPlus, 2); // Phi->K+K-
397 // for D0 -> rho0 pi+ k-
398 Int_t iRho0=113;
399 ForceParticleDecay(iRho0, kPiPlus, 2); // Rho0->pi+pi-
400 // for Lambda_c -> Delta++ K-
401 Int_t iDeltaPP = 2224;
402 Int_t productsD[2] = {kProton, kPiPlus}, multD[2] = {1, 1};
403 ForceParticleDecay(iDeltaPP, productsD, multD, 2);
404
405
406 Int_t decayP1[kNHadrons][4] =
407 {
408 {kKMinus, kPiPlus, kPiPlus, 0},
409 {kKMinus, kPiPlus, 0 , 0},
410 {kKPlus , iKstarbar0, 0 , 0},
411 {-1 , -1 , -1 , -1},
412 {kProton, iKstarbar0, 0 , 0}
413 };
414 Int_t decayP2[kNHadrons][4] =
415 {
416 {iKstarbar0, kPiPlus, 0 , 0},
417 {kKMinus , kPiPlus, kPiPlus, kPiMinus},
418 {iPhi , kPiPlus, 0 , 0},
419 {-1 , -1 , -1 , -1},
420 {iDeltaPP , kKMinus, 0 , 0}
421 };
422 Int_t decayP3[kNHadrons][4] =
423 {
424 {-1 , -1 , -1 , -1},
425 {kKMinus , kPiPlus, iRho0 , 0 },
426 {-1 , -1 , -1 , -1},
427 {-1 , -1 , -1 , -1},
428 {kProton , kKMinus, kPiPlus , 0}
429 };
430 if(optUse4Bodies==0){
431 for(Int_t iDau=0;iDau<4;iDau++){
432 decayP2[1][iDau]=-1;
433 decayP3[1][iDau]=-1;
434 }
435 }
436
437
438 Pythia8::Pythia* thePythia = TPythia8::Instance()->Pythia8();
439 Pythia8::ParticleDataTable table = thePythia->particleData;
440
441 for (Int_t ihadron = 0; ihadron < kNHadrons; ihadron++)
442 {
443 Pythia8::ParticleDataEntry* pd = table.particleDataPtr(hadron[ihadron]);
444 pd->setMayDecay(true);
445 Pythia8::DecayTable decays = pd->decay;
446
447
448 for (Int_t ic = 0; ic < decays.size(); ic++) {
449 Pythia8::DecayChannel& decCh = decays[ic];
450 if ((
451 decCh.product(0) == decayP1[ihadron][0] &&
452 decCh.product(1) == decayP1[ihadron][1] &&
453 decCh.product(2) == decayP1[ihadron][2] &&
454 decCh.product(3) == decayP1[ihadron][3] &&
455 decCh.product(4) == 0
456 )
457 || (
458 decCh.product(0) == decayP2[ihadron][0] &&
459 decCh.product(1) == decayP2[ihadron][1] &&
460 decCh.product(2) == decayP2[ihadron][2] &&
461 decCh.product(3) == decayP2[ihadron][3] &&
462 decCh.product(4) == 0
463 )
464 || (
465 decCh.product(0) == decayP3[ihadron][0] &&
466 decCh.product(1) == decayP3[ihadron][1] &&
467 decCh.product(2) == decayP3[ihadron][2] &&
468 decCh.product(3) == decayP3[ihadron][3] &&
469 decCh.product(4) == 0
470 ))
471 {
472 decCh.onMode(1);
473 } else {
474 decCh.onMode(0);
475 } // selected channel ?
476 } // decay channels
477 } // hadrons
478}
479