]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/PartCorrBase/AliCaloPID.cxx
bugfixes for CALO/PHOS in the HLT build system
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliCaloPID.cxx
index 0a5eb46f5f177ba466dcfcb44b30effcdb3de2bd..46915a006f0905874305210053765736d925dd0b 100755 (executable)
 /* $Id: AliCaloPID.cxx 21839 2007-10-29 13:49:42Z gustavo $ */
 
 //_________________________________________________________________________
-// Class for track/cluster acceptance selection
-// Selection in Central barrel, EMCAL and PHOS
-//                
+// Class for PID selection with calorimeters
+// The Output of the 2 main methods GetPdg is a PDG number identifying the cluster, 
+// being kPhoton, kElectron, kPi0 ... as defined in the header file
+//   - GetPdg(const TString calo, const Double_t * pid, const Float_t energy)
+//      Reads the PID weights array of the ESDs and depending on its magnitude identifies the particle
+//   - GetPdg(const TString calo,const TLorentzVector mom, const AliAODCaloCluster * cluster)
+//      Recalcultes PID, the bayesian or any new one to be implemented in the future
+//      Right now only the possibility to recalculate EMCAL with bayesian and simple PID.
+//      In order to recalculate Bayesian, it is necessary to load the EMCALUtils library
+//      and do SwitchOnBayesianRecalculation().
+//      To change the PID parameters from Low to High like the ones by default, use the constructor 
+//      AliCaloPID(flux)
+//      where flux is AliCaloPID::kLow or AliCaloPID::kHigh
+//      If it is necessary to change the parameters use the constructor 
+//      AliCaloPID(AliEMCALPIDUtils *utils) and set the parameters before.
+//   - SetPIDBits: Simple PID, depending on the thresholds fDispCut fTOFCut and even the
+//     result of the PID bayesian a different PID bit is set. 
+//
+//  All these methods can be called in the analysis you are interested.
+//
 //*-- Author: Gustavo Conesa (LNF-INFN) 
 //////////////////////////////////////////////////////////////////////////////
 
@@ -31,7 +48,7 @@
 #include "AliCaloPID.h"
 #include "AliAODCaloCluster.h"
 #include "AliAODPWG4Particle.h"
-
+#include "AliEMCALPIDUtils.h"
 ClassImp(AliCaloPID)
 
 
@@ -44,7 +61,8 @@ fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
 fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) , 
 fPHOSNeutralWeight(0.), fPHOSWeightFormula(0), 
 fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0),
-fDispCut(0.),fTOFCut(0.), fDebug(-1)
+fDispCut(0.),fTOFCut(0.), fDebug(-1), 
+fRecalculateBayesian(kFALSE), fParticleFlux(kLow), fEMCALPIDUtils(new AliEMCALPIDUtils)
 {
   //Ctor
   
@@ -52,6 +70,42 @@ fDispCut(0.),fTOFCut(0.), fDebug(-1)
   InitParameters();
 }
 
+//________________________________________________
+AliCaloPID::AliCaloPID(const Int_t flux) : 
+TObject(), fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),  
+fEMCALElectronWeight(0.),  fEMCALChargeWeight(0.),
+fEMCALNeutralWeight(0.),
+fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),  
+fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) , 
+fPHOSNeutralWeight(0.), fPHOSWeightFormula(0), 
+fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0),
+fDispCut(0.),fTOFCut(0.), fDebug(-1), 
+fRecalculateBayesian(kTRUE), fParticleFlux(flux), fEMCALPIDUtils(new AliEMCALPIDUtils)
+{
+       //Ctor
+       
+       //Initialize parameters
+       InitParameters();
+}
+
+//________________________________________________
+AliCaloPID::AliCaloPID(const TTask * emcalpid) : 
+TObject(), fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),  
+fEMCALElectronWeight(0.),  fEMCALChargeWeight(0.),
+fEMCALNeutralWeight(0.),
+fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),  
+fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) , 
+fPHOSNeutralWeight(0.), fPHOSWeightFormula(0), 
+fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0),
+fDispCut(0.),fTOFCut(0.), fDebug(-1), 
+fRecalculateBayesian(kTRUE), fParticleFlux(-1), fEMCALPIDUtils( (AliEMCALPIDUtils*) emcalpid)
+{
+       //Ctor
+       
+       //Initialize parameters
+       InitParameters();
+}
+
 //____________________________________________________________________________
 AliCaloPID::AliCaloPID(const AliCaloPID & pid) :   
 TObject(pid), fEMCALPhotonWeight(pid.fEMCALPhotonWeight), 
