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::SetSeed(UInt_t seed)
 547 {
 548   //
 549   // set seed in PYTHIA8
 550   // NB. 900000000 is the maximum seed (0 is not allowed)
 551   //
 552   ReadString("Random:setSeed = on");
 553   ReadString(Form("Random:seed = %d", (seed % 900000000) + 1));
 554 }
 555
 556 void AliPythia8::SetNuclei(Int_t /*a1*/, Int_t /*a2*/)
 557 {
 558 // Treat protons as inside nuclei with mass numbers a1 and a2
 559 //    The MSTP array in the PYPARS common block is used to enable and
 560 //    select the nuclear structure functions.
 561 //    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
 562 //            =1: internal PYTHIA acording to MSTP(51)
 563 //            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
 564 //    If the following mass number both not equal zero, nuclear corrections of the stf are used.
 565 //    MSTP(192) : Mass number of nucleus side 1
 566 //    MSTP(193) : Mass number of nucleus side 2
 567 //    SetMSTP(52,2);
 568 //    SetMSTP(192, a1);
 569 //    SetMSTP(193, a2);
 570 }
 571
 572
 573 AliPythia8* AliPythia8::Instance()
 574 {
 575 // Set random number generator
 576     if (fgAliPythia8) {
 577         return fgAliPythia8;
 578     } else {
 579         fgAliPythia8 = new AliPythia8();
 580         return fgAliPythia8;
 581     }
 582 }
 583
 584 void AliPythia8::PrintParticles()
 585 {
 586 // Print list of particl properties
 587     ReadString("Main:showAllParticleData");
 588 }
 589
 590 void  AliPythia8::ResetDecayTable()
 591 {
 592 //  Set default values for pythia decay switches
 593 //    Int_t i;
 594 //    for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
 595 //    for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
 596 }
 597
 598 void  AliPythia8::SetDecayTable()
 599 {
 600 //  Set default values for pythia decay switches
 601 //
 602 //    Int_t i;
 603 //    for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
 604 //    for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
 605 }
 606
 607 void  AliPythia8::Pyclus(Int_t& njet)
 608 {
 609 //  Call Pythia clustering algorithm
 610 //
 611     Bool_t ok = fClusterJet.analyze(Pythia8()->event, fYScale, fPtScale, fNJetMin, fNJetMax);
 612     njet = 0;
 613     if (ok) njet = fClusterJet.size();
 614 }
 615
 616 void  AliPythia8::Pycell(Int_t& njet)
 617 {
 618 //  Call Pythia jet reconstruction algorithm
 619 //
 620     Bool_t ok = fCellJet.analyze(Pythia8()->event, fMinEtJet, fRJet, fEtSeed);
 621     njet = 0;
 622     if (ok) njet = fCellJet.size();
 623 }
 624
 625 void AliPythia8::GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e)
 626 {
 627     // Get jet number i
 628     Float_t et = fCellJet.eT(i);
 629     px = et * TMath::Cos(fCellJet.phiWeighted(i));
 630     py = et * TMath::Sin(fCellJet.phiWeighted(i));
 631     pz = et * TMath::SinH(fCellJet.etaWeighted(i));
 632     e  = et * TMath::CosH(fCellJet.etaWeighted(i));
 633 }
 634
 635 void AliPythia8::GenerateEvent()
 636 {
 637     // Generate one event
 638     AliTPythia8::GenerateEvent();
 639 }
 640
 641 void AliPythia8::GenerateMIEvent()
 642 {
 643     // New multiple interaction scenario
 644     AliWarning("Not implemented. No event will be generated");
 645 }
 646
 647 void AliPythia8::PrintStatistics()
 648 {
 649     // End of run statistics
 650     AliTPythia8::PrintStatistics();
 651 }
 652
 653 void AliPythia8::EventListing()
 654 {
 655     // End of run statistics
 656     AliTPythia8::EventListing();
 657 }
 658
 659 Int_t AliPythia8::ProcessCode()
 660 {
 661     // Returns the subprocess code for the current event
 662     return Pythia8()->info.code();
 663 }
 664
 665 void AliPythia8::ConfigHeavyFlavor()
 666 {
 667     //
 668     // Default configuration for Heavy Flavor production
 669     //
 670     // All QCD processes
 671     //
 672     ReadString("HardQCD:all = on");
 673
 674     // No multiple interactions
 675     ReadString("PartonLevel:MI = off");
 676     ReadString("MultipleInteractions:pTmin = 0.0");
 677     ReadString("MultipleInteractions:pT0Ref = 0.0");
 678
 679     // Initial/final parton shower on (Pythia default)
 680     ReadString("PartonLevel:ISR = on");
 681     ReadString("PartonLevel:FSR = on");
 682
 683     // 2nd order alpha_s
 684     ReadString("SigmaProcess:alphaSorder = 2");
 685
 686     // QCD scales
 687     ReadString("SigmaProcess:renormScale2 = 2");
 688     ReadString("SigmaProcess:renormMultFac = 1.");
 689 }
 690
 691 void AliPythia8::AtlasTuning()
 692 {
 693     //
 694     // Configuration for the ATLAS tuning
 695     ReadString(Form("PDF:LHAPDFset = %s", AliStructFuncType::PDFsetName(kCTEQ5L).Data()));
 696     ReadString("PartonLevel:MI = on");
 697     ReadString("MultipleInteractions:pTmin = 1.9");
 698     ReadString("MultipleInteractions:pT0Ref = 1.8");
 699     ReadString("MultipleInteractions:ecmRef = 1000.");
 700     ReadString("MultipleInteractions:expPow = 0.16");
 701     ReadString("MultipleInteractions:bProfile = 2");
 702     ReadString("MultipleInteractions:coreFraction = 0.16");
 703     ReadString("MultipleInteractions:coreRadius = 0.5");
 704 //      SetPARP(85,0.33);          // Regulates gluon prod. mechanism
 705 //      SetPARP(86,0.66);          // Regulates gluon prod. mechanism
 706     ReadString("SigmaProcess:factorMultFac = 1.");
 707
 708 }
 709
 710 void AliPythia8::SetPtHardRange(Float_t ptmin, Float_t ptmax)
 711 {
 712     // Set the pt hard range
 713     ReadString(Form("PhaseSpace:pTHatMin = %13.3f", ptmin));
 714     ReadString(Form("PhaseSpace:pTHatMax = %13.3f", ptmax));
 715 }
 716
 717 void AliPythia8::SetYHardRange(Float_t /*ymin*/, Float_t /*ymax*/)
 718 {
 719     // Set the y hard range
 720     printf("YHardRange not implemented in Pythia8 !!!\n");
 721
 722 }
 723
 724
 725 void AliPythia8::SetFragmentation(Int_t flag)
 726 {
 727     // Switch fragmentation on/off
 728     if (flag) {
 729         ReadString("HadronLevel:Hadronize = on");
 730     } else {
 731         ReadString("HadronLevel:Hadronize = off");
 732     }
 733 }
 734
 735 void AliPythia8::SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2)
 736 {
 737 //  initial state radiation
 738     if (flag1) {
 739         ReadString("PartonLevel:ISR = on");
 740     } else {
 741         ReadString("PartonLevel:ISR = off");
 742     }
 743 // final state radiation
 744     if (flag2) {
 745         ReadString("PartonLevel:FSR = on");
 746     } else {
 747         ReadString("PartonLevel:FSR = off");
 748     }
 749 }
 750
 751 void AliPythia8::SetIntrinsicKt(Float_t kt)
 752 {
 753 // Set the intrinsic kt
 754         ReadString("BeamRemnants:primordialKT = on");
 755         ReadString("BeamRemnants:primordialKTsoft = 0.");
 756         ReadString(Form("BeamRemnants:primordialKThard = %13.3f", kt));
 757         ReadString("BeamRemnants:halfScaleForKT = 0.");
 758         ReadString("BeamRemnants:halfMassForKT = 0.");
 759 }
 760
 761 void AliPythia8::SwitchHFOff()
 762 {
 763     // Switch off heavy flavor
 764     // Maximum number of quark flavours used in pdf
 765     ReadString("PDFinProcess:nQuarkIn = 3");
 766     // Maximum number of flavors that can be used in showers
 767     ReadString("TimeShower:nGluonToQuark = 3");
 768     ReadString("SpaceShower:nQuarkIn = 3");
 769
 770
 771 }
 772
 773 void AliPythia8::SetPycellParameters(Float_t etaMax, Int_t nEta, Int_t nPhi,
 774                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
 775 {
 776 // Set pycell parameters
 777     fCellJet  = Pythia8::CellJet( etaMax, nEta, nPhi, 2, 0, 0., 0., thresh);
 778     fEtSeed   = etseed;
 779     fMinEtJet = minet;
 780     fRJet     = r;
 781 }
 782
 783 void AliPythia8::ModifiedSplitting()
 784 {
 785 //
 786 // We have to see how to implement this in Pythia8 !!!
 787 //
 788     // Modified splitting probability as a model for quenching
 789 //    SetPARJ(200, 0.8);
 790 //    SetMSTJ(41, 1);  // QCD radiation only
 791 //    SetMSTJ(42, 2);  // angular ordering
 792 //    SetMSTJ(44, 2);  // option to run alpha_s
 793 //    SetMSTJ(47, 0);  // No correction back to hard scattering element
 794 //    SetMSTJ(50, 0);  // No coherence in first branching
 795 //    SetPARJ(82, 1.); // Cut off for parton showers
 796 }
 797
 798
 799 void AliPythia8::InitQuenching(Float_t /*cMin*/, Float_t /*cMax*/, Float_t /*k*/, Int_t /*iECMethod*/, Float_t /*zmax*/, Int_t /*ngmax*/)
 800 {
 801     //
 802     //
 803     AliWarning("Not implemented !");
 804 }
 805
 806 void AliPythia8::SwitchHadronisationOff()
 807 {
 808     // Switch off hadronisation
 809     ReadString("HadronLevel:Hadronize = off");
 810 }
 811
 812 void AliPythia8::SwitchHadronisationOn()
 813 {
 814     // Switch on hadronisarion
 815     ReadString("HadronLevel:Hadronize = on");
 816 }
 817
 818
 819 void AliPythia8::GetXandQ(Float_t& x1, Float_t& x2, Float_t& q)
 820 {
 821     // Get x1, x2 and Q for this event
 822
 823     q  = Pythia8()->info.QFac();
 824     x1 = Pythia8()->info.x1();
 825     x2 = Pythia8()->info.x2();
 826
 827 }
 828
 829 Float_t AliPythia8::GetXSection()
 830 {
 831     // Get the total cross-section
 832     return Pythia8()->info.sigmaGen();
 833 }
 834
 835 Float_t AliPythia8::GetPtHard()
 836 {
 837     // Get the pT hard for this event
 838     return Pythia8()->info.pTHat();
 839 }
 840
 841
 842
 843
 844 AliPythia8& AliPythia8::operator=(const  AliPythia8& rhs)
 845 {
 846 // Assignment operator
 847     rhs.Copy(*this);
 848     return *this;
 849 }
 850
 851  void AliPythia8::Copy(TObject&) const
 852 {
 853     //
 854     // Copy
 855     //
 856     Fatal("Copy","Not implemented!\n");
 857 }
 858
 859 //
 860 // To be implemented
 861 //
 862 void AliPythia8::SetNumberOfParticles(Int_t /*i*/)
 863 {
 864     AliWarning("Not implemented");
 865 }
 866
 867 void AliPythia8::EditEventList(Int_t /*i*/)
 868 {
 869     AliWarning("Not implemented");
 870 }
 871
 872 void AliPythia8::Pyquen(Double_t /*a*/, Int_t /*b*/, Double_t /*c*/)
 873 {
 874     AliWarning("Cannot be used with Pythia8");
 875 }
 876
 877 void AliPythia8::HadronizeEvent()
 878 {
 879     // Needs access to HadronLevel ?
 880     AliWarning("Not yet implemented");
 881 }
 882
 883 void AliPythia8::GetQuenchingParameters(Double_t& /*xp*/, Double_t& /*yp*/, Double_t* /*z[4]*/)
 884 {
 885     AliWarning("Not yet implemented");
 886 }
 887
 888 void AliPythia8::LoadEvent(AliStack* /*stack*/, Int_t /*flag*/, Int_t /*reHadr*/)
 889 {
 890     AliWarning("Not yet implemented");
 891 }