#include <TDatabasePDG.h>
#include <TParticle.h>
#include <TPDGCode.h>
+#include <TObjArray.h>
#include <TSystem.h>
#include <TTree.h>
#include "AliConst.h"
#include "AliDecayerPythia.h"
#include "AliGenPythia.h"
+#include "AliFastGlauber.h"
#include "AliHeader.h"
#include "AliGenPythiaEventHeader.h"
#include "AliPythia.h"
AliGenPythia::AliGenPythia():
AliGenMC(),
fProcess(kPyCharm),
+ fItune(-1),
fStrucFunc(kCTEQ5L),
fKineBias(0.),
fTrials(0),
fGinit(1),
fGfinal(1),
fHadronisation(1),
+ fPatchOmegaDalitz(0),
fNpartons(0),
fReadFromFile(0),
fQuench(0),
+ fQhat(0.),
+ fLength(0.),
+ fImpact(0.),
fPtKick(1.),
fFullEvent(kTRUE),
fDecayer(new AliDecayerPythia()),
fTriggerEta(0.9),
fTriggerMultiplicity(0),
fTriggerMultiplicityEta(0),
+ fTriggerMultiplicityPtMin(0),
fCountMode(kCountAll),
fHeader(0),
fRL(0),
- fFileName(0),
+ fkFileName(0),
fFragPhotonInCalo(kFALSE),
fPi0InCalo(kFALSE) ,
fPhotonInCalo(kFALSE),
+ fEleInEMCAL(kFALSE),
fCheckEMCAL(kFALSE),
fCheckPHOS(kFALSE),
fCheckPHOSeta(kFALSE),
fFragPhotonOrPi0MinPt(0),
fPhotonMinPt(0),
+ fElectronMinPt(0),
fPHOSMinPhi(219.),
fPHOSMaxPhi(321.),
fPHOSEta(0.13),
{
// Default Constructor
fEnergyCMS = 5500.;
- SetNuclei(0,0);
if (!AliPythiaRndm::GetPythiaRandom())
AliPythiaRndm::SetPythiaRandom(GetRandom());
}
AliGenPythia::AliGenPythia(Int_t npart)
:AliGenMC(npart),
fProcess(kPyCharm),
+ fItune(-1),
fStrucFunc(kCTEQ5L),
fKineBias(0.),
fTrials(0),
fGinit(kTRUE),
fGfinal(kTRUE),
fHadronisation(kTRUE),
+ fPatchOmegaDalitz(0),
fNpartons(0),
fReadFromFile(kFALSE),
fQuench(kFALSE),
+ fQhat(0.),
+ fLength(0.),
+ fImpact(0.),
fPtKick(1.),
fFullEvent(kTRUE),
fDecayer(new AliDecayerPythia()),
fTriggerEta(0.9),
fTriggerMultiplicity(0),
fTriggerMultiplicityEta(0),
+ fTriggerMultiplicityPtMin(0),
fCountMode(kCountAll),
fHeader(0),
fRL(0),
- fFileName(0),
+ fkFileName(0),
fFragPhotonInCalo(kFALSE),
fPi0InCalo(kFALSE) ,
fPhotonInCalo(kFALSE),
+ fEleInEMCAL(kFALSE),
fCheckEMCAL(kFALSE),
fCheckPHOS(kFALSE),
fCheckPHOSeta(kFALSE),
fFragPhotonOrPi0MinPt(0),
fPhotonMinPt(0),
+ fElectronMinPt(0),
fPHOSMinPhi(219.),
fPHOSMaxPhi(321.),
fPHOSEta(0.13),
// Set random number generator
if (!AliPythiaRndm::GetPythiaRandom())
AliPythiaRndm::SetPythiaRandom(GetRandom());
- SetNuclei(0,0);
}
AliGenPythia::~AliGenPythia()
if (fReadFromFile) {
- fRL = AliRunLoader::Open(fFileName, "Partons");
+ fRL = AliRunLoader::Open(fkFileName, "Partons");
fRL->LoadKinematics();
fRL->LoadHeader();
} else {
fRL = 0x0;
}
//
- fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
+ fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
// Forward Paramters to the AliPythia object
fDecayer->SetForceDecay(fForceDecay);
// Switch off Heavy Flavors on request
fParentSelect[0] = 431;
fFlavorSelect = 4;
break;
+ case kPyLambdacppMNR:
+ fParentSelect[0] = 4122;
+ fFlavorSelect = 4;
+ break;
case kPyBeauty:
+ case kPyBeautyJets:
case kPyBeautyPbPbMNR:
case kPyBeautypPbMNR:
case kPyBeautyppMNR:
fParentSelect[0] = 443;
break;
case kPyMbDefault:
+ case kPyMbAtlasTuneMC09:
case kPyMb:
case kPyMbWithDirectPhoton:
case kPyMbNonDiffr:
Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
}
- if (fQuench) {
+ fPythia->SetPARJ(200, 0.0);
+ fPythia->SetPARJ(199, 0.0);
+ fPythia->SetPARJ(198, 0.0);
+ fPythia->SetPARJ(197, 0.0);
+
+ if (fQuench == 1) {
fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
}
- fPythia->SetPARJ(200, 0.0);
if (fQuench == 3) {
// Nestor's change of the splittings
fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
fPythia->SetMSTJ(50, 0); // No coherence in first branching
fPythia->SetPARJ(82, 1.); // Cut off for parton showers
+ } else if (fQuench == 4) {
+ // Armesto-Cunqueiro-Salgado change of the splittings.
+ AliFastGlauber* glauber = AliFastGlauber::Instance();
+ glauber->Init(2);
+ //read and store transverse almonds corresponding to differnt
+ //impact parameters.
+ glauber->SetCentralityClass(0.,0.1);
+ fPythia->SetPARJ(200, 1.);
+ fPythia->SetPARJ(198, fQhat);
+ fPythia->SetPARJ(199, fLength);
+ fPythia->SetMSTJ(42, 2); // angular ordering
+ fPythia->SetMSTJ(44, 2); // option to run alpha_s
+ fPythia->SetPARJ(82, 1.); // Cut off for parton showers
}
}
// Switch hadronisation off
//
fPythia->SetMSTJ(1, 0);
+
+ if (fQuench ==4){
+ Double_t bimp;
+ // Quenching comes through medium-modified splitting functions.
+ AliFastGlauber::Instance()->GetRandomBHard(bimp);
+ fPythia->SetPARJ(197, bimp);
+ fImpact = bimp;
+ fPythia->Qpygin0();
+ }
//
// Either produce new event or read partons from file
//
}
fNpartons = fPythia->GetN();
} else {
- printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
- fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
+ printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
+ fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
fPythia->SetN(0);
LoadEvent(fRL->Stack(), 0 , 1);
fPythia->Pyedit(21);
//
// .. and perform hadronisation
// printf("Calling hadronisation %d\n", fPythia->GetN());
- fPythia->Pyexec();
+
+ if (fPatchOmegaDalitz) {
+ fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
+ fPythia->Pyexec();
+ fPythia->DalitzDecays();
+ fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
+ }
+ fPythia->Pyexec();
}
fTrials++;
fPythia->ImportParticles(&fParticles,"All");
- Boost();
+
+ if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
+ if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
+
//
//
//
Int_t nTkbles = 0; // Trackable particles
if (fProcess != kPyMbDefault &&
fProcess != kPyMb &&
+ fProcess != kPyMbAtlasTuneMC09 &&
fProcess != kPyMbWithDirectPhoton &&
fProcess != kPyJets &&
fProcess != kPyDirectGamma &&
fProcess != kPyW &&
fProcess != kPyZ &&
fProcess != kPyCharmppMNRwmi &&
- fProcess != kPyBeautyppMNRwmi) {
+ fProcess != kPyBeautyppMNRwmi &&
+ fProcess != kPyBeautyJets) {
for (i = 0; i < np; i++) {
TParticle* iparticle = (TParticle *) fParticles.At(i);
Int_t* pParent = new Int_t[np];
for (i=0; i< np; i++) pParent[i] = -1;
- if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
+ if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
TParticle* jet1 = (TParticle *) fParticles.At(6);
TParticle* jet2 = (TParticle *) fParticles.At(7);
if (!CheckTrigger(jet1, jet2)) {
Int_t pdg = 0;
if (fFragPhotonInCalo) pdg = 22 ; // Photon
- else if (fPi0InCalo) pdg = 111 ; // Pi0
+ else if (fPi0InCalo) pdg = 111 ; // Pi0
for (i=0; i< np; i++) {
TParticle* iparticle = (TParticle *) fParticles.At(i);
if(!ok)
return 0;
}
-
+
+ // Select beauty jets with electron in EMCAL
+ if (fProcess == kPyBeautyJets && fEleInEMCAL) {
+
+ Bool_t ok = kFALSE;
+
+ Int_t pdg = 11; //electron
+
+ Float_t pt = 0.;
+ Float_t eta = 0.;
+ Float_t phi = 0.;
+ for (i=0; i< np; i++) {
+ TParticle* iparticle = (TParticle *) fParticles.At(i);
+ if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg &&
+ iparticle->Pt() > fElectronMinPt){
+ pt = iparticle->Pt();
+ phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
+ eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
+ if(IsInEMCAL(phi,eta))
+ ok =kTRUE;
+ }
+ }
+ if(!ok)
+ return 0;
+ AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
+ }
// Check for minimum multiplicity
if (fTriggerMultiplicity > 0) {
Int_t multiplicity = 0;
Int_t statusCode = iparticle->GetStatusCode();
// Initial state particle
- if (statusCode > 20)
+ if (statusCode != 1)
continue;
-
- // skip quarks and gluons
- Int_t pdgCode = TMath::Abs(iparticle->GetPdgCode());
- if (pdgCode <= 10 || pdgCode == 21)
- continue;
-
+ // eta cut
if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
continue;
-
+ // pt cut
+ if (iparticle->Pt() < fTriggerMultiplicityPtMin)
+ continue;
+
TParticlePDG* pdgPart = iparticle->GetPDG();
if (pdgPart && pdgPart->Charge() == 0)
continue;
delete [] pParent;
return 0;
}
-
- Printf("Triggered on event with multiplicity of %d > %d", multiplicity, fTriggerMultiplicity);
+ Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
}
// Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
// Check if there is a ccbar or bbbar pair with at least one of the two
// in fYMin < y < fYMax
- if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
+
+ if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
TParticle *partCheck;
TParticle *mother;
Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
fProcess == kPyZ ||
fProcess == kPyMbDefault ||
fProcess == kPyMb ||
+ fProcess == kPyMbAtlasTuneMC09 ||
fProcess == kPyMbWithDirectPhoton ||
fProcess == kPyMbNonDiffr)
&& (fCutOnChild == 1) ) {
Int_t km = iparticle->GetFirstMother();
if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
(ks != 1) ||
- (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
+ ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) {
nc++;
if (ks == 1) trackIt = 1;
Int_t ipa = iparticle->GetFirstMother()-1;
//
// Event Vertex
fHeader->SetPrimaryVertex(fVertex);
-
+ fHeader->SetInteractionTime(fEventTime);
//
// Number of primaries
fHeader->SetNProduced(fNprimaries);
//
// Jets that have triggered
- if (fProcess == kPyJets || fProcess == kPyDirectGamma)
+ //Need to store jets for b-jet studies too!
+ if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
{
Int_t ntrig, njet;
Float_t jets[4][10];
if (fQuench == 1) {
// Pythia::Quench()
fPythia->GetQuenchingParameters(xp, yp, z);
- } else {
+ } else if (fQuench == 2){
// Pyquen
Double_t r1 = PARIMP.rb1;
Double_t r2 = PARIMP.rb2;
xp = r * TMath::Cos(phi);
yp = r * TMath::Sin(phi);
+ } else if (fQuench == 4) {
+ // QPythia
+ Double_t xy[2];
+ Double_t i0i1[2];
+ AliFastGlauber::Instance()->GetSavedXY(xy);
+ AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
+ xp = xy[0];
+ yp = xy[1];
+ ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
}
+
((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
- }
+ }
//
// Store pt^hard
((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
pdg[1] = jet2->GetPdgCode();
Bool_t triggered = kFALSE;
- if (fProcess == kPyJets) {
+ if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
Int_t njets = 0;
Int_t ntrig = 0;
Float_t jets[4][10];
void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
{
-//
-// Load event into Pythia Common Block
-//
-
- Int_t npart = stack -> GetNprimary();
- Int_t n0 = 0;
+ //
+ // Load event into Pythia Common Block
+ //
+
+ Int_t npart = stack -> GetNprimary();
+ Int_t n0 = 0;
+
+ if (!flag) {
+ (fPythia->GetPyjets())->N = npart;
+ } else {
+ n0 = (fPythia->GetPyjets())->N;
+ (fPythia->GetPyjets())->N = n0 + npart;
+ }
+
+
+ for (Int_t part = 0; part < npart; part++) {
+ TParticle *mPart = stack->Particle(part);
- if (!flag) {
- (fPythia->GetPyjets())->N = npart;
- } else {
- n0 = (fPythia->GetPyjets())->N;
- (fPythia->GetPyjets())->N = n0 + npart;
+ Int_t kf = mPart->GetPdgCode();
+ Int_t ks = mPart->GetStatusCode();
+ Int_t idf = mPart->GetFirstDaughter();
+ Int_t idl = mPart->GetLastDaughter();
+
+ if (reHadr) {
+ if (ks == 11 || ks == 12) {
+ ks -= 10;
+ idf = -1;
+ idl = -1;
+ }
}
+ Float_t px = mPart->Px();
+ Float_t py = mPart->Py();
+ Float_t pz = mPart->Pz();
+ Float_t e = mPart->Energy();
+ Float_t m = mPart->GetCalcMass();
- for (Int_t part = 0; part < npart; part++) {
- TParticle *mPart = stack->Particle(part);
-
- Int_t kf = mPart->GetPdgCode();
- Int_t ks = mPart->GetStatusCode();
- Int_t idf = mPart->GetFirstDaughter();
- Int_t idl = mPart->GetLastDaughter();
-
- if (reHadr) {
+
+ (fPythia->GetPyjets())->P[0][part+n0] = px;
+ (fPythia->GetPyjets())->P[1][part+n0] = py;
+ (fPythia->GetPyjets())->P[2][part+n0] = pz;
+ (fPythia->GetPyjets())->P[3][part+n0] = e;
+ (fPythia->GetPyjets())->P[4][part+n0] = m;
+
+ (fPythia->GetPyjets())->K[1][part+n0] = kf;
+ (fPythia->GetPyjets())->K[0][part+n0] = ks;
+ (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
+ (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
+ (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
+ }
+}
+
+void AliGenPythia::LoadEvent(TObjArray* stack, Int_t flag, Int_t reHadr)
+{
+ //
+ // Load event into Pythia Common Block
+ //
+
+ Int_t npart = stack -> GetEntries();
+ Int_t n0 = 0;
+
+ if (!flag) {
+ (fPythia->GetPyjets())->N = npart;
+ } else {
+ n0 = (fPythia->GetPyjets())->N;
+ (fPythia->GetPyjets())->N = n0 + npart;
+ }
+
+
+ for (Int_t part = 0; part < npart; part++) {
+ TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
+ Int_t kf = mPart->GetPdgCode();
+ Int_t ks = mPart->GetStatusCode();
+ Int_t idf = mPart->GetFirstDaughter();
+ Int_t idl = mPart->GetLastDaughter();
+
+ if (reHadr) {
if (ks == 11 || ks == 12) {
- ks -= 10;
- idf = -1;
- idl = -1;
+ ks -= 10;
+ idf = -1;
+ idl = -1;
}
- }
-
- Float_t px = mPart->Px();
- Float_t py = mPart->Py();
- Float_t pz = mPart->Pz();
- Float_t e = mPart->Energy();
- Float_t m = mPart->GetCalcMass();
-
-
- (fPythia->GetPyjets())->P[0][part+n0] = px;
- (fPythia->GetPyjets())->P[1][part+n0] = py;
- (fPythia->GetPyjets())->P[2][part+n0] = pz;
- (fPythia->GetPyjets())->P[3][part+n0] = e;
- (fPythia->GetPyjets())->P[4][part+n0] = m;
-
- (fPythia->GetPyjets())->K[1][part+n0] = kf;
- (fPythia->GetPyjets())->K[0][part+n0] = ks;
- (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
- (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
- (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
}
+
+ Float_t px = mPart->Px();
+ Float_t py = mPart->Py();
+ Float_t pz = mPart->Pz();
+ Float_t e = mPart->Energy();
+ Float_t m = mPart->GetCalcMass();
+
+
+ (fPythia->GetPyjets())->P[0][part+n0] = px;
+ (fPythia->GetPyjets())->P[1][part+n0] = py;
+ (fPythia->GetPyjets())->P[2][part+n0] = pz;
+ (fPythia->GetPyjets())->P[3][part+n0] = e;
+ (fPythia->GetPyjets())->P[4][part+n0] = m;
+
+ (fPythia->GetPyjets())->K[1][part+n0] = kf;
+ (fPythia->GetPyjets())->K[0][part+n0] = ks;
+ (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
+ (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
+ (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
+ }
}
//loop for all particles and produce the phi rotation
Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
Double_t oldphi, newphi;
- Double_t newVx, newVy, R, Vz, time;
- Double_t newPx, newPy, pt, Pz, e;
+ Double_t newVx, newVy, r, vZ, time;
+ Double_t newPx, newPy, pt, pz, e;
for(Int_t i=0; i< np; i++) {
TParticle* iparticle = (TParticle *) fParticles.At(i);
oldphi = iparticle->Phi();
if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
- R = iparticle->R();
- newVx = R*TMath::Cos(newphi);
- newVy = R*TMath::Sin(newphi);
- Vz = iparticle->Vz(); // don't transform
+ r = iparticle->R();
+ newVx = r * TMath::Cos(newphi);
+ newVy = r * TMath::Sin(newphi);
+ vZ = iparticle->Vz(); // don't transform
time = iparticle->T(); // don't transform
pt = iparticle->Pt();
- newPx = pt*TMath::Cos(newphi);
- newPy = pt*TMath::Sin(newphi);
- Pz = iparticle->Pz(); // don't transform
- e = iparticle->Energy(); // don't transform
+ newPx = pt * TMath::Cos(newphi);
+ newPy = pt * TMath::Sin(newphi);
+ pz = iparticle->Pz(); // don't transform
+ e = iparticle->Energy(); // don't transform
// apply rotation
- iparticle->SetProductionVertex(newVx, newVy, Vz, time);
- iparticle->SetMomentum(newPx, newPy, Pz, e);
+ iparticle->SetProductionVertex(newVx, newVy, vZ, time);
+ iparticle->SetMomentum(newPx, newPy, pz, e);
} //end particle loop
}
-#ifdef never
-void AliGenPythia::Streamer(TBuffer &R__b)
-{
- // Stream an object of class AliGenPythia.
-
- if (R__b.IsReading()) {
- Version_t R__v = R__b.ReadVersion(); if (R__v) { }
- AliGenerator::Streamer(R__b);
- R__b >> (Int_t&)fProcess;
- R__b >> (Int_t&)fStrucFunc;
- R__b >> (Int_t&)fForceDecay;
- R__b >> fEnergyCMS;
- R__b >> fKineBias;
- R__b >> fTrials;
- fParentSelect.Streamer(R__b);
- fChildSelect.Streamer(R__b);
- R__b >> fXsection;
-// (AliPythia::Instance())->Streamer(R__b);
- R__b >> fPtHardMin;
- R__b >> fPtHardMax;
-// if (fDecayer) fDecayer->Streamer(R__b);
- } else {
- R__b.WriteVersion(AliGenPythia::IsA());
- AliGenerator::Streamer(R__b);
- R__b << (Int_t)fProcess;
- R__b << (Int_t)fStrucFunc;
- R__b << (Int_t)fForceDecay;
- R__b << fEnergyCMS;
- R__b << fKineBias;
- R__b << fTrials;
- fParentSelect.Streamer(R__b);
- fChildSelect.Streamer(R__b);
- R__b << fXsection;
-// R__b << fPythia;
- R__b << fPtHardMin;
- R__b << fPtHardMax;
- // fDecayer->Streamer(R__b);
- }
-}
-#endif
-
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff