#include "AliMC.h"
#include "AliLog.h"
#include "PyquenCommon.h"
-#include "AliLog.h"
ClassImp(AliGenPythia)
fPatchOmegaDalitz(0),
fNpartons(0),
fReadFromFile(0),
+ fReadLHEF(0),
fQuench(0),
fQhat(0.),
fLength(0.),
fFeedDownOpt(kTRUE),
fFragmentation(kTRUE),
fSetNuclei(kFALSE),
+ fUseNuclearPDF(kFALSE),
+ fUseLorentzBoost(kTRUE),
fNewMIS(kFALSE),
fHFoff(kFALSE),
fNucPdf(0),
fTriggerParticle(0),
fTriggerEta(0.9),
+ fTriggerMinPt(-1),
+ fTriggerMaxPt(1000),
fTriggerMultiplicity(0),
fTriggerMultiplicityEta(0),
fTriggerMultiplicityPtMin(0),
fHeader(0),
fRL(0),
fkFileName(0),
+ fkNameLHEF(0),
fFragPhotonInCalo(kFALSE),
+ fHadronInCalo(kFALSE) ,
fPi0InCalo(kFALSE) ,
- fPhotonInCalo(kFALSE),
- fEleInEMCAL(kFALSE),
+ fEtaInCalo(kFALSE) ,
+ fPhotonInCalo(kFALSE), // not in use
+ fDecayPhotonInCalo(kFALSE),
+ fForceNeutralMeson2PhotonDecay(kFALSE),
+ fEleInCalo(kFALSE),
+ fEleInEMCAL(kFALSE), // not in use
+ fCheckBarrel(kFALSE),
fCheckEMCAL(kFALSE),
fCheckPHOS(kFALSE),
fCheckPHOSeta(kFALSE),
- fFragPhotonOrPi0MinPt(0),
- fPhotonMinPt(0),
- fElectronMinPt(0),
+ fPHOSRotateCandidate(-1),
+ fTriggerParticleMinPt(0),
+ fPhotonMinPt(0), // not in use
+ fElectronMinPt(0), // not in use
fPHOSMinPhi(219.),
fPHOSMaxPhi(321.),
fPHOSEta(0.13),
fPatchOmegaDalitz(0),
fNpartons(0),
fReadFromFile(kFALSE),
+ fReadLHEF(0),
fQuench(kFALSE),
fQhat(0.),
fLength(0.),
fFeedDownOpt(kTRUE),
fFragmentation(kTRUE),
fSetNuclei(kFALSE),
+ fUseNuclearPDF(kFALSE),
+ fUseLorentzBoost(kTRUE),
fNewMIS(kFALSE),
fHFoff(kFALSE),
fNucPdf(0),
fTriggerParticle(0),
- fTriggerEta(0.9),
+ fTriggerEta(0.9),
+ fTriggerMinPt(-1),
+ fTriggerMaxPt(1000),
fTriggerMultiplicity(0),
fTriggerMultiplicityEta(0),
fTriggerMultiplicityPtMin(0),
fHeader(0),
fRL(0),
fkFileName(0),
+ fkNameLHEF(0),
fFragPhotonInCalo(kFALSE),
+ fHadronInCalo(kFALSE) ,
fPi0InCalo(kFALSE) ,
- fPhotonInCalo(kFALSE),
- fEleInEMCAL(kFALSE),
+ fEtaInCalo(kFALSE) ,
+ fPhotonInCalo(kFALSE), // not in use
+ fDecayPhotonInCalo(kFALSE),
+ fForceNeutralMeson2PhotonDecay(kFALSE),
+ fEleInCalo(kFALSE),
+ fEleInEMCAL(kFALSE), // not in use
+ fCheckBarrel(kFALSE),
fCheckEMCAL(kFALSE),
fCheckPHOS(kFALSE),
fCheckPHOSeta(kFALSE),
- fFragPhotonOrPi0MinPt(0),
- fPhotonMinPt(0),
- fElectronMinPt(0),
+ fPHOSRotateCandidate(-1),
+ fTriggerParticleMinPt(0),
+ fPhotonMinPt(0), // not in use
+ fElectronMinPt(0), // not in use
fPHOSMinPhi(219.),
fPHOSMaxPhi(321.),
fPHOSEta(0.13),
fPythia->SetCKIN(7,fYHardMin);
fPythia->SetCKIN(8,fYHardMax);
- if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
+ if (fProjectile != "p" || fTarget != "p") fPythia->SetCollisionSystem(fProjectile,fTarget);
+
+ if(fUseNuclearPDF)
+ fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
// Fragmentation?
if (fFragmentation) {
fPythia->SetMSTP(111,1);
fPythia->SetMSTP(91,0);
}
-
+ if (fReadLHEF) fPythia->OpenFortranFile(97, const_cast<char*>(fkNameLHEF));
+
if (fReadFromFile) {
fRL = AliRunLoader::Open(fkFileName, "Partons");
fRL->LoadKinematics();
case kPyCharmpPbMNR:
case kPyCharmppMNR:
case kPyCharmppMNRwmi:
+ case kPyCharmPWHG:
fParentSelect[0] = 411;
fParentSelect[1] = 421;
fParentSelect[2] = 431;
case kPyBeautypPbMNR:
case kPyBeautyppMNR:
case kPyBeautyppMNRwmi:
+ case kPyBeautyPWHG:
fParentSelect[0]= 511;
fParentSelect[1]= 521;
fParentSelect[2]= 531;
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;
}
//
//
//
AliGenMC::Init();
+
+ // Reset Lorentz boost if demanded
+ if(!fUseLorentzBoost) {
+ fDyBoost = 0;
+ Warning("Init","Demand to discard Lorentz boost.\n");
+ }
//
//
//
if (fSetNuclei) {
- fDyBoost = 0;
- Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
+ fDyBoost = 0;
+ Warning("Init","Deprecated function SetNuclei() used (nPDFs + no boost). Use SetProjectile + SetTarget + SetUseNuclearPDF + SetUseLorentzBoost instead.\n");
}
-
fPythia->SetPARJ(200, 0.0);
fPythia->SetPARJ(199, 0.0);
fPythia->SetPARJ(198, 0.0);
fPythia->SetMSTJ(44, 2); // option to run alpha_s
fPythia->SetPARJ(82, 1.); // Cut off for parton showers
}
+
+ if ( AliLog::GetDebugLevel("","AliGenPythia") >= 1 ) {
+ fPythia->Pystat(4);
+ fPythia->Pystat(5);
+ }
}
void AliGenPythia::Generate()
fProcess != kPyMbMSEL1 &&
fProcess != kPyW &&
fProcess != kPyZ &&
+ fProcess != kPyZgamma &&
fProcess != kPyCharmppMNRwmi &&
fProcess != kPyBeautyppMNRwmi &&
- fProcess != kPyBeautyJets) {
+ fProcess != kPyBeautyJets &&
+ fProcess != kPyWPWHG &&
+ fProcess != kPyJetsPWHG &&
+ fProcess != kPyCharmPWHG &&
+ fProcess != kPyBeautyPWHG) {
for (i = 0; i < np; i++) {
TParticle* iparticle = (TParticle *) fParticles.At(i);
// converts from mm/c to s
const Float_t kconv=0.001/2.999792458e8;
-
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 || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
- TParticle* jet1 = (TParticle *) fParticles.At(6);
+ if ((fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
+ && fEtMinJet > 0.) {
+ TParticle* jet1 = (TParticle *) fParticles.At(6);
TParticle* jet2 = (TParticle *) fParticles.At(7);
- if (!CheckTrigger(jet1, jet2)) {
+
+ if (!jet1 || ! jet2 || !CheckTrigger(jet1, jet2)) {
delete [] pParent;
return 0;
}
}
- // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
- if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
-
- Bool_t ok = kFALSE;
-
- Int_t pdg = 0;
- if (fFragPhotonInCalo) pdg = 22 ; // Photon
- else if (fPi0InCalo) pdg = 111 ; // Pi0
-
- for (i=0; i< np; i++) {
- TParticle* iparticle = (TParticle *) fParticles.At(i);
- if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
- iparticle->Pt() > fFragPhotonOrPi0MinPt){
- Int_t imother = iparticle->GetFirstMother() - 1;
- TParticle* pmother = (TParticle *) fParticles.At(imother);
- if(pdg == 111 ||
- (pdg == 22 && pmother->GetStatusCode() != 11)) //No photon from hadron decay
- {
- Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
- Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
- if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
- (fCheckPHOS && IsInPHOS(phi,eta)) )
- ok =kTRUE;
- }
- }
- }
- if(!ok) {
- delete[] pParent;
- return 0;
- }
+
+ // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
+ // implemented primaryly for kPyJets, but extended to any kind of process.
+ if ((fFragPhotonInCalo || fPi0InCalo || fEtaInCalo || fEleInCalo || fHadronInCalo || fDecayPhotonInCalo) &&
+ (fCheckPHOS || fCheckEMCAL || fCheckBarrel) ) {
+ Bool_t ok = TriggerOnSelectedParticles(np);
+
+ if(!ok) {
+ delete[] pParent;
+ 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) {
- delete[] pParent;
- return 0;
- }
-
- AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
- }
// Check for diffraction
if(fkTuneForDiff) {
if(fItune==320 && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
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
- if ((fProcess == kPyJets || 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);
- Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
- Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
-
- if(iparticle->GetStatusCode() == 1
- && iparticle->GetPdgCode() == pdg
- && iparticle->Pt() > fPhotonMinPt
- && eta < fPHOSEta){
-
- // first check if the photon is in PHOS phi
- if(IsInPHOS(phi,eta)){
- okd = kTRUE;
- break;
- }
- if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
-
- }
- }
-
- if(!okd && iphcand != -1) // execute rotation in phi
- RotatePhi(iphcand,okd);
-
- if(!okd) {
- delete [] pParent;
- return 0;
- }
- }
if (fTriggerParticle) {
Bool_t triggered = kFALSE;
if (kf != fTriggerParticle) continue;
if (iparticle->Pt() == 0.) continue;
if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
+ if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
triggered = kTRUE;
break;
}
}
//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 == kPyMbAtlasTuneMC09 ||
Int_t km = iparticle->GetFirstMother();
if (
(((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
- ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)
+ ((fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) && ks == 21 && km == 0 && i>1)
)
{
nc++;
{
// Treat protons as inside nuclei with mass numbers a1 and a2
- fAProjectile = a1;
- fATarget = a2;
- fNucPdf = pdfset; // 0 EKS98 1 EPS08
- fSetNuclei = kTRUE;
+ fAProjectile = a1;
+ fATarget = a2;
+ fNucPdf = pdfset; // 0 EKS98 9 EPS09LO 19 EPS09NLO
+ fUseNuclearPDF = kTRUE;
+ fSetNuclei = kTRUE;
}
// Number of primaries
fHeader->SetNProduced(fNprimaries);
//
+// Event weight
+ fHeader->SetEventWeight(fPythia->GetVINT(97));
+//
// Jets that have triggered
//Need to store jets for b-jet studies too!
- if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
+ if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
{
Int_t ntrig, njet;
Float_t jets[4][10];
((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
}
//
-// Store pt^hard
+// Store pt^hard and cross-section
((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
+ ((AliGenPythiaEventHeader*) fHeader)->SetXsection(fPythia->GetPARI(1));
+
+//
+// Store Event Weight
+ ((AliGenPythiaEventHeader*) fHeader)->SetEventWeight(fPythia->GetPARI(7));
+
//
// Pass header
//
pdg[1] = jet2->GetPdgCode();
Bool_t triggered = kFALSE;
- if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
+ if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) {
Int_t njets = 0;
Int_t ntrig = 0;
Float_t jets[4][10];
return;
}
+Bool_t AliGenPythia::IsInBarrel(Float_t eta) const
+{
+ // Is particle in Central Barrel acceptance?
+ // etamin=-etamax
+ if( eta < fTriggerEta )
+ return kTRUE;
+ else
+ return kFALSE;
+}
+
Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
{
// Is particle in EMCAL acceptance?
return kFALSE;
}
-Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta) const
+Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta, Int_t iparticle)
{
// Is particle in PHOS acceptance?
// Acceptance slightly larger considered.
// phi in degrees, etamin=-etamax
+ // iparticle is the index of the particle to be checked, for PHOS rotation case
+
if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
eta < fPHOSEta )
return kTRUE;
else
+ {
+ if( fCheckPHOSeta && eta < fPHOSEta) fPHOSRotateCandidate = iparticle;
+
return kFALSE;
+ }
}
-void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
+void AliGenPythia::RotatePhi(Bool_t& okdd)
{
+ //Rotate event in phi to enhance events in PHOS acceptance
+
+ if(fPHOSRotateCandidate < 0) return ;
+
//calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
//calculate deltaphi
- TParticle* ph = (TParticle *) fParticles.At(iphcand);
+ TParticle* ph = (TParticle *) fParticles.At(fPHOSRotateCandidate);
Double_t phphi = ph->Phi();
Double_t deltaphi = phiPHOS - phphi;
-
+
//loop for all particles and produce the phi rotation
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();
- newphi = oldphi + deltaphi;
- 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
- 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
-
- // apply rotation
- iparticle->SetProductionVertex(newVx, newVy, vZ, time);
- iparticle->SetMomentum(newPx, newPy, pz, e);
-
+ TParticle* iparticle = (TParticle *) fParticles.At(i);
+ oldphi = iparticle->Phi();
+ newphi = oldphi + deltaphi;
+ 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
+ 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
+
+ // apply rotation
+ iparticle->SetProductionVertex(newVx, newVy, vZ, time);
+ iparticle->SetMomentum(newPx, newPy, pz, e);
+
} //end particle loop
- // now let's check that we put correctly the candidate photon in PHOS
- Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
- Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
- if(IsInPHOS(phi,eta))
- okdd = kTRUE;
+ // now let's check that we put correctly the candidate photon in PHOS
+ Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
+ Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
+ if(IsInPHOS(phi,eta,-1))
+ okdd = kTRUE;
+
+ // reset the value for next event
+ fPHOSRotateCandidate = -1;
+
}
-
-
Bool_t AliGenPythia::CheckDiffraction()
{
// use this method only with Perugia-0 tune!
TParticle * part = (TParticle *) fParticles.At(iPart);
Double_t E= part->Energy();
Double_t P= part->P();
- M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
+ Double_t M2 = (fEnergyCMS-E-P)*(fEnergyCMS-E+P);
+ if(M2<0) return kFALSE;
+ M= TMath::Sqrt(M2);
}
Double_t Mmin, Mmax, wSD, wDD, wND;
}
return kFALSE;
}
+
+Bool_t AliGenPythia::CheckDetectorAcceptance(Float_t phi, Float_t eta, Int_t iparticle)
+{
+ // check the eta/phi correspond to the detectors acceptance
+ // iparticle is the index of the particle to be checked, for PHOS rotation case
+ if (fCheckPHOS && IsInPHOS (phi,eta,iparticle)) return kTRUE;
+ else if(fCheckEMCAL && IsInEMCAL (phi,eta)) return kTRUE;
+ else if(fCheckBarrel && IsInBarrel( eta)) return kTRUE;
+ else return kFALSE;
+}
+
+Bool_t AliGenPythia::TriggerOnSelectedParticles(Int_t np)
+{
+ // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
+ // implemented primaryly for kPyJets, but extended to any kind of process.
+
+ //printf("Check: frag photon %d, pi0 %d, eta %d, electron %d, hadron %d, decay %d, PHOS %d, EMCAL %d, Barrel %d \n",
+ // fFragPhotonInCalo,fPi0InCalo, fEtaInCalo,fEleInCalo,fHadronInCalo,fDecayPhotonInCalo,fCheckPHOS,fCheckEMCAL, fCheckBarrel);
+
+ Bool_t ok = kFALSE;
+ for (Int_t i=0; i< np; i++) {
+
+ TParticle* iparticle = (TParticle *) fParticles.At(i);
+
+ Int_t pdg = iparticle->GetPdgCode();
+ Int_t status = iparticle->GetStatusCode();
+ Int_t imother = iparticle->GetFirstMother() - 1;
+
+ TParticle* pmother = 0x0;
+ Int_t momStatus = -1;
+ Int_t momPdg = -1;
+ if(imother > 0 ){
+ pmother = (TParticle *) fParticles.At(imother);
+ momStatus = pmother->GetStatusCode();
+ momPdg = pmother->GetPdgCode();
+ }
+
+ ok = kFALSE;
+
+ //
+ // Check the particle type: hadron (not pi0 or eta), electron, decay photon (from pi0 or eta or any), pi0 or eta
+ //
+ // Hadron
+ if (fHadronInCalo && status == 1)
+ {
+ if(TMath::Abs(pdg) > 23 && pdg !=221 && pdg != 111) // avoid photons, electrons, muons, neutrinos and eta or pi0
+ // (in case neutral mesons were declared stable)
+ ok = kTRUE;
+ }
+ //Electron
+ else if (fEleInCalo && status == 1 && TMath::Abs(pdg) == 11)
+ {
+ ok = kTRUE;
+ }
+ //Fragmentation photon
+ else if (fFragPhotonInCalo && pdg == 22 && status == 1)
+ {
+ if(momStatus != 11) ok = kTRUE ; // No photon from hadron decay
+ }
+ // Decay photon
+ else if (fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay && pdg == 22) // pi0 status can be 1 or 11 depending on decay settings, work only for 11
+ {
+ if( momStatus == 11)
+ {
+ //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Decay photon! pdg %d, status %d, pt %2.2f, mom: pdg %d, pt %2.2f\n",
+ // pdg,status,iparticle->Pt(),momPdg,pmother->Pt());
+ ok = kTRUE ; // photon from hadron decay
+
+ //In case only decays from pi0 or eta requested
+ if(fPi0InCalo && momPdg!=111) ok = kFALSE;
+ if(fEtaInCalo && momPdg!=221) ok = kFALSE;
+ }
+
+ }
+ // Pi0 or Eta particle
+ else if ((fPi0InCalo || fEtaInCalo))
+ {
+ if(fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay ) continue ;
+
+ if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
+ {
+ //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
+ ok = kTRUE;
+ }
+ else if (fEtaInCalo && pdg == 221)
+ {
+ //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Eta! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
+ ok = kTRUE;
+ }
+
+ }// pi0 or eta
+
+ //
+ // Check that the selected particle is in the calorimeter acceptance
+ //
+ if(ok && iparticle->Pt() > fTriggerParticleMinPt)
+ {
+ //Just check if the selected particle falls in the acceptance
+ if(!fForceNeutralMeson2PhotonDecay )
+ {
+ //printf("\t Check acceptance! \n");
+ Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
+ Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
+
+ if(CheckDetectorAcceptance(phi,eta,i))
+ {
+ ok =kTRUE;
+ AliDebug(1,Form("Selected trigger pdg %d, status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
+ //printf("\t Accept \n");
+ break;
+ }
+ else ok = kFALSE;
+ }
+ //Mesons have several decay modes, select only those decaying into 2 photons
+ else if(fForceNeutralMeson2PhotonDecay && (fPi0InCalo || fEtaInCalo))
+ {
+ // In case we want the pi0/eta trigger,
+ // check the decay mode (2 photons)
+
+ //printf("\t Force decay 2 gamma\n");
+
+ Int_t ndaughters = iparticle->GetNDaughters();
+ if(ndaughters != 2){
+ ok=kFALSE;
+ continue;
+ }
+
+ TParticle*d1 = (TParticle *) fParticles.At(iparticle->GetDaughter(0)-1);
+ TParticle*d2 = (TParticle *) fParticles.At(iparticle->GetDaughter(1)-1);
+ if(!d1 || !d2) {
+ ok=kFALSE;
+ continue;
+ }
+
+ //iparticle->Print();
+ //d1->Print();
+ //d2->Print();
+
+ Int_t pdgD1 = d1->GetPdgCode();
+ Int_t pdgD2 = d2->GetPdgCode();
+ //printf("\t \t 2 daughters pdg = %d - %d\n",pdgD1,pdgD2);
+ //printf("mother %d - %d\n",d1->GetFirstMother(),d2->GetFirstMother());
+
+ if(pdgD1 != 22 || pdgD2 != 22){
+ ok = kFALSE;
+ continue;
+ }
+
+ //printf("\t accept decay\n");
+
+ //Trigger on the meson, not on the daughter
+ if(!fDecayPhotonInCalo){
+
+ Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
+ Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
+
+ if(CheckDetectorAcceptance(phi,eta,i))
+ {
+ //printf("\t Accept meson pdg %d\n",pdg);
+ ok =kTRUE;
+ AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
+ break;
+ } else {
+ ok=kFALSE;
+ continue;
+ }
+ }
+
+ //printf("Check daughters acceptance\n");
+
+ //Trigger on the meson daughters
+ //Photon 1
+ Float_t phi = d1->Phi()*180./TMath::Pi(); //Convert to degrees
+ Float_t eta =TMath::Abs(d1->Eta()); //in calos etamin=-etamax
+ if(d1->Pt() > fTriggerParticleMinPt)
+ {
+ //printf("\t Check acceptance photon 1! \n");
+ if(CheckDetectorAcceptance(phi,eta,i))
+ {
+ //printf("\t Accept Photon 1\n");
+ ok =kTRUE;
+ AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
+ break;
+ }
+ else ok = kFALSE;
+ } // pt cut
+ else ok = kFALSE;
+
+ //Photon 2
+ phi = d2->Phi()*180./TMath::Pi(); //Convert to degrees
+ eta =TMath::Abs(d2->Eta()); //in calos etamin=-etamax
+
+ if(d2->Pt() > fTriggerParticleMinPt)
+ {
+ //printf("\t Check acceptance photon 2! \n");
+ if(CheckDetectorAcceptance(phi,eta,i))
+ {
+ //printf("\t Accept Photon 2\n");
+ ok =kTRUE;
+ AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
+ break;
+ }
+ else ok = kFALSE;
+ } // pt cut
+ else ok = kFALSE;
+ } // force 2 photon daughters in pi0/eta decays
+ else ok = kFALSE;
+ } else ok = kFALSE; // check acceptance
+ } // primary loop
+
+ //
+ // If requested, rotate the particles event in phi to enhance/speed PHOS selection
+ // A particle passing all trigger conditions except phi position in PHOS, is used as reference
+ //
+ if(fCheckPHOSeta)
+ {
+ RotatePhi(ok);
+ }
+
+ return ok;
+}
+