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)
376 // Split the cluster in 2, do invariant mass, get the mass and decide
377 // if this is a photon, pi0, eta, ...
379 Int_t absId1 = -1; Int_t absId2 = -1;
380 Int_t nc = cluster->GetNCells();
381 Int_t *absIdList = new Int_t [nc];
382 Float_t *maxEList = new Float_t[nc];
383 Float_t l0 = cluster->GetM02();
389 //If too small or big E or low number of cells, or close to a bad channel skip it
390 if( l0 < fSplitM02MinCut || l0 > fSplitM02MaxCut || nc < fSplitMinNCells) return kNeutralUnknown ;
392 // Get Number of local maxima
393 nMax = caloutils->GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ;
395 //---------------------------------------------------------------------
396 // Get the 2 max indeces and do inv mass
397 //---------------------------------------------------------------------
401 absId1 = absIdList[0];
402 absId2 = absIdList[1];
407 absId1 = absIdList[0];
409 //Find second highest energy cell
411 TString calorimeter = "EMCAL";
412 if(cluster->IsPHOS()) calorimeter = "PHOS";
414 for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++)
416 Int_t absId = cluster->GetCellsAbsId()[iDigit];
417 if( absId == absId1 ) continue ;
418 Float_t endig = cells->GetCellAmplitude(absId);
419 caloutils->RecalibrateCellAmplitude(endig,calorimeter,absId);
429 // loop on maxima, find 2 highest
433 for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++)
435 Float_t endig = maxEList[iDigit];
439 absId1 = absIdList[iDigit];
441 }// first maxima loop
445 for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++)
447 if(absIdList[iDigit]==absId1) continue;
448 Float_t endig = maxEList[iDigit];
452 absId2 = absIdList[iDigit];
454 }// second maxima loop
456 } // n local maxima > 2
458 //---------------------------------------------------------------------
459 // Split the cluster energy in 2, around the highest 2 local maxima
460 //---------------------------------------------------------------------
462 AliAODCaloCluster *cluster1 = new AliAODCaloCluster(0, 0,NULL,0.,NULL,NULL,1,0);
463 AliAODCaloCluster *cluster2 = new AliAODCaloCluster(1, 0,NULL,0.,NULL,NULL,1,0);
465 caloutils->SplitEnergy(absId1,absId2,cluster, cells, cluster1, cluster2,nMax); /*absIdList, maxEList,*/
467 TLorentzVector cellMom1;
468 TLorentzVector cellMom2;
470 cluster1->GetMomentum(cellMom1,vertex);
471 cluster2->GetMomentum(cellMom2,vertex);
473 mass = (cellMom1+cellMom2).M();
474 angle = cellMom2.Angle(cellMom1.Vect());
476 if (mass < fMassPhoMax && mass > fMassPhoMin) return kPhoton;
477 else if(mass < fMassPi0Max && mass > fMassPi0Min) return kPi0;
478 else if(mass < fMassEtaMax && mass > fMassEtaMin) return kEta;
479 else return kNeutralUnknown;
483 delete [] absIdList ;
488 //_________________________________________
489 TString AliCaloPID::GetPIDParametersList()
491 //Put data member values in string to keep in output container
493 TString parList ; //this will be list of parameters used for this analysis.
494 const Int_t buffersize = 255;
495 char onePar[buffersize] ;
496 snprintf(onePar,buffersize,"--- AliCaloPID ---\n") ;
498 if(fUseBayesianWeights){
499 snprintf(onePar,buffersize,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)\n",fEMCALPhotonWeight) ;
501 snprintf(onePar,buffersize,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)\n",fEMCALPi0Weight) ;
503 snprintf(onePar,buffersize,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)\n",fEMCALElectronWeight) ;
505 snprintf(onePar,buffersize,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)\n",fEMCALChargeWeight) ;
507 snprintf(onePar,buffersize,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)\n",fEMCALNeutralWeight) ;
509 snprintf(onePar,buffersize,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)\n",fPHOSPhotonWeight) ;
511 snprintf(onePar,buffersize,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)\n",fPHOSPi0Weight) ;
513 snprintf(onePar,buffersize,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)\n",fPHOSElectronWeight) ;
515 snprintf(onePar,buffersize,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)\n",fPHOSChargeWeight) ;
517 snprintf(onePar,buffersize,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)\n",fPHOSNeutralWeight) ;
520 if(fPHOSWeightFormula){
521 snprintf(onePar,buffersize,"PHOS Photon Weight Formula: %s\n",fPHOSPhotonWeightFormulaExpression.Data() ) ;
523 snprintf(onePar,buffersize,"PHOS Pi0 Weight Formula: %s\n",fPHOSPi0WeightFormulaExpression.Data() ) ;
528 snprintf(onePar,buffersize,"EMCAL: fEMCALL0CutMin =%2.2f, fEMCALL0CutMax =%2.2f (Cut on Shower Shape) \n",fEMCALL0CutMin, fEMCALL0CutMax) ;
530 snprintf(onePar,buffersize,"EMCAL: fEMCALDEtaCut =%2.2f, fEMCALDPhiCut =%2.2f (Cut on track matching) \n",fEMCALDEtaCut, fEMCALDPhiCut) ;
532 snprintf(onePar,buffersize,"fTOFCut =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
534 snprintf(onePar,buffersize,"fPHOSRCut =%2.2f, fPHOSDispersionCut =%2.2f (Cut on Shower Shape and CPV) \n",fPHOSRCut,fPHOSDispersionCut) ;
539 if(fDoClusterSplitting)
541 snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n", fSplitM02MinCut, fSplitM02MaxCut) ;
543 snprintf(onePar,buffersize,"fMinNCells =%d \n", fSplitMinNCells) ;
545 snprintf(onePar,buffersize,"pi0 : %2.1f < m <%2.1f\n", fMassPi0Min,fMassPi0Max);
547 snprintf(onePar,buffersize,"eta : %2.1f < m <%2.1f\n", fMassEtaMin,fMassEtaMax);
549 snprintf(onePar,buffersize,"conv: %2.1f < m <%2.1f\n", fMassPhoMin,fMassPhoMax);
557 //________________________________________________
558 void AliCaloPID::Print(const Option_t * opt) const
561 //Print some relevant parameters set for the analysis
565 printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
567 if(fUseBayesianWeights)
569 printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",
570 fPHOSPhotonWeight, fPHOSPi0Weight,
571 fPHOSElectronWeight, fPHOSChargeWeight, fPHOSNeutralWeight) ;
572 printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",
573 fEMCALPhotonWeight, fEMCALPi0Weight,
574 fEMCALElectronWeight, fEMCALChargeWeight, fEMCALNeutralWeight) ;
576 printf("PHOS Parametrized weight on? = %d\n", fPHOSWeightFormula) ;
577 if(fPHOSWeightFormula)
579 printf("Photon weight formula = %s\n", fPHOSPhotonWeightFormulaExpression.Data());
580 printf("Pi0 weight formula = %s\n", fPHOSPi0WeightFormulaExpression .Data());
582 if(fRecalculateBayesian) printf(" Recalculate bayesian with Particle Flux? = %d\n",fParticleFlux);
586 printf("TOF cut = %e\n", fTOFCut);
587 printf("EMCAL Lambda0 cut min = %2.2f; max = %2.2f\n", fEMCALL0CutMin,fEMCALL0CutMax);
588 printf("EMCAL cluster-track dEta < %2.3f; dPhi < %2.3f\n", fEMCALDEtaCut, fEMCALDPhiCut);
589 printf("PHOS Treac matching cut =%2.2f, Dispersion Cut =%2.2f \n",fPHOSRCut, fPHOSDispersionCut) ;
593 if(fDoClusterSplitting)
595 printf("Min. N Cells =%d \n", fSplitMinNCells) ;
596 printf("%2.2f < lambda_0^2 <%2.2f \n",fSplitM02MinCut,fSplitM02MaxCut);
597 printf("pi0 : %2.2f<m<%2.2f \n", fMassPi0Min,fMassPi0Max);
598 printf("eta : %2.2f<m<%2.2f \n", fMassEtaMin,fMassEtaMax);
599 printf("phot: %2.2f<m<%2.2f \n", fMassPhoMin,fMassPhoMax);
606 //_________________________________________________________________
607 void AliCaloPID::PrintClusterPIDWeights(const Double_t * pid) const
609 // print PID of cluster, (AliVCluster*)cluster->GetPID()
611 printf("AliCaloPID::PrintClusterPIDWeights() \n \t ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, \n \t \
612 pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
613 pid[AliVCluster::kPhoton], pid[AliVCluster::kPi0],
614 pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon],
615 pid[AliVCluster::kPion], pid[AliVCluster::kKaon],
616 pid[AliVCluster::kProton],
617 pid[AliVCluster::kNeutron], pid[AliVCluster::kKaon0]);
621 //___________________________________________________________________________
622 void AliCaloPID::SetPIDBits(AliVCluster * cluster,
623 AliAODPWG4Particle * ph, AliCalorimeterUtils* cu,
626 //Set Bits for PID selection
629 //Double_t disp= cluster->GetDispersion() ;
630 Double_t l1 = cluster->GetM20() ;
631 Double_t l0 = cluster->GetM02() ;
632 Bool_t isDispOK = kTRUE ;
633 if(cluster->IsPHOS()){
634 if(TestPHOSDispersion(ph->Pt(),l0,l1) < fPHOSDispersionCut) isDispOK = kTRUE;
635 else isDispOK = kFALSE;
639 if(l0 > fEMCALL0CutMin && l0 < fEMCALL0CutMax) isDispOK = kTRUE;
643 ph->SetDispBit(isDispOK) ;
646 Double_t tof=cluster->GetTOF() ;
647 ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
650 Bool_t isNeutral = IsTrackMatched(cluster,cu,event);
652 ph->SetChargedBit(isNeutral);
655 ph->SetIdentifiedParticleType(GetIdentifiedParticleType(cluster));
658 printf("AliCaloPID::SetPIDBits: TOF %e, Lambda0 %2.2f, Lambda1 %2.2f\n",tof , l0, l1);
659 printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
660 ph->GetIdentifiedParticleType(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit());
664 //_________________________________________________________
665 Bool_t AliCaloPID::IsTrackMatched(AliVCluster* cluster,
666 AliCalorimeterUtils * cu,
667 AliVEvent* event) const
669 //Check if there is any track attached to this cluster
671 Int_t nMatches = cluster->GetNTracksMatched();
672 AliVTrack * track = 0;
677 //In case of ESDs, by default without match one entry with negative index, no match, reject.
678 if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
680 Int_t iESDtrack = cluster->GetTrackMatchedIndex();
681 if(iESDtrack >= 0) track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
685 printf("AliCaloPID::IsTrackMatched() - Null matched track in ESD when index is OK!\n");
690 track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0));
692 printf("AliCaloPID::IsTrackMatched() - Null matched track in AOD!\n");
697 Float_t dZ = cluster->GetTrackDz();
698 Float_t dR = cluster->GetTrackDx();
700 // if track matching was recalculated
701 if(cluster->IsEMCAL() && cu && cu->IsRecalculationOfClusterTrackMatchingOn()){
702 dR = 2000., dZ = 2000.;
703 cu->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
706 if(cluster->IsPHOS()) {
708 track->GetPxPyPz(p) ;
709 TLorentzVector trackmom(p[0],p[1],p[2],0);
710 Int_t charge = track->Charge();
711 Double_t mf = event->GetMagneticField();
712 if(TestPHOSChargedVeto(dR, dZ, trackmom.Pt(), charge, mf ) < fPHOSRCut) return kTRUE;
719 printf("AliCaloPID::IsTrackMatched - EMCAL dR %f < %f, dZ %f < %f \n",dR, fEMCALDPhiCut, dZ, fEMCALDEtaCut);
721 if(TMath::Abs(dR) < fEMCALDPhiCut &&
722 TMath::Abs(dZ) < fEMCALDEtaCut) return kTRUE;
728 } // more than 1 match, at least one track in array
733 //___________________________________________________________________________________________________
734 Float_t AliCaloPID::TestPHOSDispersion(const Double_t pt, const Double_t l1, const Double_t l2) const
736 //Check if cluster photon-like. Uses photon cluster parameterization in real pp data
737 //Returns distance in sigmas. Recommended cut 2.5
739 Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
740 Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
741 Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
742 Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
743 Double_t c =-0.35-0.550*TMath::Exp(-0.390730*pt) ;
744 Double_t r2 = 0.5* (l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
745 0.5* (l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
746 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
748 if(fDebug > 0) printf("AliCaloPID::TestPHOSDispersion() - PHOS SS R %f < %f?\n", TMath::Sqrt(r2), fPHOSDispersionCut);
750 return TMath::Sqrt(r2) ;
754 //_______________________________________________________________________________________________
755 Float_t AliCaloPID::TestPHOSChargedVeto(const Double_t dx, const Double_t dz, const Double_t pt,
756 const Int_t charge, const Double_t mf) const
758 //Checks distance to the closest track. Takes into account
759 //non-perpendicular incidence of tracks.
760 //returns distance in sigmas. Recommended cut: 2.
761 //Requires (sign) of magnetic filed. onc can find it for example as following
763 // AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
765 // mf = event->GetMagneticField(); //Positive for ++ and negative for --
770 Double_t sx = TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
771 6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+
773 Double_t sz = TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+
776 if(mf<0.){ //field --
779 meanX = TMath::Min(7.3, 3.89994*1.20679 *1.20679 /(pt*pt+1.20679*1.20679)+
780 0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
782 meanX =-TMath::Min(7.7, 3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+
783 1.23114 +4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
788 meanX =-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+
789 0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
791 meanX = TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-
792 0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
795 Double_t rz = (dz-meanZ)/sz ;
796 Double_t rx = (dx-meanX)/sx ;
799 printf("AliCaloPID::TestPHOSDispersion() - PHOS Matching R %f < %f\n",TMath::Sqrt(rx*rx+rz*rz), fPHOSRCut);
801 return TMath::Sqrt(rx*rx+rz*rz) ;