| 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 |
| 247 | fMassWidthPi0Param[3] = 0.0023; // pol1 param0 of width for E > 8 GeV |
| 248 | fMassWidthPi0Param[4] = 0.0008; // pol1 param1 of width for E > 8 GeV |
| 249 | fMassWidthPi0Param[5] = 0.130; // Mean mass value for NLM=1 |
| 250 | fMassWidthPi0Param[6] = 0.134; // Mean mass value for NLM=2 |
| 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 | |
| 356 | // Get the parametrized width of the mass |
| 357 | Float_t width = 0.009; |
| 358 | if(energy < 8) width = fMassWidthPi0Param[2]; |
| 359 | else width = fMassWidthPi0Param[3]+energy*fMassWidthPi0Param[4]; |
| 360 | |
| 361 | // Calculate the 2 sigma cut |
| 362 | Float_t minMass = meanMass-fSplitWidthSigma*width; |
| 363 | Float_t maxMass = meanMass+fSplitWidthSigma*width; |
| 364 | |
| 365 | // In case of low energy, hard cut to avoid conversions |
| 5b4c2f5b |
366 | if(energy < 10 && nlm == 1 && minMass < fMassWidthPi0Param[0] ) minMass = fMassWidthPi0Param[0]; |
| 367 | if(energy < 10 && nlm == 2 && minMass < fMassWidthPi0Param[1] ) minMass = fMassWidthPi0Param[1]; |
| 5a72d9af |
368 | |
| 369 | //printf("\t \t sigma %1.1f width %3.1f, mean Mass %3.0f, minMass %3.0f, maxMass %3.0f\n ", |
| 370 | // fSplitWidthSigma, width*1000, meanMass*1000,minMass*1000,maxMass*1000); |
| 371 | |
| 372 | if(mass < maxMass && mass > minMass) return kTRUE; |
| 373 | else return kFALSE; |
| 5a72d9af |
374 | |
| 375 | } |
| 376 | |
| 377 | //_____________________________________________________________________________________________ |
| 5b4c2f5b |
378 | Bool_t AliCaloPID::IsInPi0M02Range(const Float_t energy, const Float_t m02, const Int_t nlm) |
| 5a72d9af |
379 | { |
| 5b4c2f5b |
380 | // Select the appropriate m02 range in splitting method for pi0 |
| 5a72d9af |
381 | |
| 382 | Float_t minCut = fSplitM02MinCut; |
| 5b4c2f5b |
383 | Float_t maxCut = fSplitM02MaxCut; |
| 3c1d9afb |
384 | |
| 5a72d9af |
385 | if(!fUseSimpleM02Cut) |
| 386 | { |
| a5a3f703 |
387 | Int_t inlm = nlm-1; |
| 388 | if(nlm > 2) inlm=1; // only 2 cases defined nlm=1 and nlm>=2 |
| 5a72d9af |
389 | |
| 667432ef |
390 | minCut = fM02MinParam[inlm][0]+ |
| 5b4c2f5b |
391 | fM02MinParam[inlm][1]*energy+ |
| 392 | fM02MinParam[inlm][2]*energy*energy+ |
| 393 | fM02MinParam[inlm][3]*energy*energy*energy; |
| 394 | |
| 395 | maxCut = fM02MaxParam[inlm][0]+ |
| 396 | fM02MaxParam[inlm][1]*energy+ |
| 397 | fM02MaxParam[inlm][2]*energy*energy+ |
| 398 | fM02MaxParam[inlm][3]*energy*energy*energy; |
| 399 | |
| 667432ef |
400 | |
| 401 | // In any case and beyond validity energy range of the function, |
| 402 | // the parameter cannot be smaller than this (0.3 for nlm=1 and 0.6 for the rest) |
| 403 | if( minCut < fM02MinParam[inlm][4] || energy > fM02MinParam[inlm][5] ) minCut = fM02MinParam[inlm][4]; |
| 5b4c2f5b |
404 | if( maxCut < fM02MaxParam[inlm][4] || energy > fM02MaxParam[inlm][5] ) maxCut = fM02MaxParam[inlm][4]; |
| 405 | if( nlm>2 ) maxCut+=0.75; |
| 5a72d9af |
406 | } |
| 407 | |
| 5b4c2f5b |
408 | //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 |
409 | |
| 5b4c2f5b |
410 | if(m02 < maxCut && m02 > minCut) return kTRUE; |
| 411 | else return kFALSE; |
| 5a72d9af |
412 | |
| 49b5c49b |
413 | } |
| 414 | |
| 5a72d9af |
415 | |
| 5b4c2f5b |
416 | //_____________________________________________________________________________________________ |
| 417 | Bool_t AliCaloPID::IsInEtaM02Range(const Float_t energy, const Float_t m02, const Int_t nlm) |
| 418 | { |
| 419 | // Select the appropriate m02 range in splitting method to select eta's |
| 420 | // Use same parametrization as pi0, just shift the distributions (to be tuned) |
| 421 | |
| 422 | Float_t minCut = fSplitM02MinCut; |
| 423 | Float_t maxCut = fSplitM02MaxCut; |
| 424 | |
| 425 | if(!fUseSimpleM02Cut) |
| 426 | { |
| 427 | Int_t inlm = nlm-1; |
| 428 | if(nlm > 2) inlm=1; // only 2 cases defined nlm=1 and nlm>=2 |
| 429 | |
| 430 | Float_t shiftE = energy-20; // to be tuned |
| 431 | if(nlm==1) shiftE=energy-28; |
| 432 | |
| 433 | minCut = fM02MinParam[inlm][0]+ |
| 434 | fM02MinParam[inlm][1]*shiftE+ |
| 435 | fM02MinParam[inlm][2]*shiftE*shiftE+ |
| 436 | fM02MinParam[inlm][3]*shiftE*shiftE*shiftE; |
| 437 | |
| 438 | // In any case and beyond validity energy range of the function, |
| 439 | // the parameter cannot be smaller than this (0.3 for nlm=1 and 0.6 for the rest) |
| 440 | Float_t minE = 50; |
| 441 | if(nlm>1) minE = 35; |
| 442 | if( minCut < fM02MinParam[inlm][4]+1 || shiftE > minE ) minCut = fM02MinParam[inlm][4]+1; |
| 443 | |
| 444 | shiftE = energy+20; // to be tuned |
| 445 | |
| 446 | maxCut = 1+fM02MaxParam[inlm][0]+ |
| 447 | fM02MaxParam[inlm][1]*shiftE+ |
| 448 | fM02MaxParam[inlm][2]*shiftE*shiftE+ |
| 449 | fM02MaxParam[inlm][3]*shiftE*shiftE*shiftE; |
| 450 | |
| 451 | Float_t maxE = 50; |
| 452 | if(nlm>1) maxE = 40; |
| 453 | if( maxCut < fM02MaxParam[inlm][4]+1 || shiftE > maxE ) maxCut = fM02MaxParam[inlm][4]+1; |
| 454 | if( nlm>2 ) maxCut+=0.75; |
| 455 | } |
| 456 | |
| 457 | //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); |
| 458 | |
| 459 | if(m02 < maxCut && m02 > minCut) return kTRUE; |
| 460 | else return kFALSE; |
| 461 | |
| 462 | } |
| 463 | |
| 464 | //_____________________________________________________________________________________________ |
| 465 | Bool_t AliCaloPID::IsInConM02Range(const Float_t energy, const Float_t m02, const Int_t nlm) |
| 466 | { |
| 467 | // Select the appropriate m02 range in splitting method for converted photons |
| 468 | // Just min limit for pi0s is max for conversion. |
| 469 | |
| 470 | Float_t minCut = 0.1; |
| 471 | Float_t maxCut = 0.3; |
| 472 | |
| 473 | if(!fUseSimpleM02Cut) |
| 474 | { |
| 475 | Int_t inlm = nlm-1; |
| 476 | if(nlm > 2) inlm=1; // only 2 cases defined nlm=1 and nlm>=2 |
| 477 | |
| 478 | maxCut = fM02MinParam[inlm][0]+ |
| 479 | fM02MinParam[inlm][1]*energy+ |
| 480 | fM02MinParam[inlm][2]*energy*energy+ |
| 481 | fM02MinParam[inlm][3]*energy*energy*energy; |
| 482 | |
| 483 | // In any case and beyond validity energy range of the function, |
| 484 | // the parameter cannot be smaller than this (0.3 for nlm=1 and 0.6 for the rest) |
| 485 | if( maxCut < fM02MaxParam[inlm][4] || energy > fM02MaxParam[inlm][5] ) maxCut = fM02MaxParam[inlm][4]; |
| 486 | } |
| 487 | |
| 488 | |
| 489 | if(m02 < maxCut && m02 > minCut) return kTRUE; |
| 490 | else return kFALSE; |
| 491 | |
| 492 | } |
| 493 | |
| c5693f62 |
494 | //______________________________________________ |
| 495 | AliEMCALPIDUtils *AliCaloPID::GetEMCALPIDUtils() |
| 496 | { |
| 497 | // return pointer to AliEMCALPIDUtils, create it if needed |
| 498 | |
| 499 | if(!fEMCALPIDUtils) fEMCALPIDUtils = new AliEMCALPIDUtils ; |
| 500 | return fEMCALPIDUtils ; |
| 501 | |
| 502 | } |
| 503 | |
| 504 | |
| 49b5c49b |
505 | //______________________________________________________________________ |
| 3c1d9afb |
506 | Int_t AliCaloPID::GetIdentifiedParticleType(const AliVCluster * cluster) |
| 49b5c49b |
507 | { |
| 508 | // Returns a PDG number corresponding to the likely ID of the cluster |
| 509 | |
| 3c1d9afb |
510 | Float_t energy = cluster->E(); |
| 49b5c49b |
511 | Float_t lambda0 = cluster->GetM02(); |
| 512 | Float_t lambda1 = cluster->GetM20(); |
| 513 | |
| 514 | // --------------------- |
| 515 | // Use bayesian approach |
| 516 | // --------------------- |
| 517 | |
| 3c1d9afb |
518 | if(fUseBayesianWeights) |
| 519 | { |
| 00a38d07 |
520 | Double_t weights[AliPID::kSPECIESCN]; |
| 49b5c49b |
521 | |
| 3c1d9afb |
522 | if(cluster->IsEMCAL() && fRecalculateBayesian) |
| 523 | { |
| 49b5c49b |
524 | fEMCALPIDUtils->ComputePID(energy, lambda0); |
| 00a38d07 |
525 | for(Int_t i = 0; i < AliPID::kSPECIESCN; i++) weights[i] = fEMCALPIDUtils->GetPIDFinal(i); |
| 49b5c49b |
526 | } |
| 3c1d9afb |
527 | else |
| 528 | { |
| 00a38d07 |
529 | for(Int_t i = 0; i < AliPID::kSPECIESCN; i++) weights[i] = cluster->GetPID()[i]; |
| 49b5c49b |
530 | } |
| 531 | |
| 3c1d9afb |
532 | if(fDebug > 0) PrintClusterPIDWeights(weights); |
| 49b5c49b |
533 | |
| 3c1d9afb |
534 | return GetIdentifiedParticleTypeFromBayesWeights(cluster->IsEMCAL(), weights, energy); |
| bdd2a262 |
535 | } |
| 49b5c49b |
536 | |
| 537 | // ------------------------------------------------------- |
| 538 | // Calculate PID SS from data, do not use bayesian weights |
| 539 | // ------------------------------------------------------- |
| 540 | |
| 3c1d9afb |
541 | 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", |
| 542 | cluster->IsEMCAL(),energy,lambda0,cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(), |
| 49b5c49b |
543 | cluster->GetEmcCpvDistance(), cluster->GetDistanceToBadChannel(),cluster->GetNExMax()); |
| 544 | |
| 3c1d9afb |
545 | if(cluster->IsEMCAL()) |
| 546 | { |
| 49b5c49b |
547 | if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType() - EMCAL SS %f <%f < %f?\n",fEMCALL0CutMin, lambda0, fEMCALL0CutMax); |
| 548 | |
| 549 | if(lambda0 < fEMCALL0CutMax && lambda0 > fEMCALL0CutMin) return kPhoton ; |
| 550 | else return kNeutralUnknown ; |
| 3c1d9afb |
551 | } // EMCAL |
| 552 | else // PHOS |
| 553 | { |
| 554 | if(TestPHOSDispersion(energy,lambda0,lambda1) < fPHOSDispersionCut) return kPhoton; |
| 555 | else return kNeutralUnknown; |
| 49b5c49b |
556 | } |
| 557 | |
| 1c5acb87 |
558 | } |
| 559 | |
| 49b5c49b |
560 | //_______________________________________________________________________________ |
| 3c1d9afb |
561 | Int_t AliCaloPID::GetIdentifiedParticleTypeFromBayesWeights(const Bool_t isEMCAL, |
| 49b5c49b |
562 | const Double_t * pid, |
| 563 | const Float_t energy) |
| 564 | { |
| 565 | //Return most probable identity of the particle after bayesian weights calculated in reconstruction |
| 477d6cee |
566 | |
| 3c1d9afb |
567 | if(!pid) |
| 568 | { |
| 21a4b1c0 |
569 | printf("AliCaloPID::GetIdentifiedParticleType() - pid pointer not initialized!!!\n"); |
| 477d6cee |
570 | abort(); |
| 571 | } |
| 572 | |
| 15800db4 |
573 | Float_t wPh = fPHOSPhotonWeight ; |
| 477d6cee |
574 | Float_t wPi0 = fPHOSPi0Weight ; |
| 15800db4 |
575 | Float_t wE = fPHOSElectronWeight ; |
| 576 | Float_t wCh = fPHOSChargeWeight ; |
| 577 | Float_t wNe = fPHOSNeutralWeight ; |
| 49b5c49b |
578 | |
| 3c1d9afb |
579 | if(!isEMCAL && fPHOSWeightFormula){ |
| a5fb4114 |
580 | wPh = GetPHOSPhotonWeightFormula()->Eval(energy) ; |
| 581 | wPi0 = GetPHOSPi0WeightFormula() ->Eval(energy); |
| 582 | } |
| 3c1d9afb |
583 | else |
| 584 | { |
| 477d6cee |
585 | wPh = fEMCALPhotonWeight ; |
| 586 | wPi0 = fEMCALPi0Weight ; |
| 587 | wE = fEMCALElectronWeight ; |
| 588 | wCh = fEMCALChargeWeight ; |
| 589 | wNe = fEMCALNeutralWeight ; |
| 477d6cee |
590 | } |
| 591 | |
| 3c1d9afb |
592 | if(fDebug > 0) PrintClusterPIDWeights(pid); |
| 593 | |
| 477d6cee |
594 | Int_t pdg = kNeutralUnknown ; |
| c8fe2783 |
595 | Float_t chargedHadronWeight = pid[AliVCluster::kProton]+pid[AliVCluster::kKaon]+ |
| 49b5c49b |
596 | pid[AliVCluster::kPion]+pid[AliVCluster::kMuon]; |
| c8fe2783 |
597 | Float_t neutralHadronWeight = pid[AliVCluster::kNeutron]+pid[AliVCluster::kKaon0]; |
| 598 | Float_t allChargedWeight = pid[AliVCluster::kElectron]+pid[AliVCluster::kEleCon]+ chargedHadronWeight; |
| 599 | Float_t allNeutralWeight = pid[AliVCluster::kPhoton]+pid[AliVCluster::kPi0]+ neutralHadronWeight; |
| 477d6cee |
600 | |
| 601 | //Select most probable ID |
| 3c1d9afb |
602 | if(!isEMCAL) // PHOS |
| 603 | { |
| a5fb4114 |
604 | if(pid[AliVCluster::kPhoton] > wPh) pdg = kPhoton ; |
| 605 | else if(pid[AliVCluster::kPi0] > wPi0) pdg = kPi0 ; |
| c8fe2783 |
606 | else if(pid[AliVCluster::kElectron] > wE) pdg = kElectron ; |
| a5fb4114 |
607 | else if(pid[AliVCluster::kEleCon] > wE) pdg = kEleCon ; |
| 608 | else if(chargedHadronWeight > wCh) pdg = kChargedHadron ; |
| 609 | else if(neutralHadronWeight > wNe) pdg = kNeutralHadron ; |
| 477d6cee |
610 | else if(allChargedWeight > allNeutralWeight) |
| 611 | pdg = kChargedUnknown ; |
| 612 | else |
| 613 | pdg = kNeutralUnknown ; |
| 614 | } |
| 3c1d9afb |
615 | else //EMCAL |
| 616 | { |
| 2007809d |
617 | if(pid[AliVCluster::kPhoton] > wPh) pdg = kPhoton ; |
| 618 | else if(pid[AliVCluster::kElectron] > wE) pdg = kElectron ; |
| 619 | else if(pid[AliVCluster::kPhoton]+pid[AliVCluster::kElectron] > wPh) pdg = kPhoton ; //temporal sollution until track matching for electrons is considered |
| 620 | else if(pid[AliVCluster::kPi0] > wPi0) pdg = kPi0 ; |
| 477d6cee |
621 | else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ; |
| 622 | else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ; |
| 2007809d |
623 | else pdg = kNeutralUnknown ; |
| 477d6cee |
624 | } |
| 625 | |
| 21a4b1c0 |
626 | if(fDebug > 0)printf("AliCaloPID::GetIdentifiedParticleType:Final Pdg: %d, cluster energy %2.2f \n", pdg,energy); |
| 1c5acb87 |
627 | |
| 49b5c49b |
628 | return pdg ; |
| 9a6fa057 |
629 | |
| 1c5acb87 |
630 | } |
| 631 | |
| 3c1d9afb |
632 | //____________________________________________________________________________________________________ |
| 633 | Int_t AliCaloPID::GetIdentifiedParticleTypeFromClusterSplitting(AliVCluster* cluster, |
| 634 | AliVCaloCells* cells, |
| 635 | AliCalorimeterUtils * caloutils, |
| 636 | Double_t vertex[3], |
| 637 | Int_t & nMax, |
| bfdcf7fb |
638 | Double_t & mass, Double_t & angle, |
| 639 | Double_t & e1 , Double_t & e2 ) |
| 3c1d9afb |
640 | { |
| 641 | // Split the cluster in 2, do invariant mass, get the mass and decide |
| 642 | // if this is a photon, pi0, eta, ... |
| 643 | |
| 2bf17171 |
644 | Int_t absId1 = -1; Int_t absId2 = -1; |
| 5a72d9af |
645 | Float_t eClus = cluster->E(); |
| 646 | Float_t m02 = cluster->GetM02(); |
| 2bf17171 |
647 | const Int_t nc = cluster->GetNCells(); |
| 648 | Int_t absIdList[nc]; |
| 649 | Float_t maxEList [nc]; |
| 3c1d9afb |
650 | |
| 3c1d9afb |
651 | mass = -1.; |
| 652 | angle = -1.; |
| 667432ef |
653 | |
| cf7e2ca9 |
654 | //If too low number of cells, skip it |
| 655 | if ( nc < fSplitMinNCells) return kNeutralUnknown ; |
| 656 | |
| 657 | if(fDebug > 0) printf("\t pass nCells cut\n"); |
| 658 | |
| 3c1d9afb |
659 | // Get Number of local maxima |
| 5a72d9af |
660 | nMax = caloutils->GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ; |
| 3c1d9afb |
661 | |
| 5a72d9af |
662 | if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleTypeFromClusterSplitting() - Cluster : E %1.1f, M02 %1.2f, NLM %d, N Cells %d\n", |
| 663 | eClus,m02,nMax,nc); |
| 664 | |
| 3c1d9afb |
665 | //--------------------------------------------------------------------- |
| 666 | // Get the 2 max indeces and do inv mass |
| 667 | //--------------------------------------------------------------------- |
| 668 | |
| 669 | if ( nMax == 2 ) |
| 670 | { |
| 671 | absId1 = absIdList[0]; |
| 672 | absId2 = absIdList[1]; |
| 673 | } |
| 674 | else if( nMax == 1 ) |
| 675 | { |
| 676 | |
| 677 | absId1 = absIdList[0]; |
| 678 | |
| 679 | //Find second highest energy cell |
| 680 | |
| 681 | TString calorimeter = "EMCAL"; |
| 682 | if(cluster->IsPHOS()) calorimeter = "PHOS"; |
| 683 | Float_t enmax = 0 ; |
| 684 | for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++) |
| 685 | { |
| 686 | Int_t absId = cluster->GetCellsAbsId()[iDigit]; |
| 687 | if( absId == absId1 ) continue ; |
| 688 | Float_t endig = cells->GetCellAmplitude(absId); |
| 689 | caloutils->RecalibrateCellAmplitude(endig,calorimeter,absId); |
| 690 | if(endig > enmax) |
| 691 | { |
| 692 | enmax = endig ; |
| 693 | absId2 = absId ; |
| 694 | } |
| 695 | }// cell loop |
| 696 | }// 1 maxima |
| 697 | else |
| 698 | { // n max > 2 |
| 699 | // loop on maxima, find 2 highest |
| 700 | |
| 701 | // First max |
| 702 | Float_t enmax = 0 ; |
| 703 | for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++) |
| 704 | { |
| 705 | Float_t endig = maxEList[iDigit]; |
| 706 | if(endig > enmax) |
| 707 | { |
| 708 | enmax = endig ; |
| 709 | absId1 = absIdList[iDigit]; |
| 710 | } |
| 711 | }// first maxima loop |
| 712 | |
| 713 | // Second max |
| 714 | Float_t enmax2 = 0; |
| 715 | for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++) |
| 716 | { |
| 717 | if(absIdList[iDigit]==absId1) continue; |
| 718 | Float_t endig = maxEList[iDigit]; |
| 719 | if(endig > enmax2) |
| 720 | { |
| 721 | enmax2 = endig ; |
| 722 | absId2 = absIdList[iDigit]; |
| 723 | } |
| 724 | }// second maxima loop |
| 725 | |
| 726 | } // n local maxima > 2 |
| 727 | |
| cf7e2ca9 |
728 | if(absId2<0 || absId1<0) |
| 729 | { |
| 730 | 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", |
| 731 | nMax,absId1,absId2,eClus,nc,m02); |
| 732 | return kNeutralUnknown ; |
| 733 | } |
| 734 | |
| 3c1d9afb |
735 | //--------------------------------------------------------------------- |
| 736 | // Split the cluster energy in 2, around the highest 2 local maxima |
| 737 | //--------------------------------------------------------------------- |
| 738 | |
| 2bf17171 |
739 | AliAODCaloCluster cluster1(0, 0,NULL,0.,NULL,NULL,1,0); |
| 740 | AliAODCaloCluster cluster2(1, 0,NULL,0.,NULL,NULL,1,0); |
| 3c1d9afb |
741 | |
| 2bf17171 |
742 | caloutils->SplitEnergy(absId1,absId2,cluster, cells, &cluster1, &cluster2,nMax); /*absIdList, maxEList,*/ |
| 3c1d9afb |
743 | |
| 744 | TLorentzVector cellMom1; |
| 745 | TLorentzVector cellMom2; |
| 746 | |
| 2bf17171 |
747 | cluster1.GetMomentum(cellMom1,vertex); |
| 748 | cluster2.GetMomentum(cellMom2,vertex); |
| 3c1d9afb |
749 | |
| 750 | mass = (cellMom1+cellMom2).M(); |
| 751 | angle = cellMom2.Angle(cellMom1.Vect()); |
| bfdcf7fb |
752 | e1 = cluster1.E(); |
| 753 | e2 = cluster2.E(); |
| 5b4c2f5b |
754 | |
| 5a72d9af |
755 | // Consider clusters with splitted energy not too different to original cluster energy |
| 756 | if((e1+e2)/eClus < fSplitEFracMin) return kNeutralUnknown ; |
| 757 | |
| 758 | if(fDebug > 0) printf("\t pass Split E frac cut\n"); |
| 5b4c2f5b |
759 | |
| 667432ef |
760 | // Asymmetry of cluster |
| 761 | Float_t asy =-10; |
| 762 | if(e1+e2 > 0) asy = (e1-e2) / (e1+e2); |
| 763 | if( fUseSplitAsyCut && !IsInPi0SplitAsymmetryRange(eClus,asy,nMax) ) return kNeutralUnknown ; |
| 764 | |
| 765 | if (fDebug>0) printf("\t pass asymmetry cut\n"); |
| 5b4c2f5b |
766 | |
| 767 | Bool_t pi0OK = kFALSE; |
| 768 | Bool_t etaOK = kFALSE; |
| 769 | Bool_t conOK = kFALSE; |
| 770 | |
| 771 | //If too small or big M02, skip it |
| 772 | if (IsInPi0M02Range(eClus,m02,nMax)) pi0OK = kTRUE; |
| 773 | else if(IsInEtaM02Range(eClus,m02,nMax)) etaOK = kTRUE; |
| 774 | else if(IsInConM02Range(eClus,m02,nMax)) conOK = kTRUE; |
| 775 | |
| 5a72d9af |
776 | // Check the mass, and set an ID to the splitted cluster |
| 5b4c2f5b |
777 | if ( conOK && mass < fMassPhoMax && mass > fMassPhoMin ) { if(fDebug > 0) printf("\t Split Conv \n"); return kPhoton ; } |
| 778 | else if( etaOK && mass < fMassEtaMax && mass > fMassEtaMin ) { if(fDebug > 0) printf("\t Split Eta \n"); return kEta ; } |
| 779 | else if( pi0OK && IsInPi0SplitMassRange(cluster->E(),mass,nMax)) { if(fDebug > 0) printf("\t Split Pi0 \n"); return kPi0 ; } |
| 780 | else return kNeutralUnknown ; |
| 3c1d9afb |
781 | |
| 3c1d9afb |
782 | } |
| 783 | |
| 49b5c49b |
784 | //_________________________________________ |
| 785 | TString AliCaloPID::GetPIDParametersList() |
| 786 | { |
| 477d6cee |
787 | //Put data member values in string to keep in output container |
| 788 | |
| 789 | TString parList ; //this will be list of parameters used for this analysis. |
| 5ae09196 |
790 | const Int_t buffersize = 255; |
| 791 | char onePar[buffersize] ; |
| 792 | snprintf(onePar,buffersize,"--- AliCaloPID ---\n") ; |
| 477d6cee |
793 | parList+=onePar ; |
| 49b5c49b |
794 | if(fUseBayesianWeights){ |
| 795 | snprintf(onePar,buffersize,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)\n",fEMCALPhotonWeight) ; |
| 796 | parList+=onePar ; |
| 797 | snprintf(onePar,buffersize,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)\n",fEMCALPi0Weight) ; |
| 798 | parList+=onePar ; |
| 799 | snprintf(onePar,buffersize,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)\n",fEMCALElectronWeight) ; |
| 800 | parList+=onePar ; |
| 801 | snprintf(onePar,buffersize,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)\n",fEMCALChargeWeight) ; |
| 802 | parList+=onePar ; |
| 803 | snprintf(onePar,buffersize,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)\n",fEMCALNeutralWeight) ; |
| 804 | parList+=onePar ; |
| 805 | snprintf(onePar,buffersize,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)\n",fPHOSPhotonWeight) ; |
| 806 | parList+=onePar ; |
| 807 | snprintf(onePar,buffersize,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)\n",fPHOSPi0Weight) ; |
| 808 | parList+=onePar ; |
| 809 | snprintf(onePar,buffersize,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)\n",fPHOSElectronWeight) ; |
| 810 | parList+=onePar ; |
| 811 | snprintf(onePar,buffersize,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)\n",fPHOSChargeWeight) ; |
| 812 | parList+=onePar ; |
| 813 | snprintf(onePar,buffersize,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)\n",fPHOSNeutralWeight) ; |
| 814 | parList+=onePar ; |
| 815 | |
| 816 | if(fPHOSWeightFormula){ |
| 817 | snprintf(onePar,buffersize,"PHOS Photon Weight Formula: %s\n",fPHOSPhotonWeightFormulaExpression.Data() ) ; |
| 818 | parList+=onePar; |
| 819 | snprintf(onePar,buffersize,"PHOS Pi0 Weight Formula: %s\n",fPHOSPi0WeightFormulaExpression.Data() ) ; |
| 820 | parList+=onePar; |
| 821 | } |
| 822 | } |
| 823 | else { |
| 824 | snprintf(onePar,buffersize,"EMCAL: fEMCALL0CutMin =%2.2f, fEMCALL0CutMax =%2.2f (Cut on Shower Shape) \n",fEMCALL0CutMin, fEMCALL0CutMax) ; |
| 825 | parList+=onePar ; |
| 826 | snprintf(onePar,buffersize,"EMCAL: fEMCALDEtaCut =%2.2f, fEMCALDPhiCut =%2.2f (Cut on track matching) \n",fEMCALDEtaCut, fEMCALDPhiCut) ; |
| 827 | parList+=onePar ; |
| 828 | snprintf(onePar,buffersize,"fTOFCut =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ; |
| 829 | parList+=onePar ; |
| 830 | snprintf(onePar,buffersize,"fPHOSRCut =%2.2f, fPHOSDispersionCut =%2.2f (Cut on Shower Shape and CPV) \n",fPHOSRCut,fPHOSDispersionCut) ; |
| 831 | parList+=onePar ; |
| 832 | |
| a5fb4114 |
833 | } |
| 477d6cee |
834 | |
| 3c1d9afb |
835 | if(fDoClusterSplitting) |
| 836 | { |
| 837 | snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n", fSplitM02MinCut, fSplitM02MaxCut) ; |
| 838 | parList+=onePar ; |
| 839 | snprintf(onePar,buffersize,"fMinNCells =%d \n", fSplitMinNCells) ; |
| 840 | parList+=onePar ; |
| 841 | snprintf(onePar,buffersize,"pi0 : %2.1f < m <%2.1f\n", fMassPi0Min,fMassPi0Max); |
| 842 | parList+=onePar ; |
| 843 | snprintf(onePar,buffersize,"eta : %2.1f < m <%2.1f\n", fMassEtaMin,fMassEtaMax); |
| 844 | parList+=onePar ; |
| 845 | snprintf(onePar,buffersize,"conv: %2.1f < m <%2.1f\n", fMassPhoMin,fMassPhoMax); |
| 846 | parList+=onePar ; |
| 847 | } |
| 848 | |
| 477d6cee |
849 | return parList; |
| 850 | |
| 1c5acb87 |
851 | } |
| 852 | |
| 49b5c49b |
853 | //________________________________________________ |
| 1c5acb87 |
854 | void AliCaloPID::Print(const Option_t * opt) const |
| 855 | { |
| 477d6cee |
856 | |
| 857 | //Print some relevant parameters set for the analysis |
| 858 | if(! opt) |
| 859 | return; |
| 860 | |
| 861 | printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ; |
| 862 | |
| 3c1d9afb |
863 | if(fUseBayesianWeights) |
| 864 | { |
| 49b5c49b |
865 | printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n", |
| 3c1d9afb |
866 | fPHOSPhotonWeight, fPHOSPi0Weight, |
| 867 | fPHOSElectronWeight, fPHOSChargeWeight, fPHOSNeutralWeight) ; |
| 49b5c49b |
868 | printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n", |
| 3c1d9afb |
869 | fEMCALPhotonWeight, fEMCALPi0Weight, |
| 870 | fEMCALElectronWeight, fEMCALChargeWeight, fEMCALNeutralWeight) ; |
| 49b5c49b |
871 | |
| 872 | printf("PHOS Parametrized weight on? = %d\n", fPHOSWeightFormula) ; |
| 3c1d9afb |
873 | if(fPHOSWeightFormula) |
| 874 | { |
| 49b5c49b |
875 | printf("Photon weight formula = %s\n", fPHOSPhotonWeightFormulaExpression.Data()); |
| 876 | printf("Pi0 weight formula = %s\n", fPHOSPi0WeightFormulaExpression .Data()); |
| 877 | } |
| 878 | if(fRecalculateBayesian) printf(" Recalculate bayesian with Particle Flux? = %d\n",fParticleFlux); |
| 879 | } |
| 3c1d9afb |
880 | else |
| 881 | { |
| 882 | printf("TOF cut = %e\n", fTOFCut); |
| 883 | printf("EMCAL Lambda0 cut min = %2.2f; max = %2.2f\n", fEMCALL0CutMin,fEMCALL0CutMax); |
| 884 | printf("EMCAL cluster-track dEta < %2.3f; dPhi < %2.3f\n", fEMCALDEtaCut, fEMCALDPhiCut); |
| 885 | printf("PHOS Treac matching cut =%2.2f, Dispersion Cut =%2.2f \n",fPHOSRCut, fPHOSDispersionCut) ; |
| 49b5c49b |
886 | |
| a5fb4114 |
887 | } |
| 477d6cee |
888 | |
| 3c1d9afb |
889 | if(fDoClusterSplitting) |
| 890 | { |
| 891 | printf("Min. N Cells =%d \n", fSplitMinNCells) ; |
| 892 | printf("%2.2f < lambda_0^2 <%2.2f \n",fSplitM02MinCut,fSplitM02MaxCut); |
| 893 | printf("pi0 : %2.2f<m<%2.2f \n", fMassPi0Min,fMassPi0Max); |
| 894 | printf("eta : %2.2f<m<%2.2f \n", fMassEtaMin,fMassEtaMax); |
| 895 | printf("phot: %2.2f<m<%2.2f \n", fMassPhoMin,fMassPhoMax); |
| 896 | } |
| 897 | |
| 477d6cee |
898 | printf(" \n"); |
| 899 | |
| 1c5acb87 |
900 | } |
| 901 | |
| 3c1d9afb |
902 | //_________________________________________________________________ |
| 903 | void AliCaloPID::PrintClusterPIDWeights(const Double_t * pid) const |
| 904 | { |
| 905 | // print PID of cluster, (AliVCluster*)cluster->GetPID() |
| 906 | |
| c2791479 |
907 | printf("AliCaloPID::PrintClusterPIDWeights() \n \t ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, \n \t \ |
| 3c1d9afb |
908 | pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n", |
| 909 | pid[AliVCluster::kPhoton], pid[AliVCluster::kPi0], |
| 910 | pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon], |
| 911 | pid[AliVCluster::kPion], pid[AliVCluster::kKaon], |
| 912 | pid[AliVCluster::kProton], |
| 913 | pid[AliVCluster::kNeutron], pid[AliVCluster::kKaon0]); |
| 914 | |
| 915 | } |
| 916 | |
| 49b5c49b |
917 | //___________________________________________________________________________ |
| 3c1d9afb |
918 | void AliCaloPID::SetPIDBits(AliVCluster * cluster, |
| 49b5c49b |
919 | AliAODPWG4Particle * ph, AliCalorimeterUtils* cu, |
| 920 | AliVEvent* event) |
| 921 | { |
| 477d6cee |
922 | //Set Bits for PID selection |
| 923 | |
| 924 | //Dispersion/lambdas |
| 5ae09196 |
925 | //Double_t disp= cluster->GetDispersion() ; |
| 926 | Double_t l1 = cluster->GetM20() ; |
| 927 | Double_t l0 = cluster->GetM02() ; |
| 928 | Bool_t isDispOK = kTRUE ; |
| 9a6fa057 |
929 | if(cluster->IsPHOS()){ |
| 49b5c49b |
930 | if(TestPHOSDispersion(ph->Pt(),l0,l1) < fPHOSDispersionCut) isDispOK = kTRUE; |
| 931 | else isDispOK = kFALSE; |
| 5ae09196 |
932 | } |
| 933 | else{//EMCAL |
| 934 | |
| 49b5c49b |
935 | if(l0 > fEMCALL0CutMin && l0 < fEMCALL0CutMax) isDispOK = kTRUE; |
| 936 | |
| 5ae09196 |
937 | } |
| 938 | |
| 939 | ph->SetDispBit(isDispOK) ; |
| 477d6cee |
940 | |
| 941 | //TOF |
| 942 | Double_t tof=cluster->GetTOF() ; |
| 943 | ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ; |
| 944 | |
| 49b5c49b |
945 | //Charged |
| 946 | Bool_t isNeutral = IsTrackMatched(cluster,cu,event); |
| 5ae09196 |
947 | |
| 948 | ph->SetChargedBit(isNeutral); |
| 477d6cee |
949 | |
| 950 | //Set PID pdg |
| 3c1d9afb |
951 | ph->SetIdentifiedParticleType(GetIdentifiedParticleType(cluster)); |
| 477d6cee |
952 | |
| 5a72d9af |
953 | if(fDebug > 0) |
| 954 | { |
| 5ae09196 |
955 | printf("AliCaloPID::SetPIDBits: TOF %e, Lambda0 %2.2f, Lambda1 %2.2f\n",tof , l0, l1); |
| 477d6cee |
956 | printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n", |
| 49b5c49b |
957 | ph->GetIdentifiedParticleType(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit()); |
| 477d6cee |
958 | } |
| 1c5acb87 |
959 | } |
| 960 | |
| 09273901 |
961 | //_________________________________________________________ |
| 49b5c49b |
962 | Bool_t AliCaloPID::IsTrackMatched(AliVCluster* cluster, |
| 963 | AliCalorimeterUtils * cu, |
| 964 | AliVEvent* event) const |
| 965 | { |
| 5ae09196 |
966 | //Check if there is any track attached to this cluster |
| 967 | |
| 968 | Int_t nMatches = cluster->GetNTracksMatched(); |
| 49b5c49b |
969 | AliVTrack * track = 0; |
| 970 | Double_t p[3]; |
| 971 | |
| 44443bbd |
972 | if(nMatches > 0) |
| 973 | { |
| 49b5c49b |
974 | //In case of ESDs, by default without match one entry with negative index, no match, reject. |
| 975 | if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName()))) |
| 976 | { |
| 977 | Int_t iESDtrack = cluster->GetTrackMatchedIndex(); |
| 978 | if(iESDtrack >= 0) track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack)); |
| 979 | else return kFALSE; |
| d39cba7e |
980 | |
| 44443bbd |
981 | if (!track) |
| 982 | { |
| 983 | if(fDebug > 0) printf("AliCaloPID::IsTrackMatched() - Null matched track in ESD when index is OK!\n"); |
| 49b5c49b |
984 | return kFALSE; |
| 985 | } |
| 986 | } |
| 987 | else { // AOD |
| 988 | track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0)); |
| 44443bbd |
989 | if (!track) |
| 990 | { |
| 991 | if(fDebug > 0) printf("AliCaloPID::IsTrackMatched() - Null matched track in AOD!\n"); |
| 49b5c49b |
992 | return kFALSE; |
| c76f0089 |
993 | } |
| d39cba7e |
994 | } |
| 5ae09196 |
995 | |
| 49b5c49b |
996 | Float_t dZ = cluster->GetTrackDz(); |
| 997 | Float_t dR = cluster->GetTrackDx(); |
| 998 | |
| 999 | // if track matching was recalculated |
| 44443bbd |
1000 | if(cluster->IsEMCAL() && cu && cu->IsRecalculationOfClusterTrackMatchingOn()) |
| 1001 | { |
| 49b5c49b |
1002 | dR = 2000., dZ = 2000.; |
| 31ae6d59 |
1003 | cu->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR); |
| 49b5c49b |
1004 | } |
| 09273901 |
1005 | |
| 44443bbd |
1006 | if(cluster->IsPHOS()) |
| 1007 | { |
| 49b5c49b |
1008 | track->GetPxPyPz(p) ; |
| 1009 | TLorentzVector trackmom(p[0],p[1],p[2],0); |
| 1010 | Int_t charge = track->Charge(); |
| 1011 | Double_t mf = event->GetMagneticField(); |
| 1012 | if(TestPHOSChargedVeto(dR, dZ, trackmom.Pt(), charge, mf ) < fPHOSRCut) return kTRUE; |
| 1013 | else return kFALSE; |
| 1014 | |
| 1015 | }//PHOS |
| 5a72d9af |
1016 | else //EMCAL |
| 1017 | { |
| 1018 | if(fDebug > 1) |
| 49b5c49b |
1019 | printf("AliCaloPID::IsTrackMatched - EMCAL dR %f < %f, dZ %f < %f \n",dR, fEMCALDPhiCut, dZ, fEMCALDEtaCut); |
| 1020 | |
| 1021 | if(TMath::Abs(dR) < fEMCALDPhiCut && |
| 1022 | TMath::Abs(dZ) < fEMCALDEtaCut) return kTRUE; |
| 1023 | else return kFALSE; |
| 1024 | |
| 1025 | }//EMCAL cluster |
| 1026 | |
| 1027 | |
| 1028 | } // more than 1 match, at least one track in array |
| 1029 | else return kFALSE; |
| 1030 | |
| 1031 | } |
| 1032 | |
| 1033 | //___________________________________________________________________________________________________ |
| 1034 | Float_t AliCaloPID::TestPHOSDispersion(const Double_t pt, const Double_t l1, const Double_t l2) const |
| 1035 | { |
| 1036 | //Check if cluster photon-like. Uses photon cluster parameterization in real pp data |
| 1037 | //Returns distance in sigmas. Recommended cut 2.5 |
| 1038 | |
| 1039 | Double_t l2Mean = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ; |
| 1040 | Double_t l1Mean = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ; |
| 1041 | Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt; |
| 1042 | Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt; |
| 1043 | Double_t c =-0.35-0.550*TMath::Exp(-0.390730*pt) ; |
| f3138ecf |
1044 | Double_t r2 = 0.5* (l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma + |
| 49b5c49b |
1045 | 0.5* (l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma + |
| 1046 | 0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ; |
| 1047 | |
| f3138ecf |
1048 | if(fDebug > 0) printf("AliCaloPID::TestPHOSDispersion() - PHOS SS R %f < %f?\n", TMath::Sqrt(r2), fPHOSDispersionCut); |
| 5ae09196 |
1049 | |
| f3138ecf |
1050 | return TMath::Sqrt(r2) ; |
| 5ae09196 |
1051 | |
| 1052 | } |
| 1053 | |
| 49b5c49b |
1054 | //_______________________________________________________________________________________________ |
| 1055 | Float_t AliCaloPID::TestPHOSChargedVeto(const Double_t dx, const Double_t dz, const Double_t pt, |
| 1056 | const Int_t charge, const Double_t mf) const |
| 1057 | { |
| 1058 | //Checks distance to the closest track. Takes into account |
| 1059 | //non-perpendicular incidence of tracks. |
| 1060 | //returns distance in sigmas. Recommended cut: 2. |
| 1061 | //Requires (sign) of magnetic filed. onc can find it for example as following |
| 1062 | // Double_t mf=0. ; |
| 1063 | // AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent()); |
| 1064 | // if(event) |
| 1065 | // mf = event->GetMagneticField(); //Positive for ++ and negative for -- |
| 1066 | |
| 1067 | |
| 1068 | Double_t meanX = 0.; |
| 1069 | Double_t meanZ = 0.; |
| 1070 | Double_t sx = TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+ |
| 1071 | 6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+ |
| 1072 | 1.59219); |
| 1073 | Double_t sz = TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+ |
| 1074 | 1.60) ; |
| 1075 | |
| 1076 | if(mf<0.){ //field -- |
| 1077 | meanZ = -0.468318 ; |
| 1078 | if(charge>0) |
| 1079 | meanX = TMath::Min(7.3, 3.89994*1.20679 *1.20679 /(pt*pt+1.20679*1.20679)+ |
| 1080 | 0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ; |
| 1081 | else |
| 1082 | meanX =-TMath::Min(7.7, 3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+ |
| 1083 | 1.23114 +4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ; |
| 1084 | } |
| 1085 | else{ //Field ++ |
| 1086 | meanZ = -0.468318; |
| 1087 | if(charge>0) |
| 1088 | meanX =-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+ |
| 1089 | 0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ; |
| 1090 | else |
| 1091 | meanX = TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)- |
| 1092 | 0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ; |
| 1093 | } |
| 9a6fa057 |
1094 | |
| 49b5c49b |
1095 | Double_t rz = (dz-meanZ)/sz ; |
| 1096 | Double_t rx = (dx-meanX)/sx ; |
| 9a6fa057 |
1097 | |
| 49b5c49b |
1098 | if(fDebug > 0) |
| 1099 | printf("AliCaloPID::TestPHOSDispersion() - PHOS Matching R %f < %f\n",TMath::Sqrt(rx*rx+rz*rz), fPHOSRCut); |
| 9a6fa057 |
1100 | |
| 49b5c49b |
1101 | return TMath::Sqrt(rx*rx+rz*rz) ; |
| 9a6fa057 |
1102 | |
| 1103 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame