Introduction of the Copyright and cvs Log
[u/mrichter/AliRoot.git] / EVGEN / AliPythia.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 /*
17 $Log$
18 */
19
20
21 #include "AliPythia.h"
22 #include "AliMC.h"
23 ClassImp(AliPythia)
24
25 #ifndef WIN32
26 # define lu1ent lu1ent_
27 # define type_of_call
28 #else
29 # define lu1ent LU1ENT
30 # define type_of_call _stdcall
31 #endif
32
33 extern "C" void type_of_call 
34           lu1ent(Int_t&, Int_t&, Float_t&, Float_t&, Float_t&);
35
36
37
38 //_____________________________________________________________________________
39
40 Int_t AliPythia::fgInit=0;
41
42 AliPythia::AliPythia()
43 {
44     for (Int_t i=0; i<501; i++) {
45         fGPCode[i][0]=0; 
46         fGPCode[i][1]=0;
47     }
48 }
49
50 void  AliPythia::Lu1Ent(Int_t flag, Int_t idpart, 
51                       Float_t mom, Float_t theta,Float_t phi)
52 {
53   printf("%d %d %f %f %f\n",flag, idpart, mom, theta, phi);
54   lu1ent(flag, idpart, mom, theta, phi);
55
56 }
57 void AliPythia::DecayParticle(Int_t idpart, 
58                               Float_t mom, Float_t theta,Float_t phi)
59 {
60     Lu1Ent(0, idpart, mom, theta, phi);
61     GetPrimaries();
62 }
63
64 void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
65 {
66     fProcess = process;
67     fEcms = energy;
68     fStrucFunc = strucfunc;
69 //  don't decay p0
70     SetMDCY(LuComp(111),1,0);
71 //  select structure function 
72     SetMSTP(52,2);
73     SetMSTP(51,strucfunc);
74 //
75 // Pythia initialisation for selected processes//
76 //
77 // Make MSEL clean
78 //
79     for (Int_t i=1; i<= 200; i++) {
80         SetMSUB(i,0);
81     }
82 //  select charm production
83     switch (process) 
84     {
85     case charm:
86         SetMSEL(4);
87 //
88 //  heavy quark masses
89
90         SetPMAS(4,1,1.2);
91
92 //
93 //    primordial pT
94         SetMSTP(91,1);
95         SetPARP(91,1);
96         SetPARP(93,3);
97 //
98         break;
99     case beauty:
100         SetMSEL(5);
101         SetPMAS(5,1,4.75);
102         break;
103     case jpsi:
104         SetMSEL(0);
105 // gg->J/Psi g
106         SetMSUB(86,1);
107         break;
108     case jpsi_chi:
109         SetMSEL(0);
110 // gg->J/Psi g
111         SetMSUB(86,1);
112 // gg-> chi_0c g
113         SetMSUB(87,1);
114 // gg-> chi_1c g
115         SetMSUB(88,1);
116 // gg-> chi_2c g
117         SetMSUB(89,1);  
118     case charm_unforced:
119         SetMSEL(0);
120 // gq->qg   
121         SetMSUB(28,1);
122 // gg->qq
123         SetMSUB(53,1);
124 // gg->gg
125         SetMSUB(68,1);
126     case beauty_unforced:
127         SetMSEL(0);
128 // gq->qg   
129         SetMSUB(28,1);
130 // gg->qq
131         SetMSUB(53,1);
132 // gg->gg
133         SetMSUB(68,1);
134         break;
135     case mb:
136 // Minimum Bias pp-Collisions
137 //
138 // Tuning of parameters descibed in G. Ciapetti and A. Di Ciaccio
139 // Proc. of the LHC Workshop, Aachen 1990, Vol. II p. 155
140 //   
141 //      select Pythia min. bias model
142         SetMSEL(2);
143         SetMSUB(92,1);
144         SetMSUB(93,1);
145         SetMSUB(94,1);
146         SetMSUB(95,1);  
147 //      Multiple interactions switched on
148         SetMSTP(81,1);
149         SetMSTP(82,1);
150 //      Low-pT cut-off for hard scattering
151         SetPARP(81,1.9);
152 //      model for subsequent non-hardest interaction
153 //      90% gg->gg 10% gg->qq
154         SetPARP(86,0.9);
155 //      90% of gluon interactions have minimum string length
156         SetPARP(85,0.9);
157     }
158 //
159 //  Initialize PYTHIA
160     Initialize("CMS","p","p",fEcms);
161 }
162
163 Int_t AliPythia::CountProducts(Int_t channel, Int_t particle)
164 {
165     Int_t np=0;
166     for (Int_t i=1; i<=5; i++) {
167         if (TMath::Abs(GetKFDP(channel,i)) == particle) np++;
168     }
169     return np;
170 }
171
172 void AliPythia::AllowAllDecays()
173 {
174     Int_t i;
175     for (i=1; i<= 2000; i++) {
176         SetMDME(i,1,1);
177     }
178 //
179     for (i=0; i<501; i++){
180         fBraPart[i]=1;
181     }
182 }
183
184 void AliPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
185 {
186 //
187 //  force decay of particle into products with multiplicity mult
188
189     Int_t kc=LuComp(particle);
190     SetMDCY(kc,1,1);
191     Int_t ifirst=GetMDCY(kc,2);
192     Int_t ilast=ifirst+GetMDCY(kc,3)-1;
193     fBraPart[kc] = 1;
194 //
195 //  Loop over decay channels
196     for (Int_t channel=ifirst; channel<=ilast;channel++) {
197         if (CountProducts(channel,product) >= mult) {
198             SetMDME(channel,1,1);
199         } else {
200             SetMDME(channel,1,0);
201             fBraPart[kc]-=GetBRAT(channel);
202         }
203     }
204 }
205
206 void AliPythia::ForceDecay(Decay_t decay)
207 {
208     fDecay=decay;
209 //
210 // Make clean
211 // AllowAllDecays();
212 //
213 // select mode    
214
215     switch (decay) 
216     {
217     case semimuonic:
218         if (fProcess==charm || fProcess == charm_unforced) {
219             ForceParticleDecay(  411,13,1); // D+/-     
220             ForceParticleDecay(  421,13,1); // D0     
221             ForceParticleDecay(  431,13,1); // D_s     
222             ForceParticleDecay( 4122,13,1); // Lambda_c     
223         }
224         if (fProcess==beauty || fProcess == beauty_unforced) {
225             ForceParticleDecay(  511,13,1); // B0     
226             ForceParticleDecay(  521,13,1); // B+/-     
227             ForceParticleDecay(  531,13,1); // B_s     
228             ForceParticleDecay( 5122,13,1); // Lambda_b    
229         }
230     break;
231     case dimuon:
232         ForceParticleDecay(   41,13,2); // phi
233         ForceParticleDecay(  443,13,2); // J/Psi
234         ForceParticleDecay(30443,13,2); // Psi'
235         ForceParticleDecay(  553,13,2); // Upsilon
236         ForceParticleDecay(30553,13,2); // Upsilon'
237         break;
238     case semielectronic:
239         
240         ForceParticleDecay(  411,11,1); // D+/-     
241         ForceParticleDecay(  421,11,1); // D0     
242         ForceParticleDecay(  431,11,1); // D_s     
243         ForceParticleDecay( 4122,11,1); // Lambda_c     
244         
245         ForceParticleDecay(  511,11,1); // B0     
246         ForceParticleDecay(  521,11,1); // B+/-     
247         ForceParticleDecay(  531,11,1); // B_s     
248         ForceParticleDecay( 5122,11,1); // Lambda_b     
249         break;
250     case dielectron:
251
252         ForceParticleDecay(   41,11,2); // phi
253         ForceParticleDecay(  443,11,2); // J/Psi
254         ForceParticleDecay(30443,11,2); // Psi'
255         ForceParticleDecay(  553,11,2); // Upsilon
256         ForceParticleDecay(30553,11,2); // Upsilon'
257         break;
258     case b_jpsi_dimuon:
259         ForceParticleDecay(  511,443,1); // B0     
260         ForceParticleDecay(  521,443,1); // B+/-     
261         ForceParticleDecay(  531,443,1); // B_s     
262         ForceParticleDecay( 5122,443,1); // Lambda_b
263         ForceParticleDecay(  443,13,2);  // J/Psi    
264         break;
265     case b_psip_dimuon:
266         ForceParticleDecay(  511,30443,1); // B0     
267         ForceParticleDecay(  521,30443,1); // B+/-     
268         ForceParticleDecay(  531,30443,1); // B_s     
269         ForceParticleDecay( 5122,30443,1); // Lambda_b 
270         ForceParticleDecay(30443,13,2);    // Psi'   
271         break;
272     case b_jpsi_dielectron:
273         ForceParticleDecay(  511,443,1); // B0     
274         ForceParticleDecay(  521,443,1); // B+/-     
275         ForceParticleDecay(  531,443,1); // B_s     
276         ForceParticleDecay( 5122,443,1); // Lambda_b
277         ForceParticleDecay(  443,11,2);  // J/Psi    
278         break;
279     case b_psip_dielectron:
280         ForceParticleDecay(  511,30443,1); // B0     
281         ForceParticleDecay(  521,30443,1); // B+/-     
282         ForceParticleDecay(  531,30443,1); // B_s     
283         ForceParticleDecay( 5122,30443,1); // Lambda_b 
284         ForceParticleDecay(30443,11,2);    // Psi'   
285         break;
286     case pitomu:
287         ForceParticleDecay(211,13,1); // pi->mu     
288         break;
289     case katomu:
290         ForceParticleDecay(321,13,1); // K->mu     
291         break;
292     }
293 }
294
295
296     void AliPythia::DefineParticles()
297 {
298     if (fgInit) return;
299     fgInit=1;
300     
301     Float_t mass;
302     Float_t tlife;
303     Int_t kc, nkc, i;
304 //
305 //
306 // Some particles cloned for rare decays     
307 //
308 //  phi-> mu+mu- and phi -> e+e-
309 //  clone the original phi
310     kc  = LuComp(333);
311     nkc = 41;
312     
313     for (i=1;i<=3;i++) {
314         SetKCHG(nkc,i,GetKCHG(kc,i));
315     }
316     
317     for (i=1;i<=4;i++) {
318         SetPMAS(nkc,i,GetPMAS(kc,i));
319     }
320     SetCHAF(nkc,GetCHAF(kc));
321     fBraPart[kc]=1;
322 //    
323 //  decay
324     SetMDCY(nkc,1,1);
325     SetMDCY(nkc,2,993);
326     SetMDCY(nkc,3,2);
327 //
328 //  phi-> e+e-
329     SetMDME(993,1,1);
330     SetMDME(993,2,0);
331     SetBRAT(993,2.99e-4);
332     SetKFDP(993,1,+11);
333     SetKFDP(993,2,-11);
334     SetKFDP(993,3,0);
335     SetKFDP(993,4,0);
336     SetKFDP(993,5,0);
337 //
338 //  phi-> mu+mu-
339     SetMDME(994,1,1);
340     SetMDME(994,2,0);
341     SetBRAT(994,2.5e-4);
342     SetKFDP(994,1,+13);
343     SetKFDP(994,2,-13);
344     SetKFDP(994,3,0);
345     SetKFDP(994,4,0);
346     SetKFDP(994,5,0);
347 //
348 //          Vector mesons
349 //
350 // phi clone for dilepton decay-channel
351     kc=LuComp(41);          
352     mass =GetPMAS(kc,1);
353     tlife=GetPMAS(kc,4);
354     gMC->Gspart(113,"Phi",3,mass,0,tlife);
355     fGPCode[kc][0]=113;
356     fGPCode[kc][1]=113;
357     // J/Psi  
358     kc=LuComp(443);
359     mass =GetPMAS(kc,1);
360     tlife=GetPMAS(kc,4);
361     gMC->Gspart(114,"J/Psi",3,mass,0,tlife);
362     fGPCode[kc][0]=114;
363     fGPCode[kc][1]=114;
364     // psi prime
365     kc=LuComp(30443);
366     mass =GetPMAS(kc,1);
367     tlife=GetPMAS(kc,4);
368     gMC->Gspart(115,"Psi'",3,mass,0,tlife);
369     fGPCode[kc][0]=115;
370     fGPCode[kc][1]=115;
371     // upsilon(1s) 
372     kc=LuComp(553);     
373     mass =GetPMAS(kc,1);
374     tlife=GetPMAS(kc,4);
375     gMC->Gspart(116,"Upsilon",3,mass,0,tlife);
376     fGPCode[kc][0]=116;
377     fGPCode[kc][1]=116;
378     // upsilon(2s)      
379     kc=LuComp(30553);
380     mass =GetPMAS(kc,1);
381     tlife=GetPMAS(kc,4);
382     gMC->Gspart(117,"Upsilon'",3,mass,0,tlife);
383     fGPCode[kc][0]=117;
384     fGPCode[kc][1]=117;
385     // upsilon(3s)      
386     kc=LuComp(30553);
387     mass =GetPMAS(kc,1);
388     tlife=GetPMAS(kc,4);
389     gMC->Gspart(118,"Upsilon''",3,mass,0,tlife);
390     fGPCode[kc][0]=118;
391     fGPCode[kc][1]=118;
392 //
393 // charmed mesons
394 //
395     //  D^+/-
396     kc=LuComp(411);
397     mass =GetPMAS(kc,1);
398     tlife=GetPMAS(kc,4);
399     gMC->Gspart(119,"D^+",3,mass, 1,tlife);
400     gMC->Gspart(120,"D^-",3,mass,-1,tlife);
401     fGPCode[kc][0]=119;
402     fGPCode[kc][1]=120;
403     // D^0
404     kc=LuComp(421);
405     mass =GetPMAS(kc,1);
406     tlife=GetPMAS(kc,4);
407     gMC->Gspart(121,"D^0",3,mass,0,tlife);
408     gMC->Gspart(122,"D^0bar",3,mass,0,tlife);
409     fGPCode[kc][0]=121;
410     fGPCode[kc][1]=122;
411     // D_s
412     kc=LuComp(431);
413     mass =GetPMAS(kc,1);
414     tlife=GetPMAS(kc,4);
415     gMC->Gspart(123,"D_s^+",3,mass, 1,tlife);
416     gMC->Gspart(124,"D_s^-",3,mass,-1,tlife);
417     fGPCode[kc][0]=123;
418     fGPCode[kc][1]=124;
419     // Lambda_c
420     kc=LuComp(4122);
421     mass =GetPMAS(kc,1);
422     tlife=GetPMAS(kc,4);
423     gMC->Gspart(125,"Lambda_c+",3,mass, 1,tlife);
424     gMC->Gspart(126,"Lambda_c-",3,mass,-1,tlife);
425     fGPCode[kc][0]=125;
426     fGPCode[kc][1]=126;
427     //
428     //  beauty mesons
429     // B_0
430     kc=LuComp(511);
431     mass =GetPMAS(kc,1);
432     tlife=GetPMAS(kc,4);
433     gMC->Gspart(127,"B^0",3,mass, 0,tlife);
434     gMC->Gspart(128,"B^0bar",3,mass, 0,tlife);
435     fGPCode[kc][0]=127;
436     fGPCode[kc][1]=128;
437     // B^+-
438     kc=LuComp(521);
439     mass =GetPMAS(kc,1);
440     tlife=GetPMAS(kc,4);
441     gMC->Gspart(129,"B^+",3,mass, 1,tlife);
442     gMC->Gspart(130,"B^-",3,mass,-1,tlife);
443     fGPCode[kc][0]=129;
444     fGPCode[kc][1]=130;
445     // B_s
446     kc=LuComp(531);
447     mass =GetPMAS(kc,1);
448     tlife=GetPMAS(kc,4);
449     gMC->Gspart(131,"B_s",3,mass, 0,tlife);
450     gMC->Gspart(132,"B_s^bar",3,mass,0,tlife);
451     fGPCode[kc][0]=131;
452     fGPCode[kc][1]=132;
453     // Lambda_b
454     kc=LuComp(5122);
455     mass =GetPMAS(kc,1);
456     tlife=GetPMAS(kc,4);
457     gMC->Gspart(133,"Lambda_b",3,mass, 0,tlife);
458     gMC->Gspart(134,"Lambda_b^bar",3,mass,0,tlife);
459     fGPCode[kc][0]=133;
460     fGPCode[kc][1]=134;
461 //
462 //          set up correspondance between standard GEANT particle codes
463 //          and PYTHIA kf
464
465     kc=LuComp(22);  // gamma
466     fGPCode[kc][0]=1;
467     fGPCode[kc][1]=1;
468     
469     kc=LuComp(11);  // positron
470     fGPCode[kc][0]=2;
471     fGPCode[kc][1]=3;
472     
473     kc=LuComp(12);  // neutrino
474     fGPCode[kc][0]=4;
475     fGPCode[kc][1]=4;
476
477     kc=LuComp(13);  // muon
478     fGPCode[kc][0]=5;
479     fGPCode[kc][1]=6;
480     
481     kc=LuComp(111); // pi0
482     fGPCode[kc][0]=7;
483     fGPCode[kc][1]=7;
484
485     kc=LuComp(211); // pi+
486     fGPCode[kc][0]=8;
487     fGPCode[kc][1]=9;
488
489     kc=LuComp(130); // K0 short
490     fGPCode[kc][0]=10;
491     fGPCode[kc][1]=10;
492
493     kc=LuComp(321); // K+/-
494     fGPCode[kc][0]=11;
495     fGPCode[kc][1]=12;
496
497     kc=LuComp(2112); // neutron/anti-neutron
498     fGPCode[kc][0]=13;
499     fGPCode[kc][1]=25;
500     
501     kc=LuComp(2212); // proton/anti-proton
502     fGPCode[kc][0]=14;
503     fGPCode[kc][1]=15;
504     
505     kc=LuComp(310);  // K0 short
506     fGPCode[kc][0]=16; 
507     fGPCode[kc][1]=16;
508
509     kc=LuComp(221);  // eta
510     fGPCode[kc][0]=17;
511     fGPCode[kc][1]=17;
512
513     kc=LuComp(3122); // lambda
514     fGPCode[kc][0]=18;
515     fGPCode[kc][1]=18;
516
517     kc=LuComp(3222); // sigma+/antisigma+
518     fGPCode[kc][0]=19;
519     fGPCode[kc][1]=29;
520
521     kc=LuComp(3212); // sigma0/antisigma0
522     fGPCode[kc][0]=20;
523     fGPCode[kc][1]=28;
524
525     kc=LuComp(3112); // sigma-/antisigma-
526     fGPCode[kc][0]=21;
527     fGPCode[kc][1]=27;
528
529     kc=LuComp(3322); // xsi0-/antixsi0
530     fGPCode[kc][0]=22;
531     fGPCode[kc][1]=30;
532
533     kc=LuComp(3312); // xsi-/antixsi+
534     fGPCode[kc][0]=23;
535     fGPCode[kc][1]=31;
536
537     kc=LuComp(3334); // omega/antiomega
538     fGPCode[kc][0]=24;
539     fGPCode[kc][1]=32;
540 }
541
542
543     
544 Int_t  AliPythia::GetGeantCode(Int_t kf)
545 {
546     Int_t kc=LuComp(TMath::Abs(kf));
547     return (kf > 0) ? fGPCode[kc][0] : fGPCode[kc][1];
548 }
549     
550 Float_t  AliPythia::GetBraPart(Int_t kf)
551 {
552     Int_t kc=LuComp(TMath::Abs(kf));
553     return fBraPart[kc];
554 }
555
556     
557
558
559
560
561
562