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