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