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