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