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 /* $Id$ */
  17
  18 #include "AliPythia.h"
  19 #include "AliPythiaRndm.h"
  20
  21 ClassImp(AliPythia)
  22
  23 #ifndef WIN32
  24 # define pyclus pyclus_
  25 # define pycell pycell_
  26 # define type_of_call
  27 #else
  28 # define pyclus PYCLUS
  29 # define pycell PYCELL
  30 # define type_of_call _stdcall
  31 #endif
  32
  33 extern "C" void type_of_call pyclus(Int_t & );
  34 extern "C" void type_of_call pycell(Int_t & );
  35
  36 //_____________________________________________________________________________
  37
  38 AliPythia* AliPythia::fgAliPythia=NULL;
  39
  40 AliPythia::AliPythia()
  41 {
  42 // Default Constructor
  43 //
  44 //  Set random number
  45     if (!AliPythiaRndm::GetPythiaRandom())
  46       AliPythiaRndm::SetPythiaRandom(GetRandom());
  47
  48 }
  49
  50 void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
  51 {
  52 // Initialise the process to generate
  53     if (!AliPythiaRndm::GetPythiaRandom())
  54       AliPythiaRndm::SetPythiaRandom(GetRandom());
  55
  56     fProcess = process;
  57     fEcms = energy;
  58     fStrucFunc = strucfunc;
  59 //  don't decay p0
  60     SetMDCY(Pycomp(111),1,0);
  61 //  select structure function
  62     SetMSTP(52,2);
  63     SetMSTP(51,strucfunc);
  64 //
  65 // Pythia initialisation for selected processes//
  66 //
  67 // Make MSEL clean
  68 //
  69     for (Int_t i=1; i<= 200; i++) {
  70         SetMSUB(i,0);
  71     }
  72 //  select charm production
  73     switch (process)
  74     {
  75     case kPyCharm:
  76         SetMSEL(4);
  77 //
  78 //  heavy quark masses
  79
  80         SetPMAS(4,1,1.2);
  81         SetMSTU(16,2);
  82 //
  83 //    primordial pT
  84         SetMSTP(91,1);
  85         SetPARP(91,1.);
  86         SetPARP(93,5.);
  87 //
  88         break;
  89     case kPyBeauty:
  90         SetMSEL(5);
  91         SetPMAS(5,1,4.75);
  92         SetMSTU(16,2);
  93         break;
  94     case kPyJpsi:
  95         SetMSEL(0);
  96 // gg->J/Psi g
  97         SetMSUB(86,1);
  98         break;
  99     case kPyJpsiChi:
 100         SetMSEL(0);
 101 // gg->J/Psi g
 102         SetMSUB(86,1);
 103 // gg-> chi_0c g
 104         SetMSUB(87,1);
 105 // gg-> chi_1c g
 106         SetMSUB(88,1);
 107 // gg-> chi_2c g
 108         SetMSUB(89,1);
 109         break;
 110     case kPyCharmUnforced:
 111         SetMSEL(0);
 112 // gq->qg
 113         SetMSUB(28,1);
 114 // gg->qq
 115         SetMSUB(53,1);
 116 // gg->gg
 117         SetMSUB(68,1);
 118         break;
 119     case kPyBeautyUnforced:
 120         SetMSEL(0);
 121 // gq->qg
 122         SetMSUB(28,1);
 123 // gg->qq
 124         SetMSUB(53,1);
 125 // gg->gg
 126         SetMSUB(68,1);
 127         break;
 128     case kPyMb:
 129 // Minimum Bias pp-Collisions
 130 //
 131 //
 132 //      select Pythia min. bias model
 133         SetMSEL(0);
 134         SetMSUB(92,1);      // single diffraction AB-->XB
 135         SetMSUB(93,1);      // single diffraction AB-->AX
 136         SetMSUB(94,1);      // double diffraction
 137         SetMSUB(95,1);      // low pt production
 138         SetMSTP(81,1);      // multiple interactions switched on
 139         SetMSTP(82,3);      // model with varying impact param. & a single Gaussian
 140         SetPARP(82,3.47);   // set value pT_0  for turn-off of the cross section of
 141                             // multiple interaction at a reference energy = 14000 GeV
 142         SetPARP(89,14000.); // reference energy for the above parameter
 143         SetPARP(90,0.174);  // set exponent for energy dependence of pT_0
 144     case kPyMbNonDiffr:
 145 // Minimum Bias pp-Collisions
 146 //
 147 //
 148 //      select Pythia min. bias model
 149         SetMSEL(0);
 150         SetMSUB(95,1);      // low pt production
 151         SetMSTP(81,1);      // multiple interactions switched on
 152         SetMSTP(82,3);      // model with varying impact param. & a single Gaussian
 153         SetPARP(82,3.47);   // set value pT_0  for turn-off of the cross section of
 154                             // multiple interaction at a reference energy = 14000 GeV
 155         SetPARP(89,14000.); // reference energy for the above parameter
 156         SetPARP(90,0.174);  // set exponent for energy dependence of pT_0
 157
 158         break;
 159     case kPyJets:
 160 //
 161 //  QCD Jets
 162 //
 163         SetMSEL(1);
 164         break;
 165     case kPyDirectGamma:
 166         SetMSEL(10);
 167         break;
 168     case kPyCharmPbPbMNR:
 169     case kPyD0PbPbMNR:
 170       // Tuning of Pythia parameters aimed to get a resonable agreement
 171       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 172       // c-cbar single inclusive and double differential distributions.
 173       // This parameter settings are meant to work with Pb-Pb collisions
 174       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 175       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 176       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
 177
 178       // All QCD processes
 179       SetMSEL(1);
 180
 181       // No multiple interactions
 182       SetMSTP(81,0);
 183       SetPARP(81,0.0);
 184       SetPARP(82,0.0);
 185
 186       // Initial/final parton shower on (Pythia default)
 187       SetMSTP(61,1);
 188       SetMSTP(71,1);
 189
 190       // 2nd order alpha_s
 191       SetMSTP(2,2);
 192
 193       // QCD scales
 194       SetMSTP(32,2);
 195       SetPARP(34,1.0);
 196
 197       // Intrinsic <kT>
 198       SetMSTP(91,1);
 199       SetPARP(91,1.304);
 200       SetPARP(93,6.52);
 201
 202       // Set c-quark mass
 203       SetPMAS(4,1,1.2);
 204
 205       break;
 206     case kPyCharmpPbMNR:
 207     case kPyD0pPbMNR:
 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 p-Pb collisions
 212       // (AliGenPythia::SetNuclei) and with kCTEQ4L 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>
 236       SetMSTP(91,1);
 237       SetPARP(91,1.16);
 238       SetPARP(93,5.8);
 239
 240       // Set c-quark mass
 241       SetPMAS(4,1,1.2);
 242
 243       break;
 244     case kPyCharmppMNR:
 245     case kPyD0ppMNR:
 246       // Tuning of Pythia parameters aimed to get a resonable agreement
 247       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 248       // c-cbar single inclusive and double differential distributions.
 249       // This parameter settings are meant to work with pp collisions
 250       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 251       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 252       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
 253
 254       // All QCD processes
 255       SetMSEL(1);
 256
 257       // No multiple interactions
 258       SetMSTP(81,0);
 259       SetPARP(81,0.0);
 260       SetPARP(82,0.0);
 261
 262       // Initial/final parton shower on (Pythia default)
 263       SetMSTP(61,1);
 264       SetMSTP(71,1);
 265
 266       // 2nd order alpha_s
 267       SetMSTP(2,2);
 268
 269       // QCD scales
 270       SetMSTP(32,2);
 271       SetPARP(34,1.0);
 272
 273       // Intrinsic <kT^2>
 274       SetMSTP(91,1);
 275       SetPARP(91,1.);
 276       SetPARP(93,5.);
 277
 278       // Set c-quark mass
 279       SetPMAS(4,1,1.2);
 280
 281       break;
 282     case kPyBeautyPbPbMNR:
 283       // Tuning of Pythia parameters aimed to get a resonable agreement
 284       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 285       // b-bbar single inclusive and double differential distributions.
 286       // This parameter settings are meant to work with Pb-Pb collisions
 287       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 288       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 289       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
 290
 291       // All QCD processes
 292       SetMSEL(1);
 293
 294       // No multiple interactions
 295       SetMSTP(81,0);
 296       SetPARP(81,0.0);
 297       SetPARP(82,0.0);
 298
 299       // Initial/final parton shower on (Pythia default)
 300       SetMSTP(61,1);
 301       SetMSTP(71,1);
 302
 303       // 2nd order alpha_s
 304       SetMSTP(2,2);
 305
 306       // QCD scales
 307       SetMSTP(32,2);
 308       SetPARP(34,1.0);
 309       SetPARP(67,1.0);
 310       SetPARP(71,1.0);
 311
 312       // Intrinsic <kT>
 313       SetMSTP(91,1);
 314       SetPARP(91,2.035);
 315       SetPARP(93,10.17);
 316
 317       // Set b-quark mass
 318       SetPMAS(5,1,4.75);
 319
 320       break;
 321     case kPyBeautypPbMNR:
 322       // Tuning of Pythia parameters aimed to get a resonable agreement
 323       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 324       // b-bbar single inclusive and double differential distributions.
 325       // This parameter settings are meant to work with p-Pb collisions
 326       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 327       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 328       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
 329
 330       // All QCD processes
 331       SetMSEL(1);
 332
 333       // No multiple interactions
 334       SetMSTP(81,0);
 335       SetPARP(81,0.0);
 336       SetPARP(82,0.0);
 337
 338       // Initial/final parton shower on (Pythia default)
 339       SetMSTP(61,1);
 340       SetMSTP(71,1);
 341
 342       // 2nd order alpha_s
 343       SetMSTP(2,2);
 344
 345       // QCD scales
 346       SetMSTP(32,2);
 347       SetPARP(34,1.0);
 348       SetPARP(67,1.0);
 349       SetPARP(71,1.0);
 350
 351       // Intrinsic <kT>
 352       SetMSTP(91,1);
 353       SetPARP(91,1.60);
 354       SetPARP(93,8.00);
 355
 356       // Set b-quark mass
 357       SetPMAS(5,1,4.75);
 358
 359       break;
 360     case kPyBeautyppMNR:
 361       // Tuning of Pythia parameters aimed to get a resonable agreement
 362       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 363       // b-bbar single inclusive and double differential distributions.
 364       // This parameter settings are meant to work with pp collisions
 365       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 366       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 367       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
 368
 369       // All QCD processes
 370       SetMSEL(1);
 371
 372       // No multiple interactions
 373       SetMSTP(81,0);
 374       SetPARP(81,0.0);
 375       SetPARP(82,0.0);
 376
 377       // Initial/final parton shower on (Pythia default)
 378       SetMSTP(61,1);
 379       SetMSTP(71,1);
 380
 381       // 2nd order alpha_s
 382       SetMSTP(2,2);
 383
 384       // QCD scales
 385       SetMSTP(32,2);
 386       SetPARP(34,1.0);
 387       SetPARP(67,1.0);
 388       SetPARP(71,1.0);
 389
 390       // Intrinsic <kT>
 391       SetMSTP(91,1);
 392       SetPARP(91,1.);
 393       SetPARP(93,5.);
 394
 395       // Set b-quark mass
 396       SetPMAS(5,1,4.75);
 397
 398       break;
 399     }
 400 //
 401 //  Initialize PYTHIA
 402     SetMSTP(41,1);   // all resonance decays switched on
 403
 404     Initialize("CMS","p","p",fEcms);
 405
 406 }
 407
 408 Int_t AliPythia::CheckedLuComp(Int_t kf)
 409 {
 410 // Check Lund particle code (for debugging)
 411     Int_t kc=Pycomp(kf);
 412     printf("\n Lucomp kf,kc %d %d",kf,kc);
 413     return kc;
 414 }
 415
 416 void AliPythia::SetNuclei(Int_t a1, Int_t a2)
 417 {
 418 // Treat protons as inside nuclei with mass numbers a1 and a2
 419 //    The MSTP array in the PYPARS common block is used to enable and
 420 //    select the nuclear structure functions.
 421 //    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
 422 //            =1: internal PYTHIA acording to MSTP(51)
 423 //            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
 424 //    If the following mass number both not equal zero, nuclear corrections of the stf are used.
 425 //    MSTP(192) : Mass number of nucleus side 1
 426 //    MSTP(193) : Mass number of nucleus side 2
 427     SetMSTP(52,2);
 428     SetMSTP(192, a1);
 429     SetMSTP(193, a2);
 430 }
 431
 432
 433 AliPythia* AliPythia::Instance()
 434 {
 435 // Set random number generator
 436     if (fgAliPythia) {
 437         return fgAliPythia;
 438     } else {
 439         fgAliPythia = new AliPythia();
 440         return fgAliPythia;
 441     }
 442 }
 443
 444 void AliPythia::PrintParticles()
 445 {
 446 // Print list of particl properties
 447     Int_t np = 0;
 448     char*   name = new char[16];
 449     for (Int_t kf=0; kf<1000000; kf++) {
 450         for (Int_t c = 1;  c > -2; c-=2) {
 451             Int_t kc = Pycomp(c*kf);
 452             if (kc) {
 453                 Float_t mass  = GetPMAS(kc,1);
 454                 Float_t width = GetPMAS(kc,2);
 455                 Float_t tau   = GetPMAS(kc,4);
 456
 457                 Pyname(kf,name);
 458
 459                 np++;
 460
 461                 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
 462                        c*kf, name, mass, width, tau);
 463             }
 464         }
 465     }
 466     printf("\n Number of particles %d \n \n", np);
 467 }
 468
 469 void  AliPythia::ResetDecayTable()
 470 {
 471 //  Set default values for pythia decay switches
 472     Int_t i;
 473     for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
 474     for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
 475 }
 476
 477 void  AliPythia::SetDecayTable()
 478 {
 479 //  Set default values for pythia decay switches
 480 //
 481     Int_t i;
 482     for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
 483     for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
 484 }
 485
 486 void  AliPythia::Pyclus(Int_t& njet)
 487 {
 488 //  Call Pythia clustering algorithm
 489 //
 490     pyclus(njet);
 491 }
 492
 493 void  AliPythia::Pycell(Int_t& njet)
 494 {
 495 //  Call Pythia jet reconstruction algorithm
 496 //
 497     pycell(njet);
 498 }
 499
 500
 501