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