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     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     AliTPythia8(),
  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("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         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("BeamRemnants:primordialKT = on");
 274         ReadString("BeamRemnants:primordialKTsoft = 0.");
 275         ReadString("BeamRemnants:primordialKThard = 1.");
 276         ReadString("BeamRemnants:halfScaleForKT = 0.");
 277         ReadString("BeamRemnants: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("BeamRemnants:primordialKT = on");
 321         ReadString("BeamRemnants:primordialKTsoft = 0.");
 322         ReadString("BeamRemnants:primordialKThard = 1.304");
 323         ReadString("BeamRemnants:halfScaleForKT = 0.");
 324         ReadString("BeamRemnants: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("BeamRemnants:primordialKT = on");
 342         ReadString("BeamRemnants:primordialKTsoft = 0.");
 343         ReadString("BeamRemnants:primordialKThard = 1.16");
 344         ReadString("BeamRemnants:halfScaleForKT = 0.");
 345         ReadString("BeamRemnants: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     case kPyLambdacppMNR:
 354       // Tuning of Pythia parameters aimed to get a resonable agreement
 355       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 356       // c-cbar single inclusive and double differential distributions.
 357       // This parameter settings are meant to work with pp collisions
 358       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 359       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 360       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
 361         ConfigHeavyFlavor();
 362       // Intrinsic <kT^2>
 363         ReadString("BeamRemnants:primordialKT = on");
 364         ReadString("BeamRemnants:primordialKTsoft = 0.");
 365         ReadString("BeamRemnants:primordialKThard = 1.");
 366         ReadString("BeamRemnants:halfScaleForKT = 0.");
 367         ReadString("BeamRemnants:halfMassForKT = 0.");
 368       // Set c-quark mass
 369         ReadString("ParticleData:mcRun = 1.20");
 370       break;
 371     case kPyCharmppMNRwmi:
 372       // Tuning of Pythia parameters aimed to get a resonable agreement
 373       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 374       // c-cbar single inclusive and double differential distributions.
 375       // This parameter settings are meant to work with pp collisions
 376       // and with kCTEQ5L PDFs.
 377       // Added multiple interactions according to ATLAS tune settings.
 378       // To get a "reasonable" agreement with MNR results, events have to be
 379       // generated with the minimum ptHard (AliGenPythia::SetPtHard)
 380       // set to 2.76 GeV.
 381       // To get a "perfect" agreement with MNR results, events have to be
 382       // generated in four ptHard bins with the following relative
 383       // normalizations:
 384       // 2.76-3 GeV: 25%
 385       //    3-4 GeV: 40%
 386       //    4-8 GeV: 29%
 387       //     >8 GeV:  6%
 388         ConfigHeavyFlavor();
 389       // Intrinsic <kT^2>
 390         ReadString("BeamRemnants:primordialKT = on");
 391         ReadString("BeamRemnants:primordialKTsoft = 0.");
 392         ReadString("BeamRemnants:primordialKThard = 1.");
 393         ReadString("BeamRemnants:halfScaleForKT = 0.");
 394         ReadString("BeamRemnants:halfMassForKT = 0.");
 395       // Set c-quark mass
 396         ReadString("ParticleData:mcRun = 1.20");
 397         AtlasTuning();
 398         break;
 399     case kPyBeautyPbPbMNR:
 400       // Tuning of Pythia parameters aimed to get a resonable agreement
 401       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 402       // b-bbar single inclusive and double differential distributions.
 403       // This parameter settings are meant to work with Pb-Pb collisions
 404       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 405       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 406       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
 407         ConfigHeavyFlavor();
 408       // QCD scales
 409         ReadString("SigmaProcess:factorMultFac = 1.");
 410       // Intrinsic <kT>
 411         ReadString("BeamRemnants:primordialKT = on");
 412         ReadString("BeamRemnants:primordialKTsoft = 0.");
 413         ReadString("BeamRemnants:primordialKThard = 2.035");
 414         ReadString("BeamRemnants:halfScaleForKT = 0.");
 415         ReadString("BeamRemnants:halfMassForKT = 0.");
 416       // Set b-quark mass
 417         ReadString("ParticleData:mbRun = 4.75");
 418       break;
 419     case kPyBeautypPbMNR:
 420       // Tuning of Pythia parameters aimed to get a resonable agreement
 421       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 422       // b-bbar single inclusive and double differential distributions.
 423       // This parameter settings are meant to work with p-Pb collisions
 424       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 425       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 426       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
 427         ConfigHeavyFlavor();
 428       // QCD scales
 429         ReadString("SigmaProcess:factorMultFac = 1.");
 430       // Intrinsic <kT>
 431         ReadString("BeamRemnants:primordialKT = on");
 432         ReadString("BeamRemnants:primordialKTsoft = 0.");
 433         ReadString("BeamRemnants:primordialKThard = 1.6");
 434         ReadString("BeamRemnants:halfScaleForKT = 0.");
 435         ReadString("BeamRemnants:halfMassForKT = 0.");
 436       // Set b-quark mass
 437         ReadString("ParticleData:mbRun = 4.75");
 438       break;
 439     case kPyBeautyppMNR:
 440       // Tuning of Pythia parameters aimed to get a resonable agreement
 441       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 442       // b-bbar single inclusive and double differential distributions.
 443       // This parameter settings are meant to work with pp collisions
 444       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 445       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 446       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
 447         ConfigHeavyFlavor();
 448         // QCD scales
 449         ReadString("SigmaProcess:factorMultFac = 1.");
 450         // Intrinsic <kT>
 451         ReadString("BeamRemnants:primordialKT = on");
 452         ReadString("BeamRemnants:primordialKTsoft = 0.");
 453         ReadString("BeamRemnants:primordialKThard = 1.0");
 454         ReadString("BeamRemnants:halfScaleForKT = 0.");
 455         ReadString("BeamRemnants:halfMassForKT = 0.");
 456         // Set b-quark mass
 457         ReadString("ParticleData:mbRun = 4.75");
 458       break;
 459      case kPyBeautyppMNRwmi:
 460       // Tuning of Pythia parameters aimed to get a resonable agreement
 461       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 462       // b-bbar single inclusive and double differential distributions.
 463       // This parameter settings are meant to work with pp collisions
 464       // and with kCTEQ5L PDFs.
 465       // Added multiple interactions according to ATLAS tune settings.
 466       // To get a "reasonable" agreement with MNR results, events have to be
 467       // generated with the minimum ptHard (AliGenPythia::SetPtHard)
 468       // set to 2.76 GeV.
 469       // To get a "perfect" agreement with MNR results, events have to be
 470       // generated in four ptHard bins with the following relative
 471       // normalizations:
 472       // 2.76-4 GeV:  5%
 473       //    4-6 GeV: 31%
 474       //    6-8 GeV: 28%
 475       //     >8 GeV: 36%
 476          ConfigHeavyFlavor();
 477          // QCD scales
 478          ReadString("SigmaProcess:factorMultFac = 1.");
 479          // Intrinsic <kT>
 480         ReadString("BeamRemnants:primordialKT = on");
 481         ReadString("BeamRemnants:primordialKTsoft = 0.");
 482         ReadString("BeamRemnants:primordialKThard = 1.0");
 483         ReadString("BeamRemnants:halfScaleForKT = 0.");
 484         ReadString("BeamRemnants:halfMassForKT = 0.");
 485         // Set b-quark mass
 486         ReadString("ParticleData:mbRun = 4.75");
 487         AtlasTuning();
 488         break;
 489     case kPyW:
 490         //Inclusive production of W+/-
 491         //f fbar -> W+
 492         ReadString("WeakSingleBoson:ffbar2W = on");
 493         // Initial/final parton shower on (Pythia default)
 494         // With parton showers on we are generating "W inclusive process"
 495         ReadString("PartonLevel:ISR = on");
 496         ReadString("PartonLevel:FSR = on");
 497         break;
 498     case kPyZ:
 499         //Inclusive production of Z
 500         //f fbar -> Z/gamma
 501         ReadString("WeakSingleBoson:ffbar2gmZ = on");
 502         //only Z included, not gamma
 503         ReadString("WeakZ0:gmZmode = 2");
 504         // Initial/final parton shower on (Pythia default)
 505         // With parton showers on we are generating "Z inclusive process"
 506         ReadString("PartonLevel:ISR = on");
 507         ReadString("PartonLevel:FSR = on");
 508         break;
 509     }
 510 //
 511 //  Initialize PYTHIA
 512 //    SetMSTP(41,1);   // all resonance decays switched on
 513     Initialize(2212, 2212, fEcms);
 514 }
 515
 516 void AliPythia8::SetNuclei(Int_t /*a1*/, Int_t /*a2*/)
 517 {
 518 // Treat protons as inside nuclei with mass numbers a1 and a2
 519 //    The MSTP array in the PYPARS common block is used to enable and
 520 //    select the nuclear structure functions.
 521 //    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
 522 //            =1: internal PYTHIA acording to MSTP(51)
 523 //            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
 524 //    If the following mass number both not equal zero, nuclear corrections of the stf are used.
 525 //    MSTP(192) : Mass number of nucleus side 1
 526 //    MSTP(193) : Mass number of nucleus side 2
 527 //    SetMSTP(52,2);
 528 //    SetMSTP(192, a1);
 529 //    SetMSTP(193, a2);
 530 }
 531
 532
 533 AliPythia8* AliPythia8::Instance()
 534 {
 535 // Set random number generator
 536     if (fgAliPythia8) {
 537         return fgAliPythia8;
 538     } else {
 539         fgAliPythia8 = new AliPythia8();
 540         return fgAliPythia8;
 541     }
 542 }
 543
 544 void AliPythia8::PrintParticles()
 545 {
 546 // Print list of particl properties
 547     ReadString("Main:showAllParticleData");
 548 }
 549
 550 void  AliPythia8::ResetDecayTable()
 551 {
 552 //  Set default values for pythia decay switches
 553 //    Int_t i;
 554 //    for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
 555 //    for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
 556 }
 557
 558 void  AliPythia8::SetDecayTable()
 559 {
 560 //  Set default values for pythia decay switches
 561 //
 562 //    Int_t i;
 563 //    for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
 564 //    for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
 565 }
 566
 567 void  AliPythia8::Pyclus(Int_t& njet)
 568 {
 569 //  Call Pythia clustering algorithm
 570 //
 571     Bool_t ok = fClusterJet.analyze(Pythia8()->event, fYScale, fPtScale, fNJetMin, fNJetMax);
 572     njet = 0;
 573     if (ok) njet = fClusterJet.size();
 574 }
 575
 576 void  AliPythia8::Pycell(Int_t& njet)
 577 {
 578 //  Call Pythia jet reconstruction algorithm
 579 //
 580     Bool_t ok = fCellJet.analyze(Pythia8()->event, fMinEtJet, fRJet, fEtSeed);
 581     njet = 0;
 582     if (ok) njet = fCellJet.size();
 583 }
 584
 585 void AliPythia8::GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e)
 586 {
 587     // Get jet number i
 588     Float_t et = fCellJet.eT(i);
 589     px = et * TMath::Cos(fCellJet.phiWeighted(i));
 590     py = et * TMath::Sin(fCellJet.phiWeighted(i));
 591     pz = et * TMath::SinH(fCellJet.etaWeighted(i));
 592     e  = et * TMath::CosH(fCellJet.etaWeighted(i));
 593 }
 594
 595 void AliPythia8::GenerateEvent()
 596 {
 597     // Generate one event
 598     AliTPythia8::GenerateEvent();
 599 }
 600
 601 void AliPythia8::GenerateMIEvent()
 602 {
 603     // New multiple interaction scenario
 604     AliWarning("Not implemented. No event will be generated");
 605 }
 606
 607 void AliPythia8::PrintStatistics()
 608 {
 609     // End of run statistics
 610     AliTPythia8::PrintStatistics();
 611 }
 612
 613 void AliPythia8::EventListing()
 614 {
 615     // End of run statistics
 616     AliTPythia8::EventListing();
 617 }
 618
 619 Int_t AliPythia8::ProcessCode()
 620 {
 621     // Returns the subprocess code for the current event
 622     printf("Process %5d %5d \n", Pythia8()->info.code(),  Pythia8()->info.codeSub());
 623
 624     return Pythia8()->info.codeSub();
 625 }
 626
 627 void AliPythia8::ConfigHeavyFlavor()
 628 {
 629     //
 630     // Default configuration for Heavy Flavor production
 631     //
 632     // All QCD processes
 633     //
 634     ReadString("HardQCD:all = on");
 635
 636     // No multiple interactions
 637     ReadString("PartonLevel:MI = off");
 638     ReadString("MultipleInteractions:pTmin = 0.0");
 639     ReadString("MultipleInteractions:pT0Ref = 0.0");
 640
 641     // Initial/final parton shower on (Pythia default)
 642     ReadString("PartonLevel:ISR = on");
 643     ReadString("PartonLevel:FSR = on");
 644
 645     // 2nd order alpha_s
 646     ReadString("SigmaProcess:alphaSorder = 2");
 647
 648     // QCD scales
 649     ReadString("SigmaProcess:renormScale2 = 2");
 650     ReadString("SigmaProcess:renormMultFac = 1.");
 651 }
 652
 653 void AliPythia8::AtlasTuning()
 654 {
 655     //
 656     // Configuration for the ATLAS tuning
 657     ReadString(Form("PDF:LHAPDFset = %s", AliStructFuncType::PDFsetName(kCTEQ5L).Data()));
 658     ReadString("PartonLevel:MI = on");
 659     ReadString("MultipleInteractions:pTmin = 1.9");
 660     ReadString("MultipleInteractions:pT0Ref = 1.8");
 661     ReadString("MultipleInteractions:ecmRef = 1000.");
 662     ReadString("MultipleInteractions:expPow = 0.16");
 663     ReadString("MultipleInteractions:bProfile = 2");
 664     ReadString("MultipleInteractions:coreFraction = 0.16");
 665     ReadString("MultipleInteractions:coreRadius = 0.5");
 666 //      SetPARP(85,0.33);          // Regulates gluon prod. mechanism
 667 //      SetPARP(86,0.66);          // Regulates gluon prod. mechanism
 668     ReadString("SigmaProcess:factorMultFac = 1.");
 669 }
 670
 671 void AliPythia8::SetPtHardRange(Float_t ptmin, Float_t ptmax)
 672 {
 673     // Set the pt hard range
 674     ReadString(Form("PhaseSpace:pTHatMin = %13.3f", ptmin));
 675     ReadString(Form("PhaseSpace:pTHatMax = %13.3f", ptmax));
 676 }
 677
 678 void AliPythia8::SetYHardRange(Float_t /*ymin*/, Float_t /*ymax*/)
 679 {
 680     // Set the y hard range
 681     printf("YHardRange not implemented in Pythia8 !!!\n");
 682
 683 }
 684
 685
 686 void AliPythia8::SetFragmentation(Int_t flag)
 687 {
 688     // Switch fragmentation on/off
 689     if (flag) {
 690         ReadString("HadronLevel:Hadronize = on");
 691     } else {
 692         ReadString("HadronLevel:Hadronize = off");
 693     }
 694 }
 695
 696 void AliPythia8::SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2)
 697 {
 698 //  initial state radiation
 699     if (flag1) {
 700         ReadString("PartonLevel:ISR = on");
 701     } else {
 702         ReadString("PartonLevel:ISR = off");
 703     }
 704 // final state radiation
 705     if (flag2) {
 706         ReadString("PartonLevel:FSR = on");
 707     } else {
 708         ReadString("PartonLevel:FSR = off");
 709     }
 710 }
 711
 712 void AliPythia8::SetIntrinsicKt(Float_t kt)
 713 {
 714         ReadString("BeamRemnants:primordialKT = on");
 715         ReadString("BeamRemnants:primordialKTsoft = 0.");
 716         ReadString(Form("BeamRemnants:primordialKThard = %13.3f", kt));
 717         ReadString("BeamRemnants:halfScaleForKT = 0.");
 718         ReadString("BeamRemnants:halfMassForKT = 0.");
 719 }
 720
 721 void AliPythia8::SwitchHFOff()
 722 {
 723     // Switch off heavy flavor
 724     // Maximum number of quark flavours used in pdf
 725     ReadString("PDFinProcess:nQuarkIn = 3");
 726     // Maximum number of flavors that can be used in showers
 727     ReadString("TimeShower:nGluonToQuark = 3");
 728     ReadString("SpaceShower:nQuarkIn = 3");
 729
 730
 731 }
 732
 733 void AliPythia8::SetPycellParameters(Float_t etaMax, Int_t nEta, Int_t nPhi,
 734                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
 735 {
 736 // Set pycell parameters
 737     fCellJet  = Pythia8::CellJet( etaMax, nEta, nPhi, 2, 0, 0., 0., thresh);
 738     fEtSeed   = etseed;
 739     fMinEtJet = minet;
 740     fRJet     = r;
 741 }
 742
 743 void AliPythia8::ModifiedSplitting()
 744 {
 745 //
 746 // We have to see how to implement this in Pythia8 !!!
 747 //
 748     // Modified splitting probability as a model for quenching
 749 //    SetPARJ(200, 0.8);
 750 //    SetMSTJ(41, 1);  // QCD radiation only
 751 //    SetMSTJ(42, 2);  // angular ordering
 752 //    SetMSTJ(44, 2);  // option to run alpha_s
 753 //    SetMSTJ(47, 0);  // No correction back to hard scattering element
 754 //    SetMSTJ(50, 0);  // No coherence in first branching
 755 //    SetPARJ(82, 1.); // Cut off for parton showers
 756 }
 757
 758
 759 void AliPythia8::InitQuenching(Float_t /*cMin*/, Float_t /*cMax*/, Float_t /*k*/, Int_t /*iECMethod*/, Float_t /*zmax*/, Int_t /*ngmax*/)
 760 {
 761     //
 762     //
 763     AliWarning("Not implemented !");
 764 }
 765
 766 void AliPythia8::SwitchHadronisationOff()
 767 {
 768     // Switch off hadronisation
 769     ReadString("HadronLevel:Hadronize = off");
 770 }
 771
 772 void AliPythia8::SwitchHadronisationOn()
 773 {
 774     // Switch on hadronisarion
 775     ReadString("HadronLevel:Hadronize = on");
 776 }
 777
 778
 779 void AliPythia8::GetXandQ(Float_t& x1, Float_t& x2, Float_t& q)
 780 {
 781     // Get x1, x2 and Q for this event
 782
 783     q  = Pythia8()->info.QFac();
 784     x1 = Pythia8()->info.x1();
 785     x2 = Pythia8()->info.x2();
 786
 787 }
 788
 789 Float_t AliPythia8::GetXSection()
 790 {
 791     // Get the total cross-section
 792     return Pythia8()->info.sigmaGen();
 793 }
 794
 795 Float_t AliPythia8::GetPtHard()
 796 {
 797     // Get the pT hard for this event
 798     return Pythia8()->info.pTHat();
 799 }
 800
 801
 802
 803
 804 AliPythia8& AliPythia8::operator=(const  AliPythia8& rhs)
 805 {
 806 // Assignment operator
 807     rhs.Copy(*this);
 808     return *this;
 809 }
 810
 811  void AliPythia8::Copy(TObject&) const
 812 {
 813     //
 814     // Copy
 815     //
 816     Fatal("Copy","Not implemented!\n");
 817 }
 818
 819 //
 820 // To be implemented
 821 //
 822 void AliPythia8::SetNumberOfParticles(Int_t /*i*/)
 823 {
 824     AliWarning("Not implemented");
 825 }
 826
 827 void AliPythia8::EditEventList(Int_t /*i*/)
 828 {
 829     AliWarning("Not implemented");
 830 }
 831
 832 void AliPythia8::Pyquen(Double_t /*a*/, Int_t /*b*/, Double_t /*c*/)
 833 {
 834     AliWarning("Cannot be used with Pythia8");
 835 }
 836
 837 void AliPythia8::HadronizeEvent()
 838 {
 839     // Needs access to HadronLevel ?
 840     AliWarning("Not yet implemented");
 841 }
 842
 843 void AliPythia8::GetQuenchingParameters(Double_t& /*xp*/, Double_t& /*yp*/, Double_t* /*z[4]*/)
 844 {
 845     AliWarning("Not yet implemented");
 846 }
 847
 848 void AliPythia8::LoadEvent(AliStack* /*stack*/, Int_t /*flag*/, Int_t /*reHadr*/)
 849 {
 850     AliWarning("Not yet implemented");
 851 }