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