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