PYTHIA6/AliPythia8.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /* $Id$ */
  18 #include <TString.h>
  19 #include <TVector3.h>
  20 #include <TMath.h>
  21
  22 #include "AliPythia8.h"
  23 #include "AliLog.h"
  24 #include "AliStack.h"
  25 #include "AliPythiaRndm.h"
  26
  27
  28 ClassImp(AliPythia8)
  29
  30 // Particles produced in string fragmentation point directly to either of the two endpoints
  31 // of the string (depending in the side they were generated from).
  32 //    SetMSTU(16,2); // ????
  33 //  String drawing almost completely minimizes string length.
  34 //  Probability that an additional interaction gives two gluons
  35 //  ... with color connection to nearest neighbours
  36 //    SetPARP(85,0.9);
  37 //  ... as closed gluon loop
  38 //    SetPARP(86,0.95);
  39 // Lambda_FSR scale.
  40 //      SetPARJ(81, 0.29);
  41 // Baryon production model
  42 //      SetMSTJ(12,3);
  43 // String fragmentation
  44 //      SetMSTJ(1,1);
  45 // sea quarks can be used for baryon formation
  46 //      SetMSTP(88,2);
  47 // choice of max. virtuality for ISR
  48 //      SetMSTP(68,1);
  49 // regularisation scheme of ISR
  50 //      SetMSTP(70,2);
  51 // all resonance decays switched on
  52 //    SetMSTP(41,1);
  53 AliPythia8* AliPythia8::fgAliPythia8=NULL;
  54
  55 AliPythia8::AliPythia8():
  56     TPythia8(),
  57     AliPythiaBase(),
  58     fProcess(kPyMb),
  59     fEcms(0.),
  60     fStrucFunc(kCTEQ5L)
  61 {
  62 // Default Constructor
  63 //
  64 //  Set random number
  65     if (!AliPythiaRndm::GetPythiaRandom())
  66         AliPythiaRndm::SetPythiaRandom(GetRandom());
  67 }
  68
  69 AliPythia8::AliPythia8(const AliPythia8& pythia):
  70     TPythia8(),
  71     AliPythiaBase(),
  72     fProcess(kPyMb),
  73     fEcms(0.),
  74     fStrucFunc(kCTEQ5L)
  75 {
  76     // Copy Constructor
  77     pythia.Copy(*this);
  78 }
  79
  80 void AliPythia8::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
  81 {
  82 // Initialise the process to generate
  83     if (!AliPythiaRndm::GetPythiaRandom())
  84       AliPythiaRndm::SetPythiaRandom(GetRandom());
  85
  86     fProcess   = process;
  87     fEcms      = energy;
  88     fStrucFunc = strucfunc;
  89 //...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
  90     ReadString("111:mayDecay  = off");
  91     ReadString("310:mayDecay  = off");
  92     ReadString("3122:mayDecay = off");
  93     ReadString("3112:mayDecay = off");
  94     ReadString("3212:mayDecay = off");
  95     ReadString("3222:mayDecay = off");
  96     ReadString("3312:mayDecay = off");
  97     ReadString("3322:mayDecay = off");
  98     ReadString("3334:mayDecay = off");
  99     // Select structure function
 100
 101     ReadString("PDF:useLHAPDF = on");
 102     ReadString(Form("PDF:LHAPDFset = %s", AliStructFuncType::PDFsetName(fStrucFunc).Data()));
 103
 104     // Particles produced in string fragmentation point directly to either of the two endpoints
 105     // of the string (depending in the side they were generated from).
 106
 107 //    SetMSTU(16,2); // ????
 108
 109 //
 110 // Pythia initialisation for selected processes//
 111 //
 112     switch (process)
 113     {
 114     case kPyOldUEQ2ordered:  //Old underlying events with Q2 ordered QCD processes
 115 //        Multiple interactions on.
 116         ReadString("PartonLevel:MI = on");
 117 // Double Gaussian matter distribution.
 118         ReadString("MultipleInteractions:bProfile = 2");
 119         ReadString("MultipleInteractions:coreFraction = 0.5");
 120         ReadString("MultipleInteractions:coreRadius = 0.4");
 121 //  pT0.
 122         ReadString("MultipleInteractions:pTmin = 2.0");
 123 //  Reference energy for pT0 and energy rescaling pace.
 124         ReadString("MultipleInteractions:ecmRef = 1800.");
 125         ReadString("MultipleInteractions:ecmPow = 0.25");
 126 //  String drawing almost completely minimizes string length.
 127 //      SetPARP(85,0.9);
 128 //      SetPARP(86,0.95);
 129 // ISR and FSR activity.
 130 // Q^2 scale of the hard scattering
 131         ReadString("SigmaProcess:factorMultFac = 4.");
 132 // Lambda_FSR scale.
 133 //      SetPARJ(81, 0.29);
 134         break;
 135     case kPyOldUEQ2ordered2:
 136 // Old underlying events with Q2 ordered QCD processes
 137 // Multiple interactions on.
 138         ReadString("PartonLevel:MI = on");
 139 // Double Gaussian matter distribution.
 140         ReadString("MultipleInteractions:bProfile = 2");
 141         ReadString("MultipleInteractions:coreFraction = 0.5");
 142         ReadString("MultipleInteractions:coreRadius = 0.4");
 143 // pT0.
 144         ReadString("MultipleInteractions:pTmin = 2.0");
 145 // Reference energy for pT0 and energy rescaling pace.
 146         ReadString("MultipleInteractions:ecmRef = 1800.");
 147         ReadString("MultipleInteractions:ecmPow = 0.16");
 148 // String drawing almost completely minimizes string length.
 149 //      SetPARP(85,0.9);
 150 //      SetPARP(86,0.95);
 151 // ISR and FSR activity.
 152         ReadString("SigmaProcess:factorMultFac = 4.");
 153 // Lambda_FSR scale.
 154 //      SetPARJ(81,0.29);
 155         break;
 156     case kPyOldPopcorn:
 157 // Old production mechanism: Old Popcorn
 158         ReadString("HardQCD:all = on");
 159 //      SetMSTJ(12,3);
 160 // (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
 161 //      SetMSTP(88,2);
 162 // (D=1)see can be used to form  baryons (BARYON JUNCTION)
 163 //      SetMSTJ(1,1);
 164         AtlasTuning();
 165         break;
 166     case kPyCharm:
 167         ReadString("HardQCD:gg2ccbar = on");
 168         ReadString("HardQCD:qqbar2ccbar = on");
 169 //  heavy quark masses
 170         ReadString("ParticleData:mcRun = 1.2");
 171 //
 172 //    primordial pT
 173         ReadString("Beams:primordialKT = on");
 174         ReadString("Beams:primordialKTsoft = 0.");
 175         ReadString("Beams:primordialKThard = 1.");
 176         ReadString("Beams:halfScaleForKT = 0.");
 177         ReadString("Beams:halfMassForKT = 0.");
 178         break;
 179     case kPyBeauty:
 180         ReadString("HardQCD:gg2bbbar = on");
 181         ReadString("HardQCD:qqbar2bbbar = on");
 182         ReadString("ParticleData:mbRun = 4.75");
 183         break;
 184     case kPyJpsi:
 185 // gg->J/Psi g
 186         ReadString("Charmonium:gg2QQbar[3S1(1)]g = on");
 187         break;
 188     case kPyJpsiChi:
 189         ReadString("Charmonium:all = on");
 190         break;
 191     case kPyCharmUnforced:
 192 // gq->qg
 193         ReadString("HardQCD:gq2qg = on");
 194 // gg->qq
 195         ReadString("HardQCD:gg2qq = on");
 196 // gg->gg
 197         ReadString("HardQCD:gg2gg = on");
 198         break;
 199     case kPyBeautyUnforced:
 200 // gq->qg
 201         ReadString("HardQCD:gq2qg = on");
 202 // gg->qq
 203         ReadString("HardQCD:gg2qq = on");
 204 // gg->gg
 205         ReadString("HardQCD:gg2gg = on");
 206         break;
 207     case kPyMb:
 208 // Minimum Bias pp-Collisions
 209 //
 210 //
 211 //      select Pythia min. bias model
 212 // single diffraction AB-->XB
 213         ReadString("SoftQCD:minBias = on");
 214         ReadString("SoftQCD:singleDiffractive = on");
 215         ReadString("SoftQCD:doubleDiffractive = on");
 216         AtlasTuning();
 217         break;
 218     case kPyMbDefault:
 219 // Minimum Bias pp-Collisions
 220 //
 221 //
 222 //      select Pythia min. bias model
 223         ReadString("SoftQCD:minBias = on");
 224         ReadString("SoftQCD:singleDiffractive = on");
 225         ReadString("SoftQCD:doubleDiffractive = on");
 226         break;
 227     case kPyLhwgMb:
 228 // Les Houches Working Group 05 Minimum Bias pp-Collisions: hep-ph/0604120
 229 //  -> Pythia 6.3 or above is needed
 230 //
 231         ReadString("SoftQCD:minBias = on");
 232         ReadString("SoftQCD:singleDiffractive = on");
 233         ReadString("SoftQCD:doubleDiffractive = on");
 234         ReadString(Form("PDF:LHAPDFset = %s", AliStructFuncType::PDFsetName(kCTEQ6ll).Data()));
 235
 236 //      SetMSTP(68,1);
 237 //      SetMSTP(70,2);
 238 //      ReadString("PartonLevel:MI = on");
 239 // Double Gaussian matter distribution.
 240         ReadString("MultipleInteractions:bProfile = 2");
 241         ReadString("MultipleInteractions:coreFraction = 0.5");
 242         ReadString("MultipleInteractions:coreRadius = 0.5");
 243         ReadString("MultipleInteractions:expPow = 0.16");
 244         ReadString("MultipleInteractions:pTmin = 2.3");
 245 //      SetMSTP(88,1);
 246 //      SetPARP(85,0.9);           // Regulates gluon prod. mechanism
 247         break;
 248     case kPyMbNonDiffr:
 249 // Minimum Bias pp-Collisions
 250 //
 251 //
 252 //      select Pythia min. bias model
 253         ReadString("SoftQCD:minBias = on");
 254         AtlasTuning();
 255         break;
 256     case kPyMbMSEL1:
 257         ConfigHeavyFlavor();
 258 // Intrinsic <kT^2>
 259         ReadString("Beams:primordialKT = on");
 260         ReadString("Beams:primordialKTsoft = 0.");
 261         ReadString("Beams:primordialKThard = 1.");
 262         ReadString("Beams:halfScaleForKT = 0.");
 263         ReadString("Beams:halfMassForKT = 0.");
 264 // Set Q-quark mass
 265         ReadString("ParticleData:mcRun = 1.20");
 266         ReadString("ParticleData:mbRun = 4.78");
 267 // Atlas Tuning
 268         AtlasTuning();
 269         break;
 270     case kPyJets:
 271 //
 272 //  QCD Jets
 273 //
 274         ReadString("HardQCD:all = on");
 275 //
 276 //  Pythia Tune A (CDF)
 277 //
 278         ReadString("PartonLevel:MI = on");
 279         ReadString("MultipleInteractions:pTmin = 2.0");
 280         ReadString("MultipleInteractions:pT0Ref = 2.8");
 281         ReadString("MultipleInteractions:ecmRef = 1800.");
 282         ReadString("MultipleInteractions:expPow = 0.25");
 283         ReadString("MultipleInteractions:bProfile = 2");
 284         ReadString("MultipleInteractions:coreFraction = 0.16");
 285         ReadString("MultipleInteractions:coreRadius = 0.4");
 286         ReadString("SigmaProcess:factorMultFac = 2.5");
 287 //      SetPARP(85,0.90) ;         // Regulates gluon prod. mechanism
 288 //      SetPARP(86,0.95);          // Regulates gluon prod. mechanism
 289        break;
 290     case kPyDirectGamma:
 291         ReadString("PromptPhoton:all = on");
 292         break;
 293     case kPyCharmPbPbMNR:
 294     case kPyD0PbPbMNR:
 295     case kPyDPlusPbPbMNR:
 296     case kPyDPlusStrangePbPbMNR:
 297       // Tuning of Pythia parameters aimed to get a resonable agreement
 298       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 299       // c-cbar single inclusive and double differential distributions.
 300       // This parameter settings are meant to work with Pb-Pb collisions
 301       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 302       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 303       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
 304         ConfigHeavyFlavor();
 305       // Intrinsic <kT>
 306         ReadString("Beams:primordialKT = on");
 307         ReadString("Beams:primordialKTsoft = 0.");
 308         ReadString("Beams:primordialKThard = 1.304");
 309         ReadString("Beams:halfScaleForKT = 0.");
 310         ReadString("Beams:halfMassForKT = 0.");
 311       // Set c-quark mass
 312         ReadString("ParticleData:mcRun = 1.20");
 313         break;
 314     case kPyCharmpPbMNR:
 315     case kPyD0pPbMNR:
 316     case kPyDPluspPbMNR:
 317     case kPyDPlusStrangepPbMNR:
 318       // Tuning of Pythia parameters aimed to get a resonable agreement
 319       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 320       // c-cbar single inclusive and double differential distributions.
 321       // This parameter settings are meant to work with p-Pb collisions
 322       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 323       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 324       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
 325         ConfigHeavyFlavor();
 326       // Intrinsic <kT>
 327         ReadString("Beams:primordialKT = on");
 328         ReadString("Beams:primordialKTsoft = 0.");
 329         ReadString("Beams:primordialKThard = 1.16");
 330         ReadString("Beams:halfScaleForKT = 0.");
 331         ReadString("Beams:halfMassForKT = 0.");
 332       // Set c-quark mass
 333         ReadString("ParticleData:mcRun = 1.20");
 334       break;
 335     case kPyCharmppMNR:
 336     case kPyD0ppMNR:
 337     case kPyDPlusppMNR:
 338     case kPyDPlusStrangeppMNR:
 339       // Tuning of Pythia parameters aimed to get a resonable agreement
 340       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 341       // c-cbar single inclusive and double differential distributions.
 342       // This parameter settings are meant to work with pp collisions
 343       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 344       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 345       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
 346         ConfigHeavyFlavor();
 347       // Intrinsic <kT^2>
 348         ReadString("Beams:primordialKT = on");
 349         ReadString("Beams:primordialKTsoft = 0.");
 350         ReadString("Beams:primordialKThard = 1.");
 351         ReadString("Beams:halfScaleForKT = 0.");
 352         ReadString("Beams:halfMassForKT = 0.");
 353       // Set c-quark mass
 354         ReadString("ParticleData:mcRun = 1.20");
 355       break;
 356     case kPyCharmppMNRwmi:
 357       // Tuning of Pythia parameters aimed to get a resonable agreement
 358       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 359       // c-cbar single inclusive and double differential distributions.
 360       // This parameter settings are meant to work with pp collisions
 361       // and with kCTEQ5L PDFs.
 362       // Added multiple interactions according to ATLAS tune settings.
 363       // To get a "reasonable" agreement with MNR results, events have to be
 364       // generated with the minimum ptHard (AliGenPythia::SetPtHard)
 365       // set to 2.76 GeV.
 366       // To get a "perfect" agreement with MNR results, events have to be
 367       // generated in four ptHard bins with the following relative
 368       // normalizations:
 369       // 2.76-3 GeV: 25%
 370       //    3-4 GeV: 40%
 371       //    4-8 GeV: 29%
 372       //     >8 GeV:  6%
 373         ConfigHeavyFlavor();
 374       // Intrinsic <kT^2>
 375         ReadString("Beams:primordialKT = on");
 376         ReadString("Beams:primordialKTsoft = 0.");
 377         ReadString("Beams:primordialKThard = 1.");
 378         ReadString("Beams:halfScaleForKT = 0.");
 379         ReadString("Beams:halfMassForKT = 0.");
 380       // Set c-quark mass
 381         ReadString("ParticleData:mcRun = 1.20");
 382         AtlasTuning();
 383         break;
 384     case kPyBeautyPbPbMNR:
 385       // Tuning of Pythia parameters aimed to get a resonable agreement
 386       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 387       // b-bbar single inclusive and double differential distributions.
 388       // This parameter settings are meant to work with Pb-Pb collisions
 389       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 390       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 391       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
 392         ConfigHeavyFlavor();
 393       // QCD scales
 394         ReadString("SigmaProcess:factorMultFac = 1.");
 395       // Intrinsic <kT>
 396         ReadString("Beams:primordialKT = on");
 397         ReadString("Beams:primordialKTsoft = 0.");
 398         ReadString("Beams:primordialKThard = 2.035");
 399         ReadString("Beams:halfScaleForKT = 0.");
 400         ReadString("Beams:halfMassForKT = 0.");
 401       // Set b-quark mass
 402         ReadString("ParticleData:mbRun = 4.75");
 403       break;
 404     case kPyBeautypPbMNR:
 405       // Tuning of Pythia parameters aimed to get a resonable agreement
 406       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 407       // b-bbar single inclusive and double differential distributions.
 408       // This parameter settings are meant to work with p-Pb collisions
 409       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 410       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 411       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
 412         ConfigHeavyFlavor();
 413       // QCD scales
 414         ReadString("SigmaProcess:factorMultFac = 1.");
 415       // Intrinsic <kT>
 416         ReadString("Beams:primordialKT = on");
 417         ReadString("Beams:primordialKTsoft = 0.");
 418         ReadString("Beams:primordialKThard = 1.6");
 419         ReadString("Beams:halfScaleForKT = 0.");
 420         ReadString("Beams:halfMassForKT = 0.");
 421       // Set b-quark mass
 422         ReadString("ParticleData:mbRun = 4.75");
 423       break;
 424     case kPyBeautyppMNR:
 425       // Tuning of Pythia parameters aimed to get a resonable agreement
 426       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 427       // b-bbar single inclusive and double differential distributions.
 428       // This parameter settings are meant to work with pp collisions
 429       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 430       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 431       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
 432         ConfigHeavyFlavor();
 433         // QCD scales
 434         ReadString("SigmaProcess:factorMultFac = 1.");
 435         // Intrinsic <kT>
 436         ReadString("Beams:primordialKT = on");
 437         ReadString("Beams:primordialKTsoft = 0.");
 438         ReadString("Beams:primordialKThard = 1.0");
 439         ReadString("Beams:halfScaleForKT = 0.");
 440         ReadString("Beams:halfMassForKT = 0.");
 441         // Set b-quark mass
 442         ReadString("ParticleData:mbRun = 4.75");
 443       break;
 444      case kPyBeautyppMNRwmi:
 445       // Tuning of Pythia parameters aimed to get a resonable agreement
 446       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 447       // b-bbar single inclusive and double differential distributions.
 448       // This parameter settings are meant to work with pp collisions
 449       // and with kCTEQ5L PDFs.
 450       // Added multiple interactions according to ATLAS tune settings.
 451       // To get a "reasonable" agreement with MNR results, events have to be
 452       // generated with the minimum ptHard (AliGenPythia::SetPtHard)
 453       // set to 2.76 GeV.
 454       // To get a "perfect" agreement with MNR results, events have to be
 455       // generated in four ptHard bins with the following relative
 456       // normalizations:
 457       // 2.76-4 GeV:  5%
 458       //    4-6 GeV: 31%
 459       //    6-8 GeV: 28%
 460       //     >8 GeV: 36%
 461          ConfigHeavyFlavor();
 462          // QCD scales
 463          ReadString("SigmaProcess:factorMultFac = 1.");
 464          // Intrinsic <kT>
 465         ReadString("Beams:primordialKT = on");
 466         ReadString("Beams:primordialKTsoft = 0.");
 467         ReadString("Beams:primordialKThard = 1.0");
 468         ReadString("Beams:halfScaleForKT = 0.");
 469         ReadString("Beams:halfMassForKT = 0.");
 470         // Set b-quark mass
 471         ReadString("ParticleData:mbRun = 4.75");
 472         AtlasTuning();
 473         break;
 474     case kPyW:
 475         //Inclusive production of W+/-
 476         //f fbar -> W+
 477         ReadString("WeakSingleBoson:ffbar2W = on");
 478         // Initial/final parton shower on (Pythia default)
 479         // With parton showers on we are generating "W inclusive process"
 480         ReadString("PartonLevel:ISR = on");
 481         ReadString("PartonLevel:FSR = on");
 482         break;
 483     case kPyZ:
 484         //Inclusive production of Z
 485         //f fbar -> Z/gamma
 486         ReadString("WeakSingleBoson:ffbar2gmZ = on");
 487         //only Z included, not gamma
 488         ReadString("WeakZ0:gmZmode = 2");
 489         // Initial/final parton shower on (Pythia default)
 490         // With parton showers on we are generating "Z inclusive process"
 491         ReadString("PartonLevel:ISR = on");
 492         ReadString("PartonLevel:FSR = on");
 493         break;
 494     }
 495 //
 496 //  Initialize PYTHIA
 497 //    SetMSTP(41,1);   // all resonance decays switched on
 498     Initialize(2212, 2212, fEcms);
 499 }
 500
 501 void AliPythia8::SetNuclei(Int_t /*a1*/, Int_t /*a2*/)
 502 {
 503 // Treat protons as inside nuclei with mass numbers a1 and a2
 504 //    The MSTP array in the PYPARS common block is used to enable and
 505 //    select the nuclear structure functions.
 506 //    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
 507 //            =1: internal PYTHIA acording to MSTP(51)
 508 //            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
 509 //    If the following mass number both not equal zero, nuclear corrections of the stf are used.
 510 //    MSTP(192) : Mass number of nucleus side 1
 511 //    MSTP(193) : Mass number of nucleus side 2
 512 //    SetMSTP(52,2);
 513 //    SetMSTP(192, a1);
 514 //    SetMSTP(193, a2);
 515 }
 516
 517
 518 AliPythia8* AliPythia8::Instance()
 519 {
 520 // Set random number generator
 521     if (fgAliPythia8) {
 522         return fgAliPythia8;
 523     } else {
 524         fgAliPythia8 = new AliPythia8();
 525         return fgAliPythia8;
 526     }
 527 }
 528
 529 void AliPythia8::PrintParticles()
 530 {
 531 // Print list of particl properties
 532     ReadString("Main:showAllParticleData");
 533 }
 534
 535 void  AliPythia8::ResetDecayTable()
 536 {
 537 //  Set default values for pythia decay switches
 538 //    Int_t i;
 539 //    for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
 540 //    for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
 541 }
 542
 543 void  AliPythia8::SetDecayTable()
 544 {
 545 //  Set default values for pythia decay switches
 546 //
 547 //    Int_t i;
 548 //    for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
 549 //    for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
 550 }
 551
 552 void  AliPythia8::Pyclus(Int_t& njet)
 553 {
 554 //  Call Pythia clustering algorithm
 555 //
 556     Bool_t ok = fClusterJet.analyze(Pythia8()->event, fYScale, fPtScale, fNJetMin, fNJetMax);
 557     njet = 0;
 558     if (ok) njet = fClusterJet.size();
 559 }
 560
 561 void  AliPythia8::Pycell(Int_t& njet)
 562 {
 563 //  Call Pythia jet reconstruction algorithm
 564 //
 565     Bool_t ok = fCellJet.analyze(Pythia8()->event, fMinEtJet, fRJet, fEtSeed);
 566     njet = 0;
 567     if (ok) njet = fCellJet.size();
 568 }
 569
 570 void AliPythia8::GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e)
 571 {
 572     // Get jet number i
 573     Float_t et = fCellJet.eT(i);
 574     px = et * TMath::Cos(fCellJet.phiWeighted(i));
 575     py = et * TMath::Sin(fCellJet.phiWeighted(i));
 576     pz = et * TMath::SinH(fCellJet.etaWeighted(i));
 577     e  = et * TMath::CosH(fCellJet.etaWeighted(i));
 578 }
 579
 580 void AliPythia8::GenerateEvent()
 581 {
 582     // Generate one event
 583     TPythia8::GenerateEvent();
 584 }
 585
 586 void AliPythia8::GenerateMIEvent()
 587 {
 588     // New multiple interaction scenario
 589     AliWarning("Not implemented. No event will be generated");
 590 }
 591
 592 void AliPythia8::PrintStatistics()
 593 {
 594     // End of run statistics
 595     TPythia8::PrintStatistics();
 596 }
 597
 598 void AliPythia8::EventListing()
 599 {
 600     // End of run statistics
 601     TPythia8::EventListing();
 602 }
 603
 604 Int_t AliPythia8::ProcessCode()
 605 {
 606     // Returns the subprocess code for the current event
 607     return Pythia8()->info.codeSub();
 608 }
 609
 610 void AliPythia8::ConfigHeavyFlavor()
 611 {
 612     //
 613     // Default configuration for Heavy Flavor production
 614     //
 615     // All QCD processes
 616     //
 617     ReadString("HardQCD:all = on");
 618
 619     // No multiple interactions
 620     ReadString("PartonLevel:MI = off");
 621     ReadString("MultipleInteractions:pTmin = 0.0");
 622     ReadString("MultipleInteractions:pT0Ref = 0.0");
 623
 624     // Initial/final parton shower on (Pythia default)
 625     ReadString("PartonLevel:ISR = on");
 626     ReadString("PartonLevel:FSR = on");
 627
 628     // 2nd order alpha_s
 629     ReadString("SigmaProcess:alphaSorder = 2");
 630
 631     // QCD scales
 632     ReadString("SigmaProcess:renormScale2 = 2");
 633     ReadString("SigmaProcess:renormMultFac = 1.");
 634 }
 635
 636 void AliPythia8::AtlasTuning()
 637 {
 638     //
 639     // Configuration for the ATLAS tuning
 640     ReadString(Form("PDF:LHAPDFset = %s", AliStructFuncType::PDFsetName(kCTEQ5L).Data()));
 641     ReadString("PartonLevel:MI = on");
 642     ReadString("MultipleInteractions:pTmin = 1.9");
 643     ReadString("MultipleInteractions:pT0Ref = 1.8");
 644     ReadString("MultipleInteractions:ecmRef = 1000.");
 645     ReadString("MultipleInteractions:expPow = 0.16");
 646     ReadString("MultipleInteractions:bProfile = 2");
 647     ReadString("MultipleInteractions:coreFraction = 0.16");
 648     ReadString("MultipleInteractions:coreRadius = 0.5");
 649 //      SetPARP(85,0.33);          // Regulates gluon prod. mechanism
 650 //      SetPARP(86,0.66);          // Regulates gluon prod. mechanism
 651     ReadString("SigmaProcess:factorMultFac = 1.");
 652 }
 653
 654 void AliPythia8::SetPtHardRange(Float_t ptmin, Float_t ptmax)
 655 {
 656     // Set the pt hard range
 657     ReadString(Form("PhaseSpace:pTHatMin = %13.3f", ptmin));
 658     ReadString(Form("PhaseSpace:pTHatMax = %13.3f", ptmax));
 659 }
 660
 661 void AliPythia8::SetYHardRange(Float_t /*ymin*/, Float_t /*ymax*/)
 662 {
 663     // Set the y hard range
 664     printf("YHardRange not implemented in Pythia8 !!!\n");
 665
 666 }
 667
 668
 669 void AliPythia8::SetFragmentation(Int_t flag)
 670 {
 671     // Switch fragmentation on/off
 672     if (flag) {
 673         ReadString("HadronLevel:Hadronize = on");
 674     } else {
 675         ReadString("HadronLevel:Hadronize = off");
 676     }
 677 }
 678
 679 void AliPythia8::SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2)
 680 {
 681 //  initial state radiation
 682     if (flag1) {
 683         ReadString("PartonLevel:ISR = on");
 684     } else {
 685         ReadString("PartonLevel:ISR = off");
 686     }
 687 // final state radiation
 688     if (flag2) {
 689         ReadString("PartonLevel:FSR = on");
 690     } else {
 691         ReadString("PartonLevel:FSR = off");
 692     }
 693 }
 694
 695 void AliPythia8::SetIntrinsicKt(Float_t kt)
 696 {
 697         ReadString("Beams:primordialKT = on");
 698         ReadString("Beams:primordialKTsoft = 0.");
 699         ReadString(Form("Beams:primordialKThard = %13.3f", kt));
 700         ReadString("Beams:halfScaleForKT = 0.");
 701         ReadString("Beams:halfMassForKT = 0.");
 702 }
 703
 704 void AliPythia8::SwitchHFOff()
 705 {
 706     // Switch off heavy flavor
 707     // Maximum number of quark flavours used in pdf
 708     ReadString("PDFinProcess:nQuarkIn = 3");
 709     // Maximum number of flavors that can be used in showers
 710     ReadString("TimeShower:nGluonToQuark = 3");
 711     ReadString("SpaceShower:nQuarkIn = 3");
 712
 713
 714 }
 715
 716 void AliPythia8::SetPycellParameters(Float_t etaMax, Int_t nEta, Int_t nPhi,
 717                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
 718 {
 719 // Set pycell parameters
 720     fCellJet  = Pythia8::CellJet( etaMax, nEta, nPhi, 2, 0, 0., 0., thresh);
 721     fEtSeed   = etseed;
 722     fMinEtJet = minet;
 723     fRJet     = r;
 724 }
 725
 726 void AliPythia8::ModifiedSplitting()
 727 {
 728 //
 729 // We have to see how to implement this in Pythia8 !!!
 730 //
 731     // Modified splitting probability as a model for quenching
 732 //    SetPARJ(200, 0.8);
 733 //    SetMSTJ(41, 1);  // QCD radiation only
 734 //    SetMSTJ(42, 2);  // angular ordering
 735 //    SetMSTJ(44, 2);  // option to run alpha_s
 736 //    SetMSTJ(47, 0);  // No correction back to hard scattering element
 737 //    SetMSTJ(50, 0);  // No coherence in first branching
 738 //    SetPARJ(82, 1.); // Cut off for parton showers
 739 }
 740
 741
 742 void AliPythia8::InitQuenching(Float_t /*cMin*/, Float_t /*cMax*/, Float_t /*k*/, Int_t /*iECMethod*/, Float_t /*zmax*/, Int_t /*ngmax*/)
 743 {
 744     //
 745     //
 746     AliWarning("Not implemented !");
 747 }
 748
 749 void AliPythia8::SwitchHadronisationOff()
 750 {
 751     // Switch off hadronisation
 752     ReadString("HadronLevel:Hadronize = off");
 753 }
 754
 755 void AliPythia8::SwitchHadronisationOn()
 756 {
 757     // Switch on hadronisarion
 758     ReadString("HadronLevel:Hadronize = on");
 759 }
 760
 761
 762 void AliPythia8::GetXandQ(Float_t& x1, Float_t& x2, Float_t& q)
 763 {
 764     // Get x1, x2 and Q for this event
 765
 766     q  = Pythia8()->info.QFac();
 767     x1 = Pythia8()->info.x1();
 768     x2 = Pythia8()->info.x2();
 769
 770 }
 771
 772 Float_t AliPythia8::GetXSection()
 773 {
 774     // Get the total cross-section
 775     return Pythia8()->info.sigmaGen();
 776 }
 777
 778 Float_t AliPythia8::GetPtHard()
 779 {
 780     // Get the pT hard for this event
 781     return Pythia8()->info.pTHat();
 782 }
 783
 784
 785
 786
 787 AliPythia8& AliPythia8::operator=(const  AliPythia8& rhs)
 788 {
 789 // Assignment operator
 790     rhs.Copy(*this);
 791     return *this;
 792 }
 793
 794  void AliPythia8::Copy(TObject&) const
 795 {
 796     //
 797     // Copy
 798     //
 799     Fatal("Copy","Not implemented!\n");
 800 }
 801
 802 //
 803 // To be implemented
 804 //
 805 void AliPythia8::SetNumberOfParticles(Int_t /*i*/)
 806 {
 807     AliWarning("Not implemented");
 808 }
 809
 810 void AliPythia8::EditEventList(Int_t /*i*/)
 811 {
 812     AliWarning("Not implemented");
 813 }
 814
 815 void AliPythia8::Pyquen(Double_t /*a*/, Int_t /*b*/, Double_t /*c*/)
 816 {
 817     AliWarning("Cannot be used with Pythia8");
 818 }
 819
 820 void AliPythia8::HadronizeEvent()
 821 {
 822     // Needs access to HadronLevel ?
 823     AliWarning("Not yet implemented");
 824 }
 825
 826 void AliPythia8::GetQuenchingParameters(Double_t& /*xp*/, Double_t& /*yp*/, Double_t* /*z[4]*/)
 827 {
 828     AliWarning("Not yet implemented");
 829 }
 830
 831 void AliPythia8::LoadEvent(AliStack* /*stack*/, Int_t /*flag*/, Int_t /*reHadr*/)
 832 {
 833     AliWarning("Not yet implemented");
 834 }