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