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