EVGEN/AliPythia.cxx

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