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