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