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