[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"
#include "AliStack.h"
#include "TParticle.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), fMCGenerator("")
{
	//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),fMCGenerator(pid.fMCGenerator)
{
	// 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;
	fMCGenerator = pid.fMCGenerator;

	return *this;
	
}

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


//_________________________________________________________________________
Int_t AliCaloPID::CheckOrigin(const Int_t label, AliStack * stack) const {
  //Play with the MC stack if available
  //Check origin of the candidates, good for PYTHIA
  
  if(!stack) AliFatal("Stack is not available, check analysis settings in configuration file, STOP!!");
  
  if(label >= 0 && label <  stack->GetNtrack()){
    //Mother
    TParticle * mom = stack->Particle(label);
    Int_t mPdg = TMath::Abs(mom->GetPdgCode());
    Int_t mStatus =  mom->GetStatusCode() ;
    Int_t iParent =  mom->GetFirstMother() ;
    if(fDebug > 0 && label < 8 ) printf("AliCaloPID::CheckOrigin: Mother is parton %d\n",iParent);
    
    //GrandParent
    TParticle * parent = new TParticle ;
    Int_t pPdg = -1;
    Int_t pStatus =-1;
    if(iParent > 0){
      parent = stack->Particle(iParent);
      pPdg = TMath::Abs(parent->GetPdgCode());
      pStatus = parent->GetStatusCode();  
    }
    else if(fDebug > 0 ) printf("AliCaloPID::CheckOrigin: Parent with label %d\n",iParent);
    
    //return tag
    if(mPdg == 22){
      if(mStatus == 1){
	if(fMCGenerator == "PYTHIA"){
	  if(iParent < 8 && iParent > 5) {//outgoing partons
	    if(pPdg == 22) return kMCPrompt;
	    else  return kMCFragmentation;
	  }//Outgoing partons
	  else if(pStatus == 11){//Decay
	    if(pPdg == 111) return kMCPi0Decay ;
	    else if (pPdg == 321)  return kMCEtaDecay ;
	    else  return kMCOtherDecay ;
	  }//Decay
	  else return kMCISR; //Initial state radiation
	}//PYTHIA

	else if(fMCGenerator == "HERWIG"){	  
	  if(pStatus < 197){//Not decay
 	    while(1){
	      if(parent->GetFirstMother()<=5) break;
	      iParent = parent->GetFirstMother();
	      parent=stack->Particle(iParent);
	      pStatus= parent->GetStatusCode();
	      pPdg = parent->GetPdgCode();
	    }//Look for the parton
	    
	    if(iParent < 8 && iParent > 5) {
	      if(pPdg == 22) return kMCPrompt;
	      else  return kMCFragmentation;
	    }
	    return kMCISR;//Initial state radiation
	  }//Not decay
	  else{//Decay
	    if(pPdg == 111) return kMCPi0Decay ;
	    else if (pPdg == 321)  return kMCEtaDecay ;
	    else  return kMCOtherDecay ;
	  }//Decay
	}//HERWIG
	else return  kMCUnknown;
      }//Status 1 : Pythia generated
      else if(mStatus == 0){
	if(pPdg ==22 || pPdg ==11) return kMCConversion ;
	if(pPdg == 111) return kMCPi0Decay ;
	else if (pPdg == 221)  return kMCEtaDecay ;
	else  return kMCOtherDecay ;
      }//status 0 : geant generated
    }//Mother Photon
    else if(mPdg == 111)  return kMCPi0 ;
    else if(mPdg == 221)  return kMCEta ;
    else if(mPdg ==11){
      printf("Origin electron, pT %f\n",mom->Pt());

      if(mStatus == 0) return kMCConversion ;
      else return kMCElectron ;
    }
    else return kMCUnknown;
  }//Good label value
  else{
    if(label < 0 ) printf("AliCaloPID::CheckOrigin: *** bad label or no stack ***:  label %d \n", label);
    if(label >=  stack->GetNtrack()) printf("AliCaloPID::CheckOrigin: *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
    return kMCUnknown;
  }//Bad label
	
  return kMCUnknown;
  
}

//_______________________________________________________________
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;
	fMCGenerator = "PYTHIA";
}

//_______________________________________________________________
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("MC Generator   = %s\n",fMCGenerator.Data());

	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"
	36	#include "AliStack.h"
	37	#include "TParticle.h"
	38
	39	ClassImp(AliCaloPID)
	40
	41
	42	//________________________________________________
	43	AliCaloPID::AliCaloPID() :
	44	TObject(), fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
	45	fEMCALElectronWeight(0.), fEMCALChargeWeight(0.),
	46	fEMCALNeutralWeight(0.),
	47	fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
	48	fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) ,
	49	fPHOSNeutralWeight(0.), fPHOSWeightFormula(0),
	50	fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0),
7175a03a	51	fDispCut(0.),fTOFCut(0.), fDebug(-1), fMCGenerator("")
1c5acb87	52	{
	53	//Ctor
	54
	55	//Initialize parameters
	56	InitParameters();
	57	}
	58
	59	//____________________________________________________________________________
	60	AliCaloPID::AliCaloPID(const AliCaloPID & pid) :
	61	TObject(pid), fEMCALPhotonWeight(pid.fEMCALPhotonWeight),
	62	fEMCALPi0Weight(pid.fEMCALPi0Weight),
	63	fEMCALElectronWeight(pid.fEMCALElectronWeight),
	64	fEMCALChargeWeight(pid.fEMCALChargeWeight),
	65	fEMCALNeutralWeight(pid.fEMCALNeutralWeight),
	66	fPHOSPhotonWeight(pid.fPHOSPhotonWeight),
	67	fPHOSPi0Weight(pid.fPHOSPi0Weight),
	68	fPHOSElectronWeight(pid.fPHOSElectronWeight),
	69	fPHOSChargeWeight(pid.fPHOSChargeWeight),
	70	fPHOSNeutralWeight(pid.fPHOSNeutralWeight),
	71	fPHOSWeightFormula(pid.fPHOSWeightFormula),
	72	fPHOSPhotonWeightFormula(pid.fPHOSPhotonWeightFormula),
	73	fPHOSPi0WeightFormula(pid.fPHOSPi0WeightFormula),
	74	fDispCut(pid.fDispCut),fTOFCut(pid.fTOFCut),
