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