--- /dev/null
+
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id$ */
+#include <TString.h>
+#include <TVector3.h>
+#include <TMath.h>
+
+#include "AliPythia8.h"
+#include "AliLog.h"
+#include "AliStack.h"
+#include "AliPythiaRndm.h"
+
+
+ClassImp(AliPythia8)
+
+// Particles produced in string fragmentation point directly to either of the two endpoints
+// of the string (depending in the side they were generated from).
+// SetMSTU(16,2); // ????
+// String drawing almost completely minimizes string length.
+// Probability that an additional interaction gives two gluons
+// ... with color connection to nearest neighbours
+// SetPARP(85,0.9);
+// ... as closed gluon loop
+// SetPARP(86,0.95);
+// Lambda_FSR scale.
+// SetPARJ(81, 0.29);
+// Baryon production model
+// SetMSTJ(12,3);
+// String fragmentation
+// SetMSTJ(1,1);
+// sea quarks can be used for baryon formation
+// SetMSTP(88,2);
+// choice of max. virtuality for ISR
+// SetMSTP(68,1);
+// regularisation scheme of ISR
+// SetMSTP(70,2);
+// all resonance decays switched on
+// SetMSTP(41,1);
+AliPythia8* AliPythia8::fgAliPythia8=NULL;
+
+AliPythia8::AliPythia8():
+ TPythia8(),
+ AliPythiaBase(),
+ fProcess(kPyMb),
+ fEcms(0.),
+ fStrucFunc(kCTEQ5L),
+ fEtSeed(0.),
+ fMinEtJet(0.),
+ fRJet(0.),
+ fYScale(0.),
+ fPtScale(0.),
+ fNJetMin(0),
+ fNJetMax(0)
+{
+// Default Constructor
+//
+// Set random number
+ if (!AliPythiaRndm::GetPythiaRandom())
+ AliPythiaRndm::SetPythiaRandom(GetRandom());
+}
+
+AliPythia8::AliPythia8(const AliPythia8& pythia):
+ TPythia8(),
+ AliPythiaBase(),
+ fProcess(kPyMb),
+ fEcms(0.),
+ fStrucFunc(kCTEQ5L),
+ fEtSeed(0.),
+ fMinEtJet(0.),
+ fRJet(0.),
+ fYScale(0.),
+ fPtScale(0.),
+ fNJetMin(0),
+ fNJetMax(0)
+{
+ // Copy Constructor
+ pythia.Copy(*this);
+}
+
+void AliPythia8::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
+{
+// Initialise the process to generate
+ if (!AliPythiaRndm::GetPythiaRandom())
+ AliPythiaRndm::SetPythiaRandom(GetRandom());
+
+ fProcess = process;
+ fEcms = energy;
+ fStrucFunc = strucfunc;
+//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
+ ReadString("111:mayDecay = off");
+ ReadString("310:mayDecay = off");
+ ReadString("3122:mayDecay = off");
+ ReadString("3112:mayDecay = off");
+ ReadString("3212:mayDecay = off");
+ ReadString("3222:mayDecay = off");
+ ReadString("3312:mayDecay = off");
+ ReadString("3322:mayDecay = off");
+ ReadString("3334:mayDecay = off");
+ // Select structure function
+
+ ReadString("PDF:useLHAPDF = on");
+ ReadString(Form("PDF:LHAPDFset = %s", AliStructFuncType::PDFsetName(fStrucFunc).Data()));
+
+ // Particles produced in string fragmentation point directly to either of the two endpoints
+ // of the string (depending in the side they were generated from).
+
+// SetMSTU(16,2); // ????
+
+//
+// Pythia initialisation for selected processes//
+//
+ switch (process)
+ {
+ case kPyOldUEQ2ordered: //Old underlying events with Q2 ordered QCD processes
+// Multiple interactions on.
+ ReadString("PartonLevel:MI = on");
+// Double Gaussian matter distribution.
+ ReadString("MultipleInteractions:bProfile = 2");
+ ReadString("MultipleInteractions:coreFraction = 0.5");
+ ReadString("MultipleInteractions:coreRadius = 0.4");
+// pT0.
+ ReadString("MultipleInteractions:pTmin = 2.0");
+// Reference energy for pT0 and energy rescaling pace.
+ ReadString("MultipleInteractions:ecmRef = 1800.");
+ ReadString("MultipleInteractions:ecmPow = 0.25");
+// String drawing almost completely minimizes string length.
+// SetPARP(85,0.9);
+// SetPARP(86,0.95);
+// ISR and FSR activity.
+// Q^2 scale of the hard scattering
+ ReadString("SigmaProcess:factorMultFac = 4.");
+// Lambda_FSR scale.
+// SetPARJ(81, 0.29);
+ break;
+ case kPyOldUEQ2ordered2:
+// Old underlying events with Q2 ordered QCD processes
+// Multiple interactions on.
+ ReadString("PartonLevel:MI = on");
+// Double Gaussian matter distribution.
+ ReadString("MultipleInteractions:bProfile = 2");
+ ReadString("MultipleInteractions:coreFraction = 0.5");
+ ReadString("MultipleInteractions:coreRadius = 0.4");
+// pT0.
+ ReadString("MultipleInteractions:pTmin = 2.0");
+// Reference energy for pT0 and energy rescaling pace.
+ ReadString("MultipleInteractions:ecmRef = 1800.");
+ ReadString("MultipleInteractions:ecmPow = 0.16");
+// String drawing almost completely minimizes string length.
+// SetPARP(85,0.9);
+// SetPARP(86,0.95);
+// ISR and FSR activity.
+ ReadString("SigmaProcess:factorMultFac = 4.");
+// Lambda_FSR scale.
+// SetPARJ(81,0.29);
+ break;
+ case kPyOldPopcorn:
+// Old production mechanism: Old Popcorn
+ ReadString("HardQCD:all = on");
+// SetMSTJ(12,3);
+// (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
+// SetMSTP(88,2);
+// (D=1)see can be used to form baryons (BARYON JUNCTION)
+// SetMSTJ(1,1);
+ AtlasTuning();
+ break;
+ case kPyCharm:
+ ReadString("HardQCD:gg2ccbar = on");
+ ReadString("HardQCD:qqbar2ccbar = on");
+// heavy quark masses
+ ReadString("ParticleData:mcRun = 1.2");
+//
+// primordial pT
+ ReadString("Beams:primordialKT = on");
+ ReadString("Beams:primordialKTsoft = 0.");
+ ReadString("Beams:primordialKThard = 1.");
+ ReadString("Beams:halfScaleForKT = 0.");
+ ReadString("Beams:halfMassForKT = 0.");
+ break;
+ case kPyBeauty:
+ ReadString("HardQCD:gg2bbbar = on");
+ ReadString("HardQCD:qqbar2bbbar = on");
+ ReadString("ParticleData:mbRun = 4.75");
+ break;
+ case kPyJpsi:
+// gg->J/Psi g
+ ReadString("Charmonium:gg2QQbar[3S1(1)]g = on");
+ break;
+ case kPyJpsiChi:
+ ReadString("Charmonium:all = on");
+ break;
+ case kPyCharmUnforced:
+// gq->qg
+ ReadString("HardQCD:gq2qg = on");
+// gg->qq
+ ReadString("HardQCD:gg2qq = on");
+// gg->gg
+ ReadString("HardQCD:gg2gg = on");
+ break;
+ case kPyBeautyUnforced:
+// gq->qg
+ ReadString("HardQCD:gq2qg = on");
+// gg->qq
+ ReadString("HardQCD:gg2qq = on");
+// gg->gg
+ ReadString("HardQCD:gg2gg = on");
+ break;
+ case kPyMb:
+// Minimum Bias pp-Collisions
+//
+//
+// select Pythia min. bias model
+// single diffraction AB-->XB
+ ReadString("SoftQCD:minBias = on");
+ ReadString("SoftQCD:singleDiffractive = on");
+ ReadString("SoftQCD:doubleDiffractive = on");
+ AtlasTuning();
+ break;
+ case kPyMbDefault:
+// Minimum Bias pp-Collisions
+//
+//
+// select Pythia min. bias model
+ ReadString("SoftQCD:minBias = on");
+ ReadString("SoftQCD:singleDiffractive = on");
+ ReadString("SoftQCD:doubleDiffractive = on");
+ break;
+ case kPyLhwgMb:
+// Les Houches Working Group 05 Minimum Bias pp-Collisions: hep-ph/0604120
+// -> Pythia 6.3 or above is needed
+//
+ ReadString("SoftQCD:minBias = on");
+ ReadString("SoftQCD:singleDiffractive = on");
+ ReadString("SoftQCD:doubleDiffractive = on");
+ ReadString(Form("PDF:LHAPDFset = %s", AliStructFuncType::PDFsetName(kCTEQ6ll).Data()));
+
+// SetMSTP(68,1);
+// SetMSTP(70,2);
+// ReadString("PartonLevel:MI = on");
+// Double Gaussian matter distribution.
+ ReadString("MultipleInteractions:bProfile = 2");
+ ReadString("MultipleInteractions:coreFraction = 0.5");
+ ReadString("MultipleInteractions:coreRadius = 0.5");
+ ReadString("MultipleInteractions:expPow = 0.16");
+ ReadString("MultipleInteractions:pTmin = 2.3");
+// SetMSTP(88,1);
+// SetPARP(85,0.9); // Regulates gluon prod. mechanism
+ break;
+ case kPyMbNonDiffr:
+// Minimum Bias pp-Collisions
+//
+//
+// select Pythia min. bias model
+ ReadString("SoftQCD:minBias = on");
+ AtlasTuning();
+ break;
+ case kPyMbMSEL1:
+ ConfigHeavyFlavor();
+// Intrinsic <kT^2>
+ ReadString("Beams:primordialKT = on");
+ ReadString("Beams:primordialKTsoft = 0.");
+ ReadString("Beams:primordialKThard = 1.");
+ ReadString("Beams:halfScaleForKT = 0.");
+ ReadString("Beams:halfMassForKT = 0.");
+// Set Q-quark mass
+ ReadString("ParticleData:mcRun = 1.20");
+ ReadString("ParticleData:mbRun = 4.78");
+// Atlas Tuning
+ AtlasTuning();
+ break;
+ case kPyJets:
+//
+// QCD Jets
+//
+ ReadString("HardQCD:all = on");
+//
+// Pythia Tune A (CDF)
+//
+ ReadString("PartonLevel:MI = on");
+ ReadString("MultipleInteractions:pTmin = 2.0");
+ ReadString("MultipleInteractions:pT0Ref = 2.8");
+ ReadString("MultipleInteractions:ecmRef = 1800.");
+ ReadString("MultipleInteractions:expPow = 0.25");
+ ReadString("MultipleInteractions:bProfile = 2");
+ ReadString("MultipleInteractions:coreFraction = 0.16");
+ ReadString("MultipleInteractions:coreRadius = 0.4");
+ ReadString("SigmaProcess:factorMultFac = 2.5");
+// SetPARP(85,0.90) ; // Regulates gluon prod. mechanism
+// SetPARP(86,0.95); // Regulates gluon prod. mechanism
+ break;
+ case kPyDirectGamma:
+ ReadString("PromptPhoton:all = on");
+ break;
+ case kPyCharmPbPbMNR:
+ case kPyD0PbPbMNR:
+ case kPyDPlusPbPbMNR:
+ case kPyDPlusStrangePbPbMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // c-cbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with Pb-Pb collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
+ ConfigHeavyFlavor();
+ // Intrinsic <kT>
+ ReadString("Beams:primordialKT = on");
+ ReadString("Beams:primordialKTsoft = 0.");
+ ReadString("Beams:primordialKThard = 1.304");
+ ReadString("Beams:halfScaleForKT = 0.");
+ ReadString("Beams:halfMassForKT = 0.");
+ // Set c-quark mass
+ ReadString("ParticleData:mcRun = 1.20");
+ break;
+ case kPyCharmpPbMNR:
+ case kPyD0pPbMNR:
+ case kPyDPluspPbMNR:
+ case kPyDPlusStrangepPbMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // c-cbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with p-Pb collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
+ ConfigHeavyFlavor();
+ // Intrinsic <kT>
+ ReadString("Beams:primordialKT = on");
+ ReadString("Beams:primordialKTsoft = 0.");
+ ReadString("Beams:primordialKThard = 1.16");
+ ReadString("Beams:halfScaleForKT = 0.");
+ ReadString("Beams:halfMassForKT = 0.");
+ // Set c-quark mass
+ ReadString("ParticleData:mcRun = 1.20");
+ break;
+ case kPyCharmppMNR:
+ case kPyD0ppMNR:
+ case kPyDPlusppMNR:
+ case kPyDPlusStrangeppMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // c-cbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with pp collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
+ ConfigHeavyFlavor();
+ // Intrinsic <kT^2>
+ ReadString("Beams:primordialKT = on");
+ ReadString("Beams:primordialKTsoft = 0.");
+ ReadString("Beams:primordialKThard = 1.");
+ ReadString("Beams:halfScaleForKT = 0.");
+ ReadString("Beams:halfMassForKT = 0.");
+ // Set c-quark mass
+ ReadString("ParticleData:mcRun = 1.20");
+ break;
+ case kPyCharmppMNRwmi:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // c-cbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with pp collisions
+ // and with kCTEQ5L PDFs.
+ // Added multiple interactions according to ATLAS tune settings.
+ // To get a "reasonable" agreement with MNR results, events have to be
+ // generated with the minimum ptHard (AliGenPythia::SetPtHard)
+ // set to 2.76 GeV.
+ // To get a "perfect" agreement with MNR results, events have to be
+ // generated in four ptHard bins with the following relative
+ // normalizations:
+ // 2.76-3 GeV: 25%
+ // 3-4 GeV: 40%
+ // 4-8 GeV: 29%
+ // >8 GeV: 6%
+ ConfigHeavyFlavor();
+ // Intrinsic <kT^2>
+ ReadString("Beams:primordialKT = on");
+ ReadString("Beams:primordialKTsoft = 0.");
+ ReadString("Beams:primordialKThard = 1.");
+ ReadString("Beams:halfScaleForKT = 0.");
+ ReadString("Beams:halfMassForKT = 0.");
+ // Set c-quark mass
+ ReadString("ParticleData:mcRun = 1.20");
+ AtlasTuning();
+ break;
+ case kPyBeautyPbPbMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // b-bbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with Pb-Pb collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
+ ConfigHeavyFlavor();
+ // QCD scales
+ ReadString("SigmaProcess:factorMultFac = 1.");
+ // Intrinsic <kT>
+ ReadString("Beams:primordialKT = on");
+ ReadString("Beams:primordialKTsoft = 0.");
+ ReadString("Beams:primordialKThard = 2.035");
+ ReadString("Beams:halfScaleForKT = 0.");
+ ReadString("Beams:halfMassForKT = 0.");
+ // Set b-quark mass
+ ReadString("ParticleData:mbRun = 4.75");
+ break;
+ case kPyBeautypPbMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // b-bbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with p-Pb collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
+ ConfigHeavyFlavor();
+ // QCD scales
+ ReadString("SigmaProcess:factorMultFac = 1.");
+ // Intrinsic <kT>
+ ReadString("Beams:primordialKT = on");
+ ReadString("Beams:primordialKTsoft = 0.");
+ ReadString("Beams:primordialKThard = 1.6");
+ ReadString("Beams:halfScaleForKT = 0.");
+ ReadString("Beams:halfMassForKT = 0.");
+ // Set b-quark mass
+ ReadString("ParticleData:mbRun = 4.75");
+ break;
+ case kPyBeautyppMNR:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // b-bbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with pp collisions
+ // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+ // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+ // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
+ ConfigHeavyFlavor();
+ // QCD scales
+ ReadString("SigmaProcess:factorMultFac = 1.");
+ // Intrinsic <kT>
+ ReadString("Beams:primordialKT = on");
+ ReadString("Beams:primordialKTsoft = 0.");
+ ReadString("Beams:primordialKThard = 1.0");
+ ReadString("Beams:halfScaleForKT = 0.");
+ ReadString("Beams:halfMassForKT = 0.");
+ // Set b-quark mass
+ ReadString("ParticleData:mbRun = 4.75");
+ break;
+ case kPyBeautyppMNRwmi:
+ // Tuning of Pythia parameters aimed to get a resonable agreement
+ // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+ // b-bbar single inclusive and double differential distributions.
+ // This parameter settings are meant to work with pp collisions
+ // and with kCTEQ5L PDFs.
+ // Added multiple interactions according to ATLAS tune settings.
+ // To get a "reasonable" agreement with MNR results, events have to be
+ // generated with the minimum ptHard (AliGenPythia::SetPtHard)
+ // set to 2.76 GeV.
+ // To get a "perfect" agreement with MNR results, events have to be
+ // generated in four ptHard bins with the following relative
+ // normalizations:
+ // 2.76-4 GeV: 5%
+ // 4-6 GeV: 31%
+ // 6-8 GeV: 28%
+ // >8 GeV: 36%
+ ConfigHeavyFlavor();
+ // QCD scales
+ ReadString("SigmaProcess:factorMultFac = 1.");
+ // Intrinsic <kT>
+ ReadString("Beams:primordialKT = on");
+ ReadString("Beams:primordialKTsoft = 0.");
+ ReadString("Beams:primordialKThard = 1.0");
+ ReadString("Beams:halfScaleForKT = 0.");
+ ReadString("Beams:halfMassForKT = 0.");
+ // Set b-quark mass
+ ReadString("ParticleData:mbRun = 4.75");
+ AtlasTuning();
+ break;
+ case kPyW:
+ //Inclusive production of W+/-
+ //f fbar -> W+
+ ReadString("WeakSingleBoson:ffbar2W = on");
+ // Initial/final parton shower on (Pythia default)
+ // With parton showers on we are generating "W inclusive process"
+ ReadString("PartonLevel:ISR = on");
+ ReadString("PartonLevel:FSR = on");
+ break;
+ case kPyZ:
+ //Inclusive production of Z
+ //f fbar -> Z/gamma
+ ReadString("WeakSingleBoson:ffbar2gmZ = on");
+ //only Z included, not gamma
+ ReadString("WeakZ0:gmZmode = 2");
+ // Initial/final parton shower on (Pythia default)
+ // With parton showers on we are generating "Z inclusive process"
+ ReadString("PartonLevel:ISR = on");
+ ReadString("PartonLevel:FSR = on");
+ break;
+ }
+//
+// Initialize PYTHIA
+// SetMSTP(41,1); // all resonance decays switched on
+ Initialize(2212, 2212, fEcms);
+}
+
+void AliPythia8::SetNuclei(Int_t /*a1*/, Int_t /*a2*/)
+{
+// Treat protons as inside nuclei with mass numbers a1 and a2
+// The MSTP array in the PYPARS common block is used to enable and
+// select the nuclear structure functions.
+// MSTP(52) : (D=1) choice of proton and nuclear structure-function library
+// =1: internal PYTHIA acording to MSTP(51)
+// =2: PDFLIB proton s.f., with MSTP(51) = 1000xNGROUP+NSET
+// If the following mass number both not equal zero, nuclear corrections of the stf are used.
+// MSTP(192) : Mass number of nucleus side 1
+// MSTP(193) : Mass number of nucleus side 2
+// SetMSTP(52,2);
+// SetMSTP(192, a1);
+// SetMSTP(193, a2);
+}
+
+
+AliPythia8* AliPythia8::Instance()
+{
+// Set random number generator
+ if (fgAliPythia8) {
+ return fgAliPythia8;
+ } else {
+ fgAliPythia8 = new AliPythia8();
+ return fgAliPythia8;
+ }
+}
+
+void AliPythia8::PrintParticles()
+{
+// Print list of particl properties
+ ReadString("Main:showAllParticleData");
+}
+
+void AliPythia8::ResetDecayTable()
+{
+// Set default values for pythia decay switches
+// Int_t i;
+// for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
+// for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
+}
+
+void AliPythia8::SetDecayTable()
+{
+// Set default values for pythia decay switches
+//
+// Int_t i;
+// for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
+// for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
+}
+
+void AliPythia8::Pyclus(Int_t& njet)
+{
+// Call Pythia clustering algorithm
+//
+ Bool_t ok = fClusterJet.analyze(Pythia8()->event, fYScale, fPtScale, fNJetMin, fNJetMax);
+ njet = 0;
+ if (ok) njet = fClusterJet.size();
+}
+
+void AliPythia8::Pycell(Int_t& njet)
+{
+// Call Pythia jet reconstruction algorithm
+//
+ Bool_t ok = fCellJet.analyze(Pythia8()->event, fMinEtJet, fRJet, fEtSeed);
+ njet = 0;
+ if (ok) njet = fCellJet.size();
+}
+
+void AliPythia8::GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e)
+{
+ // Get jet number i
+ Float_t et = fCellJet.eT(i);
+ px = et * TMath::Cos(fCellJet.phiWeighted(i));
+ py = et * TMath::Sin(fCellJet.phiWeighted(i));
+ pz = et * TMath::SinH(fCellJet.etaWeighted(i));
+ e = et * TMath::CosH(fCellJet.etaWeighted(i));
+}
+
+void AliPythia8::GenerateEvent()
+{
+ // Generate one event
+ TPythia8::GenerateEvent();
+}
+
+void AliPythia8::GenerateMIEvent()
+{
+ // New multiple interaction scenario
+ AliWarning("Not implemented. No event will be generated");
+}
+
+void AliPythia8::PrintStatistics()
+{
+ // End of run statistics
+ TPythia8::PrintStatistics();
+}
+
+void AliPythia8::EventListing()
+{
+ // End of run statistics
+ TPythia8::EventListing();
+}
+
+Int_t AliPythia8::ProcessCode()
+{
+ // Returns the subprocess code for the current event
+ return Pythia8()->info.codeSub();
+}
+
+void AliPythia8::ConfigHeavyFlavor()
+{
+ //
+ // Default configuration for Heavy Flavor production
+ //
+ // All QCD processes
+ //
+ ReadString("HardQCD:all = on");
+
+ // No multiple interactions
+ ReadString("PartonLevel:MI = off");
+ ReadString("MultipleInteractions:pTmin = 0.0");
+ ReadString("MultipleInteractions:pT0Ref = 0.0");
+
+ // Initial/final parton shower on (Pythia default)
+ ReadString("PartonLevel:ISR = on");
+ ReadString("PartonLevel:FSR = on");
+
+ // 2nd order alpha_s
+ ReadString("SigmaProcess:alphaSorder = 2");
+
+ // QCD scales
+ ReadString("SigmaProcess:renormScale2 = 2");
+ ReadString("SigmaProcess:renormMultFac = 1.");
+}
+
+void AliPythia8::AtlasTuning()
+{
+ //
+ // Configuration for the ATLAS tuning
+ ReadString(Form("PDF:LHAPDFset = %s", AliStructFuncType::PDFsetName(kCTEQ5L).Data()));
+ ReadString("PartonLevel:MI = on");
+ ReadString("MultipleInteractions:pTmin = 1.9");
+ ReadString("MultipleInteractions:pT0Ref = 1.8");
+ ReadString("MultipleInteractions:ecmRef = 1000.");
+ ReadString("MultipleInteractions:expPow = 0.16");
+ ReadString("MultipleInteractions:bProfile = 2");
+ ReadString("MultipleInteractions:coreFraction = 0.16");
+ ReadString("MultipleInteractions:coreRadius = 0.5");
+// SetPARP(85,0.33); // Regulates gluon prod. mechanism
+// SetPARP(86,0.66); // Regulates gluon prod. mechanism
+ ReadString("SigmaProcess:factorMultFac = 1.");
+}
+
+void AliPythia8::SetPtHardRange(Float_t ptmin, Float_t ptmax)
+{
+ // Set the pt hard range
+ ReadString(Form("PhaseSpace:pTHatMin = %13.3f", ptmin));
+ ReadString(Form("PhaseSpace:pTHatMax = %13.3f", ptmax));
+}
+
+void AliPythia8::SetYHardRange(Float_t /*ymin*/, Float_t /*ymax*/)
+{
+ // Set the y hard range
+ printf("YHardRange not implemented in Pythia8 !!!\n");
+
+}
+
+
+void AliPythia8::SetFragmentation(Int_t flag)
+{
+ // Switch fragmentation on/off
+ if (flag) {
+ ReadString("HadronLevel:Hadronize = on");
+ } else {
+ ReadString("HadronLevel:Hadronize = off");
+ }
+}
+
+void AliPythia8::SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2)
+{
+// initial state radiation
+ if (flag1) {
+ ReadString("PartonLevel:ISR = on");
+ } else {
+ ReadString("PartonLevel:ISR = off");
+ }
+// final state radiation
+ if (flag2) {
+ ReadString("PartonLevel:FSR = on");
+ } else {
+ ReadString("PartonLevel:FSR = off");
+ }
+}
+
+void AliPythia8::SetIntrinsicKt(Float_t kt)
+{
+ ReadString("Beams:primordialKT = on");
+ ReadString("Beams:primordialKTsoft = 0.");
+ ReadString(Form("Beams:primordialKThard = %13.3f", kt));
+ ReadString("Beams:halfScaleForKT = 0.");
+ ReadString("Beams:halfMassForKT = 0.");
+}
+
+void AliPythia8::SwitchHFOff()
+{
+ // Switch off heavy flavor
+ // Maximum number of quark flavours used in pdf
+ ReadString("PDFinProcess:nQuarkIn = 3");
+ // Maximum number of flavors that can be used in showers
+ ReadString("TimeShower:nGluonToQuark = 3");
+ ReadString("SpaceShower:nQuarkIn = 3");
+
+
+}
+
+void AliPythia8::SetPycellParameters(Float_t etaMax, Int_t nEta, Int_t nPhi,
+ Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
+{
+// Set pycell parameters
+ fCellJet = Pythia8::CellJet( etaMax, nEta, nPhi, 2, 0, 0., 0., thresh);
+ fEtSeed = etseed;
+ fMinEtJet = minet;
+ fRJet = r;
+}
+
+void AliPythia8::ModifiedSplitting()
+{
+//
+// We have to see how to implement this in Pythia8 !!!
+//
+ // Modified splitting probability as a model for quenching
+// SetPARJ(200, 0.8);
+// SetMSTJ(41, 1); // QCD radiation only
+// SetMSTJ(42, 2); // angular ordering
+// SetMSTJ(44, 2); // option to run alpha_s
+// SetMSTJ(47, 0); // No correction back to hard scattering element
+// SetMSTJ(50, 0); // No coherence in first branching
+// SetPARJ(82, 1.); // Cut off for parton showers
+}
+
+
+void AliPythia8::InitQuenching(Float_t /*cMin*/, Float_t /*cMax*/, Float_t /*k*/, Int_t /*iECMethod*/, Float_t /*zmax*/, Int_t /*ngmax*/)
+{
+ //
+ //
+ AliWarning("Not implemented !");
+}
+
+void AliPythia8::SwitchHadronisationOff()
+{
+ // Switch off hadronisation
+ ReadString("HadronLevel:Hadronize = off");
+}
+
+void AliPythia8::SwitchHadronisationOn()
+{
+ // Switch on hadronisarion
+ ReadString("HadronLevel:Hadronize = on");
+}
+
+
+void AliPythia8::GetXandQ(Float_t& x1, Float_t& x2, Float_t& q)
+{
+ // Get x1, x2 and Q for this event
+
+ q = Pythia8()->info.QFac();
+ x1 = Pythia8()->info.x1();
+ x2 = Pythia8()->info.x2();
+
+}
+
+Float_t AliPythia8::GetXSection()
+{
+ // Get the total cross-section
+ return Pythia8()->info.sigmaGen();
+}
+
+Float_t AliPythia8::GetPtHard()
+{
+ // Get the pT hard for this event
+ return Pythia8()->info.pTHat();
+}
+
+
+
+
+AliPythia8& AliPythia8::operator=(const AliPythia8& rhs)
+{
+// Assignment operator
+ rhs.Copy(*this);
+ return *this;
+}
+
+ void AliPythia8::Copy(TObject&) const
+{
+ //
+ // Copy
+ //
+ Fatal("Copy","Not implemented!\n");
+}
+
+//
+// To be implemented
+//
+void AliPythia8::SetNumberOfParticles(Int_t /*i*/)
+{
+ AliWarning("Not implemented");
+}
+
+void AliPythia8::EditEventList(Int_t /*i*/)
+{
+ AliWarning("Not implemented");
+}
+
+void AliPythia8::Pyquen(Double_t /*a*/, Int_t /*b*/, Double_t /*c*/)
+{
+ AliWarning("Cannot be used with Pythia8");
+}
+
+void AliPythia8::HadronizeEvent()
+{
+ // Needs access to HadronLevel ?
+ AliWarning("Not yet implemented");
+}
+
+void AliPythia8::GetQuenchingParameters(Double_t& /*xp*/, Double_t& /*yp*/, Double_t* /*z[4]*/)
+{
+ AliWarning("Not yet implemented");
+}
+
+void AliPythia8::LoadEvent(AliStack* /*stack*/, Int_t /*flag*/, Int_t /*reHadr*/)
+{
+ AliWarning("Not yet implemented");
+}
--- /dev/null
+#ifndef ALIPYTHIA8_H
+#define ALIPYTHIA8_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id: AliPythia.h,v 1.22 2007/10/09 08:43:24 morsch Exp $ */
+
+#include <TPythia8.h>
+#include "Analysis.h"
+#include "AliPythiaBase.h"
+
+class AliStack;
+class AliPythia8 : public TPythia8, public AliPythiaBase
+{
+
+ public:
+ AliPythia8();
+ AliPythia8(const AliPythia8& pythia);
+ virtual ~AliPythia8() {;}
+ virtual Int_t Version() {return (8);}
+ // convert to compressed code and print result (for debugging only)
+ virtual Int_t CheckedLuComp(Int_t /*kf*/) {return -1;}
+ // Pythia initialisation for selected processes
+ virtual void ProcInit (Process_t process, Float_t energy, StrucFunc_t strucfunc);
+ virtual void GenerateEvent();
+ virtual void GenerateMIEvent();
+ virtual void HadronizeEvent();
+ virtual Int_t GetNumberOfParticles() {return GetN();}
+ virtual void SetNumberOfParticles(Int_t i);
+ virtual void EditEventList(Int_t i);
+ virtual void PrintStatistics();
+ virtual void EventListing();
+ virtual Int_t GetParticles(TClonesArray *particles) {return ImportParticles(particles, "All");}
+ // Treat protons as inside nuclei
+ virtual void SetNuclei(Int_t a1, Int_t a2);
+ // Print particle properties
+ virtual void PrintParticles();
+ // Reset the decay table
+ virtual void ResetDecayTable();
+ //
+ // Common Physics Configuration
+ virtual void SetPtHardRange(Float_t ptmin, Float_t ptmax);
+ virtual void SetYHardRange(Float_t ymin, Float_t ymax);
+ virtual void SetFragmentation(Int_t flag);
+ virtual void SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2);
+ virtual void SetIntrinsicKt(Float_t kt);
+ virtual void SwitchHFOff();
+ virtual void SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
+ Float_t thresh, Float_t etseed, Float_t minet, Float_t r);
+ virtual void GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e);
+ virtual void ModifiedSplitting();
+ virtual void InitQuenching(Float_t bmin, Float_t bmax, Float_t k,
+ Int_t iECMethod, Float_t zmax = 0.97,
+ Int_t ngmax = 30);
+ virtual void SwitchHadronisationOff();
+ virtual void SwitchHadronisationOn();
+ //
+ // Common Getters
+ virtual void GetXandQ(Float_t& x1, Float_t& x2, Float_t& q);
+ virtual Float_t GetXSection();
+ virtual Float_t GetPtHard();
+ //
+ //
+ virtual void SetDecayTable();
+ virtual void Pystat(Int_t /*i*/){;}
+ virtual void Pylist(Int_t /*i*/){;}
+ virtual Int_t ProcessCode();
+ virtual void Pycell(Int_t& nclus);
+ virtual void Pyclus(Int_t& nclus);
+ virtual void Pyshow(Int_t /*ip1*/, Int_t /*ip2*/, Double_t /*qmax*/) {;}
+ virtual void Pyrobo(Int_t /*imi*/, Int_t /*ima*/, Double_t /*the*/,
+ Double_t /*phi*/, Double_t /*bex*/, Double_t /*bey*/, Double_t /*bez*/) {;}
+ static AliPythia8* Instance();
+ virtual void Quench() {;}
+ // Assignment Operator
+ AliPythia8 & operator=(const AliPythia8 & rhs);
+ void Copy(TObject&) const;
+ //
+ // Not yet implemented
+ //
+ virtual void Pyquen(Double_t, Int_t, Double_t);
+ virtual void Pyexec() {;}
+ virtual void GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4]);
+ virtual void LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr);
+ virtual void ConfigHeavyFlavor();
+ virtual void AtlasTuning();
+ protected:
+ Process_t fProcess; // Process type
+ Float_t fEcms; // Centre of mass energy
+ StrucFunc_t fStrucFunc; // Structure function
+ Int_t fDefMDCY[501]; // ! Default decay switches per particle
+ Int_t fDefMDME[2001]; // ! Default decay switches per mode
+ Double_t fZQuench[4]; // ! Quenching fractions for this even
+ Pythia8::CellJet fCellJet; // ! Cell jet object
+ Float_t fEtSeed; // ! Et seed for cell jets
+ Float_t fMinEtJet; // ! Min jet et
+ Float_t fRJet; // ! Radius for cell jets
+ Pythia8::ClusterJet fClusterJet; // ! Cluster jet object
+ Float_t fYScale; // ! cut-off joining scale
+ Float_t fPtScale; // ! cut-off joining scale
+ Int_t fNJetMin; // ! min. number of jets
+ Int_t fNJetMax; // ! max. number of jets
+ static AliPythia8* fgAliPythia8; // Pointer to single instance
+
+ ClassDef(AliPythia8, 1) //ALICE UI to PYTHIA8
+};
+
+#endif
+
+
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff