Here the comment:
[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"
41
42ClassImp(AliTPCPIDResponse);
d2aa6df0 43
644666df 44const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]=
45{
46 "", //default - no name
47 "1", //all high
48 "2", //OROC only
49 "unknown"
50};
10d100d4 51
52//_________________________________________________________________________
53AliTPCPIDResponse::AliTPCPIDResponse():
644666df 54 TNamed(),
d2aa6df0 55 fMIP(50.),
644666df 56 fRes0(),
57 fResN2(),
d2aa6df0 58 fKp1(0.0283086),
59 fKp2(2.63394e+01),
60 fKp3(5.04114e-11),
61 fKp4(2.12543),
62 fKp5(4.88663),
63 fUseDatabase(kFALSE),
644666df 64 fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
65 fVoltageMap(72),
66 fLowGainIROCthreshold(-40),
67 fBadIROCthreshhold(-70),
68 fLowGainOROCthreshold(-40),
69 fBadOROCthreshhold(-40),
70 fMaxBadLengthFraction(0.5),
f84b18dd 71 fMagField(0.),
72 fhEtaCorr(0x0),
73 fhEtaSigmaPar1(0x0),
74 fSigmaPar0(0.0)
10d100d4 75{
76 //
77 // The default constructor
78 //
644666df 79 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;}
10d100d4 80}
1b45e564 81/*TODO remove?
10d100d4 82//_________________________________________________________________________
d2aa6df0 83AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
644666df 84 TNamed(),
d2aa6df0 85 fMIP(param[0]),
644666df 86 fRes0(),
87 fResN2(),
d2aa6df0 88 fKp1(0.0283086),
89 fKp2(2.63394e+01),
90 fKp3(5.04114e-11),
91 fKp4(2.12543),
92 fKp5(4.88663),
93 fUseDatabase(kFALSE),
644666df 94 fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
95 fVoltageMap(72),
96 fLowGainIROCthreshold(-40),
97 fBadIROCthreshhold(-70),
98 fLowGainOROCthreshold(-40),
99 fBadOROCthreshhold(-40),
100 fMaxBadLengthFraction(0.5),
f84b18dd 101 fMagField(0.),
102 fhEtaCorr(0x0),
103 fhEtaSigmaPar1(0x0),
104 fSigmaPar0(0.0)
10d100d4 105{
106 //
107 // The main constructor
108 //
644666df 109 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];}
110}
1b45e564 111*/
f84b18dd 112
113//_________________________________________________________________________
114AliTPCPIDResponse::~AliTPCPIDResponse()
115{
116 //
117 // Destructor
118 //
119
120 delete fhEtaCorr;
121 fhEtaCorr = 0x0;
122
123 delete fhEtaSigmaPar1;
124 fhEtaSigmaPar1 = 0x0;
125}
126
127
644666df 128//_________________________________________________________________________
129AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that):
130 TNamed(that),
131 fMIP(that.fMIP),
132 fRes0(),
133 fResN2(),
134 fKp1(that.fKp1),
135 fKp2(that.fKp2),
136 fKp3(that.fKp3),
137 fKp4(that.fKp4),
138 fKp5(that.fKp5),
139 fUseDatabase(that.fUseDatabase),
140 fResponseFunctions(that.fResponseFunctions),
141 fVoltageMap(that.fVoltageMap),
142 fLowGainIROCthreshold(that.fLowGainIROCthreshold),
143 fBadIROCthreshhold(that.fBadIROCthreshhold),
144 fLowGainOROCthreshold(that.fLowGainOROCthreshold),
145 fBadOROCthreshhold(that.fBadOROCthreshhold),
146 fMaxBadLengthFraction(that.fMaxBadLengthFraction),
f84b18dd 147 fMagField(that.fMagField),
148 fhEtaCorr(0x0),
149 fhEtaSigmaPar1(0x0),
150 fSigmaPar0(that.fSigmaPar0)
644666df 151{
152 //copy ctor
153 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
f84b18dd 154
155 // Copy eta maps
1b45e564 156 if (that.fhEtaCorr) {
3800f96f 157 fhEtaCorr = new TH2D(*(that.fhEtaCorr));
158 fhEtaCorr->SetDirectory(0);
159 }
f84b18dd 160
1b45e564 161 if (that.fhEtaSigmaPar1) {
3800f96f 162 fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
163 fhEtaSigmaPar1->SetDirectory(0);
164 }
644666df 165}
166
167//_________________________________________________________________________
168AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
169{
170 //assignment
171 if (&that==this) return *this;
172 TNamed::operator=(that);
173 fMIP=that.fMIP;
174 fKp1=that.fKp1;
175 fKp2=that.fKp2;
176 fKp3=that.fKp3;
177 fKp4=that.fKp4;
178 fKp5=that.fKp5;
179 fUseDatabase=that.fUseDatabase;
180 fResponseFunctions=that.fResponseFunctions;
181 fVoltageMap=that.fVoltageMap;
182 fLowGainIROCthreshold=that.fLowGainIROCthreshold;
183 fBadIROCthreshhold=that.fBadIROCthreshhold;
184 fLowGainOROCthreshold=that.fLowGainOROCthreshold;
185 fBadOROCthreshhold=that.fBadOROCthreshhold;
186 fMaxBadLengthFraction=that.fMaxBadLengthFraction;
187 fMagField=that.fMagField;
188 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
f84b18dd 189
190 delete fhEtaCorr;
f85a3764 191 fhEtaCorr=0x0;
1b45e564 192 if (that.fhEtaCorr) {
3800f96f 193 fhEtaCorr = new TH2D(*(that.fhEtaCorr));
194 fhEtaCorr->SetDirectory(0);
195 }
f84b18dd 196
197 delete fhEtaSigmaPar1;
f85a3764 198 fhEtaSigmaPar1=0x0;
1b45e564 199 if (that.fhEtaSigmaPar1) {
3800f96f 200 fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
201 fhEtaSigmaPar1->SetDirectory(0);
202 }
f84b18dd 203
204 fSigmaPar0 = that.fSigmaPar0;
205
644666df 206 return *this;
10d100d4 207}
208
d2aa6df0 209//_________________________________________________________________________
10d100d4 210Double_t AliTPCPIDResponse::Bethe(Double_t betaGamma) const {
211 //
212 // This is the Bethe-Bloch function normalised to 1 at the minimum
213 // WARNING
214 // Simulated and reconstructed Bethe-Bloch differs
215 // Simulated curve is the dNprim/dx
216 // Reconstructed is proportianal dNtot/dx
217 // Temporary fix for production - Simple linear correction function
218 // Future 2 Bethe Bloch formulas needed
219 // 1. for simulation
220 // 2. for reconstructed PID
221 //
d2aa6df0 222
223// const Float_t kmeanCorrection =0.1;
10d100d4 224 Double_t bb=
225 AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5);
226 return bb*fMIP;
227}
228
229//_________________________________________________________________________
230void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1,
231 Double_t kp2,
232 Double_t kp3,
233 Double_t kp4,
234 Double_t kp5) {
235 //
236 // Set the parameters of the ALEPH Bethe-Bloch formula
237 //
238 fKp1=kp1;
239 fKp2=kp2;
240 fKp3=kp3;
241 fKp4=kp4;
242 fKp5=kp5;
243}
644666df 244
10d100d4 245//_________________________________________________________________________
246void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) {
247 //
248 // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint)
249 //
644666df 250 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;}
10d100d4 251}
252
253//_________________________________________________________________________
254Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom,
d2aa6df0 255 AliPID::EParticleType n) const {
10d100d4 256 //
f84b18dd 257 // Deprecated function (for backward compatibility). Please use
258 // GetExpectedSignal(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource,
259 // Bool_t correctEta = kTRUE);
260 // instead!
261 //
262 //
10d100d4 263 // Calculates the expected PID signal as the function of
264 // the information stored in the track, for the specified particle type
265 //
266 // At the moment, these signals are just the results of calling the
267 // Bethe-Bloch formula.
268 // This can be improved. By taking into account the number of
269 // assigned clusters and/or the track dip angle, for example.
270 //
271
f84b18dd 272 //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
273 // !!! Splines for light nuclei need to be normalised to this factor !!!
274 const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3);
275
539a5a59 276 Double_t mass=AliPID::ParticleMassZ(n);
f84b18dd 277 if (!fUseDatabase) return Bethe(mom/mass) * chargeFactor;
d2aa6df0 278 //
644666df 279 const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
00a38d07 280
f84b18dd 281 if (!responseFunction) return Bethe(mom/mass) * chargeFactor;
282
00a38d07 283 return fMIP*responseFunction->Eval(mom/mass)*chargeFactor;
d2aa6df0 284
10d100d4 285}
286
287//_________________________________________________________________________
288Double_t AliTPCPIDResponse::GetExpectedSigma(const Float_t mom,
f84b18dd 289 const Int_t nPoints,
644666df 290 AliPID::EParticleType n) const {
10d100d4 291 //
f84b18dd 292 // Deprecated function (for backward compatibility). Please use
293 // GetExpectedSigma(onst AliVTrack* track, AliPID::EParticleType species,
294 // ETPCdEdxSource dedxSource, Bool_t correctEta) instead!
295 //
296 //
10d100d4 297 // Calculates the expected sigma of the PID signal as the function of
298 // the information stored in the track, for the specified particle type
299 //
644666df 300
301 if (nPoints != 0)
302 return GetExpectedSignal(mom,n)*fRes0[0]*sqrt(1. + fResN2[0]/nPoints);
303 else
304 return GetExpectedSignal(mom,n)*fRes0[0];
305}
306
307////////////////////////////////////////////////////NEW//////////////////////////////
308
309//_________________________________________________________________________
310void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario) {
311 //
312 // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint)
313 //
314 if ((Int_t)gainScenario>(Int_t)fgkNumberOfGainScenarios) return; //TODO: better handling!
315 fRes0[gainScenario]=res0;
316 fResN2[gainScenario]=resN2;
317}
318
f84b18dd 319
644666df 320//_________________________________________________________________________
f84b18dd 321Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
644666df 322 AliPID::EParticleType species,
f84b18dd 323 Double_t /*dEdx*/,
324 const TSpline3* responseFunction,
325 Bool_t correctEta) const
644666df 326{
327 // Calculates the expected PID signal as the function of
f84b18dd 328 // the information stored in the track and the given parameters,
329 // for the specified particle type
644666df 330 //
331 // At the moment, these signals are just the results of calling the
f84b18dd 332 // Bethe-Bloch formula plus, if desired, taking into account the eta dependence.
644666df 333 // This can be improved. By taking into account the number of
334 // assigned clusters and/or the track dip angle, for example.
10d100d4 335 //
336
f84b18dd 337 Double_t mom=track->GetTPCmomentum();
644666df 338 Double_t mass=AliPID::ParticleMass(species);
339
f84b18dd 340 //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
341 // !!! Splines for light nuclei need to be normalised to this factor !!!
342 const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3);
343
344 if (!responseFunction)
345 return Bethe(mom/mass) * chargeFactor;
644666df 346
f84b18dd 347 Double_t dEdxSplines = fMIP*responseFunction->Eval(mom/mass) * chargeFactor;
348
349 if (!correctEta)
350 return dEdxSplines;
351
352 //TODO Alternatively take current track dEdx
353 //return dEdxSplines * GetEtaCorrection(track, dEdx);
354 return dEdxSplines * GetEtaCorrection(track, dEdxSplines);
644666df 355}
356
f84b18dd 357
644666df 358//_________________________________________________________________________
359Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
360 AliPID::EParticleType species,
f84b18dd 361 ETPCdEdxSource dedxSource,
362 Bool_t correctEta) const
644666df 363{
364 // Calculates the expected PID signal as the function of
365 // the information stored in the track, for the specified particle type
366 //
367 // At the moment, these signals are just the results of calling the
f84b18dd 368 // Bethe-Bloch formula plus, if desired, taking into account the eta dependence.
644666df 369 // This can be improved. By taking into account the number of
370 // assigned clusters and/or the track dip angle, for example.
371 //
372
f84b18dd 373 if (!fUseDatabase) {
374 //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
375 // !!! Splines for light nuclei need to be normalised to this factor !!!
376 const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3);
644666df 377
f84b18dd 378 return Bethe(track->GetTPCmomentum() / AliPID::ParticleMass(species)) * chargeFactor;
379 }
644666df 380
f84b18dd 381 Double_t dEdx = -1;
382 Int_t nPoints = -1;
383 ETPCgainScenario gainScenario = kGainScenarioInvalid;
384 TSpline3* responseFunction = 0x0;
385
386 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) {
387 // Something is wrong with the track -> Return obviously invalid value
388 return -999;
389 }
390
391 // Charge factor already taken into account inside the following function call
392 return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta);
644666df 393}
394
395//_________________________________________________________________________
396TSpline3* AliTPCPIDResponse::GetResponseFunction( AliPID::EParticleType type,
397 AliTPCPIDResponse::ETPCgainScenario gainScenario ) const
398{
399 //get response function
400 return dynamic_cast<TSpline3*>(fResponseFunctions.At(ResponseFunctionIndex(type,gainScenario)));
401}
402
403//_________________________________________________________________________
404TSpline3* AliTPCPIDResponse::GetResponseFunction( const AliVTrack* track,
405 AliPID::EParticleType species,
f84b18dd 406 ETPCdEdxSource dedxSource) const
644666df 407{
408 //the splines are stored in an array, different scenarios
409
f84b18dd 410 Double_t dEdx = -1;
411 Int_t nPoints = -1;
412 ETPCgainScenario gainScenario = kGainScenarioInvalid;
413 TSpline3* responseFunction = 0x0;
414
415 if (ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
416 return responseFunction;
644666df 417
418 return NULL;
419}
420
421//_________________________________________________________________________
422void AliTPCPIDResponse::ResetSplines()
423{
424 //reset the array with splines
425 for (Int_t i=0;i<fResponseFunctions.GetEntriesFast();i++)
426 {
427 fResponseFunctions.AddAt(NULL,i);
428 }
429}
430//_________________________________________________________________________
431Int_t AliTPCPIDResponse::ResponseFunctionIndex( AliPID::EParticleType species,
432 ETPCgainScenario gainScenario ) const
433{
434 //get the index in fResponseFunctions given type and scenario
435 return Int_t(species)+Int_t(gainScenario)*fgkNumberOfParticleSpecies;
436}
437
438//_________________________________________________________________________
439void AliTPCPIDResponse::SetResponseFunction( TObject* o,
440 AliPID::EParticleType species,
441 ETPCgainScenario gainScenario )
442{
443 fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario));
444}
445
f84b18dd 446
644666df 447//_________________________________________________________________________
f84b18dd 448Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
449 AliPID::EParticleType species,
450 ETPCgainScenario gainScenario,
451 Double_t dEdx,
452 Int_t nPoints,
453 const TSpline3* responseFunction,
454 Bool_t correctEta) const
644666df 455{
456 // Calculates the expected sigma of the PID signal as the function of
f84b18dd 457 // the information stored in the track and the given parameters,
458 // for the specified particle type
644666df 459 //
460
f84b18dd 461 if (!responseFunction)
462 return 999;
463
464 // If no sigma map is available or if no eta correction is requested (sigma maps only for corrected eta!), use the old parametrisation
465 if (!fhEtaSigmaPar1 || !correctEta) {
466 if (nPoints != 0)
467 return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE) *
468 fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints);
469 else
470 return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE)*fRes0[gainScenario];
471 }
472
473 if (nPoints > 0) {
474 Double_t sigmaPar1 = GetSigmaPar1(track, species, dEdx, responseFunction);
475
476 return GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE)*sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints);
477 }
478 else {
479 // One should never have/take tracks with 0 dEdx clusters!
480 return 999;
481 }
644666df 482}
483
f84b18dd 484
644666df 485//_________________________________________________________________________
486Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
487 AliPID::EParticleType species,
f84b18dd 488 ETPCdEdxSource dedxSource,
489 Bool_t correctEta) const
644666df 490{
491 // Calculates the expected sigma of the PID signal as the function of
492 // the information stored in the track, for the specified particle type
493 // and dedx scenario
494 //
495
f84b18dd 496 Double_t dEdx = -1;
497 Int_t nPoints = -1;
498 ETPCgainScenario gainScenario = kGainScenarioInvalid;
499 TSpline3* responseFunction = 0x0;
500
501 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
502 return 999; //TODO: better handling!
503
504 return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta);
644666df 505}
506
f84b18dd 507
644666df 508//_________________________________________________________________________
509Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
510 AliPID::EParticleType species,
f84b18dd 511 ETPCdEdxSource dedxSource,
512 Bool_t correctEta) const
644666df 513{
f84b18dd 514 //Calculates the number of sigmas of the PID signal from the expected value
515 //for a given particle species in the presence of multiple gain scenarios
644666df 516 //inside the TPC
f84b18dd 517
518 Double_t dEdx = -1;
519 Int_t nPoints = -1;
520 ETPCgainScenario gainScenario = kGainScenarioInvalid;
521 TSpline3* responseFunction = 0x0;
522
523 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
524 return -999; //TODO: Better handling!
525
526 Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta);
527 Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta);
528 // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters
529 if (sigma >= 998)
530 return -999;
531 else
532 return (dEdx-bethe)/sigma;
644666df 533}
534
535//_________________________________________________________________________
536Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
537 AliPID::EParticleType species,
f84b18dd 538 ETPCdEdxSource dedxSource,
1b45e564 539 Double_t& dEdx,
540 Int_t& nPoints,
541 ETPCgainScenario& gainScenario,
f85a3764 542 TSpline3** responseFunction) const
644666df 543{
544 // Calculates the right parameters for PID
545 // dEdx parametrization for the proper gain scenario, dEdx
546 // and NPoints used for dEdx
547 // based on the track geometry (which chambers it crosses) for the specified particle type
548 // and preferred source of dedx.
549 // returns true on success
550
f84b18dd 551
552 if (dedxSource == kdEdxDefault) {
553 // Fast handling for default case. In addition: Keep it simple (don't call additional functions) to
554 // avoid possible bugs
f85a3764 555
556 // GetTPCsignalTunedOnData will be non-positive, if it has not been set (i.e. in case of MC NOT tuned to data).
557 // If this is the case, just take the normal signal
558 dEdx = track->GetTPCsignalTunedOnData();
559 if (dEdx <= 0) {
560 dEdx = track->GetTPCsignal();
561 }
562
f84b18dd 563 nPoints = track->GetTPCsignalN();
564 gainScenario = kDefault;
565
566 TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario));
567 *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
568
569 return kTRUE;
570 }
571
f85a3764 572 //TODO Proper handle of tuneMConData for other dEdx sources
f84b18dd 573
644666df 574 Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used)
575 Char_t ncl[3]; //same
576 Char_t nrows[3]; //same
f84b18dd 577 const AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo();
644666df 578
579 if (!dEdxInfo && dedxSource!=kdEdxDefault) //in one case its ok if we dont have the info
580 {
581 AliError("AliTPCdEdxInfo not available");
644666df 582 return kFALSE;
583 }
584
585 if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows);
586
587 //check if we cross a bad OROC in which case we reject
f85a3764 588 EChamberStatus trackOROCStatus = TrackStatus(track,2);
589 if (trackOROCStatus==kChamberOff || trackOROCStatus==kChamberLowGain)
644666df 590 {
644666df 591 return kFALSE;
592 }
593
594 switch (dedxSource)
595 {
644666df 596 case kdEdxOROC:
597 {
f85a3764 598 if (trackOROCStatus==kChamberInvalid) return kFALSE; //never reached OROC
f84b18dd 599 dEdx = signal[3];
600 nPoints = ncl[2]+ncl[1];
601 gainScenario = kOROChigh;
644666df 602 break;
603 }
604 case kdEdxHybrid:
605 {
606 //if we cross a bad IROC we use OROC dedx, if we dont we use combined
607 EChamberStatus status = TrackStatus(track,1);
608 if (status!=kChamberHighGain)
609 {
f84b18dd 610 dEdx = signal[3];
611 nPoints = ncl[2]+ncl[1];
612 gainScenario = kOROChigh;
644666df 613 }
614 else
615 {
f84b18dd 616 dEdx = track->GetTPCsignal();
617 nPoints = track->GetTPCsignalN();
618 gainScenario = kALLhigh;
644666df 619 }
620 break;
621 }
622 default:
623 {
f84b18dd 624 dEdx = 0.;
625 nPoints = 0;
626 gainScenario = kGainScenarioInvalid;
644666df 627 return kFALSE;
628 }
629 }
f84b18dd 630 TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario));
631 *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
632
644666df 633 return kTRUE;
634}
635
f84b18dd 636
637//_________________________________________________________________________
638Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, Double_t dEdxSplines) const
639{
640 //
641 // Get eta correction for the given parameters.
642 //
643
644 if (!fhEtaCorr)
645 return 1.;
646
647 Double_t tpcSignal = dEdxSplines;
648
649 if (tpcSignal < 1.)
650 return 1.;
651
652 // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()).
653 // However, this value is not available for AODs and, thus, not for AliVTrack.
654 // Fortunately, the following formula allows to approximate the local tanTheta with the
655 // global theta angle -> This is for by far most of the tracks the same, but gives at
656 // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok.
657 Double_t tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
658 Int_t binX = fhEtaCorr->GetXaxis()->FindBin(tanTheta);
659 Int_t binY = fhEtaCorr->GetYaxis()->FindBin(1. / tpcSignal);
660
661 if (binX == 0)
662 binX = 1;
663 if (binX > fhEtaCorr->GetXaxis()->GetNbins())
664 binX = fhEtaCorr->GetXaxis()->GetNbins();
665
666 if (binY == 0)
667 binY = 1;
668 if (binY > fhEtaCorr->GetYaxis()->GetNbins())
669 binY = fhEtaCorr->GetYaxis()->GetNbins();
670
671 return fhEtaCorr->GetBinContent(binX, binY);
672}
673
674
675//_________________________________________________________________________
676Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
677{
678 //
679 // Get eta correction for the given track.
680 //
681
682 if (!fhEtaCorr)
683 return 1.;
684
685 Double_t dEdx = -1;
686 Int_t nPoints = -1;
687 ETPCgainScenario gainScenario = kGainScenarioInvalid;
688 TSpline3* responseFunction = 0x0;
689
690 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
691 return 1.;
692
693 Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE);
694
695 //TODO Alternatively take current track dEdx
696 //return GetEtaCorrection(track, dEdx);
697
698 return GetEtaCorrection(track, dEdxSplines);
699}
700
701
702//_________________________________________________________________________
703Double_t AliTPCPIDResponse::GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
704{
705 //
706 // Get eta corrected dEdx for the given track. For the correction, the expected dEdx of
707 // the specified species will be used. If the species is set to AliPID::kUnknown, the
708 // dEdx of the track is used instead.
709 // WARNING: In the latter case, the eta correction might not be as good as if the
710 // expected dEdx is used, which is the way the correction factor is designed
711 // for.
712 // In any case, one has to decide carefully to which expected signal one wants to
713 // compare the corrected value - to the corrected or uncorrected.
714 // Anyhow, a safer way of looking e.g. at the n-sigma is to call
715 // the corresponding function GetNumberOfSigmas!
716 //
717
718 Double_t dEdx = -1;
719 Int_t nPoints = -1;
720 ETPCgainScenario gainScenario = kGainScenarioInvalid;
721 TSpline3* responseFunction = 0x0;
722
723 // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case
724 // it is not used anyway, so this causes no trouble.
725 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
726 return -1.;
727
728 Double_t etaCorr = 0.;
729
730 if (species < AliPID::kUnknown) {
731 Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE);
732 etaCorr = GetEtaCorrection(track, dEdxSplines);
733 }
734 else {
735 etaCorr = GetEtaCorrection(track, dEdx);
736 }
737
738 if (etaCorr <= 0)
739 return -1.;
740
741 return dEdx / etaCorr;
742}
743
744
745
746//_________________________________________________________________________
747Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, Double_t dEdx, const TSpline3* responseFunction) const
748{
749 //
750 // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track.
751 //
752
753 if (!fhEtaSigmaPar1)
754 return 999;
755
756 // The sigma maps are created with data sets that are already eta corrected and for which the
757 // splines have been re-created. Therefore, the value for the lookup needs to be the value of
758 // the splines without any additional eta correction.
759 // NOTE: This is due to the method the maps are created. The track dEdx (not the expected one!)
760 // is corrected to uniquely related a momemtum bin with an expected dEdx, where the expected dEdx
761 // equals the track dEdx for all eta (thanks to the correction and the re-fit of the splines).
762 // Consequently, looking up the uncorrected expected dEdx at a given tanTheta yields the correct
763 // sigma parameter!
764 // Also: It has to be the spline dEdx, since one wants to get the sigma for the assumption(!)
765 // of such a particle, which by assumption then has this dEdx value
766
767 Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE);
768
769 if (dEdxExpected < 1.)
770 return 999;
771
772 // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()).
773 // However, this value is not available for AODs and, thus, not or AliVTrack.
774 // Fortunately, the following formula allows to approximate the local tanTheta with the
775 // global theta angle -> This is for by far most of the tracks the same, but gives at
776 // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok.
777 Double_t tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
778 Int_t binX = fhEtaSigmaPar1->GetXaxis()->FindBin(tanTheta);
779 Int_t binY = fhEtaSigmaPar1->GetYaxis()->FindBin(1. / dEdxExpected);
780
781 if (binX == 0)
782 binX = 1;
783 if (binX > fhEtaSigmaPar1->GetXaxis()->GetNbins())
784 binX = fhEtaSigmaPar1->GetXaxis()->GetNbins();
785
786 if (binY == 0)
787 binY = 1;
788 if (binY > fhEtaSigmaPar1->GetYaxis()->GetNbins())
789 binY = fhEtaSigmaPar1->GetYaxis()->GetNbins();
790
791 return fhEtaSigmaPar1->GetBinContent(binX, binY);
792}
793
794
795//_________________________________________________________________________
796Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
797{
798 //
799 // Get eta correction for the given track.
800 //
801
802 if (!fhEtaSigmaPar1)
803 return 999;
804
805 Double_t dEdx = -1;
806 Int_t nPoints = -1;
807 ETPCgainScenario gainScenario = kGainScenarioInvalid;
808 TSpline3* responseFunction = 0x0;
809
810 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
811 return 999;
812
813 return GetSigmaPar1(track, species, dEdx, responseFunction);
814}
815
816
817//_________________________________________________________________________
818Bool_t AliTPCPIDResponse::SetEtaCorrMap(TH2D* hMap)
819{
820 //
821 // Load map for TPC eta correction (a copy is stored and will be deleted automatically).
822 // If hMap is 0x0,the eta correction will be disabled and kFALSE is returned.
823 // If the map can be set, kTRUE is returned.
824 //
825
826 delete fhEtaCorr;
827
828 if (!hMap) {
829 fhEtaCorr = 0x0;
830
831 return kFALSE;
832 }
833
834 fhEtaCorr = (TH2D*)(hMap->Clone());
835 fhEtaCorr->SetDirectory(0);
836
837 return kTRUE;
838}
839
840//_________________________________________________________________________
841Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0)
842{
843 //
844 // Load map for TPC sigma map (a copy is stored and will be deleted automatically):
845 // Parameter 1 is stored as a 2D map (1/dEdx vs. tanTheta_local) and parameter 0 is
846 // a constant. If hSigmaPar1Map is 0x0, the old sigma parametrisation will be used
847 // (and sigmaPar0 is ignored!) and kFALSE is returned.
848 // If the map can be set, sigmaPar0 is also set and kTRUE will be returned.
849 //
850
851 delete fhEtaSigmaPar1;
852
853 if (!hSigmaPar1Map) {
854 fhEtaSigmaPar1 = 0x0;
855 fSigmaPar0 = 0.0;
856
857 return kFALSE;
858 }
859
860 fhEtaSigmaPar1 = (TH2D*)(hSigmaPar1Map->Clone());
861 fhEtaSigmaPar1->SetDirectory(0);
862 fSigmaPar0 = sigmaPar0;
863
864 return kTRUE;
865}
866
867
644666df 868//_________________________________________________________________________
f85a3764 869Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner,
870 Double_t* trackPositionOuter,
871 Float_t& inphi,
644666df 872 Float_t& outphi,
f85a3764 873 Int_t& in,
1b45e564 874 Int_t& out ) const
644666df 875{
876 //calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses
877 //for OROC chamber numbers add 36
878 //returned angles are between (0,2pi)
644666df 879
880 inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]);
881 outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]);
882
883 if (inphi<0) {inphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
884 if (outphi<0) {outphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
885
886 in = sectorNumber(inphi);
887 out = sectorNumber(outphi);
888
889 //for the C side (positive z) offset by 18
890 if (trackPositionInner[2]>0.0)
891 {
892 in+=18;
893 out+=18;
894 }
895 return kTRUE;
896}
1b45e564 897
898
644666df 899//_____________________________________________________________________________
900Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const
901{
902 //calculate sector number
903 const Float_t width=TMath::TwoPi()/18.;
904 return TMath::Floor(phi/width);
905}
906
1b45e564 907
644666df 908//_____________________________________________________________________________
909void AliTPCPIDResponse::Print(Option_t* /*option*/) const
910{
911 //Print info
912 fResponseFunctions.Print();
913}
914
1b45e564 915
644666df 916//_____________________________________________________________________________
917AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const
918{
919 //status of the track: if it crosses any bad chambers return kChamberOff;
920 //IROC:layer=1, OROC:layer=2
921 if (layer<1 || layer>2) layer=1;
922 Int_t in=0;
923 Int_t out=0;
924 Float_t inphi=0.;
925 Float_t outphi=0.;
926 Float_t innerRadius = (layer==1)?83.0:133.7;
927 Float_t outerRadius = (layer==1)?133.5:247.7;
f85a3764 928
929 /////////////////////////////////////////////////////////////////////////////
930 //find out where track enters and leaves the layer.
931 //
1b45e564 932 Double_t trackPositionInner[3];
933 Double_t trackPositionOuter[3];
934
f85a3764 935 //if there is no inner param this could mean we're using the AOD track,
936 //we still can extrapolate from the vertex - so use those params.
937 const AliExternalTrackParam* ip = track->GetInnerParam();
938 if (ip) track=ip;
939
940 Bool_t trackAtInner = track->GetXYZAt(innerRadius, fMagField, trackPositionInner);
941 Bool_t trackAtOuter = track->GetXYZAt(outerRadius, fMagField, trackPositionOuter);
942
943 if (!trackAtInner)
944 {
945 //if we dont even enter inner radius we do nothing and return invalid
946 inphi=0.0;
947 outphi=0.0;
948 in=0;
949 out=0;
950 return kChamberInvalid;
951 }
952
953 if (!trackAtOuter)
954 {
955 //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position
956 Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter);
957 Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]);
958 if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius)
959 {
960 //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]);
961 }
962 else
963 {
964 inphi=0.0;
965 outphi=0.0;
966 in=0;
967 out=0;
968 return kChamberInvalid;
969 }
970 }
971
972
1b45e564 973 if (!sectorNumbersInOut(trackPositionInner,
974 trackPositionOuter,
975 inphi,
976 outphi,
977 in,
f85a3764 978 out)) return kChamberInvalid;
644666df 979
f85a3764 980 /////////////////////////////////////////////////////////////////////////////
1b45e564 981 //now we have the location of the track we can check
f85a3764 982 //if it is in a good/bad chamber
983 //
644666df 984 Bool_t sideA = kTRUE;
1b45e564 985
644666df 986 if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE;
987
988 in=in%18;
989 out=out%18;
990
991 if (TMath::Abs(in-out)>9)
992 {
993 if (TMath::Max(in,out)==out)
994 {
995 Int_t tmp=out;
996 out = in;
997 in = tmp;
998 Float_t tmpphi=outphi;
999 outphi=inphi;
1000 inphi=tmpphi;
1001 }
1002 in-=18;
1003 inphi-=TMath::TwoPi();
1004 }
10d100d4 1005 else
644666df 1006 {
1007 if (TMath::Max(in,out)==in)
1008 {
1009 Int_t tmp=out;
1010 out=in;
1011 in=tmp;
1012 Float_t tmpphi=outphi;
1013 outphi=inphi;
1014 inphi=tmpphi;
1015 }
1016 }
1017
1018 Float_t trackLengthInBad = 0.;
1019 Float_t trackLengthInLowGain = 0.;
1020 Float_t trackLengthTotal = TMath::Abs(outphi-inphi);
f85a3764 1021 Float_t lengthFractionInBadSectors = 0.;
644666df 1022
1023 const Float_t sectorWidth = TMath::TwoPi()/18.;
1b45e564 1024
644666df 1025 for (Int_t i=in; i<=out; i++)
1026 {
1027 int j=i;
1028 if (i<0) j+=18; //correct for the negative values
1029 if (!sideA) j+=18; //move to the correct side
1b45e564 1030
644666df 1031 Float_t deltaPhi = 0.;
1032 Float_t phiEdge=sectorWidth*i;
1033 if (inphi>phiEdge) {deltaPhi=phiEdge+sectorWidth-inphi;}
1034 else if ((outphi>=phiEdge) && (outphi<(phiEdge+sectorWidth))) {deltaPhi=outphi-phiEdge;}
1035 else {deltaPhi=sectorWidth;}
1b45e564 1036
644666df 1037 Float_t v = fVoltageMap[(layer==1)?(j):(j+36)];
1b45e564 1038 if (v<=fBadIROCthreshhold)
1039 {
1040 trackLengthInBad+=deltaPhi;
f85a3764 1041 lengthFractionInBadSectors=1.;
1042 }
1b45e564 1043 if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold)
f85a3764 1044 {
1045 trackLengthInLowGain+=deltaPhi;
1046 lengthFractionInBadSectors=1.;
1047 }
644666df 1048 }
1049
1050 //for now low gain and bad (off) chambers are treated equally
f85a3764 1051 if (trackLengthTotal>0)
1052 lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal;
644666df 1053
1054 //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 1055
1056 if (lengthFractionInBadSectors>fMaxBadLengthFraction)
644666df 1057 {
1058 //printf("%%%%%%%% %s kChamberLowGain\n",(layer==1)?"IROC":"OROC");
1059 return kChamberLowGain;
1060 }
1b45e564 1061
644666df 1062 return kChamberHighGain;
1063}
1064
1b45e564 1065
644666df 1066//_____________________________________________________________________________
1067Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const
1068{
1069 //return the radius of the outermost padrow containing a cluster in TPC
1070 //for the track
1071 const TBits* clusterMap=track->GetTPCClusterMapPtr();
1072 if (!clusterMap) return 0.;
1073
1074 //from AliTPCParam, radius of first IROC row
1b45e564 1075 const Float_t rfirstIROC = 8.52249984741210938e+01;
644666df 1076 const Float_t padrowHeightIROC = 0.75;
1077 const Float_t rfirstOROC0 = 1.35100006103515625e+02;
1078 const Float_t padrowHeightOROC0 = 1.0;
1079 const Float_t rfirstOROC1 = 1.99350006103515625e+02;
1080 const Float_t padrowHeightOROC1 = 1.5;
1081
1082 Int_t maxPadRow=160;
1083 while ((--maxPadRow)>0 && !clusterMap->TestBitNumber(maxPadRow)){}
1084 if (maxPadRow>126) return rfirstOROC1+(maxPadRow-126-1)*padrowHeightOROC1;
1085 if (maxPadRow>62) return rfirstOROC0+(maxPadRow-62-1)*padrowHeightOROC0;
1086 if (maxPadRow>0) return rfirstIROC+(maxPadRow-1)*padrowHeightIROC;
1087 return 0.0;
10d100d4 1088}
644666df 1089
1b45e564 1090
644666df 1091//_____________________________________________________________________________
1092Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const
1093{
1094 //calculate the coordinates of the apex of the track
1095 Double_t x[3];
1096 track->GetXYZ(x);
1097 Double_t p[3];
1098 track->GetPxPyPz(p);
1099 Double_t r = 1./track->OneOverPt()/0.0299792458/magField; //signed - will determine the direction of b
1b45e564 1100 //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 1101 //find orthogonal vector (Gram-Schmidt)
1102 Double_t alpha = (p[0]*x[0] + p[1]*x[1])/(p[0]*p[0] + p[1]*p[1]);
1103 Double_t b[2];
1b45e564 1104 b[0]=x[0]-alpha*p[0];
1105 b[1]=x[1]-alpha*p[1];
1106
644666df 1107 Double_t norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
1108 if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
1b45e564 1109 b[0]/=norm;
1110 b[1]/=norm;
1111 b[0]*=r;
1112 b[1]*=r;
1113 b[0]+=x[0];
644666df 1114 b[1]+=x[1];
1115 //printf("center: x:%.2f, y:%.2f\n",b[0],b[1]);
1b45e564 1116
644666df 1117 norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
1118 if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
1119
1120 position[0]=b[0]+b[0]*TMath::Abs(r)/norm;
1121 position[1]=b[1]+b[1]*TMath::Abs(r)/norm;
1122 position[2]=0.;
1123 return kTRUE;
f85a3764 1124}