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