@@ -68,7 +122,8 @@ fPHOSWeightFormula(pid.fPHOSWeightFormula),
 fPHOSPhotonWeightFormula(pid.fPHOSPhotonWeightFormula), 
 fPHOSPi0WeightFormula(pid.fPHOSPi0WeightFormula), 
 fDispCut(pid.fDispCut),fTOFCut(pid.fTOFCut),
-fDebug(pid.fDebug)
+fDebug(pid.fDebug), fRecalculateBayesian(pid.fRecalculateBayesian),
+fParticleFlux(pid.fParticleFlux), fEMCALPIDUtils(pid.fEMCALPIDUtils)
 {
   // cpy ctor
        
@@ -101,6 +156,10 @@ AliCaloPID & AliCaloPID::operator = (const AliCaloPID & pid)
   fTOFCut   = pid.fTOFCut;
   fDebug    = pid.fDebug;
   
+  fRecalculateBayesian = pid.fRecalculateBayesian;
+  fParticleFlux        = pid.fParticleFlux;
+  fEMCALPIDUtils       = pid.fEMCALPIDUtils;
+               
   return *this;
   
 }
@@ -111,7 +170,7 @@ AliCaloPID::~AliCaloPID() {
   
   if(fPHOSPhotonWeightFormula) delete  fPHOSPhotonWeightFormula ;
   if(fPHOSPi0WeightFormula) delete  fPHOSPi0WeightFormula ;
-  
+  if(fEMCALPIDUtils) delete fEMCALPIDUtils ;
 }
 
 
@@ -120,14 +179,14 @@ void AliCaloPID::InitParameters()
 {
   //Initialize the parameters of the PID.
   
-  fEMCALPhotonWeight   = 0.8 ;
+  fEMCALPhotonWeight   = 0.5 ;
   fEMCALPi0Weight      = 0.5 ;
-  fEMCALElectronWeight = 0.8 ;
+  fEMCALElectronWeight = 0.5 ;
   fEMCALChargeWeight   = 0.5 ;
   fEMCALNeutralWeight  = 0.5 ;
   
-  fPHOSPhotonWeight    = 0.75 ;
-  fPHOSPi0Weight       = 0.8 ;
+  fPHOSPhotonWeight    = 0.5 ;
+  fPHOSPi0Weight       = 0.5 ;
   fPHOSElectronWeight  = 0.5 ;
   fPHOSChargeWeight    = 0.5 ;
   fPHOSNeutralWeight   = 0.5 ;
@@ -142,6 +201,17 @@ void AliCaloPID::InitParameters()
   fDispCut  = 1.5;
   fTOFCut   = 5.e-9;
   fDebug = -1;
+       
+  if(fRecalculateBayesian){
+       if(fParticleFlux == kLow){
+               printf("AliCaloPID::Init() - SetLOWFluxParam\n");
+               fEMCALPIDUtils->SetLowFluxParam() ;
+       }
+       else if (fParticleFlux == kHigh){
+               printf("AliCaloPID::Init() - SetHIGHFluxParam\n");
+               fEMCALPIDUtils->SetHighFluxParam() ;
+       }
+  }
 }
 
 //_______________________________________________________________
@@ -153,11 +223,11 @@ Int_t AliCaloPID::GetPdg(const TString calo, const Double_t * pid, const Float_t
     abort();
   }
   
-  Float_t wPh =  fPHOSPhotonWeight ;
+  Float_t wPh  =  fPHOSPhotonWeight ;
   Float_t wPi0 =  fPHOSPi0Weight ;
