[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"
#ifdef __EMCALUTIL__
#include "AliEMCALPIDUtils.h"
#endif
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)
#ifdef __EMCALUTIL__
, fEMCALPIDUtils(new AliEMCALPIDUtils)
#endif
{
  //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)
#ifdef __EMCALUTIL__
, fEMCALPIDUtils(new AliEMCALPIDUtils)
#endif
{
	//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)
#ifdef __EMCALUTIL__
, fEMCALPIDUtils( (AliEMCALPIDUtils*) emcalpid)
#endif
{
	//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)
#ifdef __EMCALUTIL__
, fEMCALPIDUtils(pid.fEMCALPIDUtils)
#endif
{
  // 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;
#ifdef __EMCALUTIL__
  fEMCALPIDUtils       = pid.fEMCALPIDUtils;
#endif
	
	
  return *this;
  
}

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


//_______________________________________________________________
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){
#ifdef __EMCALUTIL__
	if(fParticleFlux == kLow){
		printf("AliCaloPID::Init() - SetLOWFluxParam\n");
		fEMCALPIDUtils->SetLowFluxParam() ;
	}
	else if (fParticleFlux == kHigh){
		printf("AliCaloPID::Init() - SetHIGHFluxParam\n");
		fEMCALPIDUtils->SetHighFluxParam() ;
	}
#endif
  }
}

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