]>
Commit | Line | Data |
---|---|---|
10d100d4 | 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 | **************************************************************************/ | |
15 | ||
16 | //----------------------------------------------------------------- | |
17 | // Implementation of the TPC PID class | |
18 | // Very naive one... Should be made better by the detector experts... | |
19 | // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch | |
20 | // With many additions and modifications suggested by | |
21 | // Alexander Kalweit, GSI, alexander.philipp.kalweit@cern.ch | |
22 | // Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de | |
644666df | 23 | // ...and some modifications by |
24 | // Mikolaj Krzewicki, GSI, mikolaj.krzewicki@cern.ch | |
f84b18dd | 25 | // ...and some modifications plus eta correction functions by |
26 | // Benjamin Hess, University of Tuebingen, bhess@cern.ch | |
10d100d4 | 27 | //----------------------------------------------------------------- |
28 | ||
d2aa6df0 | 29 | #include <TGraph.h> |
30 | #include <TObjArray.h> | |
31 | #include <TSpline.h> | |
644666df | 32 | #include <TBits.h> |
00a38d07 | 33 | #include <TMath.h> |
f84b18dd | 34 | #include <TH2D.h> |
d2aa6df0 | 35 | |
644666df | 36 | #include <AliLog.h> |
10d100d4 | 37 | #include "AliExternalTrackParam.h" |
644666df | 38 | #include "AliVTrack.h" |
d2aa6df0 | 39 | #include "AliTPCPIDResponse.h" |
644666df | 40 | #include "AliTPCdEdxInfo.h" |
912f7fdc | 41 | #include "TFile.h" |
42 | #include "TSpline.h" | |
644666df | 43 | |
44 | ClassImp(AliTPCPIDResponse); | |
d2aa6df0 | 45 | |
912f7fdc | 46 | |
47 | AliTPCPIDResponse *AliTPCPIDResponse::fgInstance =0; | |
48 | ||
644666df | 49 | const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]= |
50 | { | |
51 | "", //default - no name | |
52 | "1", //all high | |
53 | "2", //OROC only | |
54 | "unknown" | |
55 | }; | |
10d100d4 | 56 | |
57 | //_________________________________________________________________________ | |
58 | AliTPCPIDResponse::AliTPCPIDResponse(): | |
644666df | 59 | TNamed(), |
d2aa6df0 | 60 | fMIP(50.), |
644666df | 61 | fRes0(), |
62 | fResN2(), | |
d2aa6df0 | 63 | fKp1(0.0283086), |
64 | fKp2(2.63394e+01), | |
65 | fKp3(5.04114e-11), | |
66 | fKp4(2.12543), | |
67 | fKp5(4.88663), | |
68 | fUseDatabase(kFALSE), | |
644666df | 69 | fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios), |
70 | fVoltageMap(72), | |
71 | fLowGainIROCthreshold(-40), | |
72 | fBadIROCthreshhold(-70), | |
73 | fLowGainOROCthreshold(-40), | |
74 | fBadOROCthreshhold(-40), | |
75 | fMaxBadLengthFraction(0.5), | |
f84b18dd | 76 | fMagField(0.), |
77 | fhEtaCorr(0x0), | |
78 | fhEtaSigmaPar1(0x0), | |
87da0205 | 79 | fSigmaPar0(0.0), |
80 | fCurrentEventMultiplicity(0), | |
81 | fCorrFuncMultiplicity(0x0), | |
82 | fCorrFuncMultiplicityTanTheta(0x0), | |
83 | fCorrFuncSigmaMultiplicity(0x0) | |
10d100d4 | 84 | { |
85 | // | |
86 | // The default constructor | |
87 | // | |
5a9dc560 | 88 | |
89 | AliLog::SetClassDebugLevel("AliTPCPIDResponse", AliLog::kInfo); | |
90 | ||
644666df | 91 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;} |
87da0205 | 92 | |
93 | fCorrFuncMultiplicity = new TF1("fCorrFuncMultiplicity", | |
94 | "[0] + [1]*TMath::Max([4], TMath::Min(x, [3])) + [2] * TMath::Power(TMath::Max([4], TMath::Min(x, [3])), 2)", | |
95 | 0., 0.2); | |
96 | fCorrFuncMultiplicityTanTheta = new TF1("fCorrFuncMultiplicityTanTheta", "[0] * (x - [2]) + [1] * (x * x - [2] * [2])", -1.5, 1.5); | |
97 | fCorrFuncSigmaMultiplicity = new TF1("fCorrFuncSigmaMultiplicity", | |
98 | "TMath::Max(0, [0] + [1]*TMath::Min(x, [3]) + [2] * TMath::Power(TMath::Min(x, [3]), 2))", 0., 0.2); | |
99 | ||
100 | ||
101 | ResetMultiplicityCorrectionFunctions(); | |
912f7fdc | 102 | fgInstance=this; |
10d100d4 | 103 | } |
1b45e564 | 104 | /*TODO remove? |
10d100d4 | 105 | //_________________________________________________________________________ |
d2aa6df0 | 106 | AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param): |
644666df | 107 | TNamed(), |
d2aa6df0 | 108 | fMIP(param[0]), |
644666df | 109 | fRes0(), |
110 | fResN2(), | |
d2aa6df0 | 111 | fKp1(0.0283086), |
112 | fKp2(2.63394e+01), | |
113 | fKp3(5.04114e-11), | |
114 | fKp4(2.12543), | |
115 | fKp5(4.88663), | |
116 | fUseDatabase(kFALSE), | |
644666df | 117 | fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios), |
118 | fVoltageMap(72), | |
119 | fLowGainIROCthreshold(-40), | |
120 | fBadIROCthreshhold(-70), | |
121 | fLowGainOROCthreshold(-40), | |
122 | fBadOROCthreshhold(-40), | |
123 | fMaxBadLengthFraction(0.5), | |
f84b18dd | 124 | fMagField(0.), |
125 | fhEtaCorr(0x0), | |
126 | fhEtaSigmaPar1(0x0), | |
127 | fSigmaPar0(0.0) | |
10d100d4 | 128 | { |
129 | // | |
130 | // The main constructor | |
131 | // | |
644666df | 132 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];} |
133 | } | |
1b45e564 | 134 | */ |
f84b18dd | 135 | |
136 | //_________________________________________________________________________ | |
137 | AliTPCPIDResponse::~AliTPCPIDResponse() | |
138 | { | |
139 | // | |
140 | // Destructor | |
141 | // | |
142 | ||
143 | delete fhEtaCorr; | |
144 | fhEtaCorr = 0x0; | |
145 | ||
146 | delete fhEtaSigmaPar1; | |
147 | fhEtaSigmaPar1 = 0x0; | |
87da0205 | 148 | |
149 | delete fCorrFuncMultiplicity; | |
150 | fCorrFuncMultiplicity = 0x0; | |
151 | ||
152 | delete fCorrFuncMultiplicityTanTheta; | |
153 | fCorrFuncMultiplicityTanTheta = 0x0; | |
154 | ||
155 | delete fCorrFuncSigmaMultiplicity; | |
156 | fCorrFuncSigmaMultiplicity = 0x0; | |
912f7fdc | 157 | if (fgInstance==this) fgInstance=0; |
f84b18dd | 158 | } |
159 | ||
160 | ||
644666df | 161 | //_________________________________________________________________________ |
162 | AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that): | |
163 | TNamed(that), | |
164 | fMIP(that.fMIP), | |
165 | fRes0(), | |
166 | fResN2(), | |
167 | fKp1(that.fKp1), | |
168 | fKp2(that.fKp2), | |
169 | fKp3(that.fKp3), | |
170 | fKp4(that.fKp4), | |
171 | fKp5(that.fKp5), | |
172 | fUseDatabase(that.fUseDatabase), | |
173 | fResponseFunctions(that.fResponseFunctions), | |
174 | fVoltageMap(that.fVoltageMap), | |
175 | fLowGainIROCthreshold(that.fLowGainIROCthreshold), | |
176 | fBadIROCthreshhold(that.fBadIROCthreshhold), | |
177 | fLowGainOROCthreshold(that.fLowGainOROCthreshold), | |
178 | fBadOROCthreshhold(that.fBadOROCthreshhold), | |
179 | fMaxBadLengthFraction(that.fMaxBadLengthFraction), | |
f84b18dd | 180 | fMagField(that.fMagField), |
181 | fhEtaCorr(0x0), | |
182 | fhEtaSigmaPar1(0x0), | |
87da0205 | 183 | fSigmaPar0(that.fSigmaPar0), |
184 | fCurrentEventMultiplicity(that.fCurrentEventMultiplicity), | |
185 | fCorrFuncMultiplicity(0x0), | |
186 | fCorrFuncMultiplicityTanTheta(0x0), | |
187 | fCorrFuncSigmaMultiplicity(0x0) | |
644666df | 188 | { |
189 | //copy ctor | |
190 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];} | |
f84b18dd | 191 | |
192 | // Copy eta maps | |
1b45e564 | 193 | if (that.fhEtaCorr) { |
3800f96f | 194 | fhEtaCorr = new TH2D(*(that.fhEtaCorr)); |
195 | fhEtaCorr->SetDirectory(0); | |
196 | } | |
f84b18dd | 197 | |
1b45e564 | 198 | if (that.fhEtaSigmaPar1) { |
3800f96f | 199 | fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1)); |
200 | fhEtaSigmaPar1->SetDirectory(0); | |
201 | } | |
87da0205 | 202 | |
203 | // Copy multiplicity correction functions | |
204 | if (that.fCorrFuncMultiplicity) { | |
205 | fCorrFuncMultiplicity = new TF1(*(that.fCorrFuncMultiplicity)); | |
206 | } | |
207 | ||
208 | if (that.fCorrFuncMultiplicityTanTheta) { | |
209 | fCorrFuncMultiplicityTanTheta = new TF1(*(that.fCorrFuncMultiplicityTanTheta)); | |
210 | } | |
211 | ||
212 | if (that.fCorrFuncSigmaMultiplicity) { | |
213 | fCorrFuncSigmaMultiplicity = new TF1(*(that.fCorrFuncSigmaMultiplicity)); | |
214 | } | |
644666df | 215 | } |
216 | ||
217 | //_________________________________________________________________________ | |
218 | AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that) | |
219 | { | |
220 | //assignment | |
221 | if (&that==this) return *this; | |
222 | TNamed::operator=(that); | |
223 | fMIP=that.fMIP; | |
224 | fKp1=that.fKp1; | |
225 | fKp2=that.fKp2; | |
226 | fKp3=that.fKp3; | |
227 | fKp4=that.fKp4; | |
228 | fKp5=that.fKp5; | |
229 | fUseDatabase=that.fUseDatabase; | |
230 | fResponseFunctions=that.fResponseFunctions; | |
231 | fVoltageMap=that.fVoltageMap; | |
232 | fLowGainIROCthreshold=that.fLowGainIROCthreshold; | |
233 | fBadIROCthreshhold=that.fBadIROCthreshhold; | |
234 | fLowGainOROCthreshold=that.fLowGainOROCthreshold; | |
235 | fBadOROCthreshhold=that.fBadOROCthreshhold; | |
236 | fMaxBadLengthFraction=that.fMaxBadLengthFraction; | |
237 | fMagField=that.fMagField; | |
87da0205 | 238 | fCurrentEventMultiplicity=that.fCurrentEventMultiplicity; |
644666df | 239 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];} |
f84b18dd | 240 | |
241 | delete fhEtaCorr; | |
f85a3764 | 242 | fhEtaCorr=0x0; |
1b45e564 | 243 | if (that.fhEtaCorr) { |
3800f96f | 244 | fhEtaCorr = new TH2D(*(that.fhEtaCorr)); |
245 | fhEtaCorr->SetDirectory(0); | |
246 | } | |
f84b18dd | 247 | |
248 | delete fhEtaSigmaPar1; | |
f85a3764 | 249 | fhEtaSigmaPar1=0x0; |
1b45e564 | 250 | if (that.fhEtaSigmaPar1) { |
3800f96f | 251 | fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1)); |
252 | fhEtaSigmaPar1->SetDirectory(0); | |
253 | } | |
f84b18dd | 254 | |
255 | fSigmaPar0 = that.fSigmaPar0; | |
87da0205 | 256 | |
257 | delete fCorrFuncMultiplicity; | |
258 | fCorrFuncMultiplicity = 0x0; | |
259 | if (that.fCorrFuncMultiplicity) { | |
260 | fCorrFuncMultiplicity = new TF1(*(that.fCorrFuncMultiplicity)); | |
261 | } | |
262 | ||
263 | delete fCorrFuncMultiplicityTanTheta; | |
264 | fCorrFuncMultiplicityTanTheta = 0x0; | |
265 | if (that.fCorrFuncMultiplicityTanTheta) { | |
266 | fCorrFuncMultiplicityTanTheta = new TF1(*(that.fCorrFuncMultiplicityTanTheta)); | |
267 | } | |
268 | ||
269 | delete fCorrFuncSigmaMultiplicity; | |
270 | fCorrFuncSigmaMultiplicity = 0x0; | |
271 | if (that.fCorrFuncSigmaMultiplicity) { | |
272 | fCorrFuncSigmaMultiplicity = new TF1(*(that.fCorrFuncSigmaMultiplicity)); | |
273 | } | |
f84b18dd | 274 | |
644666df | 275 | return *this; |
10d100d4 | 276 | } |
277 | ||
d2aa6df0 | 278 | //_________________________________________________________________________ |
10d100d4 | 279 | Double_t AliTPCPIDResponse::Bethe(Double_t betaGamma) const { |
280 | // | |
281 | // This is the Bethe-Bloch function normalised to 1 at the minimum | |
282 | // WARNING | |
283 | // Simulated and reconstructed Bethe-Bloch differs | |
284 | // Simulated curve is the dNprim/dx | |
285 | // Reconstructed is proportianal dNtot/dx | |
286 | // Temporary fix for production - Simple linear correction function | |
287 | // Future 2 Bethe Bloch formulas needed | |
288 | // 1. for simulation | |
289 | // 2. for reconstructed PID | |
290 | // | |
d2aa6df0 | 291 | |
292 | // const Float_t kmeanCorrection =0.1; | |
10d100d4 | 293 | Double_t bb= |
294 | AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5); | |
295 | return bb*fMIP; | |
296 | } | |
297 | ||
298 | //_________________________________________________________________________ | |
299 | void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1, | |
300 | Double_t kp2, | |
301 | Double_t kp3, | |
302 | Double_t kp4, | |
303 | Double_t kp5) { | |
304 | // | |
305 | // Set the parameters of the ALEPH Bethe-Bloch formula | |
306 | // | |
307 | fKp1=kp1; | |
308 | fKp2=kp2; | |
309 | fKp3=kp3; | |
310 | fKp4=kp4; | |
311 | fKp5=kp5; | |
312 | } | |
644666df | 313 | |
10d100d4 | 314 | //_________________________________________________________________________ |
315 | void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) { | |
316 | // | |
317 | // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint) | |
318 | // | |
644666df | 319 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;} |
10d100d4 | 320 | } |
321 | ||
322 | //_________________________________________________________________________ | |
323 | Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom, | |
d2aa6df0 | 324 | AliPID::EParticleType n) const { |
10d100d4 | 325 | // |
f84b18dd | 326 | // Deprecated function (for backward compatibility). Please use |
327 | // GetExpectedSignal(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource, | |
87da0205 | 328 | // Bool_t correctEta, Bool_t correctMultiplicity); |
f84b18dd | 329 | // instead! |
330 | // | |
331 | // | |
10d100d4 | 332 | // Calculates the expected PID signal as the function of |
333 | // the information stored in the track, for the specified particle type | |
334 | // | |
335 | // At the moment, these signals are just the results of calling the | |
336 | // Bethe-Bloch formula. | |
337 | // This can be improved. By taking into account the number of | |
338 | // assigned clusters and/or the track dip angle, for example. | |
339 | // | |
340 | ||
f84b18dd | 341 | //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...) |
342 | // !!! Splines for light nuclei need to be normalised to this factor !!! | |
343 | const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3); | |
344 | ||
539a5a59 | 345 | Double_t mass=AliPID::ParticleMassZ(n); |
f84b18dd | 346 | if (!fUseDatabase) return Bethe(mom/mass) * chargeFactor; |
d2aa6df0 | 347 | // |
644666df | 348 | const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n); |
00a38d07 | 349 | |
f84b18dd | 350 | if (!responseFunction) return Bethe(mom/mass) * chargeFactor; |
351 | ||
00a38d07 | 352 | return fMIP*responseFunction->Eval(mom/mass)*chargeFactor; |
d2aa6df0 | 353 | |
10d100d4 | 354 | } |
355 | ||
356 | //_________________________________________________________________________ | |
357 | Double_t AliTPCPIDResponse::GetExpectedSigma(const Float_t mom, | |
f84b18dd | 358 | const Int_t nPoints, |
644666df | 359 | AliPID::EParticleType n) const { |
10d100d4 | 360 | // |
f84b18dd | 361 | // Deprecated function (for backward compatibility). Please use |
362 | // GetExpectedSigma(onst AliVTrack* track, AliPID::EParticleType species, | |
363 | // ETPCdEdxSource dedxSource, Bool_t correctEta) instead! | |
364 | // | |
365 | // | |
10d100d4 | 366 | // Calculates the expected sigma of the PID signal as the function of |
367 | // the information stored in the track, for the specified particle type | |
368 | // | |
644666df | 369 | |
370 | if (nPoints != 0) | |
371 | return GetExpectedSignal(mom,n)*fRes0[0]*sqrt(1. + fResN2[0]/nPoints); | |
372 | else | |
373 | return GetExpectedSignal(mom,n)*fRes0[0]; | |
374 | } | |
375 | ||
376 | ////////////////////////////////////////////////////NEW////////////////////////////// | |
377 | ||
378 | //_________________________________________________________________________ | |
379 | void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario) { | |
380 | // | |
381 | // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint) | |
382 | // | |
bd46387d | 383 | if ((Int_t)gainScenario>=(Int_t)fgkNumberOfGainScenarios) return; //TODO: better handling! |
644666df | 384 | fRes0[gainScenario]=res0; |
385 | fResN2[gainScenario]=resN2; | |
386 | } | |
387 | ||
f84b18dd | 388 | |
644666df | 389 | //_________________________________________________________________________ |
f84b18dd | 390 | Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, |
644666df | 391 | AliPID::EParticleType species, |
f84b18dd | 392 | Double_t /*dEdx*/, |
393 | const TSpline3* responseFunction, | |
87da0205 | 394 | Bool_t correctEta, |
395 | Bool_t correctMultiplicity) const | |
644666df | 396 | { |
397 | // Calculates the expected PID signal as the function of | |
f84b18dd | 398 | // the information stored in the track and the given parameters, |
399 | // for the specified particle type | |
644666df | 400 | // |
401 | // At the moment, these signals are just the results of calling the | |
87da0205 | 402 | // Bethe-Bloch formula plus, if desired, taking into account the eta dependence |
403 | // and the multiplicity dependence (for PbPb). | |
644666df | 404 | // This can be improved. By taking into account the number of |
405 | // assigned clusters and/or the track dip angle, for example. | |
10d100d4 | 406 | // |
407 | ||
f84b18dd | 408 | Double_t mom=track->GetTPCmomentum(); |
97d17530 | 409 | Double_t mass=AliPID::ParticleMassZ(species); |
644666df | 410 | |
f84b18dd | 411 | //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...) |
412 | // !!! Splines for light nuclei need to be normalised to this factor !!! | |
413 | const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3); | |
414 | ||
415 | if (!responseFunction) | |
416 | return Bethe(mom/mass) * chargeFactor; | |
644666df | 417 | |
f84b18dd | 418 | Double_t dEdxSplines = fMIP*responseFunction->Eval(mom/mass) * chargeFactor; |
87da0205 | 419 | |
420 | if (!correctEta && !correctMultiplicity) | |
f84b18dd | 421 | return dEdxSplines; |
422 | ||
87da0205 | 423 | Double_t corrFactorEta = 1.0; |
424 | Double_t corrFactorMultiplicity = 1.0; | |
425 | ||
426 | if (correctEta) { | |
ef7661fd | 427 | corrFactorEta = GetEtaCorrectionFast(track, dEdxSplines); |
87da0205 | 428 | //TODO Alternatively take current track dEdx |
ef7661fd | 429 | //corrFactorEta = GetEtaCorrectionFast(track, dEdx); |
87da0205 | 430 | } |
431 | ||
432 | if (correctMultiplicity) | |
ef7661fd | 433 | corrFactorMultiplicity = GetMultiplicityCorrectionFast(track, dEdxSplines * corrFactorEta, fCurrentEventMultiplicity); |
87da0205 | 434 | |
435 | return dEdxSplines * corrFactorEta * corrFactorMultiplicity; | |
644666df | 436 | } |
437 | ||
f84b18dd | 438 | |
644666df | 439 | //_________________________________________________________________________ |
440 | Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, | |
441 | AliPID::EParticleType species, | |
f84b18dd | 442 | ETPCdEdxSource dedxSource, |
87da0205 | 443 | Bool_t correctEta, |
444 | Bool_t correctMultiplicity) const | |
644666df | 445 | { |
446 | // Calculates the expected PID signal as the function of | |
447 | // the information stored in the track, for the specified particle type | |
448 | // | |
449 | // At the moment, these signals are just the results of calling the | |
87da0205 | 450 | // Bethe-Bloch formula plus, if desired, taking into account the eta dependence |
451 | // and the multiplicity dependence (for PbPb). | |
644666df | 452 | // This can be improved. By taking into account the number of |
453 | // assigned clusters and/or the track dip angle, for example. | |
454 | // | |
455 | ||
f84b18dd | 456 | if (!fUseDatabase) { |
457 | //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...) | |
458 | // !!! Splines for light nuclei need to be normalised to this factor !!! | |
459 | const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3); | |
644666df | 460 | |
97d17530 | 461 | return Bethe(track->GetTPCmomentum() / AliPID::ParticleMassZ(species)) * chargeFactor; |
f84b18dd | 462 | } |
644666df | 463 | |
f84b18dd | 464 | Double_t dEdx = -1; |
465 | Int_t nPoints = -1; | |
466 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
467 | TSpline3* responseFunction = 0x0; | |
468 | ||
469 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) { | |
470 | // Something is wrong with the track -> Return obviously invalid value | |
471 | return -999; | |
472 | } | |
473 | ||
474 | // Charge factor already taken into account inside the following function call | |
87da0205 | 475 | return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity); |
644666df | 476 | } |
477 | ||
478 | //_________________________________________________________________________ | |
479 | TSpline3* AliTPCPIDResponse::GetResponseFunction( AliPID::EParticleType type, | |
480 | AliTPCPIDResponse::ETPCgainScenario gainScenario ) const | |
481 | { | |
482 | //get response function | |
483 | return dynamic_cast<TSpline3*>(fResponseFunctions.At(ResponseFunctionIndex(type,gainScenario))); | |
484 | } | |
485 | ||
486 | //_________________________________________________________________________ | |
487 | TSpline3* AliTPCPIDResponse::GetResponseFunction( const AliVTrack* track, | |
488 | AliPID::EParticleType species, | |
f84b18dd | 489 | ETPCdEdxSource dedxSource) const |
644666df | 490 | { |
491 | //the splines are stored in an array, different scenarios | |
492 | ||
f84b18dd | 493 | Double_t dEdx = -1; |
494 | Int_t nPoints = -1; | |
495 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
496 | TSpline3* responseFunction = 0x0; | |
497 | ||
498 | if (ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
499 | return responseFunction; | |
644666df | 500 | |
501 | return NULL; | |
502 | } | |
503 | ||
504 | //_________________________________________________________________________ | |
505 | void AliTPCPIDResponse::ResetSplines() | |
506 | { | |
507 | //reset the array with splines | |
508 | for (Int_t i=0;i<fResponseFunctions.GetEntriesFast();i++) | |
509 | { | |
510 | fResponseFunctions.AddAt(NULL,i); | |
511 | } | |
512 | } | |
513 | //_________________________________________________________________________ | |
514 | Int_t AliTPCPIDResponse::ResponseFunctionIndex( AliPID::EParticleType species, | |
515 | ETPCgainScenario gainScenario ) const | |
516 | { | |
517 | //get the index in fResponseFunctions given type and scenario | |
518 | return Int_t(species)+Int_t(gainScenario)*fgkNumberOfParticleSpecies; | |
519 | } | |
520 | ||
521 | //_________________________________________________________________________ | |
522 | void AliTPCPIDResponse::SetResponseFunction( TObject* o, | |
523 | AliPID::EParticleType species, | |
524 | ETPCgainScenario gainScenario ) | |
525 | { | |
526 | fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario)); | |
527 | } | |
528 | ||
f84b18dd | 529 | |
644666df | 530 | //_________________________________________________________________________ |
f84b18dd | 531 | Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, |
532 | AliPID::EParticleType species, | |
533 | ETPCgainScenario gainScenario, | |
534 | Double_t dEdx, | |
535 | Int_t nPoints, | |
536 | const TSpline3* responseFunction, | |
87da0205 | 537 | Bool_t correctEta, |
538 | Bool_t correctMultiplicity) const | |
644666df | 539 | { |
540 | // Calculates the expected sigma of the PID signal as the function of | |
f84b18dd | 541 | // the information stored in the track and the given parameters, |
542 | // for the specified particle type | |
644666df | 543 | // |
544 | ||
f84b18dd | 545 | if (!responseFunction) |
546 | return 999; | |
87da0205 | 547 | |
548 | //TODO Check whether it makes sense to set correctMultiplicity to kTRUE while correctEta might be kFALSE | |
5a9dc560 | 549 | |
550 | // If eta correction (=> new sigma parametrisation) is requested, but no sigma map is available, print error message | |
551 | if (correctEta && !fhEtaSigmaPar1) { | |
552 | AliError("New sigma parametrisation requested, but sigma map not initialised (usually via AliPIDResponse). Old sigma parametrisation will be used!"); | |
553 | } | |
554 | ||
f84b18dd | 555 | // If no sigma map is available or if no eta correction is requested (sigma maps only for corrected eta!), use the old parametrisation |
556 | if (!fhEtaSigmaPar1 || !correctEta) { | |
557 | if (nPoints != 0) | |
87da0205 | 558 | return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity) * |
f84b18dd | 559 | fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints); |
560 | else | |
87da0205 | 561 | return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity)*fRes0[gainScenario]; |
f84b18dd | 562 | } |
563 | ||
564 | if (nPoints > 0) { | |
87da0205 | 565 | // Use eta correction (+ eta-dependent sigma) |
ef7661fd | 566 | Double_t sigmaPar1 = GetSigmaPar1Fast(track, species, dEdx, responseFunction); |
f84b18dd | 567 | |
87da0205 | 568 | if (correctMultiplicity) { |
569 | // In addition, take into account multiplicity dependence of mean and sigma of dEdx | |
570 | Double_t dEdxExpectedEtaCorrected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE); | |
571 | ||
572 | // GetMultiplicityCorrection and GetMultiplicitySigmaCorrection both need the eta corrected dEdxExpected | |
ef7661fd | 573 | Double_t multiplicityCorrFactor = GetMultiplicityCorrectionFast(track, dEdxExpectedEtaCorrected, fCurrentEventMultiplicity); |
574 | Double_t multiplicitySigmaCorrFactor = GetMultiplicitySigmaCorrectionFast(dEdxExpectedEtaCorrected, fCurrentEventMultiplicity); | |
87da0205 | 575 | |
576 | // multiplicityCorrFactor to correct dEdxExpected for multiplicity. In addition: Correction factor for sigma | |
577 | return (dEdxExpectedEtaCorrected * multiplicityCorrFactor) | |
578 | * (sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints) * multiplicitySigmaCorrFactor); | |
579 | } | |
580 | else { | |
581 | return GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE)* | |
582 | sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints); | |
583 | } | |
f84b18dd | 584 | } |
585 | else { | |
586 | // One should never have/take tracks with 0 dEdx clusters! | |
587 | return 999; | |
588 | } | |
644666df | 589 | } |
590 | ||
f84b18dd | 591 | |
644666df | 592 | //_________________________________________________________________________ |
593 | Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, | |
594 | AliPID::EParticleType species, | |
f84b18dd | 595 | ETPCdEdxSource dedxSource, |
87da0205 | 596 | Bool_t correctEta, |
597 | Bool_t correctMultiplicity) const | |
644666df | 598 | { |
599 | // Calculates the expected sigma of the PID signal as the function of | |
600 | // the information stored in the track, for the specified particle type | |
601 | // and dedx scenario | |
602 | // | |
603 | ||
f84b18dd | 604 | Double_t dEdx = -1; |
605 | Int_t nPoints = -1; | |
606 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
607 | TSpline3* responseFunction = 0x0; | |
608 | ||
609 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
610 | return 999; //TODO: better handling! | |
611 | ||
87da0205 | 612 | return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity); |
644666df | 613 | } |
614 | ||
f84b18dd | 615 | |
644666df | 616 | //_________________________________________________________________________ |
617 | Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track, | |
618 | AliPID::EParticleType species, | |
f84b18dd | 619 | ETPCdEdxSource dedxSource, |
87da0205 | 620 | Bool_t correctEta, |
621 | Bool_t correctMultiplicity) const | |
644666df | 622 | { |
f84b18dd | 623 | //Calculates the number of sigmas of the PID signal from the expected value |
624 | //for a given particle species in the presence of multiple gain scenarios | |
644666df | 625 | //inside the TPC |
567624b5 | 626 | |
f84b18dd | 627 | Double_t dEdx = -1; |
628 | Int_t nPoints = -1; | |
629 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
630 | TSpline3* responseFunction = 0x0; | |
631 | ||
632 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
633 | return -999; //TODO: Better handling! | |
567624b5 | 634 | |
87da0205 | 635 | Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity); |
636 | Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity); | |
f84b18dd | 637 | // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters |
638 | if (sigma >= 998) | |
639 | return -999; | |
640 | else | |
641 | return (dEdx-bethe)/sigma; | |
644666df | 642 | } |
643 | ||
567624b5 | 644 | //_________________________________________________________________________ |
645 | Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track, | |
646 | AliPID::EParticleType species, | |
647 | ETPCdEdxSource dedxSource, | |
1d59271b | 648 | Bool_t correctEta, |
87da0205 | 649 | Bool_t correctMultiplicity, |
650 | Bool_t ratio/*=kFALSE*/)const | |
567624b5 | 651 | { |
652 | //Calculates the number of sigmas of the PID signal from the expected value | |
653 | //for a given particle species in the presence of multiple gain scenarios | |
654 | //inside the TPC | |
655 | ||
656 | Double_t dEdx = -1; | |
657 | Int_t nPoints = -1; | |
658 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
659 | TSpline3* responseFunction = 0x0; | |
660 | ||
661 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
662 | return -9999.; //TODO: Better handling! | |
663 | ||
87da0205 | 664 | const Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity); |
1d59271b | 665 | |
666 | Double_t delta=-9999.; | |
667 | if (!ratio) delta=dEdx-bethe; | |
668 | else if (bethe>1.e-20) delta=dEdx/bethe; | |
87da0205 | 669 | |
1d59271b | 670 | return delta; |
567624b5 | 671 | } |
672 | ||
644666df | 673 | //_________________________________________________________________________ |
674 | Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track, | |
675 | AliPID::EParticleType species, | |
f84b18dd | 676 | ETPCdEdxSource dedxSource, |
1b45e564 | 677 | Double_t& dEdx, |
678 | Int_t& nPoints, | |
679 | ETPCgainScenario& gainScenario, | |
f85a3764 | 680 | TSpline3** responseFunction) const |
644666df | 681 | { |
682 | // Calculates the right parameters for PID | |
683 | // dEdx parametrization for the proper gain scenario, dEdx | |
684 | // and NPoints used for dEdx | |
685 | // based on the track geometry (which chambers it crosses) for the specified particle type | |
686 | // and preferred source of dedx. | |
687 | // returns true on success | |
688 | ||
f84b18dd | 689 | |
690 | if (dedxSource == kdEdxDefault) { | |
691 | // Fast handling for default case. In addition: Keep it simple (don't call additional functions) to | |
692 | // avoid possible bugs | |
f85a3764 | 693 | |
694 | // GetTPCsignalTunedOnData will be non-positive, if it has not been set (i.e. in case of MC NOT tuned to data). | |
695 | // If this is the case, just take the normal signal | |
696 | dEdx = track->GetTPCsignalTunedOnData(); | |
697 | if (dEdx <= 0) { | |
698 | dEdx = track->GetTPCsignal(); | |
699 | } | |
700 | ||
f84b18dd | 701 | nPoints = track->GetTPCsignalN(); |
702 | gainScenario = kDefault; | |
703 | ||
704 | TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario)); | |
705 | *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast? | |
706 | ||
707 | return kTRUE; | |
708 | } | |
709 | ||
f85a3764 | 710 | //TODO Proper handle of tuneMConData for other dEdx sources |
f84b18dd | 711 | |
644666df | 712 | Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used) |
713 | Char_t ncl[3]; //same | |
714 | Char_t nrows[3]; //same | |
f84b18dd | 715 | const AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo(); |
644666df | 716 | |
717 | if (!dEdxInfo && dedxSource!=kdEdxDefault) //in one case its ok if we dont have the info | |
718 | { | |
719 | AliError("AliTPCdEdxInfo not available"); | |
644666df | 720 | return kFALSE; |
721 | } | |
722 | ||
723 | if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows); | |
724 | ||
725 | //check if we cross a bad OROC in which case we reject | |
f85a3764 | 726 | EChamberStatus trackOROCStatus = TrackStatus(track,2); |
727 | if (trackOROCStatus==kChamberOff || trackOROCStatus==kChamberLowGain) | |
644666df | 728 | { |
644666df | 729 | return kFALSE; |
730 | } | |
731 | ||
732 | switch (dedxSource) | |
733 | { | |
644666df | 734 | case kdEdxOROC: |
735 | { | |
f85a3764 | 736 | if (trackOROCStatus==kChamberInvalid) return kFALSE; //never reached OROC |
f84b18dd | 737 | dEdx = signal[3]; |
738 | nPoints = ncl[2]+ncl[1]; | |
739 | gainScenario = kOROChigh; | |
644666df | 740 | break; |
741 | } | |
742 | case kdEdxHybrid: | |
743 | { | |
744 | //if we cross a bad IROC we use OROC dedx, if we dont we use combined | |
745 | EChamberStatus status = TrackStatus(track,1); | |
746 | if (status!=kChamberHighGain) | |
747 | { | |
f84b18dd | 748 | dEdx = signal[3]; |
749 | nPoints = ncl[2]+ncl[1]; | |
750 | gainScenario = kOROChigh; | |
644666df | 751 | } |
752 | else | |
753 | { | |
f84b18dd | 754 | dEdx = track->GetTPCsignal(); |
755 | nPoints = track->GetTPCsignalN(); | |
756 | gainScenario = kALLhigh; | |
644666df | 757 | } |
758 | break; | |
759 | } | |
760 | default: | |
761 | { | |
f84b18dd | 762 | dEdx = 0.; |
763 | nPoints = 0; | |
764 | gainScenario = kGainScenarioInvalid; | |
644666df | 765 | return kFALSE; |
766 | } | |
767 | } | |
f84b18dd | 768 | TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario)); |
769 | *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast? | |
770 | ||
644666df | 771 | return kTRUE; |
772 | } | |
773 | ||
f84b18dd | 774 | |
775 | //_________________________________________________________________________ | |
ef7661fd | 776 | Double_t AliTPCPIDResponse::GetEtaCorrectionFast(const AliVTrack *track, Double_t dEdxSplines) const |
f84b18dd | 777 | { |
ef7661fd | 778 | // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly. |
f84b18dd | 779 | // |
780 | // Get eta correction for the given parameters. | |
781 | // | |
782 | ||
5a9dc560 | 783 | if (!fhEtaCorr) { |
784 | // Calling this function means to request eta correction in some way. Print error message, if no map is available! | |
785 | AliError("Eta correction requested, but map not initialised (usually via AliPIDResponse). Returning eta correction factor 1!"); | |
f84b18dd | 786 | return 1.; |
5a9dc560 | 787 | } |
f84b18dd | 788 | |
789 | Double_t tpcSignal = dEdxSplines; | |
790 | ||
791 | if (tpcSignal < 1.) | |
792 | return 1.; | |
793 | ||
87da0205 | 794 | Double_t tanTheta = GetTrackTanTheta(track); |
f84b18dd | 795 | Int_t binX = fhEtaCorr->GetXaxis()->FindBin(tanTheta); |
796 | Int_t binY = fhEtaCorr->GetYaxis()->FindBin(1. / tpcSignal); | |
797 | ||
798 | if (binX == 0) | |
799 | binX = 1; | |
800 | if (binX > fhEtaCorr->GetXaxis()->GetNbins()) | |
801 | binX = fhEtaCorr->GetXaxis()->GetNbins(); | |
802 | ||
803 | if (binY == 0) | |
804 | binY = 1; | |
805 | if (binY > fhEtaCorr->GetYaxis()->GetNbins()) | |
806 | binY = fhEtaCorr->GetYaxis()->GetNbins(); | |
807 | ||
808 | return fhEtaCorr->GetBinContent(binX, binY); | |
809 | } | |
810 | ||
811 | ||
812 | //_________________________________________________________________________ | |
813 | Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
814 | { | |
815 | // | |
816 | // Get eta correction for the given track. | |
817 | // | |
818 | ||
819 | if (!fhEtaCorr) | |
820 | return 1.; | |
821 | ||
822 | Double_t dEdx = -1; | |
823 | Int_t nPoints = -1; | |
824 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
825 | TSpline3* responseFunction = 0x0; | |
826 | ||
827 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
828 | return 1.; | |
829 | ||
87da0205 | 830 | // For the eta correction, do NOT take the multiplicity corrected value of dEdx |
831 | Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE); | |
f84b18dd | 832 | |
833 | //TODO Alternatively take current track dEdx | |
ef7661fd | 834 | //return GetEtaCorrectionFast(track, dEdx); |
f84b18dd | 835 | |
ef7661fd | 836 | return GetEtaCorrectionFast(track, dEdxSplines); |
f84b18dd | 837 | } |
838 | ||
839 | ||
840 | //_________________________________________________________________________ | |
841 | Double_t AliTPCPIDResponse::GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
842 | { | |
843 | // | |
844 | // Get eta corrected dEdx for the given track. For the correction, the expected dEdx of | |
845 | // the specified species will be used. If the species is set to AliPID::kUnknown, the | |
846 | // dEdx of the track is used instead. | |
847 | // WARNING: In the latter case, the eta correction might not be as good as if the | |
848 | // expected dEdx is used, which is the way the correction factor is designed | |
849 | // for. | |
850 | // In any case, one has to decide carefully to which expected signal one wants to | |
851 | // compare the corrected value - to the corrected or uncorrected. | |
852 | // Anyhow, a safer way of looking e.g. at the n-sigma is to call | |
853 | // the corresponding function GetNumberOfSigmas! | |
854 | // | |
855 | ||
856 | Double_t dEdx = -1; | |
857 | Int_t nPoints = -1; | |
858 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
859 | TSpline3* responseFunction = 0x0; | |
860 | ||
861 | // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case | |
862 | // it is not used anyway, so this causes no trouble. | |
863 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
864 | return -1.; | |
865 | ||
866 | Double_t etaCorr = 0.; | |
867 | ||
868 | if (species < AliPID::kUnknown) { | |
87da0205 | 869 | // For the eta correction, do NOT take the multiplicity corrected value of dEdx |
870 | Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE); | |
ef7661fd | 871 | etaCorr = GetEtaCorrectionFast(track, dEdxSplines); |
f84b18dd | 872 | } |
873 | else { | |
ef7661fd | 874 | etaCorr = GetEtaCorrectionFast(track, dEdx); |
f84b18dd | 875 | } |
876 | ||
877 | if (etaCorr <= 0) | |
878 | return -1.; | |
879 | ||
880 | return dEdx / etaCorr; | |
881 | } | |
882 | ||
883 | ||
f84b18dd | 884 | //_________________________________________________________________________ |
ef7661fd | 885 | Double_t AliTPCPIDResponse::GetSigmaPar1Fast(const AliVTrack *track, AliPID::EParticleType species, Double_t dEdx, |
886 | const TSpline3* responseFunction) const | |
f84b18dd | 887 | { |
ef7661fd | 888 | // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly. |
f84b18dd | 889 | // |
890 | // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track. | |
891 | // | |
892 | ||
5a9dc560 | 893 | if (!fhEtaSigmaPar1) { |
894 | // Calling this function means to request new sigma parametrisation in some way. Print error message, if no map is available! | |
895 | AliError("New sigma parametrisation requested, but sigma map not initialised (usually via AliPIDResponse). Returning error value for sigma parameter1 = 999!"); | |
f84b18dd | 896 | return 999; |
5a9dc560 | 897 | } |
f84b18dd | 898 | |
899 | // The sigma maps are created with data sets that are already eta corrected and for which the | |
900 | // splines have been re-created. Therefore, the value for the lookup needs to be the value of | |
901 | // the splines without any additional eta correction. | |
902 | // NOTE: This is due to the method the maps are created. The track dEdx (not the expected one!) | |
903 | // is corrected to uniquely related a momemtum bin with an expected dEdx, where the expected dEdx | |
904 | // equals the track dEdx for all eta (thanks to the correction and the re-fit of the splines). | |
905 | // Consequently, looking up the uncorrected expected dEdx at a given tanTheta yields the correct | |
906 | // sigma parameter! | |
907 | // Also: It has to be the spline dEdx, since one wants to get the sigma for the assumption(!) | |
908 | // of such a particle, which by assumption then has this dEdx value | |
87da0205 | 909 | |
910 | // For the eta correction, do NOT take the multiplicity corrected value of dEdx | |
911 | Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE); | |
f84b18dd | 912 | |
913 | if (dEdxExpected < 1.) | |
914 | return 999; | |
915 | ||
ef7661fd | 916 | Double_t tanTheta = GetTrackTanTheta(track); |
f84b18dd | 917 | Int_t binX = fhEtaSigmaPar1->GetXaxis()->FindBin(tanTheta); |
918 | Int_t binY = fhEtaSigmaPar1->GetYaxis()->FindBin(1. / dEdxExpected); | |
919 | ||
920 | if (binX == 0) | |
921 | binX = 1; | |
922 | if (binX > fhEtaSigmaPar1->GetXaxis()->GetNbins()) | |
923 | binX = fhEtaSigmaPar1->GetXaxis()->GetNbins(); | |
924 | ||
925 | if (binY == 0) | |
926 | binY = 1; | |
927 | if (binY > fhEtaSigmaPar1->GetYaxis()->GetNbins()) | |
928 | binY = fhEtaSigmaPar1->GetYaxis()->GetNbins(); | |
929 | ||
930 | return fhEtaSigmaPar1->GetBinContent(binX, binY); | |
931 | } | |
932 | ||
933 | ||
934 | //_________________________________________________________________________ | |
935 | Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
936 | { | |
937 | // | |
87da0205 | 938 | // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track. |
f84b18dd | 939 | // |
940 | ||
941 | if (!fhEtaSigmaPar1) | |
942 | return 999; | |
943 | ||
944 | Double_t dEdx = -1; | |
945 | Int_t nPoints = -1; | |
946 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
947 | TSpline3* responseFunction = 0x0; | |
948 | ||
949 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
950 | return 999; | |
951 | ||
ef7661fd | 952 | return GetSigmaPar1Fast(track, species, dEdx, responseFunction); |
f84b18dd | 953 | } |
954 | ||
955 | ||
956 | //_________________________________________________________________________ | |
957 | Bool_t AliTPCPIDResponse::SetEtaCorrMap(TH2D* hMap) | |
958 | { | |
959 | // | |
960 | // Load map for TPC eta correction (a copy is stored and will be deleted automatically). | |
961 | // If hMap is 0x0,the eta correction will be disabled and kFALSE is returned. | |
962 | // If the map can be set, kTRUE is returned. | |
963 | // | |
964 | ||
965 | delete fhEtaCorr; | |
966 | ||
967 | if (!hMap) { | |
968 | fhEtaCorr = 0x0; | |
969 | ||
970 | return kFALSE; | |
971 | } | |
972 | ||
973 | fhEtaCorr = (TH2D*)(hMap->Clone()); | |
974 | fhEtaCorr->SetDirectory(0); | |
975 | ||
976 | return kTRUE; | |
977 | } | |
978 | ||
87da0205 | 979 | |
f84b18dd | 980 | //_________________________________________________________________________ |
981 | Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0) | |
982 | { | |
983 | // | |
984 | // Load map for TPC sigma map (a copy is stored and will be deleted automatically): | |
985 | // Parameter 1 is stored as a 2D map (1/dEdx vs. tanTheta_local) and parameter 0 is | |
986 | // a constant. If hSigmaPar1Map is 0x0, the old sigma parametrisation will be used | |
987 | // (and sigmaPar0 is ignored!) and kFALSE is returned. | |
988 | // If the map can be set, sigmaPar0 is also set and kTRUE will be returned. | |
989 | // | |
990 | ||
991 | delete fhEtaSigmaPar1; | |
992 | ||
993 | if (!hSigmaPar1Map) { | |
994 | fhEtaSigmaPar1 = 0x0; | |
995 | fSigmaPar0 = 0.0; | |
996 | ||
997 | return kFALSE; | |
998 | } | |
999 | ||
1000 | fhEtaSigmaPar1 = (TH2D*)(hSigmaPar1Map->Clone()); | |
1001 | fhEtaSigmaPar1->SetDirectory(0); | |
1002 | fSigmaPar0 = sigmaPar0; | |
1003 | ||
1004 | return kTRUE; | |
1005 | } | |
1006 | ||
1007 | ||
87da0205 | 1008 | //_________________________________________________________________________ |
1009 | Double_t AliTPCPIDResponse::GetTrackTanTheta(const AliVTrack *track) const | |
1010 | { | |
1011 | // Extract the tanTheta from the information available in the AliVTrack | |
1012 | ||
1013 | // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()). | |
1014 | // However, this value is not available for AODs and, thus, not for AliVTrack. | |
1015 | // Fortunately, the following formula allows to approximate the local tanTheta with the | |
1016 | // global theta angle -> This is for by far most of the tracks the same, but gives at | |
1017 | // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok. | |
1018 | ||
ef7661fd | 1019 | /* |
1020 | const AliExternalTrackParam* innerParam = track->GetInnerParam(); | |
1021 | Double_t tanTheta = 0; | |
1022 | if (innerParam) | |
1023 | tanTheta = innerParam->GetTgl(); | |
1024 | else | |
1025 | tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0); | |
1026 | ||
1027 | // Constant in formula for B in kGauss (factor 0.1 to convert B from Tesla to kGauss), | |
1028 | // pT in GeV/c (factor c*1e-9 to arrive at GeV/c) and curvature in 1/cm (factor 0.01 to get from m to cm) | |
1029 | const Double_t constant = TMath::C()* 1e-9 * 0.1 * 0.01; | |
1030 | const Double_t curvature = fMagField * constant / track->Pt(); // in 1./cm | |
1031 | ||
1032 | Double_t averageddzdr = 0.; | |
1033 | Int_t nParts = 0; | |
1034 | ||
1035 | for (Double_t r = 85; r < 245; r++) { | |
1036 | Double_t sinPhiLocal = TMath::Abs(r*curvature*0.5); | |
1037 | ||
1038 | // Cut on |sin(phi)| as during reco | |
1039 | if (TMath::Abs(sinPhiLocal) <= 0.95) { | |
1040 | const Double_t phiLocal = TMath::ASin(sinPhiLocal); | |
1041 | const Double_t tanPhiLocal = TMath::Tan(phiLocal); | |
1042 | ||
1043 | averageddzdr += tanTheta * TMath::Sqrt(1. + tanPhiLocal * tanPhiLocal); | |
1044 | nParts++; | |
1045 | } | |
1046 | } | |
1047 | ||
1048 | if (nParts > 0) | |
1049 | averageddzdr /= nParts; | |
1050 | else { | |
1051 | AliError("Problems during determination of dz/dr. Returning pure tanTheta as best estimate!"); | |
1052 | return tanTheta; | |
1053 | } | |
1054 | ||
1055 | //printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\ntanThetaGlobalFromTheta/tanTheta/Averageddzdr: %f / %f / %f\n\n", | |
1056 | // track->Pt(), constant, fMagField, 1./curvature, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, averageddzdr); | |
1057 | ||
1058 | return averageddzdr; | |
1059 | */ | |
1060 | ||
87da0205 | 1061 | |
1062 | // Alternatively (in average, the effect is found to be negligable!): | |
1063 | // Take local tanTheta from TPC inner wall, if available (currently only for ESDs available) | |
ef7661fd | 1064 | //const AliExternalTrackParam* innerParam = track->GetInnerParam(); |
1065 | //if (innerParam) { | |
1066 | // return innerParam->GetTgl(); | |
1067 | //} | |
87da0205 | 1068 | |
1069 | return TMath::Tan(-track->Theta() + TMath::Pi() / 2.0); | |
1070 | } | |
1071 | ||
1072 | ||
1073 | //_________________________________________________________________________ | |
ef7661fd | 1074 | Double_t AliTPCPIDResponse::GetMultiplicityCorrectionFast(const AliVTrack *track, const Double_t dEdxExpected, const Int_t multiplicity) const |
87da0205 | 1075 | { |
ef7661fd | 1076 | // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly. |
1077 | // | |
87da0205 | 1078 | // Calculate the multiplicity correction factor for this track for the given multiplicity. |
1079 | // The parameter dEdxExpected should take into account the eta correction already! | |
1080 | ||
1081 | // Multiplicity depends on pure dEdx. Therefore, correction factor depends indirectly on eta | |
1082 | // => Use eta corrected value for dEdexExpected. | |
1083 | ||
1084 | if (dEdxExpected <= 0 || multiplicity <= 0) | |
1085 | return 1.0; | |
1086 | ||
1087 | const Double_t dEdxExpectedInv = 1. / dEdxExpected; | |
1088 | Double_t relSlope = fCorrFuncMultiplicity->Eval(dEdxExpectedInv); | |
1089 | ||
1090 | const Double_t tanTheta = GetTrackTanTheta(track); | |
1091 | relSlope += fCorrFuncMultiplicityTanTheta->Eval(tanTheta); | |
1092 | ||
1093 | return (1. + relSlope * multiplicity); | |
1094 | } | |
1095 | ||
1096 | ||
1097 | //_________________________________________________________________________ | |
1098 | Double_t AliTPCPIDResponse::GetMultiplicityCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
1099 | { | |
1100 | // | |
1101 | // Get multiplicity correction for the given track (w.r.t. the mulitplicity of the current event) | |
1102 | // | |
1103 | ||
1104 | //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse | |
1105 | ||
1106 | // No correction in case of multiplicity <= 0 | |
1107 | if (fCurrentEventMultiplicity <= 0) | |
1108 | return 1.; | |
1109 | ||
1110 | Double_t dEdx = -1; | |
1111 | Int_t nPoints = -1; | |
1112 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
1113 | TSpline3* responseFunction = 0x0; | |
1114 | ||
1115 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
1116 | return 1.; | |
1117 | ||
1118 | //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid? | |
1119 | // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course) | |
1120 | Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE); | |
1121 | ||
ef7661fd | 1122 | return GetMultiplicityCorrectionFast(track, dEdxExpected, fCurrentEventMultiplicity); |
87da0205 | 1123 | } |
1124 | ||
1125 | ||
1126 | //_________________________________________________________________________ | |
1127 | Double_t AliTPCPIDResponse::GetMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
1128 | { | |
1129 | // | |
1130 | // Get multiplicity corrected dEdx for the given track. For the correction, the expected dEdx of | |
1131 | // the specified species will be used. If the species is set to AliPID::kUnknown, the | |
1132 | // dEdx of the track is used instead. | |
1133 | // WARNING: In the latter case, the correction might not be as good as if the | |
1134 | // expected dEdx is used, which is the way the correction factor is designed | |
1135 | // for. | |
1136 | // In any case, one has to decide carefully to which expected signal one wants to | |
1137 | // compare the corrected value - to the corrected or uncorrected. | |
1138 | // Anyhow, a safer way of looking e.g. at the n-sigma is to call | |
1139 | // the corresponding function GetNumberOfSigmas! | |
1140 | // | |
1141 | ||
1142 | //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse | |
1143 | ||
1144 | Double_t dEdx = -1; | |
1145 | Int_t nPoints = -1; | |
1146 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
1147 | TSpline3* responseFunction = 0x0; | |
1148 | ||
1149 | // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case | |
1150 | // it is not used anyway, so this causes no trouble. | |
1151 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
1152 | return -1.; | |
1153 | ||
1154 | ||
1155 | // No correction in case of multiplicity <= 0 | |
1156 | if (fCurrentEventMultiplicity <= 0) | |
1157 | return dEdx; | |
1158 | ||
1159 | Double_t multiplicityCorr = 0.; | |
1160 | ||
1161 | // TODO Normally, one should use the eta corrected values in BOTH of the following cases. Does it make sense to use the uncorrected dEdx values? | |
1162 | if (species < AliPID::kUnknown) { | |
1163 | // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course). | |
1164 | // However, one needs the eta corrected value! | |
1165 | Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE); | |
ef7661fd | 1166 | multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdxSplines, fCurrentEventMultiplicity); |
87da0205 | 1167 | } |
1168 | else { | |
1169 | // One needs the eta corrected value to determine the multiplicity correction factor! | |
ef7661fd | 1170 | Double_t etaCorr = GetEtaCorrectionFast(track, dEdx); |
1171 | multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdx * etaCorr, fCurrentEventMultiplicity); | |
87da0205 | 1172 | } |
1173 | ||
1174 | if (multiplicityCorr <= 0) | |
1175 | return -1.; | |
1176 | ||
1177 | return dEdx / multiplicityCorr; | |
1178 | } | |
1179 | ||
1180 | ||
1181 | //_________________________________________________________________________ | |
1182 | Double_t AliTPCPIDResponse::GetEtaAndMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, | |
1183 | ETPCdEdxSource dedxSource) const | |
1184 | { | |
1185 | // | |
1186 | // Get multiplicity and eta corrected dEdx for the given track. For the correction, | |
1187 | // the expected dEdx of the specified species will be used. If the species is set | |
1188 | // to AliPID::kUnknown, the dEdx of the track is used instead. | |
1189 | // WARNING: In the latter case, the correction might not be as good as if the | |
1190 | // expected dEdx is used, which is the way the correction factor is designed | |
1191 | // for. | |
1192 | // In any case, one has to decide carefully to which expected signal one wants to | |
1193 | // compare the corrected value - to the corrected or uncorrected. | |
1194 | // Anyhow, a safer way of looking e.g. at the n-sigma is to call | |
1195 | // the corresponding function GetNumberOfSigmas! | |
1196 | // | |
1197 | ||
1198 | Double_t dEdx = -1; | |
1199 | Int_t nPoints = -1; | |
1200 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
1201 | TSpline3* responseFunction = 0x0; | |
1202 | ||
1203 | // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case | |
1204 | // it is not used anyway, so this causes no trouble. | |
1205 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
1206 | return -1.; | |
1207 | ||
1208 | Double_t multiplicityCorr = 0.; | |
1209 | Double_t etaCorr = 0.; | |
1210 | ||
1211 | if (species < AliPID::kUnknown) { | |
1212 | // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course) | |
1213 | Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE); | |
ef7661fd | 1214 | etaCorr = GetEtaCorrectionFast(track, dEdxSplines); |
1215 | multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdxSplines * etaCorr, fCurrentEventMultiplicity); | |
87da0205 | 1216 | } |
1217 | else { | |
ef7661fd | 1218 | etaCorr = GetEtaCorrectionFast(track, dEdx); |
1219 | multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdx * etaCorr, fCurrentEventMultiplicity); | |
87da0205 | 1220 | } |
1221 | ||
1222 | if (multiplicityCorr <= 0 || etaCorr <= 0) | |
1223 | return -1.; | |
1224 | ||
1225 | return dEdx / multiplicityCorr / etaCorr; | |
1226 | } | |
1227 | ||
1228 | ||
1229 | //_________________________________________________________________________ | |
ef7661fd | 1230 | Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrectionFast(const Double_t dEdxExpected, const Int_t multiplicity) const |
87da0205 | 1231 | { |
ef7661fd | 1232 | // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly. |
1233 | // | |
87da0205 | 1234 | // Calculate the multiplicity sigma correction factor for the corresponding expected dEdx and for the given multiplicity. |
1235 | // The parameter dEdxExpected should take into account the eta correction already! | |
1236 | ||
1237 | // Multiplicity dependence of sigma depends on the real dEdx at zero multiplicity, | |
1238 | // i.e. the eta (only) corrected dEdxExpected value has to be used | |
1239 | // since all maps etc. have been created for ~zero multiplicity | |
1240 | ||
1241 | if (dEdxExpected <= 0 || multiplicity <= 0) | |
1242 | return 1.0; | |
1243 | ||
1244 | Double_t relSigmaSlope = fCorrFuncSigmaMultiplicity->Eval(1. / dEdxExpected); | |
1245 | ||
1246 | return (1. + relSigmaSlope * multiplicity); | |
1247 | } | |
1248 | ||
1249 | ||
1250 | //_________________________________________________________________________ | |
1251 | Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
1252 | { | |
1253 | // | |
1254 | // Get multiplicity sigma correction for the given track (w.r.t. the mulitplicity of the current event) | |
1255 | // | |
1256 | ||
1257 | //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse | |
1258 | ||
1259 | // No correction in case of multiplicity <= 0 | |
1260 | if (fCurrentEventMultiplicity <= 0) | |
1261 | return 1.; | |
1262 | ||
1263 | Double_t dEdx = -1; | |
1264 | Int_t nPoints = -1; | |
1265 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
1266 | TSpline3* responseFunction = 0x0; | |
1267 | ||
1268 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
1269 | return 1.; | |
1270 | ||
1271 | //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid? | |
1272 | // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course) | |
1273 | Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE); | |
1274 | ||
ef7661fd | 1275 | return GetMultiplicitySigmaCorrectionFast(dEdxExpected, fCurrentEventMultiplicity); |
87da0205 | 1276 | } |
1277 | ||
1278 | ||
1279 | //_________________________________________________________________________ | |
1280 | void AliTPCPIDResponse::ResetMultiplicityCorrectionFunctions() | |
1281 | { | |
1282 | // Default values: No correction, i.e. overall correction factor should be one | |
1283 | ||
1284 | // This function should always return 0.0, if multcorr disabled | |
1285 | fCorrFuncMultiplicity->SetParameter(0, 0.); | |
1286 | fCorrFuncMultiplicity->SetParameter(1, 0.); | |
1287 | fCorrFuncMultiplicity->SetParameter(2, 0.); | |
1288 | fCorrFuncMultiplicity->SetParameter(3, 0.); | |
1289 | fCorrFuncMultiplicity->SetParameter(4, 0.); | |
1290 | ||
1291 | // This function should always return 0., if multCorr disabled | |
1292 | fCorrFuncMultiplicityTanTheta->SetParameter(0, 0.); | |
1293 | fCorrFuncMultiplicityTanTheta->SetParameter(1, 0.); | |
1294 | fCorrFuncMultiplicityTanTheta->SetParameter(2, 0.); | |
1295 | ||
1296 | // This function should always return 0.0, if mutlcorr disabled | |
1297 | fCorrFuncSigmaMultiplicity->SetParameter(0, 0.); | |
1298 | fCorrFuncSigmaMultiplicity->SetParameter(1, 0.); | |
1299 | fCorrFuncSigmaMultiplicity->SetParameter(2, 0.); | |
1300 | fCorrFuncSigmaMultiplicity->SetParameter(3, 0.); | |
1301 | } | |
1302 | ||
1303 | ||
644666df | 1304 | //_________________________________________________________________________ |
f85a3764 | 1305 | Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner, |
1306 | Double_t* trackPositionOuter, | |
1307 | Float_t& inphi, | |
644666df | 1308 | Float_t& outphi, |
f85a3764 | 1309 | Int_t& in, |
1b45e564 | 1310 | Int_t& out ) const |
644666df | 1311 | { |
1312 | //calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses | |
1313 | //for OROC chamber numbers add 36 | |
1314 | //returned angles are between (0,2pi) | |
644666df | 1315 | |
1316 | inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]); | |
1317 | outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]); | |
1318 | ||
1319 | if (inphi<0) {inphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi | |
1320 | if (outphi<0) {outphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi | |
1321 | ||
1322 | in = sectorNumber(inphi); | |
1323 | out = sectorNumber(outphi); | |
1324 | ||
1325 | //for the C side (positive z) offset by 18 | |
1326 | if (trackPositionInner[2]>0.0) | |
1327 | { | |
1328 | in+=18; | |
1329 | out+=18; | |
1330 | } | |
1331 | return kTRUE; | |
1332 | } | |
1b45e564 | 1333 | |
1334 | ||
644666df | 1335 | //_____________________________________________________________________________ |
1336 | Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const | |
1337 | { | |
1338 | //calculate sector number | |
1339 | const Float_t width=TMath::TwoPi()/18.; | |
1340 | return TMath::Floor(phi/width); | |
1341 | } | |
1342 | ||
1b45e564 | 1343 | |
644666df | 1344 | //_____________________________________________________________________________ |
1345 | void AliTPCPIDResponse::Print(Option_t* /*option*/) const | |
1346 | { | |
1347 | //Print info | |
1348 | fResponseFunctions.Print(); | |
1349 | } | |
1350 | ||
1b45e564 | 1351 | |
644666df | 1352 | //_____________________________________________________________________________ |
1353 | AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const | |
1354 | { | |
1355 | //status of the track: if it crosses any bad chambers return kChamberOff; | |
1356 | //IROC:layer=1, OROC:layer=2 | |
1357 | if (layer<1 || layer>2) layer=1; | |
1358 | Int_t in=0; | |
1359 | Int_t out=0; | |
1360 | Float_t inphi=0.; | |
1361 | Float_t outphi=0.; | |
1362 | Float_t innerRadius = (layer==1)?83.0:133.7; | |
1363 | Float_t outerRadius = (layer==1)?133.5:247.7; | |
f85a3764 | 1364 | |
1365 | ///////////////////////////////////////////////////////////////////////////// | |
1366 | //find out where track enters and leaves the layer. | |
1367 | // | |
1b45e564 | 1368 | Double_t trackPositionInner[3]; |
1369 | Double_t trackPositionOuter[3]; | |
1370 | ||
f85a3764 | 1371 | //if there is no inner param this could mean we're using the AOD track, |
1372 | //we still can extrapolate from the vertex - so use those params. | |
1373 | const AliExternalTrackParam* ip = track->GetInnerParam(); | |
1374 | if (ip) track=ip; | |
1375 | ||
1376 | Bool_t trackAtInner = track->GetXYZAt(innerRadius, fMagField, trackPositionInner); | |
1377 | Bool_t trackAtOuter = track->GetXYZAt(outerRadius, fMagField, trackPositionOuter); | |
1378 | ||
1379 | if (!trackAtInner) | |
1380 | { | |
1381 | //if we dont even enter inner radius we do nothing and return invalid | |
1382 | inphi=0.0; | |
1383 | outphi=0.0; | |
1384 | in=0; | |
1385 | out=0; | |
1386 | return kChamberInvalid; | |
1387 | } | |
1388 | ||
1389 | if (!trackAtOuter) | |
1390 | { | |
1391 | //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position | |
1392 | Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter); | |
1393 | Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]); | |
1394 | if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius) | |
1395 | { | |
1396 | //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]); | |
1397 | } | |
1398 | else | |
1399 | { | |
1400 | inphi=0.0; | |
1401 | outphi=0.0; | |
1402 | in=0; | |
1403 | out=0; | |
1404 | return kChamberInvalid; | |
1405 | } | |
1406 | } | |
1407 | ||
1408 | ||
1b45e564 | 1409 | if (!sectorNumbersInOut(trackPositionInner, |
1410 | trackPositionOuter, | |
1411 | inphi, | |
1412 | outphi, | |
1413 | in, | |
f85a3764 | 1414 | out)) return kChamberInvalid; |
644666df | 1415 | |
f85a3764 | 1416 | ///////////////////////////////////////////////////////////////////////////// |
1b45e564 | 1417 | //now we have the location of the track we can check |
f85a3764 | 1418 | //if it is in a good/bad chamber |
1419 | // | |
644666df | 1420 | Bool_t sideA = kTRUE; |
1b45e564 | 1421 | |
644666df | 1422 | if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE; |
1423 | ||
1424 | in=in%18; | |
1425 | out=out%18; | |
1426 | ||
1427 | if (TMath::Abs(in-out)>9) | |
1428 | { | |
1429 | if (TMath::Max(in,out)==out) | |
1430 | { | |
1431 | Int_t tmp=out; | |
1432 | out = in; | |
1433 | in = tmp; | |
1434 | Float_t tmpphi=outphi; | |
1435 | outphi=inphi; | |
1436 | inphi=tmpphi; | |
1437 | } | |
1438 | in-=18; | |
1439 | inphi-=TMath::TwoPi(); | |
1440 | } | |
10d100d4 | 1441 | else |
644666df | 1442 | { |
1443 | if (TMath::Max(in,out)==in) | |
1444 | { | |
1445 | Int_t tmp=out; | |
1446 | out=in; | |
1447 | in=tmp; | |
1448 | Float_t tmpphi=outphi; | |
1449 | outphi=inphi; | |
1450 | inphi=tmpphi; | |
1451 | } | |
1452 | } | |
1453 | ||
1454 | Float_t trackLengthInBad = 0.; | |
1455 | Float_t trackLengthInLowGain = 0.; | |
1456 | Float_t trackLengthTotal = TMath::Abs(outphi-inphi); | |
f85a3764 | 1457 | Float_t lengthFractionInBadSectors = 0.; |
644666df | 1458 | |
1459 | const Float_t sectorWidth = TMath::TwoPi()/18.; | |
1b45e564 | 1460 | |
644666df | 1461 | for (Int_t i=in; i<=out; i++) |
1462 | { | |
1463 | int j=i; | |
1464 | if (i<0) j+=18; //correct for the negative values | |
1465 | if (!sideA) j+=18; //move to the correct side | |
1b45e564 | 1466 | |
644666df | 1467 | Float_t deltaPhi = 0.; |
1468 | Float_t phiEdge=sectorWidth*i; | |
1469 | if (inphi>phiEdge) {deltaPhi=phiEdge+sectorWidth-inphi;} | |
1470 | else if ((outphi>=phiEdge) && (outphi<(phiEdge+sectorWidth))) {deltaPhi=outphi-phiEdge;} | |
1471 | else {deltaPhi=sectorWidth;} | |
1b45e564 | 1472 | |
644666df | 1473 | Float_t v = fVoltageMap[(layer==1)?(j):(j+36)]; |
1b45e564 | 1474 | if (v<=fBadIROCthreshhold) |
1475 | { | |
1476 | trackLengthInBad+=deltaPhi; | |
f85a3764 | 1477 | lengthFractionInBadSectors=1.; |
1478 | } | |
1b45e564 | 1479 | if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold) |
f85a3764 | 1480 | { |
1481 | trackLengthInLowGain+=deltaPhi; | |
1482 | lengthFractionInBadSectors=1.; | |
1483 | } | |
644666df | 1484 | } |
1485 | ||
1486 | //for now low gain and bad (off) chambers are treated equally | |
f85a3764 | 1487 | if (trackLengthTotal>0) |
1488 | lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal; | |
644666df | 1489 | |
1490 | //printf("### side: %s, pt: %.2f, pz: %.2f, in: %i, out: %i, phiIN: %.2f, phiOUT: %.2f, rIN: %.2f, rOUT: %.2f\n",(sideA)?"A":"C",track->Pt(),track->Pz(),in,out,inphi,outphi,innerRadius,outerRadius); | |
1b45e564 | 1491 | |
1492 | if (lengthFractionInBadSectors>fMaxBadLengthFraction) | |
644666df | 1493 | { |
1494 | //printf("%%%%%%%% %s kChamberLowGain\n",(layer==1)?"IROC":"OROC"); | |
1495 | return kChamberLowGain; | |
1496 | } | |
1b45e564 | 1497 | |
644666df | 1498 | return kChamberHighGain; |
1499 | } | |
1500 | ||
1b45e564 | 1501 | |
644666df | 1502 | //_____________________________________________________________________________ |
1503 | Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const | |
1504 | { | |
1505 | //return the radius of the outermost padrow containing a cluster in TPC | |
1506 | //for the track | |
1507 | const TBits* clusterMap=track->GetTPCClusterMapPtr(); | |
1508 | if (!clusterMap) return 0.; | |
1509 | ||
1510 | //from AliTPCParam, radius of first IROC row | |
1b45e564 | 1511 | const Float_t rfirstIROC = 8.52249984741210938e+01; |
644666df | 1512 | const Float_t padrowHeightIROC = 0.75; |
1513 | const Float_t rfirstOROC0 = 1.35100006103515625e+02; | |
1514 | const Float_t padrowHeightOROC0 = 1.0; | |
1515 | const Float_t rfirstOROC1 = 1.99350006103515625e+02; | |
1516 | const Float_t padrowHeightOROC1 = 1.5; | |
1517 | ||
1518 | Int_t maxPadRow=160; | |
1519 | while ((--maxPadRow)>0 && !clusterMap->TestBitNumber(maxPadRow)){} | |
1520 | if (maxPadRow>126) return rfirstOROC1+(maxPadRow-126-1)*padrowHeightOROC1; | |
1521 | if (maxPadRow>62) return rfirstOROC0+(maxPadRow-62-1)*padrowHeightOROC0; | |
1522 | if (maxPadRow>0) return rfirstIROC+(maxPadRow-1)*padrowHeightIROC; | |
1523 | return 0.0; | |
10d100d4 | 1524 | } |
644666df | 1525 | |
1b45e564 | 1526 | |
644666df | 1527 | //_____________________________________________________________________________ |
1528 | Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const | |
1529 | { | |
1530 | //calculate the coordinates of the apex of the track | |
1531 | Double_t x[3]; | |
1532 | track->GetXYZ(x); | |
1533 | Double_t p[3]; | |
1534 | track->GetPxPyPz(p); | |
1535 | Double_t r = 1./track->OneOverPt()/0.0299792458/magField; //signed - will determine the direction of b | |
1b45e564 | 1536 | //printf("b: %.2f, x:%.2f, y:%.2f, pt: %.2f, px:%.2f, py%.2f, r: %.2f\n",magField, x[0],x[1],track->Pt(), p[0],p[1],r); |
644666df | 1537 | //find orthogonal vector (Gram-Schmidt) |
1538 | Double_t alpha = (p[0]*x[0] + p[1]*x[1])/(p[0]*p[0] + p[1]*p[1]); | |
1539 | Double_t b[2]; | |
1b45e564 | 1540 | b[0]=x[0]-alpha*p[0]; |
1541 | b[1]=x[1]-alpha*p[1]; | |
1542 | ||
644666df | 1543 | Double_t norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); |
1544 | if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE; | |
1b45e564 | 1545 | b[0]/=norm; |
1546 | b[1]/=norm; | |
1547 | b[0]*=r; | |
1548 | b[1]*=r; | |
1549 | b[0]+=x[0]; | |
644666df | 1550 | b[1]+=x[1]; |
1551 | //printf("center: x:%.2f, y:%.2f\n",b[0],b[1]); | |
1b45e564 | 1552 | |
644666df | 1553 | norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); |
1554 | if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE; | |
1555 | ||
1556 | position[0]=b[0]+b[0]*TMath::Abs(r)/norm; | |
1557 | position[1]=b[1]+b[1]*TMath::Abs(r)/norm; | |
1558 | position[2]=0.; | |
1559 | return kTRUE; | |
f85a3764 | 1560 | } |
912f7fdc | 1561 | |
1562 | Double_t AliTPCPIDResponse::EvaldEdxSpline(Double_t bg,Int_t entry){ | |
1563 | // | |
1564 | // Evaluate the dEdx response for given entry | |
1565 | // | |
1566 | TSpline * spline = (TSpline*)fSplineArray.At(entry); | |
1567 | if (spline) return spline->Eval(bg); | |
1568 | return 0; | |
1569 | } | |
1570 | ||
1571 | ||
1572 | Bool_t AliTPCPIDResponse::RegisterSpline(const char * name, Int_t index){ | |
1573 | // | |
1574 | // register spline to be used for drawing comparisons | |
1575 | // | |
1576 | TFile * fTPCBB = TFile::Open("$ALICE_ROOT/OADB/COMMON/PID/data/TPCPIDResponse.root"); | |
1577 | TObjArray *arrayTPCPID= (TObjArray*) fTPCBB->Get("TPCPIDResponse"); | |
1578 | if (fSplineArray.GetEntriesFast()<index) fSplineArray.Expand(index*2); | |
1579 | TSpline3 *spline=0; | |
1580 | if (arrayTPCPID){ | |
1581 | spline = (TSpline3*)arrayTPCPID->FindObject(name); | |
1582 | if (spline) fSplineArray.AddAt(spline->Clone(),index); | |
1583 | } | |
1584 | delete arrayTPCPID; | |
1585 | delete fTPCBB; | |
1586 | return (spline!=0); | |
1587 | } |