]> 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
Remove unused files and add ifdef at the beginning of .h file
[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$
18*/
19
75c6d54e 20
fe4da5cc 21#include "AliPythia.h"
22#include "AliMC.h"
23ClassImp(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
33extern "C" void type_of_call
34 lu1ent(Int_t&, Int_t&, Float_t&, Float_t&, Float_t&);
35
36
37
38//_____________________________________________________________________________
39
40Int_t AliPythia::fgInit=0;
41
75c6d54e 42AliPythia::AliPythia()
43{
44 for (Int_t i=0; i<501; i++) {
45 fGPCode[i][0]=0;
46 fGPCode[i][1]=0;
47 }
48}
49
fe4da5cc 50void 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}
57void 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
64void 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;
75c6d54e 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);
fe4da5cc 157 }
158//
159// Initialize PYTHIA
160 Initialize("CMS","p","p",fEcms);
161}
162
163Int_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
172void 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
184void AliPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
185{
186//
187// force decay of particle into products with multiplicity mult
75c6d54e 188
fe4da5cc 189 Int_t kc=LuComp(particle);
75c6d54e 190 SetMDCY(kc,1,1);
fe4da5cc 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
206void 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:
75c6d54e 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 }
fe4da5cc 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;
75c6d54e 286 case pitomu:
287 ForceParticleDecay(211,13,1); // pi->mu
288 break;
289 case katomu:
290 ForceParticleDecay(321,13,1); // K->mu
291 break;
fe4da5cc 292 }
293}
294
295
296 void AliPythia::DefineParticles()
297{
298 if (fgInit) return;
299 fgInit=1;
fe4da5cc 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//
75c6d54e 350// phi clone for dilepton decay-channel
fe4da5cc 351 kc=LuComp(41);
352 mass =GetPMAS(kc,1);
353 tlife=GetPMAS(kc,4);
cfce8870 354 gMC->Gspart(113,"Phi",3,mass,0,tlife);
fe4da5cc 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);
cfce8870 361 gMC->Gspart(114,"J/Psi",3,mass,0,tlife);
fe4da5cc 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);
cfce8870 368 gMC->Gspart(115,"Psi'",3,mass,0,tlife);
fe4da5cc 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);
cfce8870 375 gMC->Gspart(116,"Upsilon",3,mass,0,tlife);
fe4da5cc 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);
cfce8870 382 gMC->Gspart(117,"Upsilon'",3,mass,0,tlife);
fe4da5cc 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);
cfce8870 389 gMC->Gspart(118,"Upsilon''",3,mass,0,tlife);
fe4da5cc 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);
cfce8870 399 gMC->Gspart(119,"D^+",3,mass, 1,tlife);
400 gMC->Gspart(120,"D^-",3,mass,-1,tlife);
fe4da5cc 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);
cfce8870 407 gMC->Gspart(121,"D^0",3,mass,0,tlife);
408 gMC->Gspart(122,"D^0bar",3,mass,0,tlife);
fe4da5cc 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);
cfce8870 415 gMC->Gspart(123,"D_s^+",3,mass, 1,tlife);
416 gMC->Gspart(124,"D_s^-",3,mass,-1,tlife);
fe4da5cc 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);
cfce8870 423 gMC->Gspart(125,"Lambda_c+",3,mass, 1,tlife);
424 gMC->Gspart(126,"Lambda_c-",3,mass,-1,tlife);
fe4da5cc 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);
cfce8870 433 gMC->Gspart(127,"B^0",3,mass, 0,tlife);
434 gMC->Gspart(128,"B^0bar",3,mass, 0,tlife);
fe4da5cc 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);
cfce8870 441 gMC->Gspart(129,"B^+",3,mass, 1,tlife);
442 gMC->Gspart(130,"B^-",3,mass,-1,tlife);
fe4da5cc 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);
cfce8870 449 gMC->Gspart(131,"B_s",3,mass, 0,tlife);
450 gMC->Gspart(132,"B_s^bar",3,mass,0,tlife);
fe4da5cc 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);
cfce8870 457 gMC->Gspart(133,"Lambda_b",3,mass, 0,tlife);
458 gMC->Gspart(134,"Lambda_b^bar",3,mass,0,tlife);
fe4da5cc 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
544Int_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
550Float_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