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