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