1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //_________________________________________________________________________
17 // Class for PID selection with calorimeters
18 // The Output of the main method GetIdentifiedParticleType is a PDG number identifying the cluster,
19 // being kPhoton, kElectron, kPi0 ... as defined in the header file
20 // - GetIdentifiedParticleType(const AliVCluster * cluster)
21 // Assignes a PID tag to the cluster, right now there is the possibility to : use bayesian weights from reco,
22 // recalculate them (EMCAL) or use other procedures not used in reco.
23 // In order to recalculate Bayesian, it is necessary to load the EMCALUtils library
24 // and do SwitchOnBayesianRecalculation().
25 // To change the PID parameters from Low to High like the ones by default, use the constructor
27 // where flux is AliCaloPID::kLow or AliCaloPID::kHigh
28 // If it is necessary to change the parameters use the constructor
29 // AliCaloPID(AliEMCALPIDUtils *utils) and set the parameters before.
31 // - GetGetIdentifiedParticleTypeFromBayesian(const Double_t * pid, const Float_t energy)
32 // Reads the PID weights array of the ESDs and depending on its magnitude identifies the particle,
33 // executed when bayesian is ON by GetIdentifiedParticleType(const AliVCluster * cluster)
34 // - SetPIDBits: Simple PID, depending on the thresholds fLOCut fTOFCut and even the
35 // result of the PID bayesian a different PID bit is set.
37 // - IsTrackMatched(): Independent method that needs to be combined with GetIdentifiedParticleType to know if the cluster was matched
39 //*-- Author: Gustavo Conesa (INFN-LNF)
40 //////////////////////////////////////////////////////////////////////////////
43 // --- ROOT system ---
48 // ---- ANALYSIS system ----
49 #include "AliCaloPID.h"
50 #include "AliAODCaloCluster.h"
51 #include "AliVCaloCells.h"
52 #include "AliVTrack.h"
53 #include "AliAODPWG4Particle.h"
54 #include "AliCalorimeterUtils.h"
55 #include "AliVEvent.h"
58 #include "AliEMCALPIDUtils.h"
63 //________________________
64 AliCaloPID::AliCaloPID() :
65 TObject(), fDebug(-1), fParticleFlux(kLow),
67 fEMCALPIDUtils(), fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
68 fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
69 fEMCALElectronWeight(0.), fEMCALChargeWeight(0.), fEMCALNeutralWeight(0.),
70 fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
71 fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) , fPHOSNeutralWeight(0.),
72 fPHOSWeightFormula(0), fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
73 fPHOSPhotonWeightFormulaExpression(""),
74 fPHOSPi0WeightFormulaExpression(""),
76 fEMCALL0CutMax(100.), fEMCALL0CutMin(0),
77 fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.),
79 fPHOSDispersionCut(1000), fPHOSRCut(1000),
80 fDoClusterSplitting(kFALSE),
81 fSplitM02MaxCut(0), fSplitM02MinCut(0), fSplitMinNCells(0),
82 fMassEtaMin(0), fMassEtaMax(0),
83 fMassPi0Min(0), fMassPi0Max(0),
84 fMassPhoMin(0), fMassPhoMax(0)
88 //Initialize parameters
92 //________________________________________
93 AliCaloPID::AliCaloPID(const Int_t flux) :
94 TObject(), fDebug(-1), fParticleFlux(flux),
96 fEMCALPIDUtils(), fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
97 fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
98 fEMCALElectronWeight(0.), fEMCALChargeWeight(0.), fEMCALNeutralWeight(0.),
99 fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
100 fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) , fPHOSNeutralWeight(0.),
101 fPHOSWeightFormula(0), fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
102 fPHOSPhotonWeightFormulaExpression(""),
103 fPHOSPi0WeightFormulaExpression(""),
105 fEMCALL0CutMax(100.), fEMCALL0CutMin(0),
106 fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.),
108 fPHOSDispersionCut(1000), fPHOSRCut(1000),
109 fDoClusterSplitting(kFALSE),
110 fSplitM02MaxCut(0), fSplitM02MinCut(0), fSplitMinNCells(0),
111 fMassEtaMin(0), fMassEtaMax(0),
112 fMassPi0Min(0), fMassPi0Max(0),
113 fMassPhoMin(0), fMassPhoMax(0)
117 //Initialize parameters
122 //_______________________________________________
123 AliCaloPID::AliCaloPID(const TNamed * emcalpid) :
124 TObject(), fDebug(-1), fParticleFlux(kLow),
126 fEMCALPIDUtils((AliEMCALPIDUtils*)emcalpid),
127 fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
128 fEMCALPhotonWeight(0.), fEMCALPi0Weight(0.),
129 fEMCALElectronWeight(0.), fEMCALChargeWeight(0.), fEMCALNeutralWeight(0.),
130 fPHOSPhotonWeight(0.), fPHOSPi0Weight(0.),
131 fPHOSElectronWeight(0.), fPHOSChargeWeight(0.) , fPHOSNeutralWeight(0.),
132 fPHOSWeightFormula(0), fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
133 fPHOSPhotonWeightFormulaExpression(""),
134 fPHOSPi0WeightFormulaExpression(""),
136 fEMCALL0CutMax(100.), fEMCALL0CutMin(0),
137 fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.),
139 fPHOSDispersionCut(1000), fPHOSRCut(1000),
140 fDoClusterSplitting(kFALSE),
141 fSplitM02MaxCut(0), fSplitM02MinCut(0), fSplitMinNCells(0),
142 fMassEtaMin(0), fMassEtaMax(0),
143 fMassPi0Min(0), fMassPi0Max(0),
144 fMassPhoMin(0), fMassPhoMax(0)
148 //Initialize parameters
152 //_______________________
153 AliCaloPID::~AliCaloPID()
157 delete fPHOSPhotonWeightFormula ;
158 delete fPHOSPi0WeightFormula ;
159 delete fEMCALPIDUtils ;
163 //_______________________________
164 void AliCaloPID::InitParameters()
166 //Initialize the parameters of the PID.
169 fEMCALPhotonWeight = 0.6 ;
170 fEMCALPi0Weight = 0.6 ;
171 fEMCALElectronWeight = 0.6 ;
172 fEMCALChargeWeight = 0.6 ;
173 fEMCALNeutralWeight = 0.6 ;
175 fPHOSPhotonWeight = 0.6 ;
176 fPHOSPi0Weight = 0.6 ;
177 fPHOSElectronWeight = 0.6 ;
178 fPHOSChargeWeight = 0.6 ;
179 fPHOSNeutralWeight = 0.6 ;
181 //Formula to set the PID weight threshold for photon or pi0
182 fPHOSWeightFormula = kFALSE;
183 fPHOSPhotonWeightFormulaExpression = "0.98*(x<40)+ 0.68*(x>=100)+(x>=40 && x<100)*(0.98+x*(6e-3)-x*x*(2e-04)+x*x*x*(1.1e-06))";
184 fPHOSPi0WeightFormulaExpression = "0.98*(x<65)+ 0.915*(x>=100)+(x>=65 && x-x*(1.95e-3)-x*x*(4.31e-05)+x*x*x*(3.61e-07))" ;
186 if(fRecalculateBayesian){
187 if(fParticleFlux == kLow){
188 printf("AliCaloPID::Init() - SetLOWFluxParam\n");
189 fEMCALPIDUtils->SetLowFluxParam() ;
191 else if (fParticleFlux == kHigh){
192 printf("AliCaloPID::Init() - SetHIGHFluxParam\n");
193 fEMCALPIDUtils->SetHighFluxParam() ;
197 //PID recalculation, not bayesian
200 fEMCALL0CutMax = 0.3 ;
201 fEMCALL0CutMin = 0.01;
203 fEMCALDPhiCut = 0.05; // Same cut as in AliEMCALRecoUtils
204 fEMCALDEtaCut = 0.025;// Same cut as in AliEMCALRecoUtils
206 // PHOS / EMCAL, not used
211 fPHOSDispersionCut = 2.5;
215 fSplitM02MinCut = 0.26 ;
216 fSplitM02MaxCut = 100 ;
217 fSplitMinNCells = 4 ;
230 //______________________________________________
231 AliEMCALPIDUtils *AliCaloPID::GetEMCALPIDUtils()
233 // return pointer to AliEMCALPIDUtils, create it if needed
235 if(!fEMCALPIDUtils) fEMCALPIDUtils = new AliEMCALPIDUtils ;
236 return fEMCALPIDUtils ;
241 //______________________________________________________________________
242 Int_t AliCaloPID::GetIdentifiedParticleType(const AliVCluster * cluster)
244 // Returns a PDG number corresponding to the likely ID of the cluster
246 Float_t energy = cluster->E();
247 Float_t lambda0 = cluster->GetM02();
248 Float_t lambda1 = cluster->GetM20();
250 // ---------------------
251 // Use bayesian approach
252 // ---------------------
254 if(fUseBayesianWeights)
256 Double_t weights[AliPID::kSPECIESN];
258 if(cluster->IsEMCAL() && fRecalculateBayesian)
260 fEMCALPIDUtils->ComputePID(energy, lambda0);
261 for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = fEMCALPIDUtils->GetPIDFinal(i);
265 for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = cluster->GetPID()[i];
268 if(fDebug > 0) PrintClusterPIDWeights(weights);
270 return GetIdentifiedParticleTypeFromBayesWeights(cluster->IsEMCAL(), weights, energy);
273 // -------------------------------------------------------
274 // Calculate PID SS from data, do not use bayesian weights
275 // -------------------------------------------------------
277 if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType: EMCAL %d?, E %3.2f, l0 %3.2f, l1 %3.2f, disp %3.2f, tof %1.11f, distCPV %3.2f, distToBC %1.1f, NMax %d\n",
278 cluster->IsEMCAL(),energy,lambda0,cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(),
279 cluster->GetEmcCpvDistance(), cluster->GetDistanceToBadChannel(),cluster->GetNExMax());
281 if(cluster->IsEMCAL())
283 if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType() - EMCAL SS %f <%f < %f?\n",fEMCALL0CutMin, lambda0, fEMCALL0CutMax);
285 if(lambda0 < fEMCALL0CutMax && lambda0 > fEMCALL0CutMin) return kPhoton ;
286 else return kNeutralUnknown ;
290 if(TestPHOSDispersion(energy,lambda0,lambda1) < fPHOSDispersionCut) return kPhoton;
291 else return kNeutralUnknown;
296 //_______________________________________________________________________________
297 Int_t AliCaloPID::GetIdentifiedParticleTypeFromBayesWeights(const Bool_t isEMCAL,
298 const Double_t * pid,
299 const Float_t energy)
301 //Return most probable identity of the particle after bayesian weights calculated in reconstruction
305 printf("AliCaloPID::GetIdentifiedParticleType() - pid pointer not initialized!!!\n");
309 Float_t wPh = fPHOSPhotonWeight ;
310 Float_t wPi0 = fPHOSPi0Weight ;
311 Float_t wE = fPHOSElectronWeight ;
312 Float_t wCh = fPHOSChargeWeight ;
313 Float_t wNe = fPHOSNeutralWeight ;
315 if(!isEMCAL && fPHOSWeightFormula){
316 wPh = GetPHOSPhotonWeightFormula()->Eval(energy) ;
317 wPi0 = GetPHOSPi0WeightFormula() ->Eval(energy);
321 wPh = fEMCALPhotonWeight ;
322 wPi0 = fEMCALPi0Weight ;
323 wE = fEMCALElectronWeight ;
324 wCh = fEMCALChargeWeight ;
325 wNe = fEMCALNeutralWeight ;
328 if(fDebug > 0) PrintClusterPIDWeights(pid);
330 Int_t pdg = kNeutralUnknown ;
331 Float_t chargedHadronWeight = pid[AliVCluster::kProton]+pid[AliVCluster::kKaon]+
332 pid[AliVCluster::kPion]+pid[AliVCluster::kMuon];
333 Float_t neutralHadronWeight = pid[AliVCluster::kNeutron]+pid[AliVCluster::kKaon0];
334 Float_t allChargedWeight = pid[AliVCluster::kElectron]+pid[AliVCluster::kEleCon]+ chargedHadronWeight;
335 Float_t allNeutralWeight = pid[AliVCluster::kPhoton]+pid[AliVCluster::kPi0]+ neutralHadronWeight;
337 //Select most probable ID
340 if(pid[AliVCluster::kPhoton] > wPh) pdg = kPhoton ;
341 else if(pid[AliVCluster::kPi0] > wPi0) pdg = kPi0 ;
342 else if(pid[AliVCluster::kElectron] > wE) pdg = kElectron ;
343 else if(pid[AliVCluster::kEleCon] > wE) pdg = kEleCon ;
344 else if(chargedHadronWeight > wCh) pdg = kChargedHadron ;
345 else if(neutralHadronWeight > wNe) pdg = kNeutralHadron ;
346 else if(allChargedWeight > allNeutralWeight)
347 pdg = kChargedUnknown ;
349 pdg = kNeutralUnknown ;
353 if(pid[AliVCluster::kPhoton] > wPh) pdg = kPhoton ;
354 else if(pid[AliVCluster::kElectron] > wE) pdg = kElectron ;
355 else if(pid[AliVCluster::kPhoton]+pid[AliVCluster::kElectron] > wPh) pdg = kPhoton ; //temporal sollution until track matching for electrons is considered
356 else if(pid[AliVCluster::kPi0] > wPi0) pdg = kPi0 ;
357 else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;
358 else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ;
359 else pdg = kNeutralUnknown ;
362 if(fDebug > 0)printf("AliCaloPID::GetIdentifiedParticleType:Final Pdg: %d, cluster energy %2.2f \n", pdg,energy);
368 //____________________________________________________________________________________________________
369 Int_t AliCaloPID::GetIdentifiedParticleTypeFromClusterSplitting(AliVCluster* cluster,
370 AliVCaloCells* cells,
371 AliCalorimeterUtils * caloutils,
374 Double_t & mass, Double_t & angle,
375 Double_t & e1 , Double_t & e2 )
377 // Split the cluster in 2, do invariant mass, get the mass and decide
378 // if this is a photon, pi0, eta, ...
380 Int_t absId1 = -1; Int_t absId2 = -1;
381 Float_t l0 = cluster->GetM02();
382 const Int_t nc = cluster->GetNCells();
384 Float_t maxEList [nc];
390 //If too small or big E or low number of cells, or close to a bad channel skip it
391 if( l0 < fSplitM02MinCut || l0 > fSplitM02MaxCut || nc < fSplitMinNCells)
393 return kNeutralUnknown ;
396 // Get Number of local maxima
397 nMax = caloutils->GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ;
399 //---------------------------------------------------------------------
400 // Get the 2 max indeces and do inv mass
401 //---------------------------------------------------------------------
405 absId1 = absIdList[0];
406 absId2 = absIdList[1];
411 absId1 = absIdList[0];
413 //Find second highest energy cell
415 TString calorimeter = "EMCAL";
416 if(cluster->IsPHOS()) calorimeter = "PHOS";
418 for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++)
420 Int_t absId = cluster->GetCellsAbsId()[iDigit];
421 if( absId == absId1 ) continue ;
422 Float_t endig = cells->GetCellAmplitude(absId);
423 caloutils->RecalibrateCellAmplitude(endig,calorimeter,absId);
433 // loop on maxima, find 2 highest
437 for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++)
439 Float_t endig = maxEList[iDigit];
443 absId1 = absIdList[iDigit];
445 }// first maxima loop
449 for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++)
451 if(absIdList[iDigit]==absId1) continue;
452 Float_t endig = maxEList[iDigit];
456 absId2 = absIdList[iDigit];
458 }// second maxima loop
460 } // n local maxima > 2
462 //---------------------------------------------------------------------
463 // Split the cluster energy in 2, around the highest 2 local maxima
464 //---------------------------------------------------------------------
466 AliAODCaloCluster cluster1(0, 0,NULL,0.,NULL,NULL,1,0);
467 AliAODCaloCluster cluster2(1, 0,NULL,0.,NULL,NULL,1,0);
469 caloutils->SplitEnergy(absId1,absId2,cluster, cells, &cluster1, &cluster2,nMax); /*absIdList, maxEList,*/
471 TLorentzVector cellMom1;
472 TLorentzVector cellMom2;
474 cluster1.GetMomentum(cellMom1,vertex);
475 cluster2.GetMomentum(cellMom2,vertex);
477 mass = (cellMom1+cellMom2).M();
478 angle = cellMom2.Angle(cellMom1.Vect());
482 if (mass < fMassPhoMax && mass > fMassPhoMin) return kPhoton;
483 else if(mass < fMassPi0Max && mass > fMassPi0Min) return kPi0;
484 else if(mass < fMassEtaMax && mass > fMassEtaMin) return kEta;
485 else return kNeutralUnknown;
489 //_________________________________________
490 TString AliCaloPID::GetPIDParametersList()
492 //Put data member values in string to keep in output container
494 TString parList ; //this will be list of parameters used for this analysis.
495 const Int_t buffersize = 255;
496 char onePar[buffersize] ;
497 snprintf(onePar,buffersize,"--- AliCaloPID ---\n") ;
499 if(fUseBayesianWeights){
500 snprintf(onePar,buffersize,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)\n",fEMCALPhotonWeight) ;
502 snprintf(onePar,buffersize,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)\n",fEMCALPi0Weight) ;
504 snprintf(onePar,buffersize,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)\n",fEMCALElectronWeight) ;
506 snprintf(onePar,buffersize,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)\n",fEMCALChargeWeight) ;
508 snprintf(onePar,buffersize,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)\n",fEMCALNeutralWeight) ;
510 snprintf(onePar,buffersize,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)\n",fPHOSPhotonWeight) ;
512 snprintf(onePar,buffersize,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)\n",fPHOSPi0Weight) ;
514 snprintf(onePar,buffersize,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)\n",fPHOSElectronWeight) ;
516 snprintf(onePar,buffersize,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)\n",fPHOSChargeWeight) ;
518 snprintf(onePar,buffersize,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)\n",fPHOSNeutralWeight) ;
521 if(fPHOSWeightFormula){
522 snprintf(onePar,buffersize,"PHOS Photon Weight Formula: %s\n",fPHOSPhotonWeightFormulaExpression.Data() ) ;
524 snprintf(onePar,buffersize,"PHOS Pi0 Weight Formula: %s\n",fPHOSPi0WeightFormulaExpression.Data() ) ;
529 snprintf(onePar,buffersize,"EMCAL: fEMCALL0CutMin =%2.2f, fEMCALL0CutMax =%2.2f (Cut on Shower Shape) \n",fEMCALL0CutMin, fEMCALL0CutMax) ;
531 snprintf(onePar,buffersize,"EMCAL: fEMCALDEtaCut =%2.2f, fEMCALDPhiCut =%2.2f (Cut on track matching) \n",fEMCALDEtaCut, fEMCALDPhiCut) ;
533 snprintf(onePar,buffersize,"fTOFCut =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
535 snprintf(onePar,buffersize,"fPHOSRCut =%2.2f, fPHOSDispersionCut =%2.2f (Cut on Shower Shape and CPV) \n",fPHOSRCut,fPHOSDispersionCut) ;
540 if(fDoClusterSplitting)
542 snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n", fSplitM02MinCut, fSplitM02MaxCut) ;
544 snprintf(onePar,buffersize,"fMinNCells =%d \n", fSplitMinNCells) ;
546 snprintf(onePar,buffersize,"pi0 : %2.1f < m <%2.1f\n", fMassPi0Min,fMassPi0Max);
548 snprintf(onePar,buffersize,"eta : %2.1f < m <%2.1f\n", fMassEtaMin,fMassEtaMax);
550 snprintf(onePar,buffersize,"conv: %2.1f < m <%2.1f\n", fMassPhoMin,fMassPhoMax);
558 //________________________________________________
559 void AliCaloPID::Print(const Option_t * opt) const
562 //Print some relevant parameters set for the analysis
566 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
568 if(fUseBayesianWeights)
570 printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",
571 fPHOSPhotonWeight, fPHOSPi0Weight,
572 fPHOSElectronWeight, fPHOSChargeWeight, fPHOSNeutralWeight) ;
573 printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",
574 fEMCALPhotonWeight, fEMCALPi0Weight,
575 fEMCALElectronWeight, fEMCALChargeWeight, fEMCALNeutralWeight) ;
577 printf("PHOS Parametrized weight on? = %d\n", fPHOSWeightFormula) ;
578 if(fPHOSWeightFormula)
580 printf("Photon weight formula = %s\n", fPHOSPhotonWeightFormulaExpression.Data());
581 printf("Pi0 weight formula = %s\n", fPHOSPi0WeightFormulaExpression .Data());
583 if(fRecalculateBayesian) printf(" Recalculate bayesian with Particle Flux? = %d\n",fParticleFlux);
587 printf("TOF cut = %e\n", fTOFCut);
588 printf("EMCAL Lambda0 cut min = %2.2f; max = %2.2f\n", fEMCALL0CutMin,fEMCALL0CutMax);
589 printf("EMCAL cluster-track dEta < %2.3f; dPhi < %2.3f\n", fEMCALDEtaCut, fEMCALDPhiCut);
590 printf("PHOS Treac matching cut =%2.2f, Dispersion Cut =%2.2f \n",fPHOSRCut, fPHOSDispersionCut) ;
594 if(fDoClusterSplitting)
596 printf("Min. N Cells =%d \n", fSplitMinNCells) ;
597 printf("%2.2f < lambda_0^2 <%2.2f \n",fSplitM02MinCut,fSplitM02MaxCut);
598 printf("pi0 : %2.2f<m<%2.2f \n", fMassPi0Min,fMassPi0Max);
599 printf("eta : %2.2f<m<%2.2f \n", fMassEtaMin,fMassEtaMax);
600 printf("phot: %2.2f<m<%2.2f \n", fMassPhoMin,fMassPhoMax);
607 //_________________________________________________________________
608 void AliCaloPID::PrintClusterPIDWeights(const Double_t * pid) const
610 // print PID of cluster, (AliVCluster*)cluster->GetPID()
612 printf("AliCaloPID::PrintClusterPIDWeights() \n \t ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, \n \t \
613 pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
614 pid[AliVCluster::kPhoton], pid[AliVCluster::kPi0],
615 pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon],
616 pid[AliVCluster::kPion], pid[AliVCluster::kKaon],
617 pid[AliVCluster::kProton],
618 pid[AliVCluster::kNeutron], pid[AliVCluster::kKaon0]);
622 //___________________________________________________________________________
623 void AliCaloPID::SetPIDBits(AliVCluster * cluster,
624 AliAODPWG4Particle * ph, AliCalorimeterUtils* cu,
627 //Set Bits for PID selection
630 //Double_t disp= cluster->GetDispersion() ;
631 Double_t l1 = cluster->GetM20() ;
632 Double_t l0 = cluster->GetM02() ;
633 Bool_t isDispOK = kTRUE ;
634 if(cluster->IsPHOS()){
635 if(TestPHOSDispersion(ph->Pt(),l0,l1) < fPHOSDispersionCut) isDispOK = kTRUE;
636 else isDispOK = kFALSE;
640 if(l0 > fEMCALL0CutMin && l0 < fEMCALL0CutMax) isDispOK = kTRUE;
644 ph->SetDispBit(isDispOK) ;
647 Double_t tof=cluster->GetTOF() ;
648 ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
651 Bool_t isNeutral = IsTrackMatched(cluster,cu,event);
653 ph->SetChargedBit(isNeutral);
656 ph->SetIdentifiedParticleType(GetIdentifiedParticleType(cluster));
659 printf("AliCaloPID::SetPIDBits: TOF %e, Lambda0 %2.2f, Lambda1 %2.2f\n",tof , l0, l1);
660 printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
661 ph->GetIdentifiedParticleType(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit());
665 //_________________________________________________________
666 Bool_t AliCaloPID::IsTrackMatched(AliVCluster* cluster,
667 AliCalorimeterUtils * cu,
668 AliVEvent* event) const
670 //Check if there is any track attached to this cluster
672 Int_t nMatches = cluster->GetNTracksMatched();
673 AliVTrack * track = 0;
678 //In case of ESDs, by default without match one entry with negative index, no match, reject.
679 if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
681 Int_t iESDtrack = cluster->GetTrackMatchedIndex();
682 if(iESDtrack >= 0) track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
687 if(fDebug > 0) printf("AliCaloPID::IsTrackMatched() - Null matched track in ESD when index is OK!\n");
692 track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0));
695 if(fDebug > 0) printf("AliCaloPID::IsTrackMatched() - Null matched track in AOD!\n");
700 Float_t dZ = cluster->GetTrackDz();
701 Float_t dR = cluster->GetTrackDx();
703 // if track matching was recalculated
704 if(cluster->IsEMCAL() && cu && cu->IsRecalculationOfClusterTrackMatchingOn())
706 dR = 2000., dZ = 2000.;
707 cu->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
710 if(cluster->IsPHOS())
712 track->GetPxPyPz(p) ;
713 TLorentzVector trackmom(p[0],p[1],p[2],0);
714 Int_t charge = track->Charge();
715 Double_t mf = event->GetMagneticField();
716 if(TestPHOSChargedVeto(dR, dZ, trackmom.Pt(), charge, mf ) < fPHOSRCut) return kTRUE;
723 printf("AliCaloPID::IsTrackMatched - EMCAL dR %f < %f, dZ %f < %f \n",dR, fEMCALDPhiCut, dZ, fEMCALDEtaCut);
725 if(TMath::Abs(dR) < fEMCALDPhiCut &&
726 TMath::Abs(dZ) < fEMCALDEtaCut) return kTRUE;
732 } // more than 1 match, at least one track in array
737 //___________________________________________________________________________________________________
738 Float_t AliCaloPID::TestPHOSDispersion(const Double_t pt, const Double_t l1, const Double_t l2) const
740 //Check if cluster photon-like. Uses photon cluster parameterization in real pp data
741 //Returns distance in sigmas. Recommended cut 2.5
743 Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
744 Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
745 Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
746 Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
747 Double_t c =-0.35-0.550*TMath::Exp(-0.390730*pt) ;
748 Double_t r2 = 0.5* (l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
749 0.5* (l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
750 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
752 if(fDebug > 0) printf("AliCaloPID::TestPHOSDispersion() - PHOS SS R %f < %f?\n", TMath::Sqrt(r2), fPHOSDispersionCut);
754 return TMath::Sqrt(r2) ;
758 //_______________________________________________________________________________________________
759 Float_t AliCaloPID::TestPHOSChargedVeto(const Double_t dx, const Double_t dz, const Double_t pt,
760 const Int_t charge, const Double_t mf) const
762 //Checks distance to the closest track. Takes into account
763 //non-perpendicular incidence of tracks.
764 //returns distance in sigmas. Recommended cut: 2.
765 //Requires (sign) of magnetic filed. onc can find it for example as following
767 // AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
769 // mf = event->GetMagneticField(); //Positive for ++ and negative for --
774 Double_t sx = TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
775 6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+
777 Double_t sz = TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+
780 if(mf<0.){ //field --
783 meanX = TMath::Min(7.3, 3.89994*1.20679 *1.20679 /(pt*pt+1.20679*1.20679)+
784 0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
786 meanX =-TMath::Min(7.7, 3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+
787 1.23114 +4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
792 meanX =-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+
793 0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
795 meanX = TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-
796 0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
799 Double_t rz = (dz-meanZ)/sz ;
800 Double_t rx = (dx-meanX)/sx ;
803 printf("AliCaloPID::TestPHOSDispersion() - PHOS Matching R %f < %f\n",TMath::Sqrt(rx*rx+rz*rz), fPHOSRCut);
805 return TMath::Sqrt(rx*rx+rz*rz) ;