EVGEN/AliPythia.cxx

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