Line data Source code
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
23 : // Jens Wiechula, IKF, Jens.Wiechula@cern.ch
24 : // ...and some modifications by
25 : // Mikolaj Krzewicki, GSI, mikolaj.krzewicki@cern.ch
26 : // ...and some modifications plus eta correction functions by
27 : // Benjamin Hess, University of Tuebingen, bhess@cern.ch
28 : //-----------------------------------------------------------------
29 :
30 : #include <TGraph.h>
31 : #include <TError.h>
32 : #include <TObjArray.h>
33 : #include <TSpline.h>
34 : #include <TBits.h>
35 : #include <TMath.h>
36 : #include <TH2D.h>
37 : #include <TSystem.h>
38 : #include <TMD5.h>
39 :
40 : #include <AliLog.h>
41 : #include "AliExternalTrackParam.h"
42 : #include "AliVTrack.h"
43 : #include "AliTPCPIDResponse.h"
44 : #include "AliTPCdEdxInfo.h"
45 : #include "AliOADBContainer.h"
46 : #include "TFile.h"
47 : #include "TSpline.h"
48 :
49 176 : ClassImp(AliTPCPIDResponse);
50 :
51 :
52 : AliTPCPIDResponse *AliTPCPIDResponse::fgInstance =0;
53 :
54 : const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]=
55 : {
56 : "", //default - no name
57 : "1", //all high
58 : "2", //OROC only
59 : "unknown"
60 : };
61 :
62 : //_________________________________________________________________________
63 : AliTPCPIDResponse::AliTPCPIDResponse():
64 10 : TNamed(),
65 10 : fMIP(50.),
66 10 : fRes0(),
67 10 : fResN2(),
68 10 : fKp1(0.0283086),
69 10 : fKp2(2.63394e+01),
70 10 : fKp3(5.04114e-11),
71 10 : fKp4(2.12543),
72 10 : fKp5(4.88663),
73 10 : fUseDatabase(kFALSE),
74 10 : fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
75 10 : fOADBContainer(0x0),
76 10 : fVoltageMap(72),
77 10 : fLowGainIROCthreshold(-40),
78 10 : fBadIROCthreshhold(-70),
79 10 : fLowGainOROCthreshold(-40),
80 10 : fBadOROCthreshhold(-40),
81 10 : fMaxBadLengthFraction(0.5),
82 10 : fMagField(0.),
83 10 : fhEtaCorr(0x0),
84 10 : fhEtaSigmaPar1(0x0),
85 10 : fSigmaPar0(0.0),
86 10 : fCurrentEventMultiplicity(0),
87 10 : fCorrFuncMultiplicity(0x0),
88 10 : fCorrFuncMultiplicityTanTheta(0x0),
89 10 : fCorrFuncSigmaMultiplicity(0x0),
90 10 : fdEdxType(kdEdxTrack),
91 10 : fdEdxChargeType(0),
92 10 : fdEdxWeightType(0),
93 10 : fIROCweight(1.),
94 10 : fOROCmedWeight(1.),
95 10 : fOROClongWeight(1.),
96 10 : fSplineArray()
97 50 : {
98 : //
99 : // The default constructor
100 : //
101 :
102 10 : AliLog::SetClassDebugLevel("AliTPCPIDResponse", AliLog::kInfo);
103 :
104 80 : for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;}
105 :
106 30 : fCorrFuncMultiplicity = new TF1("fCorrFuncMultiplicity",
107 : "[0] + [1]*TMath::Max([4], TMath::Min(x, [3])) + [2] * TMath::Power(TMath::Max([4], TMath::Min(x, [3])), 2)",
108 : 0., 0.2);
109 30 : fCorrFuncMultiplicityTanTheta = new TF1("fCorrFuncMultiplicityTanTheta", "[0] * (x - [2]) + [1] * (x * x - [2] * [2])", -1.5, 1.5);
110 30 : fCorrFuncSigmaMultiplicity = new TF1("fCorrFuncSigmaMultiplicity",
111 : "TMath::Max(0.0, [0] + [1]*TMath::Min(x, [3]) + [2] * TMath::Power(TMath::Min(x, [3]), 2))", 0., 0.2);
112 :
113 :
114 10 : ResetMultiplicityCorrectionFunctions();
115 10 : fgInstance=this;
116 20 : }
117 : /*TODO remove?
118 : //_________________________________________________________________________
119 : AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
120 : TNamed(),
121 : fMIP(param[0]),
122 : fRes0(),
123 : fResN2(),
124 : fKp1(0.0283086),
125 : fKp2(2.63394e+01),
126 : fKp3(5.04114e-11),
127 : fKp4(2.12543),
128 : fKp5(4.88663),
129 : fUseDatabase(kFALSE),
130 : fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
131 : fVoltageMap(72),
132 : fLowGainIROCthreshold(-40),
133 : fBadIROCthreshhold(-70),
134 : fLowGainOROCthreshold(-40),
135 : fBadOROCthreshhold(-40),
136 : fMaxBadLengthFraction(0.5),
137 : fMagField(0.),
138 : fhEtaCorr(0x0),
139 : fhEtaSigmaPar1(0x0),
140 : fSigmaPar0(0.0)
141 : {
142 : //
143 : // The main constructor
144 : //
145 : for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];}
146 : }
147 : */
148 :
149 : //_________________________________________________________________________
150 : AliTPCPIDResponse::~AliTPCPIDResponse()
151 40 : {
152 : //
153 : // Destructor
154 : //
155 :
156 10 : delete fhEtaCorr;
157 10 : fhEtaCorr = 0x0;
158 :
159 10 : delete fhEtaSigmaPar1;
160 10 : fhEtaSigmaPar1 = 0x0;
161 :
162 20 : delete fCorrFuncMultiplicity;
163 10 : fCorrFuncMultiplicity = 0x0;
164 :
165 20 : delete fCorrFuncMultiplicityTanTheta;
166 10 : fCorrFuncMultiplicityTanTheta = 0x0;
167 :
168 20 : delete fCorrFuncSigmaMultiplicity;
169 10 : fCorrFuncSigmaMultiplicity = 0x0;
170 20 : if (fgInstance==this) fgInstance=0;
171 :
172 10 : delete fOADBContainer;
173 20 : }
174 :
175 :
176 : //_________________________________________________________________________
177 : AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that):
178 0 : TNamed(that),
179 0 : fMIP(that.fMIP),
180 0 : fRes0(),
181 0 : fResN2(),
182 0 : fKp1(that.fKp1),
183 0 : fKp2(that.fKp2),
184 0 : fKp3(that.fKp3),
185 0 : fKp4(that.fKp4),
186 0 : fKp5(that.fKp5),
187 0 : fUseDatabase(that.fUseDatabase),
188 0 : fResponseFunctions(that.fResponseFunctions),
189 0 : fOADBContainer(0x0),
190 0 : fVoltageMap(that.fVoltageMap),
191 0 : fLowGainIROCthreshold(that.fLowGainIROCthreshold),
192 0 : fBadIROCthreshhold(that.fBadIROCthreshhold),
193 0 : fLowGainOROCthreshold(that.fLowGainOROCthreshold),
194 0 : fBadOROCthreshhold(that.fBadOROCthreshhold),
195 0 : fMaxBadLengthFraction(that.fMaxBadLengthFraction),
196 0 : fMagField(that.fMagField),
197 0 : fhEtaCorr(0x0),
198 0 : fhEtaSigmaPar1(0x0),
199 0 : fSigmaPar0(that.fSigmaPar0),
200 0 : fCurrentEventMultiplicity(that.fCurrentEventMultiplicity),
201 0 : fCorrFuncMultiplicity(0x0),
202 0 : fCorrFuncMultiplicityTanTheta(0x0),
203 0 : fCorrFuncSigmaMultiplicity(0x0),
204 0 : fdEdxType(kdEdxTrack),
205 0 : fdEdxChargeType(that.fdEdxChargeType),
206 0 : fdEdxWeightType(that.fdEdxWeightType),
207 0 : fIROCweight(that.fIROCweight),
208 0 : fOROCmedWeight(that.fOROCmedWeight),
209 0 : fOROClongWeight(that.fOROClongWeight),
210 0 : fSplineArray()
211 0 : {
212 : //copy ctor
213 0 : for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
214 :
215 : // Copy eta maps
216 0 : if (that.fhEtaCorr) {
217 0 : fhEtaCorr = new TH2D(*(that.fhEtaCorr));
218 0 : fhEtaCorr->SetDirectory(0);
219 : }
220 :
221 0 : if (that.fhEtaSigmaPar1) {
222 0 : fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
223 0 : fhEtaSigmaPar1->SetDirectory(0);
224 : }
225 :
226 : // Copy multiplicity correction functions
227 0 : if (that.fCorrFuncMultiplicity) {
228 0 : fCorrFuncMultiplicity = new TF1(*(that.fCorrFuncMultiplicity));
229 0 : }
230 :
231 0 : if (that.fCorrFuncMultiplicityTanTheta) {
232 0 : fCorrFuncMultiplicityTanTheta = new TF1(*(that.fCorrFuncMultiplicityTanTheta));
233 0 : }
234 :
235 0 : if (that.fCorrFuncSigmaMultiplicity) {
236 0 : fCorrFuncSigmaMultiplicity = new TF1(*(that.fCorrFuncSigmaMultiplicity));
237 0 : }
238 0 : }
239 :
240 : //_________________________________________________________________________
241 : AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
242 : {
243 : //assignment
244 0 : if (&that==this) return *this;
245 0 : TNamed::operator=(that);
246 0 : fMIP=that.fMIP;
247 0 : fKp1=that.fKp1;
248 0 : fKp2=that.fKp2;
249 0 : fKp3=that.fKp3;
250 0 : fKp4=that.fKp4;
251 0 : fKp5=that.fKp5;
252 0 : fUseDatabase=that.fUseDatabase;
253 0 : fResponseFunctions=that.fResponseFunctions;
254 0 : fOADBContainer=0x0;
255 0 : fVoltageMap=that.fVoltageMap;
256 0 : fLowGainIROCthreshold=that.fLowGainIROCthreshold;
257 0 : fBadIROCthreshhold=that.fBadIROCthreshhold;
258 0 : fLowGainOROCthreshold=that.fLowGainOROCthreshold;
259 0 : fBadOROCthreshhold=that.fBadOROCthreshhold;
260 0 : fMaxBadLengthFraction=that.fMaxBadLengthFraction;
261 0 : fMagField=that.fMagField;
262 0 : fCurrentEventMultiplicity=that.fCurrentEventMultiplicity;
263 0 : for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
264 :
265 0 : delete fhEtaCorr;
266 0 : fhEtaCorr=0x0;
267 0 : if (that.fhEtaCorr) {
268 0 : fhEtaCorr = new TH2D(*(that.fhEtaCorr));
269 0 : fhEtaCorr->SetDirectory(0);
270 0 : }
271 :
272 0 : delete fhEtaSigmaPar1;
273 0 : fhEtaSigmaPar1=0x0;
274 0 : if (that.fhEtaSigmaPar1) {
275 0 : fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
276 0 : fhEtaSigmaPar1->SetDirectory(0);
277 0 : }
278 :
279 0 : fSigmaPar0 = that.fSigmaPar0;
280 :
281 0 : delete fCorrFuncMultiplicity;
282 0 : fCorrFuncMultiplicity = 0x0;
283 0 : if (that.fCorrFuncMultiplicity) {
284 0 : fCorrFuncMultiplicity = new TF1(*(that.fCorrFuncMultiplicity));
285 0 : }
286 :
287 0 : delete fCorrFuncMultiplicityTanTheta;
288 0 : fCorrFuncMultiplicityTanTheta = 0x0;
289 0 : if (that.fCorrFuncMultiplicityTanTheta) {
290 0 : fCorrFuncMultiplicityTanTheta = new TF1(*(that.fCorrFuncMultiplicityTanTheta));
291 0 : }
292 :
293 0 : delete fCorrFuncSigmaMultiplicity;
294 0 : fCorrFuncSigmaMultiplicity = 0x0;
295 0 : if (that.fCorrFuncSigmaMultiplicity) {
296 0 : fCorrFuncSigmaMultiplicity = new TF1(*(that.fCorrFuncSigmaMultiplicity));
297 0 : }
298 :
299 0 : fdEdxType =that.fdEdxType;
300 0 : fdEdxChargeType=that.fdEdxChargeType;
301 0 : fdEdxWeightType=that.fdEdxWeightType;
302 0 : fIROCweight =that.fIROCweight;
303 0 : fOROCmedWeight =that.fOROCmedWeight;
304 0 : fOROClongWeight=that.fOROClongWeight;
305 :
306 0 : return *this;
307 0 : }
308 :
309 : //_________________________________________________________________________
310 : Double_t AliTPCPIDResponse::Bethe(Double_t betaGamma) const {
311 : //
312 : // This is the Bethe-Bloch function normalised to 1 at the minimum
313 : // WARNING
314 : // Simulated and reconstructed Bethe-Bloch differs
315 : // Simulated curve is the dNprim/dx
316 : // Reconstructed is proportianal dNtot/dx
317 : // Temporary fix for production - Simple linear correction function
318 : // Future 2 Bethe Bloch formulas needed
319 : // 1. for simulation
320 : // 2. for reconstructed PID
321 : //
322 :
323 : // const Float_t kmeanCorrection =0.1;
324 : Double_t bb=
325 4824 : AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5);
326 2412 : return bb*fMIP;
327 : }
328 :
329 : //_________________________________________________________________________
330 : void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1,
331 : Double_t kp2,
332 : Double_t kp3,
333 : Double_t kp4,
334 : Double_t kp5) {
335 : //
336 : // Set the parameters of the ALEPH Bethe-Bloch formula
337 : //
338 16 : fKp1=kp1;
339 8 : fKp2=kp2;
340 8 : fKp3=kp3;
341 8 : fKp4=kp4;
342 8 : fKp5=kp5;
343 8 : }
344 :
345 : //_________________________________________________________________________
346 : void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) {
347 : //
348 : // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint)
349 : //
350 0 : for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;}
351 0 : }
352 :
353 : //_________________________________________________________________________
354 : Double_t AliTPCPIDResponse::GetExpectedSignal(Float_t mom,
355 : AliPID::EParticleType n) const {
356 : //
357 : // Deprecated function (for backward compatibility). Please use
358 : // GetExpectedSignal(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource,
359 : // Bool_t correctEta, Bool_t correctMultiplicity);
360 : // instead!
361 : //
362 : //
363 : // Calculates the expected PID signal as the function of
364 : // the information stored in the track, for the specified particle type
365 : //
366 : // At the moment, these signals are just the results of calling the
367 : // Bethe-Bloch formula.
368 : // This can be improved. By taking into account the number of
369 : // assigned clusters and/or the track dip angle, for example.
370 : //
371 :
372 : //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
373 : // !!! Splines for light nuclei need to be normalised to this factor !!!
374 0 : const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3);
375 :
376 0 : Double_t mass=AliPID::ParticleMassZ(n);
377 0 : if (!fUseDatabase) return Bethe(mom/mass) * chargeFactor;
378 : //
379 0 : const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
380 :
381 0 : if (!responseFunction) return Bethe(mom/mass) * chargeFactor;
382 :
383 0 : return fMIP*responseFunction->Eval(mom/mass)*chargeFactor;
384 :
385 0 : }
386 :
387 : //_________________________________________________________________________
388 : Double_t AliTPCPIDResponse::GetExpectedSigma(Float_t mom,
389 : Int_t nPoints,
390 : AliPID::EParticleType n) const {
391 : //
392 : // Deprecated function (for backward compatibility). Please use
393 : // GetExpectedSigma(onst AliVTrack* track, AliPID::EParticleType species,
394 : // ETPCdEdxSource dedxSource, Bool_t correctEta) instead!
395 : //
396 : //
397 : // Calculates the expected sigma of the PID signal as the function of
398 : // the information stored in the track, for the specified particle type
399 : //
400 :
401 0 : if (nPoints != 0)
402 0 : return GetExpectedSignal(mom,n)*fRes0[0]*sqrt(1. + fResN2[0]/nPoints);
403 : else
404 0 : return GetExpectedSignal(mom,n)*fRes0[0];
405 0 : }
406 :
407 : ////////////////////////////////////////////////////NEW//////////////////////////////
408 :
409 : //_________________________________________________________________________
410 : void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario) {
411 : //
412 : // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint)
413 : //
414 0 : if ((Int_t)gainScenario>=(Int_t)fgkNumberOfGainScenarios) return; //TODO: better handling!
415 0 : fRes0[gainScenario]=res0;
416 0 : fResN2[gainScenario]=resN2;
417 0 : }
418 :
419 :
420 : //_________________________________________________________________________
421 : Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
422 : AliPID::EParticleType species,
423 : Double_t /*dEdx*/,
424 : const TSpline3* responseFunction,
425 : Bool_t correctEta,
426 : Bool_t correctMultiplicity) const
427 : {
428 : // Calculates the expected PID signal as the function of
429 : // the information stored in the track and the given parameters,
430 : // for the specified particle type
431 : //
432 : // At the moment, these signals are just the results of calling the
433 : // Bethe-Bloch formula plus, if desired, taking into account the eta dependence
434 : // and the multiplicity dependence (for PbPb).
435 : // This can be improved. By taking into account the number of
436 : // assigned clusters and/or the track dip angle, for example.
437 : //
438 :
439 1206 : Double_t mom=track->GetTPCmomentum();
440 1206 : Double_t mass=AliPID::ParticleMassZ(species);
441 :
442 : //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
443 : // !!! Splines for light nuclei need to be normalised to this factor !!!
444 1206 : const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3);
445 :
446 1206 : if (!responseFunction)
447 1206 : return Bethe(mom/mass) * chargeFactor;
448 :
449 0 : Double_t dEdxSplines = fMIP*responseFunction->Eval(mom/mass) * chargeFactor;
450 :
451 0 : if (!correctEta && !correctMultiplicity)
452 0 : return dEdxSplines;
453 :
454 : Double_t corrFactorEta = 1.0;
455 : Double_t corrFactorMultiplicity = 1.0;
456 :
457 0 : if (correctEta) {
458 0 : corrFactorEta = GetEtaCorrectionFast(track, dEdxSplines);
459 : //TODO Alternatively take current track dEdx
460 : //corrFactorEta = GetEtaCorrectionFast(track, dEdx);
461 0 : }
462 :
463 0 : if (correctMultiplicity)
464 0 : corrFactorMultiplicity = GetMultiplicityCorrectionFast(track, dEdxSplines * corrFactorEta, fCurrentEventMultiplicity);
465 :
466 0 : return dEdxSplines * corrFactorEta * corrFactorMultiplicity;
467 1206 : }
468 :
469 :
470 : //_________________________________________________________________________
471 : Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
472 : AliPID::EParticleType species,
473 : ETPCdEdxSource dedxSource,
474 : Bool_t correctEta,
475 : Bool_t correctMultiplicity) const
476 : {
477 : // Calculates the expected PID signal as the function of
478 : // the information stored in the track, for the specified particle type
479 : //
480 : // At the moment, these signals are just the results of calling the
481 : // Bethe-Bloch formula plus, if desired, taking into account the eta dependence
482 : // and the multiplicity dependence (for PbPb).
483 : // This can be improved. By taking into account the number of
484 : // assigned clusters and/or the track dip angle, for example.
485 : //
486 :
487 1206 : if (!fUseDatabase) {
488 : //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
489 : // !!! Splines for light nuclei need to be normalised to this factor !!!
490 1206 : const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3);
491 :
492 1206 : return Bethe(track->GetTPCmomentum() / AliPID::ParticleMassZ(species)) * chargeFactor;
493 : }
494 :
495 0 : Double_t dEdx = -1;
496 0 : Int_t nPoints = -1;
497 0 : ETPCgainScenario gainScenario = kGainScenarioInvalid;
498 0 : TSpline3* responseFunction = 0x0;
499 :
500 0 : if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) {
501 : // Something is wrong with the track -> Return obviously invalid value
502 0 : return -999;
503 : }
504 :
505 : // Charge factor already taken into account inside the following function call
506 0 : return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity);
507 1206 : }
508 :
509 : //_________________________________________________________________________
510 : TSpline3* AliTPCPIDResponse::GetResponseFunction( AliPID::EParticleType type,
511 : AliTPCPIDResponse::ETPCgainScenario gainScenario ) const
512 : {
513 : //get response function
514 0 : return dynamic_cast<TSpline3*>(fResponseFunctions.At(ResponseFunctionIndex(type,gainScenario)));
515 : }
516 :
517 : //_________________________________________________________________________
518 : TSpline3* AliTPCPIDResponse::GetResponseFunction( const AliVTrack* track,
519 : AliPID::EParticleType species,
520 : ETPCdEdxSource dedxSource) const
521 : {
522 : //the splines are stored in an array, different scenarios
523 :
524 0 : Double_t dEdx = -1;
525 0 : Int_t nPoints = -1;
526 0 : ETPCgainScenario gainScenario = kGainScenarioInvalid;
527 0 : TSpline3* responseFunction = 0x0;
528 :
529 0 : if (ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
530 0 : return responseFunction;
531 :
532 0 : return NULL;
533 0 : }
534 :
535 : //_________________________________________________________________________
536 : void AliTPCPIDResponse::ResetSplines()
537 : {
538 : //reset the array with splines
539 0 : for (Int_t i=0;i<fResponseFunctions.GetEntriesFast();i++)
540 : {
541 0 : fResponseFunctions.AddAt(NULL,i);
542 : }
543 0 : }
544 : //_________________________________________________________________________
545 : Int_t AliTPCPIDResponse::ResponseFunctionIndex( AliPID::EParticleType species,
546 : ETPCgainScenario gainScenario ) const
547 : {
548 : //get the index in fResponseFunctions given type and scenario
549 2412 : return Int_t(species)+Int_t(gainScenario)*fgkNumberOfParticleSpecies;
550 : }
551 :
552 : //_________________________________________________________________________
553 : void AliTPCPIDResponse::SetResponseFunction( TObject* o,
554 : AliPID::EParticleType species,
555 : ETPCgainScenario gainScenario )
556 : {
557 0 : fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario));
558 0 : }
559 :
560 :
561 : //_________________________________________________________________________
562 : Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
563 : AliPID::EParticleType species,
564 : ETPCgainScenario gainScenario,
565 : Double_t dEdx,
566 : Int_t nPoints,
567 : const TSpline3* responseFunction,
568 : Bool_t correctEta,
569 : Bool_t correctMultiplicity) const
570 : {
571 : // Calculates the expected sigma of the PID signal as the function of
572 : // the information stored in the track and the given parameters,
573 : // for the specified particle type
574 : //
575 :
576 : //if (!responseFunction)
577 : //return 999;
578 :
579 : //TODO Check whether it makes sense to set correctMultiplicity to kTRUE while correctEta might be kFALSE
580 :
581 : // If eta correction (=> new sigma parametrisation) is requested, but no sigma map is available, print error message
582 1206 : if (correctEta && !fhEtaSigmaPar1) {
583 0 : AliError("New sigma parametrisation requested, but sigma map not initialised (usually via AliPIDResponse). Old sigma parametrisation will be used!");
584 0 : }
585 :
586 : // If no sigma map is available or if no eta correction is requested (sigma maps only for corrected eta!), use the old parametrisation
587 1206 : if (!fhEtaSigmaPar1 || !correctEta) {
588 2412 : if (nPoints != 0)
589 3618 : return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity) *
590 2412 : fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints);
591 : else
592 0 : return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity)*fRes0[gainScenario];
593 : }
594 :
595 0 : if (nPoints > 0) {
596 : // Use eta correction (+ eta-dependent sigma)
597 0 : Double_t sigmaPar1 = GetSigmaPar1Fast(track, species, dEdx, responseFunction);
598 :
599 0 : if (correctMultiplicity) {
600 : // In addition, take into account multiplicity dependence of mean and sigma of dEdx
601 0 : Double_t dEdxExpectedEtaCorrected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
602 :
603 : // GetMultiplicityCorrection and GetMultiplicitySigmaCorrection both need the eta corrected dEdxExpected
604 0 : Double_t multiplicityCorrFactor = GetMultiplicityCorrectionFast(track, dEdxExpectedEtaCorrected, fCurrentEventMultiplicity);
605 0 : Double_t multiplicitySigmaCorrFactor = GetMultiplicitySigmaCorrectionFast(dEdxExpectedEtaCorrected, fCurrentEventMultiplicity);
606 :
607 : // multiplicityCorrFactor to correct dEdxExpected for multiplicity. In addition: Correction factor for sigma
608 0 : return (dEdxExpectedEtaCorrected * multiplicityCorrFactor)
609 0 : * (sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints) * multiplicitySigmaCorrFactor);
610 : }
611 : else {
612 0 : return GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE)*
613 0 : sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints);
614 : }
615 : }
616 : else {
617 : // One should never have/take tracks with 0 dEdx clusters!
618 0 : return 999;
619 : }
620 1206 : }
621 :
622 :
623 : //_________________________________________________________________________
624 : Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
625 : AliPID::EParticleType species,
626 : ETPCdEdxSource dedxSource,
627 : Bool_t correctEta,
628 : Bool_t correctMultiplicity) const
629 : {
630 : // Calculates the expected sigma of the PID signal as the function of
631 : // the information stored in the track, for the specified particle type
632 : // and dedx scenario
633 : //
634 :
635 1206 : Double_t dEdx = -1;
636 1206 : Int_t nPoints = -1;
637 1206 : ETPCgainScenario gainScenario = kGainScenarioInvalid;
638 1206 : TSpline3* responseFunction = 0x0;
639 :
640 1206 : if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
641 0 : return 999; //TODO: better handling!
642 :
643 1206 : return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity);
644 1206 : }
645 :
646 :
647 : //_________________________________________________________________________
648 : Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
649 : AliPID::EParticleType species,
650 : ETPCdEdxSource dedxSource,
651 : Bool_t correctEta,
652 : Bool_t correctMultiplicity) const
653 : {
654 : //Calculates the number of sigmas of the PID signal from the expected value
655 : //for a given particle species in the presence of multiple gain scenarios
656 : //inside the TPC
657 :
658 0 : Double_t dEdx = -1;
659 0 : Int_t nPoints = -1;
660 0 : ETPCgainScenario gainScenario = kGainScenarioInvalid;
661 0 : TSpline3* responseFunction = 0x0;
662 :
663 0 : if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
664 0 : return -999; //TODO: Better handling!
665 :
666 0 : Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity);
667 0 : Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity);
668 : // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters
669 0 : if (sigma >= 998)
670 0 : return -999;
671 : else
672 0 : return (dEdx-bethe)/sigma;
673 0 : }
674 :
675 : //_________________________________________________________________________
676 : Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track,
677 : AliPID::EParticleType species,
678 : ETPCdEdxSource dedxSource,
679 : Bool_t correctEta,
680 : Bool_t correctMultiplicity,
681 : Bool_t ratio/*=kFALSE*/)const
682 : {
683 : //Calculates the number of sigmas of the PID signal from the expected value
684 : //for a given particle species in the presence of multiple gain scenarios
685 : //inside the TPC
686 :
687 0 : Double_t dEdx = -1;
688 0 : Int_t nPoints = -1;
689 0 : ETPCgainScenario gainScenario = kGainScenarioInvalid;
690 0 : TSpline3* responseFunction = 0x0;
691 :
692 0 : if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
693 0 : return -9999.; //TODO: Better handling!
694 :
695 0 : const Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity);
696 :
697 : Double_t delta=-9999.;
698 0 : if (!ratio) delta=dEdx-bethe;
699 0 : else if (bethe>1.e-20) delta=dEdx/bethe;
700 :
701 0 : return delta;
702 0 : }
703 :
704 : //_________________________________________________________________________
705 : Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
706 : AliPID::EParticleType species,
707 : ETPCdEdxSource dedxSource,
708 : Double_t& dEdx,
709 : Int_t& nPoints,
710 : ETPCgainScenario& gainScenario,
711 : TSpline3** responseFunction) const
712 : {
713 : // Calculates the right parameters for PID
714 : // dEdx parametrization for the proper gain scenario, dEdx
715 : // and NPoints used for dEdx
716 : // based on the track geometry (which chambers it crosses) for the specified particle type
717 : // and preferred source of dedx.
718 : // returns true on success
719 :
720 :
721 2412 : if (dedxSource == kdEdxDefault) {
722 : // Fast handling for default case. In addition: Keep it simple (don't call additional functions) to
723 : // avoid possible bugs
724 :
725 : // GetTPCsignalTunedOnData will be non-positive, if it has not been set (i.e. in case of MC NOT tuned to data).
726 : // If this is the case, just take the normal signal
727 1206 : dEdx = track->GetTPCsignalTunedOnData();
728 1206 : if (dEdx <= 0) {
729 : // dEdx = track->GetTPCsignal();
730 1206 : dEdx = GetTrackdEdx(track);
731 1206 : }
732 :
733 1206 : nPoints = track->GetTPCsignalN();
734 1206 : gainScenario = kDefault;
735 :
736 1206 : TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario));
737 2412 : *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
738 :
739 : return kTRUE;
740 : }
741 :
742 : //TODO Proper handle of tuneMConData for other dEdx sources
743 :
744 0 : Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used)
745 0 : Char_t ncl[3]; //same
746 0 : Char_t nrows[3]; //same
747 0 : const AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo();
748 :
749 0 : if (!dEdxInfo && dedxSource!=kdEdxDefault) //in one case its ok if we dont have the info
750 : {
751 0 : AliError("AliTPCdEdxInfo not available");
752 0 : return kFALSE;
753 : }
754 :
755 0 : if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows);
756 :
757 : //check if we cross a bad OROC in which case we reject
758 0 : EChamberStatus trackOROCStatus = TrackStatus(track,2);
759 0 : if (trackOROCStatus==kChamberOff || trackOROCStatus==kChamberLowGain)
760 : {
761 0 : return kFALSE;
762 : }
763 :
764 0 : switch (dedxSource)
765 : {
766 : case kdEdxOROC:
767 : {
768 0 : if (trackOROCStatus==kChamberInvalid) return kFALSE; //never reached OROC
769 0 : dEdx = signal[3];
770 0 : nPoints = ncl[2]+ncl[1];
771 0 : gainScenario = kOROChigh;
772 0 : break;
773 : }
774 : case kdEdxHybrid:
775 : {
776 : //if we cross a bad IROC we use OROC dedx, if we dont we use combined
777 0 : EChamberStatus status = TrackStatus(track,1);
778 0 : if (status!=kChamberHighGain)
779 : {
780 0 : dEdx = signal[3];
781 0 : nPoints = ncl[2]+ncl[1];
782 0 : gainScenario = kOROChigh;
783 0 : }
784 : else
785 : {
786 : // dEdx = track->GetTPCsignal();
787 0 : dEdx = GetTrackdEdx(track);
788 0 : nPoints = track->GetTPCsignalN();
789 0 : gainScenario = kALLhigh;
790 : }
791 : break;
792 : }
793 : default:
794 : {
795 0 : dEdx = 0.;
796 0 : nPoints = 0;
797 0 : gainScenario = kGainScenarioInvalid;
798 0 : return kFALSE;
799 : }
800 : }
801 0 : TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario));
802 0 : *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
803 :
804 : return kTRUE;
805 1206 : }
806 :
807 :
808 : //_________________________________________________________________________
809 : Double_t AliTPCPIDResponse::GetEtaCorrectionFast(const AliVTrack *track, Double_t dEdxSplines) const
810 : {
811 : // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
812 : //
813 : // Get eta correction for the given parameters.
814 : //
815 :
816 0 : if (!fhEtaCorr) {
817 : // Calling this function means to request eta correction in some way. Print error message, if no map is available!
818 0 : AliError("Eta correction requested, but map not initialised (usually via AliPIDResponse). Returning eta correction factor 1!");
819 0 : return 1.;
820 : }
821 :
822 : Double_t tpcSignal = dEdxSplines;
823 :
824 0 : if (tpcSignal < 1.)
825 0 : return 1.;
826 :
827 0 : Double_t tanTheta = GetTrackTanTheta(track);
828 0 : Int_t binX = fhEtaCorr->GetXaxis()->FindFixBin(tanTheta);
829 0 : Int_t binY = fhEtaCorr->GetYaxis()->FindFixBin(1. / tpcSignal);
830 :
831 0 : if (binX == 0)
832 0 : binX = 1;
833 0 : if (binX > fhEtaCorr->GetXaxis()->GetNbins())
834 0 : binX = fhEtaCorr->GetXaxis()->GetNbins();
835 :
836 0 : if (binY == 0)
837 0 : binY = 1;
838 0 : if (binY > fhEtaCorr->GetYaxis()->GetNbins())
839 0 : binY = fhEtaCorr->GetYaxis()->GetNbins();
840 :
841 0 : return fhEtaCorr->GetBinContent(binX, binY);
842 0 : }
843 :
844 :
845 : //_________________________________________________________________________
846 : Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
847 : {
848 : //
849 : // Get eta correction for the given track.
850 : //
851 :
852 0 : if (!fhEtaCorr)
853 0 : return 1.;
854 :
855 0 : Double_t dEdx = -1;
856 0 : Int_t nPoints = -1;
857 0 : ETPCgainScenario gainScenario = kGainScenarioInvalid;
858 0 : TSpline3* responseFunction = 0x0;
859 :
860 0 : if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
861 0 : return 1.;
862 :
863 : // For the eta correction, do NOT take the multiplicity corrected value of dEdx
864 0 : Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
865 :
866 : //TODO Alternatively take current track dEdx
867 : //return GetEtaCorrectionFast(track, dEdx);
868 :
869 0 : return GetEtaCorrectionFast(track, dEdxSplines);
870 0 : }
871 :
872 :
873 : //_________________________________________________________________________
874 : Double_t AliTPCPIDResponse::GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
875 : {
876 : //
877 : // Get eta corrected dEdx for the given track. For the correction, the expected dEdx of
878 : // the specified species will be used. If the species is set to AliPID::kUnknown, the
879 : // dEdx of the track is used instead.
880 : // WARNING: In the latter case, the eta correction might not be as good as if the
881 : // expected dEdx is used, which is the way the correction factor is designed
882 : // for.
883 : // In any case, one has to decide carefully to which expected signal one wants to
884 : // compare the corrected value - to the corrected or uncorrected.
885 : // Anyhow, a safer way of looking e.g. at the n-sigma is to call
886 : // the corresponding function GetNumberOfSigmas!
887 : //
888 :
889 0 : Double_t dEdx = -1;
890 0 : Int_t nPoints = -1;
891 0 : ETPCgainScenario gainScenario = kGainScenarioInvalid;
892 0 : TSpline3* responseFunction = 0x0;
893 :
894 : // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case
895 : // it is not used anyway, so this causes no trouble.
896 0 : if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
897 0 : return -1.;
898 :
899 : Double_t etaCorr = 0.;
900 :
901 0 : if (species < AliPID::kUnknown) {
902 : // For the eta correction, do NOT take the multiplicity corrected value of dEdx
903 0 : Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
904 0 : etaCorr = GetEtaCorrectionFast(track, dEdxSplines);
905 0 : }
906 : else {
907 0 : etaCorr = GetEtaCorrectionFast(track, dEdx);
908 : }
909 :
910 0 : if (etaCorr <= 0)
911 0 : return -1.;
912 :
913 0 : return dEdx / etaCorr;
914 0 : }
915 :
916 :
917 : //_________________________________________________________________________
918 : Double_t AliTPCPIDResponse::GetSigmaPar1Fast(const AliVTrack *track, AliPID::EParticleType species, Double_t dEdx,
919 : const TSpline3* responseFunction) const
920 : {
921 : // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
922 : //
923 : // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track.
924 : //
925 :
926 0 : if (!fhEtaSigmaPar1) {
927 : // Calling this function means to request new sigma parametrisation in some way. Print error message, if no map is available!
928 0 : AliError("New sigma parametrisation requested, but sigma map not initialised (usually via AliPIDResponse). Returning error value for sigma parameter1 = 999!");
929 0 : return 999;
930 : }
931 :
932 : // The sigma maps are created with data sets that are already eta corrected and for which the
933 : // splines have been re-created. Therefore, the value for the lookup needs to be the value of
934 : // the splines without any additional eta correction.
935 : // NOTE: This is due to the method the maps are created. The track dEdx (not the expected one!)
936 : // is corrected to uniquely related a momemtum bin with an expected dEdx, where the expected dEdx
937 : // equals the track dEdx for all eta (thanks to the correction and the re-fit of the splines).
938 : // Consequently, looking up the uncorrected expected dEdx at a given tanTheta yields the correct
939 : // sigma parameter!
940 : // Also: It has to be the spline dEdx, since one wants to get the sigma for the assumption(!)
941 : // of such a particle, which by assumption then has this dEdx value
942 :
943 : // For the eta correction, do NOT take the multiplicity corrected value of dEdx
944 0 : Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
945 :
946 0 : if (dEdxExpected < 1.)
947 0 : return 999;
948 :
949 0 : Double_t tanTheta = GetTrackTanTheta(track);
950 0 : Int_t binX = fhEtaSigmaPar1->GetXaxis()->FindFixBin(tanTheta);
951 0 : Int_t binY = fhEtaSigmaPar1->GetYaxis()->FindFixBin(1. / dEdxExpected);
952 :
953 0 : if (binX == 0)
954 0 : binX = 1;
955 0 : if (binX > fhEtaSigmaPar1->GetXaxis()->GetNbins())
956 0 : binX = fhEtaSigmaPar1->GetXaxis()->GetNbins();
957 :
958 0 : if (binY == 0)
959 0 : binY = 1;
960 0 : if (binY > fhEtaSigmaPar1->GetYaxis()->GetNbins())
961 0 : binY = fhEtaSigmaPar1->GetYaxis()->GetNbins();
962 :
963 0 : return fhEtaSigmaPar1->GetBinContent(binX, binY);
964 0 : }
965 :
966 :
967 : //_________________________________________________________________________
968 : Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
969 : {
970 : //
971 : // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track.
972 : //
973 :
974 0 : if (!fhEtaSigmaPar1)
975 0 : return 999;
976 :
977 0 : Double_t dEdx = -1;
978 0 : Int_t nPoints = -1;
979 0 : ETPCgainScenario gainScenario = kGainScenarioInvalid;
980 0 : TSpline3* responseFunction = 0x0;
981 :
982 0 : if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
983 0 : return 999;
984 :
985 0 : return GetSigmaPar1Fast(track, species, dEdx, responseFunction);
986 0 : }
987 :
988 :
989 : //_________________________________________________________________________
990 : Bool_t AliTPCPIDResponse::SetEtaCorrMap(TH2D* hMap)
991 : {
992 : //
993 : // Load map for TPC eta correction (a copy is stored and will be deleted automatically).
994 : // If hMap is 0x0,the eta correction will be disabled and kFALSE is returned.
995 : // If the map can be set, kTRUE is returned.
996 : //
997 :
998 0 : delete fhEtaCorr;
999 :
1000 0 : if (!hMap) {
1001 0 : fhEtaCorr = 0x0;
1002 :
1003 0 : return kFALSE;
1004 : }
1005 :
1006 0 : fhEtaCorr = (TH2D*)(hMap->Clone());
1007 0 : fhEtaCorr->SetDirectory(0);
1008 :
1009 0 : return kTRUE;
1010 0 : }
1011 :
1012 :
1013 : //_________________________________________________________________________
1014 : Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0)
1015 : {
1016 : //
1017 : // Load map for TPC sigma map (a copy is stored and will be deleted automatically):
1018 : // Parameter 1 is stored as a 2D map (1/dEdx vs. tanTheta_local) and parameter 0 is
1019 : // a constant. If hSigmaPar1Map is 0x0, the old sigma parametrisation will be used
1020 : // (and sigmaPar0 is ignored!) and kFALSE is returned.
1021 : // If the map can be set, sigmaPar0 is also set and kTRUE will be returned.
1022 : //
1023 :
1024 0 : delete fhEtaSigmaPar1;
1025 :
1026 0 : if (!hSigmaPar1Map) {
1027 0 : fhEtaSigmaPar1 = 0x0;
1028 0 : fSigmaPar0 = 0.0;
1029 :
1030 0 : return kFALSE;
1031 : }
1032 :
1033 0 : fhEtaSigmaPar1 = (TH2D*)(hSigmaPar1Map->Clone());
1034 0 : fhEtaSigmaPar1->SetDirectory(0);
1035 0 : fSigmaPar0 = sigmaPar0;
1036 :
1037 0 : return kTRUE;
1038 0 : }
1039 :
1040 :
1041 : //_________________________________________________________________________
1042 : Double_t AliTPCPIDResponse::GetTrackTanTheta(const AliVTrack *track) const
1043 : {
1044 : // Extract the tanTheta from the information available in the AliVTrack
1045 :
1046 : // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()).
1047 : // However, this value is not available for AODs and, thus, not for AliVTrack.
1048 : // Fortunately, the following formula allows to approximate the local tanTheta with the
1049 : // global theta angle -> This is for by far most of the tracks the same, but gives at
1050 : // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok.
1051 :
1052 : /*
1053 : const AliExternalTrackParam* innerParam = track->GetInnerParam();
1054 : Double_t tanTheta = 0;
1055 : if (innerParam)
1056 : tanTheta = innerParam->GetTgl();
1057 : else
1058 : tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
1059 :
1060 : // Constant in formula for B in kGauss (factor 0.1 to convert B from Tesla to kGauss),
1061 : // 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)
1062 : const Double_t constant = TMath::C()* 1e-9 * 0.1 * 0.01;
1063 : const Double_t curvature = fMagField * constant / track->Pt(); // in 1./cm
1064 :
1065 : Double_t averageddzdr = 0.;
1066 : Int_t nParts = 0;
1067 :
1068 : for (Double_t r = 85; r < 245; r++) {
1069 : Double_t sinPhiLocal = TMath::Abs(r*curvature*0.5);
1070 :
1071 : // Cut on |sin(phi)| as during reco
1072 : if (TMath::Abs(sinPhiLocal) <= 0.95) {
1073 : const Double_t phiLocal = TMath::ASin(sinPhiLocal);
1074 : const Double_t tanPhiLocal = TMath::Tan(phiLocal);
1075 :
1076 : averageddzdr += tanTheta * TMath::Sqrt(1. + tanPhiLocal * tanPhiLocal);
1077 : nParts++;
1078 : }
1079 : }
1080 :
1081 : if (nParts > 0)
1082 : averageddzdr /= nParts;
1083 : else {
1084 : AliError("Problems during determination of dz/dr. Returning pure tanTheta as best estimate!");
1085 : return tanTheta;
1086 : }
1087 :
1088 : //printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\ntanThetaGlobalFromTheta/tanTheta/Averageddzdr: %f / %f / %f\n\n",
1089 : // track->Pt(), constant, fMagField, 1./curvature, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, averageddzdr);
1090 :
1091 : return averageddzdr;
1092 : */
1093 :
1094 :
1095 : // Alternatively (in average, the effect is found to be negligable!):
1096 : // Take local tanTheta from TPC inner wall, if available (currently only for ESDs available)
1097 : //const AliExternalTrackParam* innerParam = track->GetInnerParam();
1098 : //if (innerParam) {
1099 : // return innerParam->GetTgl();
1100 : //}
1101 :
1102 0 : return TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
1103 : }
1104 :
1105 :
1106 : //_________________________________________________________________________
1107 : Double_t AliTPCPIDResponse::GetMultiplicityCorrectionFast(const AliVTrack *track, Double_t dEdxExpected, Int_t multiplicity) const
1108 : {
1109 : // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
1110 : //
1111 : // Calculate the multiplicity correction factor for this track for the given multiplicity.
1112 : // The parameter dEdxExpected should take into account the eta correction already!
1113 :
1114 : // Multiplicity depends on pure dEdx. Therefore, correction factor depends indirectly on eta
1115 : // => Use eta corrected value for dEdexExpected.
1116 :
1117 0 : if (dEdxExpected <= 0 || multiplicity <= 0)
1118 0 : return 1.0;
1119 :
1120 0 : const Double_t dEdxExpectedInv = 1. / dEdxExpected;
1121 0 : Double_t relSlope = fCorrFuncMultiplicity->Eval(dEdxExpectedInv);
1122 :
1123 0 : const Double_t tanTheta = GetTrackTanTheta(track);
1124 0 : relSlope += fCorrFuncMultiplicityTanTheta->Eval(tanTheta);
1125 :
1126 0 : return (1. + relSlope * multiplicity);
1127 0 : }
1128 :
1129 :
1130 : //_________________________________________________________________________
1131 : Double_t AliTPCPIDResponse::GetMultiplicityCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
1132 : {
1133 : //
1134 : // Get multiplicity correction for the given track (w.r.t. the multiplicity of the current event)
1135 : //
1136 :
1137 : //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
1138 :
1139 : // No correction in case of multiplicity <= 0
1140 0 : if (fCurrentEventMultiplicity <= 0)
1141 0 : return 1.;
1142 :
1143 0 : Double_t dEdx = -1;
1144 0 : Int_t nPoints = -1;
1145 0 : ETPCgainScenario gainScenario = kGainScenarioInvalid;
1146 0 : TSpline3* responseFunction = 0x0;
1147 :
1148 0 : if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
1149 0 : return 1.;
1150 :
1151 : //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid?
1152 : // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course)
1153 0 : Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
1154 :
1155 0 : return GetMultiplicityCorrectionFast(track, dEdxExpected, fCurrentEventMultiplicity);
1156 0 : }
1157 :
1158 :
1159 : //_________________________________________________________________________
1160 : Double_t AliTPCPIDResponse::GetMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
1161 : {
1162 : //
1163 : // Get multiplicity corrected dEdx for the given track. For the correction, the expected dEdx of
1164 : // the specified species will be used. If the species is set to AliPID::kUnknown, the
1165 : // dEdx of the track is used instead.
1166 : // WARNING: In the latter case, the correction might not be as good as if the
1167 : // expected dEdx is used, which is the way the correction factor is designed
1168 : // for.
1169 : // In any case, one has to decide carefully to which expected signal one wants to
1170 : // compare the corrected value - to the corrected or uncorrected.
1171 : // Anyhow, a safer way of looking e.g. at the n-sigma is to call
1172 : // the corresponding function GetNumberOfSigmas!
1173 : //
1174 :
1175 : //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
1176 :
1177 0 : Double_t dEdx = -1;
1178 0 : Int_t nPoints = -1;
1179 0 : ETPCgainScenario gainScenario = kGainScenarioInvalid;
1180 0 : TSpline3* responseFunction = 0x0;
1181 :
1182 : // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case
1183 : // it is not used anyway, so this causes no trouble.
1184 0 : if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
1185 0 : return -1.;
1186 :
1187 :
1188 : // No correction in case of multiplicity <= 0
1189 0 : if (fCurrentEventMultiplicity <= 0)
1190 0 : return dEdx;
1191 :
1192 : Double_t multiplicityCorr = 0.;
1193 :
1194 : // 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?
1195 0 : if (species < AliPID::kUnknown) {
1196 : // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course).
1197 : // However, one needs the eta corrected value!
1198 0 : Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
1199 0 : multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdxSplines, fCurrentEventMultiplicity);
1200 0 : }
1201 : else {
1202 : // One needs the eta corrected value to determine the multiplicity correction factor!
1203 0 : Double_t etaCorr = GetEtaCorrectionFast(track, dEdx);
1204 0 : multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdx * etaCorr, fCurrentEventMultiplicity);
1205 : }
1206 :
1207 0 : if (multiplicityCorr <= 0)
1208 0 : return -1.;
1209 :
1210 0 : return dEdx / multiplicityCorr;
1211 0 : }
1212 :
1213 :
1214 : //_________________________________________________________________________
1215 : Double_t AliTPCPIDResponse::GetEtaAndMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species,
1216 : ETPCdEdxSource dedxSource) const
1217 : {
1218 : //
1219 : // Get multiplicity and eta corrected dEdx for the given track. For the correction,
1220 : // the expected dEdx of the specified species will be used. If the species is set
1221 : // to AliPID::kUnknown, the dEdx of the track is used instead.
1222 : // WARNING: In the latter case, the correction might not be as good as if the
1223 : // expected dEdx is used, which is the way the correction factor is designed
1224 : // for.
1225 : // In any case, one has to decide carefully to which expected signal one wants to
1226 : // compare the corrected value - to the corrected or uncorrected.
1227 : // Anyhow, a safer way of looking e.g. at the n-sigma is to call
1228 : // the corresponding function GetNumberOfSigmas!
1229 : //
1230 :
1231 0 : Double_t dEdx = -1;
1232 0 : Int_t nPoints = -1;
1233 0 : ETPCgainScenario gainScenario = kGainScenarioInvalid;
1234 0 : TSpline3* responseFunction = 0x0;
1235 :
1236 : // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case
1237 : // it is not used anyway, so this causes no trouble.
1238 0 : if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
1239 0 : return -1.;
1240 :
1241 : Double_t multiplicityCorr = 0.;
1242 : Double_t etaCorr = 0.;
1243 :
1244 0 : if (species < AliPID::kUnknown) {
1245 : // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course)
1246 0 : Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
1247 0 : etaCorr = GetEtaCorrectionFast(track, dEdxSplines);
1248 0 : multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdxSplines * etaCorr, fCurrentEventMultiplicity);
1249 0 : }
1250 : else {
1251 0 : etaCorr = GetEtaCorrectionFast(track, dEdx);
1252 0 : multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdx * etaCorr, fCurrentEventMultiplicity);
1253 : }
1254 :
1255 0 : if (multiplicityCorr <= 0 || etaCorr <= 0)
1256 0 : return -1.;
1257 :
1258 0 : return dEdx / multiplicityCorr / etaCorr;
1259 0 : }
1260 :
1261 :
1262 : //_________________________________________________________________________
1263 : Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrectionFast(Double_t dEdxExpected, Int_t multiplicity) const
1264 : {
1265 : // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
1266 : //
1267 : // Calculate the multiplicity sigma correction factor for the corresponding expected dEdx and for the given multiplicity.
1268 : // The parameter dEdxExpected should take into account the eta correction already!
1269 :
1270 : // Multiplicity dependence of sigma depends on the real dEdx at zero multiplicity,
1271 : // i.e. the eta (only) corrected dEdxExpected value has to be used
1272 : // since all maps etc. have been created for ~zero multiplicity
1273 :
1274 0 : if (dEdxExpected <= 0 || multiplicity <= 0)
1275 0 : return 1.0;
1276 :
1277 0 : Double_t relSigmaSlope = fCorrFuncSigmaMultiplicity->Eval(1. / dEdxExpected);
1278 :
1279 0 : return (1. + relSigmaSlope * multiplicity);
1280 0 : }
1281 :
1282 :
1283 : //_________________________________________________________________________
1284 : Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
1285 : {
1286 : //
1287 : // Get multiplicity sigma correction for the given track (w.r.t. the multiplicity of the current event)
1288 : //
1289 :
1290 : //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
1291 :
1292 : // No correction in case of multiplicity <= 0
1293 0 : if (fCurrentEventMultiplicity <= 0)
1294 0 : return 1.;
1295 :
1296 0 : Double_t dEdx = -1;
1297 0 : Int_t nPoints = -1;
1298 0 : ETPCgainScenario gainScenario = kGainScenarioInvalid;
1299 0 : TSpline3* responseFunction = 0x0;
1300 :
1301 0 : if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
1302 0 : return 1.;
1303 :
1304 : //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid?
1305 : // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course)
1306 0 : Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
1307 :
1308 0 : return GetMultiplicitySigmaCorrectionFast(dEdxExpected, fCurrentEventMultiplicity);
1309 0 : }
1310 :
1311 : //_____________________________________________________________________________
1312 : Double_t AliTPCPIDResponse::GetTrackdEdx(const AliVTrack* track) const
1313 : {
1314 : // get dEdx of the track depending on the selected dEdx type
1315 2412 : if (fdEdxType == kdEdxInfo) {
1316 0 : const AliTPCdEdxInfo *dEdxInfo = track->GetTPCdEdxInfo();
1317 0 : if (dEdxInfo) {
1318 0 : return dEdxInfo->GetWeightedMean(fdEdxChargeType, fdEdxWeightType, fIROCweight, fOROCmedWeight, fOROClongWeight);
1319 : } else {
1320 0 : AliError("Could not find dEdx info, using default signal");
1321 : }
1322 0 : }
1323 :
1324 1206 : return track->GetTPCsignal();
1325 1206 : }
1326 :
1327 :
1328 : //_________________________________________________________________________
1329 : void AliTPCPIDResponse::ResetMultiplicityCorrectionFunctions()
1330 : {
1331 : // Default values: No correction, i.e. overall correction factor should be one
1332 :
1333 : // This function should always return 0.0, if multcorr disabled
1334 20 : fCorrFuncMultiplicity->SetParameter(0, 0.);
1335 10 : fCorrFuncMultiplicity->SetParameter(1, 0.);
1336 10 : fCorrFuncMultiplicity->SetParameter(2, 0.);
1337 10 : fCorrFuncMultiplicity->SetParameter(3, 0.);
1338 10 : fCorrFuncMultiplicity->SetParameter(4, 0.);
1339 :
1340 : // This function should always return 0., if multCorr disabled
1341 10 : fCorrFuncMultiplicityTanTheta->SetParameter(0, 0.);
1342 10 : fCorrFuncMultiplicityTanTheta->SetParameter(1, 0.);
1343 10 : fCorrFuncMultiplicityTanTheta->SetParameter(2, 0.);
1344 :
1345 : // This function should always return 0.0, if mutlcorr disabled
1346 10 : fCorrFuncSigmaMultiplicity->SetParameter(0, 0.);
1347 10 : fCorrFuncSigmaMultiplicity->SetParameter(1, 0.);
1348 10 : fCorrFuncSigmaMultiplicity->SetParameter(2, 0.);
1349 10 : fCorrFuncSigmaMultiplicity->SetParameter(3, 0.);
1350 10 : }
1351 :
1352 :
1353 : //_________________________________________________________________________
1354 : Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner,
1355 : Double_t* trackPositionOuter,
1356 : Float_t& inphi,
1357 : Float_t& outphi,
1358 : Int_t& in,
1359 : Int_t& out ) const
1360 : {
1361 : //calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses
1362 : //for OROC chamber numbers add 36
1363 : //returned angles are between (0,2pi)
1364 :
1365 0 : inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]);
1366 0 : outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]);
1367 :
1368 0 : if (inphi<0) {inphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
1369 0 : if (outphi<0) {outphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
1370 :
1371 0 : in = sectorNumber(inphi);
1372 0 : out = sectorNumber(outphi);
1373 :
1374 : //for the C side (positive z) offset by 18
1375 0 : if (trackPositionInner[2]>0.0)
1376 : {
1377 0 : in+=18;
1378 0 : out+=18;
1379 0 : }
1380 0 : return kTRUE;
1381 : }
1382 :
1383 :
1384 : //_____________________________________________________________________________
1385 : Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const
1386 : {
1387 : //calculate sector number
1388 0 : const Float_t width=TMath::TwoPi()/18.;
1389 0 : return TMath::Floor(phi/width);
1390 : }
1391 :
1392 :
1393 : //_____________________________________________________________________________
1394 : void AliTPCPIDResponse::Print(Option_t* /*option*/) const
1395 : {
1396 : //Print info
1397 0 : fResponseFunctions.Print();
1398 0 : }
1399 :
1400 :
1401 : //_____________________________________________________________________________
1402 : AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const
1403 : {
1404 : //status of the track: if it crosses any bad chambers return kChamberOff;
1405 : //IROC:layer=1, OROC:layer=2
1406 0 : if (layer<1 || layer>2) layer=1;
1407 0 : Int_t in=0;
1408 0 : Int_t out=0;
1409 0 : Float_t inphi=0.;
1410 0 : Float_t outphi=0.;
1411 0 : Float_t innerRadius = (layer==1)?83.0:133.7;
1412 0 : Float_t outerRadius = (layer==1)?133.5:247.7;
1413 :
1414 : /////////////////////////////////////////////////////////////////////////////
1415 : //find out where track enters and leaves the layer.
1416 : //
1417 0 : Double_t trackPositionInner[3];
1418 0 : Double_t trackPositionOuter[3];
1419 :
1420 : //if there is no inner param this could mean we're using the AOD track,
1421 : //we still can extrapolate from the vertex - so use those params.
1422 0 : const AliExternalTrackParam* ip = track->GetInnerParam();
1423 0 : if (ip) track=dynamic_cast<const AliVTrack*>(ip);
1424 :
1425 0 : Bool_t trackAtInner = track->GetXYZAt(innerRadius, fMagField, trackPositionInner);
1426 0 : Bool_t trackAtOuter = track->GetXYZAt(outerRadius, fMagField, trackPositionOuter);
1427 :
1428 0 : if (!trackAtInner)
1429 : {
1430 : //if we dont even enter inner radius we do nothing and return invalid
1431 0 : inphi=0.0;
1432 0 : outphi=0.0;
1433 0 : in=0;
1434 0 : out=0;
1435 0 : return kChamberInvalid;
1436 : }
1437 :
1438 0 : if (!trackAtOuter)
1439 : {
1440 : //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position
1441 0 : Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter);
1442 0 : Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]);
1443 0 : if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius)
1444 : {
1445 : //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]);
1446 : }
1447 : else
1448 : {
1449 0 : inphi=0.0;
1450 0 : outphi=0.0;
1451 0 : in=0;
1452 0 : out=0;
1453 0 : return kChamberInvalid;
1454 : }
1455 0 : }
1456 :
1457 :
1458 0 : if (!sectorNumbersInOut(trackPositionInner,
1459 : trackPositionOuter,
1460 : inphi,
1461 : outphi,
1462 : in,
1463 0 : out)) return kChamberInvalid;
1464 :
1465 : /////////////////////////////////////////////////////////////////////////////
1466 : //now we have the location of the track we can check
1467 : //if it is in a good/bad chamber
1468 : //
1469 : Bool_t sideA = kTRUE;
1470 :
1471 0 : if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE;
1472 :
1473 0 : in=in%18;
1474 0 : out=out%18;
1475 :
1476 0 : if (TMath::Abs(in-out)>9)
1477 : {
1478 0 : if (TMath::Max(in,out)==out)
1479 : {
1480 : Int_t tmp=out;
1481 0 : out = in;
1482 0 : in = tmp;
1483 0 : Float_t tmpphi=outphi;
1484 0 : outphi=inphi;
1485 0 : inphi=tmpphi;
1486 0 : }
1487 0 : in-=18;
1488 0 : inphi-=TMath::TwoPi();
1489 0 : }
1490 : else
1491 : {
1492 0 : if (TMath::Max(in,out)==in)
1493 : {
1494 0 : Int_t tmp=out;
1495 0 : out=in;
1496 0 : in=tmp;
1497 0 : Float_t tmpphi=outphi;
1498 0 : outphi=inphi;
1499 0 : inphi=tmpphi;
1500 0 : }
1501 : }
1502 :
1503 : Float_t trackLengthInBad = 0.;
1504 : Float_t trackLengthInLowGain = 0.;
1505 0 : Float_t trackLengthTotal = TMath::Abs(outphi-inphi);
1506 : Float_t lengthFractionInBadSectors = 0.;
1507 :
1508 0 : const Float_t sectorWidth = TMath::TwoPi()/18.;
1509 :
1510 0 : for (Int_t i=in; i<=out; i++)
1511 : {
1512 : int j=i;
1513 0 : if (i<0) j+=18; //correct for the negative values
1514 0 : if (!sideA) j+=18; //move to the correct side
1515 :
1516 : Float_t deltaPhi = 0.;
1517 0 : Float_t phiEdge=sectorWidth*i;
1518 0 : if (inphi>phiEdge) {deltaPhi=phiEdge+sectorWidth-inphi;}
1519 0 : else if ((outphi>=phiEdge) && (outphi<(phiEdge+sectorWidth))) {deltaPhi=outphi-phiEdge;}
1520 : else {deltaPhi=sectorWidth;}
1521 :
1522 0 : Float_t v = fVoltageMap[(layer==1)?(j):(j+36)];
1523 0 : if (v<=fBadIROCthreshhold)
1524 : {
1525 0 : trackLengthInBad+=deltaPhi;
1526 : lengthFractionInBadSectors=1.;
1527 0 : }
1528 0 : if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold)
1529 : {
1530 0 : trackLengthInLowGain+=deltaPhi;
1531 : lengthFractionInBadSectors=1.;
1532 0 : }
1533 : }
1534 :
1535 : //for now low gain and bad (off) chambers are treated equally
1536 0 : if (trackLengthTotal>0)
1537 0 : lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal;
1538 :
1539 : //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);
1540 :
1541 0 : if (lengthFractionInBadSectors>fMaxBadLengthFraction)
1542 : {
1543 : //printf("%%%%%%%% %s kChamberLowGain\n",(layer==1)?"IROC":"OROC");
1544 0 : return kChamberLowGain;
1545 : }
1546 :
1547 0 : return kChamberHighGain;
1548 0 : }
1549 :
1550 :
1551 : //_____________________________________________________________________________
1552 : Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const
1553 : {
1554 : //return the radius of the outermost padrow containing a cluster in TPC
1555 : //for the track
1556 0 : const TBits* clusterMap=track->GetTPCClusterMapPtr();
1557 0 : if (!clusterMap) return 0.;
1558 :
1559 : //from AliTPCParam, radius of first IROC row
1560 : const Float_t rfirstIROC = 8.52249984741210938e+01;
1561 : const Float_t padrowHeightIROC = 0.75;
1562 : const Float_t rfirstOROC0 = 1.35100006103515625e+02;
1563 : const Float_t padrowHeightOROC0 = 1.0;
1564 : const Float_t rfirstOROC1 = 1.99350006103515625e+02;
1565 : const Float_t padrowHeightOROC1 = 1.5;
1566 :
1567 : Int_t maxPadRow=160;
1568 0 : while ((--maxPadRow)>0 && !clusterMap->TestBitNumber(maxPadRow)){}
1569 0 : if (maxPadRow>126) return rfirstOROC1+(maxPadRow-126-1)*padrowHeightOROC1;
1570 0 : if (maxPadRow>62) return rfirstOROC0+(maxPadRow-62-1)*padrowHeightOROC0;
1571 0 : if (maxPadRow>0) return rfirstIROC+(maxPadRow-1)*padrowHeightIROC;
1572 0 : return 0.0;
1573 0 : }
1574 :
1575 :
1576 : //_____________________________________________________________________________
1577 : Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const
1578 : {
1579 : //calculate the coordinates of the apex of the track
1580 0 : Double_t x[3];
1581 0 : track->GetXYZ(x);
1582 0 : Double_t p[3];
1583 0 : track->GetPxPyPz(p);
1584 0 : Double_t r = 1./track->OneOverPt()/0.0299792458/magField; //signed - will determine the direction of b
1585 : //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);
1586 : //find orthogonal vector (Gram-Schmidt)
1587 0 : Double_t alpha = (p[0]*x[0] + p[1]*x[1])/(p[0]*p[0] + p[1]*p[1]);
1588 : Double_t b[2];
1589 0 : b[0]=x[0]-alpha*p[0];
1590 0 : b[1]=x[1]-alpha*p[1];
1591 :
1592 0 : Double_t norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
1593 0 : if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
1594 0 : b[0]/=norm;
1595 0 : b[1]/=norm;
1596 0 : b[0]*=r;
1597 0 : b[1]*=r;
1598 0 : b[0]+=x[0];
1599 0 : b[1]+=x[1];
1600 : //printf("center: x:%.2f, y:%.2f\n",b[0],b[1]);
1601 :
1602 0 : norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
1603 0 : if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
1604 :
1605 0 : position[0]=b[0]+b[0]*TMath::Abs(r)/norm;
1606 0 : position[1]=b[1]+b[1]*TMath::Abs(r)/norm;
1607 0 : position[2]=0.;
1608 0 : return kTRUE;
1609 0 : }
1610 :
1611 : Double_t AliTPCPIDResponse::EvaldEdxSpline(Double_t bg,Int_t entry){
1612 : //
1613 : // Evaluate the dEdx response for given entry
1614 : //
1615 0 : TSpline * spline = (TSpline*)fSplineArray.At(entry);
1616 0 : if (spline) return spline->Eval(bg);
1617 0 : return 0;
1618 0 : }
1619 :
1620 :
1621 : Bool_t AliTPCPIDResponse::RegisterSpline(const char * name, Int_t index){
1622 : //
1623 : // register spline to be used for drawing comparisons
1624 : //
1625 0 : TFile * fTPCBB = TFile::Open("$ALICE_PHYSICS/OADB/COMMON/PID/data/TPCPIDResponse.root");
1626 0 : TObjArray *arrayTPCPID= (TObjArray*) fTPCBB->Get("TPCPIDResponse");
1627 0 : if (fSplineArray.GetEntriesFast()<index) fSplineArray.Expand(index*2);
1628 : TSpline3 *spline=0;
1629 0 : if (arrayTPCPID){
1630 0 : spline = (TSpline3*)arrayTPCPID->FindObject(name);
1631 0 : if (spline) fSplineArray.AddAt(spline->Clone(),index);
1632 : }
1633 0 : delete arrayTPCPID;
1634 0 : delete fTPCBB;
1635 0 : return (spline!=0);
1636 : }
1637 :
1638 : //_____________________________________________________________________________
1639 : Bool_t AliTPCPIDResponse::InitFromOADB(const Int_t run, const char* pass, const char* oadbFile/*="$ALICE_PHYSICS/OADB/COMMON/PID/data/TPCPIDResponseOADB.root"*/, Bool_t initMultiplicityCorrection/*=kTRUE*/)
1640 : {
1641 : //
1642 : //
1643 : //
1644 :
1645 0 : AliInfo( "----------------------| Initialisation TPC PID Response from OADB |----------------------");
1646 0 : AliInfo(TString::Format("----------------------| Run: %d, pass: %-19s |----------------------", run, pass));
1647 0 : if (!fOADBContainer) {
1648 0 : fOADBContainer = new AliOADBContainer("TPCSplines");
1649 0 : fOADBContainer->InitFromFile(oadbFile,"TPCSplines");
1650 0 : }
1651 :
1652 0 : const TString spass(pass);
1653 0 : const Int_t passOrig=spass.Atoi();
1654 : Int_t passIter=passOrig;
1655 :
1656 0 : const TObjArray *arr=dynamic_cast<TObjArray*>(fOADBContainer->GetObject(run,"",spass.Data()));
1657 0 : while (!arr && --passIter > 0) {
1658 0 : TString passCurrent=TString::Format("%d", passIter);
1659 0 : arr=dynamic_cast<TObjArray*>(fOADBContainer->GetObject(run,"",passCurrent.Data()));
1660 0 : }
1661 :
1662 0 : if (!arr) {
1663 0 : AliError("***** Risk for unreliable TPC PID detected: ********");
1664 0 : AliError(Form(" could not find a valid OADB object for run %d pass %s", run, pass));
1665 0 : AliError( " Most probably this is because the PID response for this pass is not available, yet");
1666 0 : AliError( " Please also check https://twiki.cern.ch/twiki/bin/view/ALICE/TPCSplines");
1667 0 : return kFALSE;
1668 : }
1669 :
1670 0 : if (passIter<passOrig) {
1671 0 : AliWarning( "***** Risk for unreliable TPC PID detected: ********");
1672 0 : AliWarning(Form(" Using splines from a previous pass: %d<%d", passIter, passOrig));
1673 0 : AliWarning( " Most probably this is because the PID response for this pass has not been extracted, yet");
1674 0 : AliWarning( " Please also check https://twiki.cern.ch/twiki/bin/view/ALICE/TPCSplines");
1675 : }
1676 :
1677 : //===| Set up of splines |====================================================
1678 : // clear response functions and reset spline usage
1679 0 : fResponseFunctions.Clear();
1680 0 : SetUseDatabase(kFALSE);
1681 :
1682 0 : const TObjArray *arrSplines = static_cast<TObjArray*>(arr->FindObject("Splines"));
1683 0 : if (!arrSplines) {
1684 0 : AliError("***** Risk for unreliable TPC PID detected: ********");
1685 0 : AliError(Form(" could not find array of splines for run %d", run));
1686 0 : AliError( " This should not happen, plese report");
1687 0 : return kFALSE;
1688 : }
1689 0 : SetSplinesFromArray(arrSplines);
1690 :
1691 : //===| Set up multiplicity correction |=======================================
1692 0 : if (initMultiplicityCorrection) {
1693 0 : const TObject *multiplicityCorrection=arr->FindObject("MultiplicityCorrection");
1694 0 : if (multiplicityCorrection) {
1695 0 : const TString multiplicityData(multiplicityCorrection->GetTitle());
1696 0 : const Bool_t res=SetMultiplicityCorrectionFromString(multiplicityData);
1697 0 : if (!res) {
1698 0 : AliError("Problem setting up multiplicity correction for TPC PID");
1699 : }
1700 0 : }
1701 0 : } else {
1702 0 : AliInfo("Multiplicity correction explicitly disabled");
1703 : }
1704 :
1705 : //===| Set up of dEdx type |==================================================
1706 0 : const TNamed *dEdxType=static_cast<TNamed*>(arr->FindObject("dEdxType"));
1707 0 : if (dEdxType) {
1708 0 : const TString dEdxTypeSet(dEdxType->GetTitle());
1709 0 : SetdEdxTypeFromString(dEdxTypeSet);
1710 0 : }
1711 :
1712 : //===| resolution parametrisation |===========================================
1713 0 : const TNamed *dEdxResolution=static_cast<TNamed*>(arr->FindObject("dEdxResolution"));
1714 0 : if (dEdxResolution) {
1715 0 : const TString dEdxResolutionString(dEdxResolution->GetTitle());
1716 0 : SetdEdxResolutionFromString(dEdxResolutionString);
1717 0 : }
1718 :
1719 : return kTRUE;
1720 0 : }
1721 :
1722 : //______________________________________________________________________________
1723 : Bool_t AliTPCPIDResponse::SetSplinesFromArray(const TObjArray* arrSplines)
1724 : {
1725 : // Set up internal spline array from order array of splines 'arrSplines'
1726 : // arrSplines is assumes to have the splines for the single particles inc
1727 : // the numbered position as given by AliPID::EParticleType
1728 :
1729 : //---| get default splines for missing species |------------------------------
1730 0 : TObject *protonSpline = arrSplines->At(AliPID::kProton );
1731 0 : TObject *pionSpline = arrSplines->At(AliPID::kPion );
1732 : // TObject *allSpline = arrSplines->At(AliPID::kSPECIESC);
1733 :
1734 : //---| set up the spline array |----------------------------------------------
1735 0 : for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie) {
1736 0 : TSpline3 *responseFunction = static_cast<TSpline3*>(arrSplines->At(ispecie));
1737 :
1738 : //---| spline replacement |-------------------------------------------------
1739 : // usually no splines are extracted for muons and light nuclei
1740 : // in case there is no dedicated spline, use the
1741 : // * pion spline for muons and the
1742 : // * proton spline for light nuclei
1743 0 : if (!responseFunction) {
1744 0 : if (ispecie==Int_t(AliPID::kMuon)) {
1745 0 : responseFunction=static_cast<TSpline3*>(pionSpline);
1746 0 : } else if (ispecie>Int_t(AliPID::kProton)) {
1747 0 : responseFunction=static_cast<TSpline3*>(protonSpline);
1748 0 : }
1749 : }
1750 :
1751 0 : if (!responseFunction) {
1752 0 : AliError("***** Risk for unreliable TPC PID detected: ********");
1753 0 : AliError(Form(" No spline found for %s, this should not happen, please report", AliPID::ParticleName(ispecie) ));
1754 0 : continue;
1755 : }
1756 0 : SetUseDatabase(kTRUE);
1757 0 : SetResponseFunction((AliPID::EParticleType)ispecie, responseFunction);
1758 0 : AliInfo(Form("Adding spline: %d - %s (MD5(spline) = %s)",ispecie,responseFunction->GetName(),
1759 : GetChecksum(responseFunction).Data()));
1760 :
1761 0 : }
1762 :
1763 0 : return kTRUE;
1764 0 : }
1765 :
1766 : //______________________________________________________________________________
1767 : Bool_t AliTPCPIDResponse::SetMultiplicityCorrectionFromString(const TString& multiplicityData)
1768 : {
1769 : //
1770 : //
1771 : //
1772 :
1773 : const TObjArray *arrParameters=0x0;
1774 0 : if (!(arrParameters=GetMultiplicityCorrectionArrayFromString(multiplicityData))) {
1775 0 : return kFALSE;
1776 : }
1777 :
1778 0 : TString log("Setting multiplicity correction parameters for mean; tan theta; sigma: ");
1779 0 : log.Append(multiplicityData);
1780 0 : AliInfo(log.Data());
1781 :
1782 0 : const TObjArray *arrPar1 = static_cast<TObjArray*>(arrParameters->At(0));
1783 0 : for (Int_t ipar=0; ipar<arrPar1->GetEntriesFast(); ++ipar) {
1784 0 : SetParameterMultiplicityCorrection(ipar, static_cast<TObjString*>(arrPar1->At(ipar))->String().Atof());
1785 : }
1786 :
1787 0 : const TObjArray *arrPar2 = static_cast<TObjArray*>(arrParameters->At(1));
1788 0 : for (Int_t ipar=0; ipar<arrPar2->GetEntriesFast(); ++ipar) {
1789 0 : SetParameterMultiplicityCorrectionTanTheta(ipar, static_cast<TObjString*>(arrPar2->At(ipar))->String().Atof());
1790 : }
1791 :
1792 0 : const TObjArray *arrPar3 = static_cast<TObjArray*>(arrParameters->At(2));
1793 0 : for (Int_t ipar=0; ipar<arrPar3->GetEntriesFast(); ++ipar) {
1794 0 : SetParameterMultiplicitySigmaCorrection(ipar, static_cast<TObjString*>(arrPar3->At(ipar))->String().Atof());
1795 : }
1796 :
1797 0 : delete arrParameters;
1798 : return kTRUE;
1799 0 : }
1800 :
1801 : //______________________________________________________________________________
1802 : TObjArray* AliTPCPIDResponse::GetMultiplicityCorrectionArrayFromString(const TString& corrections)
1803 : {
1804 : // the corrections string is supposed to be in the format
1805 : // parMC_0,parMC_1,parMC_2,parMC_3, parMC_4; parMCTT_0,parMCTT_1,parMCTT_2; parMSC_0,parMSC_1,parMSC_2, parMSC_3
1806 : // where parMC are the parameters for AliTPCPIDResponse::SetParameterMultiplicityCorrection
1807 : // parMCTT are the parameters for AliTPCPIDResponse::SetParameterMultiplicityCorrectionTanTheta
1808 : // parMSC are the parameters for AliTPCPIDResponse::SetParameterMultiplicitySigmaCorrection
1809 :
1810 : const Int_t nSets=3;
1811 : const Int_t nPars[nSets]={5,3,4};
1812 :
1813 0 : AliTPCPIDResponse temp;
1814 0 : TObjArray *arrCorrectionSets = corrections.Tokenize(";");
1815 0 : if (arrCorrectionSets->GetEntriesFast() != nSets){
1816 0 : temp.Error("AliTPCPIDResponse::CheckMultiplicityCorrectionString","Number of parameter sets not equal to %d. Please read documentation of this function", nSets);
1817 0 : delete arrCorrectionSets;
1818 0 : return 0x0;
1819 : }
1820 :
1821 0 : for (Int_t iset=0 ;iset<nSets; ++iset) {
1822 0 : TObjString *string=(TObjString*)arrCorrectionSets->RemoveAt(iset);
1823 0 : TObjArray *arrTmp=string->String().Tokenize(",");
1824 0 : delete string;
1825 : string=0x0;
1826 :
1827 0 : if (arrTmp->GetEntriesFast() != nPars[iset]){
1828 0 : temp.Error("AliTPCPIDResponse::CheckMultiplicityCorrectionString","Number of parameters in set %d not equal to %d. Please read documentation of this function", iset, nPars[iset]);
1829 0 : delete arrCorrectionSets;
1830 0 : delete arrTmp;
1831 0 : return 0x0;
1832 : }
1833 :
1834 0 : arrCorrectionSets->AddAt(arrTmp, iset);
1835 0 : }
1836 :
1837 0 : return arrCorrectionSets;
1838 0 : }
1839 :
1840 : //______________________________________________________________________________
1841 : Bool_t AliTPCPIDResponse::SetdEdxTypeFromString(const TString& dEdxTypeSet)
1842 : {
1843 : // Set up the dEdx usage parsing dEdxTypeSet
1844 : // 6 comma separated values are expected.
1845 : //
1846 : // The format assumes is:
1847 : // 0: dEdxType (AliTPCPIDResponse::ETPCdEdxType)
1848 : // 1: dEdxChargeType (qTot, qMax)
1849 : // 2: dEdxWeightType (0-measured pid clusters, 1-measured pid clusters + subThreshold clusters)
1850 : // 3: dEdxIROCweight (specific weight for IROC)
1851 : // 4: dEdxOROCmedWeight (specific weight for OROC medium pads)
1852 : // 5: dEdxOROClongWeight (specific weight for OROC long pads)
1853 0 : TObjArray *arrParams = dEdxTypeSet.Tokenize(",");
1854 :
1855 : Bool_t retVal=kFALSE;
1856 :
1857 0 : if (arrParams->GetEntriesFast() == 6) {
1858 0 : const Int_t dEdxType = ((TObjString*)arrParams->At(0))->String().Atoi();
1859 0 : const Int_t dEdxChargeType = ((TObjString*)arrParams->At(1))->String().Atoi();
1860 0 : const Int_t dEdxWeightType = ((TObjString*)arrParams->At(2))->String().Atoi();
1861 0 : const Double_t dEdxIROCweight = ((TObjString*)arrParams->At(3))->String().Atof();
1862 0 : const Double_t dEdxOROCmedWeight = ((TObjString*)arrParams->At(4))->String().Atof();
1863 0 : const Double_t dEdxOROClongWeight = ((TObjString*)arrParams->At(5))->String().Atof();
1864 0 : SetdEdxType((AliTPCPIDResponse::ETPCdEdxType)dEdxType, dEdxChargeType, dEdxWeightType, dEdxIROCweight, dEdxOROCmedWeight, dEdxOROClongWeight);
1865 0 : AliInfo(TString::Format("Setting custom TPC dEdxType: %d, %d, %d, %.2f, %.2f, %.2f",
1866 : dEdxType, dEdxChargeType, dEdxWeightType, dEdxIROCweight, dEdxOROCmedWeight, dEdxOROClongWeight));
1867 : retVal=kTRUE;
1868 0 : } else {
1869 0 : AliError("Wrong number of parameters for custom TPC dEdxType");
1870 : }
1871 0 : delete arrParams;
1872 :
1873 0 : return retVal;
1874 0 : }
1875 :
1876 : //______________________________________________________________________________
1877 : Bool_t AliTPCPIDResponse::SetdEdxResolutionFromString(const TString& dEdxResolutionString)
1878 : {
1879 : // Set up the dEdx resolution parsing dEdxResolutionString
1880 : // 2 comma separated values are expected to describe the relative resolution
1881 : // if resN2=0 no dependence on the number of clusters will be described
1882 : // sigma_rel = res0 * sqrt(1+resN2/npoint)
1883 : //
1884 : // The format assumes is:
1885 : // 0: res0
1886 : // 1: resN2
1887 0 : TObjArray *arrParams = dEdxResolutionString.Tokenize(",");
1888 :
1889 : Bool_t retVal=kFALSE;
1890 :
1891 0 : if (arrParams->GetEntriesFast() == 2) {
1892 0 : const Float_t res0 = ((TObjString*)arrParams->At(0))->String().Atoi();
1893 0 : const Float_t resN2 = ((TObjString*)arrParams->At(1))->String().Atoi();
1894 0 : AliInfo(TString::Format("Setting custom TPC dEdxResolution: %.2f, %.2f",
1895 : fRes0[0], fResN2[0]));
1896 0 : fRes0 [0]=res0;
1897 0 : fResN2[0]=resN2;
1898 : retVal=kTRUE;
1899 0 : } else {
1900 0 : AliError("Wrong number of parameters for custom TPC dEdxResolution");
1901 : }
1902 0 : delete arrParams;
1903 :
1904 0 : return retVal;
1905 0 : }
1906 :
1907 : //______________________________________________________________________________
1908 : TString AliTPCPIDResponse::GetChecksum(const TObject* obj)
1909 : {
1910 : // Return the checksum for an object obj (tested to work properly at least for histograms and TSplines).
1911 :
1912 0 : TString fileName = Form("tempChecksum.C"); // File name must be fixed for data type "TSpline3", since the file name will end up in the file content!
1913 :
1914 : // For parallel processing, a unique file pathname is required. Uniqueness can be guaranteed by using a unique directory name
1915 : UInt_t index = 0;
1916 0 : TString uniquePathName = Form("tempChecksum_%u", index);
1917 :
1918 : // To get a unique path name, increase the index until no directory
1919 : // of such a name exists.
1920 : // NOTE: gSystem->AccessPathName(...) returns kTRUE, if the access FAILED!
1921 0 : while (!gSystem->AccessPathName(uniquePathName.Data()))
1922 0 : uniquePathName = Form("tempChecksum_%u", ++index);
1923 :
1924 0 : AliTPCPIDResponse temp;
1925 0 : if (gSystem->mkdir(uniquePathName.Data()) < 0) {
1926 0 : temp.Error("AliTPCPIDResponse::GetChecksum","Could not create temporary directory to store temp file for checksum determination!");
1927 0 : return "ERROR";
1928 : }
1929 :
1930 0 : TString option = "";
1931 :
1932 : // Save object as a macro, which will be deleted immediately after the checksum has been computed
1933 : // (does not work for desired data types if saved as *.root for some reason) - one only wants to compare the content, not
1934 : // the modification time etc. ...
1935 0 : if (dynamic_cast<const TH1*>(obj)) {
1936 0 : option = "colz"; // Histos need this option, since w/o this option, a counter is added to the filename
1937 : }
1938 :
1939 :
1940 : // SaveAs must be called with the fixed fileName only, since the first argument goes into the file content
1941 : // for some object types. Thus, change the directory, save the file and then go back
1942 0 : TString oldDir = gSystem->pwd();
1943 0 : gSystem->cd(uniquePathName.Data());
1944 0 : obj->SaveAs(fileName.Data(), option.Data());
1945 0 : gSystem->cd(oldDir.Data());
1946 :
1947 : // Use the file to calculate the MD5 checksum
1948 0 : TMD5* md5 = TMD5::FileChecksum(Form("%s/%s", uniquePathName.Data(), fileName.Data()));
1949 0 : TString checksum = md5->AsString();
1950 :
1951 : // Clean up
1952 0 : delete md5;
1953 0 : gSystem->Exec(Form("rm -rf %s", uniquePathName.Data()));
1954 :
1955 0 : return checksum;
1956 0 : }
|
