PYTHIA8/AliPythia8.cxx

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