21cbe653ca6ff637959328621ecd8e24efd8ebed
[u/mrichter/AliRoot.git] / PWG / CaloTrackCorrBase / 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
16 //_________________________________________________________________________
17 // Class for PID selection with calorimeters
18 // The Output of the main method GetIdentifiedParticleType is a PDG number identifying the cluster, 
19 // being kPhoton, kElectron, kPi0 ... as defined in the header file
20 //   - GetIdentifiedParticleType(const AliVCluster * cluster) 
21 //      Assignes a PID tag to the cluster, right now there is the possibility to : use bayesian weights from reco, 
22 //      recalculate them (EMCAL) or use other procedures not used in reco.
23 //      In order to recalculate Bayesian, it is necessary to load the EMCALUtils library
24 //      and do SwitchOnBayesianRecalculation().
25 //      To change the PID parameters from Low to High like the ones by default, use the constructor 
26 //      AliCaloPID(flux)
27 //      where flux is AliCaloPID::kLow or AliCaloPID::kHigh
28 //      If it is necessary to change the parameters use the constructor 
29 //      AliCaloPID(AliEMCALPIDUtils *utils) and set the parameters before.
30
31 //   - GetGetIdentifiedParticleTypeFromBayesian(const Double_t * pid, const Float_t energy)
32 //      Reads the PID weights array of the ESDs and depending on its magnitude identifies the particle, 
33 //      executed when bayesian is ON by GetIdentifiedParticleType(const AliVCluster * cluster) 
34 //   - SetPIDBits: Simple PID, depending on the thresholds fLOCut fTOFCut and even the
35 //     result of the PID bayesian a different PID bit is set. 
36 //
37 //   - IsTrackMatched(): Independent method that needs to be combined with GetIdentifiedParticleType to know if the cluster was matched
38 //
39 //*-- Author: Gustavo Conesa (INFN-LNF)
40 //////////////////////////////////////////////////////////////////////////////
41
42
43 // --- ROOT system ---
44 #include <TMath.h>
45 #include <TString.h>
46 #include <TList.h>
47
48 // ---- ANALYSIS system ----
49 #include "AliCaloPID.h"
50 #include "AliAODCaloCluster.h"
51 #include "AliVCaloCells.h"
52 #include "AliVTrack.h"
53 #include "AliAODPWG4Particle.h"
54 #include "AliCalorimeterUtils.h"
55 #include "AliVEvent.h"
56
57 // ---- Detector ----
58 #include "AliEMCALPIDUtils.h"
59
60 ClassImp(AliCaloPID)
61
62
63 //________________________
64 AliCaloPID::AliCaloPID() : 
65 TObject(),                fDebug(-1),                  fParticleFlux(kLow),
66 //Bayesian
67 fEMCALPIDUtils(),         fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
68 fEMCALPhotonWeight(0.),   fEMCALPi0Weight(0.),  
69 fEMCALElectronWeight(0.), fEMCALChargeWeight(0.),      fEMCALNeutralWeight(0.),
70 fPHOSPhotonWeight(0.),    fPHOSPi0Weight(0.),  
71 fPHOSElectronWeight(0.),  fPHOSChargeWeight(0.) ,      fPHOSNeutralWeight(0.), 
72 fPHOSWeightFormula(0),    fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
73 fPHOSPhotonWeightFormulaExpression(""), 
74 fPHOSPi0WeightFormulaExpression(""),
75 //PID calculation
76 fEMCALL0CutMax(100.),     fEMCALL0CutMin(0),           
77 fEMCALDEtaCut(2000.),     fEMCALDPhiCut(2000.),
78 fTOFCut(0.), 
79 fPHOSDispersionCut(1000), fPHOSRCut(1000),
80 fDoClusterSplitting(kFALSE),
81 fSplitM02MaxCut(0),       fSplitM02MinCut(0),          fSplitMinNCells(0),
82 fMassEtaMin(0),           fMassEtaMax(0),
83 fMassPi0Min(0),           fMassPi0Max(0),
84 fMassPhoMin(0),           fMassPhoMax(0)
85 {
86   //Ctor
87   
88   //Initialize parameters
89   InitParameters();
90 }
91
92 //________________________________________
93 AliCaloPID::AliCaloPID(const Int_t flux) : 
94 TObject(),                fDebug(-1),                  fParticleFlux(flux),
95 //Bayesian
96 fEMCALPIDUtils(),         fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
97 fEMCALPhotonWeight(0.),   fEMCALPi0Weight(0.),  
98 fEMCALElectronWeight(0.), fEMCALChargeWeight(0.),      fEMCALNeutralWeight(0.),
99 fPHOSPhotonWeight(0.),    fPHOSPi0Weight(0.),  
100 fPHOSElectronWeight(0.),  fPHOSChargeWeight(0.) ,      fPHOSNeutralWeight(0.), 
101 fPHOSWeightFormula(0),    fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
102 fPHOSPhotonWeightFormulaExpression(""), 
103 fPHOSPi0WeightFormulaExpression(""),
104 //PID calculation
105 fEMCALL0CutMax(100.),     fEMCALL0CutMin(0),           
106 fEMCALDEtaCut(2000.),     fEMCALDPhiCut(2000.),
107 fTOFCut(0.), 
108 fPHOSDispersionCut(1000), fPHOSRCut(1000),
109 fDoClusterSplitting(kFALSE),
110 fSplitM02MaxCut(0),       fSplitM02MinCut(0),          fSplitMinNCells(0),
111 fMassEtaMin(0),           fMassEtaMax(0),
112 fMassPi0Min(0),           fMassPi0Max(0),
113 fMassPhoMin(0),           fMassPhoMax(0)
114 {
115   //Ctor
116         
117   //Initialize parameters
118   InitParameters();
119   
120 }
121
122 //_______________________________________________
123 AliCaloPID::AliCaloPID(const TNamed * emcalpid) : 
124 TObject(),                   fDebug(-1),                  fParticleFlux(kLow),
125 //Bayesian
126 fEMCALPIDUtils((AliEMCALPIDUtils*)emcalpid),         
127 fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
128 fEMCALPhotonWeight(0.),      fEMCALPi0Weight(0.),  
129 fEMCALElectronWeight(0.),    fEMCALChargeWeight(0.),      fEMCALNeutralWeight(0.),
130 fPHOSPhotonWeight(0.),       fPHOSPi0Weight(0.),  
131 fPHOSElectronWeight(0.),     fPHOSChargeWeight(0.) ,      fPHOSNeutralWeight(0.), 
132 fPHOSWeightFormula(0),       fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
133 fPHOSPhotonWeightFormulaExpression(""), 
134 fPHOSPi0WeightFormulaExpression(""),
135 //PID calculation
136 fEMCALL0CutMax(100.),        fEMCALL0CutMin(0),   
137 fEMCALDEtaCut(2000.),        fEMCALDPhiCut(2000.),
138 fTOFCut(0.), 
139 fPHOSDispersionCut(1000),    fPHOSRCut(1000),
140 fDoClusterSplitting(kFALSE),
141 fSplitM02MaxCut(0),       fSplitM02MinCut(0),          fSplitMinNCells(0),
142 fMassEtaMin(0),           fMassEtaMax(0),
143 fMassPi0Min(0),           fMassPi0Max(0),
144 fMassPhoMin(0),           fMassPhoMax(0)
145 {
146   //Ctor
147   
148   //Initialize parameters
149   InitParameters();
150 }
151
152 //_______________________
153 AliCaloPID::~AliCaloPID() 
154 {
155   //Dtor
156   
157   delete fPHOSPhotonWeightFormula ;
158   delete fPHOSPi0WeightFormula ;
159   delete fEMCALPIDUtils ;
160   
161 }
162
163 //_______________________________
164 void AliCaloPID::InitParameters()
165 {
166   //Initialize the parameters of the PID.
167   
168   // Bayesian
169   fEMCALPhotonWeight   = 0.6 ;
170   fEMCALPi0Weight      = 0.6 ;
171   fEMCALElectronWeight = 0.6 ;
172   fEMCALChargeWeight   = 0.6 ;
173   fEMCALNeutralWeight  = 0.6 ;
174   
175   fPHOSPhotonWeight    = 0.6 ;
176   fPHOSPi0Weight       = 0.6 ;
177   fPHOSElectronWeight  = 0.6 ;
178   fPHOSChargeWeight    = 0.6 ;
179   fPHOSNeutralWeight   = 0.6 ;
180   
181   //Formula to set the PID weight threshold for photon or pi0
182   fPHOSWeightFormula       = kFALSE;
183   fPHOSPhotonWeightFormulaExpression = "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))";
184   fPHOSPi0WeightFormulaExpression    = "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))"   ;
185   
186   if(fRecalculateBayesian){
187     if(fParticleFlux == kLow){
188       printf("AliCaloPID::Init() - SetLOWFluxParam\n");
189       fEMCALPIDUtils->SetLowFluxParam() ;
190     }
191     else if (fParticleFlux == kHigh){
192       printf("AliCaloPID::Init() - SetHIGHFluxParam\n");
193       fEMCALPIDUtils->SetHighFluxParam() ;
194     }
195   }
196   
197   //PID recalculation, not bayesian
198   
199   //EMCAL
200   fEMCALL0CutMax = 0.3 ;
201   fEMCALL0CutMin = 0.01;
202   
203   fEMCALDPhiCut  = 0.05; // Same cut as in AliEMCALRecoUtils
204   fEMCALDEtaCut  = 0.025;// Same cut as in AliEMCALRecoUtils
205
206   // PHOS / EMCAL, not used
207   fTOFCut        = 1.e-6;
208   
209   //PHOS
210   fPHOSRCut          = 2. ; 
211   fPHOSDispersionCut = 2.5;
212   
213   // Cluster splitting
214   
215   fSplitM02MinCut = 0.26 ;
216   fSplitM02MaxCut = 100  ;
217   fSplitMinNCells = 4    ;
218
219   fMassEtaMin  = 0.4;
220   fMassEtaMax  = 0.6;
221   
222   fMassPi0Min  = 0.08;
223   fMassPi0Max  = 0.20;
224   
225   fMassPhoMin  = 0.0;
226   fMassPhoMax  = 0.05;
227   
228 }
229
230 //______________________________________________
231 AliEMCALPIDUtils *AliCaloPID::GetEMCALPIDUtils() 
232 {
233   // return pointer to AliEMCALPIDUtils, create it if needed
234   
235   if(!fEMCALPIDUtils) fEMCALPIDUtils = new AliEMCALPIDUtils ; 
236   return fEMCALPIDUtils ; 
237   
238 }
239
240
241 //______________________________________________________________________
242 Int_t AliCaloPID::GetIdentifiedParticleType(const AliVCluster * cluster) 
243 {
244   // Returns a PDG number corresponding to the likely ID of the cluster
245   
246   Float_t energy  = cluster->E();       
247   Float_t lambda0 = cluster->GetM02();
248   Float_t lambda1 = cluster->GetM20();
249   
250   // ---------------------
251   // Use bayesian approach
252   // ---------------------
253   
254   if(fUseBayesianWeights)
255   {
256     Double_t weights[AliPID::kSPECIESN];
257     
258     if(cluster->IsEMCAL() && fRecalculateBayesian)
259     {           
260       fEMCALPIDUtils->ComputePID(energy, lambda0);
261       for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = fEMCALPIDUtils->GetPIDFinal(i);
262     }
263     else 
264     {
265       for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = cluster->GetPID()[i];
266     }
267
268     if(fDebug > 0)  PrintClusterPIDWeights(weights);
269     
270     return GetIdentifiedParticleTypeFromBayesWeights(cluster->IsEMCAL(), weights, energy);
271   }
272   
273   // -------------------------------------------------------
274   // Calculate PID SS from data, do not use bayesian weights
275   // -------------------------------------------------------
276   
277   if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType: EMCAL %d?, E %3.2f, l0 %3.2f, l1 %3.2f, disp %3.2f, tof %1.11f, distCPV %3.2f, distToBC %1.1f, NMax %d\n",
278                         cluster->IsEMCAL(),energy,lambda0,cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(), 
279                         cluster->GetEmcCpvDistance(), cluster->GetDistanceToBadChannel(),cluster->GetNExMax());
280   
281   if(cluster->IsEMCAL())
282   {
283     if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType() - EMCAL SS %f <%f < %f?\n",fEMCALL0CutMin, lambda0, fEMCALL0CutMax);
284     
285     if(lambda0 < fEMCALL0CutMax && lambda0 > fEMCALL0CutMin) return kPhoton ;
286     else                                                     return kNeutralUnknown ; 
287   }    // EMCAL
288   else // PHOS
289   {    
290     if(TestPHOSDispersion(energy,lambda0,lambda1) < fPHOSDispersionCut) return kPhoton;
291     else                                                                return kNeutralUnknown;
292   }
293   
294 }
295
296 //_______________________________________________________________________________
297 Int_t AliCaloPID::GetIdentifiedParticleTypeFromBayesWeights(const Bool_t isEMCAL, 
298                                                             const Double_t * pid, 
299                                                             const Float_t energy) 
300 {
301   //Return most probable identity of the particle after bayesian weights calculated in reconstruction
302   
303   if(!pid)
304   { 
305     printf("AliCaloPID::GetIdentifiedParticleType() - pid pointer not initialized!!!\n");
306     abort();
307   }
308   
309   Float_t wPh  =  fPHOSPhotonWeight ;
310   Float_t wPi0 =  fPHOSPi0Weight ;
311   Float_t wE   =  fPHOSElectronWeight ;
312   Float_t wCh  =  fPHOSChargeWeight ;
313   Float_t wNe  =  fPHOSNeutralWeight ;
314   
315   if(!isEMCAL && fPHOSWeightFormula){
316     wPh  = GetPHOSPhotonWeightFormula()->Eval(energy) ;
317     wPi0 = GetPHOSPi0WeightFormula()   ->Eval(energy);
318   }
319   else
320   {
321     wPh  =  fEMCALPhotonWeight ;
322     wPi0 =  fEMCALPi0Weight ;
323     wE   =  fEMCALElectronWeight ;
324     wCh  =  fEMCALChargeWeight ;
325     wNe  =  fEMCALNeutralWeight ;
326   }
327   
328   if(fDebug > 0) PrintClusterPIDWeights(pid);
329     
330   Int_t pdg = kNeutralUnknown ;
331   Float_t chargedHadronWeight = pid[AliVCluster::kProton]+pid[AliVCluster::kKaon]+
332   pid[AliVCluster::kPion]+pid[AliVCluster::kMuon];
333   Float_t neutralHadronWeight = pid[AliVCluster::kNeutron]+pid[AliVCluster::kKaon0];
334   Float_t allChargedWeight    = pid[AliVCluster::kElectron]+pid[AliVCluster::kEleCon]+ chargedHadronWeight;
335   Float_t allNeutralWeight    = pid[AliVCluster::kPhoton]+pid[AliVCluster::kPi0]+ neutralHadronWeight;
336   
337   //Select most probable ID
338   if(!isEMCAL) // PHOS
339   {
340     if(pid[AliVCluster::kPhoton] > wPh)        pdg = kPhoton ;
341     else if(pid[AliVCluster::kPi0] > wPi0)     pdg = kPi0 ; 
342     else if(pid[AliVCluster::kElectron] > wE)  pdg = kElectron ;
343     else if(pid[AliVCluster::kEleCon] >  wE)   pdg = kEleCon ;
344     else if(chargedHadronWeight > wCh)         pdg = kChargedHadron ;  
345     else if(neutralHadronWeight > wNe)         pdg = kNeutralHadron ; 
346     else if(allChargedWeight >  allNeutralWeight)
347       pdg = kChargedUnknown ; 
348     else 
349       pdg = kNeutralUnknown ;
350   }
351   else //EMCAL
352   {
353     if(pid[AliVCluster::kPhoton]  > wPh)                     pdg = kPhoton ;
354     else if(pid[AliVCluster::kElectron]  > wE)               pdg = kElectron ;
355     else if(pid[AliVCluster::kPhoton]+pid[AliVCluster::kElectron]  > wPh) pdg = kPhoton ; //temporal sollution until track matching for electrons is considered
356     else if(pid[AliVCluster::kPi0] > wPi0)                   pdg = kPi0 ; 
357     else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;  
358     else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ; 
359     else                                                     pdg = kNeutralUnknown ;
360   }
361   
362   if(fDebug > 0)printf("AliCaloPID::GetIdentifiedParticleType:Final Pdg: %d, cluster energy %2.2f \n", pdg,energy);
363
364   return pdg ;
365   
366 }
367
368 //____________________________________________________________________________________________________
369 Int_t AliCaloPID::GetIdentifiedParticleTypeFromClusterSplitting(AliVCluster* cluster, 
370                                                                 AliVCaloCells* cells,
371                                                                 AliCalorimeterUtils * caloutils,
372                                                                 Double_t   vertex[3],
373                                                                 Int_t    & nMax,
374                                                                 Double_t & mass, Double_t & angle,
375                                                                 Double_t & e1  , Double_t & e2    ) 
376 {
377   // Split the cluster in 2, do invariant mass, get the mass and decide 
378   // if this is a photon, pi0, eta, ...
379   
380   Int_t   absId1 = -1; Int_t absId2 = -1;
381   Float_t l0     = cluster->GetM02();
382   const Int_t nc = cluster->GetNCells();
383   Int_t   absIdList[nc]; 
384   Float_t maxEList [nc]; 
385   
386   nMax  = -1 ;
387   mass  = -1.;
388   angle = -1.;
389   
390   //If too small or big E or low number of cells, or close to a bad channel skip it
391   if( l0 < fSplitM02MinCut || l0 > fSplitM02MaxCut || nc < fSplitMinNCells) 
392   {
393     return kNeutralUnknown ; 
394   }
395   
396   // Get Number of local maxima
397   nMax      = caloutils->GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ;  
398   
399   //---------------------------------------------------------------------
400   // Get the 2 max indeces and do inv mass
401   //---------------------------------------------------------------------
402   
403   if     ( nMax == 2 ) 
404   {
405     absId1 = absIdList[0];
406     absId2 = absIdList[1];
407   }
408   else if( nMax == 1 )
409   {
410     
411     absId1 = absIdList[0];
412     
413     //Find second highest energy cell
414     
415     TString  calorimeter = "EMCAL";
416     if(cluster->IsPHOS()) calorimeter = "PHOS";
417     Float_t enmax = 0 ;
418     for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++)
419     {
420       Int_t absId = cluster->GetCellsAbsId()[iDigit];
421       if( absId == absId1 ) continue ; 
422       Float_t endig = cells->GetCellAmplitude(absId);
423       caloutils->RecalibrateCellAmplitude(endig,calorimeter,absId); 
424       if(endig > enmax) 
425       {
426         enmax  = endig ;
427         absId2 = absId ;
428       }
429     }// cell loop
430   }// 1 maxima 
431   else
432   {  // n max > 2
433     // loop on maxima, find 2 highest
434     
435     // First max
436     Float_t enmax = 0 ;
437     for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++)
438     {
439       Float_t endig = maxEList[iDigit];
440       if(endig > enmax) 
441       {
442         enmax  = endig ;
443         absId1 = absIdList[iDigit];
444       }
445     }// first maxima loop
446     
447     // Second max 
448     Float_t enmax2 = 0;
449     for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++)
450     {
451       if(absIdList[iDigit]==absId1) continue;
452       Float_t endig = maxEList[iDigit];
453       if(endig > enmax2) 
454       {
455         enmax2  = endig ;
456         absId2 = absIdList[iDigit];
457       }
458     }// second maxima loop
459     
460   } // n local maxima > 2
461   
462   //---------------------------------------------------------------------
463   // Split the cluster energy in 2, around the highest 2 local maxima
464   //---------------------------------------------------------------------  
465   
466   AliAODCaloCluster cluster1(0, 0,NULL,0.,NULL,NULL,1,0);
467   AliAODCaloCluster cluster2(1, 0,NULL,0.,NULL,NULL,1,0);
468   
469   caloutils->SplitEnergy(absId1,absId2,cluster, cells, &cluster1, &cluster2,nMax); /*absIdList, maxEList,*/
470   
471   TLorentzVector cellMom1; 
472   TLorentzVector cellMom2;  
473   
474   cluster1.GetMomentum(cellMom1,vertex);
475   cluster2.GetMomentum(cellMom2,vertex);
476   
477   mass  = (cellMom1+cellMom2).M();
478   angle = cellMom2.Angle(cellMom1.Vect());
479   e1    = cluster1.E();
480   e2    = cluster2.E();
481   
482   if     (mass < fMassPhoMax && mass > fMassPhoMin) return kPhoton;
483   else if(mass < fMassPi0Max && mass > fMassPi0Min) return kPi0;
484   else if(mass < fMassEtaMax && mass > fMassEtaMin) return kEta;
485   else                                              return kNeutralUnknown;
486   
487 }
488
489 //_________________________________________
490 TString  AliCaloPID::GetPIDParametersList()  
491 {
492   //Put data member values in string to keep in output container
493   
494   TString parList ; //this will be list of parameters used for this analysis.
495   const Int_t buffersize = 255;
496   char onePar[buffersize] ;
497   snprintf(onePar,buffersize,"--- AliCaloPID ---\n") ;
498   parList+=onePar ;     
499   if(fUseBayesianWeights){
500     snprintf(onePar,buffersize,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)\n",fEMCALPhotonWeight) ;
501     parList+=onePar ;
502     snprintf(onePar,buffersize,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)\n",fEMCALPi0Weight) ;
503     parList+=onePar ;
504     snprintf(onePar,buffersize,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)\n",fEMCALElectronWeight) ;
505     parList+=onePar ;
506     snprintf(onePar,buffersize,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)\n",fEMCALChargeWeight) ;
507     parList+=onePar ;
508     snprintf(onePar,buffersize,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)\n",fEMCALNeutralWeight) ;
509     parList+=onePar ;
510     snprintf(onePar,buffersize,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)\n",fPHOSPhotonWeight) ;
511     parList+=onePar ;
512     snprintf(onePar,buffersize,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)\n",fPHOSPi0Weight) ;
513     parList+=onePar ;
514     snprintf(onePar,buffersize,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)\n",fPHOSElectronWeight) ;
515     parList+=onePar ;
516     snprintf(onePar,buffersize,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)\n",fPHOSChargeWeight) ;
517     parList+=onePar ;
518     snprintf(onePar,buffersize,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)\n",fPHOSNeutralWeight) ;
519     parList+=onePar ;
520     
521     if(fPHOSWeightFormula){
522       snprintf(onePar,buffersize,"PHOS Photon Weight Formula: %s\n",fPHOSPhotonWeightFormulaExpression.Data() ) ;
523       parList+=onePar;
524       snprintf(onePar,buffersize,"PHOS Pi0    Weight Formula: %s\n",fPHOSPi0WeightFormulaExpression.Data()    ) ;
525       parList+=onePar;    
526     }
527   }
528   else {
529     snprintf(onePar,buffersize,"EMCAL: fEMCALL0CutMin =%2.2f, fEMCALL0CutMax =%2.2f  (Cut on Shower Shape) \n",fEMCALL0CutMin, fEMCALL0CutMax) ;
530     parList+=onePar ;
531     snprintf(onePar,buffersize,"EMCAL: fEMCALDEtaCut =%2.2f, fEMCALDPhiCut =%2.2f  (Cut on track matching) \n",fEMCALDEtaCut, fEMCALDPhiCut) ;
532     parList+=onePar ;
533     snprintf(onePar,buffersize,"fTOFCut  =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
534     parList+=onePar ;   
535     snprintf(onePar,buffersize,"fPHOSRCut =%2.2f, fPHOSDispersionCut =%2.2f  (Cut on Shower Shape and CPV) \n",fPHOSRCut,fPHOSDispersionCut) ;
536     parList+=onePar ;
537     
538   }
539   
540   if(fDoClusterSplitting)
541   {
542     snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n",    fSplitM02MinCut, fSplitM02MaxCut) ;
543     parList+=onePar ;
544     snprintf(onePar,buffersize,"fMinNCells =%d \n",        fSplitMinNCells) ;
545     parList+=onePar ;    
546     snprintf(onePar,buffersize,"pi0 : %2.1f < m <%2.1f\n", fMassPi0Min,fMassPi0Max);
547     parList+=onePar ;
548     snprintf(onePar,buffersize,"eta : %2.1f < m <%2.1f\n", fMassEtaMin,fMassEtaMax);
549     parList+=onePar ;
550     snprintf(onePar,buffersize,"conv: %2.1f < m <%2.1f\n", fMassPhoMin,fMassPhoMax);
551     parList+=onePar ;
552   }
553   
554   return parList; 
555   
556 }
557
558 //________________________________________________
559 void AliCaloPID::Print(const Option_t * opt) const
560 {
561   
562   //Print some relevant parameters set for the analysis
563   if(! opt)
564     return;
565   
566   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
567   
568   if(fUseBayesianWeights)
569   {
570     printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",  
571            fPHOSPhotonWeight,   fPHOSPi0Weight, 
572            fPHOSElectronWeight, fPHOSChargeWeight, fPHOSNeutralWeight) ; 
573     printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",   
574            fEMCALPhotonWeight,   fEMCALPi0Weight, 
575            fEMCALElectronWeight, fEMCALChargeWeight, fEMCALNeutralWeight) ; 
576     
577     printf("PHOS Parametrized weight on?  =     %d\n",  fPHOSWeightFormula) ; 
578     if(fPHOSWeightFormula)
579     {
580       printf("Photon weight formula = %s\n", fPHOSPhotonWeightFormulaExpression.Data());
581       printf("Pi0    weight formula = %s\n", fPHOSPi0WeightFormulaExpression   .Data());
582     }
583     if(fRecalculateBayesian) printf(" Recalculate bayesian with Particle Flux?    = %d\n",fParticleFlux);
584   }
585   else 
586   {
587     printf("TOF cut        = %e\n",                                   fTOFCut);
588     printf("EMCAL Lambda0 cut min = %2.2f; max = %2.2f\n",            fEMCALL0CutMin,fEMCALL0CutMax);
589     printf("EMCAL cluster-track dEta < %2.3f; dPhi < %2.3f\n",        fEMCALDEtaCut, fEMCALDPhiCut);
590     printf("PHOS Treac matching cut =%2.2f, Dispersion Cut =%2.2f \n",fPHOSRCut,     fPHOSDispersionCut) ;
591     
592   }
593   
594   if(fDoClusterSplitting)
595   {
596     printf("Min. N Cells =%d \n",         fSplitMinNCells) ;
597     printf("%2.2f < lambda_0^2 <%2.2f \n",fSplitM02MinCut,fSplitM02MaxCut);
598     printf("pi0 : %2.2f<m<%2.2f \n",      fMassPi0Min,fMassPi0Max);
599     printf("eta : %2.2f<m<%2.2f \n",      fMassEtaMin,fMassEtaMax);
600     printf("phot: %2.2f<m<%2.2f \n",      fMassPhoMin,fMassPhoMax);
601   }
602   
603   printf(" \n");
604   
605
606
607 //_________________________________________________________________
608 void AliCaloPID::PrintClusterPIDWeights(const Double_t * pid) const
609 {
610   // print PID of cluster, (AliVCluster*)cluster->GetPID()
611   
612   printf("AliCaloPID::PrintClusterPIDWeights() \n \t ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, \n \t \
613          pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
614          pid[AliVCluster::kPhoton],    pid[AliVCluster::kPi0],
615          pid[AliVCluster::kElectron],  pid[AliVCluster::kEleCon],
616          pid[AliVCluster::kPion],      pid[AliVCluster::kKaon], 
617          pid[AliVCluster::kProton],
618          pid[AliVCluster::kNeutron],   pid[AliVCluster::kKaon0]);
619   
620 }
621
622 //___________________________________________________________________________
623 void AliCaloPID::SetPIDBits(AliVCluster * cluster, 
624                             AliAODPWG4Particle * ph, AliCalorimeterUtils* cu, 
625                             AliVEvent* event) 
626 {
627   //Set Bits for PID selection
628   
629   //Dispersion/lambdas
630   //Double_t disp= cluster->GetDispersion()  ;
631   Double_t l1  = cluster->GetM20() ;
632   Double_t l0  = cluster->GetM02() ; 
633   Bool_t isDispOK = kTRUE ;
634   if(cluster->IsPHOS()){ 
635     if(TestPHOSDispersion(ph->Pt(),l0,l1) < fPHOSDispersionCut) isDispOK = kTRUE;
636     else                                                        isDispOK = kFALSE; 
637   }
638   else{//EMCAL
639     
640     if(l0 > fEMCALL0CutMin && l0 < fEMCALL0CutMax) isDispOK = kTRUE;
641
642   }
643   
644   ph->SetDispBit(isDispOK) ;
645   
646   //TOF
647   Double_t tof=cluster->GetTOF()  ;
648   ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ; 
649   
650   //Charged 
651   Bool_t isNeutral = IsTrackMatched(cluster,cu,event);
652   
653   ph->SetChargedBit(isNeutral);
654   
655   //Set PID pdg
656   ph->SetIdentifiedParticleType(GetIdentifiedParticleType(cluster));
657   
658   if(fDebug > 0){ 
659     printf("AliCaloPID::SetPIDBits: TOF %e, Lambda0 %2.2f, Lambda1 %2.2f\n",tof , l0, l1);      
660     printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
661            ph->GetIdentifiedParticleType(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit()); 
662   }
663 }
664
665 //_________________________________________________________
666 Bool_t AliCaloPID::IsTrackMatched(AliVCluster* cluster,
667                                   AliCalorimeterUtils * cu, 
668                                   AliVEvent* event) const 
669 {
670   //Check if there is any track attached to this cluster
671   
672   Int_t nMatches = cluster->GetNTracksMatched();
673   AliVTrack * track = 0;
674   Double_t p[3];
675
676   if(nMatches > 0)
677   {
678     //In case of ESDs, by default without match one entry with negative index, no match, reject.
679     if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
680     {    
681       Int_t iESDtrack = cluster->GetTrackMatchedIndex();
682       if(iESDtrack >= 0) track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
683       else return kFALSE;
684       
685       if (!track)
686       {
687         if(fDebug > 0) printf("AliCaloPID::IsTrackMatched() - Null matched track in ESD when index is OK!\n");
688         return kFALSE;
689       }
690     }      
691     else { // AOD
692       track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0));
693       if (!track)
694       {
695         if(fDebug > 0) printf("AliCaloPID::IsTrackMatched() - Null matched track in AOD!\n");
696         return kFALSE;
697       }
698     }
699     
700     Float_t dZ  = cluster->GetTrackDz();
701     Float_t dR  = cluster->GetTrackDx();
702     
703     // if track matching was recalculated
704     if(cluster->IsEMCAL() && cu && cu->IsRecalculationOfClusterTrackMatchingOn())
705     {
706       dR = 2000., dZ = 2000.;
707       cu->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
708     }
709         
710     if(cluster->IsPHOS()) 
711     {
712       track->GetPxPyPz(p) ;
713       TLorentzVector trackmom(p[0],p[1],p[2],0);
714       Int_t charge = track->Charge();
715       Double_t mf  = event->GetMagneticField();
716       if(TestPHOSChargedVeto(dR, dZ, trackmom.Pt(), charge, mf ) < fPHOSRCut) return kTRUE;
717       else                                                                    return kFALSE;
718       
719     }//PHOS
720     else {//EMCAL
721       
722     if(fDebug > 0) 
723         printf("AliCaloPID::IsTrackMatched - EMCAL dR %f < %f, dZ %f < %f \n",dR, fEMCALDPhiCut, dZ, fEMCALDEtaCut);
724       
725       if(TMath::Abs(dR) < fEMCALDPhiCut && 
726          TMath::Abs(dZ) < fEMCALDEtaCut)   return kTRUE;
727       else                                 return kFALSE;
728       
729     }//EMCAL cluster 
730     
731     
732   } // more than 1 match, at least one track in array
733   else return kFALSE;
734     
735 }
736
737 //___________________________________________________________________________________________________
738 Float_t AliCaloPID::TestPHOSDispersion(const Double_t pt, const Double_t l1, const Double_t l2) const 
739 {
740   //Check if cluster photon-like. Uses photon cluster parameterization in real pp data 
741   //Returns distance in sigmas. Recommended cut 2.5
742   
743   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
744   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
745   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
746   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
747   Double_t c       =-0.35-0.550*TMath::Exp(-0.390730*pt) ;
748   Double_t r2      = 0.5*  (l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma + 
749                      0.5*  (l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
750                      0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
751   
752   if(fDebug > 0) printf("AliCaloPID::TestPHOSDispersion() - PHOS SS R %f < %f?\n", TMath::Sqrt(r2), fPHOSDispersionCut);
753   
754   return TMath::Sqrt(r2) ; 
755   
756 }
757
758 //_______________________________________________________________________________________________
759 Float_t AliCaloPID::TestPHOSChargedVeto(const Double_t dx,  const Double_t dz, const Double_t pt, 
760                                         const Int_t charge, const Double_t mf) const 
761 {
762   //Checks distance to the closest track. Takes into account 
763   //non-perpendicular incidence of tracks.
764   //returns distance in sigmas. Recommended cut: 2.
765   //Requires (sign) of magnetic filed. onc can find it for example as following
766   //  Double_t mf=0. ;
767   //  AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
768   //  if(event)
769   //    mf = event->GetMagneticField(); //Positive for ++ and negative for --
770   
771   
772   Double_t meanX = 0.;
773   Double_t meanZ = 0.;
774   Double_t sx = TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
775                            6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+
776                            1.59219);
777   Double_t sz = TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+
778                            1.60) ;
779   
780   if(mf<0.){ //field --
781     meanZ = -0.468318 ;
782     if(charge>0)
783       meanX = TMath::Min(7.3, 3.89994*1.20679 *1.20679 /(pt*pt+1.20679*1.20679)+  
784                          0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
785     else
786       meanX =-TMath::Min(7.7, 3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+
787                          1.23114 +4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
788   }
789   else{ //Field ++
790     meanZ = -0.468318;
791     if(charge>0)
792       meanX =-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+
793                          0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
794     else
795       meanX = TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-
796                          0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;     
797   }
798   
799   Double_t rz = (dz-meanZ)/sz ;
800   Double_t rx = (dx-meanX)/sx ;
801   
802   if(fDebug > 0) 
803     printf("AliCaloPID::TestPHOSDispersion() - PHOS Matching R %f < %f\n",TMath::Sqrt(rx*rx+rz*rz), fPHOSRCut);
804   
805   return TMath::Sqrt(rx*rx+rz*rz) ;
806   
807 }