[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(new TFormula(*pid.fPHOSPhotonWeightFormula)), 
//fPHOSPi0WeightFormula   (new TFormula(*pid.fPHOSPi0WeightFormula)), 
fDispCut(pid.fDispCut),fTOFCut(pid.fTOFCut),
fDebug(pid.fDebug), fRecalculateBayesian(pid.fRecalculateBayesian),
fParticleFlux(pid.fParticleFlux), 
fEMCALPIDUtils(new AliEMCALPIDUtils(*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; 
//  delete fPHOSPhotonWeightFormula;  fPHOSPhotonWeightFormula = new TFormula(*pid.fPHOSPhotonWeightFormula); 
//  delete fPHOSPi0WeightFormula;     fPHOSPi0WeightFormula    = new TFormula(*pid.fPHOSPi0WeightFormula);
//  
//  fDispCut  = pid.fDispCut;
//  fTOFCut   = pid.fTOFCut;
//  fDebug    = pid.fDebug;
//  
//  fRecalculateBayesian = pid.fRecalculateBayesian;
//  fParticleFlux        = pid.fParticleFlux;
//	
//  delete fEMCALPIDUtils; fEMCALPIDUtils       = new AliEMCALPIDUtils(*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;
  //if(fPHOSPhotonWeightFormula) delete fPHOSPhotonWeightFormula;
  //if(fPHOSPi0WeightFormula)    delete fPHOSPi0WeightFormula;
  //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[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]);
  
  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;
  
  //Select most probable ID
  if(calo=="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 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{//EMCAL

    if(pid[AliVCluster::kPhoton]+pid[AliVCluster::kElectron]  > wPh) pdg = kPhoton ; //temporal sollution until track matching for electrons is considered
    //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(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;  
    else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ; 
    else pdg =  kNeutralUnknown ;
    
  }
  
  if(fDebug > 0)printf("AliCaloPID::GetPdg:Final Pdg: %d, cluster energy %2.2f \n", pdg,energy);
   
  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->GetDistanceToBadChannel(),cluster->GetNExMax());

  if(calo == "EMCAL") {
	  //Recalculate Bayesian
	  if(fRecalculateBayesian){	  
		  if(fDebug > 0)  {
			  const Double_t  *pid0 = cluster->GetPID();
			  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[AliVCluster::kPhoton], pid0[AliVCluster::kPi0],
								 pid0[AliVCluster::kElectron], pid0[AliVCluster::kEleCon],
								 pid0[AliVCluster::kPion], pid0[AliVCluster::kKaon], pid0[AliVCluster::kProton],
								 pid0[AliVCluster::kNeutron], pid0[AliVCluster::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->GetPID(),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>
a33b76ed	45	//#include <TFormula.h>
1c5acb87	46
1c5acb87	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.) ,
a33b76ed	62	fPHOSNeutralWeight(0.), //fPHOSWeightFormula(0),
a33b76ed	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.) ,
a33b76ed	80	fPHOSNeutralWeight(0.), //fPHOSWeightFormula(0),
a33b76ed	81	//fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0),
bdd2a262	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.) ,
a33b76ed	98	fPHOSNeutralWeight(0.), //fPHOSWeightFormula(0),
a33b76ed	99	//fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0),
bdd2a262	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),
a33b76ed	121	//fPHOSWeightFormula(pid.fPHOSWeightFormula),
	122	//fPHOSPhotonWeightFormula(new TFormula(*pid.fPHOSPhotonWeightFormula)),
	123	//fPHOSPi0WeightFormula (new TFormula(*pid.fPHOSPi0WeightFormula)),
1c5acb87	124	fDispCut(pid.fDispCut),fTOFCut(pid.fTOFCut),
bdd2a262	125	fDebug(pid.fDebug), fRecalculateBayesian(pid.fRecalculateBayesian),
fb9841b8	126	fParticleFlux(pid.fParticleFlux),
fb9841b8	127	fEMCALPIDUtils(new AliEMCALPIDUtils(*pid.fEMCALPIDUtils))
1c5acb87	128	{
477d6cee	129	// cpy ctor
1c5acb87	130
	131	}
	132
	133	//_________________________________________________________________________
614701c6	134	//AliCaloPID & AliCaloPID::operator = (const AliCaloPID & pid)
	135	//{
	136	// // assignment operator
	137	//
	138	// if(&pid == this) return *this;
	139	//
	140	// fEMCALPhotonWeight = pid. fEMCALPhotonWeight ;
	141	// fEMCALPi0Weight = pid.fEMCALPi0Weight ;
	142	// fEMCALElectronWeight = pid.fEMCALElectronWeight;
	143	// fEMCALChargeWeight = pid.fEMCALChargeWeight;
	144	// fEMCALNeutralWeight = pid.fEMCALNeutralWeight;
	145	//
	146	// fPHOSPhotonWeight = pid.fPHOSPhotonWeight ;
	147	// fPHOSPi0Weight = pid.fPHOSPi0Weight ;
	148	// fPHOSElectronWeight = pid.fPHOSElectronWeight;
	149	// fPHOSChargeWeight = pid.fPHOSChargeWeight;
	150	// fPHOSNeutralWeight = pid.fPHOSNeutralWeight;
	151	//
	152	// fPHOSWeightFormula = pid.fPHOSWeightFormula;
	153	// delete fPHOSPhotonWeightFormula; fPHOSPhotonWeightFormula = new TFormula(*pid.fPHOSPhotonWeightFormula);
	154	// delete fPHOSPi0WeightFormula; fPHOSPi0WeightFormula = new TFormula(*pid.fPHOSPi0WeightFormula);
	155	//
	156	// fDispCut = pid.fDispCut;
	157	// fTOFCut = pid.fTOFCut;
	158	// fDebug = pid.fDebug;
	159	//
	160	// fRecalculateBayesian = pid.fRecalculateBayesian;
	161	// fParticleFlux = pid.fParticleFlux;
	162	//
	163	// delete fEMCALPIDUtils; fEMCALPIDUtils = new AliEMCALPIDUtils(*pid.fEMCALPIDUtils);
	164	//
	165	// return *this;
	166	//
	167	//}
1c5acb87	168
	169	//_________________________________
	170	AliCaloPID::~AliCaloPID() {
477d6cee	171	//Dtor
477d6cee	172
a33b76ed	173	// if(fPHOSPhotonWeightFormula) delete fPHOSPhotonWeightFormula ;
a33b76ed	174	// if(fPHOSPi0WeightFormula) delete fPHOSPi0WeightFormula ;
fb9841b8	175	if(fEMCALPIDUtils) delete fEMCALPIDUtils ;
1c5acb87	176	}
	177
	178
1c5acb87	179	//_______________________________________________________________
	180	void AliCaloPID::InitParameters()
	181	{
477d6cee	182	//Initialize the parameters of the PID.
477d6cee	183
691c43c9	184	fEMCALPhotonWeight = 0.5 ;
477d6cee	185	fEMCALPi0Weight = 0.5 ;
691c43c9	186	fEMCALElectronWeight = 0.5 ;
477d6cee	187	fEMCALChargeWeight = 0.5 ;
	188	fEMCALNeutralWeight = 0.5 ;
	189
691c43c9	190	fPHOSPhotonWeight = 0.5 ;
691c43c9	191	fPHOSPi0Weight = 0.5 ;
477d6cee	192	fPHOSElectronWeight = 0.5 ;
	193	fPHOSChargeWeight = 0.5 ;
	194	fPHOSNeutralWeight = 0.5 ;
	195
	196	//Formula to set the PID weight threshold for photon or pi0
a33b76ed	197	//fPHOSWeightFormula = kTRUE;
	198	//if(fPHOSPhotonWeightFormula) delete fPHOSPhotonWeightFormula;
	199	//if(fPHOSPi0WeightFormula) delete fPHOSPi0WeightFormula;
	200	//fPHOSPhotonWeightFormula =
	201	//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))");
	202	//fPHOSPi0WeightFormula =
	203	//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	204
	205	fDispCut = 1.5;
	206	fTOFCut = 5.e-9;
	207	fDebug = -1;
bdd2a262	208
bdd2a262	209	if(fRecalculateBayesian){
bdd2a262	210	if(fParticleFlux == kLow){
	211	printf("AliCaloPID::Init() - SetLOWFluxParam\n");
	212	fEMCALPIDUtils->SetLowFluxParam() ;
	213	}
	214	else if (fParticleFlux == kHigh){
	215	printf("AliCaloPID::Init() - SetHIGHFluxParam\n");
	216	fEMCALPIDUtils->SetHighFluxParam() ;
	217	}
bdd2a262	218	}
1c5acb87	219	}
	220
	221	//_______________________________________________________________
	222	Int_t AliCaloPID::GetPdg(const TString calo, const Double_t * pid, const Float_t energy) const {
477d6cee	223	//Return most probable identity of the particle.
	224
	225	if(!pid){
	226	printf("AliCaloPID::GetPdg() - pid pointer not initialized!!!\n");
	227	abort();
	228	}
	229
15800db4	230	Float_t wPh = fPHOSPhotonWeight ;
477d6cee	231	Float_t wPi0 = fPHOSPi0Weight ;
15800db4	232	Float_t wE = fPHOSElectronWeight ;
	233	Float_t wCh = fPHOSChargeWeight ;
	234	Float_t wNe = fPHOSNeutralWeight ;
477d6cee	235
477d6cee	236
a33b76ed	237	// if(calo == "PHOS" && fPHOSWeightFormula){
	238	// wPh = fPHOSPhotonWeightFormula->Eval(energy) ;
	239	// wPi0 = fPHOSPi0WeightFormula->Eval(energy);
	240	// }
477d6cee	241
	242	if(calo == "EMCAL"){
	243
	244	wPh = fEMCALPhotonWeight ;
	245	wPi0 = fEMCALPi0Weight ;
	246	wE = fEMCALElectronWeight ;
	247	wCh = fEMCALChargeWeight ;
	248	wNe = fEMCALNeutralWeight ;
	249
	250	}
	251
	252	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",
c8fe2783	253	calo.Data(),pid[AliVCluster::kPhoton], pid[AliVCluster::kPi0],
	254	pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon],
	255	pid[AliVCluster::kPion], pid[AliVCluster::kKaon], pid[AliVCluster::kProton],
	256	pid[AliVCluster::kNeutron], pid[AliVCluster::kKaon0]);
477d6cee	257
477d6cee	258	Int_t pdg = kNeutralUnknown ;
c8fe2783	259	Float_t chargedHadronWeight = pid[AliVCluster::kProton]+pid[AliVCluster::kKaon]+
	260	pid[AliVCluster::kPion]+pid[AliVCluster::kMuon];
	261	Float_t neutralHadronWeight = pid[AliVCluster::kNeutron]+pid[AliVCluster::kKaon0];
	262	Float_t allChargedWeight = pid[AliVCluster::kElectron]+pid[AliVCluster::kEleCon]+ chargedHadronWeight;
	263	Float_t allNeutralWeight = pid[AliVCluster::kPhoton]+pid[AliVCluster::kPi0]+ neutralHadronWeight;
477d6cee	264
	265	//Select most probable ID
	266	if(calo=="PHOS"){
c8fe2783	267	if(pid[AliVCluster::kPhoton] > wPh) pdg = kPhoton ;
	268	else if(pid[AliVCluster::kPi0] > wPi0) pdg = kPi0 ;
	269	else if(pid[AliVCluster::kElectron] > wE) pdg = kElectron ;
	270	else if(pid[AliVCluster::kEleCon] > wE) pdg = kEleCon ;
477d6cee	271	else if(chargedHadronWeight > wCh) pdg = kChargedHadron ;
	272	else if(neutralHadronWeight > wNe) pdg = kNeutralHadron ;
	273	else if(allChargedWeight > allNeutralWeight)
	274	pdg = kChargedUnknown ;
	275	else
	276	pdg = kNeutralUnknown ;
	277	}
	278	else{//EMCAL
bdd2a262	279
c8fe2783	280	if(pid[AliVCluster::kPhoton]+pid[AliVCluster::kElectron] > wPh) pdg = kPhoton ; //temporal sollution until track matching for electrons is considered
	281	//if(pid[AliVCluster::kPhoton] > wPh) pdg = kPhoton ;
	282	else if(pid[AliVCluster::kPi0] > wPi0) pdg = kPi0 ;
	283	//else if(pid[AliVCluster::kElectron] > wE) pdg = kElectron ;
477d6cee	284	else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;
	285	else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ;
	286	else pdg = kNeutralUnknown ;
	287
	288	}
	289
ebc881ec	290	if(fDebug > 0)printf("AliCaloPID::GetPdg:Final Pdg: %d, cluster energy %2.2f \n", pdg,energy);
477d6cee	291
	292	return pdg ;
	293
1c5acb87	294	}
	295
	296	//_______________________________________________________________
	297	Int_t AliCaloPID::GetPdg(const TString calo,const TLorentzVector mom, const AliAODCaloCluster * cluster) const {
477d6cee	298	//Recalculated PID with all parameters
bdd2a262	299	Float_t lambda0 = cluster->GetM02();
	300	Float_t energy = mom.E();
	301
	302	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",
	303	calo.Data(),energy,lambda0,cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(),
c8fe2783	304	cluster->GetEmcCpvDistance(), cluster->GetDistanceToBadChannel(),cluster->GetNExMax());
bdd2a262	305
477d6cee	306	if(calo == "EMCAL") {
bdd2a262	307	//Recalculate Bayesian
bdd2a262	308	if(fRecalculateBayesian){
bdd2a262	309	if(fDebug > 0) {
c8fe2783	310	const Double_t *pid0 = cluster->GetPID();
bdd2a262	311	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",
c8fe2783	312	calo.Data(),pid0[AliVCluster::kPhoton], pid0[AliVCluster::kPi0],
	313	pid0[AliVCluster::kElectron], pid0[AliVCluster::kEleCon],
	314	pid0[AliVCluster::kPion], pid0[AliVCluster::kKaon], pid0[AliVCluster::kProton],
	315	pid0[AliVCluster::kNeutron], pid0[AliVCluster::kKaon0]);
bdd2a262	316	}
	317
	318	fEMCALPIDUtils->ComputePID(energy, lambda0);
	319	Double_t pid[AliPID::kSPECIESN];
	320	for(Int_t i = 0; i < AliPID::kSPECIESN; i++) pid[i] = fEMCALPIDUtils->GetPIDFinal(i);
	321	return GetPdg(calo, pid, energy);
	322
bdd2a262	323
	324	}
	325
