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