#include <TMath.h>
#include <TString.h>
#include <TList.h>
#include "AliCaloPID.h"
#include "AliAODCaloCluster.h"
#include "AliVCaloCells.h"
#include "AliVTrack.h"
#include "AliAODPWG4Particle.h"
#include "AliCalorimeterUtils.h"
#include "AliVEvent.h"
#include "AliLog.h"
#include "AliEMCALPIDUtils.h"
ClassImp(AliCaloPID)
AliCaloPID::AliCaloPID() :
TObject(), fDebug(-1), fParticleFlux(kLow),
fEMCALPIDUtils(), fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
fEMCALElectronWeight(0.), fEMCALChargeWeight(0.), fEMCALNeutralWeight(0.),
fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) , fPHOSNeutralWeight(0.),
fPHOSWeightFormula(0), fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
fPHOSPhotonWeightFormulaExpression(""),
fPHOSPi0WeightFormulaExpression(""),
fEMCALL0CutMax(100.), fEMCALL0CutMin(0),
fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.),
fTOFCut(0.),
fPHOSDispersionCut(1000), fPHOSRCut(1000),
fUseSimpleMassCut(kFALSE),
fUseSimpleM02Cut(kFALSE),
fUseSplitAsyCut(kFALSE),
fUseSplitSSCut(kTRUE),
fSplitM02MaxCut(0), fSplitM02MinCut(0), fSplitMinNCells(0),
fMassEtaMin(0), fMassEtaMax(0),
fMassPi0Min(0), fMassPi0Max(0),
fMassPhoMin(0), fMassPhoMax(0),
fM02MaxParamShiftNLMN(0),
fSplitWidthSigma(0), fMassShiftHighECell(0)
{
InitParameters();
}
AliCaloPID::AliCaloPID(Int_t flux) :
TObject(), fDebug(-1), fParticleFlux(flux),
fEMCALPIDUtils(), fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
fEMCALElectronWeight(0.), fEMCALChargeWeight(0.), fEMCALNeutralWeight(0.),
fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) , fPHOSNeutralWeight(0.),
fPHOSWeightFormula(0), fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
fPHOSPhotonWeightFormulaExpression(""),
fPHOSPi0WeightFormulaExpression(""),
fEMCALL0CutMax(100.), fEMCALL0CutMin(0),
fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.),
fTOFCut(0.),
fPHOSDispersionCut(1000), fPHOSRCut(1000),
fUseSimpleMassCut(kFALSE),
fUseSimpleM02Cut(kFALSE),
fUseSplitAsyCut(kFALSE),
fUseSplitSSCut(kTRUE),
fSplitM02MaxCut(0), fSplitM02MinCut(0), fSplitMinNCells(0),
fMassEtaMin(0), fMassEtaMax(0),
fMassPi0Min(0), fMassPi0Max(0),
fMassPhoMin(0), fMassPhoMax(0),
fM02MaxParamShiftNLMN(0),
fSplitWidthSigma(0), fMassShiftHighECell(0)
{
InitParameters();
}
AliCaloPID::AliCaloPID(const TNamed * emcalpid) :
TObject(), fDebug(-1), fParticleFlux(kLow),
fEMCALPIDUtils((AliEMCALPIDUtils*)emcalpid),
fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
fEMCALElectronWeight(0.), fEMCALChargeWeight(0.), fEMCALNeutralWeight(0.),
fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) , fPHOSNeutralWeight(0.),
fPHOSWeightFormula(0), fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
fPHOSPhotonWeightFormulaExpression(""),
fPHOSPi0WeightFormulaExpression(""),
fEMCALL0CutMax(100.), fEMCALL0CutMin(0),
fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.),
fTOFCut(0.),
fPHOSDispersionCut(1000), fPHOSRCut(1000),
fUseSimpleMassCut(kFALSE),
fUseSimpleM02Cut(kFALSE),
fUseSplitAsyCut(kFALSE),
fUseSplitSSCut(kTRUE),
fSplitM02MaxCut(0), fSplitM02MinCut(0), fSplitMinNCells(0),
fMassEtaMin(0), fMassEtaMax(0),
fMassPi0Min(0), fMassPi0Max(0),
fMassPhoMin(0), fMassPhoMax(0),
fM02MaxParamShiftNLMN(0),
fSplitWidthSigma(0), fMassShiftHighECell(0)
{
InitParameters();
}
AliCaloPID::~AliCaloPID()
{
delete fPHOSPhotonWeightFormula ;
delete fPHOSPi0WeightFormula ;
delete fEMCALPIDUtils ;
}
void AliCaloPID::InitParameters()
{
fEMCALPhotonWeight = 0.6 ;
fEMCALPi0Weight = 0.6 ;
fEMCALElectronWeight = 0.6 ;
fEMCALChargeWeight = 0.6 ;
fEMCALNeutralWeight = 0.6 ;
fPHOSPhotonWeight = 0.6 ;
fPHOSPi0Weight = 0.6 ;
fPHOSElectronWeight = 0.6 ;
fPHOSChargeWeight = 0.6 ;
fPHOSNeutralWeight = 0.6 ;
fPHOSWeightFormula = kFALSE;
fPHOSPhotonWeightFormulaExpression = "0.98*(x<40)+ 0.68*(x>=100)+(x>=40 && x<100)*(0.98+x*(6e-3)-x*x*(2e-04)+x*x*x*(1.1e-06))";
fPHOSPi0WeightFormulaExpression = "0.98*(x<65)+ 0.915*(x>=100)+(x>=65 && x-x*(1.95e-3)-x*x*(4.31e-05)+x*x*x*(3.61e-07))" ;
if(fRecalculateBayesian)
{
if(fParticleFlux == kLow)
{
AliInfo("SetLOWFluxParam");
fEMCALPIDUtils->SetLowFluxParam() ;
}
else if (fParticleFlux == kHigh)
{
AliInfo("SetHighFluxParam");
fEMCALPIDUtils->SetHighFluxParam() ;
}
}
fEMCALL0CutMax = 0.3 ;
fEMCALL0CutMin = 0.01;
fEMCALDPhiCut = 0.05;
fEMCALDEtaCut = 0.025;
fTOFCut = 1.e-6;
fPHOSRCut = 2. ;
fPHOSDispersionCut = 2.5;
fSplitM02MinCut = 0.3 ;
fSplitM02MaxCut = 5 ;
fSplitMinNCells = 4 ;
fMassEtaMin = 0.4;
fMassEtaMax = 0.6;
fMassPi0Min = 0.11;
fMassPi0Max = 0.18;
fMassPhoMin = 0.0;
fMassPhoMax = 0.08;
fMassPi0Param[0][0] = 0 ;
fMassPi0Param[0][1] = 0 ;
fMassPi0Param[0][2] = 0 ;
fMassPi0Param[0][3] = 0.044 ;
fMassPi0Param[0][4] = 0.0049;
fMassPi0Param[0][5] = 0.070 ;
fMassPi0Param[1][0] = 0.115 ;
fMassPi0Param[1][1] = 0.00096;
fMassPi0Param[1][2] = 21 ;
fMassPi0Param[1][3] = 0.10 ;
fMassPi0Param[1][4] = 0.0017;
fMassPi0Param[1][5] = 0.070 ;
fWidthPi0Param[0][0] = 0.012 ;
fWidthPi0Param[0][1] = 0.0 ;
fWidthPi0Param[0][2] = 19 ;
fWidthPi0Param[0][3] = 0.0012;
fWidthPi0Param[0][4] = 0.0006;
fWidthPi0Param[0][5] = 0.0 ;
fWidthPi0Param[1][0] = 0.009 ;
fWidthPi0Param[1][1] = 0.000 ;
fWidthPi0Param[1][2] = 10 ;
fWidthPi0Param[1][3] = 0.0023 ;
fWidthPi0Param[1][4] = 0.00067;
fWidthPi0Param[1][5] = 0.000 ;
fMassShiftHighECell = 0;
fM02MinParam[0][0] = 2.135 ;
fM02MinParam[0][1] =-0.245 ;
fM02MinParam[0][2] = 0.0 ;
fM02MinParam[0][3] = 0.0 ;
fM02MinParam[0][4] = 0.0 ;
fM02MinParam[1][0] = 2.135 ;
fM02MinParam[1][1] =-0.245 ;
fM02MinParam[1][2] = 0.0 ;
fM02MinParam[1][3] = 0.0 ;
fM02MinParam[1][4] = 0.0 ;
fM02MaxParam[0][0] = 0.0662 ;
fM02MaxParam[0][1] =-0.0201 ;
fM02MaxParam[0][2] =-0.0955 ;
fM02MaxParam[0][3] = 0.00186;
fM02MaxParam[0][4] = 9.91 ;
fM02MaxParam[1][0] = 0.353 ;
fM02MaxParam[1][1] =-0.0264 ;
fM02MaxParam[1][2] =-0.524 ;
fM02MaxParam[1][3] = 0.00559;
fM02MaxParam[1][4] = 21.9 ;
fM02MaxParamShiftNLMN = 0.75;
fAsyMinParam[0][0] = 0.96 ;
fAsyMinParam[0][1] = 0 ;
fAsyMinParam[0][2] =-879 ;
fAsyMinParam[0][3] = 0.96 ;
fAsyMinParam[1][0] = 0.95 ;
fAsyMinParam[1][1] = 0.0015;
fAsyMinParam[1][2] =-233 ;
fAsyMinParam[1][3] = 1.0 ;
fSplitEFracMin[0] = 0.0 ;
fSplitEFracMin[1] = 0.0 ;
fSplitEFracMin[2] = 0.0 ;
fSubClusterEMin[0] = 0.0;
fSubClusterEMin[1] = 0.0;
fSubClusterEMin[2] = 0.0;
fSplitWidthSigma = 3. ;
}
Bool_t AliCaloPID::IsInPi0SplitAsymmetryRange(Float_t energy, Float_t asy, Int_t nlm) const
{
if(!fUseSplitAsyCut) return kTRUE ;
Float_t abasy = TMath::Abs(asy);
Int_t inlm = nlm-1;
if(nlm > 2) inlm=1;
Float_t cut = fAsyMinParam[inlm][0] + fAsyMinParam[inlm][1]*energy + fAsyMinParam[inlm][2]/energy/energy/energy ;
if( cut > fAsyMinParam[inlm][3] ) cut = fAsyMinParam[inlm][3];
if(abasy < cut) return kTRUE;
else return kFALSE;
}
Bool_t AliCaloPID::IsInPi0SplitMassRange(Float_t energy, Float_t mass, Int_t nlm) const
{
if(fUseSimpleMassCut)
{
if(mass < fMassPi0Max && mass > fMassPi0Min) return kTRUE;
else return kFALSE;
}
Int_t inlm = nlm-1;
if(nlm > 2) inlm=1;
Float_t meanMass = energy * fMassPi0Param[inlm][1] + fMassPi0Param[inlm][0];
if(energy > fMassPi0Param[inlm][2]) meanMass = energy * fMassPi0Param[inlm][4] + fMassPi0Param[inlm][3];
meanMass -= fMassShiftHighECell;
Float_t width = 0.009;
if (energy > 8 && energy < fWidthPi0Param[inlm][2])
width = energy * fWidthPi0Param[inlm][1] + fWidthPi0Param[inlm][0];
else if( energy > fWidthPi0Param[inlm][2])
width = energy * energy * fWidthPi0Param[inlm][5] + energy * fWidthPi0Param[inlm][4] + fWidthPi0Param[inlm][3];
Float_t minMass = meanMass-fSplitWidthSigma*width;
Float_t maxMass = meanMass+fSplitWidthSigma*width;
if(energy < 10 && minMass < fMassPi0Param[inlm][5] ) minMass = fMassPi0Param[inlm][5];
if(mass < maxMass && mass > minMass) return kTRUE;
else return kFALSE;
}
Bool_t AliCaloPID::IsInM02Range(Float_t m02) const
{
Float_t minCut = fSplitM02MinCut;
Float_t maxCut = fSplitM02MaxCut;
if(m02 < maxCut && m02 > minCut) return kTRUE;
else return kFALSE;
}
Bool_t AliCaloPID::IsInPi0M02Range(Float_t energy, Float_t m02, Int_t nlm) const
{
if(!fUseSplitSSCut) return kTRUE ;
if(!IsInM02Range(m02)) return kFALSE ;
else if(!fUseSimpleM02Cut)
{
Int_t inlm = nlm-1;
if(nlm > 2) inlm=1;
Float_t minCut = fSplitM02MinCut;
Float_t maxCut = fSplitM02MaxCut;
if(energy > 1) minCut = TMath::Exp( fM02MinParam[inlm][0] + fM02MinParam[inlm][1]*energy ) +
fM02MinParam[inlm][2] + fM02MinParam[inlm][3]*energy + fM02MinParam[inlm][4]/energy;
if(energy > 1) maxCut = TMath::Exp( fM02MaxParam[inlm][0] + fM02MaxParam[inlm][1]*energy ) +
fM02MaxParam[inlm][2] + fM02MaxParam[inlm][3]*energy + fM02MaxParam[inlm][4]/energy;
if( minCut < fSplitM02MinCut) minCut = fSplitM02MinCut;
if( maxCut > fSplitM02MaxCut) maxCut = fSplitM02MaxCut;
if( nlm > 2 ) maxCut+=fM02MaxParamShiftNLMN;
if(m02 < maxCut && m02 > minCut) return kTRUE;
else return kFALSE;
}
else return kTRUE;
}
Bool_t AliCaloPID::IsInEtaM02Range(Float_t energy, Float_t m02, Int_t nlm) const
{
if(!fUseSplitSSCut) return kTRUE ;
if(!IsInM02Range(m02)) return kFALSE ;
else if(!fUseSimpleM02Cut)
{
Int_t inlm = nlm-1;
if(nlm > 2) inlm=1;
Float_t minCut = fSplitM02MinCut;
Float_t maxCut = fSplitM02MaxCut;
Float_t shiftE = energy-20;
if(nlm==1) shiftE=energy-28;
if(shiftE > 1) minCut = TMath::Exp( fM02MinParam[inlm][0] + fM02MinParam[inlm][1]*shiftE ) +
fM02MinParam[inlm][2] + fM02MinParam[inlm][3]*shiftE + fM02MinParam[inlm][4]/shiftE;
if( minCut < fSplitM02MinCut) minCut = fSplitM02MinCut;
shiftE = energy+20;
if(shiftE > 1) maxCut = 1 + TMath::Exp( fM02MaxParam[inlm][0] + fM02MaxParam[inlm][1]*shiftE ) +
fM02MaxParam[inlm][2] + fM02MaxParam[inlm][3]*shiftE + fM02MaxParam[inlm][4]/shiftE;
if( maxCut > fSplitM02MaxCut) maxCut = fSplitM02MaxCut;
if( nlm > 2 ) maxCut+=fM02MaxParamShiftNLMN;
if(m02 < maxCut && m02 > minCut) return kTRUE;
else return kFALSE;
}
else return kTRUE;
}
Bool_t AliCaloPID::IsInConM02Range(Float_t energy, Float_t m02, Int_t nlm) const
{
if(!fUseSplitSSCut) return kTRUE ;
Float_t minCut = 0.1;
Float_t maxCut = fSplitM02MinCut;
if(!fUseSimpleM02Cut)
{
Int_t inlm = nlm-1;
if(nlm > 2) inlm=1;
if(energy > 1) maxCut = TMath::Exp( fM02MinParam[inlm][0] + fM02MinParam[inlm][1]*energy ) +
fM02MinParam[inlm][2] + fM02MinParam[inlm][3]*energy + fM02MinParam[inlm][4]/energy;
if( maxCut < fSplitM02MinCut) maxCut = fSplitM02MinCut;
}
if(m02 < maxCut && m02 > minCut) return kTRUE;
else return kFALSE;
}
AliEMCALPIDUtils *AliCaloPID::GetEMCALPIDUtils()
{
if(!fEMCALPIDUtils) fEMCALPIDUtils = new AliEMCALPIDUtils ;
return fEMCALPIDUtils ;
}
Int_t AliCaloPID::GetIdentifiedParticleType(AliVCluster * cluster)
{
Float_t energy = cluster->E();
Float_t lambda0 = cluster->GetM02();
Float_t lambda1 = cluster->GetM20();
if(fUseBayesianWeights)
{
Double_t weights[AliPID::kSPECIESCN];
if(cluster->IsEMCAL() && fRecalculateBayesian)
{
fEMCALPIDUtils->ComputePID(energy, lambda0);
for(Int_t i = 0; i < AliPID::kSPECIESCN; i++) weights[i] = fEMCALPIDUtils->GetPIDFinal(i);
}
else
{
for(Int_t i = 0; i < AliPID::kSPECIESCN; i++) weights[i] = cluster->GetPID()[i];
}
if(fDebug > 0) PrintClusterPIDWeights(weights);
return GetIdentifiedParticleTypeFromBayesWeights(cluster->IsEMCAL(), weights, energy);
}
AliDebug(1,Form("EMCAL %d?, E %3.2f, l0 %3.2f, l1 %3.2f, disp %3.2f, tof %1.11f, distCPV %3.2f, distToBC %1.1f, NMax %d",
cluster->IsEMCAL(),energy,lambda0,cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(),
cluster->GetEmcCpvDistance(), cluster->GetDistanceToBadChannel(),cluster->GetNExMax()));
if(cluster->IsEMCAL())
{
AliDebug(1,Form("EMCAL SS %f <%f < %f?",fEMCALL0CutMin, lambda0, fEMCALL0CutMax));
if(lambda0 < fEMCALL0CutMax && lambda0 > fEMCALL0CutMin) return kPhoton ;
else return kNeutralUnknown ;
}
else
{
if(TestPHOSDispersion(energy,lambda0,lambda1) < fPHOSDispersionCut) return kPhoton;
else return kNeutralUnknown;
}
}
Int_t AliCaloPID::GetIdentifiedParticleTypeFromBayesWeights(Bool_t isEMCAL, Double_t * pid, Float_t energy)
{
if(!pid)
{
AliFatal("pid pointer not initialized!!!");
return kNeutralUnknown;
}
Float_t wPh = fPHOSPhotonWeight ;
Float_t wPi0 = fPHOSPi0Weight ;
Float_t wE = fPHOSElectronWeight ;
Float_t wCh = fPHOSChargeWeight ;
Float_t wNe = fPHOSNeutralWeight ;
if(!isEMCAL && fPHOSWeightFormula){
wPh = GetPHOSPhotonWeightFormula()->Eval(energy) ;
wPi0 = GetPHOSPi0WeightFormula() ->Eval(energy);
}
else
{
wPh = fEMCALPhotonWeight ;
wPi0 = fEMCALPi0Weight ;
wE = fEMCALElectronWeight ;
wCh = fEMCALChargeWeight ;
wNe = fEMCALNeutralWeight ;
}
if(fDebug > 0) PrintClusterPIDWeights(pid);
Int_t pdg = kNeutralUnknown ;
Float_t chargedHadronWeight = pid[AliVCluster::kProton]+pid[AliVCluster::kKaon]+
pid[AliVCluster::kPion]+pid[AliVCluster::kMuon];
Float_t neutralHadronWeight = pid[AliVCluster::kNeutron]+pid[AliVCluster::kKaon0];
Float_t allChargedWeight = pid[AliVCluster::kElectron]+pid[AliVCluster::kEleCon]+ chargedHadronWeight;
Float_t allNeutralWeight = pid[AliVCluster::kPhoton]+pid[AliVCluster::kPi0]+ neutralHadronWeight;
if(!isEMCAL)
{
if(pid[AliVCluster::kPhoton] > wPh) pdg = kPhoton ;
else if(pid[AliVCluster::kPi0] > wPi0) pdg = kPi0 ;
else if(pid[AliVCluster::kElectron] > wE) pdg = kElectron ;
else if(pid[AliVCluster::kEleCon] > wE) pdg = kEleCon ;
else if(chargedHadronWeight > wCh) pdg = kChargedHadron ;
else if(neutralHadronWeight > wNe) pdg = kNeutralHadron ;
else if(allChargedWeight > allNeutralWeight)
pdg = kChargedUnknown ;
else
pdg = kNeutralUnknown ;
}
else
{
if(pid[AliVCluster::kPhoton] > wPh) pdg = kPhoton ;
else if(pid[AliVCluster::kElectron] > wE) pdg = kElectron ;
else if(pid[AliVCluster::kPhoton]+pid[AliVCluster::kElectron] > wPh) pdg = kPhoton ;
else if(pid[AliVCluster::kPi0] > wPi0) pdg = kPi0 ;
else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;
else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ;
else pdg = kNeutralUnknown ;
}
AliDebug(1,Form("Final Pdg: %d, cluster energy %2.2f", pdg,energy));
return pdg ;
}
Int_t AliCaloPID::GetIdentifiedParticleTypeFromClusterSplitting(AliVCluster* cluster,
AliVCaloCells* cells,
AliCalorimeterUtils * caloutils,
Double_t vertex[3],
Int_t & nMax,
Double_t & mass, Double_t & angle,
TLorentzVector & l1, TLorentzVector & l2,
Int_t & absId1, Int_t & absId2,
Float_t & distbad1, Float_t & distbad2,
Bool_t & fidcut1, Bool_t & fidcut2 ) const
{
Float_t eClus = cluster->E();
Float_t m02 = cluster->GetM02();
const Int_t nc = cluster->GetNCells();
Int_t absIdList[nc];
Float_t maxEList [nc];
mass = -1.;
angle = -1.;
if ( nc < fSplitMinNCells) return kNeutralUnknown ;
AliDebug(2,"\t pass nCells cut");
nMax = caloutils->GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ;
AliDebug(1,Form("Cluster : E %1.1f, M02 %1.2f, NLM %d, N Cells %d",eClus,m02,nMax,nc));
Int_t calorimeter = AliCalorimeterUtils::kEMCAL;
if(cluster->IsPHOS()) calorimeter = AliCalorimeterUtils::kPHOS;
if ( nMax == 2 )
{
absId1 = absIdList[0];
absId2 = absIdList[1];
Float_t en1 = cells->GetCellAmplitude(absId1);
caloutils->RecalibrateCellAmplitude(en1,calorimeter,absId1);
Float_t en2 = cells->GetCellAmplitude(absId2);
caloutils->RecalibrateCellAmplitude(en2,calorimeter,absId2);
if(en1 < en2)
{
absId2 = absIdList[0];
absId1 = absIdList[1];
}
}
else if( nMax == 1 )
{
absId1 = absIdList[0];
Float_t enmax = 0 ;
for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++)
{
Int_t absId = cluster->GetCellsAbsId()[iDigit];
if( absId == absId1 ) continue ;
Float_t endig = cells->GetCellAmplitude(absId);
caloutils->RecalibrateCellAmplitude(endig,calorimeter,absId);
if(endig > enmax)
{
enmax = endig ;
absId2 = absId ;
}
}
}
else
{
Float_t enmax = 0 ;
for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++)
{
Float_t endig = maxEList[iDigit];
if(endig > enmax)
{
enmax = endig ;
absId1 = absIdList[iDigit];
}
}
Float_t enmax2 = 0;
for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++)
{
if(absIdList[iDigit]==absId1) continue;
Float_t endig = maxEList[iDigit];
if(endig > enmax2)
{
enmax2 = endig ;
absId2 = absIdList[iDigit];
}
}
}
if(absId2<0 || absId1<0)
{
AliDebug(1,Form("Bad index for local maxima : N max %d, i1 %d, i2 %d, cluster E %2.2f, ncells %d, m02 %2.2f",
nMax,absId1,absId2,eClus,nc,m02));
return kNeutralUnknown ;
}
AliAODCaloCluster cluster1(0, 0,NULL,0.,NULL,NULL,1,0);
AliAODCaloCluster cluster2(1, 0,NULL,0.,NULL,NULL,1,0);
caloutils->SplitEnergy(absId1,absId2,cluster, cells, &cluster1, &cluster2,nMax);
fidcut1 = caloutils->GetEMCALRecoUtils()->CheckCellFiducialRegion(caloutils->GetEMCALGeometry(), &cluster1,cells);
fidcut2 = caloutils->GetEMCALRecoUtils()->CheckCellFiducialRegion(caloutils->GetEMCALGeometry(), &cluster2,cells);
caloutils->GetEMCALRecoUtils()->RecalculateClusterDistanceToBadChannel(caloutils->GetEMCALGeometry(),cells,&cluster1);
caloutils->GetEMCALRecoUtils()->RecalculateClusterDistanceToBadChannel(caloutils->GetEMCALGeometry(),cells,&cluster2);
distbad1 = cluster1.GetDistanceToBadChannel();
distbad2 = cluster2.GetDistanceToBadChannel();
cluster1.GetMomentum(l1,vertex);
cluster2.GetMomentum(l2,vertex);
mass = (l1+l2).M();
angle = l2.Angle(l1.Vect());
Float_t e1 = cluster1.E();
Float_t e2 = cluster2.E();
Float_t splitFracCut = 0;
if(nMax < 3) splitFracCut = fSplitEFracMin[nMax-1];
else splitFracCut = fSplitEFracMin[2];
if((e1+e2)/eClus < splitFracCut) return kNeutralUnknown ;
AliDebug(2,"\t pass Split E frac cut");
Float_t minECut = fSubClusterEMin[2];
if (nMax == 2) minECut = fSubClusterEMin[1];
else if(nMax == 1) minECut = fSubClusterEMin[0];
if(e1 < minECut || e2 < minECut)
{
return kNeutralUnknown ;
}
AliDebug(2,"\t pass min sub-cluster E cut");
Float_t asy =-10;
if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
if( !IsInPi0SplitAsymmetryRange(eClus,asy,nMax) ) return kNeutralUnknown ;
AliDebug(2,"\t pass asymmetry cut");
Bool_t pi0OK = kFALSE;
Bool_t etaOK = kFALSE;
Bool_t conOK = kFALSE;
if (IsInPi0M02Range(eClus,m02,nMax)) pi0OK = kTRUE;
else if(IsInEtaM02Range(eClus,m02,nMax)) etaOK = kTRUE;
else if(IsInConM02Range(eClus,m02,nMax)) conOK = kTRUE;
Float_t energy = eClus;
if(nMax > 2) energy = e1+e2;
if ( conOK && mass < fMassPhoMax && mass > fMassPhoMin ) { AliDebug(2,"\t Split Conv"); return kPhoton ; }
else if( etaOK && mass < fMassEtaMax && mass > fMassEtaMin ) { AliDebug(2,"\t Split Eta" ); return kEta ; }
else if( pi0OK && IsInPi0SplitMassRange(energy,mass,nMax) ) { AliDebug(2,"\t Split Pi0" ); return kPi0 ; }
else return kNeutralUnknown ;
}
TString AliCaloPID::GetPIDParametersList()
{
TString parList ;
const Int_t buffersize = 255;
char onePar[buffersize] ;
snprintf(onePar,buffersize,"--- AliCaloPID ---") ;
parList+=onePar ;
if(fUseBayesianWeights)
{
snprintf(onePar,buffersize,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)",fEMCALPhotonWeight) ;
parList+=onePar ;
snprintf(onePar,buffersize,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)",fEMCALPi0Weight) ;
parList+=onePar ;
snprintf(onePar,buffersize,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)",fEMCALElectronWeight) ;
parList+=onePar ;
snprintf(onePar,buffersize,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)",fEMCALChargeWeight) ;
parList+=onePar ;
snprintf(onePar,buffersize,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)",fEMCALNeutralWeight) ;
parList+=onePar ;
snprintf(onePar,buffersize,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)",fPHOSPhotonWeight) ;
parList+=onePar ;
snprintf(onePar,buffersize,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)",fPHOSPi0Weight) ;
parList+=onePar ;
snprintf(onePar,buffersize,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)",fPHOSElectronWeight) ;
parList+=onePar ;
snprintf(onePar,buffersize,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)",fPHOSChargeWeight) ;
parList+=onePar ;
snprintf(onePar,buffersize,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)",fPHOSNeutralWeight) ;
parList+=onePar ;
if(fPHOSWeightFormula)
{
snprintf(onePar,buffersize,"PHOS Photon Weight Formula: %s",fPHOSPhotonWeightFormulaExpression.Data() ) ;
parList+=onePar;
snprintf(onePar,buffersize,"PHOS Pi0 Weight Formula: %s",fPHOSPi0WeightFormulaExpression.Data() ) ;
parList+=onePar;
}
}
else
{
snprintf(onePar,buffersize,"EMCAL: fEMCALL0CutMin =%2.2f, fEMCALL0CutMax =%2.2f (Cut on Shower Shape)",fEMCALL0CutMin, fEMCALL0CutMax) ;
parList+=onePar ;
snprintf(onePar,buffersize,"EMCAL: fEMCALDEtaCut =%2.2f, fEMCALDPhiCut =%2.2f (Cut on track matching)",fEMCALDEtaCut, fEMCALDPhiCut) ;
parList+=onePar ;
snprintf(onePar,buffersize,"fTOFCut =%e (Cut on TOF, used in PID evaluation)",fTOFCut) ;
parList+=onePar ;
snprintf(onePar,buffersize,"fPHOSRCut =%2.2f, fPHOSDispersionCut =%2.2f (Cut on Shower Shape and CPV)",fPHOSRCut,fPHOSDispersionCut) ;
parList+=onePar ;
}
if(fUseSimpleM02Cut)
{
snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f", fSplitM02MinCut, fSplitM02MaxCut) ;
parList+=onePar ;
}
snprintf(onePar,buffersize,"fMinNCells =%d", fSplitMinNCells) ;
parList+=onePar ;
if(fUseSimpleMassCut)
{
snprintf(onePar,buffersize,"pi0 : %2.1f < m <%2.1f", fMassPi0Min,fMassPi0Max);
parList+=onePar ;
}
snprintf(onePar,buffersize,"eta : %2.1f < m <%2.1f", fMassEtaMin,fMassEtaMax);
parList+=onePar ;
snprintf(onePar,buffersize,"conv: %2.1f < m <%2.1f", fMassPhoMin,fMassPhoMax);
parList+=onePar ;
return parList;
}
void AliCaloPID::Print(const Option_t * opt) const
{
if(! opt)
return;
printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
if(fUseBayesianWeights)
{
printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",
fPHOSPhotonWeight, fPHOSPi0Weight,
fPHOSElectronWeight, fPHOSChargeWeight, fPHOSNeutralWeight) ;
printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",
fEMCALPhotonWeight, fEMCALPi0Weight,
fEMCALElectronWeight, fEMCALChargeWeight, fEMCALNeutralWeight) ;
printf("PHOS Parametrized weight on? = %d\n", fPHOSWeightFormula) ;
if(fPHOSWeightFormula)
{
printf("Photon weight formula = %s\n", fPHOSPhotonWeightFormulaExpression.Data());
printf("Pi0 weight formula = %s\n", fPHOSPi0WeightFormulaExpression .Data());
}
if(fRecalculateBayesian) printf(" Recalculate bayesian with Particle Flux? = %d\n",fParticleFlux);
}
else
{
printf("TOF cut = %e\n", fTOFCut);
printf("EMCAL Lambda0 cut min = %2.2f; max = %2.2f\n", fEMCALL0CutMin,fEMCALL0CutMax);
printf("EMCAL cluster-track dEta < %2.3f; dPhi < %2.3f\n", fEMCALDEtaCut, fEMCALDPhiCut);
printf("PHOS Treac matching cut =%2.2f, Dispersion Cut =%2.2f \n",fPHOSRCut, fPHOSDispersionCut) ;
}
printf("Min. N Cells =%d \n", fSplitMinNCells) ;
if(fUseSimpleM02Cut) printf("%2.2f < lambda_0^2 <%2.2f \n",fSplitM02MinCut,fSplitM02MaxCut);
if(fUseSimpleMassCut)printf("pi0 : %2.2f<m<%2.2f \n", fMassPi0Min,fMassPi0Max);
printf("eta : %2.2f<m<%2.2f \n", fMassEtaMin,fMassEtaMax);
printf("phot: %2.2f<m<%2.2f \n", fMassPhoMin,fMassPhoMax);
printf(" \n");
}
void AliCaloPID::PrintClusterPIDWeights(const Double_t * pid) const
{
printf("AliCaloPID::PrintClusterPIDWeights() \n \t ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, \n \t \
pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
pid[AliVCluster::kPhoton], pid[AliVCluster::kPi0],
pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon],
pid[AliVCluster::kPion], pid[AliVCluster::kKaon],
pid[AliVCluster::kProton],
pid[AliVCluster::kNeutron], pid[AliVCluster::kKaon0]);
}
void AliCaloPID::SetPIDBits(AliVCluster * cluster,
AliAODPWG4Particle * ph, AliCalorimeterUtils* cu,
AliVEvent* event)
{
Double_t l1 = cluster->GetM20() ;
Double_t l0 = cluster->GetM02() ;
Bool_t isDispOK = kTRUE ;
if(cluster->IsPHOS()){
if(TestPHOSDispersion(ph->Pt(),l0,l1) < fPHOSDispersionCut) isDispOK = kTRUE;
else isDispOK = kFALSE;
}
else{
if(l0 > fEMCALL0CutMin && l0 < fEMCALL0CutMax) isDispOK = kTRUE;
}
ph->SetDispBit(isDispOK) ;
Double_t tof=cluster->GetTOF() ;
ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
Bool_t isNeutral = IsTrackMatched(cluster,cu,event);
ph->SetChargedBit(isNeutral);
ph->SetIdentifiedParticleType(GetIdentifiedParticleType(cluster));
AliDebug(1,Form("TOF %e, Lambda0 %2.2f, Lambda1 %2.2f",tof , l0, l1));
AliDebug(1,Form("pdg %d, bits: TOF %d, Dispersion %d, Charge %d",
ph->GetIdentifiedParticleType(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit()));
}
Bool_t AliCaloPID::IsTrackMatched(AliVCluster* cluster,
AliCalorimeterUtils * cu,
AliVEvent* event) const
{
Int_t nMatches = cluster->GetNTracksMatched();
AliVTrack * track = 0;
if(nMatches > 0)
{
if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
{
Int_t iESDtrack = cluster->GetTrackMatchedIndex();
if(iESDtrack >= 0) track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
else return kFALSE;
if (!track)
{
AliWarning("Null matched track in ESD when index is OK!");
return kFALSE;
}
}
else
{
track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0));
if (!track)
{
AliWarning("Null matched track in AOD!");
return kFALSE;
}
}
Float_t dZ = cluster->GetTrackDz();
Float_t dR = cluster->GetTrackDx();
if(cluster->IsEMCAL() && cu && cu->IsRecalculationOfClusterTrackMatchingOn())
{
dR = 2000., dZ = 2000.;
cu->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
}
if(cluster->IsPHOS())
{
Int_t charge = track->Charge();
Double_t mf = event->GetMagneticField();
if(TestPHOSChargedVeto(dR, dZ, track->Pt(), charge, mf ) < fPHOSRCut) return kTRUE;
else return kFALSE;
}
else
{
AliDebug(1,Form("EMCAL dR %f < %f, dZ %f < %f ",dR, fEMCALDPhiCut, dZ, fEMCALDEtaCut));
if(TMath::Abs(dR) < fEMCALDPhiCut &&
TMath::Abs(dZ) < fEMCALDEtaCut) return kTRUE;
else return kFALSE;
}
}
else return kFALSE;
}
Float_t AliCaloPID::TestPHOSDispersion(const Double_t pt, const Double_t l1, const Double_t l2) const
{
Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
Double_t c =-0.35-0.550*TMath::Exp(-0.390730*pt) ;
Double_t r2 = 0.5* (l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
0.5* (l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
AliDebug(1,Form("PHOS SS R %f < %f?", TMath::Sqrt(r2), fPHOSDispersionCut));
return TMath::Sqrt(r2) ;
}
Float_t AliCaloPID::TestPHOSChargedVeto(const Double_t dx, const Double_t dz, const Double_t pt,
const Int_t charge, const Double_t mf) const
{
Double_t meanX = 0.;
Double_t meanZ = 0.;
Double_t sx = TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+
1.59219);
Double_t sz = TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+
1.60) ;
if(mf<0.){
meanZ = -0.468318 ;
if(charge>0)
meanX = TMath::Min(7.3, 3.89994*1.20679 *1.20679 /(pt*pt+1.20679*1.20679)+
0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
else
meanX =-TMath::Min(7.7, 3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+
1.23114 +4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
}
else{
meanZ = -0.468318;
if(charge>0)
meanX =-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+
0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
else
meanX = TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-
0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
}
Double_t rz = (dz-meanZ)/sz ;
Double_t rx = (dx-meanX)/sx ;
AliDebug(1,Form("PHOS Matching R %f < %f",TMath::Sqrt(rx*rx+rz*rz), fPHOSRCut));
return TMath::Sqrt(rx*rx+rz*rz) ;
}