[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 track/cluster acceptance selection
// Selection in Central barrel, EMCAL and PHOS
//                
//*-- Author: Gustavo Conesa (LNF-INFN) 
//////////////////////////////////////////////////////////////////////////////


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

//---- ANALYSIS system ----
#include "AliLog.h"
#include "AliCaloPID.h"
#include "AliAODCaloCluster.h"
#include "AliAODPWG4Particle.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)
{
	//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)
{
	// 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;
	
	return *this;
	
}

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


//_______________________________________________________________
void AliCaloPID::InitParameters()
{
	//Initialize the parameters of the PID.
	
	fEMCALPhotonWeight   = 0.8 ;
	fEMCALPi0Weight      = 0.5 ;
	fEMCALElectronWeight = 0.8 ;
	fEMCALChargeWeight   = 0.5 ;
	fEMCALNeutralWeight  = 0.5 ;
	
	fPHOSPhotonWeight    = 0.75 ;
	fPHOSPi0Weight       = 0.8 ;
	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;
}

//_______________________________________________________________
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) AliFatal("pid pointer not initialized!!!");
	
	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]  > 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
	if(fDebug > 0)printf("AliCaloPID::GetPdg: Calorimeter %s, E %3.2f, l0 %3.2f, l1 %3.2f, disp %3.2f, tof %1.11f, distCPV %3.2f, distToBC %1.1f, NMax %d\n",
						 calo.Data(),mom.E(),cluster->GetM02(),cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(), 
						 cluster->GetEmcCpvDistance(), cluster->GetDistToBadChannel(),cluster->GetNExMax());
	
	if(calo == "EMCAL") {
		if(cluster->GetM02()< 0.25) return kPhoton ;
		else return  kNeutralHadron ; 
	}
	
	//   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){
		parList+="PHOS Photon Weight Formula: "+(fPHOSPhotonWeightFormula->GetExpFormula("p"));
		parList+="PHOS Pi0    Weight Formula: "+(fPHOSPi0WeightFormula->GetExpFormula("p"));
	}
	
	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(" \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	//_________________________________________________________________________
	18	// Class for track/cluster acceptance selection
	19	// Selection in Central barrel, EMCAL and PHOS
	20	//
	21	//*-- Author: Gustavo Conesa (LNF-INFN)
	22	//////////////////////////////////////////////////////////////////////////////
	23
	24
	25	// --- ROOT system ---
	26	#include <TMath.h>
	27	#include <TLorentzVector.h>
	28	#include <TString.h>
	29	#include <TFormula.h>
	30
	31	//---- ANALYSIS system ----
	32	#include "AliLog.h"
	33	#include "AliCaloPID.h"
	34	#include "AliAODCaloCluster.h"
	35	#include "AliAODPWG4Particle.h"
1c5acb87	36
	37	ClassImp(AliCaloPID)
	38
	39
	40	//________________________________________________
	41	AliCaloPID::AliCaloPID() :
	42	TObject(), fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
	43	fEMCALElectronWeight(0.), fEMCALChargeWeight(0.),
	44	fEMCALNeutralWeight(0.),
	45	fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
	46	fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) ,
	47	fPHOSNeutralWeight(0.), fPHOSWeightFormula(0),
	48	fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0),
6639984f	49	fDispCut(0.),fTOFCut(0.), fDebug(-1)
1c5acb87	50	{
	51	//Ctor
	52
	53	//Initialize parameters
	54	InitParameters();
	55	}
	56
	57	//____________________________________________________________________________
	58	AliCaloPID::AliCaloPID(const AliCaloPID & pid) :
	59	TObject(pid), fEMCALPhotonWeight(pid.fEMCALPhotonWeight),
	60	fEMCALPi0Weight(pid.fEMCALPi0Weight),
	61	fEMCALElectronWeight(pid.fEMCALElectronWeight),
	62	fEMCALChargeWeight(pid.fEMCALChargeWeight),
	63	fEMCALNeutralWeight(pid.fEMCALNeutralWeight),
	64	fPHOSPhotonWeight(pid.fPHOSPhotonWeight),
	65	fPHOSPi0Weight(pid.fPHOSPi0Weight),
	66	fPHOSElectronWeight(pid.fPHOSElectronWeight),
	67	fPHOSChargeWeight(pid.fPHOSChargeWeight),
	68	fPHOSNeutralWeight(pid.fPHOSNeutralWeight),
	69	fPHOSWeightFormula(pid.fPHOSWeightFormula),
	70	fPHOSPhotonWeightFormula(pid.fPHOSPhotonWeightFormula),
	71	fPHOSPi0WeightFormula(pid.fPHOSPi0WeightFormula),
	72	fDispCut(pid.fDispCut),fTOFCut(pid.fTOFCut),
