[u/mrichter/AliRoot.git] / PWG4 / 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 "AliAODCluster.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) 
{
  //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) 


{
  // 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;

  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))");
}


//________________________________________________________________
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<<<<<<<<<<<<\n");
//     fPHOSPhotonWeightFormula->Print();
//     printf(">>>>>>>>>>> Pi0    weight formula<<<<<<<<<<<<\n");
//     fPHOSPhotonWeightFormula->Print();
//   }
  printf(" \n");

} 

//_______________________________________________________________
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 ;
    
  }
  
//   printf("PID: 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
    //Temporal solution, electrons and photons not differenciated
    if(pid[AliAODCluster::kPhoton] + pid[AliAODCluster::kElectron]  > wPh) pdg = kPhoton ;
    else if(pid[AliAODCluster::kPi0] > wPi0) pdg = kPi0 ; 
    else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;  
    else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ; 
    else pdg =  kNeutralUnknown ;

  }


  //printf("Final Pdg: %d \n", pdg);
  

  return pdg ;
  
}

//_______________________________________________________________
Int_t AliCaloPID::GetPdg(const TString calo, const TLorentzVector mom, const Double_t l0, 
			 const Double_t l1, const Double_t disp, const Double_t tof, 
			 const Double_t distCPV) const {
  //Recalculated PID with all parameters
  AliDebug(2,Form("Calorimeter %s, E %3.2f, l0 %3.2f, l1 %3.2f, disp %3.2f, tof %3.2f, distCPV %3.2f",
		  calo.Data(),mom.E(),l0,l1,disp,tof,distCPV));

  if(calo == "EMCAL") {
    if(l0 < 0.25) return kPhoton ;
    else return  kNeutralHadron ; 
  }
  
  return  kNeutralHadron ; 
  
}
Commit	Line	Data
d92b41ad	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
d92b41ad	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 "AliAODCluster.h"
	35
	36	ClassImp(AliCaloPID)
	37
	38
	39	//________________________________________________
	40	AliCaloPID::AliCaloPID() :
	41	TObject(), fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
	42	fEMCALElectronWeight(0.), fEMCALChargeWeight(0.),
	43	fEMCALNeutralWeight(0.),
	44	fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
	45	fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) ,
	46	fPHOSNeutralWeight(0.), fPHOSWeightFormula(0),
	47	fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0)
	48	{
	49	//Ctor
	50
	51	//Initialize parameters
	52	InitParameters();
	53	}
	54
	55	//____________________________________________________________________________
	56	AliCaloPID::AliCaloPID(const AliCaloPID & pid) :
	57	TObject(pid), fEMCALPhotonWeight(pid.fEMCALPhotonWeight),
	58	fEMCALPi0Weight(pid.fEMCALPi0Weight),
	59	fEMCALElectronWeight(pid.fEMCALElectronWeight),
	60	fEMCALChargeWeight(pid.fEMCALChargeWeight),
	61	fEMCALNeutralWeight(pid.fEMCALNeutralWeight),
	62	fPHOSPhotonWeight(pid.fPHOSPhotonWeight),
	63	fPHOSPi0Weight(pid.fPHOSPi0Weight),
	64	fPHOSElectronWeight(pid.fPHOSElectronWeight),
	65	fPHOSChargeWeight(pid.fPHOSChargeWeight),
	66	fPHOSNeutralWeight(pid.fPHOSNeutralWeight),
	67	fPHOSWeightFormula(pid.fPHOSWeightFormula),
	68	fPHOSPhotonWeightFormula(pid.fPHOSPhotonWeightFormula),
	69	fPHOSPi0WeightFormula(pid.fPHOSPi0WeightFormula)
	70
	71
	72	{
	73	// cpy ctor
	74
	75	}
	76
	77	//_________________________________________________________________________
	78	AliCaloPID & AliCaloPID::operator = (const AliCaloPID & pid)
	79	{
	80	// assignment operator
81
82	if(&pid == this) return *this;
83
84	fEMCALPhotonWeight = pid. fEMCALPhotonWeight ;
85	fEMCALPi0Weight = pid.fEMCALPi0Weight ;
86	fEMCALElectronWeight = pid.fEMCALElectronWeight;
87	fEMCALChargeWeight = pid.fEMCALChargeWeight;
88	fEMCALNeutralWeight = pid.fEMCALNeutralWeight;
89
90	fPHOSPhotonWeight = pid.fPHOSPhotonWeight ;
91	fPHOSPi0Weight = pid.fPHOSPi0Weight ;
92	fPHOSElectronWeight = pid.fPHOSElectronWeight;
93	fPHOSChargeWeight = pid.fPHOSChargeWeight;
94	fPHOSNeutralWeight = pid.fPHOSNeutralWeight;
95
96	fPHOSWeightFormula = pid.fPHOSWeightFormula;
97	fPHOSPhotonWeightFormula = pid.fPHOSPhotonWeightFormula;
98	fPHOSPi0WeightFormula = pid.fPHOSPi0WeightFormula;
99
100	return *this;
101
102	}
103
104	//_________________________________
105	AliCaloPID::~AliCaloPID() {
106	//Dtor
107
108	if(fPHOSPhotonWeightFormula) delete fPHOSPhotonWeightFormula ;
109	if(fPHOSPi0WeightFormula) delete fPHOSPi0WeightFormula ;
110
111	}
112
113
114	//_______________________________________________________________
115	void AliCaloPID::InitParameters()
116	{
117	//Initialize the parameters of the PID.
118
119	fEMCALPhotonWeight = 0.8 ;
120	fEMCALPi0Weight = 0.5 ;
121	fEMCALElectronWeight = 0.8 ;
122	fEMCALChargeWeight = 0.5 ;
123	fEMCALNeutralWeight = 0.5 ;
124
125	fPHOSPhotonWeight = 0.75 ;
126	fPHOSPi0Weight = 0.8 ;
127	fPHOSElectronWeight = 0.5 ;
128	fPHOSChargeWeight = 0.5 ;
129	fPHOSNeutralWeight = 0.5 ;
130
131	//Formula to set the PID weight threshold for photon or pi0
132	fPHOSWeightFormula = kTRUE;
133	fPHOSPhotonWeightFormula =
134	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))");
135	fPHOSPi0WeightFormula =
136	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))");
137	}
138
139
140	//________________________________________________________________
141	void AliCaloPID::Print(const Option_t * opt) const
142	{
143
144	//Print some relevant parameters set for the analysis
145	if(! opt)
146	return;
147
148	printf("*** Print: %s %s ****\n", GetName(), GetTitle() ) ;
149
150	printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",
151	fPHOSPhotonWeight, fPHOSPi0Weight,
152	fPHOSElectronWeight, fPHOSChargeWeight, fPHOSNeutralWeight) ;
153	printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",
154	fEMCALPhotonWeight, fEMCALPi0Weight,
155	fEMCALElectronWeight, fEMCALChargeWeight, fEMCALNeutralWeight) ;
156
157	printf("PHOS Parametrized weight on? = %d\n", fPHOSWeightFormula) ;
158	// if(fPHOSWeightFormula){
159	// printf(">>>>>>>>>>> Photon weight formula<<<<<<<<<<<<\n");
160	// fPHOSPhotonWeightFormula->Print();
161	// printf(">>>>>>>>>>> Pi0 weight formula<<<<<<<<<<<<\n");
162	// fPHOSPhotonWeightFormula->Print();
163	// }
164	printf(" \n");
165
166	}
167
168	//_______________________________________________________________
169	Int_t AliCaloPID::GetPdg(const TString calo, const Double_t * pid, const Float_t energy) const {
170	//Return most probable identity of the particle.
171
172	if(!pid) AliFatal("pid pointer not initialized!!!");
173
174	Float_t wPh = fPHOSPhotonWeight ;
175	Float_t wPi0 = fPHOSPi0Weight ;
176	Float_t wE = fPHOSElectronWeight ;
177	Float_t wCh = fPHOSChargeWeight ;
178	Float_t wNe = fPHOSNeutralWeight ;
179
180
181	if(calo == "PHOS" && fPHOSWeightFormula){
182	wPh = fPHOSPhotonWeightFormula->Eval(energy) ;
183	wPi0 = fPHOSPi0WeightFormula->Eval(energy);
184	}
185
186	if(calo == "EMCAL"){
187
188	wPh = fEMCALPhotonWeight ;
189	wPi0 = fEMCALPi0Weight ;
190	wE = fEMCALElectronWeight ;
191	wCh = fEMCALChargeWeight ;
192	wNe = fEMCALNeutralWeight ;
193
194	}
195
196	// printf("PID: 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",
197	// calo.Data(),pid[AliAODCluster::kPhoton], pid[AliAODCluster::kPi0],
198	// pid[AliAODCluster::kElectron], pid[AliAODCluster::kEleCon],
199	// pid[AliAODCluster::kPion], pid[AliAODCluster::kKaon], pid[AliAODCluster::kProton],
200	// pid[AliAODCluster::kNeutron], pid[AliAODCluster::kKaon0]);
201
202	Int_t pdg = kNeutralUnknown ;
203	Float_t chargedHadronWeight = pid[AliAODCluster::kProton]+pid[AliAODCluster::kKaon]+
204	pid[AliAODCluster::kPion]+pid[AliAODCluster::kMuon];
205	Float_t neutralHadronWeight = pid[AliAODCluster::kNeutron]+pid[AliAODCluster::kKaon0];
206	Float_t allChargedWeight = pid[AliAODCluster::kElectron]+pid[AliAODCluster::kEleCon]+ chargedHadronWeight;
207	Float_t allNeutralWeight = pid[AliAODCluster::kPhoton]+pid[AliAODCluster::kPi0]+ neutralHadronWeight;
208
209	//Select most probable ID
210	if(calo=="PHOS"){
211	if(pid[AliAODCluster::kPhoton] > wPh) pdg = kPhoton ;
212	else if(pid[AliAODCluster::kPi0] > wPi0) pdg = kPi0 ;
213	else if(pid[AliAODCluster::kElectron] > wE) pdg = kElectron ;
214	else if(pid[AliAODCluster::kEleCon] > wE) pdg = kEleCon ;
215	else if(chargedHadronWeight > wCh) pdg = kChargedHadron ;
216	else if(neutralHadronWeight > wNe) pdg = kNeutralHadron ;
217	else if(allChargedWeight > allNeutralWeight)
218	pdg = kChargedUnknown ;
219	else
220	pdg = kNeutralUnknown ;
221	}
222	else{//EMCAL
223	//Temporal solution, electrons and photons not differenciated
224	if(pid[AliAODCluster::kPhoton] + pid[AliAODCluster::kElectron] > wPh) pdg = kPhoton ;
225	else if(pid[AliAODCluster::kPi0] > wPi0) pdg = kPi0 ;
226	else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;
227	else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ;
228	else pdg = kNeutralUnknown ;
229
230	}
231
232
233	//printf("Final Pdg: %d \n", pdg);
234
235
236
237	return pdg ;
238
239	}
240
241	//_______________________________________________________________
242	Int_t AliCaloPID::GetPdg(const TString calo, const TLorentzVector mom, const Double_t l0,
243	const Double_t l1, const Double_t disp, const Double_t tof,
244	const Double_t distCPV) const {
245	//Recalculated PID with all parameters
2cc4014d	246	AliDebug(2,Form("Calorimeter %s, E %3.2f, l0 %3.2f, l1 %3.2f, disp %3.2f, tof %3.2f, distCPV %3.2f",
	247	calo.Data(),mom.E(),l0,l1,disp,tof,distCPV));
	248
d92b41ad	249	if(calo == "EMCAL") {
	250	if(l0 < 0.25) return kPhoton ;
	251	else return kNeutralHadron ;
	252	}
	253
	254	return kNeutralHadron ;
	255
	256	}
	257