[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliCaloPID.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
/* $Id: AliCaloPID.cxx 21839 2007-10-29 13:49:42Z gustavo $ */

//_________________________________________________________________________
// 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) 
//////////////////////////////////////////////////////////////////////////////


// --- ROOT system ---
#include <TMath.h>
#include <TString.h>
#include <TFormula.h>

//---- ANALYSIS system ----
#include "AliCaloPID.h"
#include "AliAODCaloCluster.h"
#include "AliAODPWG4Particle.h"
#include "AliEMCALPIDUtils.h"
ClassImp(AliCaloPID)


//________________________________________________
AliCaloPID::AliCaloPID() : 
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(kFALSE), fParticleFlux(kLow), fEMCALPIDUtils(new AliEMCALPIDUtils)
{
  //Ctor
  
  //Initialize parameters
  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), 
fEMCALPi0Weight(pid.fEMCALPi0Weight), 
fEMCALElectronWeight(pid.fEMCALElectronWeight), 
fEMCALChargeWeight(pid.fEMCALChargeWeight), 
fEMCALNeutralWeight(pid.fEMCALNeutralWeight), 
fPHOSPhotonWeight(pid.fPHOSPhotonWeight),
fPHOSPi0Weight(pid.fPHOSPi0Weight),
fPHOSElectronWeight(pid.fPHOSElectronWeight), 
fPHOSChargeWeight(pid.fPHOSChargeWeight),
fPHOSNeutralWeight(pid.fPHOSNeutralWeight),
fPHOSWeightFormula(pid.fPHOSWeightFormula), 
fPHOSPhotonWeightFormula(pid.fPHOSPhotonWeightFormula), 
fPHOSPi0WeightFormula(pid.fPHOSPi0WeightFormula), 
fDispCut(pid.fDispCut),fTOFCut(pid.fTOFCut),
fDebug(pid.fDebug), fRecalculateBayesian(pid.fRecalculateBayesian),
fParticleFlux(pid.fParticleFlux), fEMCALPIDUtils(pid.fEMCALPIDUtils)
{
  // cpy ctor
	
}

//_________________________________________________________________________
AliCaloPID & AliCaloPID::operator = (const AliCaloPID & pid)
{
  // assignment operator
  
  if(&pid == this) return *this;
  
  fEMCALPhotonWeight = pid. fEMCALPhotonWeight ;
  fEMCALPi0Weight = pid.fEMCALPi0Weight ;
  fEMCALElectronWeight = pid.fEMCALElectronWeight; 
  fEMCALChargeWeight = pid.fEMCALChargeWeight;
  fEMCALNeutralWeight = pid.fEMCALNeutralWeight;
  
  fPHOSPhotonWeight = pid.fPHOSPhotonWeight ;
  fPHOSPi0Weight = pid.fPHOSPi0Weight ;
  fPHOSElectronWeight = pid.fPHOSElectronWeight; 
  fPHOSChargeWeight = pid.fPHOSChargeWeight;
  fPHOSNeutralWeight = pid.fPHOSNeutralWeight;
  
  fPHOSWeightFormula       = pid.fPHOSWeightFormula; 
  fPHOSPhotonWeightFormula = pid.fPHOSPhotonWeightFormula; 
  fPHOSPi0WeightFormula    = pid.fPHOSPi0WeightFormula;
  
  fDispCut  = pid.fDispCut;
  fTOFCut   = pid.fTOFCut;
  fDebug    = pid.fDebug;
  
  fRecalculateBayesian = pid.fRecalculateBayesian;
  fParticleFlux        = pid.fParticleFlux;
  fEMCALPIDUtils       = pid.fEMCALPIDUtils;
		
  return *this;
  
}

//_________________________________
AliCaloPID::~AliCaloPID() {
  //Dtor
  
  if(fPHOSPhotonWeightFormula) delete  fPHOSPhotonWeightFormula ;
  if(fPHOSPi0WeightFormula) delete  fPHOSPi0WeightFormula ;
  if(fEMCALPIDUtils) delete fEMCALPIDUtils ;
}


//_______________________________________________________________
void AliCaloPID::InitParameters()
{
  //Initialize the parameters of the PID.
  
  fEMCALPhotonWeight   = 0.5 ;
  fEMCALPi0Weight      = 0.5 ;
  fEMCALElectronWeight = 0.5 ;
  fEMCALChargeWeight   = 0.5 ;
  fEMCALNeutralWeight  = 0.5 ;
  
  fPHOSPhotonWeight    = 0.5 ;
  fPHOSPi0Weight       = 0.5 ;
  fPHOSElectronWeight  = 0.5 ;
  fPHOSChargeWeight    = 0.5 ;
  fPHOSNeutralWeight   = 0.5 ;
  
  //Formula to set the PID weight threshold for photon or pi0
  fPHOSWeightFormula = kTRUE;
  fPHOSPhotonWeightFormula = 
    new TFormula("photonWeight","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))");
  fPHOSPi0WeightFormula = 
    new TFormula("pi0Weight","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))");
  
  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() ;
	}
  }
}

//_______________________________________________________________
Int_t AliCaloPID::GetPdg(const TString calo, const Double_t * pid, const Float_t energy) const {
  //Return most probable identity of the particle.
  
  if(!pid){ 
    printf("AliCaloPID::GetPdg() - pid pointer not initialized!!!\n");
    abort();
  }
  
  Float_t wPh  =  fPHOSPhotonWeight ;
  Float_t wPi0 =  fPHOSPi0Weight ;
  Float_t wE   =  fPHOSElectronWeight ;
  Float_t wCh  =  fPHOSChargeWeight ;
  Float_t wNe  =  fPHOSNeutralWeight ;
  
  
  if(calo == "PHOS" && fPHOSWeightFormula){
    wPh  = fPHOSPhotonWeightFormula->Eval(energy) ;
    wPi0 = fPHOSPi0WeightFormula->Eval(energy);
  }
  
  if(calo == "EMCAL"){
    
    wPh  =  fEMCALPhotonWeight ;
    wPi0 =  fEMCALPi0Weight ;
    wE   =  fEMCALElectronWeight ;
    wCh  =  fEMCALChargeWeight ;
    wNe  =  fEMCALNeutralWeight ;
    
  }
  
  if(fDebug > 0)  printf("AliCaloPID::GetPdg: 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[AliAODCluster::kPhoton], pid[AliAODCluster::kPi0],
			 pid[AliAODCluster::kElectron], pid[AliAODCluster::kEleCon],
			 pid[AliAODCluster::kPion], pid[AliAODCluster::kKaon], pid[AliAODCluster::kProton],
			 pid[AliAODCluster::kNeutron], pid[AliAODCluster::kKaon0]);
  
  Int_t pdg = kNeutralUnknown ;
  Float_t chargedHadronWeight = pid[AliAODCluster::kProton]+pid[AliAODCluster::kKaon]+
    pid[AliAODCluster::kPion]+pid[AliAODCluster::kMuon];
  Float_t neutralHadronWeight = pid[AliAODCluster::kNeutron]+pid[AliAODCluster::kKaon0];
  Float_t allChargedWeight    = pid[AliAODCluster::kElectron]+pid[AliAODCluster::kEleCon]+ chargedHadronWeight;
  Float_t allNeutralWeight    = pid[AliAODCluster::kPhoton]+pid[AliAODCluster::kPi0]+ neutralHadronWeight;
  
  //Select most probable ID
  if(calo=="PHOS"){
    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::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{//EMCAL

    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(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;  
    else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ; 
    else pdg =  kNeutralUnknown ;
    
  }
  
  if(fDebug > 0)printf("AliCaloPID::GetPdg:Final Pdg: %d \n", pdg);
   
  return pdg ;
  
}

//_______________________________________________________________
Int_t AliCaloPID::GetPdg(const TString calo,const TLorentzVector mom, const AliAODCaloCluster * cluster) const {
  //Recalculated PID with all parameters
  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") {
	  //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") {
  //    if(cluster->GetM02()< 0.25) return kPhoton ;
  //    else return  kNeutralHadron ; 
  //  }
  
  return  kNeutralHadron ; 
  
}

//__________________________________________________
TString  AliCaloPID::GetPIDParametersList()  {
  //Put data member values in string to keep in output container
  
  TString parList ; //this will be list of parameters used for this analysis.
  char onePar[255] ;
  sprintf(onePar,"--- AliCaloPID ---\n") ;
  parList+=onePar ;	
  sprintf(onePar,"fDispCut =%2.2f (Cut on dispersion, used in PID evaluation) \n",fDispCut) ;
  parList+=onePar ;
  sprintf(onePar,"fTOFCut  =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
  parList+=onePar ;
  sprintf(onePar,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)\n",fEMCALPhotonWeight) ;
  parList+=onePar ;
  sprintf(onePar,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)\n",fEMCALPi0Weight) ;
  parList+=onePar ;
  sprintf(onePar,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)\n",fEMCALElectronWeight) ;
  parList+=onePar ;
  sprintf(onePar,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)\n",fEMCALChargeWeight) ;
  parList+=onePar ;
  sprintf(onePar,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)\n",fEMCALNeutralWeight) ;
  parList+=onePar ;
  sprintf(onePar,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)\n",fPHOSPhotonWeight) ;
  parList+=onePar ;
  sprintf(onePar,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)\n",fPHOSPi0Weight) ;
  parList+=onePar ;
  sprintf(onePar,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)\n",fPHOSElectronWeight) ;
  parList+=onePar ;
  sprintf(onePar,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)\n",fPHOSChargeWeight) ;
  parList+=onePar ;
  sprintf(onePar,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)\n",fPHOSNeutralWeight) ;
  parList+=onePar ;
  
  if(fPHOSWeightFormula){
	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; 
  
}

//________________________________________________________________
void AliCaloPID::Print(const Option_t * opt) const
{
  
  //Print some relevant parameters set for the analysis
  if(! opt)
    return;
  
  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
  
  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", (fPHOSPhotonWeightFormula->GetExpFormula("p")).Data());
    printf("Pi0    weight formula = %s\n", (fPHOSPhotonWeightFormula->GetExpFormula("p")).Data());
  }
  
  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");
  
} 

//_______________________________________________________________
void AliCaloPID::SetPIDBits(const TString calo, const AliAODCaloCluster * cluster, AliAODPWG4Particle * ph) {
  //Set Bits for PID selection
  
  //Dispersion/lambdas
  Double_t disp=cluster->GetDispersion()  ;
	//    Double_t m20=calo->GetM20() ;
	//    Double_t m02=calo->GetM02() ; 
  ph->SetDispBit(disp<fDispCut) ;  
  
  //TOF
  Double_t tof=cluster->GetTOF()  ;
  ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ; 
  
  //Charged veto
  //    Double_t cpvR=calo->GetEmcCpvDistance() ; 
  Int_t ntr=cluster->GetNTracksMatched();  //number of track matched
  ph->SetChargedBit(ntr>0) ;  //Temporary cut, should we evaluate distance?
  
  //Set PID pdg
  ph->SetPdg(GetPdg(calo,cluster->PID(),ph->E()));
  
  if(fDebug > 0){ 
    printf("AliCaloPID::SetPIDBits: TOF %e, Dispersion %2.2f, NTracks %d\n",tof , disp, ntr); 	
    printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
	     ph->GetPdg(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit()); 
  }
}
Commit	Line	Data
1c5acb87	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15	/* $Id: AliCaloPID.cxx 21839 2007-10-29 13:49:42Z gustavo $ */
	16
	17	//_________________________________________________________________________
bdd2a262	18	// Class for PID selection with calorimeters
	19	// The Output of the 2 main methods GetPdg is a PDG number identifying the cluster,
	20	// being kPhoton, kElectron, kPi0 ... as defined in the header file
	21	// - GetPdg(const TString calo, const Double_t * pid, const Float_t energy)
	22	// Reads the PID weights array of the ESDs and depending on its magnitude identifies the particle
	23	// - GetPdg(const TString calo,const TLorentzVector mom, const AliAODCaloCluster * cluster)
	24	// Recalcultes PID, the bayesian or any new one to be implemented in the future
	25	// Right now only the possibility to recalculate EMCAL with bayesian and simple PID.
	26	// In order to recalculate Bayesian, it is necessary to load the EMCALUtils library
	27	// and do SwitchOnBayesianRecalculation().
	28	// To change the PID parameters from Low to High like the ones by default, use the constructor
	29	// AliCaloPID(flux)
	30	// where flux is AliCaloPID::kLow or AliCaloPID::kHigh
	31	// If it is necessary to change the parameters use the constructor
	32	// AliCaloPID(AliEMCALPIDUtils *utils) and set the parameters before.
	33	// - SetPIDBits: Simple PID, depending on the thresholds fDispCut fTOFCut and even the
	34	// result of the PID bayesian a different PID bit is set.
	35	//
	36	// All these methods can be called in the analysis you are interested.
	37	//
1c5acb87	38	//*-- Author: Gustavo Conesa (LNF-INFN)
	39	//////////////////////////////////////////////////////////////////////////////
	40
	41
	42	// --- ROOT system ---
	43	#include <TMath.h>
1c5acb87	44	#include <TString.h>
	45	#include <TFormula.h>
	46
	47	//---- ANALYSIS system ----
1c5acb87	48	#include "AliCaloPID.h"
	49	#include "AliAODCaloCluster.h"
	50	#include "AliAODPWG4Particle.h"
bdd2a262	51	#include "AliEMCALPIDUtils.h"
1c5acb87	52	ClassImp(AliCaloPID)
	53
	54
	55	//________________________________________________
	56	AliCaloPID::AliCaloPID() :
	57	TObject(), fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
	58	fEMCALElectronWeight(0.), fEMCALChargeWeight(0.),
	59	fEMCALNeutralWeight(0.),
	60	fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
	61	fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) ,
	62	fPHOSNeutralWeight(0.), fPHOSWeightFormula(0),
	63	fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0),
bdd2a262	64	fDispCut(0.),fTOFCut(0.), fDebug(-1),
e5dbdaf0	65	fRecalculateBayesian(kFALSE), fParticleFlux(kLow), fEMCALPIDUtils(new AliEMCALPIDUtils)
1c5acb87	66	{
477d6cee	67	//Ctor
	68
	69	//Initialize parameters
	70	InitParameters();
1c5acb87	71	}
1c5acb87	72
bdd2a262	73	//________________________________________________
	74	AliCaloPID::AliCaloPID(const Int_t flux) :
	75	TObject(), fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
	76	fEMCALElectronWeight(0.), fEMCALChargeWeight(0.),
	77	fEMCALNeutralWeight(0.),
	78	fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
	79	fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) ,
	80	fPHOSNeutralWeight(0.), fPHOSWeightFormula(0),
	81	fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0),
	82	fDispCut(0.),fTOFCut(0.), fDebug(-1),
