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