[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 $ */

/* History of cvs commits:
 *
 * $Log$
 *
 *
 */

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