e5dbdaf0	83	fRecalculateBayesian(kTRUE), fParticleFlux(flux), fEMCALPIDUtils(new AliEMCALPIDUtils)
bdd2a262	84	{
	85	//Ctor
	86
	87	//Initialize parameters
	88	InitParameters();
	89	}
	90
	91	//________________________________________________
	92	AliCaloPID::AliCaloPID(const TTask * emcalpid) :
	93	TObject(), fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
	94	fEMCALElectronWeight(0.), fEMCALChargeWeight(0.),
	95	fEMCALNeutralWeight(0.),
	96	fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
	97	fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) ,
	98	fPHOSNeutralWeight(0.), fPHOSWeightFormula(0),
	99	fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0),
	100	fDispCut(0.),fTOFCut(0.), fDebug(-1),
e5dbdaf0	101	fRecalculateBayesian(kTRUE), fParticleFlux(-1), fEMCALPIDUtils( (AliEMCALPIDUtils*) emcalpid)
bdd2a262	102	{
	103	//Ctor
	104
	105	//Initialize parameters
	106	InitParameters();
	107	}
	108
1c5acb87	109	//____________________________________________________________________________
	110	AliCaloPID::AliCaloPID(const AliCaloPID & pid) :
	111	TObject(pid), fEMCALPhotonWeight(pid.fEMCALPhotonWeight),
	112	fEMCALPi0Weight(pid.fEMCALPi0Weight),
	113	fEMCALElectronWeight(pid.fEMCALElectronWeight),
	114	fEMCALChargeWeight(pid.fEMCALChargeWeight),
	115	fEMCALNeutralWeight(pid.fEMCALNeutralWeight),
	116	fPHOSPhotonWeight(pid.fPHOSPhotonWeight),
	117	fPHOSPi0Weight(pid.fPHOSPi0Weight),
	118	fPHOSElectronWeight(pid.fPHOSElectronWeight),
	119	fPHOSChargeWeight(pid.fPHOSChargeWeight),
	120	fPHOSNeutralWeight(pid.fPHOSNeutralWeight),
	121	fPHOSWeightFormula(pid.fPHOSWeightFormula),
	122	fPHOSPhotonWeightFormula(pid.fPHOSPhotonWeightFormula),
	123	fPHOSPi0WeightFormula(pid.fPHOSPi0WeightFormula),
	124	fDispCut(pid.fDispCut),fTOFCut(pid.fTOFCut),
