]>
Commit | Line | Data |
---|---|---|
1c5acb87 | 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 | **************************************************************************/ | |
1c5acb87 | 15 | |
16 | //_________________________________________________________________________ | |
bdd2a262 | 17 | // Class for PID selection with calorimeters |
49b5c49b | 18 | // The Output of the main method GetIdentifiedParticleType is a PDG number identifying the cluster, |
bdd2a262 | 19 | // being kPhoton, kElectron, kPi0 ... as defined in the header file |
3c1d9afb | 20 | // - GetIdentifiedParticleType(const AliVCluster * cluster) |
49b5c49b | 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. | |
bdd2a262 | 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 | |
26 | // AliCaloPID(flux) | |
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. | |
49b5c49b | 30 | |
3c1d9afb | 31 | // - GetGetIdentifiedParticleTypeFromBayesian(const Double_t * pid, const Float_t energy) |
49b5c49b | 32 | // Reads the PID weights array of the ESDs and depending on its magnitude identifies the particle, |
3c1d9afb | 33 | // executed when bayesian is ON by GetIdentifiedParticleType(const AliVCluster * cluster) |
49b5c49b | 34 | // - SetPIDBits: Simple PID, depending on the thresholds fLOCut fTOFCut and even the |
bdd2a262 | 35 | // result of the PID bayesian a different PID bit is set. |
36 | // | |
49b5c49b | 37 | // - IsTrackMatched(): Independent method that needs to be combined with GetIdentifiedParticleType to know if the cluster was matched |
bdd2a262 | 38 | // |
49b5c49b | 39 | //*-- Author: Gustavo Conesa (INFN-LNF) |
1c5acb87 | 40 | ////////////////////////////////////////////////////////////////////////////// |
41 | ||
42 | ||
43 | // --- ROOT system --- | |
44 | #include <TMath.h> | |
1c5acb87 | 45 | #include <TString.h> |
f21fc003 | 46 | #include <TList.h> |
1c5acb87 | 47 | |
c5693f62 | 48 | // ---- ANALYSIS system ---- |
1c5acb87 | 49 | #include "AliCaloPID.h" |
3c1d9afb | 50 | #include "AliAODCaloCluster.h" |
51 | #include "AliVCaloCells.h" | |
d39cba7e | 52 | #include "AliVTrack.h" |
1c5acb87 | 53 | #include "AliAODPWG4Particle.h" |
f2ccb5b8 | 54 | #include "AliCalorimeterUtils.h" |
49b5c49b | 55 | #include "AliVEvent.h" |
f2ccb5b8 | 56 | |
c5693f62 | 57 | // ---- Detector ---- |
58 | #include "AliEMCALPIDUtils.h" | |
59 | ||
1c5acb87 | 60 | ClassImp(AliCaloPID) |
61 | ||
62 | ||
49b5c49b | 63 | //________________________ |
1c5acb87 | 64 | AliCaloPID::AliCaloPID() : |
49b5c49b | 65 | TObject(), fDebug(-1), fParticleFlux(kLow), |
66 | //Bayesian | |
67 | fEMCALPIDUtils(), fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE), | |
a5fb4114 | 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), | |
49b5c49b | 73 | fPHOSPhotonWeightFormulaExpression(""), |
74 | fPHOSPi0WeightFormulaExpression(""), | |
75 | //PID calculation | |
76 | fEMCALL0CutMax(100.), fEMCALL0CutMin(0), | |
77 | fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.), | |
78 | fTOFCut(0.), | |
3c1d9afb | 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) | |
1c5acb87 | 85 | { |
477d6cee | 86 | //Ctor |
87 | ||
88 | //Initialize parameters | |
89 | InitParameters(); | |
1c5acb87 | 90 | } |
91 | ||
49b5c49b | 92 | //________________________________________ |
bdd2a262 | 93 | AliCaloPID::AliCaloPID(const Int_t flux) : |
49b5c49b | 94 | TObject(), fDebug(-1), fParticleFlux(flux), |
95 | //Bayesian | |
96 | fEMCALPIDUtils(), fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE), | |
a5fb4114 | 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), | |
49b5c49b | 102 | fPHOSPhotonWeightFormulaExpression(""), |
103 | fPHOSPi0WeightFormulaExpression(""), | |
104 | //PID calculation | |
105 | fEMCALL0CutMax(100.), fEMCALL0CutMin(0), | |
106 | fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.), | |
107 | fTOFCut(0.), | |
3c1d9afb | 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) | |
bdd2a262 | 114 | { |
9a6fa057 | 115 | //Ctor |
bdd2a262 | 116 | |
9a6fa057 | 117 | //Initialize parameters |
118 | InitParameters(); | |
49b5c49b | 119 | |
bdd2a262 | 120 | } |
121 | ||
49b5c49b | 122 | //_______________________________________________ |
f21fc003 | 123 | AliCaloPID::AliCaloPID(const TNamed * emcalpid) : |
49b5c49b | 124 | TObject(), fDebug(-1), fParticleFlux(kLow), |
125 | //Bayesian | |
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(""), | |
135 | //PID calculation | |
136 | fEMCALL0CutMax(100.), fEMCALL0CutMin(0), | |
137 | fEMCALDEtaCut(2000.), fEMCALDPhiCut(2000.), | |
138 | fTOFCut(0.), | |
3c1d9afb | 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) | |
bdd2a262 | 145 | { |
9a6fa057 | 146 | //Ctor |
49b5c49b | 147 | |
9a6fa057 | 148 | //Initialize parameters |
149 | InitParameters(); | |
bdd2a262 | 150 | } |
151 | ||
49b5c49b | 152 | //_______________________ |
153 | AliCaloPID::~AliCaloPID() | |
154 | { | |
477d6cee | 155 | //Dtor |
156 | ||
a5fb4114 | 157 | delete fPHOSPhotonWeightFormula ; |
158 | delete fPHOSPi0WeightFormula ; | |
159 | delete fEMCALPIDUtils ; | |
49b5c49b | 160 | |
9a6fa057 | 161 | } |
1c5acb87 | 162 | |
49b5c49b | 163 | //_______________________________ |
1c5acb87 | 164 | void AliCaloPID::InitParameters() |
165 | { | |
477d6cee | 166 | //Initialize the parameters of the PID. |
167 | ||
49b5c49b | 168 | // Bayesian |
2007809d | 169 | fEMCALPhotonWeight = 0.6 ; |
170 | fEMCALPi0Weight = 0.6 ; | |
171 | fEMCALElectronWeight = 0.6 ; | |
172 | fEMCALChargeWeight = 0.6 ; | |
173 | fEMCALNeutralWeight = 0.6 ; | |
477d6cee | 174 | |
2007809d | 175 | fPHOSPhotonWeight = 0.6 ; |
176 | fPHOSPi0Weight = 0.6 ; | |
177 | fPHOSElectronWeight = 0.6 ; | |
178 | fPHOSChargeWeight = 0.6 ; | |
179 | fPHOSNeutralWeight = 0.6 ; | |
477d6cee | 180 | |
181 | //Formula to set the PID weight threshold for photon or pi0 | |
a5fb4114 | 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))" ; | |
185 | ||
bdd2a262 | 186 | if(fRecalculateBayesian){ |
49b5c49b | 187 | if(fParticleFlux == kLow){ |
188 | printf("AliCaloPID::Init() - SetLOWFluxParam\n"); | |
189 | fEMCALPIDUtils->SetLowFluxParam() ; | |
190 | } | |
191 | else if (fParticleFlux == kHigh){ | |
192 | printf("AliCaloPID::Init() - SetHIGHFluxParam\n"); | |
193 | fEMCALPIDUtils->SetHighFluxParam() ; | |
194 | } | |
195 | } | |
196 | ||
197 | //PID recalculation, not bayesian | |
198 | ||
199 | //EMCAL | |
200 | fEMCALL0CutMax = 0.3 ; | |
201 | fEMCALL0CutMin = 0.01; | |
202 | ||
203 | fEMCALDPhiCut = 0.05; // Same cut as in AliEMCALRecoUtils | |
204 | fEMCALDEtaCut = 0.025;// Same cut as in AliEMCALRecoUtils | |
205 | ||
206 | // PHOS / EMCAL, not used | |
207 | fTOFCut = 1.e-6; | |
208 | ||
209 | //PHOS | |
210 | fPHOSRCut = 2. ; | |
211 | fPHOSDispersionCut = 2.5; | |
212 | ||
3c1d9afb | 213 | // Cluster splitting |
214 | ||
215 | fSplitM02MinCut = 0.26 ; | |
216 | fSplitM02MaxCut = 100 ; | |
217 | fSplitMinNCells = 4 ; | |
218 | ||
219 | fMassEtaMin = 0.4; | |
220 | fMassEtaMax = 0.6; | |
221 | ||
222 | fMassPi0Min = 0.08; | |
223 | fMassPi0Max = 0.20; | |
224 | ||
225 | fMassPhoMin = 0.0; | |
226 | fMassPhoMax = 0.05; | |
227 | ||
49b5c49b | 228 | } |
229 | ||
c5693f62 | 230 | //______________________________________________ |
231 | AliEMCALPIDUtils *AliCaloPID::GetEMCALPIDUtils() | |
232 | { | |
233 | // return pointer to AliEMCALPIDUtils, create it if needed | |
234 | ||
235 | if(!fEMCALPIDUtils) fEMCALPIDUtils = new AliEMCALPIDUtils ; | |
236 | return fEMCALPIDUtils ; | |
237 | ||
238 | } | |
239 | ||
240 | ||
49b5c49b | 241 | //______________________________________________________________________ |
3c1d9afb | 242 | Int_t AliCaloPID::GetIdentifiedParticleType(const AliVCluster * cluster) |
49b5c49b | 243 | { |
244 | // Returns a PDG number corresponding to the likely ID of the cluster | |
245 | ||
3c1d9afb | 246 | Float_t energy = cluster->E(); |
49b5c49b | 247 | Float_t lambda0 = cluster->GetM02(); |
248 | Float_t lambda1 = cluster->GetM20(); | |
249 | ||
250 | // --------------------- | |
251 | // Use bayesian approach | |
252 | // --------------------- | |
253 | ||
3c1d9afb | 254 | if(fUseBayesianWeights) |
255 | { | |
49b5c49b | 256 | Double_t weights[AliPID::kSPECIESN]; |
257 | ||
3c1d9afb | 258 | if(cluster->IsEMCAL() && fRecalculateBayesian) |
259 | { | |
49b5c49b | 260 | fEMCALPIDUtils->ComputePID(energy, lambda0); |
261 | for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = fEMCALPIDUtils->GetPIDFinal(i); | |
262 | } | |
3c1d9afb | 263 | else |
264 | { | |
49b5c49b | 265 | for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = cluster->GetPID()[i]; |
266 | } | |
267 | ||
3c1d9afb | 268 | if(fDebug > 0) PrintClusterPIDWeights(weights); |
49b5c49b | 269 | |
3c1d9afb | 270 | return GetIdentifiedParticleTypeFromBayesWeights(cluster->IsEMCAL(), weights, energy); |
bdd2a262 | 271 | } |
49b5c49b | 272 | |
273 | // ------------------------------------------------------- | |
274 | // Calculate PID SS from data, do not use bayesian weights | |
275 | // ------------------------------------------------------- | |
276 | ||
3c1d9afb | 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(), | |
49b5c49b | 279 | cluster->GetEmcCpvDistance(), cluster->GetDistanceToBadChannel(),cluster->GetNExMax()); |
280 | ||
3c1d9afb | 281 | if(cluster->IsEMCAL()) |
282 | { | |
49b5c49b | 283 | if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType() - EMCAL SS %f <%f < %f?\n",fEMCALL0CutMin, lambda0, fEMCALL0CutMax); |
284 | ||
285 | if(lambda0 < fEMCALL0CutMax && lambda0 > fEMCALL0CutMin) return kPhoton ; | |
286 | else return kNeutralUnknown ; | |
3c1d9afb | 287 | } // EMCAL |
288 | else // PHOS | |
289 | { | |
290 | if(TestPHOSDispersion(energy,lambda0,lambda1) < fPHOSDispersionCut) return kPhoton; | |
291 | else return kNeutralUnknown; | |
49b5c49b | 292 | } |
293 | ||
1c5acb87 | 294 | } |
295 | ||
49b5c49b | 296 | //_______________________________________________________________________________ |
3c1d9afb | 297 | Int_t AliCaloPID::GetIdentifiedParticleTypeFromBayesWeights(const Bool_t isEMCAL, |
49b5c49b | 298 | const Double_t * pid, |
299 | const Float_t energy) | |
300 | { | |
301 | //Return most probable identity of the particle after bayesian weights calculated in reconstruction | |
477d6cee | 302 | |
3c1d9afb | 303 | if(!pid) |
304 | { | |
21a4b1c0 | 305 | printf("AliCaloPID::GetIdentifiedParticleType() - pid pointer not initialized!!!\n"); |
477d6cee | 306 | abort(); |
307 | } | |
308 | ||
15800db4 | 309 | Float_t wPh = fPHOSPhotonWeight ; |
477d6cee | 310 | Float_t wPi0 = fPHOSPi0Weight ; |
15800db4 | 311 | Float_t wE = fPHOSElectronWeight ; |
312 | Float_t wCh = fPHOSChargeWeight ; | |
313 | Float_t wNe = fPHOSNeutralWeight ; | |
49b5c49b | 314 | |
3c1d9afb | 315 | if(!isEMCAL && fPHOSWeightFormula){ |
a5fb4114 | 316 | wPh = GetPHOSPhotonWeightFormula()->Eval(energy) ; |
317 | wPi0 = GetPHOSPi0WeightFormula() ->Eval(energy); | |
318 | } | |
3c1d9afb | 319 | else |
320 | { | |
477d6cee | 321 | wPh = fEMCALPhotonWeight ; |
322 | wPi0 = fEMCALPi0Weight ; | |
323 | wE = fEMCALElectronWeight ; | |
324 | wCh = fEMCALChargeWeight ; | |
325 | wNe = fEMCALNeutralWeight ; | |
477d6cee | 326 | } |
327 | ||
3c1d9afb | 328 | if(fDebug > 0) PrintClusterPIDWeights(pid); |
329 | ||
477d6cee | 330 | Int_t pdg = kNeutralUnknown ; |
c8fe2783 | 331 | Float_t chargedHadronWeight = pid[AliVCluster::kProton]+pid[AliVCluster::kKaon]+ |
49b5c49b | 332 | pid[AliVCluster::kPion]+pid[AliVCluster::kMuon]; |
c8fe2783 | 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; | |
477d6cee | 336 | |
337 | //Select most probable ID | |
3c1d9afb | 338 | if(!isEMCAL) // PHOS |
339 | { | |
a5fb4114 | 340 | if(pid[AliVCluster::kPhoton] > wPh) pdg = kPhoton ; |
341 | else if(pid[AliVCluster::kPi0] > wPi0) pdg = kPi0 ; | |
c8fe2783 | 342 | else if(pid[AliVCluster::kElectron] > wE) pdg = kElectron ; |
a5fb4114 | 343 | else if(pid[AliVCluster::kEleCon] > wE) pdg = kEleCon ; |
344 | else if(chargedHadronWeight > wCh) pdg = kChargedHadron ; | |
345 | else if(neutralHadronWeight > wNe) pdg = kNeutralHadron ; | |
477d6cee | 346 | else if(allChargedWeight > allNeutralWeight) |
347 | pdg = kChargedUnknown ; | |
348 | else | |
349 | pdg = kNeutralUnknown ; | |
350 | } | |
3c1d9afb | 351 | else //EMCAL |
352 | { | |
2007809d | 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 ; | |
477d6cee | 357 | else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ; |
358 | else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ; | |
2007809d | 359 | else pdg = kNeutralUnknown ; |
477d6cee | 360 | } |
361 | ||
21a4b1c0 | 362 | if(fDebug > 0)printf("AliCaloPID::GetIdentifiedParticleType:Final Pdg: %d, cluster energy %2.2f \n", pdg,energy); |
1c5acb87 | 363 | |
49b5c49b | 364 | return pdg ; |
9a6fa057 | 365 | |
1c5acb87 | 366 | } |
367 | ||
3c1d9afb | 368 | //____________________________________________________________________________________________________ |
369 | Int_t AliCaloPID::GetIdentifiedParticleTypeFromClusterSplitting(AliVCluster* cluster, | |
370 | AliVCaloCells* cells, | |
371 | AliCalorimeterUtils * caloutils, | |
372 | Double_t vertex[3], | |
373 | Int_t & nMax, | |
374 | Double_t & mass, Double_t & angle) | |
375 | { | |
376 | // Split the cluster in 2, do invariant mass, get the mass and decide | |
377 | // if this is a photon, pi0, eta, ... | |
378 | ||
2bf17171 | 379 | Int_t absId1 = -1; Int_t absId2 = -1; |
380 | Float_t l0 = cluster->GetM02(); | |
381 | const Int_t nc = cluster->GetNCells(); | |
382 | Int_t absIdList[nc]; | |
383 | Float_t maxEList [nc]; | |
3c1d9afb | 384 | |
385 | nMax = -1 ; | |
386 | mass = -1.; | |
387 | angle = -1.; | |
388 | ||
389 | //If too small or big E or low number of cells, or close to a bad channel skip it | |
55d66f31 | 390 | if( l0 < fSplitM02MinCut || l0 > fSplitM02MaxCut || nc < fSplitMinNCells) |
391 | { | |
392 | return kNeutralUnknown ; | |
393 | } | |
3c1d9afb | 394 | |
395 | // Get Number of local maxima | |
396 | nMax = caloutils->GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ; | |
397 | ||
398 | //--------------------------------------------------------------------- | |
399 | // Get the 2 max indeces and do inv mass | |
400 | //--------------------------------------------------------------------- | |
401 | ||
402 | if ( nMax == 2 ) | |
403 | { | |
404 | absId1 = absIdList[0]; | |
405 | absId2 = absIdList[1]; | |
406 | } | |
407 | else if( nMax == 1 ) | |
408 | { | |
409 | ||
410 | absId1 = absIdList[0]; | |
411 | ||
412 | //Find second highest energy cell | |
413 | ||
414 | TString calorimeter = "EMCAL"; | |
415 | if(cluster->IsPHOS()) calorimeter = "PHOS"; | |
416 | Float_t enmax = 0 ; | |
417 | for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++) | |
418 | { | |
419 | Int_t absId = cluster->GetCellsAbsId()[iDigit]; | |
420 | if( absId == absId1 ) continue ; | |
421 | Float_t endig = cells->GetCellAmplitude(absId); | |
422 | caloutils->RecalibrateCellAmplitude(endig,calorimeter,absId); | |
423 | if(endig > enmax) | |
424 | { | |
425 | enmax = endig ; | |
426 | absId2 = absId ; | |
427 | } | |
428 | }// cell loop | |
429 | }// 1 maxima | |
430 | else | |
431 | { // n max > 2 | |
432 | // loop on maxima, find 2 highest | |
433 | ||
434 | // First max | |
435 | Float_t enmax = 0 ; | |
436 | for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++) | |
437 | { | |
438 | Float_t endig = maxEList[iDigit]; | |
439 | if(endig > enmax) | |
440 | { | |
441 | enmax = endig ; | |
442 | absId1 = absIdList[iDigit]; | |
443 | } | |
444 | }// first maxima loop | |
445 | ||
446 | // Second max | |
447 | Float_t enmax2 = 0; | |
448 | for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++) | |
449 | { | |
450 | if(absIdList[iDigit]==absId1) continue; | |
451 | Float_t endig = maxEList[iDigit]; | |
452 | if(endig > enmax2) | |
453 | { | |
454 | enmax2 = endig ; | |
455 | absId2 = absIdList[iDigit]; | |
456 | } | |
457 | }// second maxima loop | |
458 | ||
459 | } // n local maxima > 2 | |
460 | ||
461 | //--------------------------------------------------------------------- | |
462 | // Split the cluster energy in 2, around the highest 2 local maxima | |
463 | //--------------------------------------------------------------------- | |
464 | ||
2bf17171 | 465 | AliAODCaloCluster cluster1(0, 0,NULL,0.,NULL,NULL,1,0); |
466 | AliAODCaloCluster cluster2(1, 0,NULL,0.,NULL,NULL,1,0); | |
3c1d9afb | 467 | |
2bf17171 | 468 | caloutils->SplitEnergy(absId1,absId2,cluster, cells, &cluster1, &cluster2,nMax); /*absIdList, maxEList,*/ |
3c1d9afb | 469 | |
470 | TLorentzVector cellMom1; | |
471 | TLorentzVector cellMom2; | |
472 | ||
2bf17171 | 473 | cluster1.GetMomentum(cellMom1,vertex); |
474 | cluster2.GetMomentum(cellMom2,vertex); | |
3c1d9afb | 475 | |
476 | mass = (cellMom1+cellMom2).M(); | |
477 | angle = cellMom2.Angle(cellMom1.Vect()); | |
478 | ||
479 | if (mass < fMassPhoMax && mass > fMassPhoMin) return kPhoton; | |
480 | else if(mass < fMassPi0Max && mass > fMassPi0Min) return kPi0; | |
481 | else if(mass < fMassEtaMax && mass > fMassEtaMin) return kEta; | |
482 | else return kNeutralUnknown; | |
483 | ||
3c1d9afb | 484 | } |
485 | ||
49b5c49b | 486 | //_________________________________________ |
487 | TString AliCaloPID::GetPIDParametersList() | |
488 | { | |
477d6cee | 489 | //Put data member values in string to keep in output container |
490 | ||
491 | TString parList ; //this will be list of parameters used for this analysis. | |
5ae09196 | 492 | const Int_t buffersize = 255; |
493 | char onePar[buffersize] ; | |
494 | snprintf(onePar,buffersize,"--- AliCaloPID ---\n") ; | |
477d6cee | 495 | parList+=onePar ; |
49b5c49b | 496 | if(fUseBayesianWeights){ |
497 | snprintf(onePar,buffersize,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)\n",fEMCALPhotonWeight) ; | |
498 | parList+=onePar ; | |
499 | snprintf(onePar,buffersize,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)\n",fEMCALPi0Weight) ; | |
500 | parList+=onePar ; | |
501 | snprintf(onePar,buffersize,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)\n",fEMCALElectronWeight) ; | |
502 | parList+=onePar ; | |
503 | snprintf(onePar,buffersize,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)\n",fEMCALChargeWeight) ; | |
504 | parList+=onePar ; | |
505 | snprintf(onePar,buffersize,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)\n",fEMCALNeutralWeight) ; | |
506 | parList+=onePar ; | |
507 | snprintf(onePar,buffersize,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)\n",fPHOSPhotonWeight) ; | |
508 | parList+=onePar ; | |
509 | snprintf(onePar,buffersize,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)\n",fPHOSPi0Weight) ; | |
510 | parList+=onePar ; | |
511 | snprintf(onePar,buffersize,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)\n",fPHOSElectronWeight) ; | |
512 | parList+=onePar ; | |
513 | snprintf(onePar,buffersize,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)\n",fPHOSChargeWeight) ; | |
514 | parList+=onePar ; | |
515 | snprintf(onePar,buffersize,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)\n",fPHOSNeutralWeight) ; | |
516 | parList+=onePar ; | |
517 | ||
518 | if(fPHOSWeightFormula){ | |
519 | snprintf(onePar,buffersize,"PHOS Photon Weight Formula: %s\n",fPHOSPhotonWeightFormulaExpression.Data() ) ; | |
520 | parList+=onePar; | |
521 | snprintf(onePar,buffersize,"PHOS Pi0 Weight Formula: %s\n",fPHOSPi0WeightFormulaExpression.Data() ) ; | |
522 | parList+=onePar; | |
523 | } | |
524 | } | |
525 | else { | |
526 | snprintf(onePar,buffersize,"EMCAL: fEMCALL0CutMin =%2.2f, fEMCALL0CutMax =%2.2f (Cut on Shower Shape) \n",fEMCALL0CutMin, fEMCALL0CutMax) ; | |
527 | parList+=onePar ; | |
528 | snprintf(onePar,buffersize,"EMCAL: fEMCALDEtaCut =%2.2f, fEMCALDPhiCut =%2.2f (Cut on track matching) \n",fEMCALDEtaCut, fEMCALDPhiCut) ; | |
529 | parList+=onePar ; | |
530 | snprintf(onePar,buffersize,"fTOFCut =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ; | |
531 | parList+=onePar ; | |
532 | snprintf(onePar,buffersize,"fPHOSRCut =%2.2f, fPHOSDispersionCut =%2.2f (Cut on Shower Shape and CPV) \n",fPHOSRCut,fPHOSDispersionCut) ; | |
533 | parList+=onePar ; | |
534 | ||
a5fb4114 | 535 | } |
477d6cee | 536 | |
3c1d9afb | 537 | if(fDoClusterSplitting) |
538 | { | |
539 | snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n", fSplitM02MinCut, fSplitM02MaxCut) ; | |
540 | parList+=onePar ; | |
541 | snprintf(onePar,buffersize,"fMinNCells =%d \n", fSplitMinNCells) ; | |
542 | parList+=onePar ; | |
543 | snprintf(onePar,buffersize,"pi0 : %2.1f < m <%2.1f\n", fMassPi0Min,fMassPi0Max); | |
544 | parList+=onePar ; | |
545 | snprintf(onePar,buffersize,"eta : %2.1f < m <%2.1f\n", fMassEtaMin,fMassEtaMax); | |
546 | parList+=onePar ; | |
547 | snprintf(onePar,buffersize,"conv: %2.1f < m <%2.1f\n", fMassPhoMin,fMassPhoMax); | |
548 | parList+=onePar ; | |
549 | } | |
550 | ||
477d6cee | 551 | return parList; |
552 | ||
1c5acb87 | 553 | } |
554 | ||
49b5c49b | 555 | //________________________________________________ |
1c5acb87 | 556 | void AliCaloPID::Print(const Option_t * opt) const |
557 | { | |
477d6cee | 558 | |
559 | //Print some relevant parameters set for the analysis | |
560 | if(! opt) | |
561 | return; | |
562 | ||
563 | printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ; | |
564 | ||
3c1d9afb | 565 | if(fUseBayesianWeights) |
566 | { | |
49b5c49b | 567 | printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n", |
3c1d9afb | 568 | fPHOSPhotonWeight, fPHOSPi0Weight, |
569 | fPHOSElectronWeight, fPHOSChargeWeight, fPHOSNeutralWeight) ; | |
49b5c49b | 570 | printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n", |
3c1d9afb | 571 | fEMCALPhotonWeight, fEMCALPi0Weight, |
572 | fEMCALElectronWeight, fEMCALChargeWeight, fEMCALNeutralWeight) ; | |
49b5c49b | 573 | |
574 | printf("PHOS Parametrized weight on? = %d\n", fPHOSWeightFormula) ; | |
3c1d9afb | 575 | if(fPHOSWeightFormula) |
576 | { | |
49b5c49b | 577 | printf("Photon weight formula = %s\n", fPHOSPhotonWeightFormulaExpression.Data()); |
578 | printf("Pi0 weight formula = %s\n", fPHOSPi0WeightFormulaExpression .Data()); | |
579 | } | |
580 | if(fRecalculateBayesian) printf(" Recalculate bayesian with Particle Flux? = %d\n",fParticleFlux); | |
581 | } | |
3c1d9afb | 582 | else |
583 | { | |
584 | printf("TOF cut = %e\n", fTOFCut); | |
585 | printf("EMCAL Lambda0 cut min = %2.2f; max = %2.2f\n", fEMCALL0CutMin,fEMCALL0CutMax); | |
586 | printf("EMCAL cluster-track dEta < %2.3f; dPhi < %2.3f\n", fEMCALDEtaCut, fEMCALDPhiCut); | |
587 | printf("PHOS Treac matching cut =%2.2f, Dispersion Cut =%2.2f \n",fPHOSRCut, fPHOSDispersionCut) ; | |
49b5c49b | 588 | |
a5fb4114 | 589 | } |
477d6cee | 590 | |
3c1d9afb | 591 | if(fDoClusterSplitting) |
592 | { | |
593 | printf("Min. N Cells =%d \n", fSplitMinNCells) ; | |
594 | printf("%2.2f < lambda_0^2 <%2.2f \n",fSplitM02MinCut,fSplitM02MaxCut); | |
595 | printf("pi0 : %2.2f<m<%2.2f \n", fMassPi0Min,fMassPi0Max); | |
596 | printf("eta : %2.2f<m<%2.2f \n", fMassEtaMin,fMassEtaMax); | |
597 | printf("phot: %2.2f<m<%2.2f \n", fMassPhoMin,fMassPhoMax); | |
598 | } | |
599 | ||
477d6cee | 600 | printf(" \n"); |
601 | ||
1c5acb87 | 602 | } |
603 | ||
3c1d9afb | 604 | //_________________________________________________________________ |
605 | void AliCaloPID::PrintClusterPIDWeights(const Double_t * pid) const | |
606 | { | |
607 | // print PID of cluster, (AliVCluster*)cluster->GetPID() | |
608 | ||
c2791479 | 609 | printf("AliCaloPID::PrintClusterPIDWeights() \n \t ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, \n \t \ |
3c1d9afb | 610 | pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n", |
611 | pid[AliVCluster::kPhoton], pid[AliVCluster::kPi0], | |
612 | pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon], | |
613 | pid[AliVCluster::kPion], pid[AliVCluster::kKaon], | |
614 | pid[AliVCluster::kProton], | |
615 | pid[AliVCluster::kNeutron], pid[AliVCluster::kKaon0]); | |
616 | ||
617 | } | |
618 | ||
49b5c49b | 619 | //___________________________________________________________________________ |
3c1d9afb | 620 | void AliCaloPID::SetPIDBits(AliVCluster * cluster, |
49b5c49b | 621 | AliAODPWG4Particle * ph, AliCalorimeterUtils* cu, |
622 | AliVEvent* event) | |
623 | { | |
477d6cee | 624 | //Set Bits for PID selection |
625 | ||
626 | //Dispersion/lambdas | |
5ae09196 | 627 | //Double_t disp= cluster->GetDispersion() ; |
628 | Double_t l1 = cluster->GetM20() ; | |
629 | Double_t l0 = cluster->GetM02() ; | |
630 | Bool_t isDispOK = kTRUE ; | |
9a6fa057 | 631 | if(cluster->IsPHOS()){ |
49b5c49b | 632 | if(TestPHOSDispersion(ph->Pt(),l0,l1) < fPHOSDispersionCut) isDispOK = kTRUE; |
633 | else isDispOK = kFALSE; | |
5ae09196 | 634 | } |
635 | else{//EMCAL | |
636 | ||
49b5c49b | 637 | if(l0 > fEMCALL0CutMin && l0 < fEMCALL0CutMax) isDispOK = kTRUE; |
638 | ||
5ae09196 | 639 | } |
640 | ||
641 | ph->SetDispBit(isDispOK) ; | |
477d6cee | 642 | |
643 | //TOF | |
644 | Double_t tof=cluster->GetTOF() ; | |
645 | ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ; | |
646 | ||
49b5c49b | 647 | //Charged |
648 | Bool_t isNeutral = IsTrackMatched(cluster,cu,event); | |
5ae09196 | 649 | |
650 | ph->SetChargedBit(isNeutral); | |
477d6cee | 651 | |
652 | //Set PID pdg | |
3c1d9afb | 653 | ph->SetIdentifiedParticleType(GetIdentifiedParticleType(cluster)); |
477d6cee | 654 | |
655 | if(fDebug > 0){ | |
5ae09196 | 656 | printf("AliCaloPID::SetPIDBits: TOF %e, Lambda0 %2.2f, Lambda1 %2.2f\n",tof , l0, l1); |
477d6cee | 657 | printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n", |
49b5c49b | 658 | ph->GetIdentifiedParticleType(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit()); |
477d6cee | 659 | } |
1c5acb87 | 660 | } |
661 | ||
09273901 | 662 | //_________________________________________________________ |
49b5c49b | 663 | Bool_t AliCaloPID::IsTrackMatched(AliVCluster* cluster, |
664 | AliCalorimeterUtils * cu, | |
665 | AliVEvent* event) const | |
666 | { | |
5ae09196 | 667 | //Check if there is any track attached to this cluster |
668 | ||
669 | Int_t nMatches = cluster->GetNTracksMatched(); | |
49b5c49b | 670 | AliVTrack * track = 0; |
671 | Double_t p[3]; | |
672 | ||
673 | if(nMatches > 0){ | |
ae182e60 | 674 | |
49b5c49b | 675 | //In case of ESDs, by default without match one entry with negative index, no match, reject. |
676 | if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName()))) | |
677 | { | |
678 | Int_t iESDtrack = cluster->GetTrackMatchedIndex(); | |
679 | if(iESDtrack >= 0) track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack)); | |
680 | else return kFALSE; | |
d39cba7e | 681 | |
49b5c49b | 682 | if (!track){ |
683 | printf("AliCaloPID::IsTrackMatched() - Null matched track in ESD when index is OK!\n"); | |
684 | return kFALSE; | |
685 | } | |
686 | } | |
687 | else { // AOD | |
688 | track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0)); | |
689 | if (!track){ | |
690 | printf("AliCaloPID::IsTrackMatched() - Null matched track in AOD!\n"); | |
691 | return kFALSE; | |
c76f0089 | 692 | } |
d39cba7e | 693 | } |
5ae09196 | 694 | |
49b5c49b | 695 | Float_t dZ = cluster->GetTrackDz(); |
696 | Float_t dR = cluster->GetTrackDx(); | |
697 | ||
698 | // if track matching was recalculated | |
09273901 | 699 | if(cluster->IsEMCAL() && cu && cu->IsRecalculationOfClusterTrackMatchingOn()){ |
49b5c49b | 700 | dR = 2000., dZ = 2000.; |
31ae6d59 | 701 | cu->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR); |
49b5c49b | 702 | } |
09273901 | 703 | |
49b5c49b | 704 | if(cluster->IsPHOS()) { |
705 | ||
706 | track->GetPxPyPz(p) ; | |
707 | TLorentzVector trackmom(p[0],p[1],p[2],0); | |
708 | Int_t charge = track->Charge(); | |
709 | Double_t mf = event->GetMagneticField(); | |
710 | if(TestPHOSChargedVeto(dR, dZ, trackmom.Pt(), charge, mf ) < fPHOSRCut) return kTRUE; | |
711 | else return kFALSE; | |
712 | ||
713 | }//PHOS | |
714 | else {//EMCAL | |
715 | ||
09273901 | 716 | if(fDebug > 0) |
49b5c49b | 717 | printf("AliCaloPID::IsTrackMatched - EMCAL dR %f < %f, dZ %f < %f \n",dR, fEMCALDPhiCut, dZ, fEMCALDEtaCut); |
718 | ||
719 | if(TMath::Abs(dR) < fEMCALDPhiCut && | |
720 | TMath::Abs(dZ) < fEMCALDEtaCut) return kTRUE; | |
721 | else return kFALSE; | |
722 | ||
723 | }//EMCAL cluster | |
724 | ||
725 | ||
726 | } // more than 1 match, at least one track in array | |
727 | else return kFALSE; | |
728 | ||
729 | } | |
730 | ||
731 | //___________________________________________________________________________________________________ | |
732 | Float_t AliCaloPID::TestPHOSDispersion(const Double_t pt, const Double_t l1, const Double_t l2) const | |
733 | { | |
734 | //Check if cluster photon-like. Uses photon cluster parameterization in real pp data | |
735 | //Returns distance in sigmas. Recommended cut 2.5 | |
736 | ||
737 | Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ; | |
738 | Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ; | |
739 | Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt; | |
740 | Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt; | |
741 | Double_t c =-0.35-0.550*TMath::Exp(-0.390730*pt) ; | |
f3138ecf | 742 | Double_t r2 = 0.5* (l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma + |
49b5c49b | 743 | 0.5* (l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma + |
744 | 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ; | |
745 | ||
f3138ecf | 746 | if(fDebug > 0) printf("AliCaloPID::TestPHOSDispersion() - PHOS SS R %f < %f?\n", TMath::Sqrt(r2), fPHOSDispersionCut); |
5ae09196 | 747 | |
f3138ecf | 748 | return TMath::Sqrt(r2) ; |
5ae09196 | 749 | |
750 | } | |
751 | ||
49b5c49b | 752 | //_______________________________________________________________________________________________ |
753 | Float_t AliCaloPID::TestPHOSChargedVeto(const Double_t dx, const Double_t dz, const Double_t pt, | |
754 | const Int_t charge, const Double_t mf) const | |
755 | { | |
756 | //Checks distance to the closest track. Takes into account | |
757 | //non-perpendicular incidence of tracks. | |
758 | //returns distance in sigmas. Recommended cut: 2. | |
759 | //Requires (sign) of magnetic filed. onc can find it for example as following | |
760 | // Double_t mf=0. ; | |
761 | // AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent()); | |
762 | // if(event) | |
763 | // mf = event->GetMagneticField(); //Positive for ++ and negative for -- | |
764 | ||
765 | ||
766 | Double_t meanX = 0.; | |
767 | Double_t meanZ = 0.; | |
768 | Double_t sx = TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+ | |
769 | 6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+ | |
770 | 1.59219); | |
771 | Double_t sz = TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+ | |
772 | 1.60) ; | |
773 | ||
774 | if(mf<0.){ //field -- | |
775 | meanZ = -0.468318 ; | |
776 | if(charge>0) | |
777 | meanX = TMath::Min(7.3, 3.89994*1.20679 *1.20679 /(pt*pt+1.20679*1.20679)+ | |
778 | 0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ; | |
779 | else | |
780 | meanX =-TMath::Min(7.7, 3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+ | |
781 | 1.23114 +4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ; | |
782 | } | |
783 | else{ //Field ++ | |
784 | meanZ = -0.468318; | |
785 | if(charge>0) | |
786 | meanX =-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+ | |
787 | 0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ; | |
788 | else | |
789 | meanX = TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)- | |
790 | 0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ; | |
791 | } | |
9a6fa057 | 792 | |
49b5c49b | 793 | Double_t rz = (dz-meanZ)/sz ; |
794 | Double_t rx = (dx-meanX)/sx ; | |
9a6fa057 | 795 | |
49b5c49b | 796 | if(fDebug > 0) |
797 | printf("AliCaloPID::TestPHOSDispersion() - PHOS Matching R %f < %f\n",TMath::Sqrt(rx*rx+rz*rz), fPHOSRCut); | |
9a6fa057 | 798 | |
49b5c49b | 799 | return TMath::Sqrt(rx*rx+rz*rz) ; |
9a6fa057 | 800 | |
801 | } |