fPythia(0),
fProcess(kPyCharm),
fStrucFunc(kCTEQ5L),
- fEnergyCMS(5500.),
fKineBias(0.),
fTrials(0),
fTrialsRun(0),
fPHOSEta(0.13),
fEMCALMinPhi(79.),
fEMCALMaxPhi(191.),
- fEMCALEta(0.71)
-
+ fEMCALEta(0.71),
+ fItune(-1),
+ fInfo(1)
{
// Default Constructor
+ fEnergyCMS = 5500.;
SetNuclei(0,0);
if (!AliPythiaRndm::GetPythiaRandom())
AliPythiaRndm::SetPythiaRandom(GetRandom());
fPythia(pythia),
fProcess(kPyCharm),
fStrucFunc(kCTEQ5L),
- fEnergyCMS(5500.),
fKineBias(0.),
fTrials(0),
fTrialsRun(0),
fPHOSEta(0.13),
fEMCALMinPhi(79.),
fEMCALMaxPhi(191.),
- fEMCALEta(0.71)
+ fEMCALEta(0.71),
+ fItune(-1),
+ fInfo(1)
{
// default charm production at 5. 5 TeV
// semimuonic decay
// structure function GRVHO
//
+ fEnergyCMS = 5500.;
fName = "Pythia";
fTitle= "Particle Generator using PYTHIA";
SetForceDecay();
// Set random number generator
if (!AliPythiaRndm::GetPythiaRandom())
AliPythiaRndm::SetPythiaRandom(GetRandom());
- fParticles = new TClonesArray("TParticle",1000);
SetNuclei(0,0);
}
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
case kPyCharmpPbMNR:
case kPyCharmppMNR:
case kPyCharmppMNRwmi:
+ case kPyCharmPWHG:
fParentSelect[0] = 411;
fParentSelect[1] = 421;
fParentSelect[2] = 431;
fParentSelect[3] = 4122;
+ fParentSelect[4] = 4232;
+ fParentSelect[5] = 4132;
+ fParentSelect[6] = 4332;
fFlavorSelect = 4;
break;
case kPyD0PbPbMNR:
fParentSelect[0] = 431;
fFlavorSelect = 4;
break;
+ case kPyLambdacppMNR:
+ fParentSelect[0] = 4122;
+ fFlavorSelect = 4;
+ break;
case kPyBeauty:
+ case kPyBeautyJets:
case kPyBeautyPbPbMNR:
case kPyBeautypPbMNR:
case kPyBeautyppMNR:
case kPyBeautyppMNRwmi:
+ case kPyBeautyPWHG:
fParentSelect[0]= 511;
fParentSelect[1]= 521;
fParentSelect[2]= 531;
case kPyJpsi:
fParentSelect[0] = 443;
break;
+ case kPyMbAtlasTuneMC09:
case kPyMbDefault:
case kPyMb:
+ case kPyMbWithDirectPhoton:
case kPyMbNonDiffr:
case kPyMbMSEL1:
case kPyJets:
+ case kPyJetsPWHG:
case kPyDirectGamma:
case kPyLhwgMb:
break;
+ case kPyWPWHG:
case kPyW:
case kPyZ:
+ case kPyZgamma:
+ case kPyMBRSingleDiffraction:
+ case kPyMBRDoubleDiffraction:
+ case kPyMBRCentralDiffraction:
break;
}
//
fDecayer->ForceDecay();
- Float_t polar[3] = {0,0,0};
- Float_t origin[3] = {0,0,0};
- Float_t p[4];
+ Double_t polar[3] = {0,0,0};
+ Double_t origin[3] = {0,0,0};
+ Double_t p[4];
// converts from mm/c to s
const Float_t kconv=0.001/2.999792458e8;
//
}
fNpartons = fPythia->GetNumberOfParticles();
} 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->SetNumberOfParticles(0);
fPythia->LoadEvent(fRL->Stack(), 0 , 1);
fPythia->EditEventList(21);
// printf("Calling hadronisation %d\n", fPythia->GetN());
// fPythia->HadronizeEvent();
fTrials++;
- fPythia->GetParticles(fParticles);
+ fPythia->GetParticles(&fParticles);
Boost();
//
//
Int_t i;
fNprimaries = 0;
- Int_t np = fParticles->GetEntriesFast();
+ Int_t np = fParticles.GetEntriesFast();
if (np == 0) continue;
//
Int_t nTkbles = 0; // Trackable particles
if (fProcess != kPyMbDefault &&
fProcess != kPyMb &&
+ fProcess != kPyMbWithDirectPhoton &&
fProcess != kPyJets &&
fProcess != kPyDirectGamma &&
fProcess != kPyMbNonDiffr &&
fProcess != kPyMbMSEL1 &&
fProcess != kPyW &&
fProcess != kPyZ &&
+ fProcess != kPyZgamma &&
fProcess != kPyCharmppMNRwmi &&
- fProcess != kPyBeautyppMNRwmi) {
+ fProcess != kPyBeautyppMNRwmi &&
+ fProcess != kPyWPWHG &&
+ fProcess != kPyJetsPWHG &&
+ fProcess != kPyCharmPWHG &&
+ fProcess != kPyBeautyPWHG) {
for (i = 0; i < np; i++) {
- TParticle* iparticle = (TParticle *) fParticles->At(i);
+ TParticle* iparticle = (TParticle *) fParticles.At(i);
Int_t ks = iparticle->GetStatusCode();
kf = CheckPDGCode(iparticle->GetPdgCode());
// No initial state partons
if (kf >= fFlavorSelect && kf <= 6) {
Int_t idau = (fPythia->Version() == 6) ? (iparticle->GetFirstDaughter() - 1) :(iparticle->GetFirstDaughter());
if (idau > -1) {
- TParticle* daughter = (TParticle *) fParticles->At(idau);
+ TParticle* daughter = (TParticle *) fParticles.At(idau);
Int_t pdgD = daughter->GetPdgCode();
if (pdgD == 91 || pdgD == 92) {
Int_t jmin = (fPythia->Version() == 6) ? (daughter->GetFirstDaughter() - 1) : (daughter->GetFirstDaughter());
Int_t jmax = (fPythia->Version() == 6) ? (daughter->GetLastDaughter() - 1) : (daughter->GetLastDaughter());
- for (Int_t j = jmin; j <= jmax; j++)
- ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
+ for (Int_t jp = jmin; jp <= jmax; jp++)
+ ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
} // is string or cluster
} // has daughter
} // heavy quark
if (ipa > -1) {
- TParticle * mother = (TParticle *) fParticles->At(ipa);
+ TParticle * mother = (TParticle *) fParticles.At(ipa);
kfMo = TMath::Abs(mother->GetPdgCode());
}
if (nc > 0) {
for (i = 0; i < np; i++) {
if (!pSelected[i]) continue;
- TParticle * iparticle = (TParticle *) fParticles->At(i);
+ TParticle * iparticle = (TParticle *) fParticles.At(i);
kf = CheckPDGCode(iparticle->GetPdgCode());
Int_t ks = iparticle->GetStatusCode();
p[0] = iparticle->Px();
origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
- Float_t tof = kconv*iparticle->T();
+ Float_t tof = fTime + kconv*iparticle->T();
Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
}
if (jev >= fNpart || fNpart == -1) {
fKineBias=Float_t(fNpart)/Float_t(fTrials);
- fPythia->GetXandQ(fQ, fX1, fX2);
+ if (fInfo) fPythia->GetXandQ(fX1, fX2, fQ);
fTrialsRun += fTrials;
fNev++;
MakeHeader();
// converts from mm/c to s
const Float_t kconv = 0.001 / 2.999792458e8;
- Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
+ Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
Int_t* pParent = new Int_t[np];
for (i=0; i< np; i++) pParent[i] = -1;
- if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
- TParticle* jet1 = (TParticle *) fParticles->At(6);
- TParticle* jet2 = (TParticle *) fParticles->At(7);
+ if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG ) {
+ TParticle* jet1 = (TParticle *) fParticles.At(6);
+ TParticle* jet2 = (TParticle *) fParticles.At(7);
if (!CheckTrigger(jet1, jet2)) {
delete [] pParent;
return 0;
}
// Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
- if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
+ if ((fProcess == kPyJets || fProcess == kPyJetsPWHG) && (fFragPhotonInCalo || fPi0InCalo) ) {
Bool_t ok = kFALSE;
else if (fPi0InCalo) pdg = 111 ; // Pi0
for (i=0; i< np; i++) {
- TParticle* iparticle = (TParticle *) fParticles->At(i);
+ TParticle* iparticle = (TParticle *) fParticles.At(i);
if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
iparticle->Pt() > fFragPhotonOrPi0MinPt){
Int_t imother = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
- TParticle* pmother = (TParticle *) fParticles->At(imother);
+ TParticle* pmother = (TParticle *) fParticles.At(imother);
if(pdg == 111 ||
(pdg == 22 && pmother->GetStatusCode() != 11))//No photon from hadron decay
{
}
}
}
- if(!ok)
- return 0;
+ if(!ok){
+ delete [] pParent;
+ return 0;
+ }
}
// Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
- if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
+ if ((fProcess == kPyJets || fProcess == kPyJetsPWHG || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
Bool_t okd = kFALSE;
Int_t pdg = 22;
Int_t iphcand = -1;
for (i=0; i< np; i++) {
- TParticle* iparticle = (TParticle *) fParticles->At(i);
+ TParticle* iparticle = (TParticle *) fParticles.At(i);
Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
if(!okd && iphcand != -1) // execute rotation in phi
RotatePhi(iphcand,okd);
- if(!okd)
- return 0;
+ if(!okd) {
+ delete[] pParent;
+ return 0;
+ }
}
if (fTriggerParticle) {
Bool_t triggered = kFALSE;
for (i = 0; i < np; i++) {
- TParticle * iparticle = (TParticle *) fParticles->At(i);
+ TParticle * iparticle = (TParticle *) fParticles.At(i);
kf = CheckPDGCode(iparticle->GetPdgCode());
if (kf != fTriggerParticle) continue;
if (iparticle->Pt() == 0.) continue;
// 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) {
- TParticle *hvq;
+ TParticle *partCheck;
+ TParticle *mother;
Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
- Float_t yQ;
- Int_t pdgQ;
+ Bool_t theChild=kFALSE;
+ Float_t y;
+ Int_t pdg,mpdg,mpdgUpperFamily;
for(i=0; i<np; i++) {
- hvq = (TParticle*)fParticles->At(i);
- pdgQ = hvq->GetPdgCode();
- if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
- if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
- yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
- (hvq->Energy()-hvq->Pz()+1.e-13));
- if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
+ partCheck = (TParticle*)fParticles.At(i);
+ pdg = partCheck->GetPdgCode();
+ if(TMath::Abs(pdg) == fFlavorSelect) { // quark
+ if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
+ y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
+ (partCheck->Energy()-partCheck->Pz()+1.e-13));
+ if(y>fYMin && y<fYMax) inYcut=kTRUE;
+ }
+
+ if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
+ Int_t mi = partCheck->GetFirstMother() - 1;
+ if(mi<0) continue;
+ mother = (TParticle*)fParticles.At(mi);
+ mpdg=TMath::Abs(mother->GetPdgCode());
+ mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
+ if ( ParentSelected(mpdg) ||
+ (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
+ if (KinematicSelection(partCheck,1)) {
+ theChild=kTRUE;
+ }
+ }
+ }
}
- if (!theQ || !theQbar || !inYcut) {
+ if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
delete[] pParent;
return 0;
}
+ if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
+ delete[] pParent;
+ return 0;
+ }
+
}
//Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
- if ( (fProcess == kPyW ||
+ if ( (
+ fProcess == kPyWPWHG ||
+ fProcess == kPyW ||
fProcess == kPyZ ||
+ fProcess == kPyZgamma ||
fProcess == kPyMbDefault ||
fProcess == kPyMb ||
+ fProcess == kPyMbWithDirectPhoton ||
fProcess == kPyMbNonDiffr)
&& (fCutOnChild == 1) ) {
if ( !CheckKinematicsOnChild() ) {
for (i = 0; i < np; i++) {
Int_t trackIt = 0;
- TParticle * iparticle = (TParticle *) fParticles->At(i);
+ TParticle * iparticle = (TParticle *) fParticles.At(i);
kf = CheckPDGCode(iparticle->GetPdgCode());
Int_t ks = iparticle->GetStatusCode();
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 == kPyJetsPWHG) && ks == 21 && km == 0 && i>1)) {
nc++;
if (ks == 1) trackIt = 1;
origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
- Float_t tof = fEventTime + kconv * iparticle->T();
+ Float_t tof = fTime + fEventTime + kconv * iparticle->T();
PushTrack(fTrackIt*trackIt, iparent, kf,
p[0], p[1], p[2], p[3],
//
// Special Treatment to store color-flow
//
+ /*
if (ks == 3 || ks == 13 || ks == 14) {
TParticle* particle = 0;
if (fStack) {
particle = fStack->Particle(nt);
} else {
- particle = gAlice->Stack()->Particle(nt);
+ particle = AliRunLoader::Instance()->Stack()->Particle(nt);
}
-// particle->SetFirstDaughter(fPythia->GetK(2, i));
-// particle->SetLastDaughter(fPythia->GetK(3, i));
+ particle->SetFirstDaughter(fPythia->GetK(2, i));
+ particle->SetLastDaughter(fPythia->GetK(3, i));
}
-
+ */
KeepTrack(nt);
pParent[i] = nt;
SetHighWaterMark(nt);
//
// Event Vertex
fHeader->SetPrimaryVertex(fVertex);
-
+ fHeader->SetInteractionTime(fTime+fEventTime);
//
// Number of primaries
fHeader->SetNProduced(fNprimaries);
//
// Jets that have triggered
- if (fProcess == kPyJets)
+ if (fProcess == kPyJets || fProcess == kPyJetsPWHG)
{
Int_t ntrig, njet;
Float_t jets[4][10];
// Store quenching parameters
//
if (fQuench){
- Double_t z[4];
- Double_t xp, yp;
+ Double_t z[4] = {0.};
+ Double_t xp = 0.;
+ Double_t yp = 0.;
if (fQuench == 1) {
// Pythia::Quench()
fPythia->GetQuenchingParameters(xp, yp, z);
fHeader = 0x0;
}
-Bool_t AliGenPythiaPlus::CheckTrigger(TParticle* jet1, TParticle* jet2)
+Bool_t AliGenPythiaPlus::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
{
// Check the kinematic trigger condition
//
pdg[1] = jet2->GetPdgCode();
Bool_t triggered = kFALSE;
- if (fProcess == kPyJets) {
+ if (fProcess == kPyJets || fProcess == kPyJetsPWHG) {
Int_t njets = 0;
Int_t ntrig = 0;
Float_t jets[4][10];
Int_t nPartAcc = 0; //number of particles in the acceptance range
Int_t numberOfAcceptedParticles = 1;
if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
- Int_t npart = fParticles->GetEntriesFast();
+ Int_t npart = fParticles.GetEntriesFast();
for (j = 0; j<npart; j++) {
- TParticle * jparticle = (TParticle *) fParticles->At(j);
+ TParticle * jparticle = (TParticle *) fParticles.At(j);
kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
ks = jparticle->GetStatusCode();
km = jparticle->GetFirstMother();
-Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta)
+Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta) const
{
// Is particle in EMCAL acceptance?
// phi in degrees, etamin=-etamax
return kFALSE;
}
-Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta)
+Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta) const
{
// Is particle in PHOS acceptance?
// Acceptance slightly larger considered.
Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
//calculate deltaphi
- TParticle* ph = (TParticle *) fParticles->At(iphcand);
+ TParticle* ph = (TParticle *) fParticles.At(iphcand);
Double_t phphi = ph->Phi();
Double_t deltaphi = phiPHOS - phphi;
//loop for all particles and produce the phi rotation
- Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
+ 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;
for(Int_t i=0; i< np; i++) {
- TParticle* iparticle = (TParticle *) fParticles->At(i);
+ TParticle* iparticle = (TParticle *) fParticles.At(i);
oldphi = iparticle->Phi();
newphi = oldphi + deltaphi;
if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle