Decayer with alice decay options using Pythia8
[u/mrichter/AliRoot.git] / PYTHIA8 / AliDecayerPythia8.cxx
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
26 ClassImp(AliDecayerPythia8)
27
28 AliDecayerPythia8::AliDecayerPythia8():
29     TPythia8Decayer(),
30     fDecay(kAll),
31     fHeavyFlavour(kTRUE)
32 {
33     // Constructor
34 }
35
36 void 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
267 Float_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
291 Int_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
303 void 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
328 void 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
359 Float_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
368 void  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
381 void 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