EVGEN/AliPythia.cxx

   1
   2 #include "AliPythia.h"
   3 #include "AliMC.h"
   4 ClassImp(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
  14 extern "C" void type_of_call
  15           lu1ent(Int_t&, Int_t&, Float_t&, Float_t&, Float_t&);
  16
  17
  18
  19 //_____________________________________________________________________________
  20
  21 Int_t AliPythia::fgInit=0;
  22
  23 AliPythia::AliPythia()
  24 {
  25     for (Int_t i=0; i<501; i++) {
  26         fGPCode[i][0]=0;
  27         fGPCode[i][1]=0;
  28     }
  29 }
  30
  31 void  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 }
  38 void 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
  45 void 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;
 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);
 138     }
 139 //
 140 //  Initialize PYTHIA
 141     Initialize("CMS","p","p",fEcms);
 142 }
 143
 144 Int_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
 153 void 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
 165 void AliPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
 166 {
 167 //
 168 //  force decay of particle into products with multiplicity mult
 169
 170     Int_t kc=LuComp(particle);
 171     SetMDCY(kc,1,1);
 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
 187 void 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:
 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         }
 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;
 267     case pitomu:
 268         ForceParticleDecay(211,13,1); // pi->mu
 269         break;
 270     case katomu:
 271         ForceParticleDecay(321,13,1); // K->mu
 272         break;
 273     }
 274 }
 275
 276
 277     void AliPythia::DefineParticles()
 278 {
 279     if (fgInit) return;
 280     fgInit=1;
 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 //
 331 // phi clone for dilepton decay-channel
 332     kc=LuComp(41);
 333     mass =GetPMAS(kc,1);
 334     tlife=GetPMAS(kc,4);
 335     gMC->Gspart(113,"Phi",3,mass,0,tlife);
 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);
 342     gMC->Gspart(114,"J/Psi",3,mass,0,tlife);
 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);
 349     gMC->Gspart(115,"Psi'",3,mass,0,tlife);
 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);
 356     gMC->Gspart(116,"Upsilon",3,mass,0,tlife);
 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);
 363     gMC->Gspart(117,"Upsilon'",3,mass,0,tlife);
 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);
 370     gMC->Gspart(118,"Upsilon''",3,mass,0,tlife);
 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);
 380     gMC->Gspart(119,"D^+",3,mass, 1,tlife);
 381     gMC->Gspart(120,"D^-",3,mass,-1,tlife);
 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);
 388     gMC->Gspart(121,"D^0",3,mass,0,tlife);
 389     gMC->Gspart(122,"D^0bar",3,mass,0,tlife);
 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);
 396     gMC->Gspart(123,"D_s^+",3,mass, 1,tlife);
 397     gMC->Gspart(124,"D_s^-",3,mass,-1,tlife);
 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);
 404     gMC->Gspart(125,"Lambda_c+",3,mass, 1,tlife);
 405     gMC->Gspart(126,"Lambda_c-",3,mass,-1,tlife);
 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);
 414     gMC->Gspart(127,"B^0",3,mass, 0,tlife);
 415     gMC->Gspart(128,"B^0bar",3,mass, 0,tlife);
 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);
 422     gMC->Gspart(129,"B^+",3,mass, 1,tlife);
 423     gMC->Gspart(130,"B^-",3,mass,-1,tlife);
 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);
 430     gMC->Gspart(131,"B_s",3,mass, 0,tlife);
 431     gMC->Gspart(132,"B_s^bar",3,mass,0,tlife);
 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);
 438     gMC->Gspart(133,"Lambda_b",3,mass, 0,tlife);
 439     gMC->Gspart(134,"Lambda_b^bar",3,mass,0,tlife);
 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
 525 Int_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
 531 Float_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