c8fe2783	326	// If no use of bayesian, simple PID
	327	if(lambda0 < 0.25) return kPhoton ;
	328	//else return kNeutralHadron ;
	329	else return kPi0 ;
477d6cee	330	}
	331
	332	// if(calo == "PHOS") {
	333	// if(cluster->GetM02()< 0.25) return kPhoton ;
	334	// else return kNeutralHadron ;
	335	// }
	336
	337	return kNeutralHadron ;
	338
1c5acb87	339	}
	340
	341	//__________________________________________________
	342	TString AliCaloPID::GetPIDParametersList() {
477d6cee	343	//Put data member values in string to keep in output container
	344
	345	TString parList ; //this will be list of parameters used for this analysis.
	346	char onePar[255] ;
	347	sprintf(onePar,"--- AliCaloPID ---\n") ;
	348	parList+=onePar ;
	349	sprintf(onePar,"fDispCut =%2.2f (Cut on dispersion, used in PID evaluation) \n",fDispCut) ;
	350	parList+=onePar ;
	351	sprintf(onePar,"fTOFCut =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
	352	parList+=onePar ;
	353	sprintf(onePar,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)\n",fEMCALPhotonWeight) ;
	354	parList+=onePar ;
	355	sprintf(onePar,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)\n",fEMCALPi0Weight) ;
	356	parList+=onePar ;
	357	sprintf(onePar,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)\n",fEMCALElectronWeight) ;
	358	parList+=onePar ;
	359	sprintf(onePar,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)\n",fEMCALChargeWeight) ;
	360	parList+=onePar ;
	361	sprintf(onePar,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)\n",fEMCALNeutralWeight) ;
	362	parList+=onePar ;
	363	sprintf(onePar,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)\n",fPHOSPhotonWeight) ;
	364	parList+=onePar ;
	365	sprintf(onePar,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)\n",fPHOSPi0Weight) ;
	366	parList+=onePar ;
	367	sprintf(onePar,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)\n",fPHOSElectronWeight) ;
	368	parList+=onePar ;
	369	sprintf(onePar,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)\n",fPHOSChargeWeight) ;
	370	parList+=onePar ;
	371	sprintf(onePar,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)\n",fPHOSNeutralWeight) ;
	372	parList+=onePar ;
	373
a33b76ed	374	// if(fPHOSWeightFormula){
	375	// sprintf(onePar,"PHOS Photon Weight Formula: %s\n",(fPHOSPhotonWeightFormula->GetExpFormula("p")).Data()) ;
	376	// parList+=onePar;
	377	// sprintf(onePar,"PHOS Pi0 Weight Formula: %s\n",(fPHOSPi0WeightFormula->GetExpFormula("p")).Data()) ;
	378	// parList+=onePar;
	379	// }
477d6cee	380
	381	return parList;
	382
1c5acb87	383	}
	384
	385	//________________________________________________________________
	386	void AliCaloPID::Print(const Option_t * opt) const
	387	{
477d6cee	388
	389	//Print some relevant parameters set for the analysis
	390	if(! opt)
	391	return;
	392
	393	printf("*** Print: %s %s ****\n", GetName(), GetTitle() ) ;
	394
	395	printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",
	396	fPHOSPhotonWeight, fPHOSPi0Weight,
	397	fPHOSElectronWeight, fPHOSChargeWeight, fPHOSNeutralWeight) ;
	398	printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",
	399	fEMCALPhotonWeight, fEMCALPi0Weight,
	400	fEMCALElectronWeight, fEMCALChargeWeight, fEMCALNeutralWeight) ;
	401
