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