6639984f	73	fDebug(pid.fDebug)
1c5acb87	74	{
	75	// cpy ctor
	76
	77	}
	78
	79	//_________________________________________________________________________
	80	AliCaloPID & AliCaloPID::operator = (const AliCaloPID & pid)
	81	{
	82	// assignment operator
	83
	84	if(&pid == this) return *this;
	85
	86	fEMCALPhotonWeight = pid. fEMCALPhotonWeight ;
	87	fEMCALPi0Weight = pid.fEMCALPi0Weight ;
	88	fEMCALElectronWeight = pid.fEMCALElectronWeight;
	89	fEMCALChargeWeight = pid.fEMCALChargeWeight;
	90	fEMCALNeutralWeight = pid.fEMCALNeutralWeight;
	91
	92	fPHOSPhotonWeight = pid.fPHOSPhotonWeight ;
	93	fPHOSPi0Weight = pid.fPHOSPi0Weight ;
	94	fPHOSElectronWeight = pid.fPHOSElectronWeight;
	95	fPHOSChargeWeight = pid.fPHOSChargeWeight;
	96	fPHOSNeutralWeight = pid.fPHOSNeutralWeight;
	97
	98	fPHOSWeightFormula = pid.fPHOSWeightFormula;
	99	fPHOSPhotonWeightFormula = pid.fPHOSPhotonWeightFormula;
	100	fPHOSPi0WeightFormula = pid.fPHOSPi0WeightFormula;
	101
	102	fDispCut = pid.fDispCut;
	103	fTOFCut = pid.fTOFCut;
	104	fDebug = pid.fDebug;
6639984f	105
1c5acb87	106	return *this;
	107
	108	}
	109
	110	//_________________________________
	111	AliCaloPID::~AliCaloPID() {
	112	//Dtor
	113
	114	if(fPHOSPhotonWeightFormula) delete fPHOSPhotonWeightFormula ;
	115	if(fPHOSPi0WeightFormula) delete fPHOSPi0WeightFormula ;
	116
	117	}
	118
	119
1c5acb87	120	//_______________________________________________________________
	121	void AliCaloPID::InitParameters()
	122	{
	123	//Initialize the parameters of the PID.
	124
	125	fEMCALPhotonWeight = 0.8 ;
	126	fEMCALPi0Weight = 0.5 ;
	127	fEMCALElectronWeight = 0.8 ;
	128	fEMCALChargeWeight = 0.5 ;
	129	fEMCALNeutralWeight = 0.5 ;
	130
	131	fPHOSPhotonWeight = 0.75 ;
	132	fPHOSPi0Weight = 0.8 ;
	133	fPHOSElectronWeight = 0.5 ;
	134	fPHOSChargeWeight = 0.5 ;
	135	fPHOSNeutralWeight = 0.5 ;
	136
	137	//Formula to set the PID weight threshold for photon or pi0
	138	fPHOSWeightFormula = kTRUE;
	139	fPHOSPhotonWeightFormula =
	140	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))");
	141	fPHOSPi0WeightFormula =
	142	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))");
	143
	144	fDispCut = 1.5;
	145	fTOFCut = 5.e-9;
	146	fDebug = -1;
	147	}
	148
	149	//_______________________________________________________________
	150	Int_t AliCaloPID::GetPdg(const TString calo, const Double_t * pid, const Float_t energy) const {
	151	//Return most probable identity of the particle.
	152
	153	if(!pid) AliFatal("pid pointer not initialized!!!");
	154
	155	Float_t wPh = fPHOSPhotonWeight ;
	156	Float_t wPi0 = fPHOSPi0Weight ;
	157	Float_t wE = fPHOSElectronWeight ;
	158	Float_t wCh = fPHOSChargeWeight ;
	159	Float_t wNe = fPHOSNeutralWeight ;
	160
	161
	162	if(calo == "PHOS" && fPHOSWeightFormula){
	163	wPh = fPHOSPhotonWeightFormula->Eval(energy) ;
	164	wPi0 = fPHOSPi0WeightFormula->Eval(energy);
	165	}
	166
	167	if(calo == "EMCAL"){
	168
	169	wPh = fEMCALPhotonWeight ;
	170	wPi0 = fEMCALPi0Weight ;
	171	wE = fEMCALElectronWeight ;
	172	wCh = fEMCALChargeWeight ;
	173	wNe = fEMCALNeutralWeight ;
	174
	175	}
	176
	177	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",
	178	calo.Data(),pid[AliAODCluster::kPhoton], pid[AliAODCluster::kPi0],
	179	pid[AliAODCluster::kElectron], pid[AliAODCluster::kEleCon],
	180	pid[AliAODCluster::kPion], pid[AliAODCluster::kKaon], pid[AliAODCluster::kProton],
	181	pid[AliAODCluster::kNeutron], pid[AliAODCluster::kKaon0]);
	182
	183	Int_t pdg = kNeutralUnknown ;
184	Float_t chargedHadronWeight = pid[AliAODCluster::kProton]+pid[AliAODCluster::kKaon]+
185	pid[AliAODCluster::kPion]+pid[AliAODCluster::kMuon];
186	Float_t neutralHadronWeight = pid[AliAODCluster::kNeutron]+pid[AliAODCluster::kKaon0];
187	Float_t allChargedWeight = pid[AliAODCluster::kElectron]+pid[AliAODCluster::kEleCon]+ chargedHadronWeight;
188	Float_t allNeutralWeight = pid[AliAODCluster::kPhoton]+pid[AliAODCluster::kPi0]+ neutralHadronWeight;
189
190	//Select most probable ID
191	if(calo=="PHOS"){
192	if(pid[AliAODCluster::kPhoton] > wPh) pdg = kPhoton ;
193	else if(pid[AliAODCluster::kPi0] > wPi0) pdg = kPi0 ;
194	else if(pid[AliAODCluster::kElectron] > wE) pdg = kElectron ;
195	else if(pid[AliAODCluster::kEleCon] > wE) pdg = kEleCon ;
196	else if(chargedHadronWeight > wCh) pdg = kChargedHadron ;
197	else if(neutralHadronWeight > wNe) pdg = kNeutralHadron ;
198	else if(allChargedWeight > allNeutralWeight)
199	pdg = kChargedUnknown ;
200	else
201	pdg = kNeutralUnknown ;
202	}
203	else{//EMCAL
7175a03a	204	if(pid[AliAODCluster::kPhoton] > wPh) pdg = kPhoton ;
1c5acb87	205	else if(pid[AliAODCluster::kPi0] > wPi0) pdg = kPi0 ;
7175a03a	206	else if(pid[AliAODCluster::kElectron] > wE) pdg = kElectron ;
1c5acb87	207	else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;
	208	else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ;
	209	else pdg = kNeutralUnknown ;
	210
	211	}
	212
	213
	214	if(fDebug > 0)printf("AliCaloPID::GetPdg:Final Pdg: %d \n", pdg);
	215
	216
	217
	218	return pdg ;
	219
	220	}
	221
	222	//_______________________________________________________________
	223	Int_t AliCaloPID::GetPdg(const TString calo,const TLorentzVector mom, const AliAODCaloCluster * cluster) const {
	224	//Recalculated PID with all parameters
	225	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",
	226	calo.Data(),mom.E(),cluster->GetM02(),cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(),
	227	cluster->GetEmcCpvDistance(), cluster->GetDistToBadChannel(),cluster->GetNExMax());
	228
	229	if(calo == "EMCAL") {
	230	if(cluster->GetM02()< 0.25) return kPhoton ;
	231	else return kNeutralHadron ;
	232	}
	233
	234	// if(calo == "PHOS") {
	235	// if(cluster->GetM02()< 0.25) return kPhoton ;
	236	// else return kNeutralHadron ;
	237	// }
	238
	239	return kNeutralHadron ;
	240
	241	}
	242
	243	//__________________________________________________
	244	TString AliCaloPID::GetPIDParametersList() {
	245	//Put data member values in string to keep in output container
	246
	247	TString parList ; //this will be list of parameters used for this analysis.
	248	char onePar[255] ;
	249	sprintf(onePar,"--- AliCaloPID ---\n") ;
	250	parList+=onePar ;
	251	sprintf(onePar,"fDispCut =%2.2f (Cut on dispersion, used in PID evaluation) \n",fDispCut) ;
	252	parList+=onePar ;
	253	sprintf(onePar,"fTOFCut =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
	254	parList+=onePar ;
	255	sprintf(onePar,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)\n",fEMCALPhotonWeight) ;
	256	parList+=onePar ;
	257	sprintf(onePar,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)\n",fEMCALPi0Weight) ;
	258	parList+=onePar ;
	259	sprintf(onePar,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)\n",fEMCALElectronWeight) ;
	260	parList+=onePar ;
	261	sprintf(onePar,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)\n",fEMCALChargeWeight) ;
	262	parList+=onePar ;
	263	sprintf(onePar,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)\n",fEMCALNeutralWeight) ;
	264	parList+=onePar ;
	265	sprintf(onePar,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)\n",fPHOSPhotonWeight) ;
	266	parList+=onePar ;
	267	sprintf(onePar,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)\n",fPHOSPi0Weight) ;
	268	parList+=onePar ;
	269	sprintf(onePar,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)\n",fPHOSElectronWeight) ;
	270	parList+=onePar ;
271	sprintf(onePar,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)\n",fPHOSChargeWeight) ;
272	parList+=onePar ;
273	sprintf(onePar,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)\n",fPHOSNeutralWeight) ;
274	parList+=onePar ;
275
276	if(fPHOSWeightFormula){
277	parList+="PHOS Photon Weight Formula: "+(fPHOSPhotonWeightFormula->GetExpFormula("p"));
278	parList+="PHOS Pi0 Weight Formula: "+(fPHOSPi0WeightFormula->GetExpFormula("p"));
279	}
280
281	return parList;
282
283	}
284
285	//________________________________________________________________
286	void AliCaloPID::Print(const Option_t * opt) const
287	{
288
289	//Print some relevant parameters set for the analysis
290	if(! opt)
291	return;
292
293	printf("*** Print: %s %s ****\n", GetName(), GetTitle() ) ;
294
295	printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",
296	fPHOSPhotonWeight, fPHOSPi0Weight,
297	fPHOSElectronWeight, fPHOSChargeWeight, fPHOSNeutralWeight) ;
298	printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",
299	fEMCALPhotonWeight, fEMCALPi0Weight,
300	fEMCALElectronWeight, fEMCALChargeWeight, fEMCALNeutralWeight) ;
301
302	printf("PHOS Parametrized weight on? = %d\n", fPHOSWeightFormula) ;
303	if(fPHOSWeightFormula){
304	printf("Photon weight formula = %s\n", (fPHOSPhotonWeightFormula->GetExpFormula("p")).Data());
305	printf("Pi0 weight formula = %s\n", (fPHOSPhotonWeightFormula->GetExpFormula("p")).Data());
306	}
307
308	printf("TOF cut = %e\n",fTOFCut);
309	printf("Dispersion cut = %2.2f\n",fDispCut);
310	printf("Debug level = %d\n",fDebug);
7175a03a	311
1c5acb87	312	printf(" \n");
	313
	314	}
	315
	316	//_______________________________________________________________
	317	void AliCaloPID::SetPIDBits(const TString calo, const AliAODCaloCluster * cluster, AliAODPWG4Particle * ph) {
	318	//Set Bits for PID selection
	319
	320	//Dispersion/lambdas
	321	Double_t disp=cluster->GetDispersion() ;
	322	// Double_t m20=calo->GetM20() ;
	323	// Double_t m02=calo->GetM02() ;
	324	ph->SetDispBit(disp<fDispCut) ;
	325
	326	//TOF
	327	Double_t tof=cluster->GetTOF() ;
	328	ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
	329
	330	//Charged veto
	331	// Double_t cpvR=calo->GetEmcCpvDistance() ;
	332	Int_t ntr=cluster->GetNTracksMatched(); //number of track matched
	333	ph->SetChargedBit(ntr>0) ; //Temporary cut, should we evaluate distance?
	334
	335	//Set PID pdg
6639984f	336	ph->SetPdg(GetPdg(calo,cluster->PID(),ph->E()));
	337
	338	if(fDebug > 0){
	339	printf("AliCaloPID::SetPIDBits: TOF %e, Dispersion %2.2f, NTracks %d\n",tof , disp, ntr);
	340	printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
	341	ph->GetPdg(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit());
	342	}
1c5acb87	343	}
	344
	345