bdd2a262	125	fDebug(pid.fDebug), fRecalculateBayesian(pid.fRecalculateBayesian),
e5dbdaf0	126	fParticleFlux(pid.fParticleFlux), fEMCALPIDUtils(pid.fEMCALPIDUtils)
1c5acb87	127	{
477d6cee	128	// cpy ctor
1c5acb87	129
	130	}
	131
	132	//_________________________________________________________________________
	133	AliCaloPID & AliCaloPID::operator = (const AliCaloPID & pid)
	134	{
477d6cee	135	// assignment operator
	136
	137	if(&pid == this) return *this;
	138
	139	fEMCALPhotonWeight = pid. fEMCALPhotonWeight ;
	140	fEMCALPi0Weight = pid.fEMCALPi0Weight ;
	141	fEMCALElectronWeight = pid.fEMCALElectronWeight;
	142	fEMCALChargeWeight = pid.fEMCALChargeWeight;
	143	fEMCALNeutralWeight = pid.fEMCALNeutralWeight;
	144
	145	fPHOSPhotonWeight = pid.fPHOSPhotonWeight ;
	146	fPHOSPi0Weight = pid.fPHOSPi0Weight ;
	147	fPHOSElectronWeight = pid.fPHOSElectronWeight;
	148	fPHOSChargeWeight = pid.fPHOSChargeWeight;
	149	fPHOSNeutralWeight = pid.fPHOSNeutralWeight;
	150
	151	fPHOSWeightFormula = pid.fPHOSWeightFormula;
	152	fPHOSPhotonWeightFormula = pid.fPHOSPhotonWeightFormula;
	153	fPHOSPi0WeightFormula = pid.fPHOSPi0WeightFormula;
	154
	155	fDispCut = pid.fDispCut;
	156	fTOFCut = pid.fTOFCut;
	157	fDebug = pid.fDebug;
	158
bdd2a262	159	fRecalculateBayesian = pid.fRecalculateBayesian;
bdd2a262	160	fParticleFlux = pid.fParticleFlux;
bdd2a262	161	fEMCALPIDUtils = pid.fEMCALPIDUtils;
e5dbdaf0	162
477d6cee	163	return *this;
477d6cee	164
1c5acb87	165	}
	166
	167	//_________________________________
	168	AliCaloPID::~AliCaloPID() {
477d6cee	169	//Dtor
	170
	171	if(fPHOSPhotonWeightFormula) delete fPHOSPhotonWeightFormula ;
	172	if(fPHOSPi0WeightFormula) delete fPHOSPi0WeightFormula ;
e5dbdaf0	173	if(fEMCALPIDUtils) delete fEMCALPIDUtils ;
1c5acb87	174	}
	175
	176
1c5acb87	177	//_______________________________________________________________
	178	void AliCaloPID::InitParameters()
	179	{
477d6cee	180	//Initialize the parameters of the PID.
477d6cee	181
691c43c9	182	fEMCALPhotonWeight = 0.5 ;
477d6cee	183	fEMCALPi0Weight = 0.5 ;
691c43c9	184	fEMCALElectronWeight = 0.5 ;
477d6cee	185	fEMCALChargeWeight = 0.5 ;
	186	fEMCALNeutralWeight = 0.5 ;
	187
691c43c9	188	fPHOSPhotonWeight = 0.5 ;
691c43c9	189	fPHOSPi0Weight = 0.5 ;
477d6cee	190	fPHOSElectronWeight = 0.5 ;
	191	fPHOSChargeWeight = 0.5 ;
	192	fPHOSNeutralWeight = 0.5 ;
	193
	194	//Formula to set the PID weight threshold for photon or pi0
	195	fPHOSWeightFormula = kTRUE;
	196	fPHOSPhotonWeightFormula =
1c5acb87	197	new TFormula("photonWeight","0.98(x<40)+ 0.68(x>=100)+(x>=40 && x<100)(0.98+x(6e-3)-xx(2e-04)+xxx*(1.1e-06))");
477d6cee	198	fPHOSPi0WeightFormula =
1c5acb87	199	new TFormula("pi0Weight","0.98(x<65)+ 0.915(x>=100)+(x>=65 && x-x(1.95e-3)-xx(4.31e-05)+xxx(3.61e-07))");
477d6cee	200
	201	fDispCut = 1.5;
	202	fTOFCut = 5.e-9;
	203	fDebug = -1;
bdd2a262	204
bdd2a262	205	if(fRecalculateBayesian){
bdd2a262	206	if(fParticleFlux == kLow){
	207	printf("AliCaloPID::Init() - SetLOWFluxParam\n");
	208	fEMCALPIDUtils->SetLowFluxParam() ;
	209	}
	210	else if (fParticleFlux == kHigh){
	211	printf("AliCaloPID::Init() - SetHIGHFluxParam\n");
	212	fEMCALPIDUtils->SetHighFluxParam() ;
	213	}
bdd2a262	214	}
1c5acb87	215	}
	216
	217	//_______________________________________________________________
	218	Int_t AliCaloPID::GetPdg(const TString calo, const Double_t * pid, const Float_t energy) const {
477d6cee	219	//Return most probable identity of the particle.
	220
	221	if(!pid){
	222	printf("AliCaloPID::GetPdg() - pid pointer not initialized!!!\n");
	223	abort();
	224	}
	225
15800db4	226	Float_t wPh = fPHOSPhotonWeight ;
477d6cee	227	Float_t wPi0 = fPHOSPi0Weight ;
15800db4	228	Float_t wE = fPHOSElectronWeight ;
	229	Float_t wCh = fPHOSChargeWeight ;
	230	Float_t wNe = fPHOSNeutralWeight ;
477d6cee	231
	232
	233	if(calo == "PHOS" && fPHOSWeightFormula){
	234	wPh = fPHOSPhotonWeightFormula->Eval(energy) ;
	235	wPi0 = fPHOSPi0WeightFormula->Eval(energy);
	236	}
	237
	238	if(calo == "EMCAL"){
	239
	240	wPh = fEMCALPhotonWeight ;
	241	wPi0 = fEMCALPi0Weight ;
	242	wE = fEMCALElectronWeight ;
	243	wCh = fEMCALChargeWeight ;
	244	wNe = fEMCALNeutralWeight ;
	245
	246	}
	247
	248	if(fDebug > 0) printf("AliCaloPID::GetPdg: 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",
	249	calo.Data(),pid[AliAODCluster::kPhoton], pid[AliAODCluster::kPi0],
	250	pid[AliAODCluster::kElectron], pid[AliAODCluster::kEleCon],
	251	pid[AliAODCluster::kPion], pid[AliAODCluster::kKaon], pid[AliAODCluster::kProton],
	252	pid[AliAODCluster::kNeutron], pid[AliAODCluster::kKaon0]);
	253
	254	Int_t pdg = kNeutralUnknown ;
	255	Float_t chargedHadronWeight = pid[AliAODCluster::kProton]+pid[AliAODCluster::kKaon]+
1c5acb87	256	pid[AliAODCluster::kPion]+pid[AliAODCluster::kMuon];
477d6cee	257	Float_t neutralHadronWeight = pid[AliAODCluster::kNeutron]+pid[AliAODCluster::kKaon0];
	258	Float_t allChargedWeight = pid[AliAODCluster::kElectron]+pid[AliAODCluster::kEleCon]+ chargedHadronWeight;
	259	Float_t allNeutralWeight = pid[AliAODCluster::kPhoton]+pid[AliAODCluster::kPi0]+ neutralHadronWeight;
	260
	261	//Select most probable ID
	262	if(calo=="PHOS"){
	263	if(pid[AliAODCluster::kPhoton] > wPh) pdg = kPhoton ;
	264	else if(pid[AliAODCluster::kPi0] > wPi0) pdg = kPi0 ;
	265	else if(pid[AliAODCluster::kElectron] > wE) pdg = kElectron ;
	266	else if(pid[AliAODCluster::kEleCon] > wE) pdg = kEleCon ;
	267	else if(chargedHadronWeight > wCh) pdg = kChargedHadron ;
	268	else if(neutralHadronWeight > wNe) pdg = kNeutralHadron ;
	269	else if(allChargedWeight > allNeutralWeight)
	270	pdg = kChargedUnknown ;
	271	else
	272	pdg = kNeutralUnknown ;
	273	}
	274	else{//EMCAL
bdd2a262	275
15800db4	276	if(pid[AliAODCluster::kPhoton]+pid[AliAODCluster::kElectron] > wPh) pdg = kPhoton ; //temporal sollution until track matching for electrons is considered
15800db4	277	//if(pid[AliAODCluster::kPhoton] > wPh) pdg = kPhoton ;
477d6cee	278	else if(pid[AliAODCluster::kPi0] > wPi0) pdg = kPi0 ;
15800db4	279	//else if(pid[AliAODCluster::kElectron] > wE) pdg = kElectron ;
477d6cee	280	else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;
	281	else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ;
	282	else pdg = kNeutralUnknown ;
	283
	284	}
	285
	286	if(fDebug > 0)printf("AliCaloPID::GetPdg:Final Pdg: %d \n", pdg);
	287
	288	return pdg ;
	289
1c5acb87	290	}
	291
	292	//_______________________________________________________________
	293	Int_t AliCaloPID::GetPdg(const TString calo,const TLorentzVector mom, const AliAODCaloCluster * cluster) const {
477d6cee	294	//Recalculated PID with all parameters
bdd2a262	295	Float_t lambda0 = cluster->GetM02();
	296	Float_t energy = mom.E();
	297
	298	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",
	299	calo.Data(),energy,lambda0,cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(),
	300	cluster->GetEmcCpvDistance(), cluster->GetDistToBadChannel(),cluster->GetNExMax());
	301
477d6cee	302	if(calo == "EMCAL") {
bdd2a262	303	//Recalculate Bayesian
bdd2a262	304	if(fRecalculateBayesian){
bdd2a262	305	if(fDebug > 0) {
	306	const Double_t *pid0 = cluster->PID();
	307	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",
	308	calo.Data(),pid0[AliAODCluster::kPhoton], pid0[AliAODCluster::kPi0],
	309	pid0[AliAODCluster::kElectron], pid0[AliAODCluster::kEleCon],
	310	pid0[AliAODCluster::kPion], pid0[AliAODCluster::kKaon], pid0[AliAODCluster::kProton],
	311	pid0[AliAODCluster::kNeutron], pid0[AliAODCluster::kKaon0]);
	312	}
	313
	314	fEMCALPIDUtils->ComputePID(energy, lambda0);
	315	Double_t pid[AliPID::kSPECIESN];
	316	for(Int_t i = 0; i < AliPID::kSPECIESN; i++) pid[i] = fEMCALPIDUtils->GetPIDFinal(i);
	317	return GetPdg(calo, pid, energy);
	318
bdd2a262	319
	320	}
	321
	322	// If no use of bayesian, simple PID
	323	if(lambda0 < 0.25) return kPhoton ;
bb5c272c	324	//else return kNeutralHadron ;
bb5c272c	325	else return kPi0 ;
477d6cee	326	}
	327
	328	// if(calo == "PHOS") {
	329	// if(cluster->GetM02()< 0.25) return kPhoton ;
	330	// else return kNeutralHadron ;
	331	// }
	332
	333	return kNeutralHadron ;
	334
