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