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