1c5acb87	335	}
	336
	337	//__________________________________________________
	338	TString AliCaloPID::GetPIDParametersList() {
477d6cee	339	//Put data member values in string to keep in output container
	340
	341	TString parList ; //this will be list of parameters used for this analysis.
	342	char onePar[255] ;
	343	sprintf(onePar,"--- AliCaloPID ---\n") ;
	344	parList+=onePar ;
	345	sprintf(onePar,"fDispCut =%2.2f (Cut on dispersion, used in PID evaluation) \n",fDispCut) ;
	346	parList+=onePar ;
	347	sprintf(onePar,"fTOFCut =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
	348	parList+=onePar ;
	349	sprintf(onePar,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)\n",fEMCALPhotonWeight) ;
	350	parList+=onePar ;
	351	sprintf(onePar,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)\n",fEMCALPi0Weight) ;
	352	parList+=onePar ;
	353	sprintf(onePar,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)\n",fEMCALElectronWeight) ;
	354	parList+=onePar ;
	355	sprintf(onePar,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)\n",fEMCALChargeWeight) ;
	356	parList+=onePar ;
	357	sprintf(onePar,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)\n",fEMCALNeutralWeight) ;
	358	parList+=onePar ;
	359	sprintf(onePar,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)\n",fPHOSPhotonWeight) ;
	360	parList+=onePar ;
	361	sprintf(onePar,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)\n",fPHOSPi0Weight) ;
	362	parList+=onePar ;
	363	sprintf(onePar,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)\n",fPHOSElectronWeight) ;
	364	parList+=onePar ;
	365	sprintf(onePar,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)\n",fPHOSChargeWeight) ;
	366	parList+=onePar ;
	367	sprintf(onePar,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)\n",fPHOSNeutralWeight) ;
	368	parList+=onePar ;
	369
	370	if(fPHOSWeightFormula){
62ccac92	371	sprintf(onePar,"PHOS Photon Weight Formula: %s\n",(fPHOSPhotonWeightFormula->GetExpFormula("p")).Data()) ;
	372	parList+=onePar;
	373	sprintf(onePar,"PHOS Pi0 Weight Formula: %s\n",(fPHOSPi0WeightFormula->GetExpFormula("p")).Data()) ;
	374	parList+=onePar;
477d6cee	375	}
	376
	377	return parList;
	378
1c5acb87	379	}
	380
	381	//________________________________________________________________
	382	void AliCaloPID::Print(const Option_t * opt) const
	383	{
477d6cee	384
	385	//Print some relevant parameters set for the analysis
	386	if(! opt)
	387	return;
	388
	389	printf("*** Print: %s %s ****\n", GetName(), GetTitle() ) ;
	390
	391	printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",
	392	fPHOSPhotonWeight, fPHOSPi0Weight,
	393	fPHOSElectronWeight, fPHOSChargeWeight, fPHOSNeutralWeight) ;
	394	printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",
	395	fEMCALPhotonWeight, fEMCALPi0Weight,
	396	fEMCALElectronWeight, fEMCALChargeWeight, fEMCALNeutralWeight) ;
	397
	398	printf("PHOS Parametrized weight on? = %d\n", fPHOSWeightFormula) ;
	399	if(fPHOSWeightFormula){
	400	printf("Photon weight formula = %s\n", (fPHOSPhotonWeightFormula->GetExpFormula("p")).Data());
	401	printf("Pi0 weight formula = %s\n", (fPHOSPhotonWeightFormula->GetExpFormula("p")).Data());
	402	}
	403
	404	printf("TOF cut = %e\n",fTOFCut);
	405	printf("Dispersion cut = %2.2f\n",fDispCut);
	406	printf("Debug level = %d\n",fDebug);