7175a03a	75	fDebug(pid.fDebug),fMCGenerator(pid.fMCGenerator)
1c5acb87	76	{
	77	// cpy ctor
	78
	79	}
	80
	81	//_________________________________________________________________________
	82	AliCaloPID & AliCaloPID::operator = (const AliCaloPID & pid)
	83	{
	84	// assignment operator
	85
	86	if(&pid == this) return *this;
	87
	88	fEMCALPhotonWeight = pid. fEMCALPhotonWeight ;
	89	fEMCALPi0Weight = pid.fEMCALPi0Weight ;
	90	fEMCALElectronWeight = pid.fEMCALElectronWeight;
	91	fEMCALChargeWeight = pid.fEMCALChargeWeight;
	92	fEMCALNeutralWeight = pid.fEMCALNeutralWeight;
	93
	94	fPHOSPhotonWeight = pid.fPHOSPhotonWeight ;
	95	fPHOSPi0Weight = pid.fPHOSPi0Weight ;
	96	fPHOSElectronWeight = pid.fPHOSElectronWeight;
	97	fPHOSChargeWeight = pid.fPHOSChargeWeight;
	98	fPHOSNeutralWeight = pid.fPHOSNeutralWeight;
	99
	100	fPHOSWeightFormula = pid.fPHOSWeightFormula;
	101	fPHOSPhotonWeightFormula = pid.fPHOSPhotonWeightFormula;
	102	fPHOSPi0WeightFormula = pid.fPHOSPi0WeightFormula;
	103
	104	fDispCut = pid.fDispCut;
	105	fTOFCut = pid.fTOFCut;
	106	fDebug = pid.fDebug;
7175a03a	107	fMCGenerator = pid.fMCGenerator;
7175a03a	108
1c5acb87	109	return *this;
	110
	111	}
	112
	113	//_________________________________
	114	AliCaloPID::~AliCaloPID() {
	115	//Dtor
	116
	117	if(fPHOSPhotonWeightFormula) delete fPHOSPhotonWeightFormula ;
	118	if(fPHOSPi0WeightFormula) delete fPHOSPi0WeightFormula ;
	119
	120	}
	121
	122
	123	//_________________________________________________________________________
	124	Int_t AliCaloPID::CheckOrigin(const Int_t label, AliStack * stack) const {
7175a03a	125	//Play with the MC stack if available
	126	//Check origin of the candidates, good for PYTHIA
	127
	128	if(!stack) AliFatal("Stack is not available, check analysis settings in configuration file, STOP!!");
	129
	130	if(label >= 0 && label < stack->GetNtrack()){
	131	//Mother
	132	TParticle * mom = stack->Particle(label);
	133	Int_t mPdg = TMath::Abs(mom->GetPdgCode());
	134	Int_t mStatus = mom->GetStatusCode() ;
	135	Int_t iParent = mom->GetFirstMother() ;
	136	if(fDebug > 0 && label < 8 ) printf("AliCaloPID::CheckOrigin: Mother is parton %d\n",iParent);
	137
	138	//GrandParent
	139	TParticle * parent = new TParticle ;
	140	Int_t pPdg = -1;
	141	Int_t pStatus =-1;
	142	if(iParent > 0){
	143	parent = stack->Particle(iParent);
	144	pPdg = TMath::Abs(parent->GetPdgCode());
	145	pStatus = parent->GetStatusCode();
	146	}
	147	else if(fDebug > 0 ) printf("AliCaloPID::CheckOrigin: Parent with label %d\n",iParent);
	148
	149	//return tag
	150	if(mPdg == 22){
	151	if(mStatus == 1){
	152	if(fMCGenerator == "PYTHIA"){
	153	if(iParent < 8 && iParent > 5) {//outgoing partons
	154	if(pPdg == 22) return kMCPrompt;
	155	else return kMCFragmentation;
	156	}//Outgoing partons
	157	else if(pStatus == 11){//Decay
	158	if(pPdg == 111) return kMCPi0Decay ;
	159	else if (pPdg == 321) return kMCEtaDecay ;
	160	else return kMCOtherDecay ;
	161	}//Decay
	162	else return kMCISR; //Initial state radiation
	163	}//PYTHIA
	164
	165	else if(fMCGenerator == "HERWIG"){
	166	if(pStatus < 197){//Not decay
	167	while(1){
	168	if(parent->GetFirstMother()<=5) break;
	169	iParent = parent->GetFirstMother();
	170	parent=stack->Particle(iParent);
	171	pStatus= parent->GetStatusCode();
	172	pPdg = parent->GetPdgCode();
	173	}//Look for the parton
	174
	175	if(iParent < 8 && iParent > 5) {
	176	if(pPdg == 22) return kMCPrompt;
	177	else return kMCFragmentation;
	178	}
	179	return kMCISR;//Initial state radiation
	180	}//Not decay
	181	else{//Decay
	182	if(pPdg == 111) return kMCPi0Decay ;
	183	else if (pPdg == 321) return kMCEtaDecay ;
	184	else return kMCOtherDecay ;
	185	}//Decay
	186	}//HERWIG
	187	else return kMCUnknown;
	188	}//Status 1 : Pythia generated
189	else if(mStatus == 0){
190	if(pPdg ==22 \|\| pPdg ==11) return kMCConversion ;
191	if(pPdg == 111) return kMCPi0Decay ;
192	else if (pPdg == 221) return kMCEtaDecay ;
193	else return kMCOtherDecay ;
194	}//status 0 : geant generated
195	}//Mother Photon
196	else if(mPdg == 111) return kMCPi0 ;
197	else if(mPdg == 221) return kMCEta ;
198	else if(mPdg ==11){
199	printf("Origin electron, pT %f\n",mom->Pt());
200
201	if(mStatus == 0) return kMCConversion ;
202	else return kMCElectron ;
203	}
204	else return kMCUnknown;
205	}//Good label value
206	else{
207	if(label < 0 ) printf("AliCaloPID::CheckOrigin: * bad label or no stack *: label %d \n", label);
208	if(label >= stack->GetNtrack()) printf("AliCaloPID::CheckOrigin: * large label *: label %d, n tracks %d \n", label, stack->GetNtrack());
209	return kMCUnknown;
210	}//Bad label
211
212	return kMCUnknown;
213
1c5acb87	214	}
	215
	216	//_______________________________________________________________
	217	void AliCaloPID::InitParameters()
	218	{
	219	//Initialize the parameters of the PID.
	220
	221	fEMCALPhotonWeight = 0.8 ;
	222	fEMCALPi0Weight = 0.5 ;
	223	fEMCALElectronWeight = 0.8 ;
	224	fEMCALChargeWeight = 0.5 ;
	225	fEMCALNeutralWeight = 0.5 ;
	226
	227	fPHOSPhotonWeight = 0.75 ;
	228	fPHOSPi0Weight = 0.8 ;
	229	fPHOSElectronWeight = 0.5 ;
	230	fPHOSChargeWeight = 0.5 ;
	231	fPHOSNeutralWeight = 0.5 ;
	232
	233	//Formula to set the PID weight threshold for photon or pi0
	234	fPHOSWeightFormula = kTRUE;
	235	fPHOSPhotonWeightFormula =
	236	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))");
	237	fPHOSPi0WeightFormula =
	238	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))");
	239
	240	fDispCut = 1.5;
	241	fTOFCut = 5.e-9;
	242	fDebug = -1;
7175a03a	243	fMCGenerator = "PYTHIA";
1c5acb87	244	}
	245
	246	//_______________________________________________________________
	247	Int_t AliCaloPID::GetPdg(const TString calo, const Double_t * pid, const Float_t energy) const {
	248	//Return most probable identity of the particle.
	249
	250	if(!pid) AliFatal("pid pointer not initialized!!!");
	251
	252	Float_t wPh = fPHOSPhotonWeight ;
	253	Float_t wPi0 = fPHOSPi0Weight ;
	254	Float_t wE = fPHOSElectronWeight ;
	255	Float_t wCh = fPHOSChargeWeight ;
	256	Float_t wNe = fPHOSNeutralWeight ;
	257
	258
	259	if(calo == "PHOS" && fPHOSWeightFormula){
	260	wPh = fPHOSPhotonWeightFormula->Eval(energy) ;
	261	wPi0 = fPHOSPi0WeightFormula->Eval(energy);
	262	}
	263
	264	if(calo == "EMCAL"){
	265
	266	wPh = fEMCALPhotonWeight ;
	267	wPi0 = fEMCALPi0Weight ;
	268	wE = fEMCALElectronWeight ;
	269	wCh = fEMCALChargeWeight ;
	270	wNe = fEMCALNeutralWeight ;
	271
	272	}
	273
	274	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",
	275	calo.Data(),pid[AliAODCluster::kPhoton], pid[AliAODCluster::kPi0],
	276	pid[AliAODCluster::kElectron], pid[AliAODCluster::kEleCon],
	277	pid[AliAODCluster::kPion], pid[AliAODCluster::kKaon], pid[AliAODCluster::kProton],
	278	pid[AliAODCluster::kNeutron], pid[AliAODCluster::kKaon0]);
	279
	280	Int_t pdg = kNeutralUnknown ;
	281	Float_t chargedHadronWeight = pid[AliAODCluster::kProton]+pid[AliAODCluster::kKaon]+
	282	pid[AliAODCluster::kPion]+pid[AliAODCluster::kMuon];
	283	Float_t neutralHadronWeight = pid[AliAODCluster::kNeutron]+pid[AliAODCluster::kKaon0];
	284	Float_t allChargedWeight = pid[AliAODCluster::kElectron]+pid[AliAODCluster::kEleCon]+ chargedHadronWeight;
	285	Float_t allNeutralWeight = pid[AliAODCluster::kPhoton]+pid[AliAODCluster::kPi0]+ neutralHadronWeight;
	286
	287	//Select most probable ID
	288	if(calo=="PHOS"){
	289	if(pid[AliAODCluster::kPhoton] > wPh) pdg = kPhoton ;
	290	else if(pid[AliAODCluster::kPi0] > wPi0) pdg = kPi0 ;
	291	else if(pid[AliAODCluster::kElectron] > wE) pdg = kElectron ;
	292	else if(pid[AliAODCluster::kEleCon] > wE) pdg = kEleCon ;
	293	else if(chargedHadronWeight > wCh) pdg = kChargedHadron ;
	294	else if(neutralHadronWeight > wNe) pdg = kNeutralHadron ;
	295	else if(allChargedWeight > allNeutralWeight)
	296	pdg = kChargedUnknown ;
	297	else
	298	pdg = kNeutralUnknown ;
	299	}
	300	else{//EMCAL
7175a03a	301	if(pid[AliAODCluster::kPhoton] > wPh) pdg = kPhoton ;
1c5acb87	302	else if(pid[AliAODCluster::kPi0] > wPi0) pdg = kPi0 ;
7175a03a	303	else if(pid[AliAODCluster::kElectron] > wE) pdg = kElectron ;
1c5acb87	304	else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;
	305	else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ;
	306	else pdg = kNeutralUnknown ;
	307
	308	}
	309
	310
	311	if(fDebug > 0)printf("AliCaloPID::GetPdg:Final Pdg: %d \n", pdg);
	312
	313
	314
	315	return pdg ;
	316
	317	}
	318
	319	//_______________________________________________________________
	320	Int_t AliCaloPID::GetPdg(const TString calo,const TLorentzVector mom, const AliAODCaloCluster * cluster) const {
	321	//Recalculated PID with all parameters
	322	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",
	323	calo.Data(),mom.E(),cluster->GetM02(),cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(),
	324	cluster->GetEmcCpvDistance(), cluster->GetDistToBadChannel(),cluster->GetNExMax());
	325
	326	if(calo == "EMCAL") {
	327	if(cluster->GetM02()< 0.25) return kPhoton ;
	328	else return kNeutralHadron ;
	329	}
	330
	331	// if(calo == "PHOS") {
	332	// if(cluster->GetM02()< 0.25) return kPhoton ;
	333	// else return kNeutralHadron ;
	334	// }
	335
	336	return kNeutralHadron ;
	337
	338	}
	339
	340	//__________________________________________________
	341	TString AliCaloPID::GetPIDParametersList() {
	342	//Put data member values in string to keep in output container
	343
	344	TString parList ; //this will be list of parameters used for this analysis.
	345	char onePar[255] ;
	346	sprintf(onePar,"--- AliCaloPID ---\n") ;
	347	parList+=onePar ;
	348	sprintf(onePar,"fDispCut =%2.2f (Cut on dispersion, used in PID evaluation) \n",fDispCut) ;
	349	parList+=onePar ;
	350	sprintf(onePar,"fTOFCut =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
	351	parList+=onePar ;
	352	sprintf(onePar,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)\n",fEMCALPhotonWeight) ;
	353	parList+=onePar ;
	354	sprintf(onePar,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)\n",fEMCALPi0Weight) ;
	355	parList+=onePar ;
	356	sprintf(onePar,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)\n",fEMCALElectronWeight) ;
	357	parList+=onePar ;
	358	sprintf(onePar,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)\n",fEMCALChargeWeight) ;
	359	parList+=onePar ;
	360	sprintf(onePar,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)\n",fEMCALNeutralWeight) ;
	361	parList+=onePar ;
	362	sprintf(onePar,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)\n",fPHOSPhotonWeight) ;
	363	parList+=onePar ;
	364	sprintf(onePar,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)\n",fPHOSPi0Weight) ;
	365	parList+=onePar ;
	366	sprintf(onePar,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)\n",fPHOSElectronWeight) ;
	367	parList+=onePar ;
368	sprintf(onePar,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)\n",fPHOSChargeWeight) ;
369	parList+=onePar ;
370	sprintf(onePar,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)\n",fPHOSNeutralWeight) ;
371	parList+=onePar ;
372
373	if(fPHOSWeightFormula){
374	parList+="PHOS Photon Weight Formula: "+(fPHOSPhotonWeightFormula->GetExpFormula("p"));
375	parList+="PHOS Pi0 Weight Formula: "+(fPHOSPi0WeightFormula->GetExpFormula("p"));
376	}
377
378	return parList;
379
380	}
381
382	//________________________________________________________________
383	void AliCaloPID::Print(const Option_t * opt) const
384	{
385
386	//Print some relevant parameters set for the analysis
387	if(! opt)
388	return;
389
390	printf("*** Print: %s %s ****\n", GetName(), GetTitle() ) ;
391
392	printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",
393	fPHOSPhotonWeight, fPHOSPi0Weight,
394	fPHOSElectronWeight, fPHOSChargeWeight, fPHOSNeutralWeight) ;
395	printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",
396	fEMCALPhotonWeight, fEMCALPi0Weight,
397	fEMCALElectronWeight, fEMCALChargeWeight, fEMCALNeutralWeight) ;
398
399	printf("PHOS Parametrized weight on? = %d\n", fPHOSWeightFormula) ;
400	if(fPHOSWeightFormula){
401	printf("Photon weight formula = %s\n", (fPHOSPhotonWeightFormula->GetExpFormula("p")).Data());
402	printf("Pi0 weight formula = %s\n", (fPHOSPhotonWeightFormula->GetExpFormula("p")).Data());
403	}
404
405	printf("TOF cut = %e\n",fTOFCut);
406	printf("Dispersion cut = %2.2f\n",fDispCut);
407	printf("Debug level = %d\n",fDebug);
7175a03a	408	printf("MC Generator = %s\n",fMCGenerator.Data());
7175a03a	409
1c5acb87	410	printf(" \n");
	411
	412	}
	413
	414	//_______________________________________________________________
	415	void AliCaloPID::SetPIDBits(const TString calo, const AliAODCaloCluster * cluster, AliAODPWG4Particle * ph) {
	416	//Set Bits for PID selection
	417
	418	//Dispersion/lambdas
	419	Double_t disp=cluster->GetDispersion() ;
	420	// Double_t m20=calo->GetM20() ;
	421	// Double_t m02=calo->GetM02() ;
	422	ph->SetDispBit(disp<fDispCut) ;
	423
	424	//TOF
	425	Double_t tof=cluster->GetTOF() ;
	426	ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
	427
	428	//Charged veto
	429	// Double_t cpvR=calo->GetEmcCpvDistance() ;
	430	Int_t ntr=cluster->GetNTracksMatched(); //number of track matched
	431	ph->SetChargedBit(ntr>0) ; //Temporary cut, should we evaluate distance?
	432
	433	//Set PID pdg
	434	ph->SetPdg(GetPdg(calo,cluster->PID(),ph->E()));
	435
	436	if(fDebug > 0){
	437	printf("AliCaloPID::SetPIDBits: TOF %e, Dispersion %2.2f, NTracks %d\n",tof , disp, ntr);
	438	printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
	439	ph->GetPdg(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit());
	440	}
	441	}
	442
	443