// Class for PID selection with calorimeters
// The Output of the main method GetIdentifiedParticleType is a PDG number identifying the cluster,
// being kPhoton, kElectron, kPi0 ... as defined in the header file
-// - GetIdentifiedParticleType(const TString calo, const TLorentzVector mom, const AliVCluster * cluster)
+// - GetIdentifiedParticleType(const AliVCluster * cluster)
// Assignes a PID tag to the cluster, right now there is the possibility to : use bayesian weights from reco,
// recalculate them (EMCAL) or use other procedures not used in reco.
// In order to recalculate Bayesian, it is necessary to load the EMCALUtils library
// If it is necessary to change the parameters use the constructor
// AliCaloPID(AliEMCALPIDUtils *utils) and set the parameters before.
-// - GetGetIdentifiedParticleTypeFromBayesian(const TString calo, const Double_t * pid, const Float_t energy)
+// - GetGetIdentifiedParticleTypeFromBayesian(const Double_t * pid, const Float_t energy)
// Reads the PID weights array of the ESDs and depending on its magnitude identifies the particle,
-// executed when bayesian is ON by GetIdentifiedParticleType(const TString calo, const TLorentzVector mom, const AliVCluster * cluster)
+// executed when bayesian is ON by GetIdentifiedParticleType(const AliVCluster * cluster)
// - SetPIDBits: Simple PID, depending on the thresholds fLOCut fTOFCut and even the
// result of the PID bayesian a different PID bit is set.
//
// ---- ANALYSIS system ----
#include "AliCaloPID.h"
-#include "AliVCluster.h"
+#include "AliAODCaloCluster.h"
+#include "AliVCaloCells.h"
#include "AliVTrack.h"
#include "AliAODPWG4Particle.h"
#include "AliCalorimeterUtils.h"
fEMCALL0CutMax(100.), fEMCALL0CutMin(0),
fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.),
fTOFCut(0.),
-fPHOSDispersionCut(1000), fPHOSRCut(1000)
+fPHOSDispersionCut(1000), fPHOSRCut(1000),
+fDoClusterSplitting(kFALSE),
+fSplitM02MaxCut(0), fSplitM02MinCut(0), fSplitMinNCells(0),
+fMassEtaMin(0), fMassEtaMax(0),
+fMassPi0Min(0), fMassPi0Max(0),
+fMassPhoMin(0), fMassPhoMax(0)
{
//Ctor
fEMCALL0CutMax(100.), fEMCALL0CutMin(0),
fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.),
fTOFCut(0.),
-fPHOSDispersionCut(1000), fPHOSRCut(1000)
+fPHOSDispersionCut(1000), fPHOSRCut(1000),
+fDoClusterSplitting(kFALSE),
+fSplitM02MaxCut(0), fSplitM02MinCut(0), fSplitMinNCells(0),
+fMassEtaMin(0), fMassEtaMax(0),
+fMassPi0Min(0), fMassPi0Max(0),
+fMassPhoMin(0), fMassPhoMax(0)
{
//Ctor
fEMCALL0CutMax(100.), fEMCALL0CutMin(0),
fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.),
fTOFCut(0.),
-fPHOSDispersionCut(1000), fPHOSRCut(1000)
+fPHOSDispersionCut(1000), fPHOSRCut(1000),
+fDoClusterSplitting(kFALSE),
+fSplitM02MaxCut(0), fSplitM02MinCut(0), fSplitMinNCells(0),
+fMassEtaMin(0), fMassEtaMax(0),
+fMassPi0Min(0), fMassPi0Max(0),
+fMassPhoMin(0), fMassPhoMax(0)
{
//Ctor
fPHOSRCut = 2. ;
fPHOSDispersionCut = 2.5;
+ // Cluster splitting
+
+ fSplitM02MinCut = 0.26 ;
+ fSplitM02MaxCut = 100 ;
+ fSplitMinNCells = 4 ;
+
+ fMassEtaMin = 0.4;
+ fMassEtaMax = 0.6;
+
+ fMassPi0Min = 0.08;
+ fMassPi0Max = 0.20;
+
+ fMassPhoMin = 0.0;
+ fMassPhoMax = 0.05;
+
}
//______________________________________________
//______________________________________________________________________
-Int_t AliCaloPID::GetIdentifiedParticleType(const TString calo,
- const TLorentzVector mom,
- const AliVCluster * cluster)
+Int_t AliCaloPID::GetIdentifiedParticleType(const AliVCluster * cluster)
{
// Returns a PDG number corresponding to the likely ID of the cluster
- Float_t energy = mom.E();
+ Float_t energy = cluster->E();
Float_t lambda0 = cluster->GetM02();
Float_t lambda1 = cluster->GetM20();
// Use bayesian approach
// ---------------------
- if(fUseBayesianWeights){
-
+ if(fUseBayesianWeights)
+ {
Double_t weights[AliPID::kSPECIESN];
- if(calo == "EMCAL"&& fRecalculateBayesian){
+ if(cluster->IsEMCAL() && fRecalculateBayesian)
+ {
fEMCALPIDUtils->ComputePID(energy, lambda0);
for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = fEMCALPIDUtils->GetPIDFinal(i);
}
- else {
+ else
+ {
for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = cluster->GetPID()[i];
}
- if(fDebug > 0) {
- printf("AliCaloPID::GetIdentifiedParticleType: BEFORE calo %s, ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, hadrons: pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
- calo.Data(),
- weights[AliVCluster::kPhoton], weights[AliVCluster::kPi0],
- weights[AliVCluster::kElectron], weights[AliVCluster::kEleCon],
- weights[AliVCluster::kPion], weights[AliVCluster::kKaon],
- weights[AliVCluster::kProton],
- weights[AliVCluster::kNeutron], weights[AliVCluster::kKaon0]);
- }
+ if(fDebug > 0) PrintClusterPIDWeights(weights);
- return GetIdentifiedParticleTypeFromBayesWeights(calo, weights, energy);
+ return GetIdentifiedParticleTypeFromBayesWeights(cluster->IsEMCAL(), weights, energy);
}
// -------------------------------------------------------
// Calculate PID SS from data, do not use bayesian weights
// -------------------------------------------------------
- if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType: Calorimeter %s, E %3.2f, l0 %3.2f, l1 %3.2f, disp %3.2f, tof %1.11f, distCPV %3.2f, distToBC %1.1f, NMax %d\n",
- calo.Data(),energy,lambda0,cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(),
+ if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType: 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\n",
+ cluster->IsEMCAL(),energy,lambda0,cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(),
cluster->GetEmcCpvDistance(), cluster->GetDistanceToBadChannel(),cluster->GetNExMax());
- if(cluster->IsEMCAL()){
-
+ if(cluster->IsEMCAL())
+ {
if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType() - EMCAL SS %f <%f < %f?\n",fEMCALL0CutMin, lambda0, fEMCALL0CutMax);
if(lambda0 < fEMCALL0CutMax && lambda0 > fEMCALL0CutMin) return kPhoton ;
else return kNeutralUnknown ;
- }//EMCAL
- else {//PHOS
- if(TestPHOSDispersion(mom.Pt(),lambda0,lambda1) < fPHOSDispersionCut) return kPhoton;
- else return kNeutralUnknown;
+ } // EMCAL
+ else // PHOS
+ {
+ if(TestPHOSDispersion(energy,lambda0,lambda1) < fPHOSDispersionCut) return kPhoton;
+ else return kNeutralUnknown;
}
}
//_______________________________________________________________________________
-Int_t AliCaloPID::GetIdentifiedParticleTypeFromBayesWeights(const TString calo,
+Int_t AliCaloPID::GetIdentifiedParticleTypeFromBayesWeights(const Bool_t isEMCAL,
const Double_t * pid,
const Float_t energy)
{
//Return most probable identity of the particle after bayesian weights calculated in reconstruction
- if(!pid){
+ if(!pid)
+ {
printf("AliCaloPID::GetIdentifiedParticleType() - pid pointer not initialized!!!\n");
abort();
}
Float_t wCh = fPHOSChargeWeight ;
Float_t wNe = fPHOSNeutralWeight ;
- if(calo == "PHOS" && fPHOSWeightFormula){
+ if(!isEMCAL && fPHOSWeightFormula){
wPh = GetPHOSPhotonWeightFormula()->Eval(energy) ;
wPi0 = GetPHOSPi0WeightFormula() ->Eval(energy);
}
-
- if(calo == "EMCAL"){
-
+ else
+ {
wPh = fEMCALPhotonWeight ;
wPi0 = fEMCALPi0Weight ;
wE = fEMCALElectronWeight ;
wCh = fEMCALChargeWeight ;
wNe = fEMCALNeutralWeight ;
-
}
- if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType: calo %s, ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, hadrons: pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
- calo.Data(),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]);
-
+ 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 allNeutralWeight = pid[AliVCluster::kPhoton]+pid[AliVCluster::kPi0]+ neutralHadronWeight;
//Select most probable ID
- if(calo=="PHOS"){
+ if(!isEMCAL) // PHOS
+ {
if(pid[AliVCluster::kPhoton] > wPh) pdg = kPhoton ;
else if(pid[AliVCluster::kPi0] > wPi0) pdg = kPi0 ;
else if(pid[AliVCluster::kElectron] > wE) pdg = kElectron ;
else
pdg = kNeutralUnknown ;
}
- else{//EMCAL
-
+ else //EMCAL
+ {
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 ; //temporal sollution until track matching for electrons is considered
}
+//____________________________________________________________________________________________________
+Int_t AliCaloPID::GetIdentifiedParticleTypeFromClusterSplitting(AliVCluster* cluster,
+ AliVCaloCells* cells,
+ AliCalorimeterUtils * caloutils,
+ Double_t vertex[3],
+ Int_t & nMax,
+ Double_t & mass, Double_t & angle)
+{
+ // Split the cluster in 2, do invariant mass, get the mass and decide
+ // if this is a photon, pi0, eta, ...
+
+ Int_t absId1 = -1; Int_t absId2 = -1;
+ Int_t nc = cluster->GetNCells();
+ Int_t *absIdList = new Int_t [nc];
+ Float_t *maxEList = new Float_t[nc];
+ Float_t l0 = cluster->GetM02();
+
+ nMax = -1 ;
+ mass = -1.;
+ angle = -1.;
+
+ //If too small or big E or low number of cells, or close to a bad channel skip it
+ if( l0 < fSplitM02MinCut || l0 > fSplitM02MaxCut || nc < fSplitMinNCells) return kNeutralUnknown ;
+
+ // Get Number of local maxima
+ nMax = caloutils->GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ;
+
+ //---------------------------------------------------------------------
+ // Get the 2 max indeces and do inv mass
+ //---------------------------------------------------------------------
+
+ if ( nMax == 2 )
+ {
+ absId1 = absIdList[0];
+ absId2 = absIdList[1];
+ }
+ else if( nMax == 1 )
+ {
+
+ absId1 = absIdList[0];
+
+ //Find second highest energy cell
+
+ TString calorimeter = "EMCAL";
+ if(cluster->IsPHOS()) calorimeter = "PHOS";
+ 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 ;
+ }
+ }// cell loop
+ }// 1 maxima
+ else
+ { // n max > 2
+ // loop on maxima, find 2 highest
+
+ // First max
+ 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];
+ }
+ }// first maxima loop
+
+ // Second max
+ 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];
+ }
+ }// second maxima loop
+
+ } // n local maxima > 2
+
+ //---------------------------------------------------------------------
+ // Split the cluster energy in 2, around the highest 2 local maxima
+ //---------------------------------------------------------------------
+
+ AliAODCaloCluster *cluster1 = new AliAODCaloCluster(0, 0,NULL,0.,NULL,NULL,1,0);
+ AliAODCaloCluster *cluster2 = new AliAODCaloCluster(1, 0,NULL,0.,NULL,NULL,1,0);
+
+ caloutils->SplitEnergy(absId1,absId2,cluster, cells, cluster1, cluster2,nMax); /*absIdList, maxEList,*/
+
+ TLorentzVector cellMom1;
+ TLorentzVector cellMom2;
+
+ cluster1->GetMomentum(cellMom1,vertex);
+ cluster2->GetMomentum(cellMom2,vertex);
+
+ mass = (cellMom1+cellMom2).M();
+ angle = cellMom2.Angle(cellMom1.Vect());
+
+ if (mass < fMassPhoMax && mass > fMassPhoMin) return kPhoton;
+ else if(mass < fMassPi0Max && mass > fMassPi0Min) return kPi0;
+ else if(mass < fMassEtaMax && mass > fMassEtaMin) return kEta;
+ else return kNeutralUnknown;
+
+ delete cluster1;
+ delete cluster2;
+ delete [] absIdList ;
+ delete [] maxEList ;
+
+}
+
//_________________________________________
TString AliCaloPID::GetPIDParametersList()
{
}
+ if(fDoClusterSplitting)
+ {
+ snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n", fSplitM02MinCut, fSplitM02MaxCut) ;
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"fMinNCells =%d \n", fSplitMinNCells) ;
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"pi0 : %2.1f < m <%2.1f\n", fMassPi0Min,fMassPi0Max);
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"eta : %2.1f < m <%2.1f\n", fMassEtaMin,fMassEtaMax);
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"conv: %2.1f < m <%2.1f\n", fMassPhoMin,fMassPhoMax);
+ parList+=onePar ;
+ }
+
return parList;
}
printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
- if(fUseBayesianWeights){
+ 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) ;
+ 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) ;
+ fEMCALPhotonWeight, fEMCALPi0Weight,
+ fEMCALElectronWeight, fEMCALChargeWeight, fEMCALNeutralWeight) ;
printf("PHOS Parametrized weight on? = %d\n", fPHOSWeightFormula) ;
- if(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) ;
+ 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) ;
}
+ if(fDoClusterSplitting)
+ {
+ printf("Min. N Cells =%d \n", fSplitMinNCells) ;
+ printf("%2.2f < lambda_0^2 <%2.2f \n",fSplitM02MinCut,fSplitM02MaxCut);
+ 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
+{
+ // print PID of cluster, (AliVCluster*)cluster->GetPID()
+
+ 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(const TString calo, AliVCluster * cluster,
+void AliCaloPID::SetPIDBits(AliVCluster * cluster,
AliAODPWG4Particle * ph, AliCalorimeterUtils* cu,
AliVEvent* event)
{
ph->SetChargedBit(isNeutral);
//Set PID pdg
- ph->SetIdentifiedParticleType(GetIdentifiedParticleType(calo,*ph->Momentum(),cluster));
+ ph->SetIdentifiedParticleType(GetIdentifiedParticleType(cluster));
if(fDebug > 0){
printf("AliCaloPID::SetPIDBits: TOF %e, Lambda0 %2.2f, Lambda1 %2.2f\n",tof , l0, l1);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff