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