a33b76ed	402	// printf("PHOS Parametrized weight on? = %d\n", fPHOSWeightFormula) ;
	403	// if(fPHOSWeightFormula){
	404	// printf("Photon weight formula = %s\n", (fPHOSPhotonWeightFormula->GetExpFormula("p")).Data());
	405	// printf("Pi0 weight formula = %s\n", (fPHOSPhotonWeightFormula->GetExpFormula("p")).Data());
	406	// }
477d6cee	407
	408	printf("TOF cut = %e\n",fTOFCut);
	409	printf("Dispersion cut = %2.2f\n",fDispCut);
	410	printf("Debug level = %d\n",fDebug);
bdd2a262	411	printf("Recalculate Bayesian? = %d\n",fRecalculateBayesian);
bdd2a262	412	if(fRecalculateBayesian) printf("Particle Flux? = %d\n",fParticleFlux);
477d6cee	413	printf(" \n");
477d6cee	414
1c5acb87	415	}
	416
	417	//_______________________________________________________________
	418	void AliCaloPID::SetPIDBits(const TString calo, const AliAODCaloCluster * cluster, AliAODPWG4Particle * ph) {
477d6cee	419	//Set Bits for PID selection
	420
	421	//Dispersion/lambdas
	422	Double_t disp=cluster->GetDispersion() ;
1c5acb87	423	// Double_t m20=calo->GetM20() ;
1c5acb87	424	// Double_t m02=calo->GetM02() ;
477d6cee	425	ph->SetDispBit(disp<fDispCut) ;
	426
	427	//TOF
	428	Double_t tof=cluster->GetTOF() ;
	429	ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
	430
	431	//Charged veto
	432	// Double_t cpvR=calo->GetEmcCpvDistance() ;
	433	Int_t ntr=cluster->GetNTracksMatched(); //number of track matched
	434	ph->SetChargedBit(ntr>0) ; //Temporary cut, should we evaluate distance?
	435
	436	//Set PID pdg
c8fe2783	437	ph->SetPdg(GetPdg(calo,cluster->GetPID(),ph->E()));
477d6cee	438
	439	if(fDebug > 0){
	440	printf("AliCaloPID::SetPIDBits: TOF %e, Dispersion %2.2f, NTracks %d\n",tof , disp, ntr);
	441	printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
6639984f	442	ph->GetPdg(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit());
477d6cee	443	}
1c5acb87	444	}
	445
	446