AliPhysics  3337bb0 (3337bb0)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliCaloPID.cxx
Go to the documentation of this file.
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 // --- ROOT system ---
17 #include <TMath.h>
18 #include <TString.h>
19 #include <TList.h>
20 
21 // ---- ANALYSIS system ----
22 #include "AliCaloPID.h"
23 #include "AliAODCaloCluster.h"
24 #include "AliESDCaloCluster.h"
25 #include "AliVCaloCells.h"
26 #include "AliVTrack.h"
27 #include "AliCaloTrackParticle.h"
28 #include "AliCalorimeterUtils.h"
29 #include "AliFiducialCut.h" // detector enum definition
30 #include "AliVEvent.h"
31 #include "AliLog.h"
32 
33 // ---- Detector ----
34 #include "AliEMCALPIDUtils.h"
35 
37 ClassImp(AliCaloPID) ;
39 
40 //________________________
43 //________________________
45 TObject(), fDebug(-1), fParticleFlux(kLow),
46 //Bayesian
47 fEMCALPIDUtils(), fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
48 fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
49 fEMCALElectronWeight(0.), fEMCALChargeWeight(0.), fEMCALNeutralWeight(0.),
50 fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
51 fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) , fPHOSNeutralWeight(0.),
52 fPHOSWeightFormula(0), fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
53 fPHOSPhotonWeightFormulaExpression(""),
54 fPHOSPi0WeightFormulaExpression(""),
55 //PID calculation
56 fEMCALL0CutMax(100.), fEMCALL0CutMin(0),
57 fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.),
58 fEMCALUseTrackPtDepMatchingCut(0),
59 fEMCALFuncTrackPtDepDEta(0), fEMCALFuncTrackPtDepDPhi(0),
60 fEMCALFuncTrackPtDepDEtaString(""), fEMCALFuncTrackPtDepDPhiString(""),
61 fEMCALFuncTrackPtDepDEtaNParam(0) , fEMCALFuncTrackPtDepDPhiNParam(0),
62 fEMCALFuncTrackPtDepDEtaParam (0) , fEMCALFuncTrackPtDepDPhiParam (0),
63 fTOFCut(0.),
64 fPHOSDispersionCut(1000), fPHOSRCut(1000),
65 //Split
66 fUseSimpleMassCut(kFALSE),
67 fUseSimpleM02Cut(kFALSE),
68 fUseSplitAsyCut(kFALSE),
69 fUseSplitSSCut(kTRUE),
70 fSplitM02MaxCut(0), fSplitM02MinCut(0), fSplitMinNCells(0),
71 fMassEtaMin(0), fMassEtaMax(0),
72 fMassPi0Min(0), fMassPi0Max(0),
73 fMassPhoMin(0), fMassPhoMax(0),
74 fM02MaxParamShiftNLMN(0),
75 fSplitWidthSigma(0), fMassShiftHighECell(0)
76 {
78 }
79 
80 //________________________________________
85 //________________________________________
87 TObject(), fDebug(-1), fParticleFlux(flux),
88 //Bayesian
89 fEMCALPIDUtils(), fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
90 fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
91 fEMCALElectronWeight(0.), fEMCALChargeWeight(0.), fEMCALNeutralWeight(0.),
92 fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
93 fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) , fPHOSNeutralWeight(0.),
94 fPHOSWeightFormula(0), fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
95 fPHOSPhotonWeightFormulaExpression(""),
96 fPHOSPi0WeightFormulaExpression(""),
97 //PID calculation
98 fEMCALL0CutMax(100.), fEMCALL0CutMin(0),
99 fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.),
100 fEMCALUseTrackPtDepMatchingCut(0),
101 fEMCALFuncTrackPtDepDEta(0), fEMCALFuncTrackPtDepDPhi(0),
102 fEMCALFuncTrackPtDepDEtaString(""), fEMCALFuncTrackPtDepDPhiString(""),
103 fEMCALFuncTrackPtDepDEtaNParam(0) , fEMCALFuncTrackPtDepDPhiNParam(0),
104 fEMCALFuncTrackPtDepDEtaParam (0) , fEMCALFuncTrackPtDepDPhiParam (0),
105 fTOFCut(0.),
106 fPHOSDispersionCut(1000), fPHOSRCut(1000),
107 //Split
108 fUseSimpleMassCut(kFALSE),
109 fUseSimpleM02Cut(kFALSE),
110 fUseSplitAsyCut(kFALSE),
111 fUseSplitSSCut(kTRUE),
112 fSplitM02MaxCut(0), fSplitM02MinCut(0), fSplitMinNCells(0),
113 fMassEtaMin(0), fMassEtaMax(0),
114 fMassPi0Min(0), fMassPi0Max(0),
115 fMassPhoMin(0), fMassPhoMax(0),
116 fM02MaxParamShiftNLMN(0),
117 fSplitWidthSigma(0), fMassShiftHighECell(0)
118 {
119  InitParameters();
120 }
121 
122 //_______________________________________________
127 //_______________________________________________
128 AliCaloPID::AliCaloPID(const TNamed * emcalpid) :
129 TObject(), fDebug(-1), fParticleFlux(kLow),
130 //Bayesian
131 fEMCALPIDUtils((AliEMCALPIDUtils*)emcalpid),
132 fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
133 fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
134 fEMCALElectronWeight(0.), fEMCALChargeWeight(0.), fEMCALNeutralWeight(0.),
135 fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
136 fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) , fPHOSNeutralWeight(0.),
137 fPHOSWeightFormula(0), fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
138 fPHOSPhotonWeightFormulaExpression(""),
139 fPHOSPi0WeightFormulaExpression(""),
140 //PID calculation
141 fEMCALL0CutMax(100.), fEMCALL0CutMin(0),
142 fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.),
143 fEMCALUseTrackPtDepMatchingCut(0),
144 fEMCALFuncTrackPtDepDEta(0), fEMCALFuncTrackPtDepDPhi(0),
145 fEMCALFuncTrackPtDepDEtaString(""), fEMCALFuncTrackPtDepDPhiString(""),
146 fEMCALFuncTrackPtDepDEtaNParam(0) , fEMCALFuncTrackPtDepDPhiNParam(0),
147 fEMCALFuncTrackPtDepDEtaParam (0) , fEMCALFuncTrackPtDepDPhiParam (0),
148 fTOFCut(0.),
149 fPHOSDispersionCut(1000), fPHOSRCut(1000),
150 //Split
151 fUseSimpleMassCut(kFALSE),
152 fUseSimpleM02Cut(kFALSE),
153 fUseSplitAsyCut(kFALSE),
154 fUseSplitSSCut(kTRUE),
155 fSplitM02MaxCut(0), fSplitM02MinCut(0), fSplitMinNCells(0),
156 fMassEtaMin(0), fMassEtaMax(0),
157 fMassPi0Min(0), fMassPi0Max(0),
158 fMassPhoMin(0), fMassPhoMax(0),
159 fM02MaxParamShiftNLMN(0),
160 fSplitWidthSigma(0), fMassShiftHighECell(0)
161 {
162  InitParameters();
163 }
164 
165 //_______________________
166 // Destructor.
167 //_______________________
169 {
170  delete fPHOSPhotonWeightFormula ;
171  delete fPHOSPi0WeightFormula ;
172  delete fEMCALPIDUtils ;
177 }
178 
179 //_______________________________
180 // Initialize the parameters of the PID.
181 //_______________________________
183 {
184  // Bayesian
185  fEMCALPhotonWeight = 0.6 ;
186  fEMCALPi0Weight = 0.6 ;
187  fEMCALElectronWeight = 0.6 ;
188  fEMCALChargeWeight = 0.6 ;
189  fEMCALNeutralWeight = 0.6 ;
190 
191  fPHOSPhotonWeight = 0.6 ;
192  fPHOSPi0Weight = 0.6 ;
193  fPHOSElectronWeight = 0.6 ;
194  fPHOSChargeWeight = 0.6 ;
195  fPHOSNeutralWeight = 0.6 ;
196 
197  //Formula to set the PID weight threshold for photon or pi0
198  fPHOSWeightFormula = kFALSE;
199  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))";
200  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))" ;
201 
203  {
204  if(fParticleFlux == kLow)
205  {
206  AliInfo("SetLOWFluxParam");
207  fEMCALPIDUtils->SetLowFluxParam() ;
208  }
209  else if (fParticleFlux == kHigh)
210  {
211  AliInfo("SetHighFluxParam");
212  fEMCALPIDUtils->SetHighFluxParam() ;
213  }
214  }
215 
216  //PID recalculation, not bayesian
217 
218  //EMCAL
219  fEMCALL0CutMax = 0.3 ;
220  fEMCALL0CutMin = 0.01;
221 
222  // Fix Track Matching
223  fEMCALDPhiCut = 0.05; // Same cut as in AliEMCALRecoUtils
224  fEMCALDEtaCut = 0.025;// Same cut as in AliEMCALRecoUtils
225 
226  // Pt dependent track matching
228 
229  // PHOS / EMCAL, not used
230  fTOFCut = 1.e-6;
231 
232  //PHOS
233  fPHOSRCut = 2. ;
234  fPHOSDispersionCut = 2.5;
235 
236  // Cluster splitting
237 
238  fSplitM02MinCut = 0.3 ;
239  fSplitM02MaxCut = 5 ;
240  fSplitMinNCells = 4 ;
241 
242  fMassEtaMin = 0.4;
243  fMassEtaMax = 0.6;
244 
245  fMassPi0Min = 0.11;
246  fMassPi0Max = 0.18;
247 
248  fMassPhoMin = 0.0;
249  fMassPhoMax = 0.08;
250 
251  fMassPi0Param[0][0] = 0 ; // Constant term on mass dependence
252  fMassPi0Param[0][1] = 0 ; // slope term on mass dependence
253  fMassPi0Param[0][2] = 0 ; // E function change
254  fMassPi0Param[0][3] = 0.044 ; // constant term on mass dependence
255  fMassPi0Param[0][4] = 0.0049; // slope term on mass dependence
256  fMassPi0Param[0][5] = 0.070 ; // Absolute low mass cut
257 
258  fMassPi0Param[1][0] = 0.115 ; // Constant term below 21 GeV
259  fMassPi0Param[1][1] = 0.00096; // slope term below 21 GeV
260  fMassPi0Param[1][2] = 21 ; // E function change
261  fMassPi0Param[1][3] = 0.10 ; // constant term on mass dependence
262  fMassPi0Param[1][4] = 0.0017; // slope term on mass dependence
263  fMassPi0Param[1][5] = 0.070 ; // Absolute low mass cut
264 
265  fWidthPi0Param[0][0] = 0.012 ; // Constant term on width dependence
266  fWidthPi0Param[0][1] = 0.0 ; // Slope term on width dependence
267  fWidthPi0Param[0][2] = 19 ; // E function change
268  fWidthPi0Param[0][3] = 0.0012; // Constant term on width dependence
269  fWidthPi0Param[0][4] = 0.0006; // Slope term on width dependence
270  fWidthPi0Param[0][5] = 0.0 ; // xx term
271 
272  fWidthPi0Param[1][0] = 0.009 ; // Constant term on width dependence
273  fWidthPi0Param[1][1] = 0.000 ; // Slope term on width dependence
274  fWidthPi0Param[1][2] = 10 ; // E function change
275  fWidthPi0Param[1][3] = 0.0023 ; // Constant term on width dependence
276  fWidthPi0Param[1][4] = 0.00067; // Slope term on width dependence
277  fWidthPi0Param[1][5] = 0.000 ;// xx term
278 
279  fMassShiftHighECell = 0; // Shift of cuts in case of higher energy threshold in cells, 5 MeV when Ecell>150 MeV
280 
281  //TF1 *lM02MinNLM1 = new TF1("M02MinNLM1","exp(2.135-0.245*x)",6,13.6);
282  fM02MinParam[0][0] = 2.135 ;
283  fM02MinParam[0][1] =-0.245 ;
284  fM02MinParam[0][2] = 0.0 ;
285  fM02MinParam[0][3] = 0.0 ;
286  fM02MinParam[0][4] = 0.0 ;
287 
288  // Same as NLM=1 for NLM=2
289  fM02MinParam[1][0] = 2.135 ;
290  fM02MinParam[1][1] =-0.245 ;
291  fM02MinParam[1][2] = 0.0 ;
292  fM02MinParam[1][3] = 0.0 ;
293  fM02MinParam[1][4] = 0.0 ;
294 
295  //TF1 *lM02MaxNLM1 = new TF1("M02MaxNLM1","exp(0.0662-0.0201*x)-0.0955+0.00186*x[0]+9.91/x[0]",6,100);
296  fM02MaxParam[0][0] = 0.0662 ;
297  fM02MaxParam[0][1] =-0.0201 ;
298  fM02MaxParam[0][2] =-0.0955 ;
299  fM02MaxParam[0][3] = 0.00186;
300  fM02MaxParam[0][4] = 9.91 ;
301 
302  //TF1 *lM02MaxNLM2 = new TF1("M02MaxNLM2","exp(0.353-0.0264*x)-0.524+0.00559*x[0]+21.9/x[0]",6,100);
303  fM02MaxParam[1][0] = 0.353 ;
304  fM02MaxParam[1][1] =-0.0264 ;
305  fM02MaxParam[1][2] =-0.524 ;
306  fM02MaxParam[1][3] = 0.00559;
307  fM02MaxParam[1][4] = 21.9 ;
308 
309  fM02MaxParamShiftNLMN = 0.75;
310 
311  //TF1 *lAsyNLM1 = new TF1("lAsyNLM1","0.96-879/(x*x*x)",5,100);
312  fAsyMinParam[0][0] = 0.96 ;
313  fAsyMinParam[0][1] = 0 ;
314  fAsyMinParam[0][2] =-879 ;
315  fAsyMinParam[0][3] = 0.96 ; // Absolute max
316 
317  //TF1 *lAsyNLM2 = new TF1("lAsyNLM2","0.95+0.0015*x-233/(x*x*x)",5,100);
318  fAsyMinParam[1][0] = 0.95 ;
319  fAsyMinParam[1][1] = 0.0015;
320  fAsyMinParam[1][2] =-233 ;
321  fAsyMinParam[1][3] = 1.0 ; // Absolute max
322 
323  fSplitEFracMin[0] = 0.0 ; // 0.96
324  fSplitEFracMin[1] = 0.0 ; // 0.96
325  fSplitEFracMin[2] = 0.0 ; // 0.7
326 
327  fSubClusterEMin[0] = 0.0; // 3 GeV
328  fSubClusterEMin[1] = 0.0; // 1 GeV
329  fSubClusterEMin[2] = 0.0; // 1 GeV
330 
331  fSplitWidthSigma = 3. ;
332 }
333 
334 
335 //_______________________________
342 //_______________________________
344 {
345  fEMCALFuncTrackPtDepDEtaString = "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])" ;
346  fEMCALFuncTrackPtDepDPhiString = "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])" ;
347 
350 
353 
356 
363 }
364 
365 //_________________________________________________________________________________________
371 //_________________________________________________________________________________________
373 {
374  if(!fUseSplitAsyCut) return kTRUE ;
375 
376  Float_t abasy = TMath::Abs(asy);
377 
378  Int_t inlm = nlm-1;
379  if(nlm > 2) inlm=1; // only 2 cases defined nlm=1 and nlm>=2
380 
381  // Get the parametrized min cut of asymmetry for NLM=2 up to 11 GeV
382 
383  Float_t cut = fAsyMinParam[inlm][0] + fAsyMinParam[inlm][1]*energy + fAsyMinParam[inlm][2]/energy/energy/energy ;
384 
385  // In any case and beyond validity energy range of the function,
386  // the parameter cannot be smaller than 1
387  if( cut > fAsyMinParam[inlm][3] ) cut = fAsyMinParam[inlm][3];
388 
389  //printf("energy %2.2f - nlm: %d (%d)- p0 %f, p1 %f, p2 %f, p3 %f ; cut: %2.2f\n",energy,nlm,inlm,
390  // fAsyMinParam[inlm][0],fAsyMinParam[inlm][1],fAsyMinParam[inlm][2],fAsyMinParam[inlm][3],cut);
391 
392  if(abasy < cut) return kTRUE;
393  else return kFALSE;
394 }
395 
396 //______________________________________________________________________________________
402 //______________________________________________________________________________________
404 {
406  {
407  if(mass < fMassPi0Max && mass > fMassPi0Min) return kTRUE;
408  else return kFALSE;
409  }
410 
411  // Get the selected mean value as reference for the mass
412  Int_t inlm = nlm-1;
413  if(nlm > 2) inlm=1; // only 2 cases defined nlm=1 and nlm>=2
414 
415  Float_t meanMass = energy * fMassPi0Param[inlm][1] + fMassPi0Param[inlm][0];
416  if(energy > fMassPi0Param[inlm][2]) meanMass = energy * fMassPi0Param[inlm][4] + fMassPi0Param[inlm][3];
417 
418  // In case of higher energy cell cut than 50 MeV, smaller mean mass 0-5 MeV, not really necessary
419  meanMass -= fMassShiftHighECell;
420 
421  // Get the parametrized width of the mass
422  Float_t width = 0.009;
423  if (energy > 8 && energy < fWidthPi0Param[inlm][2])
424  width = energy * fWidthPi0Param[inlm][1] + fWidthPi0Param[inlm][0];
425  else if( energy > fWidthPi0Param[inlm][2])
426  width = energy * energy * fWidthPi0Param[inlm][5] + energy * fWidthPi0Param[inlm][4] + fWidthPi0Param[inlm][3];
427 
428  // Calculate the 2 sigma cut
429  Float_t minMass = meanMass-fSplitWidthSigma*width;
430  Float_t maxMass = meanMass+fSplitWidthSigma*width;
431 
432  // In case of low energy, hard cut to avoid conversions
433  if(energy < 10 && minMass < fMassPi0Param[inlm][5] ) minMass = fMassPi0Param[inlm][5];
434 
435  //printf("E %2.2f, mass %1.1f, nlm %d: sigma %1.1f width %3.1f, mean Mass %3.0f, minMass %3.0f, maxMass %3.0f\n ",
436  // energy,mass *1000, inlm, fSplitWidthSigma, width*1000, meanMass*1000,minMass*1000,maxMass*1000);
437 
438  if(mass < maxMass && mass > minMass) return kTRUE;
439  else return kFALSE;
440 }
441 
442 //________________________________________________
446 //________________________________________________
448 {
449  Float_t minCut = fSplitM02MinCut;
450  Float_t maxCut = fSplitM02MaxCut;
451 
452  if(m02 < maxCut && m02 > minCut) return kTRUE;
453  else return kFALSE;
454 }
455 
456 //_______________________________________________________________________________
462 //_______________________________________________________________________________
464 {
465  if(!fUseSplitSSCut) return kTRUE ;
466 
467  //First check the absolute minimum and maximum
468  if(!IsInM02Range(m02)) return kFALSE ;
469 
470  //If requested, check the E dependent cuts
471  else if(!fUseSimpleM02Cut)
472  {
473  Int_t inlm = nlm-1;
474  if(nlm > 2) inlm=1; // only 2 cases defined nlm=1 and nlm>=2
475 
476  Float_t minCut = fSplitM02MinCut;
477  Float_t maxCut = fSplitM02MaxCut;
478 
479  //e^{a+bx} + c + dx + e/x
480  if(energy > 1) minCut = TMath::Exp( fM02MinParam[inlm][0] + fM02MinParam[inlm][1]*energy ) +
481  fM02MinParam[inlm][2] + fM02MinParam[inlm][3]*energy + fM02MinParam[inlm][4]/energy;
482 
483  if(energy > 1) maxCut = TMath::Exp( fM02MaxParam[inlm][0] + fM02MaxParam[inlm][1]*energy ) +
484  fM02MaxParam[inlm][2] + fM02MaxParam[inlm][3]*energy + fM02MaxParam[inlm][4]/energy;
485 
486  // In any case and beyond validity energy range of the function,
487  // the parameter cannot be smaller than 0.3 or larger than 4-5
488  if( minCut < fSplitM02MinCut) minCut = fSplitM02MinCut;
489  if( maxCut > fSplitM02MaxCut) maxCut = fSplitM02MaxCut;
490  if( nlm > 2 ) maxCut+=fM02MaxParamShiftNLMN;
491 
492  //if(energy > 7) printf("\t \t E %2.2f, nlm %d, m02 %2.2f, minM02 %2.2f, maxM02 %2.2f\n",energy, nlm, m02,minCut,maxCut);
493 
494  if(m02 < maxCut && m02 > minCut) return kTRUE;
495  else return kFALSE;
496 
497  }
498 
499  else return kTRUE;
500 }
501 
502 
503 //______________________________________________________________________________
504 // Select the appropriate shower shape range in splitting method to select eta's
505 // Use same parametrization as pi0, just shift the distributions (to be tuned)
510 //______________________________________________________________________________
512 {
513  if(!fUseSplitSSCut) return kTRUE ;
514 
515  //First check the absolute minimum and maximum
516  if(!IsInM02Range(m02)) return kFALSE ;
517 
518  //DO NOT USE, study parametrization
519 
520  //If requested, check the E dependent cuts
521  else if(!fUseSimpleM02Cut)
522  {
523  Int_t inlm = nlm-1;
524  if(nlm > 2) inlm=1; // only 2 cases defined nlm=1 and nlm>=2
525 
526  Float_t minCut = fSplitM02MinCut;
527  Float_t maxCut = fSplitM02MaxCut;
528 
529  Float_t shiftE = energy-20; // to be tuned
530  if(nlm==1) shiftE=energy-28;
531 
532  //e^{a+bx} + c + dx + e/x
533  if(shiftE > 1) minCut = TMath::Exp( fM02MinParam[inlm][0] + fM02MinParam[inlm][1]*shiftE ) +
534  fM02MinParam[inlm][2] + fM02MinParam[inlm][3]*shiftE + fM02MinParam[inlm][4]/shiftE;
535 
536  // In any case the parameter cannot be smaller than 0.3
537  if( minCut < fSplitM02MinCut) minCut = fSplitM02MinCut;
538 
539  shiftE = energy+20; // to be tuned
540 
541  if(shiftE > 1) maxCut = 1 + TMath::Exp( fM02MaxParam[inlm][0] + fM02MaxParam[inlm][1]*shiftE ) +
542  fM02MaxParam[inlm][2] + fM02MaxParam[inlm][3]*shiftE + fM02MaxParam[inlm][4]/shiftE;
543 
544  // In any case the parameter cannot be smaller than 4-5
545  if( maxCut > fSplitM02MaxCut) maxCut = fSplitM02MaxCut;
546  if( nlm > 2 ) maxCut+=fM02MaxParamShiftNLMN;
547 
548  //if(energy>6)printf("\t \t E %2.2f, nlm %d, m02 %2.2f, minM02 %2.2f, maxM02 %2.2f\n",energy, nlm, m02,minCut,maxCut);
549 
550  if(m02 < maxCut && m02 > minCut) return kTRUE;
551  else return kFALSE;
552 
553  }
554 
555  else return kTRUE;
556 }
557 
558 //______________________________________________________________________________
565 //______________________________________________________________________________
567 {
568  if(!fUseSplitSSCut) return kTRUE ;
569 
570  Float_t minCut = 0.1;
571  Float_t maxCut = fSplitM02MinCut;
572 
573  if(!fUseSimpleM02Cut)
574  {
575  Int_t inlm = nlm-1;
576  if(nlm > 2) inlm=1; // only 2 cases defined nlm=1 and nlm>=2
577 
578  //e^{a+bx} + c + dx + e/x
579  if(energy > 1) maxCut = TMath::Exp( fM02MinParam[inlm][0] + fM02MinParam[inlm][1]*energy ) +
580  fM02MinParam[inlm][2] + fM02MinParam[inlm][3]*energy + fM02MinParam[inlm][4]/energy;
581 
582  if( maxCut < fSplitM02MinCut) maxCut = fSplitM02MinCut;
583  }
584 
585  if(m02 < maxCut && m02 > minCut) return kTRUE;
586  else return kFALSE;
587 
588 }
589 
590 //______________________________________________
592 //______________________________________________
593 AliEMCALPIDUtils *AliCaloPID::GetEMCALPIDUtils()
594 {
595  if(!fEMCALPIDUtils) fEMCALPIDUtils = new AliEMCALPIDUtils ;
596 
597  return fEMCALPIDUtils ;
598 }
599 
600 
601 //________________________________________________________________
604 //________________________________________________________________
606 {
607  Float_t energy = cluster->E();
608  Float_t lambda0 = cluster->GetM02();
609  Float_t lambda1 = cluster->GetM20();
610 
611  // ---------------------
612  // Use bayesian approach
613  // ---------------------
614 
616  {
617  Double_t weights[AliPID::kSPECIESCN];
618 
619  if(cluster->IsEMCAL() && fRecalculateBayesian)
620  {
621  fEMCALPIDUtils->ComputePID(energy, lambda0);
622  for(Int_t i = 0; i < AliPID::kSPECIESCN; i++) weights[i] = fEMCALPIDUtils->GetPIDFinal(i);
623  }
624  else
625  {
626  for(Int_t i = 0; i < AliPID::kSPECIESCN; i++) weights[i] = cluster->GetPID()[i];
627  }
628 
629  if(fDebug > 0) PrintClusterPIDWeights(weights);
630 
631  return GetIdentifiedParticleTypeFromBayesWeights(cluster->IsEMCAL(), weights, energy);
632  }
633 
634  // -------------------------------------------------------
635  // Calculate PID SS from data, do not use bayesian weights
636  // -------------------------------------------------------
637 
638  AliDebug(1,Form("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",
639  cluster->IsEMCAL(),energy,lambda0,cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(),
640  cluster->GetEmcCpvDistance(), cluster->GetDistanceToBadChannel(),cluster->GetNExMax()));
641 
642  if(cluster->IsEMCAL())
643  {
644  AliDebug(1,Form("EMCAL SS %f <%f < %f?",fEMCALL0CutMin, lambda0, fEMCALL0CutMax));
645 
646  if(lambda0 < fEMCALL0CutMax && lambda0 > fEMCALL0CutMin) return kPhoton ;
647  else return kNeutralUnknown ;
648  } // EMCAL
649  else // PHOS
650  {
651  if(TestPHOSDispersion(energy,lambda0,lambda1) < fPHOSDispersionCut) return kPhoton;
652  else return kNeutralUnknown;
653  }
654 }
655 
656 //_________________________________________________________________________________________________________
661 //_________________________________________________________________________________________________________
663 {
664  if(!pid)
665  {
666  AliFatal("pid pointer not initialized!!!");
667  return kNeutralUnknown; // not needed, added to content coverity
668  }
669 
670  Float_t wPh = fPHOSPhotonWeight ;
671  Float_t wPi0 = fPHOSPi0Weight ;
673  Float_t wCh = fPHOSChargeWeight ;
675 
676  if(!isEMCAL && fPHOSWeightFormula){
677  wPh = GetPHOSPhotonWeightFormula()->Eval(energy) ;
678  wPi0 = GetPHOSPi0WeightFormula() ->Eval(energy);
679  }
680  else
681  {
682  wPh = fEMCALPhotonWeight ;
683  wPi0 = fEMCALPi0Weight ;
684  wE = fEMCALElectronWeight ;
685  wCh = fEMCALChargeWeight ;
686  wNe = fEMCALNeutralWeight ;
687  }
688 
689  if(fDebug > 0) PrintClusterPIDWeights(pid);
690 
692  Float_t chargedHadronWeight = pid[AliVCluster::kProton]+pid[AliVCluster::kKaon]+
693  pid[AliVCluster::kPion]+pid[AliVCluster::kMuon];
694  Float_t neutralHadronWeight = pid[AliVCluster::kNeutron]+pid[AliVCluster::kKaon0];
695  Float_t allChargedWeight = pid[AliVCluster::kElectron]+pid[AliVCluster::kEleCon]+ chargedHadronWeight;
696  Float_t allNeutralWeight = pid[AliVCluster::kPhoton]+pid[AliVCluster::kPi0]+ neutralHadronWeight;
697 
698  //Select most probable ID
699  if(!isEMCAL) // PHOS
700  {
701  if(pid[AliVCluster::kPhoton] > wPh) pdg = kPhoton ;
702  else if(pid[AliVCluster::kPi0] > wPi0) pdg = kPi0 ;
703  else if(pid[AliVCluster::kElectron] > wE) pdg = kElectron ;
704  else if(pid[AliVCluster::kEleCon] > wE) pdg = kEleCon ;
705  else if(chargedHadronWeight > wCh) pdg = kChargedHadron ;
706  else if(neutralHadronWeight > wNe) pdg = kNeutralHadron ;
707  else if(allChargedWeight > allNeutralWeight)
708  pdg = kChargedUnknown ;
709  else
710  pdg = kNeutralUnknown ;
711  }
712  else //EMCAL
713  {
714  if(pid[AliVCluster::kPhoton] > wPh) pdg = kPhoton ;
715  else if(pid[AliVCluster::kElectron] > wE) pdg = kElectron ;
716  else if(pid[AliVCluster::kPhoton]+pid[AliVCluster::kElectron] > wPh) pdg = kPhoton ; //temporal sollution until track matching for electrons is considered
717  else if(pid[AliVCluster::kPi0] > wPi0) pdg = kPi0 ;
718  else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;
719  else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ;
720  else pdg = kNeutralUnknown ;
721  }
722 
723  AliDebug(1,Form("Final Pdg: %d, cluster energy %2.2f", pdg,energy));
724 
725  return pdg ;
726 
727 }
728 
729 //____________________________________________________________________________________________________________
748 //____________________________________________________________________________________________________________
750  AliVCaloCells* cells,
751  AliCalorimeterUtils * caloutils,
752  Double_t vertex[3],
753  Int_t & nMax,
754  Double_t & mass, Double_t & angle,
755  TLorentzVector & l1, TLorentzVector & l2,
756  Int_t & absId1, Int_t & absId2,
757  Float_t & distbad1, Float_t & distbad2,
758  Bool_t & fidcut1, Bool_t & fidcut2 ) const
759 {
760  Float_t eClus = cluster->E();
761  Float_t m02 = cluster->GetM02();
762  const Int_t nc = cluster->GetNCells();
763  Int_t absIdList[nc];
764  Float_t maxEList [nc];
765 
766  mass = -1.;
767  angle = -1.;
768 
769  //If too low number of cells, skip it
770  if ( nc < fSplitMinNCells) return kNeutralUnknown ;
771 
772  AliDebug(2,"\t pass nCells cut");
773 
774  // Get Number of local maxima
775  nMax = caloutils->GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ;
776 
777  AliDebug(1,Form("Cluster : E %1.1f, M02 %1.2f, NLM %d, N Cells %d",eClus,m02,nMax,nc));
778 
779  //---------------------------------------------------------------------
780  // Get the 2 max indeces and do inv mass
781  //---------------------------------------------------------------------
782 
784  if(cluster->IsPHOS()) calorimeter = AliFiducialCut::kPHOS;
785 
786  if ( nMax == 2 )
787  {
788  absId1 = absIdList[0];
789  absId2 = absIdList[1];
790 
791  //Order in energy
792  Float_t en1 = cells->GetCellAmplitude(absId1);
793  caloutils->RecalibrateCellAmplitude(en1,calorimeter,absId1);
794  Float_t en2 = cells->GetCellAmplitude(absId2);
795  caloutils->RecalibrateCellAmplitude(en2,calorimeter,absId2);
796  if(en1 < en2)
797  {
798  absId2 = absIdList[0];
799  absId1 = absIdList[1];
800  }
801  }
802  else if( nMax == 1 )
803  {
804 
805  absId1 = absIdList[0];
806 
807  //Find second highest energy cell
808 
809  Float_t enmax = 0 ;
810  for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++)
811  {
812  Int_t absId = cluster->GetCellsAbsId()[iDigit];
813  if( absId == absId1 ) continue ;
814  Float_t endig = cells->GetCellAmplitude(absId);
815  caloutils->RecalibrateCellAmplitude(endig,calorimeter,absId);
816  if(endig > enmax)
817  {
818  enmax = endig ;
819  absId2 = absId ;
820  }
821  }// cell loop
822  }// 1 maxima
823  else
824  { // n max > 2
825  // loop on maxima, find 2 highest
826 
827  // First max
828  Float_t enmax = 0 ;
829  for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++)
830  {
831  Float_t endig = maxEList[iDigit];
832  if(endig > enmax)
833  {
834  enmax = endig ;
835  absId1 = absIdList[iDigit];
836  }
837  }// first maxima loop
838 
839  // Second max
840  Float_t enmax2 = 0;
841  for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++)
842  {
843  if(absIdList[iDigit]==absId1) continue;
844  Float_t endig = maxEList[iDigit];
845  if(endig > enmax2)
846  {
847  enmax2 = endig ;
848  absId2 = absIdList[iDigit];
849  }
850  }// second maxima loop
851 
852  } // n local maxima > 2
853 
854  if(absId2<0 || absId1<0)
855  {
856  AliDebug(1,Form("Bad index for local maxima : N max %d, i1 %d, i2 %d, cluster E %2.2f, ncells %d, m02 %2.2f",
857  nMax,absId1,absId2,eClus,nc,m02));
858  return kNeutralUnknown ;
859  }
860 
861  //---------------------------------------------------------------------
862  // Split the cluster energy in 2, around the highest 2 local maxima
863  //---------------------------------------------------------------------
864 
865  AliAODCaloCluster cluster1(0, 0,NULL,0.,NULL,NULL,1,0);
866  AliAODCaloCluster cluster2(1, 0,NULL,0.,NULL,NULL,1,0);
867 
868  caloutils->SplitEnergy(absId1,absId2,cluster, cells, &cluster1, &cluster2,nMax); /*absIdList, maxEList,*/
869 
870  fidcut1 = caloutils->GetEMCALRecoUtils()->CheckCellFiducialRegion(caloutils->GetEMCALGeometry(), &cluster1,cells);
871  fidcut2 = caloutils->GetEMCALRecoUtils()->CheckCellFiducialRegion(caloutils->GetEMCALGeometry(), &cluster2,cells);
872 
873  caloutils->GetEMCALRecoUtils()->RecalculateClusterDistanceToBadChannel(caloutils->GetEMCALGeometry(),cells,&cluster1);
874  caloutils->GetEMCALRecoUtils()->RecalculateClusterDistanceToBadChannel(caloutils->GetEMCALGeometry(),cells,&cluster2);
875 
876  distbad1 = cluster1.GetDistanceToBadChannel();
877  distbad2 = cluster2.GetDistanceToBadChannel();
878 // if(!fidcut2 || !fidcut1 || distbad1 < 2 || distbad2 < 2)
879 // printf("*** Dist to bad channel cl %f, cl1 %f, cl2 %f; fid cut cl %d, cl1 %d, cl2 %d \n",
880 // cluster->GetDistanceToBadChannel(),distbad1,distbad2,
881 // caloutils->GetEMCALRecoUtils()->CheckCellFiducialRegion(caloutils->GetEMCALGeometry(), cluster,cells),fidcut1,fidcut2);
882 
883  cluster1.GetMomentum(l1,vertex);
884  cluster2.GetMomentum(l2,vertex);
885 
886  mass = (l1+l2).M();
887  angle = l2.Angle(l1.Vect());
888  Float_t e1 = cluster1.E();
889  Float_t e2 = cluster2.E();
890 
891  // Consider clusters with splitted energy not too different to original cluster energy
892  Float_t splitFracCut = 0;
893  if(nMax < 3) splitFracCut = fSplitEFracMin[nMax-1];
894  else splitFracCut = fSplitEFracMin[2];
895  if((e1+e2)/eClus < splitFracCut) return kNeutralUnknown ;
896 
897  AliDebug(2,"\t pass Split E frac cut");
898 
899  // Consider sub-clusters with minimum energy
900  Float_t minECut = fSubClusterEMin[2];
901  if (nMax == 2) minECut = fSubClusterEMin[1];
902  else if(nMax == 1) minECut = fSubClusterEMin[0];
903  if(e1 < minECut || e2 < minECut)
904  {
905  //printf("Reject: e1 %2.1f, e2 %2.1f, cut %2.1f\n",e1,e2,minECut);
906  return kNeutralUnknown ;
907  }
908 
909  AliDebug(2,"\t pass min sub-cluster E cut");
910 
911  // Asymmetry of cluster
912  Float_t asy =-10;
913  if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
914 
915  if( !IsInPi0SplitAsymmetryRange(eClus,asy,nMax) ) return kNeutralUnknown ;
916 
917 
918  AliDebug(2,"\t pass asymmetry cut");
919 
920  Bool_t pi0OK = kFALSE;
921  Bool_t etaOK = kFALSE;
922  Bool_t conOK = kFALSE;
923 
924  //If too small or big M02, skip it
925  if (IsInPi0M02Range(eClus,m02,nMax)) pi0OK = kTRUE;
926  else if(IsInEtaM02Range(eClus,m02,nMax)) etaOK = kTRUE;
927  else if(IsInConM02Range(eClus,m02,nMax)) conOK = kTRUE;
928 
929  Float_t energy = eClus;
930  if(nMax > 2) energy = e1+e2; // In case of NLM>2 use mass cut for NLM=2 but for the split sum not the cluster energy that is not the pi0 E.
931 
932  // Check the mass, and set an ID to the splitted cluster
933  if ( conOK && mass < fMassPhoMax && mass > fMassPhoMin ) { AliDebug(2,"\t Split Conv"); return kPhoton ; }
934  else if( etaOK && mass < fMassEtaMax && mass > fMassEtaMin ) { AliDebug(2,"\t Split Eta" ); return kEta ; }
935  else if( pi0OK && IsInPi0SplitMassRange(energy,mass,nMax) ) { AliDebug(2,"\t Split Pi0" ); return kPi0 ; }
936  else return kNeutralUnknown ;
937 
938 }
939 
940 //_________________________________________
942 //_________________________________________
944 {
945  TString parList ; //this will be list of parameters used for this analysis.
946  const Int_t buffersize = 255;
947  char onePar[buffersize] ;
948  snprintf(onePar,buffersize,"--- AliCaloPID ---") ;
949  parList+=onePar ;
951  {
952  snprintf(onePar,buffersize,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)",fEMCALPhotonWeight) ;
953  parList+=onePar ;
954  snprintf(onePar,buffersize,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)",fEMCALPi0Weight) ;
955  parList+=onePar ;
956  snprintf(onePar,buffersize,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)",fEMCALElectronWeight) ;
957  parList+=onePar ;
958  snprintf(onePar,buffersize,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)",fEMCALChargeWeight) ;
959  parList+=onePar ;
960  snprintf(onePar,buffersize,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)",fEMCALNeutralWeight) ;
961  parList+=onePar ;
962  snprintf(onePar,buffersize,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)",fPHOSPhotonWeight) ;
963  parList+=onePar ;
964  snprintf(onePar,buffersize,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)",fPHOSPi0Weight) ;
965  parList+=onePar ;
966  snprintf(onePar,buffersize,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)",fPHOSElectronWeight) ;
967  parList+=onePar ;
968  snprintf(onePar,buffersize,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)",fPHOSChargeWeight) ;
969  parList+=onePar ;
970  snprintf(onePar,buffersize,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)",fPHOSNeutralWeight) ;
971  parList+=onePar ;
972 
974  {
975  snprintf(onePar,buffersize,"PHOS Photon Weight Formula: %s",fPHOSPhotonWeightFormulaExpression.Data() ) ;
976  parList+=onePar;
977  snprintf(onePar,buffersize,"PHOS Pi0 Weight Formula: %s",fPHOSPi0WeightFormulaExpression.Data() ) ;
978  parList+=onePar;
979  }
980  }
981  else
982  {
983  snprintf(onePar,buffersize,"EMCAL: fEMCALL0CutMin =%2.2f, fEMCALL0CutMax =%2.2f (Cut on Shower Shape)",fEMCALL0CutMin, fEMCALL0CutMax) ;
984  parList+=onePar ;
985  snprintf(onePar,buffersize,"EMCAL: fEMCALDEtaCut =%2.2f, fEMCALDPhiCut =%2.2f (Cut on track matching)",fEMCALDEtaCut, fEMCALDPhiCut) ;
986  parList+=onePar ;
987  snprintf(onePar,buffersize,"fTOFCut =%e (Cut on TOF, used in PID evaluation)",fTOFCut) ;
988  parList+=onePar ;
989  snprintf(onePar,buffersize,"fPHOSRCut =%2.2f, fPHOSDispersionCut =%2.2f (Cut on Shower Shape and CPV)",fPHOSRCut,fPHOSDispersionCut) ;
990  parList+=onePar ;
991 
992  }
993 
994  if(fUseSimpleM02Cut)
995  {
996  snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f", fSplitM02MinCut, fSplitM02MaxCut) ;
997  parList+=onePar ;
998  }
999  snprintf(onePar,buffersize,"fMinNCells =%d", fSplitMinNCells) ;
1000  parList+=onePar ;
1001  if(fUseSimpleMassCut)
1002  {
1003  snprintf(onePar,buffersize,"pi0 : %2.1f < m <%2.1f", fMassPi0Min,fMassPi0Max);
1004  parList+=onePar ;
1005  }
1006  snprintf(onePar,buffersize,"eta : %2.1f < m <%2.1f", fMassEtaMin,fMassEtaMax);
1007  parList+=onePar ;
1008  snprintf(onePar,buffersize,"conv: %2.1f < m <%2.1f", fMassPhoMin,fMassPhoMax);
1009  parList+=onePar ;
1010 
1011 
1012  return parList;
1013 }
1014 
1015 //________________________________________________
1017 //________________________________________________
1018 void AliCaloPID::Print(const Option_t * opt) const
1019 {
1020  if(! opt)
1021  return;
1022 
1023  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1024 
1026  {
1027  printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",
1030  printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",
1033 
1034  printf("PHOS Parametrized weight on? = %d\n", fPHOSWeightFormula) ;
1035  if(fPHOSWeightFormula)
1036  {
1037  printf("Photon weight formula = %s\n", fPHOSPhotonWeightFormulaExpression.Data());
1038  printf("Pi0 weight formula = %s\n", fPHOSPi0WeightFormulaExpression .Data());
1039  }
1040  if(fRecalculateBayesian) printf(" Recalculate bayesian with Particle Flux? = %d\n",fParticleFlux);
1041  }
1042  else
1043  {
1044  printf("TOF cut = %e\n", fTOFCut);
1045  printf("EMCAL Lambda0 cut min = %2.2f; max = %2.2f\n", fEMCALL0CutMin,fEMCALL0CutMax);
1046  printf("EMCAL cluster-track dEta < %2.3f; dPhi < %2.3f\n", fEMCALDEtaCut, fEMCALDPhiCut);
1047  printf("PHOS Treac matching cut =%2.2f, Dispersion Cut =%2.2f \n",fPHOSRCut, fPHOSDispersionCut) ;
1048 
1049  }
1050 
1051  printf("Min. N Cells =%d \n", fSplitMinNCells) ;
1052  if(fUseSimpleM02Cut) printf("%2.2f < lambda_0^2 <%2.2f \n",fSplitM02MinCut,fSplitM02MaxCut);
1053  if(fUseSimpleMassCut)printf("pi0 : %2.2f<m<%2.2f \n", fMassPi0Min,fMassPi0Max);
1054  printf("eta : %2.2f<m<%2.2f \n", fMassEtaMin,fMassEtaMax);
1055  printf("phot: %2.2f<m<%2.2f \n", fMassPhoMin,fMassPhoMax);
1056 
1057  printf(" \n");
1058 }
1059 
1060 //_________________________________________________________________
1061 // Print PID of cluster, (AliVCluster*)cluster->GetPID()
1062 //_________________________________________________________________
1064 {
1065  printf("AliCaloPID::PrintClusterPIDWeights() \n \t ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, \n \t \
1066  pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
1067  pid[AliVCluster::kPhoton], pid[AliVCluster::kPi0],
1068  pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon],
1069  pid[AliVCluster::kPion], pid[AliVCluster::kKaon],
1070  pid[AliVCluster::kProton],
1071  pid[AliVCluster::kNeutron], pid[AliVCluster::kKaon0]);
1072 }
1073 
1074 //___________________________________________________________________________
1076 //___________________________________________________________________________
1077 void AliCaloPID::SetPIDBits(AliVCluster * cluster,
1079  AliVEvent* event)
1080 {
1081  // Dispersion/lambdas
1082  //Double_t disp= cluster->GetDispersion() ;
1083  Double_t l1 = cluster->GetM20() ;
1084  Double_t l0 = cluster->GetM02() ;
1085  Bool_t isDispOK = kTRUE ;
1086  if(cluster->IsPHOS()){
1087  if(TestPHOSDispersion(ph->Pt(),l0,l1) < fPHOSDispersionCut) isDispOK = kTRUE;
1088  else isDispOK = kFALSE;
1089  }
1090  else{//EMCAL
1091 
1092  if(l0 > fEMCALL0CutMin && l0 < fEMCALL0CutMax) isDispOK = kTRUE;
1093 
1094  }
1095 
1096  ph->SetDispBit(isDispOK) ;
1097 
1098  //TOF
1099  Double_t tof=cluster->GetTOF() ;
1100  ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
1101 
1102  //Charged
1103  Bool_t isNeutral = IsTrackMatched(cluster,cu,event);
1104 
1105  ph->SetChargedBit(isNeutral);
1106 
1107  //Set PID pdg
1109 
1110  AliDebug(1,Form("TOF %e, Lambda0 %2.2f, Lambda1 %2.2f",tof , l0, l1));
1111  AliDebug(1,Form("pdg %d, bits: TOF %d, Dispersion %d, Charge %d",
1112  ph->GetIdentifiedParticleType(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit()));
1113 }
1114 
1115 //_________________________________________________________
1122 //_________________________________________________________
1124 {
1125  if( !fEMCALFuncTrackPtDepDEta )
1126  {
1128 
1129  for(Int_t iparam = 0; iparam < fEMCALFuncTrackPtDepDEtaNParam; iparam++)
1130  fEMCALFuncTrackPtDepDEta->SetParameter(iparam,fEMCALFuncTrackPtDepDEtaParam[iparam]);
1131  }
1132 
1133  return fEMCALFuncTrackPtDepDEta ;
1134 }
1135 
1136 //_________________________________________________________
1143 //_________________________________________________________
1145 {
1146  if ( !fEMCALFuncTrackPtDepDPhi )
1147  {
1149 
1150  for(Int_t iparam = 0; iparam < fEMCALFuncTrackPtDepDPhiNParam; iparam++)
1151  fEMCALFuncTrackPtDepDPhi->SetParameter(iparam,fEMCALFuncTrackPtDepDPhiParam[iparam]);
1152  }
1153 
1154  return fEMCALFuncTrackPtDepDPhi;
1155 }
1156 
1157 //_________________________________________________________
1163 //_________________________________________________________
1164 Bool_t AliCaloPID::IsTrackMatched(AliVCluster* cluster,
1165  AliCalorimeterUtils * cu,
1166  AliVEvent* event)
1167 {
1168  Int_t nMatches = cluster->GetNTracksMatched();
1169  AliVTrack * track = 0;
1170 
1171  // At least one match
1172  //
1173  if(nMatches <= 0) return kFALSE;
1174 
1175  // Select the track, depending on ESD or AODs
1176  //
1177  //In case of ESDs,
1178  //by default without match one entry with negative index, no match, reject.
1179  //
1180  if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
1181  {
1182  Int_t iESDtrack = ((AliESDCaloCluster*)cluster)->GetTracksMatched()->At(0); //cluster->GetTrackMatchedIndex();
1183 
1184  if(iESDtrack >= 0) track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
1185  else return kFALSE;
1186 
1187  if (!track)
1188  {
1189  AliWarning(Form("Null matched track in ESD for index %d",iESDtrack));
1190  return kFALSE;
1191  }
1192  } // ESDs
1193  else
1194  { // AODs
1195  track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0));
1196  if (!track)
1197  {
1198  AliWarning("Null matched track in AOD!");
1199  return kFALSE;
1200  }
1201  } // AODs
1202 
1203  Float_t dEta = cluster->GetTrackDz();
1204  Float_t dPhi = cluster->GetTrackDx();
1205 
1206  // Comment out, new value already set in AliCalorimeterUtils::RecalculateClusterTrackMatching()
1207  // when executed in the reader.
1208  // // if track matching was recalculated
1209  // if(cluster->IsEMCAL() && cu && cu->IsRecalculationOfClusterTrackMatchingOn())
1210  // {
1211  // dR = 2000., dZ = 2000.;
1212  // cu->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1213  // //AliDebug(2,"Residuals, (Old, New): z (%2.4f,%2.4f), x (%2.4f,%2.4f)\n", cluster->GetTrackDz(),dZ,cluster->GetTrackDx(),dR));
1214  // }
1215 
1216  if(cluster->IsPHOS())
1217  {
1218  Int_t charge = track->Charge();
1219  Double_t mf = event->GetMagneticField();
1220  if(TestPHOSChargedVeto(dPhi, dEta, track->Pt(), charge, mf ) < fPHOSRCut)
1221  return kTRUE;
1222  else
1223  return kFALSE;
1224 
1225  } // PHOS
1226  else // EMCAL
1227  {
1228  AliDebug(1,Form("EMCAL dPhi %f < %f, dEta %f < %f ",dPhi, fEMCALDPhiCut, dEta, fEMCALDEtaCut));
1229 
1231  {
1232  if(TMath::Abs(dPhi) < fEMCALDPhiCut &&
1233  TMath::Abs(dEta) < fEMCALDEtaCut) return kTRUE;
1234  else return kFALSE;
1235  }
1236  else
1237  {
1238  Float_t trackPt = track->Pt();
1239 
1240  Bool_t matchDEta = kFALSE;
1241  if( TMath::Abs(dEta) < GetEMCALFuncTrackPtDepDEta()->Eval(trackPt))
1242  matchDEta = kTRUE;
1243  else
1244  matchDEta = kFALSE;
1245 
1246  Bool_t matchDPhi = kFALSE;
1247  if( TMath::Abs(dPhi) < GetEMCALFuncTrackPtDepDPhi()->Eval(trackPt))
1248  matchDPhi = kTRUE;
1249  else
1250  matchDPhi = kFALSE;
1251 
1252 // printf("Cluster E %2.2f, track pT %2.2f, dEta %2.2f, dPhi %2.2f, cut eta %2.2f, cut phi %2.2f, match eta %d, match phi %d\n",
1253 // cluster->E(),trackPt,dEta,dPhi,
1254 // GetEMCALFuncTrackPtDepDEta()->Eval(trackPt), GetEMCALFuncTrackPtDepDPhi()->Eval(trackPt),
1255 // matchDEta, matchDPhi);
1256 
1257  if(matchDPhi && matchDEta) return kTRUE ;
1258  else return kFALSE;
1259 
1260  }
1261  }// EMCAL cluster
1262 }
1263 
1264 //___________________________________________________________________________________________________
1272 //___________________________________________________________________________________________________
1274 {
1275  Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1276  Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1277  Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1278  Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1279  Double_t c =-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1280  Double_t r2 = 0.5* (l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1281  0.5* (l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1282  0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1283 
1284  AliDebug(1,Form("PHOS SS R %f < %f?", TMath::Sqrt(r2), fPHOSDispersionCut));
1285 
1286  return TMath::Sqrt(r2) ;
1287 }
1288 
1289 //_______________________________________________________________________________________________
1304 //_______________________________________________________________________________________________
1306  Int_t charge, Double_t mf) const
1307 {
1308  Double_t meanX = 0.;
1309  Double_t meanZ = 0.;
1310  Double_t sx = TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
1311  6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+
1312  1.59219);
1313  Double_t sz = TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+
1314  1.60) ;
1315 
1316  if(mf<0.){ //field --
1317  meanZ = -0.468318 ;
1318  if(charge>0)
1319  meanX = TMath::Min(7.3, 3.89994*1.20679 *1.20679 /(pt*pt+1.20679*1.20679)+
1320  0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
1321  else
1322  meanX =-TMath::Min(7.7, 3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+
1323  1.23114 +4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
1324  }
1325  else{ //Field ++
1326  meanZ = -0.468318;
1327  if(charge>0)
1328  meanX =-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+
1329  0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
1330  else
1331  meanX = TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-
1332  0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
1333  }
1334 
1335  Double_t rz = (dz-meanZ)/sz ;
1336  Double_t rx = (dx-meanX)/sx ;
1337 
1338  AliDebug(1,Form("PHOS Matching R %f < %f",TMath::Sqrt(rx*rx+rz*rz), fPHOSRCut));
1339 
1340  return TMath::Sqrt(rx*rx+rz*rz) ;
1341 }
1342 
Int_t charge
void InitParamTrackMatchPtDependent()
Definition: AliCaloPID.cxx:343
virtual void SetTOFBit(Bool_t tof)
Int_t pdg
TString fPHOSPhotonWeightFormulaExpression
Photon weight formula in string.
Definition: AliCaloPID.h:333
Bool_t IsInEtaM02Range(Float_t energy, Float_t m02, Int_t nlm) const
Definition: AliCaloPID.cxx:511
virtual Double_t Pt() const
Float_t fTOFCut
Cut on TOF, used in PID evaluation.
Definition: AliCaloPID.h:358
Float_t fMassPhoMin
Min Photon mass.
Definition: AliCaloPID.h:376
void SetPIDBits(AliVCluster *cluster, AliCaloTrackParticle *aodph, AliCalorimeterUtils *cu, AliVEvent *event)
Set Bits for PID selection.
TFormula * GetPHOSPi0WeightFormula()
Definition: AliCaloPID.h:167
double Double_t
Definition: External.C:58
Int_t fEMCALFuncTrackPtDepDEtaNParam
number of formula parameters for matching eta residual pT track dependent
Definition: AliCaloPID.h:348
Float_t fPHOSNeutralWeight
Bayesian PID weight for neutral hadrons in PHOS.
Definition: AliCaloPID.h:328
Bool_t fUseBayesianWeights
Select clusters based on weights calculated in reconstruction.
Definition: AliCaloPID.h:316
Bool_t IsInPi0SplitAsymmetryRange(Float_t energy, Float_t asy, Int_t nlm) const
Definition: AliCaloPID.cxx:372
Float_t fMassPhoMax
Min Photon mass.
Definition: AliCaloPID.h:377
Float_t fMassPi0Max
Min Pi0 mass, simple cut case.
Definition: AliCaloPID.h:375
Float_t fEMCALElectronWeight
Bayesian PID weight for electrons in EMCAL.
Definition: AliCaloPID.h:321
virtual ~AliCaloPID()
Definition: AliCaloPID.cxx:168
AliEMCALRecoUtils * GetEMCALRecoUtils() const
Float_t fPHOSPhotonWeight
Bayesian PID weight for photons in PHOS.
Definition: AliCaloPID.h:324
Float_t fSplitEFracMin[3]
Definition: AliCaloPID.h:384
Float_t fM02MaxParamShiftNLMN
shift of max M02 for NLM>2.
Definition: AliCaloPID.h:382
Float_t fWidthPi0Param[2][6]
Width param, 2 regions in energy.
Definition: AliCaloPID.h:379
Double_t mass
energy
Definition: HFPtSpectrum.C:44
Bool_t fUseSplitAsyCut
Remove splitted clusters with too large asymmetry.
Definition: AliCaloPID.h:367
TF1 * fEMCALFuncTrackPtDepDEta
TF1 for track pT dependent cut in matching eta residual.
Definition: AliCaloPID.h:344
Float_t fEMCALPi0Weight
Bayesian PID weight for pi0 in EMCAL.
Definition: AliCaloPID.h:320
Float_t fEMCALPhotonWeight
Bayesian PID weight for photons in EMCAL.
Definition: AliCaloPID.h:319
Float_t fSubClusterEMin[3]
Do not use sub-clusters with too low energy depeding on NLM.
Definition: AliCaloPID.h:386
Bool_t IsTrackMatched(AliVCluster *cluster, AliCalorimeterUtils *cu, AliVEvent *event)
Float_t fMassEtaMax
Max Eta mass.
Definition: AliCaloPID.h:373
TCanvas * c
Definition: TestFitELoss.C:172
TString fEMCALFuncTrackPtDepDPhiString
TF1 for track pT dependent cut in matching phi residual, formula string.
Definition: AliCaloPID.h:347
TString fPHOSPi0WeightFormulaExpression
Pi0 weight formula in string.
Definition: AliCaloPID.h:334
TFormula * GetPHOSPhotonWeightFormula()
Definition: AliCaloPID.h:161
Float_t fEMCALDEtaCut
Track matching fixed cut on eta residual.
Definition: AliCaloPID.h:340
const TString calorimeter
Definition: anaM.C:36
virtual void SetIdentifiedParticleType(Int_t pdg)
Float_t fMassPi0Min
Min Pi0 mass, simple cut case.
Definition: AliCaloPID.h:374
Float_t fPHOSChargeWeight
Bayesian PID weight for charged hadrons in PHOS.
Definition: AliCaloPID.h:327
virtual Bool_t GetDispBit() const
Float_t fMassPi0Param[2][6]
Mean mass param, 2 regions in energy.
Definition: AliCaloPID.h:378
TString GetPIDParametersList()
Put data member values in string to keep in output container.
Definition: AliCaloPID.cxx:943
void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
Float_t fM02MaxParam[2][5]
5 param for expo + pol fit on M02 maximum for pi0 selection.
Definition: AliCaloPID.h:381
Bool_t fRecalculateBayesian
Recalculate PID bayesian or use simple PID?
Definition: AliCaloPID.h:317
TF1 * GetEMCALFuncTrackPtDepDEta()
virtual Bool_t GetChargedBit() const
AliEMCALGeometry * GetEMCALGeometry() const
int Int_t
Definition: External.C:63
Float_t fSplitWidthSigma
Cut on mass+-width*fSplitWidthSigma.
Definition: AliCaloPID.h:387
Container for input particle information on CaloTrackCorr package.
Int_t fDebug
Debug level.
Definition: AliCaloPID.h:310
float Float_t
Definition: External.C:68
Bool_t IsInPi0M02Range(Float_t energy, Float_t m02, Int_t nlm) const
Definition: AliCaloPID.cxx:463
Int_t fEMCALFuncTrackPtDepDPhiNParam
number of formula parameters for matching eta residual pT track dependent
Definition: AliCaloPID.h:349
Float_t fSplitM02MaxCut
Study clusters with l0 smaller than cut.
Definition: AliCaloPID.h:369
Float_t * fEMCALFuncTrackPtDepDPhiParam
Formula parameters for track matching eta residual pT track dependent.
Definition: AliCaloPID.h:355
Float_t TestPHOSDispersion(Double_t pt, Double_t m20, Double_t m02) const
Bool_t fPHOSWeightFormula
Use parametrized weight threshold, function of energy.
Definition: AliCaloPID.h:330
Float_t TestPHOSChargedVeto(Double_t dx, Double_t dz, Double_t ptTrack, Int_t chargeTrack, Double_t mf) const
Float_t fEMCALNeutralWeight
Bayesian PID weight for neutral hadrons in EMCAL.
Definition: AliCaloPID.h:323
Int_t GetNumberOfLocalMaxima(AliVCluster *cluster, AliVCaloCells *cells)
Find the number of local maxima in cluster.
Float_t fPHOSRCut
Track-Cluster distance cut for track matching in PHOS.
Definition: AliCaloPID.h:361
Int_t GetIdentifiedParticleType(AliVCluster *cluster)
Definition: AliCaloPID.cxx:605
virtual void SetChargedBit(Bool_t ch)
TF1 * GetEMCALFuncTrackPtDepDPhi()
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
TString fEMCALFuncTrackPtDepDEtaString
TF1 for track pT dependent cut in matching eta residual, formula string.
Definition: AliCaloPID.h:346
Bool_t fUseSplitSSCut
Remove splitted clusters out of shower shape band.
Definition: AliCaloPID.h:368
Bool_t fEMCALUseTrackPtDepMatchingCut
Activate the matching selection, pT dependent.
Definition: AliCaloPID.h:343
Float_t fMassShiftHighECell
Shift cuts 5 MeV for Ecell > 150 MeV, default Ecell > 50 MeV.
Definition: AliCaloPID.h:388
Bool_t IsInPi0SplitMassRange(Float_t energy, Float_t mass, Int_t nlm) const
Definition: AliCaloPID.cxx:403
Float_t fPHOSElectronWeight
Bayesian PID weight for electrons in PHOS.
Definition: AliCaloPID.h:326
Float_t fEMCALDPhiCut
Track matching fixed cut on phi residual.
Definition: AliCaloPID.h:341
virtual Bool_t GetTOFBit() const
Float_t fEMCALL0CutMax
Max Cut on shower shape lambda0, used in PID evaluation, only EMCAL.
Definition: AliCaloPID.h:338
Int_t GetIdentifiedParticleTypeFromClusterSplitting(AliVCluster *cluster, AliVCaloCells *cells, AliCalorimeterUtils *caloutils, Double_t vertex[3], Int_t &nLocMax, Double_t &mass, Double_t &angle, TLorentzVector &l1, TLorentzVector &l2, Int_t &absId1, Int_t &absId2, Float_t &distbad1, Float_t &distbad2, Bool_t &fidcut1, Bool_t &fidcut2) const
Definition: AliCaloPID.cxx:749
Float_t * fEMCALFuncTrackPtDepDEtaParam
Formula parameters for track matching eta residual pT track dependent.
Definition: AliCaloPID.h:352
TFormula * fPHOSPhotonWeightFormula
Formula for photon weight.
Definition: AliCaloPID.h:331
Double_t minMass
Bool_t IsInConM02Range(Float_t energy, Float_t m02, Int_t nlm) const
Definition: AliCaloPID.cxx:566
Float_t fPHOSDispersionCut
Shower shape elipse radious cut.
Definition: AliCaloPID.h:360
Int_t GetIdentifiedParticleTypeFromBayesWeights(Bool_t isEMCAL, Double_t *pid, Float_t energy)
Definition: AliCaloPID.cxx:662
void RecalibrateCellAmplitude(Float_t &amp, Int_t calo, Int_t absId) const
Recalculate cell energy if recalibration factor.
virtual void SetDispBit(Bool_t disp)
const char Option_t
Definition: External.C:48
Bool_t fUseSimpleM02Cut
Use simple min-max M02 cut.
Definition: AliCaloPID.h:366
Class for PID selection with calorimeters.
Definition: AliCaloPID.h:53
Float_t fEMCALL0CutMin
Min Cut on shower shape lambda0, used in PID evaluation, only EMCAL.
Definition: AliCaloPID.h:339
AliEMCALPIDUtils * GetEMCALPIDUtils()
Definition: AliCaloPID.cxx:593
Double_t maxMass
Float_t fMassEtaMin
Min Eta mass.
Definition: AliCaloPID.h:372
bool Bool_t
Definition: External.C:53
Class with utils specific to calorimeter clusters/cells.
Float_t fEMCALChargeWeight
Bayesian PID weight for charged hadrons in EMCAL.
Definition: AliCaloPID.h:322
AliEMCALPIDUtils * fEMCALPIDUtils
Pointer to EMCALPID to redo the PID Bayesian calculation.
Definition: AliCaloPID.h:315
void InitParameters()
Definition: AliCaloPID.cxx:182
Int_t fParticleFlux
Particle flux for setting PID parameters.
Definition: AliCaloPID.h:311
Float_t fM02MinParam[2][5]
5 param for expo + pol fit on M02 minimum for pi0 selection (maximum for conversions).
Definition: AliCaloPID.h:380
void PrintClusterPIDWeights(const Double_t *pid) const
Bool_t fUseSimpleMassCut
Use simple min-max pi0 mass cut.
Definition: AliCaloPID.h:365
Float_t fPHOSPi0Weight
Bayesian PID weight for pi0 in PHOS.
Definition: AliCaloPID.h:325
virtual Int_t GetIdentifiedParticleType() const
TFormula * fPHOSPi0WeightFormula
Formula for pi0 weight.
Definition: AliCaloPID.h:332
Bool_t IsInM02Range(Float_t m02) const
Definition: AliCaloPID.cxx:447
Int_t fSplitMinNCells
Study clusters with ncells larger than cut.
Definition: AliCaloPID.h:371
void SplitEnergy(Int_t absId1, Int_t absId2, AliVCluster *cluster, AliVCaloCells *cells, AliAODCaloCluster *cluster1, AliAODCaloCluster *cluster2, Int_t nMax, Int_t eventNumber=0)
Float_t fSplitM02MinCut
Study clusters with l0 larger than cut, simple case.
Definition: AliCaloPID.h:370
Float_t fAsyMinParam[2][4]
4 param for fit on asymmetry minimum, for 2 cases, NLM=1 and NLM>=2.
Definition: AliCaloPID.h:383
TF1 * fEMCALFuncTrackPtDepDPhi
TF1 for track pT dependent cut in matching phi residual.
Definition: AliCaloPID.h:345