// --- ROOT system ---
-#include <TLorentzVector.h>
#include <TObjArray.h>
// --- AliRoot system ---
#include "AliCaloTrackReader.h"
#include "AliMixedEvent.h"
#include "AliCaloPID.h"
+#include "AliLog.h"
ClassImp(AliIsolationCut)
fConeSize(0.),
fPtThreshold(0.),
fPtThresholdMax(10000.),
-fSumPtThreshold(0.),
-fPtFraction(0.),
+fSumPtThreshold(0.),
+fSumPtThresholdMax(10000.),
+fPtFraction(0.),
fICMethod(0),
fPartInCone(0),
-fDebug(-1),
-fFracIsThresh(1)
+fDebug(0),
+fFracIsThresh(1),
+fMomentum(),
+fTrackVector()
{
//default ctor
if(coneA > excessA) return coneA / (coneA-excessA);
else
{
- printf("AliIsolationCut::CalculateExcessAreaFraction() - Please Check : Excess Track %2.3f, coneA %2.2f, excessA %2.2f, angle %2.2f,factor %2.2f\n",
- excess,coneA, excessA, angle*TMath::RadToDeg(), coneA / (coneA-excessA));
+ AliWarning(Form("Please Check : Excess Track %2.3f, coneA %2.2f, excessA %2.2f, angle %2.2f,factor %2.2f",
+ excess,coneA, excessA, angle*TMath::RadToDeg(), coneA / (coneA-excessA)));
return 1;
}
}
if(phiC<0) phiC+=TMath::TwoPi();
Float_t etaC = pCandidate->Eta() ;
- if(pCandidate->GetDetector()=="EMCAL")
+ if(pCandidate->GetDetectorTag() == AliCaloTrackReader::kEMCAL)
{
AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
AliCalorimeterUtils *cu = reader->GetCaloUtils();
{
//Get absolute (col,row) of candidate
Int_t iEta=-1, iPhi=-1, iRCU = -1;
- Int_t nSupMod = cu->GetModuleNumberCellIndexes(absId, pCandidate->GetDetector(), iEta, iPhi, iRCU);
+ Int_t nSupMod = cu->GetModuleNumberCellIndexes(absId, pCandidate->GetDetectorTag(), iEta, iPhi, iRCU);
Int_t colC = iEta;
if (nSupMod % 2) colC = AliEMCALGeoParams::fgkEMCALCols + iEta ;
}
}//end of cells loop
}
-
- else if(fDebug>0) printf("cluster with bad (eta,phi) in EMCal for energy density calculation\n");
+ else AliWarning("Cluster with bad (eta,phi) in EMCal for energy density calculation");
if (coneCells > 0.)
{
if(phiC<0) phiC+=TMath::TwoPi();
Float_t etaC = pCandidate->Eta() ;
- if(pCandidate->GetDetector()=="EMCAL")
+ if(pCandidate->GetDetectorTag() == AliCaloTrackReader::kEMCAL)
{
AliEMCALGeometry* eGeom = AliEMCALGeometry::GetInstance();
AliCalorimeterUtils *cu = reader->GetCaloUtils();
{
//Get absolute (col,row) of candidate
Int_t iEta=-1, iPhi=-1, iRCU = -1;
- Int_t nSupMod = cu->GetModuleNumberCellIndexes(absId, pCandidate->GetDetector(),
+ Int_t nSupMod = cu->GetModuleNumberCellIndexes(absId, pCandidate->GetDetectorTag(),
iEta, iPhi, iRCU);
Int_t colC = iEta;
}
}//end of cells loop
}
-
- else if(fDebug > 0) printf("cluster with bad (eta,phi) in EMCal for energy density coeff calculation\n");
+ else AliWarning("Cluster with bad (eta,phi) in EMCal for energy density coeff calculation");
if (coneCells > 0.)
{
parList+=onePar ;
snprintf(onePar,buffersize,"fConeSize: (isolation cone size) %1.2f\n",fConeSize) ;
parList+=onePar ;
- snprintf(onePar,buffersize,"fPtThreshold =%1.2f (isolation pt threshold) \n",fPtThreshold) ;
+ snprintf(onePar,buffersize,"fPtThreshold >%2.2f;<%2.2f (isolation pt threshold) \n",fPtThreshold,fPtThresholdMax) ;
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"fSumPtThreshold >%2.2f;<%2.2f (isolation sum pt threshold) \n",fSumPtThreshold,fSumPtThresholdMax) ;
parList+=onePar ;
- snprintf(onePar,buffersize,"fPtFraction=%1.2f (isolation pt threshold fraction ) \n",fPtFraction) ;
+ snprintf(onePar,buffersize,"fPtFraction=%2.2f (isolation pt threshold fraction) \n",fPtFraction) ;
parList+=onePar ;
snprintf(onePar,buffersize,"fICMethod=%d (isolation cut case) \n",fICMethod) ;
parList+=onePar ;
snprintf(onePar,buffersize,"fPartInCone=%d \n",fPartInCone) ;
parList+=onePar ;
- snprintf(onePar,buffersize,"fFracIsThresh=%i \n",fFracIsThresh) ;
+ snprintf(onePar,buffersize,"fFracIsThresh=%i \n",fFracIsThresh) ;
parList+=onePar ;
return parList;
//Initialize the parameters of the analysis.
fConeSize = 0.4 ;
- fPtThreshold = 1. ;
+ fPtThreshold = 0.5 ;
fPtThresholdMax = 10000. ;
- fSumPtThreshold = 0.5 ;
- fPtFraction = 0.1 ;
- fPartInCone = kOnlyCharged;
- fICMethod = kSumPtFracIC; // 0 pt threshol method, 1 cone pt sum method
+ fSumPtThreshold = 1.0 ;
+ fSumPtThresholdMax = 10000. ;
+ fPtFraction = 0.1 ;
+ fPartInCone = kNeutralAndCharged;
+ fICMethod = kSumPtIC; // 0 pt threshol method, 1 cone pt sum method
fFracIsThresh = 1;
}
TString aodArrayRefName,
Int_t & n,
Int_t & nfrac,
- Float_t & coneptsum,
- Bool_t & isolated) const
+ Float_t & coneptsum, Float_t & ptLead,
+ Bool_t & isolated)
{
//Search in cone around a candidate particle if it is isolated
Float_t ptC = pCandidate->Pt() ;
if(phiC<0) phiC+=TMath::TwoPi();
Float_t etaC = pCandidate->Eta() ;
- Float_t pt = -100. ;
- Float_t eta = -100. ;
- Float_t phi = -100. ;
- Float_t rad = -100. ;
+ Float_t pt = -100. ;
+ Float_t eta = -100. ;
+ Float_t phi = -100. ;
+ Float_t rad = -100. ;
Float_t coneptsumCluster = 0;
Float_t coneptsumTrack = 0;
nfrac = 0 ;
isolated = kFALSE;
- if(fDebug>0)
- {
- printf("AliIsolationCut::MakeIsolationCut() - Cadidate pT %2.2f, eta %2.2f, phi %2.2f, cone %1.2f, thres %2.2f, Fill AOD? %d",
- pCandidate->Pt(), pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), fConeSize,fPtThreshold,bFillAOD);
- if(plCTS) printf(", nTracks %d" ,plCTS->GetEntriesFast());
- if(plNe) printf(", nClusters %d",plNe ->GetEntriesFast());
-
- printf("\n");
- }
-
+ AliDebug(1,Form("Candidate pT %2.2f, eta %2.2f, phi %2.2f, cone %1.2f, thres %2.2f, Fill AOD? %d",
+ pCandidate->Pt(), pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), fConeSize,fPtThreshold,bFillAOD));
+
//Initialize the array with refrences
TObjArray * refclusters = 0x0;
TObjArray * reftracks = 0x0;
if(plCTS &&
(fPartInCone==kOnlyCharged || fPartInCone==kNeutralAndCharged))
{
- TVector3 p3;
for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ )
{
AliVTrack* track = dynamic_cast<AliVTrack*>(plCTS->At(ipr)) ;
if(track->GetID() == pCandidate->GetTrackLabel(0) || track->GetID() == pCandidate->GetTrackLabel(1) ||
track->GetID() == pCandidate->GetTrackLabel(2) || track->GetID() == pCandidate->GetTrackLabel(3) ) continue ;
- p3.SetXYZ(track->Px(),track->Py(),track->Pz());
- pt = p3.Pt();
- eta = p3.Eta();
- phi = p3.Phi() ;
+ fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
+ pt = fTrackVector.Pt();
+ eta = fTrackVector.Eta();
+ phi = fTrackVector.Phi() ;
}
else
{// Mixed event stored in AliAODPWG4Particles
AliAODPWG4Particle * trackmix = dynamic_cast<AliAODPWG4Particle*>(plCTS->At(ipr)) ;
if(!trackmix)
{
- printf("AliIsolationCut::MakeIsolationCut() - Wrong track data type, continue\n");
+ AliWarning("Wrong track data type, continue");
continue;
}
eta = trackmix->Eta();
phi = trackmix->Phi() ;
}
-
+
if( phi < 0 ) phi+=TMath::TwoPi();
rad = Radius(etaC, phiC, eta, phi);
// Only loop the particle at the same side of candidate
if(TMath::Abs(phi-phiC) > TMath::PiOver2()) continue ;
- // If at the same side has particle larger than candidate,
- // then candidate can not be the leading, skip such events
- if(pt > ptC)
- {
- n = -1;
- nfrac = -1;
- coneptsumTrack = -1;
- isolated = kFALSE;
-
- pCandidate->SetLeadingParticle(kFALSE);
-
- if(bFillAOD && reftracks)
- {
- reftracks->Clear();
- delete reftracks;
- }
-
- return ;
- }
+// // If at the same side has particle larger than candidate,
+// // then candidate can not be the leading, skip such events
+// if(pt > ptC)
+// {
+// n = -1;
+// nfrac = -1;
+// coneptsumTrack = -1;
+// isolated = kFALSE;
+//
+// pCandidate->SetLeadingParticle(kFALSE);
+//
+// if(bFillAOD && reftracks)
+// {
+// reftracks->Clear();
+// delete reftracks;
+// }
+//
+// return ;
+// }
- //Check if there is any particle inside cone with pt larger than fPtThreshold
+ // // Check if there is any particle inside cone with pt larger than fPtThreshold
+ // Check if the leading particule inside the cone has a ptLead larger than fPtThreshold
- if( fDebug > 0 )
- printf("\t track %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
+ AliDebug(2,Form("\t Track %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad));
if(rad < fConeSize)
{
- if(fDebug > 0) printf(" - inside candidate cone");
+ AliDebug(2,"Inside candidate cone");
if(bFillAOD)
{
}
coneptsumTrack+=pt;
- if(pt > fPtThreshold && pt < fPtThresholdMax) n++;
- if(pt > fPtFraction*ptC ) nfrac++;
-
+
+ if( ptLead < pt ) ptLead = pt;
+
+// // *Before*, count particles in cone
+// if(pt > fPtThreshold && pt < fPtThresholdMax) n++;
+//
+// //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
+// if(fFracIsThresh)
+// {
+// if( fPtFraction*ptC < fPtThreshold )
+// {
+// if( pt > fPtThreshold ) nfrac++ ;
+// }
+// else
+// {
+// if( pt > fPtFraction*ptC ) nfrac++;
+// }
+// }
+// else
+// {
+// if( pt > fPtFraction*ptC ) nfrac++;
+// }
+
} // Inside cone
- if(fDebug>0) printf("\n");
-
}// charged particle loop
-
}//Tracks
if(plNe &&
(fPartInCone==kOnlyNeutral || fPartInCone==kNeutralAndCharged))
{
- TLorentzVector mom ;
for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ )
{
pid->IsTrackMatched(calo,reader->GetCaloUtils(),reader->GetInputEvent()) ) continue ;
//Assume that come from vertex in straight line
- calo->GetMomentum(mom,reader->GetVertex(evtIndex)) ;
+ calo->GetMomentum(fMomentum,reader->GetVertex(evtIndex)) ;
- pt = mom.Pt() ;
- eta = mom.Eta() ;
- phi = mom.Phi() ;
+ pt = fMomentum.Pt() ;
+ eta = fMomentum.Eta() ;
+ phi = fMomentum.Phi() ;
}
else
{// Mixed event stored in AliAODPWG4Particles
AliAODPWG4Particle * calomix = dynamic_cast<AliAODPWG4Particle*>(plNe->At(ipr)) ;
if(!calomix)
{
- printf("AliIsolationCut::MakeIsolationCut() - Wrong calo data type, continue\n");
+ AliWarning("Wrong calo data type, continue");
continue;
}
// Only loop the particle at the same side of candidate
if(TMath::Abs(phi-phiC)>TMath::PiOver2()) continue ;
- // If at the same side has particle larger than candidate,
- // then candidate can not be the leading, skip such events
- if(pt > ptC)
- {
- n = -1;
- nfrac = -1;
- coneptsumCluster = -1;
- isolated = kFALSE;
-
- pCandidate->SetLeadingParticle(kFALSE);
-
- if(bFillAOD)
- {
- if(reftracks)
- {
- reftracks ->Clear();
- delete reftracks;
- }
-
- if(refclusters)
- {
- refclusters->Clear();
- delete refclusters;
- }
- }
- return ;
- }
+// // If at the same side has particle larger than candidate,
+// // then candidate can not be the leading, skip such events
+// if(pt > ptC)
+// {
+// n = -1;
+// nfrac = -1;
+// coneptsumCluster = -1;
+// isolated = kFALSE;
+//
+// pCandidate->SetLeadingParticle(kFALSE);
+//
+// if(bFillAOD)
+// {
+// if(reftracks)
+// {
+// reftracks ->Clear();
+// delete reftracks;
+// }
+//
+// if(refclusters)
+// {
+// refclusters->Clear();
+// delete refclusters;
+// }
+// }
+// return ;
+// }
//Check if there is any particle inside cone with pt larger than fPtThreshold
- if(fDebug > 0 )
- printf("\t cluster %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad);
+ AliDebug(2,Form("\t Cluster %d, pT %2.2f, eta %1.2f, phi %2.2f, R candidate %2.2f", ipr,pt,eta,phi,rad));
if(rad < fConeSize)
{
- if(fDebug > 0 ) printf(" - inside candidate cone");
+ AliDebug(2,"Inside candidate cone");
if(bFillAOD)
{
}
coneptsumCluster+=pt;
- if(pt > fPtThreshold && pt < fPtThresholdMax) n++;
- //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
- if(fFracIsThresh)
- {
- if( fPtFraction*ptC<fPtThreshold)
- {
- if(pt>fPtThreshold) nfrac++ ;
- }
- else
- {
- if(pt>fPtFraction*ptC) nfrac++;
- }
- }
- else
- {
- if(pt>fPtFraction*ptC) nfrac++;
- }
+
+ if( ptLead < pt ) ptLead = pt;
+
+// // *Before*, count particles in cone
+// if(pt > fPtThreshold && pt < fPtThresholdMax) n++;
+//
+// //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
+// if(fFracIsThresh)
+// {
+// if( fPtFraction*ptC < fPtThreshold )
+// {
+// if( pt > fPtThreshold ) nfrac++ ;
+// }
+// else
+// {
+// if( pt > fPtFraction*ptC ) nfrac++;
+// }
+// }
+// else
+// {
+// if( pt > fPtFraction*ptC ) nfrac++;
+// }
}//in cone
- if(fDebug>0) printf("\n");
-
}// neutral particle loop
}//neutrals
-
//Add reference arrays to AOD when filling AODs only
- if(bFillAOD)
+ if(bFillAOD)
{
if(refclusters) pCandidate->AddObjArray(refclusters);
if(reftracks) pCandidate->AddObjArray(reftracks);
}
- coneptsum = coneptsumCluster+coneptsumTrack;
+ coneptsum = coneptsumCluster + coneptsumTrack;
+
+ // *Now*, just check the leading particle in the cone if the threshold is passed
+ if(ptLead > fPtThreshold && ptLead < fPtThresholdMax) n = 1;
+
+ //if fPtFraction*ptC<fPtThreshold then consider the fPtThreshold directly
+ if(fFracIsThresh)
+ {
+ if( fPtFraction*ptC < fPtThreshold )
+ {
+ if( ptLead > fPtThreshold ) nfrac = 1 ;
+ }
+ else
+ {
+ if( ptLead > fPtFraction*ptC ) nfrac = 1;
+ }
+ }
+ else
+ {
+ if( ptLead > fPtFraction*ptC ) nfrac = 1;
+ }
+
+ //-------------------------------------------------------------------
+ //Check isolation, depending on selected isolation criteria requested
- //Check isolation, depending on selected isolation criteria
if( fICMethod == kPtThresIC)
{
- if(n==0) isolated = kTRUE ;
+ if( n == 0 ) isolated = kTRUE ;
+
+ AliDebug(1,Form("pT Cand %2.2f, pT Lead %2.2f, %2.2f<pT Lead< %2.2f, isolated %d",
+ ptC,ptLead,fPtThreshold,fPtThresholdMax,isolated));
}
- else if( fICMethod == kSumPtIC)
+ else if( fICMethod == kSumPtIC )
{
- if(coneptsum < fSumPtThreshold)
- isolated = kTRUE ;
+ if( coneptsum > fSumPtThreshold &&
+ coneptsum < fSumPtThresholdMax )
+ isolated = kFALSE ;
+ else
+ isolated = kTRUE ;
+
+ AliDebug(1,Form("pT Cand %2.2f, SumPt %2.2f, %2.2f<Sum pT< %2.2f, isolated %d",
+ ptC,ptLead,fSumPtThreshold,fSumPtThresholdMax,isolated));
}
- else if( fICMethod == kPtFracIC)
+ else if( fICMethod == kPtFracIC )
{
- if(nfrac==0) isolated = kTRUE ;
+ if(nfrac == 0 ) isolated = kTRUE ;
}
- else if( fICMethod == kSumPtFracIC)
+ else if( fICMethod == kSumPtFracIC )
{
//when the fPtFraction*ptC < fSumPtThreshold then consider the later case
- // printf("photon analysis IsDataMC() ?%i\n",IsDataMC());
- if(fFracIsThresh )
+ // printf("photon analysis IsDataMC() ?%i\n",IsDataMC());
+ if( fFracIsThresh )
{
- if( fPtFraction*ptC < fSumPtThreshold && coneptsum < fSumPtThreshold) isolated = kTRUE ;
- if( fPtFraction*ptC > fSumPtThreshold && coneptsum < fPtFraction*ptC) isolated = kTRUE ;
+ if( fPtFraction*ptC < fSumPtThreshold && coneptsum < fSumPtThreshold ) isolated = kTRUE ;
+ if( fPtFraction*ptC > fSumPtThreshold && coneptsum < fPtFraction*ptC ) isolated = kTRUE ;
}
else
{
- if(coneptsum < fPtFraction*ptC) isolated = kTRUE ;
+ if( coneptsum < fPtFraction*ptC ) isolated = kTRUE ;
}
}
- else if( fICMethod == kSumDensityIC)
+ else if( fICMethod == kSumDensityIC )
{
// Get good cell density (number of active cells over all cells in cone)
// and correct energy in cone
Float_t cellDensity = GetCellDensity(pCandidate,reader);
- if(coneptsum < fSumPtThreshold*cellDensity)
+ if( coneptsum < fSumPtThreshold*cellDensity )
isolated = kTRUE;
}
- else if( fICMethod == kSumBkgSubIC)
+ else if( fICMethod == kSumBkgSubIC )
{
Double_t coneptsumBkg = 0.;
Float_t etaBandPtSumTrackNorm = 0;
coneptsum = coneptsumCluster+coneptsumTrack;
coneptsum -= coneptsumBkg;
- if(coneptsum < fSumPtThreshold)
- isolated = kTRUE ;
- }
+
+ if( coneptsum > fSumPtThreshold && coneptsum < fSumPtThresholdMax )
+ isolated = kFALSE ;
+ else
+ isolated = kTRUE ;
+
+ }
}
printf("IC method = %d\n", fICMethod ) ;
printf("Cone Size = %1.2f\n", fConeSize ) ;
- printf("pT threshold = %2.1f\n", fPtThreshold) ;
+ printf("pT threshold = >%2.1f;<%2.1f\n", fPtThreshold , fPtThresholdMax) ;
+ printf("Sum pT threshold = >%2.1f;<%2.1f\n", fSumPtThreshold,fSumPtThresholdMax) ;
printf("pT fraction = %3.1f\n", fPtFraction ) ;
printf("particle type in cone = %d\n", fPartInCone ) ;
printf("using fraction for high pt leading instead of frac ? %i\n",fFracIsThresh);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff