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.6  1999/11/09 07:38:48  fca
  19 Changes for compatibility with version 2.23 of ROOT
  20
  21 Revision 1.5  1999/11/03 17:43:20  fca
  22 New version from G.Martinez & A.Morsch
  23
  24 Revision 1.4  1999/09/29 09:24:14  fca
  25 Introduction of the Copyright and cvs Log
  26
  27 */
  28
  29
  30 #include "AliPythia.h"
  31 #include "AliRun.h"
  32
  33 ClassImp(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
  43 extern "C" void type_of_call
  44           lu1ent(Int_t&, Int_t&, Float_t&, Float_t&, Float_t&);
  45
  46
  47
  48 //_____________________________________________________________________________
  49
  50 Int_t AliPythia::fgInit=0;
  51
  52 AliPythia::AliPythia()
  53 {
  54 // Constructor
  55     for (Int_t i=0; i<501; i++) {
  56         fGPCode[i][0]=0;
  57         fGPCode[i][1]=0;
  58     }
  59 }
  60
  61 void  AliPythia::Lu1Ent(Int_t flag, Int_t idpart,
  62                       Float_t mom, Float_t theta,Float_t phi)
  63 {
  64 // Wrap of Pythia lu1ent subroutine
  65   printf("%d %d %f %f %f\n",flag, idpart, mom, theta, phi);
  66   lu1ent(flag, idpart, mom, theta, phi);
  67
  68 }
  69 void AliPythia::DecayParticle(Int_t idpart,
  70                               Float_t mom, Float_t theta,Float_t phi)
  71 {
  72 // Decay a particle
  73     Lu1Ent(0, idpart, mom, theta, phi);
  74     GetPrimaries();
  75 }
  76
  77 void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
  78 {
  79 // Initialise the process to generate
  80     fProcess = process;
  81     fEcms = energy;
  82     fStrucFunc = strucfunc;
  83 //  don't decay p0
  84     SetMDCY(Lucomp(111),1,0);
  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;
 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);
 171     }
 172 //
 173 //  Initialize PYTHIA
 174     Initialize("CMS","p","p",fEcms);
 175 }
 176
 177 Int_t AliPythia::CountProducts(Int_t channel, Int_t particle)
 178 {
 179 // Count number of decay products
 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
 187 void AliPythia::AllowAllDecays()
 188 {
 189 // Reset decay flags
 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
 200 void AliPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
 201 {
 202 //
 203 //  force decay of particle into products with multiplicity mult
 204
 205     Int_t kc=Lucomp(particle);
 206     SetMDCY(kc,1,1);
 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
 222 void AliPythia::ForceDecay(Decay_t decay)
 223 {
 224 // Force a particle decay mode
 225     fDecay=decay;
 226 //
 227 // Make clean
 228 // AllowAllDecays();
 229 //
 230 // select mode
 231
 232     switch (decay)
 233     {
 234     case semimuonic:
 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         }
 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;
 303     case pitomu:
 304         ForceParticleDecay(211,13,1); // pi->mu
 305         break;
 306     case katomu:
 307         ForceParticleDecay(321,13,1); // K->mu
 308         break;
 309     case all:
 310     case nodecay:
 311         break;
 312     }
 313 }
 314
 315
 316     void AliPythia::DefineParticles()
 317 {
 318 // Define new particles
 319     if (fgInit) return;
 320     fgInit=1;
 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
 331     kc  = Lucomp(333);
 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 //
 371 // phi clone for dilepton decay-channel
 372     kc=Lucomp(41);
 373     mass =GetPMAS(kc,1);
 374     tlife=GetPMAS(kc,4);
 375     gMC->Gspart(113,"Phi",3,mass,0,tlife);
 376     fGPCode[kc][0]=113;
 377     fGPCode[kc][1]=113;
 378     // J/Psi
 379     kc=Lucomp(443);
 380     mass =GetPMAS(kc,1);
 381     tlife=GetPMAS(kc,4);
 382     gMC->Gspart(114,"J/Psi",3,mass,0,tlife);
 383     fGPCode[kc][0]=114;
 384     fGPCode[kc][1]=114;
 385     // psi prime
 386     kc=Lucomp(30443);
 387     mass =GetPMAS(kc,1);
 388     tlife=GetPMAS(kc,4);
 389     gMC->Gspart(115,"Psi'",3,mass,0,tlife);
 390     fGPCode[kc][0]=115;
 391     fGPCode[kc][1]=115;
 392     // upsilon(1s)
 393     kc=Lucomp(553);
 394     mass =GetPMAS(kc,1);
 395     tlife=GetPMAS(kc,4);
 396     gMC->Gspart(116,"Upsilon",3,mass,0,tlife);
 397     fGPCode[kc][0]=116;
 398     fGPCode[kc][1]=116;
 399     // upsilon(2s)
 400     kc=Lucomp(30553);
 401     mass =GetPMAS(kc,1);
 402     tlife=GetPMAS(kc,4);
 403     gMC->Gspart(117,"Upsilon'",3,mass,0,tlife);
 404     fGPCode[kc][0]=117;
 405     fGPCode[kc][1]=117;
 406     // upsilon(3s)
 407     kc=Lucomp(30553);
 408     mass =GetPMAS(kc,1);
 409     tlife=GetPMAS(kc,4);
 410     gMC->Gspart(118,"Upsilon''",3,mass,0,tlife);
 411     fGPCode[kc][0]=118;
 412     fGPCode[kc][1]=118;
 413 //
 414 // charmed mesons
 415 //
 416     //  D^+/-
 417     kc=Lucomp(411);
 418     mass =GetPMAS(kc,1);
 419     tlife=GetPMAS(kc,4);
 420     gMC->Gspart(119,"D^+",3,mass, 1,tlife);
 421     gMC->Gspart(120,"D^-",3,mass,-1,tlife);
 422     fGPCode[kc][0]=119;
 423     fGPCode[kc][1]=120;
 424     // D^0
 425     kc=Lucomp(421);
 426     mass =GetPMAS(kc,1);
 427     tlife=GetPMAS(kc,4);
 428     gMC->Gspart(121,"D^0",3,mass,0,tlife);
 429     gMC->Gspart(122,"D^0bar",3,mass,0,tlife);
 430     fGPCode[kc][0]=121;
 431     fGPCode[kc][1]=122;
 432     // D_s
 433     kc=Lucomp(431);
 434     mass =GetPMAS(kc,1);
 435     tlife=GetPMAS(kc,4);
 436     gMC->Gspart(123,"D_s^+",3,mass, 1,tlife);
 437     gMC->Gspart(124,"D_s^-",3,mass,-1,tlife);
 438     fGPCode[kc][0]=123;
 439     fGPCode[kc][1]=124;
 440     // Lambda_c
 441     kc=Lucomp(4122);
 442     mass =GetPMAS(kc,1);
 443     tlife=GetPMAS(kc,4);
 444     gMC->Gspart(125,"Lambda_c+",3,mass, 1,tlife);
 445     gMC->Gspart(126,"Lambda_c-",3,mass,-1,tlife);
 446     fGPCode[kc][0]=125;
 447     fGPCode[kc][1]=126;
 448     //
 449     //  beauty mesons
 450     // B_0
 451     kc=Lucomp(511);
 452     mass =GetPMAS(kc,1);
 453     tlife=GetPMAS(kc,4);
 454     gMC->Gspart(127,"B^0",3,mass, 0,tlife);
 455     gMC->Gspart(128,"B^0bar",3,mass, 0,tlife);
 456     fGPCode[kc][0]=127;
 457     fGPCode[kc][1]=128;
 458     // B^+-
 459     kc=Lucomp(521);
 460     mass =GetPMAS(kc,1);
 461     tlife=GetPMAS(kc,4);
 462     gMC->Gspart(129,"B^+",3,mass, 1,tlife);
 463     gMC->Gspart(130,"B^-",3,mass,-1,tlife);
 464     fGPCode[kc][0]=129;
 465     fGPCode[kc][1]=130;
 466     // B_s
 467     kc=Lucomp(531);
 468     mass =GetPMAS(kc,1);
 469     tlife=GetPMAS(kc,4);
 470     gMC->Gspart(131,"B_s",3,mass, 0,tlife);
 471     gMC->Gspart(132,"B_s^bar",3,mass,0,tlife);
 472     fGPCode[kc][0]=131;
 473     fGPCode[kc][1]=132;
 474     // Lambda_b
 475     kc=Lucomp(5122);
 476     mass =GetPMAS(kc,1);
 477     tlife=GetPMAS(kc,4);
 478     gMC->Gspart(133,"Lambda_b",3,mass, 0,tlife);
 479     gMC->Gspart(134,"Lambda_b^bar",3,mass,0,tlife);
 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
 486     kc=Lucomp(22);  // gamma
 487     fGPCode[kc][0]=1;
 488     fGPCode[kc][1]=1;
 489
 490     kc=Lucomp(11);  // positron
 491     fGPCode[kc][0]=2;
 492     fGPCode[kc][1]=3;
 493
 494     kc=Lucomp(12);  // neutrino
 495     fGPCode[kc][0]=4;
 496     fGPCode[kc][1]=4;
 497
 498     kc=Lucomp(13);  // muon
 499     fGPCode[kc][0]=5;
 500     fGPCode[kc][1]=6;
 501
 502     kc=Lucomp(111); // pi0
 503     fGPCode[kc][0]=7;
 504     fGPCode[kc][1]=7;
 505
 506     kc=Lucomp(211); // pi+
 507     fGPCode[kc][0]=8;
 508     fGPCode[kc][1]=9;
 509
 510     kc=Lucomp(130); // K0 short
 511     fGPCode[kc][0]=10;
 512     fGPCode[kc][1]=10;
 513
 514     kc=Lucomp(321); // K+/-
 515     fGPCode[kc][0]=11;
 516     fGPCode[kc][1]=12;
 517
 518     kc=Lucomp(2112); // neutron/anti-neutron
 519     fGPCode[kc][0]=13;
 520     fGPCode[kc][1]=25;
 521
 522     kc=Lucomp(2212); // proton/anti-proton
 523     fGPCode[kc][0]=14;
 524     fGPCode[kc][1]=15;
 525
 526     kc=Lucomp(310);  // K0 short
 527     fGPCode[kc][0]=16;
 528     fGPCode[kc][1]=16;
 529
 530     kc=Lucomp(221);  // eta
 531     fGPCode[kc][0]=17;
 532     fGPCode[kc][1]=17;
 533
 534     kc=Lucomp(3122); // lambda
 535     fGPCode[kc][0]=18;
 536     fGPCode[kc][1]=18;
 537
 538     kc=Lucomp(3222); // sigma+/antisigma+
 539     fGPCode[kc][0]=19;
 540     fGPCode[kc][1]=29;
 541
 542     kc=Lucomp(3212); // sigma0/antisigma0
 543     fGPCode[kc][0]=20;
 544     fGPCode[kc][1]=28;
 545
 546     kc=Lucomp(3112); // sigma-/antisigma-
 547     fGPCode[kc][0]=21;
 548     fGPCode[kc][1]=27;
 549
 550     kc=Lucomp(3322); // xsi0-/antixsi0
 551     fGPCode[kc][0]=22;
 552     fGPCode[kc][1]=30;
 553
 554     kc=Lucomp(3312); // xsi-/antixsi+
 555     fGPCode[kc][0]=23;
 556     fGPCode[kc][1]=31;
 557
 558     kc=Lucomp(3334); // omega/antiomega
 559     fGPCode[kc][0]=24;
 560     fGPCode[kc][1]=32;
 561 }
 562
 563
 564
 565 Int_t  AliPythia::GetGeantCode(Int_t kf)
 566 {
 567 // Get geant code
 568     Int_t kc=Lucomp(TMath::Abs(kf));
 569     return (kf > 0) ? fGPCode[kc][0] : fGPCode[kc][1];
 570 }
 571
 572 Float_t  AliPythia::GetBraPart(Int_t kf)
 573 {
 574 // Get branching ratio
 575     Int_t kc=Lucomp(TMath::Abs(kf));
 576     return fBraPart[kc];
 577 }
 578
 579 Int_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 }
 586
 587
 588
 589
 590
 591