bdd2a262	407	printf("Recalculate Bayesian? = %d\n",fRecalculateBayesian);
bdd2a262	408	if(fRecalculateBayesian) printf("Particle Flux? = %d\n",fParticleFlux);
477d6cee	409	printf(" \n");
477d6cee	410
1c5acb87	411	}
	412
	413	//_______________________________________________________________
	414	void AliCaloPID::SetPIDBits(const TString calo, const AliAODCaloCluster * cluster, AliAODPWG4Particle * ph) {
477d6cee	415	//Set Bits for PID selection
	416
	417	//Dispersion/lambdas
	418	Double_t disp=cluster->GetDispersion() ;
1c5acb87	419	// Double_t m20=calo->GetM20() ;
1c5acb87	420	// Double_t m02=calo->GetM02() ;
477d6cee	421	ph->SetDispBit(disp<fDispCut) ;
	422
	423	//TOF
	424	Double_t tof=cluster->GetTOF() ;
	425	ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
	426
	427	//Charged veto
	428	// Double_t cpvR=calo->GetEmcCpvDistance() ;
	429	Int_t ntr=cluster->GetNTracksMatched(); //number of track matched
	430	ph->SetChargedBit(ntr>0) ; //Temporary cut, should we evaluate distance?
	431
	432	//Set PID pdg
	433	ph->SetPdg(GetPdg(calo,cluster->PID(),ph->E()));
	434
	435	if(fDebug > 0){
	436	printf("AliCaloPID::SetPIDBits: TOF %e, Dispersion %2.2f, NTracks %d\n",tof , disp, ntr);
	437	printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
6639984f	438	ph->GetPdg(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit());
477d6cee	439	}
1c5acb87	440	}
	441
	442