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