#include "AliMC.h"
#include "AliLog.h"
#include "PyquenCommon.h"
-#include "AliLog.h"
ClassImp(AliGenPythia)
fFlavorSelect(0),
fXsection(0.),
fPythia(0),
+ fWeightPower(0.),
fPtHardMin(0.),
fPtHardMax(1.e4),
fYHardMin(-1.e10),
fYHardMax(1.e10),
fGinit(1),
fGfinal(1),
+ fCRoff(0),
fHadronisation(1),
fPatchOmegaDalitz(0),
+ fDecayerExodus(0),
fNpartons(0),
fReadFromFile(0),
+ fReadLHEF(0),
fQuench(0),
fQhat(0.),
fLength(0.),
fEtaMaxGamma(20.),
fPhiMinGamma(0.),
fPhiMaxGamma(2. * TMath::Pi()),
+ fUseYCutHQ(kFALSE),
+ fYMinHQ(-20.),
+ fYMaxHQ(20.),
fPycellEtaMax(2.),
fPycellNEta(274),
fPycellNPhi(432),
fFeedDownOpt(kTRUE),
fFragmentation(kTRUE),
fSetNuclei(kFALSE),
+ fUseNuclearPDF(kFALSE),
+ fUseLorentzBoost(kTRUE),
fNewMIS(kFALSE),
fHFoff(kFALSE),
fNucPdf(0),
fTriggerParticle(0),
fTriggerEta(0.9),
+ fTriggerY(999.),
+ fTriggerEtaMin(0.9),
fTriggerMinPt(-1),
fTriggerMaxPt(1000),
fTriggerMultiplicity(0),
fTriggerMultiplicityEta(0),
+ fTriggerMultiplicityEtaMin(0),
+ fTriggerMultiplicityEtaMax(0),
fTriggerMultiplicityPtMin(0),
fCountMode(kCountAll),
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),
+ fPHOSRotateCandidate(-1),
fTriggerParticleMinPt(0),
- fPhotonMinPt(0),
- fElectronMinPt(0),
+ fPhotonMinPt(0), // not in use
+ fElectronMinPt(0), // not in use
fPHOSMinPhi(219.),
fPHOSMaxPhi(321.),
fPHOSEta(0.13),
fFlavorSelect(0),
fXsection(0.),
fPythia(0),
+ fWeightPower(0.),
fPtHardMin(0.),
fPtHardMax(1.e4),
fYHardMin(-1.e10),
fYHardMax(1.e10),
fGinit(kTRUE),
fGfinal(kTRUE),
+ fCRoff(kFALSE),
fHadronisation(kTRUE),
- fPatchOmegaDalitz(0),
+ fPatchOmegaDalitz(0),
+ fDecayerExodus(0),
fNpartons(0),
fReadFromFile(kFALSE),
+ fReadLHEF(0),
fQuench(kFALSE),
fQhat(0.),
fLength(0.),
fEtaMaxGamma(20.),
fPhiMinGamma(0.),
fPhiMaxGamma(2. * TMath::Pi()),
+ fUseYCutHQ(kFALSE),
+ fYMinHQ(-20.),
+ fYMaxHQ(20.),
fPycellEtaMax(2.),
fPycellNEta(274),
fPycellNPhi(432),
fFeedDownOpt(kTRUE),
fFragmentation(kTRUE),
fSetNuclei(kFALSE),
+ fUseNuclearPDF(kFALSE),
+ fUseLorentzBoost(kTRUE),
fNewMIS(kFALSE),
fHFoff(kFALSE),
fNucPdf(0),
fTriggerParticle(0),
fTriggerEta(0.9),
+ fTriggerY(999.),
+ fTriggerEtaMin(0.9),
fTriggerMinPt(-1),
fTriggerMaxPt(1000),
fTriggerMultiplicity(0),
fTriggerMultiplicityEta(0),
+ fTriggerMultiplicityEtaMin(0),
+ fTriggerMultiplicityEtaMax(0),
fTriggerMultiplicityPtMin(0),
fCountMode(kCountAll),
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),
+ fPHOSRotateCandidate(-1),
fTriggerParticleMinPt(0),
- fPhotonMinPt(0),
- fElectronMinPt(0),
+ fPhotonMinPt(0), // not in use
+ fElectronMinPt(0), // not in use
fPHOSMinPhi(219.),
fPHOSMaxPhi(321.),
fPHOSEta(0.13),
//
fParentWeight=1./Float_t(fNpart);
//
-
-
+ if (fWeightPower != 0)
+ fPythia->SetWeightPower(fWeightPower);
fPythia->SetCKIN(3,fPtHardMin);
fPythia->SetCKIN(4,fPtHardMax);
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(61,fGinit);
// final state radiation
fPythia->SetMSTP(71,fGfinal);
+ //color reconnection strength
+ if(fCRoff==1)fPythia->SetMSTP(95,0);
// pt - kick
if (fPtKick > 0.) {
fPythia->SetMSTP(91,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()
fPythia->DalitzDecays();
fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
}
- fPythia->Pyexec();
- }
+
+ else if (fDecayerExodus) {
+
+ fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
+ fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
+ fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 0);
+ fPythia->Pyexec();
+ fPythia->OmegaDalitz();
+ fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
+ fPythia->PizeroDalitz();
+ fPythia->PhiDalitz();
+ fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 1);
+ fPythia->EtaDalitz();
+ fPythia->EtaprimeDalitz();
+ fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
+ fPythia->RhoDirect();
+ fPythia->OmegaDirect();
+ fPythia->PhiDirect();
+ fPythia->JPsiDirect();
+ }
+
+ fPythia->Pyexec();
+ }
fTrials++;
fPythia->ImportParticles(&fParticles,"All");
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 or hadrons going to PHOS or EMCAL
- if ( fProcess == kPyJets &&
- (fFragPhotonInCalo || fPi0InCalo || fHadronInCalo) &&
- (fCheckPHOS || fCheckEMCAL) ) {
-
- Bool_t ok = kFALSE;
-
- for (i=0; i< np; i++) {
-
- TParticle* iparticle = (TParticle *) fParticles.At(i);
-
- Int_t pdg = iparticle->GetPdgCode();
- Int_t status = iparticle->GetStatusCode();
- ok = kFALSE;
-
- if (fFragPhotonInCalo && pdg == 22 && status == 1)
- {
- Int_t imother = iparticle->GetFirstMother() - 1;
- TParticle* pmother = (TParticle *) fParticles.At(imother);
-
- if(pmother->GetStatusCode() != 11) ok = kTRUE ; // No photon from hadron decay
- }
- else if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
- {
- //printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
- ok = kTRUE;
- }
- else 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;
- }
-
- if(ok && iparticle->Pt() > fTriggerParticleMinPt)
- {
- 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;
- 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));
- break;
- }
- }
- else ok = kFALSE;
- }
+
+ // 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;
}
}
- // 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)) ) {
// eta cut
if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
continue;
+ //multiplicity check for a given eta range
+ if ((fTriggerMultiplicityEtaMin != fTriggerMultiplicityEtaMax) &&
+ (iparticle->Eta() < fTriggerMultiplicityEtaMin || iparticle->Eta() > fTriggerMultiplicityEtaMax))
+ continue;
// pt cut
if (iparticle->Pt() < fTriggerMultiplicityPtMin)
continue;
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;
kf = CheckPDGCode(iparticle->GetPdgCode());
if (kf != fTriggerParticle) continue;
if (iparticle->Pt() == 0.) continue;
- if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
+ if (TMath::Abs(iparticle->Y()) > fTriggerY) continue;
+ if (fTriggerEtaMin == fTriggerEta) {
+ if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
+ } else {
+ if (iparticle->Eta() > fTriggerEta || iparticle->Eta() < fTriggerEtaMin) continue;
+ }
if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
triggered = kTRUE;
break;
Bool_t theChild=kFALSE;
Bool_t triggered=kFALSE;
Float_t y;
- Int_t pdg,mpdg,mpdgUpperFamily;
- for(i=0; i<np; i++) {
+ Int_t pdg, mpdg, mpdgUpperFamily;
+ for(i = 0; i < np; i++) {
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(fUseYCutHQ && y>fYMinHQ && y<fYMaxHQ) inYcut=kTRUE;
+ if(!fUseYCutHQ && y>fYMin && y<fYMax) inYcut=kTRUE;
}
if(fTriggerParticle) {
if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
}
//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 ||
return 0;
}
}
-
+
for (i = 0; i < np; i++) {
Int_t trackIt = 0;
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;
}
// Builds the event header, to be called after each event
if (fHeader) delete fHeader;
fHeader = new AliGenPythiaEventHeader("Pythia");
+ fHeader->SetTitle(GetTitle());
//
// Event type
if(fProcDiff>0){
// 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)*fPythia->GetPARI(10));
+ // PARI(7) is 1 or -1, for weighted generation with accept/reject, e.g. POWHEG
+ // PARI(10) is a weight associated with reweighted generation, using Pyevwt
//
// 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;
+}
+