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.23  2002/05/06 07:17:29  morsch
  19 Pyr gives random number r in interval 0 < r < 1.
  20
  21 Revision 1.22  2002/04/26 10:28:48  morsch
  22 Option kPyBeautyPbMNR added (N. Carrer).
  23
  24 Revision 1.21  2002/03/25 14:46:16  morsch
  25 Case  kPyD0PbMNR added (N. Carrer).
  26
  27 Revision 1.20  2002/03/03 13:48:50  morsch
  28 Option  kPyCharmPbMNR added. Produce charm pairs in agreement with MNR
  29 NLO calculations (Nicola Carrer).
  30
  31 Revision 1.19  2002/02/20 08:52:20  morsch
  32 Correct documentation of SetNuclei method.
  33
  34 Revision 1.18  2002/02/07 10:43:06  morsch
  35 Tuned pp-min.bias settings (M.Monteno, R.Ugoccioni and N.Carrer)
  36
  37 Revision 1.17  2001/12/19 15:40:43  morsch
  38 For kPyJets enforce simple jet topology, i.e no initial or final state
  39 gluon radiation and no primordial pT.
  40
  41 Revision 1.16  2001/10/12 11:13:59  morsch
  42 Missing break statements added (thanks to Nicola Carrer)
  43
  44 Revision 1.15  2001/03/27 10:54:50  morsch
  45 Add ResetDecayTable() and SsetDecayTable() methods.
  46
  47 Revision 1.14  2001/03/09 13:03:40  morsch
  48 Process_t and Struc_Func_t moved to AliPythia.h
  49
  50 Revision 1.13  2000/12/18 08:55:35  morsch
  51 Make AliPythia dependent generartors work with new scheme of random number generation
  52
  53 Revision 1.12  2000/11/30 07:12:50  alibrary
  54 Introducing new Rndm and QA classes
  55
  56 Revision 1.11  2000/10/20 06:30:06  fca
  57 Use version 0 to avoid streamer generation
  58
  59 Revision 1.10  2000/10/06 14:18:44  morsch
  60 Upper cut of prim. pT distribution set to 5. GeV
  61
  62 Revision 1.9  2000/09/18 10:41:35  morsch
  63 Add possibility to use nuclear structure functions from PDF library V8.
  64
  65 Revision 1.8  2000/09/06 14:26:24  morsch
  66 Decayer functionality of AliPythia has been moved to AliDecayerPythia.
  67 Class is now a singleton.
  68
  69 Revision 1.7  2000/06/09 20:34:50  morsch
  70 All coding rule violations except RS3 corrected
  71
  72 Revision 1.6  1999/11/09 07:38:48  fca
  73 Changes for compatibility with version 2.23 of ROOT
  74
  75 Revision 1.5  1999/11/03 17:43:20  fca
  76 New version from G.Martinez & A.Morsch
  77
  78 Revision 1.4  1999/09/29 09:24:14  fca
  79 Introduction of the Copyright and cvs Log
  80
  81 */
  82
  83
  84 #include "AliPythia.h"
  85
  86 ClassImp(AliPythia)
  87
  88 //_____________________________________________________________________________
  89
  90 AliPythia* AliPythia::fgAliPythia=NULL;
  91
  92 AliPythia::AliPythia()
  93 {
  94 // Default Constructor
  95 //
  96 //  Set random number
  97     if (!sRandom) sRandom=fRandom;
  98
  99 }
 100
 101 void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
 102 {
 103 // Initialise the process to generate
 104     fProcess = process;
 105     fEcms = energy;
 106     fStrucFunc = strucfunc;
 107 //  don't decay p0
 108     SetMDCY(Pycomp(111),1,0);
 109 //  select structure function
 110     SetMSTP(52,2);
 111     SetMSTP(51,strucfunc);
 112 //
 113 // Pythia initialisation for selected processes//
 114 //
 115 // Make MSEL clean
 116 //
 117     for (Int_t i=1; i<= 200; i++) {
 118         SetMSUB(i,0);
 119     }
 120 //  select charm production
 121     switch (process)
 122     {
 123     case kPyCharm:
 124         SetMSEL(4);
 125 //
 126 //  heavy quark masses
 127
 128         SetPMAS(4,1,1.2);
 129         SetMSTU(16,2);
 130 //
 131 //    primordial pT
 132         SetMSTP(91,1);
 133         SetPARP(91,1.);
 134         SetPARP(93,5.);
 135 //
 136         break;
 137     case kPyBeauty:
 138         SetMSEL(5);
 139         SetPMAS(5,1,4.75);
 140         SetMSTU(16,2);
 141         break;
 142     case kPyJpsi:
 143         SetMSEL(0);
 144 // gg->J/Psi g
 145         SetMSUB(86,1);
 146         break;
 147     case kPyJpsiChi:
 148         SetMSEL(0);
 149 // gg->J/Psi g
 150         SetMSUB(86,1);
 151 // gg-> chi_0c g
 152         SetMSUB(87,1);
 153 // gg-> chi_1c g
 154         SetMSUB(88,1);
 155 // gg-> chi_2c g
 156         SetMSUB(89,1);
 157         break;
 158     case kPyCharmUnforced:
 159         SetMSEL(0);
 160 // gq->qg
 161         SetMSUB(28,1);
 162 // gg->qq
 163         SetMSUB(53,1);
 164 // gg->gg
 165         SetMSUB(68,1);
 166         break;
 167     case kPyBeautyUnforced:
 168         SetMSEL(0);
 169 // gq->qg
 170         SetMSUB(28,1);
 171 // gg->qq
 172         SetMSUB(53,1);
 173 // gg->gg
 174         SetMSUB(68,1);
 175         break;
 176     case kPyMb:
 177 // Minimum Bias pp-Collisions
 178 //
 179 //
 180 //      select Pythia min. bias model
 181         SetMSEL(0);
 182         SetMSUB(92,1);      // single diffraction AB-->XB
 183         SetMSUB(93,1);      // single diffraction AB-->AX
 184         SetMSUB(94,1);      // double diffraction
 185         SetMSUB(95,1);      // low pt production
 186         SetMSTP(81,1);      // multiple interactions switched on
 187         SetMSTP(82,3);      // model with varying impact param. & a single Gaussian
 188         SetPARP(82,3.47);   // set value pT_0  for turn-off of the cross section of
 189                             // multiple interaction at a reference energy = 14000 GeV
 190         SetPARP(89,14000.); // reference energy for the above parameter
 191         SetPARP(90,0.174);  // set exponent for energy dependence of pT_0
 192     case kPyMbNonDiffr:
 193 // Minimum Bias pp-Collisions
 194 //
 195 //
 196 //      select Pythia min. bias model
 197         SetMSEL(0);
 198         SetMSUB(95,1);      // low pt production
 199         SetMSTP(81,1);      // multiple interactions switched on
 200         SetMSTP(82,3);      // model with varying impact param. & a single Gaussian
 201         SetPARP(82,3.47);   // set value pT_0  for turn-off of the cross section of
 202                             // multiple interaction at a reference energy = 14000 GeV
 203         SetPARP(89,14000.); // reference energy for the above parameter
 204         SetPARP(90,0.174);  // set exponent for energy dependence of pT_0
 205
 206         break;
 207     case kPyJets:
 208         SetMSEL(1);
 209 // no initial state radiation
 210         SetMSTP(61,0);
 211 // no final state radiation
 212         SetMSTP(71,0);
 213 // no primordial pT
 214         SetMSTP(91,0);
 215 //      SetMSTP(111,0);
 216         SetMSTU(16,1);
 217         SetMSTJ(1,1);
 218
 219         break;
 220     case kPyDirectGamma:
 221         SetMSEL(10);
 222         break;
 223     case kPyCharmPbMNR:
 224     case kPyD0PbMNR:
 225       // Tuning of Pythia parameters aimed to get a resonable agreement
 226       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 227       // c-cbar single inclusive and double differential distributions.
 228       // This parameter settings are meant to work with Pb-Pb collisions
 229       // (AliGenPythia::SetNuclei) and with kCTEQ_4L PDFs.
 230       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 231       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
 232
 233       // All QCD processes
 234       SetMSEL(1);
 235
 236       // No multiple interactions
 237       SetMSTP(81,0);
 238       SetPARP(81,0.0);
 239       SetPARP(82,0.0);
 240
 241       // Initial/final parton shower on (Pythia default)
 242       SetMSTP(61,1);
 243       SetMSTP(71,1);
 244
 245       // 2nd order alpha_s
 246       SetMSTP(2,2);
 247
 248       // QCD scales
 249       SetMSTP(32,2);
 250       SetPARP(34,1.0);
 251
 252       // Intrinsic <kT^2>
 253       SetMSTP(91,1);
 254       SetPARP(91,1.304);
 255       SetPARP(93,6.52);
 256
 257       // Set c-quark mass
 258       SetPMAS(4,1,1.2);
 259
 260       break;
 261     case kPyBeautyPbMNR:
 262       // Tuning of Pythia parameters aimed to get a resonable agreement
 263       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 264       // b-bbar single inclusive and double differential distributions.
 265       // This parameter settings are meant to work with Pb-Pb collisions
 266       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 267       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 268       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
 269
 270       // All QCD processes
 271       SetMSEL(1);
 272
 273       // No multiple interactions
 274       SetMSTP(81,0);
 275       SetPARP(81,0.0);
 276       SetPARP(82,0.0);
 277
 278       // Initial/final parton shower on (Pythia default)
 279       SetMSTP(61,1);
 280       SetMSTP(71,1);
 281
 282       // 2nd order alpha_s
 283       SetMSTP(2,2);
 284
 285       // QCD scales
 286       SetMSTP(32,2);
 287       SetPARP(34,1.0);
 288       SetPARP(67,1.0);
 289       SetPARP(71,1.0);
 290
 291       // Intrinsic <kT^2>
 292       SetMSTP(91,1);
 293       SetPARP(91,2.035);
 294       SetPARP(93,10.17);
 295
 296       // Set b-quark mass
 297       SetPMAS(5,1,4.75);
 298
 299       break;
 300     }
 301 //
 302 //  Initialize PYTHIA
 303     SetMSTP(41,1);   // all resonance decays switched on
 304
 305     Initialize("CMS","p","p",fEcms);
 306
 307 }
 308
 309 Int_t AliPythia::CheckedLuComp(Int_t kf)
 310 {
 311 // Check Lund particle code (for debugging)
 312     Int_t kc=Pycomp(kf);
 313     printf("\n Lucomp kf,kc %d %d",kf,kc);
 314     return kc;
 315 }
 316
 317 void AliPythia::SetNuclei(Int_t a1, Int_t a2)
 318 {
 319 // Treat protons as inside nuclei with mass numbers a1 and a2
 320 //    The MSTP array in the PYPARS common block is used to enable and
 321 //    select the nuclear structure functions.
 322 //    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
 323 //            =1: internal PYTHIA acording to MSTP(51)
 324 //            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
 325 //    If the following mass number both not equal zero, nuclear corrections of the stf are used.
 326 //    MSTP(192) : Mass number of nucleus side 1
 327 //    MSTP(193) : Mass number of nucleus side 2
 328     SetMSTP(52,2);
 329     SetMSTP(192, a1);
 330     SetMSTP(193, a2);
 331 }
 332
 333
 334 AliPythia* AliPythia::Instance()
 335 {
 336 // Set random number generator
 337     if (fgAliPythia) {
 338         return fgAliPythia;
 339     } else {
 340         fgAliPythia = new AliPythia();
 341         return fgAliPythia;
 342     }
 343 }
 344
 345 void AliPythia::PrintParticles()
 346 {
 347 // Print list of particl properties
 348     Int_t np = 0;
 349
 350     for (Int_t kf=0; kf<1000000; kf++) {
 351         for (Int_t c = 1;  c > -2; c-=2) {
 352
 353             Int_t kc = Pycomp(c*kf);
 354             if (kc) {
 355                 Float_t mass  = GetPMAS(kc,1);
 356                 Float_t width = GetPMAS(kc,2);
 357                 Float_t tau   = GetPMAS(kc,4);
 358
 359                 char*   name = new char[8];
 360                 Pyname(kf,name);
 361
 362                 np++;
 363
 364                 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
 365                        c*kf, name, mass, width, tau);
 366             }
 367         }
 368     }
 369     printf("\n Number of particles %d \n \n", np);
 370 }
 371
 372 void  AliPythia::ResetDecayTable()
 373 {
 374 //  Set default values for pythia decay switches
 375     Int_t i;
 376     for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
 377     for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
 378 }
 379
 380 void  AliPythia::SetDecayTable()
 381 {
 382 //  Set default values for pythia decay switches
 383 //
 384     Int_t i;
 385     for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
 386     for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
 387 }
 388
 389
 390 #ifndef WIN32
 391 #define pyr    pyr_
 392 #define pyrset pyrset_
 393 #define pyrget pyrget_
 394 #else
 395 #define pyr    PYR
 396 #define pyrset PYRSET
 397 #define pyrget PYRGET
 398 #endif
 399
 400 extern "C" {
 401   Double_t pyr(Int_t*)
 402 {
 403       Float_t r;
 404       do r=sRandom->Rndm(); while(0 >= r || r >= 1);
 405       return r;
 406 }
 407   void pyrset(Int_t*,Int_t*) {}
 408   void pyrget(Int_t*,Int_t*) {}
 409 }
 410
 411
 412
 413