-  Float_t wE =  fPHOSElectronWeight ;
-  Float_t wCh =  fPHOSChargeWeight ;
-  Float_t wNe =  fPHOSNeutralWeight ;
+  Float_t wE   =  fPHOSElectronWeight ;
+  Float_t wCh  =  fPHOSChargeWeight ;
+  Float_t wNe  =  fPHOSNeutralWeight ;
   
   
   if(calo == "PHOS" && fPHOSWeightFormula){
@@ -202,9 +272,11 @@ Int_t AliCaloPID::GetPdg(const TString calo, const Double_t * pid, const Float_t
       pdg = kNeutralUnknown ;
   }
   else{//EMCAL
-    if(pid[AliAODCluster::kPhoton]  > wPh) pdg = kPhoton ;
+
+    if(pid[AliAODCluster::kPhoton]+pid[AliAODCluster::kElectron]  > wPh) pdg = kPhoton ; //temporal sollution until track matching for electrons is considered
+    //if(pid[AliAODCluster::kPhoton]  > wPh) pdg = kPhoton ;
     else if(pid[AliAODCluster::kPi0] > wPi0) pdg = kPi0 ; 
-    else if(pid[AliAODCluster::kElectron]  > wE) pdg = kElectron ;
+    //else if(pid[AliAODCluster::kElectron]  > wE) pdg = kElectron ;
     else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;  
     else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ; 
     else pdg =  kNeutralUnknown ;
@@ -220,13 +292,37 @@ Int_t AliCaloPID::GetPdg(const TString calo, const Double_t * pid, const Float_t
 //_______________________________________________________________
 Int_t AliCaloPID::GetPdg(const TString calo,const TLorentzVector mom, const AliAODCaloCluster * cluster) const {
   //Recalculated PID with all parameters
-  if(fDebug > 0)printf("AliCaloPID::GetPdg: 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(),mom.E(),cluster->GetM02(),cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(), 
-                      cluster->GetEmcCpvDistance(), cluster->GetDistToBadChannel(),cluster->GetNExMax());
-  
+  Float_t lambda0 = cluster->GetM02();
+  Float_t energy = mom.E();    
+
+  if(fDebug > 0) printf("AliCaloPID::GetPdg: 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(), 
+                                               cluster->GetEmcCpvDistance(), cluster->GetDistToBadChannel(),cluster->GetNExMax());
+
   if(calo == "EMCAL") {
-    if(cluster->GetM02()< 0.25) return kPhoton ;
-    else return  kNeutralHadron ; 
+         //Recalculate Bayesian
+         if(fRecalculateBayesian){       
+                 if(fDebug > 0)  {
+                         const Double_t  *pid0 = cluster->PID();
+                         printf("AliCaloPID::GetPdg: 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(),pid0[AliAODCluster::kPhoton], pid0[AliAODCluster::kPi0],
+                                                                pid0[AliAODCluster::kElectron], pid0[AliAODCluster::kEleCon],
+                                                                pid0[AliAODCluster::kPion], pid0[AliAODCluster::kKaon], pid0[AliAODCluster::kProton],
+                                                                pid0[AliAODCluster::kNeutron], pid0[AliAODCluster::kKaon0]);
+                 }
+                 
+                fEMCALPIDUtils->ComputePID(energy, lambda0);
+                Double_t pid[AliPID::kSPECIESN];
+                for(Int_t i = 0; i < AliPID::kSPECIESN; i++) pid[i] = fEMCALPIDUtils->GetPIDFinal(i);
+                return GetPdg(calo, pid, energy);
+                 
+         
+       }
+         
+       // If no use of bayesian, simple PID  
+    if(lambda0 < 0.25) return kPhoton ;
+    //else return  kNeutralHadron ; 
+       else return  kPi0 ;
   }
   
   //   if(calo == "PHOS") {
@@ -272,8 +368,10 @@ TString  AliCaloPID::GetPIDParametersList()  {
   parList+=onePar ;
   
   if(fPHOSWeightFormula){
-    parList+="PHOS Photon Weight Formula: "+(fPHOSPhotonWeightFormula->GetExpFormula("p"));
-    parList+="PHOS Pi0    Weight Formula: "+(fPHOSPi0WeightFormula->GetExpFormula("p"));
+       sprintf(onePar,"PHOS Photon Weight Formula: %s\n",(fPHOSPhotonWeightFormula->GetExpFormula("p")).Data()) ;
+    parList+=onePar;
+       sprintf(onePar,"PHOS Pi0    Weight Formula: %s\n",(fPHOSPi0WeightFormula->GetExpFormula("p")).Data()) ;
+       parList+=onePar;          
   }
   
   return parList; 
@@ -306,7 +404,8 @@ void AliCaloPID::Print(const Option_t * opt) const
   printf("TOF cut        = %e\n",fTOFCut);
   printf("Dispersion cut = %2.2f\n",fDispCut);
   printf("Debug level    = %d\n",fDebug);
-  
+  printf("Recalculate Bayesian?    = %d\n",fRecalculateBayesian);
+ if(fRecalculateBayesian) printf("Particle Flux?    = %d\n",fParticleFlux);
   printf(" \n");
   
 }