add common abstract interface classes for flat and fat ESDs
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTPCPIDResponse.cxx
CommitLineData
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
44ClassImp(AliTPCPIDResponse);
d2aa6df0 45
912f7fdc 46
47AliTPCPIDResponse *AliTPCPIDResponse::fgInstance =0;
48
644666df 49const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]=
50{
51 "", //default - no name
52 "1", //all high
53 "2", //OROC only
54 "unknown"
55};
10d100d4 56
57//_________________________________________________________________________
58AliTPCPIDResponse::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 106AliTPCPIDResponse::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//_________________________________________________________________________
137AliTPCPIDResponse::~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//_________________________________________________________________________
162AliTPCPIDResponse::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//_________________________________________________________________________
218AliTPCPIDResponse& 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 279Double_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//_________________________________________________________________________
299void 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//_________________________________________________________________________
315void 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//_________________________________________________________________________
205d83d8 323Double_t AliTPCPIDResponse::GetExpectedSignal(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//_________________________________________________________________________
205d83d8 357Double_t AliTPCPIDResponse::GetExpectedSigma(Float_t mom,
358 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//_________________________________________________________________________
379void 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 390Double_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//_________________________________________________________________________
440Double_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//_________________________________________________________________________
479TSpline3* 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//_________________________________________________________________________
487TSpline3* 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//_________________________________________________________________________
505void 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//_________________________________________________________________________
514Int_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//_________________________________________________________________________
522void 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 531Double_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//_________________________________________________________________________
593Double_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//_________________________________________________________________________
617Float_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
644//_________________________________________________________________________
567624b5 645Float_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
673//_________________________________________________________________________
644666df 674Bool_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 776Double_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);
7294ba92 795 Int_t binX = fhEtaCorr->GetXaxis()->FindFixBin(tanTheta);
796 Int_t binY = fhEtaCorr->GetYaxis()->FindFixBin(1. / tpcSignal);
f84b18dd 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//_________________________________________________________________________
813Double_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//_________________________________________________________________________
841Double_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 885Double_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);
7294ba92 917 Int_t binX = fhEtaSigmaPar1->GetXaxis()->FindFixBin(tanTheta);
918 Int_t binY = fhEtaSigmaPar1->GetYaxis()->FindFixBin(1. / dEdxExpected);
f84b18dd 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//_________________________________________________________________________
935Double_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//_________________________________________________________________________
957Bool_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//_________________________________________________________________________
981Bool_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
644666df 1008//_________________________________________________________________________
87da0205 1009Double_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//_________________________________________________________________________
205d83d8 1074Double_t AliTPCPIDResponse::GetMultiplicityCorrectionFast(const AliVTrack *track, Double_t dEdxExpected, 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//_________________________________________________________________________
1098Double_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//_________________________________________________________________________
1127Double_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//_________________________________________________________________________
1182Double_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//_________________________________________________________________________
205d83d8 1230Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrectionFast(Double_t dEdxExpected, 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//_________________________________________________________________________
1251Double_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//_________________________________________________________________________
1280void 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
1304//_________________________________________________________________________
f85a3764 1305Bool_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//_____________________________________________________________________________
1336Int_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//_____________________________________________________________________________
1345void AliTPCPIDResponse::Print(Option_t* /*option*/) const
1346{
1347 //Print info
1348 fResponseFunctions.Print();
1349}
1350
1b45e564 1351
644666df 1352//_____________________________________________________________________________
1353AliTPCPIDResponse::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.
50ff8c36 1373 const AliVVtrack* ip = track->GetInnerParam();
1374 if (ip) track=dynamic_cast<const AliVTrack*>(ip);
f85a3764 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//_____________________________________________________________________________
1503Float_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//_____________________________________________________________________________
644666df 1528Bool_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
1562Double_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
1572Bool_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}