Missing initialisation of headers added,
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliCaloPID.cxx
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 <TString.h>
28 #include <TFormula.h>
29
30 //---- ANALYSIS system ----
31 #include "AliCaloPID.h"
32 #include "AliAODCaloCluster.h"
33 #include "AliAODPWG4Particle.h"
34
35 ClassImp(AliCaloPID)
36
37
38 //________________________________________________
39 AliCaloPID::AliCaloPID() : 
40 TObject(), fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),  
41 fEMCALElectronWeight(0.),  fEMCALChargeWeight(0.),
42 fEMCALNeutralWeight(0.),
43 fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),  
44 fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) , 
45 fPHOSNeutralWeight(0.), fPHOSWeightFormula(0), 
46 fPHOSPhotonWeightFormula(0x0), fPHOSPi0WeightFormula(0x0),
47 fDispCut(0.),fTOFCut(0.), fDebug(-1)
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 fDispCut(pid.fDispCut),fTOFCut(pid.fTOFCut),
71 fDebug(pid.fDebug)
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   fDispCut  = pid.fDispCut;
101   fTOFCut   = pid.fTOFCut;
102   fDebug    = pid.fDebug;
103   
104   return *this;
105   
106 }
107
108 //_________________________________
109 AliCaloPID::~AliCaloPID() {
110   //Dtor
111   
112   if(fPHOSPhotonWeightFormula) delete  fPHOSPhotonWeightFormula ;
113   if(fPHOSPi0WeightFormula) delete  fPHOSPi0WeightFormula ;
114   
115 }
116
117
118 //_______________________________________________________________
119 void AliCaloPID::InitParameters()
120 {
121   //Initialize the parameters of the PID.
122   
123   fEMCALPhotonWeight   = 0.8 ;
124   fEMCALPi0Weight      = 0.5 ;
125   fEMCALElectronWeight = 0.8 ;
126   fEMCALChargeWeight   = 0.5 ;
127   fEMCALNeutralWeight  = 0.5 ;
128   
129   fPHOSPhotonWeight    = 0.75 ;
130   fPHOSPi0Weight       = 0.8 ;
131   fPHOSElectronWeight  = 0.5 ;
132   fPHOSChargeWeight    = 0.5 ;
133   fPHOSNeutralWeight   = 0.5 ;
134   
135   //Formula to set the PID weight threshold for photon or pi0
136   fPHOSWeightFormula = kTRUE;
137   fPHOSPhotonWeightFormula = 
138     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))");
139   fPHOSPi0WeightFormula = 
140     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))");
141   
142   fDispCut  = 1.5;
143   fTOFCut   = 5.e-9;
144   fDebug = -1;
145 }
146
147 //_______________________________________________________________
148 Int_t AliCaloPID::GetPdg(const TString calo, const Double_t * pid, const Float_t energy) const {
149   //Return most probable identity of the particle.
150   
151   if(!pid){ 
152     printf("AliCaloPID::GetPdg() - pid pointer not initialized!!!\n");
153     abort();
154   }
155   
156   Float_t wPh  =  fPHOSPhotonWeight ;
157   Float_t wPi0 =  fPHOSPi0Weight ;
158   Float_t wE   =  fPHOSElectronWeight ;
159   Float_t wCh  =  fPHOSChargeWeight ;
160   Float_t wNe  =  fPHOSNeutralWeight ;
161   
162   
163   if(calo == "PHOS" && fPHOSWeightFormula){
164     wPh  = fPHOSPhotonWeightFormula->Eval(energy) ;
165     wPi0 = fPHOSPi0WeightFormula->Eval(energy);
166   }
167   
168   if(calo == "EMCAL"){
169     
170     wPh  =  fEMCALPhotonWeight ;
171     wPi0 =  fEMCALPi0Weight ;
172     wE   =  fEMCALElectronWeight ;
173     wCh  =  fEMCALChargeWeight ;
174     wNe  =  fEMCALNeutralWeight ;
175     
176   }
177   
178   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",
179                          calo.Data(),pid[AliAODCluster::kPhoton], pid[AliAODCluster::kPi0],
180                          pid[AliAODCluster::kElectron], pid[AliAODCluster::kEleCon],
181                          pid[AliAODCluster::kPion], pid[AliAODCluster::kKaon], pid[AliAODCluster::kProton],
182                          pid[AliAODCluster::kNeutron], pid[AliAODCluster::kKaon0]);
183   
184   Int_t pdg = kNeutralUnknown ;
185   Float_t chargedHadronWeight = pid[AliAODCluster::kProton]+pid[AliAODCluster::kKaon]+
186     pid[AliAODCluster::kPion]+pid[AliAODCluster::kMuon];
187   Float_t neutralHadronWeight = pid[AliAODCluster::kNeutron]+pid[AliAODCluster::kKaon0];
188   Float_t allChargedWeight    = pid[AliAODCluster::kElectron]+pid[AliAODCluster::kEleCon]+ chargedHadronWeight;
189   Float_t allNeutralWeight    = pid[AliAODCluster::kPhoton]+pid[AliAODCluster::kPi0]+ neutralHadronWeight;
190   
191   //Select most probable ID
192   if(calo=="PHOS"){
193     if(pid[AliAODCluster::kPhoton] > wPh) pdg = kPhoton ;
194     else if(pid[AliAODCluster::kPi0] > wPi0) pdg = kPi0 ; 
195     else if(pid[AliAODCluster::kElectron] > wE)  pdg = kElectron ;
196     else if(pid[AliAODCluster::kEleCon] >  wE) pdg = kEleCon ;
197     else if(chargedHadronWeight > wCh) pdg = kChargedHadron ;  
198     else if(neutralHadronWeight > wNe) pdg = kNeutralHadron ; 
199     else if(allChargedWeight >  allNeutralWeight)
200       pdg = kChargedUnknown ; 
201     else 
202       pdg = kNeutralUnknown ;
203   }
204   else{//EMCAL
205     if(pid[AliAODCluster::kPhoton]+pid[AliAODCluster::kElectron]  > wPh) pdg = kPhoton ; //temporal sollution until track matching for electrons is considered
206     //if(pid[AliAODCluster::kPhoton]  > wPh) pdg = kPhoton ;
207     else if(pid[AliAODCluster::kPi0] > wPi0) pdg = kPi0 ; 
208     //else if(pid[AliAODCluster::kElectron]  > wE) pdg = kElectron ;
209     else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;  
210     else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ; 
211     else pdg =  kNeutralUnknown ;
212     
213   }
214   
215   if(fDebug > 0)printf("AliCaloPID::GetPdg:Final Pdg: %d \n", pdg);
216    
217   return pdg ;
218   
219 }
220
221 //_______________________________________________________________
222 Int_t AliCaloPID::GetPdg(const TString calo,const TLorentzVector mom, const AliAODCaloCluster * cluster) const {
223   //Recalculated PID with all parameters
224   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",
225                        calo.Data(),mom.E(),cluster->GetM02(),cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(), 
226                        cluster->GetEmcCpvDistance(), cluster->GetDistToBadChannel(),cluster->GetNExMax());
227   
228   if(calo == "EMCAL") {
229     if(cluster->GetM02() < 0.25) return kPhoton ;
230     //else return  kNeutralHadron ; 
231         else return  kPi0 ;
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         sprintf(onePar,"PHOS Photon Weight Formula: %s\n",(fPHOSPhotonWeightFormula->GetExpFormula("p")).Data()) ;
278     parList+=onePar;
279         sprintf(onePar,"PHOS Pi0    Weight Formula: %s\n",(fPHOSPi0WeightFormula->GetExpFormula("p")).Data()) ;
280         parList+=onePar;          
281   }
282   
283   return parList; 
284   
285 }
286
287 //________________________________________________________________
288 void AliCaloPID::Print(const Option_t * opt) const
289 {
290   
291   //Print some relevant parameters set for the analysis
292   if(! opt)
293     return;
294   
295   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
296   
297   printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",  
298          fPHOSPhotonWeight,  fPHOSPi0Weight, 
299          fPHOSElectronWeight,  fPHOSChargeWeight,   fPHOSNeutralWeight) ; 
300   printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",   
301          fEMCALPhotonWeight,  fEMCALPi0Weight, 
302          fEMCALElectronWeight,  fEMCALChargeWeight,  fEMCALNeutralWeight) ; 
303   
304   printf("PHOS Parametrized weight on?  =     %d\n",  fPHOSWeightFormula) ; 
305   if(fPHOSWeightFormula){
306     printf("Photon weight formula = %s\n", (fPHOSPhotonWeightFormula->GetExpFormula("p")).Data());
307     printf("Pi0    weight formula = %s\n", (fPHOSPhotonWeightFormula->GetExpFormula("p")).Data());
308   }
309   
310   printf("TOF cut        = %e\n",fTOFCut);
311   printf("Dispersion cut = %2.2f\n",fDispCut);
312   printf("Debug level    = %d\n",fDebug);
313   
314   printf(" \n");
315   
316
317
318 //_______________________________________________________________
319 void AliCaloPID::SetPIDBits(const TString calo, const AliAODCaloCluster * cluster, AliAODPWG4Particle * ph) {
320   //Set Bits for PID selection
321   
322   //Dispersion/lambdas
323   Double_t disp=cluster->GetDispersion()  ;
324         //    Double_t m20=calo->GetM20() ;
325         //    Double_t m02=calo->GetM02() ; 
326   ph->SetDispBit(disp<fDispCut) ;  
327   
328   //TOF
329   Double_t tof=cluster->GetTOF()  ;
330   ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ; 
331   
332   //Charged veto
333   //    Double_t cpvR=calo->GetEmcCpvDistance() ; 
334   Int_t ntr=cluster->GetNTracksMatched();  //number of track matched
335   ph->SetChargedBit(ntr>0) ;  //Temporary cut, should we evaluate distance?
336   
337   //Set PID pdg
338   ph->SetPdg(GetPdg(calo,cluster->PID(),ph->E()));
339   
340   if(fDebug > 0){ 
341     printf("AliCaloPID::SetPIDBits: TOF %e, Dispersion %2.2f, NTracks %d\n",tof , disp, ntr);   
342     printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
343              ph->GetPdg(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit()); 
344   }
345 }
346
347