1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-----------------------------------------------------------------
17 // Implementation of the TPC PID class
18 // Very naive one... Should be made better by the detector experts...
19 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
20 // With many additions and modifications suggested by
21 // Alexander Kalweit, GSI, alexander.philipp.kalweit@cern.ch
22 // Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de
23 // ...and some modifications by
24 // Mikolaj Krzewicki, GSI, mikolaj.krzewicki@cern.ch
25 // ...and some modifications plus eta correction functions by
26 // Benjamin Hess, University of Tuebingen, bhess@cern.ch
27 //-----------------------------------------------------------------
30 #include <TObjArray.h>
37 #include "AliExternalTrackParam.h"
38 #include "AliVTrack.h"
39 #include "AliTPCPIDResponse.h"
40 #include "AliTPCdEdxInfo.h"
44 ClassImp(AliTPCPIDResponse);
47 AliTPCPIDResponse *AliTPCPIDResponse::fgInstance =0;
49 const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]=
51 "", //default - no name
57 //_________________________________________________________________________
58 AliTPCPIDResponse::AliTPCPIDResponse():
69 fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
71 fLowGainIROCthreshold(-40),
72 fBadIROCthreshhold(-70),
73 fLowGainOROCthreshold(-40),
74 fBadOROCthreshhold(-40),
75 fMaxBadLengthFraction(0.5),
80 fCurrentEventMultiplicity(0),
81 fCorrFuncMultiplicity(0x0),
82 fCorrFuncMultiplicityTanTheta(0x0),
83 fCorrFuncSigmaMultiplicity(0x0)
86 // The default constructor
89 AliLog::SetClassDebugLevel("AliTPCPIDResponse", AliLog::kInfo);
91 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;}
93 fCorrFuncMultiplicity = new TF1("fCorrFuncMultiplicity",
94 "[0] + [1]*TMath::Max([4], TMath::Min(x, [3])) + [2] * TMath::Power(TMath::Max([4], TMath::Min(x, [3])), 2)",
96 fCorrFuncMultiplicityTanTheta = new TF1("fCorrFuncMultiplicityTanTheta", "[0] * (x - [2]) + [1] * (x * x - [2] * [2])", -1.5, 1.5);
97 fCorrFuncSigmaMultiplicity = new TF1("fCorrFuncSigmaMultiplicity",
98 "TMath::Max(0, [0] + [1]*TMath::Min(x, [3]) + [2] * TMath::Power(TMath::Min(x, [3]), 2))", 0., 0.2);
101 ResetMultiplicityCorrectionFunctions();
105 //_________________________________________________________________________
106 AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
116 fUseDatabase(kFALSE),
117 fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
119 fLowGainIROCthreshold(-40),
120 fBadIROCthreshhold(-70),
121 fLowGainOROCthreshold(-40),
122 fBadOROCthreshhold(-40),
123 fMaxBadLengthFraction(0.5),
130 // The main constructor
132 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];}
136 //_________________________________________________________________________
137 AliTPCPIDResponse::~AliTPCPIDResponse()
146 delete fhEtaSigmaPar1;
147 fhEtaSigmaPar1 = 0x0;
149 delete fCorrFuncMultiplicity;
150 fCorrFuncMultiplicity = 0x0;
152 delete fCorrFuncMultiplicityTanTheta;
153 fCorrFuncMultiplicityTanTheta = 0x0;
155 delete fCorrFuncSigmaMultiplicity;
156 fCorrFuncSigmaMultiplicity = 0x0;
157 if (fgInstance==this) fgInstance=0;
161 //_________________________________________________________________________
162 AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that):
172 fUseDatabase(that.fUseDatabase),
173 fResponseFunctions(that.fResponseFunctions),
174 fVoltageMap(that.fVoltageMap),
175 fLowGainIROCthreshold(that.fLowGainIROCthreshold),
176 fBadIROCthreshhold(that.fBadIROCthreshhold),
177 fLowGainOROCthreshold(that.fLowGainOROCthreshold),
178 fBadOROCthreshhold(that.fBadOROCthreshhold),
179 fMaxBadLengthFraction(that.fMaxBadLengthFraction),
180 fMagField(that.fMagField),
183 fSigmaPar0(that.fSigmaPar0),
184 fCurrentEventMultiplicity(that.fCurrentEventMultiplicity),
185 fCorrFuncMultiplicity(0x0),
186 fCorrFuncMultiplicityTanTheta(0x0),
187 fCorrFuncSigmaMultiplicity(0x0)
190 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
193 if (that.fhEtaCorr) {
194 fhEtaCorr = new TH2D(*(that.fhEtaCorr));
195 fhEtaCorr->SetDirectory(0);
198 if (that.fhEtaSigmaPar1) {
199 fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
200 fhEtaSigmaPar1->SetDirectory(0);
203 // Copy multiplicity correction functions
204 if (that.fCorrFuncMultiplicity) {
205 fCorrFuncMultiplicity = new TF1(*(that.fCorrFuncMultiplicity));
208 if (that.fCorrFuncMultiplicityTanTheta) {
209 fCorrFuncMultiplicityTanTheta = new TF1(*(that.fCorrFuncMultiplicityTanTheta));
212 if (that.fCorrFuncSigmaMultiplicity) {
213 fCorrFuncSigmaMultiplicity = new TF1(*(that.fCorrFuncSigmaMultiplicity));
217 //_________________________________________________________________________
218 AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
221 if (&that==this) return *this;
222 TNamed::operator=(that);
229 fUseDatabase=that.fUseDatabase;
230 fResponseFunctions=that.fResponseFunctions;
231 fVoltageMap=that.fVoltageMap;
232 fLowGainIROCthreshold=that.fLowGainIROCthreshold;
233 fBadIROCthreshhold=that.fBadIROCthreshhold;
234 fLowGainOROCthreshold=that.fLowGainOROCthreshold;
235 fBadOROCthreshhold=that.fBadOROCthreshhold;
236 fMaxBadLengthFraction=that.fMaxBadLengthFraction;
237 fMagField=that.fMagField;
238 fCurrentEventMultiplicity=that.fCurrentEventMultiplicity;
239 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
243 if (that.fhEtaCorr) {
244 fhEtaCorr = new TH2D(*(that.fhEtaCorr));
245 fhEtaCorr->SetDirectory(0);
248 delete fhEtaSigmaPar1;
250 if (that.fhEtaSigmaPar1) {
251 fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
252 fhEtaSigmaPar1->SetDirectory(0);
255 fSigmaPar0 = that.fSigmaPar0;
257 delete fCorrFuncMultiplicity;
258 fCorrFuncMultiplicity = 0x0;
259 if (that.fCorrFuncMultiplicity) {
260 fCorrFuncMultiplicity = new TF1(*(that.fCorrFuncMultiplicity));
263 delete fCorrFuncMultiplicityTanTheta;
264 fCorrFuncMultiplicityTanTheta = 0x0;
265 if (that.fCorrFuncMultiplicityTanTheta) {
266 fCorrFuncMultiplicityTanTheta = new TF1(*(that.fCorrFuncMultiplicityTanTheta));
269 delete fCorrFuncSigmaMultiplicity;
270 fCorrFuncSigmaMultiplicity = 0x0;
271 if (that.fCorrFuncSigmaMultiplicity) {
272 fCorrFuncSigmaMultiplicity = new TF1(*(that.fCorrFuncSigmaMultiplicity));
278 //_________________________________________________________________________
279 Double_t AliTPCPIDResponse::Bethe(Double_t betaGamma) const {
281 // This is the Bethe-Bloch function normalised to 1 at the minimum
283 // Simulated and reconstructed Bethe-Bloch differs
284 // Simulated curve is the dNprim/dx
285 // Reconstructed is proportianal dNtot/dx
286 // Temporary fix for production - Simple linear correction function
287 // Future 2 Bethe Bloch formulas needed
289 // 2. for reconstructed PID
292 // const Float_t kmeanCorrection =0.1;
294 AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5);
298 //_________________________________________________________________________
299 void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1,
305 // Set the parameters of the ALEPH Bethe-Bloch formula
314 //_________________________________________________________________________
315 void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) {
317 // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint)
319 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;}
322 //_________________________________________________________________________
323 Double_t AliTPCPIDResponse::GetExpectedSignal(Float_t mom,
324 AliPID::EParticleType n) const {
326 // Deprecated function (for backward compatibility). Please use
327 // GetExpectedSignal(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource,
328 // Bool_t correctEta, Bool_t correctMultiplicity);
332 // Calculates the expected PID signal as the function of
333 // the information stored in the track, for the specified particle type
335 // At the moment, these signals are just the results of calling the
336 // Bethe-Bloch formula.
337 // This can be improved. By taking into account the number of
338 // assigned clusters and/or the track dip angle, for example.
341 //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
342 // !!! Splines for light nuclei need to be normalised to this factor !!!
343 const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3);
345 Double_t mass=AliPID::ParticleMassZ(n);
346 if (!fUseDatabase) return Bethe(mom/mass) * chargeFactor;
348 const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
350 if (!responseFunction) return Bethe(mom/mass) * chargeFactor;
352 return fMIP*responseFunction->Eval(mom/mass)*chargeFactor;
356 //_________________________________________________________________________
357 Double_t AliTPCPIDResponse::GetExpectedSigma(Float_t mom,
359 AliPID::EParticleType n) const {
361 // Deprecated function (for backward compatibility). Please use
362 // GetExpectedSigma(onst AliVTrack* track, AliPID::EParticleType species,
363 // ETPCdEdxSource dedxSource, Bool_t correctEta) instead!
366 // Calculates the expected sigma of the PID signal as the function of
367 // the information stored in the track, for the specified particle type
371 return GetExpectedSignal(mom,n)*fRes0[0]*sqrt(1. + fResN2[0]/nPoints);
373 return GetExpectedSignal(mom,n)*fRes0[0];
376 ////////////////////////////////////////////////////NEW//////////////////////////////
378 //_________________________________________________________________________
379 void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario) {
381 // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint)
383 if ((Int_t)gainScenario>=(Int_t)fgkNumberOfGainScenarios) return; //TODO: better handling!
384 fRes0[gainScenario]=res0;
385 fResN2[gainScenario]=resN2;
389 //_________________________________________________________________________
390 Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
391 AliPID::EParticleType species,
393 const TSpline3* responseFunction,
395 Bool_t correctMultiplicity) const
397 // Calculates the expected PID signal as the function of
398 // the information stored in the track and the given parameters,
399 // for the specified particle type
401 // At the moment, these signals are just the results of calling the
402 // Bethe-Bloch formula plus, if desired, taking into account the eta dependence
403 // and the multiplicity dependence (for PbPb).
404 // This can be improved. By taking into account the number of
405 // assigned clusters and/or the track dip angle, for example.
408 Double_t mom=track->GetTPCmomentum();
409 Double_t mass=AliPID::ParticleMassZ(species);
411 //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
412 // !!! Splines for light nuclei need to be normalised to this factor !!!
413 const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3);
415 if (!responseFunction)
416 return Bethe(mom/mass) * chargeFactor;
418 Double_t dEdxSplines = fMIP*responseFunction->Eval(mom/mass) * chargeFactor;
420 if (!correctEta && !correctMultiplicity)
423 Double_t corrFactorEta = 1.0;
424 Double_t corrFactorMultiplicity = 1.0;
427 corrFactorEta = GetEtaCorrectionFast(track, dEdxSplines);
428 //TODO Alternatively take current track dEdx
429 //corrFactorEta = GetEtaCorrectionFast(track, dEdx);
432 if (correctMultiplicity)
433 corrFactorMultiplicity = GetMultiplicityCorrectionFast(track, dEdxSplines * corrFactorEta, fCurrentEventMultiplicity);
435 return dEdxSplines * corrFactorEta * corrFactorMultiplicity;
439 //_________________________________________________________________________
440 Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
441 AliPID::EParticleType species,
442 ETPCdEdxSource dedxSource,
444 Bool_t correctMultiplicity) const
446 // Calculates the expected PID signal as the function of
447 // the information stored in the track, for the specified particle type
449 // At the moment, these signals are just the results of calling the
450 // Bethe-Bloch formula plus, if desired, taking into account the eta dependence
451 // and the multiplicity dependence (for PbPb).
452 // This can be improved. By taking into account the number of
453 // assigned clusters and/or the track dip angle, for example.
457 //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
458 // !!! Splines for light nuclei need to be normalised to this factor !!!
459 const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3);
461 return Bethe(track->GetTPCmomentum() / AliPID::ParticleMassZ(species)) * chargeFactor;
466 ETPCgainScenario gainScenario = kGainScenarioInvalid;
467 TSpline3* responseFunction = 0x0;
469 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) {
470 // Something is wrong with the track -> Return obviously invalid value
474 // Charge factor already taken into account inside the following function call
475 return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity);
478 //_________________________________________________________________________
479 TSpline3* AliTPCPIDResponse::GetResponseFunction( AliPID::EParticleType type,
480 AliTPCPIDResponse::ETPCgainScenario gainScenario ) const
482 //get response function
483 return dynamic_cast<TSpline3*>(fResponseFunctions.At(ResponseFunctionIndex(type,gainScenario)));
486 //_________________________________________________________________________
487 TSpline3* AliTPCPIDResponse::GetResponseFunction( const AliVTrack* track,
488 AliPID::EParticleType species,
489 ETPCdEdxSource dedxSource) const
491 //the splines are stored in an array, different scenarios
495 ETPCgainScenario gainScenario = kGainScenarioInvalid;
496 TSpline3* responseFunction = 0x0;
498 if (ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
499 return responseFunction;
504 //_________________________________________________________________________
505 void AliTPCPIDResponse::ResetSplines()
507 //reset the array with splines
508 for (Int_t i=0;i<fResponseFunctions.GetEntriesFast();i++)
510 fResponseFunctions.AddAt(NULL,i);
513 //_________________________________________________________________________
514 Int_t AliTPCPIDResponse::ResponseFunctionIndex( AliPID::EParticleType species,
515 ETPCgainScenario gainScenario ) const
517 //get the index in fResponseFunctions given type and scenario
518 return Int_t(species)+Int_t(gainScenario)*fgkNumberOfParticleSpecies;
521 //_________________________________________________________________________
522 void AliTPCPIDResponse::SetResponseFunction( TObject* o,
523 AliPID::EParticleType species,
524 ETPCgainScenario gainScenario )
526 fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario));
530 //_________________________________________________________________________
531 Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
532 AliPID::EParticleType species,
533 ETPCgainScenario gainScenario,
536 const TSpline3* responseFunction,
538 Bool_t correctMultiplicity) const
540 // Calculates the expected sigma of the PID signal as the function of
541 // the information stored in the track and the given parameters,
542 // for the specified particle type
545 if (!responseFunction)
548 //TODO Check whether it makes sense to set correctMultiplicity to kTRUE while correctEta might be kFALSE
550 // If eta correction (=> new sigma parametrisation) is requested, but no sigma map is available, print error message
551 if (correctEta && !fhEtaSigmaPar1) {
552 AliError("New sigma parametrisation requested, but sigma map not initialised (usually via AliPIDResponse). Old sigma parametrisation will be used!");
555 // If no sigma map is available or if no eta correction is requested (sigma maps only for corrected eta!), use the old parametrisation
556 if (!fhEtaSigmaPar1 || !correctEta) {
558 return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity) *
559 fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints);
561 return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity)*fRes0[gainScenario];
565 // Use eta correction (+ eta-dependent sigma)
566 Double_t sigmaPar1 = GetSigmaPar1Fast(track, species, dEdx, responseFunction);
568 if (correctMultiplicity) {
569 // In addition, take into account multiplicity dependence of mean and sigma of dEdx
570 Double_t dEdxExpectedEtaCorrected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
572 // GetMultiplicityCorrection and GetMultiplicitySigmaCorrection both need the eta corrected dEdxExpected
573 Double_t multiplicityCorrFactor = GetMultiplicityCorrectionFast(track, dEdxExpectedEtaCorrected, fCurrentEventMultiplicity);
574 Double_t multiplicitySigmaCorrFactor = GetMultiplicitySigmaCorrectionFast(dEdxExpectedEtaCorrected, fCurrentEventMultiplicity);
576 // multiplicityCorrFactor to correct dEdxExpected for multiplicity. In addition: Correction factor for sigma
577 return (dEdxExpectedEtaCorrected * multiplicityCorrFactor)
578 * (sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints) * multiplicitySigmaCorrFactor);
581 return GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE)*
582 sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints);
586 // One should never have/take tracks with 0 dEdx clusters!
592 //_________________________________________________________________________
593 Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
594 AliPID::EParticleType species,
595 ETPCdEdxSource dedxSource,
597 Bool_t correctMultiplicity) const
599 // Calculates the expected sigma of the PID signal as the function of
600 // the information stored in the track, for the specified particle type
606 ETPCgainScenario gainScenario = kGainScenarioInvalid;
607 TSpline3* responseFunction = 0x0;
609 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
610 return 999; //TODO: better handling!
612 return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity);
616 //_________________________________________________________________________
617 Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
618 AliPID::EParticleType species,
619 ETPCdEdxSource dedxSource,
621 Bool_t correctMultiplicity) const
623 //Calculates the number of sigmas of the PID signal from the expected value
624 //for a given particle species in the presence of multiple gain scenarios
629 ETPCgainScenario gainScenario = kGainScenarioInvalid;
630 TSpline3* responseFunction = 0x0;
632 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
633 return -999; //TODO: Better handling!
635 Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity);
636 Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity);
637 // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters
641 return (dEdx-bethe)/sigma;
644 //_________________________________________________________________________
645 Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track,
646 AliPID::EParticleType species,
647 ETPCdEdxSource dedxSource,
649 Bool_t correctMultiplicity,
650 Bool_t ratio/*=kFALSE*/)const
652 //Calculates the number of sigmas of the PID signal from the expected value
653 //for a given particle species in the presence of multiple gain scenarios
658 ETPCgainScenario gainScenario = kGainScenarioInvalid;
659 TSpline3* responseFunction = 0x0;
661 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
662 return -9999.; //TODO: Better handling!
664 const Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity);
666 Double_t delta=-9999.;
667 if (!ratio) delta=dEdx-bethe;
668 else if (bethe>1.e-20) delta=dEdx/bethe;
673 //_________________________________________________________________________
674 Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
675 AliPID::EParticleType species,
676 ETPCdEdxSource dedxSource,
679 ETPCgainScenario& gainScenario,
680 TSpline3** responseFunction) const
682 // Calculates the right parameters for PID
683 // dEdx parametrization for the proper gain scenario, dEdx
684 // and NPoints used for dEdx
685 // based on the track geometry (which chambers it crosses) for the specified particle type
686 // and preferred source of dedx.
687 // returns true on success
690 if (dedxSource == kdEdxDefault) {
691 // Fast handling for default case. In addition: Keep it simple (don't call additional functions) to
692 // avoid possible bugs
694 // GetTPCsignalTunedOnData will be non-positive, if it has not been set (i.e. in case of MC NOT tuned to data).
695 // If this is the case, just take the normal signal
696 dEdx = track->GetTPCsignalTunedOnData();
698 dEdx = track->GetTPCsignal();
701 nPoints = track->GetTPCsignalN();
702 gainScenario = kDefault;
704 TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario));
705 *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
710 //TODO Proper handle of tuneMConData for other dEdx sources
712 Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used)
713 Char_t ncl[3]; //same
714 Char_t nrows[3]; //same
715 const AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo();
717 if (!dEdxInfo && dedxSource!=kdEdxDefault) //in one case its ok if we dont have the info
719 AliError("AliTPCdEdxInfo not available");
723 if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows);
725 //check if we cross a bad OROC in which case we reject
726 EChamberStatus trackOROCStatus = TrackStatus(track,2);
727 if (trackOROCStatus==kChamberOff || trackOROCStatus==kChamberLowGain)
736 if (trackOROCStatus==kChamberInvalid) return kFALSE; //never reached OROC
738 nPoints = ncl[2]+ncl[1];
739 gainScenario = kOROChigh;
744 //if we cross a bad IROC we use OROC dedx, if we dont we use combined
745 EChamberStatus status = TrackStatus(track,1);
746 if (status!=kChamberHighGain)
749 nPoints = ncl[2]+ncl[1];
750 gainScenario = kOROChigh;
754 dEdx = track->GetTPCsignal();
755 nPoints = track->GetTPCsignalN();
756 gainScenario = kALLhigh;
764 gainScenario = kGainScenarioInvalid;
768 TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario));
769 *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
775 //_________________________________________________________________________
776 Double_t AliTPCPIDResponse::GetEtaCorrectionFast(const AliVTrack *track, Double_t dEdxSplines) const
778 // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
780 // Get eta correction for the given parameters.
784 // Calling this function means to request eta correction in some way. Print error message, if no map is available!
785 AliError("Eta correction requested, but map not initialised (usually via AliPIDResponse). Returning eta correction factor 1!");
789 Double_t tpcSignal = dEdxSplines;
794 Double_t tanTheta = GetTrackTanTheta(track);
795 Int_t binX = fhEtaCorr->GetXaxis()->FindFixBin(tanTheta);
796 Int_t binY = fhEtaCorr->GetYaxis()->FindFixBin(1. / tpcSignal);
800 if (binX > fhEtaCorr->GetXaxis()->GetNbins())
801 binX = fhEtaCorr->GetXaxis()->GetNbins();
805 if (binY > fhEtaCorr->GetYaxis()->GetNbins())
806 binY = fhEtaCorr->GetYaxis()->GetNbins();
808 return fhEtaCorr->GetBinContent(binX, binY);
812 //_________________________________________________________________________
813 Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
816 // Get eta correction for the given track.
824 ETPCgainScenario gainScenario = kGainScenarioInvalid;
825 TSpline3* responseFunction = 0x0;
827 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
830 // For the eta correction, do NOT take the multiplicity corrected value of dEdx
831 Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
833 //TODO Alternatively take current track dEdx
834 //return GetEtaCorrectionFast(track, dEdx);
836 return GetEtaCorrectionFast(track, dEdxSplines);
840 //_________________________________________________________________________
841 Double_t AliTPCPIDResponse::GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
844 // Get eta corrected dEdx for the given track. For the correction, the expected dEdx of
845 // the specified species will be used. If the species is set to AliPID::kUnknown, the
846 // dEdx of the track is used instead.
847 // WARNING: In the latter case, the eta correction might not be as good as if the
848 // expected dEdx is used, which is the way the correction factor is designed
850 // In any case, one has to decide carefully to which expected signal one wants to
851 // compare the corrected value - to the corrected or uncorrected.
852 // Anyhow, a safer way of looking e.g. at the n-sigma is to call
853 // the corresponding function GetNumberOfSigmas!
858 ETPCgainScenario gainScenario = kGainScenarioInvalid;
859 TSpline3* responseFunction = 0x0;
861 // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case
862 // it is not used anyway, so this causes no trouble.
863 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
866 Double_t etaCorr = 0.;
868 if (species < AliPID::kUnknown) {
869 // For the eta correction, do NOT take the multiplicity corrected value of dEdx
870 Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
871 etaCorr = GetEtaCorrectionFast(track, dEdxSplines);
874 etaCorr = GetEtaCorrectionFast(track, dEdx);
880 return dEdx / etaCorr;
884 //_________________________________________________________________________
885 Double_t AliTPCPIDResponse::GetSigmaPar1Fast(const AliVTrack *track, AliPID::EParticleType species, Double_t dEdx,
886 const TSpline3* responseFunction) const
888 // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
890 // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track.
893 if (!fhEtaSigmaPar1) {
894 // Calling this function means to request new sigma parametrisation in some way. Print error message, if no map is available!
895 AliError("New sigma parametrisation requested, but sigma map not initialised (usually via AliPIDResponse). Returning error value for sigma parameter1 = 999!");
899 // The sigma maps are created with data sets that are already eta corrected and for which the
900 // splines have been re-created. Therefore, the value for the lookup needs to be the value of
901 // the splines without any additional eta correction.
902 // NOTE: This is due to the method the maps are created. The track dEdx (not the expected one!)
903 // is corrected to uniquely related a momemtum bin with an expected dEdx, where the expected dEdx
904 // equals the track dEdx for all eta (thanks to the correction and the re-fit of the splines).
905 // Consequently, looking up the uncorrected expected dEdx at a given tanTheta yields the correct
907 // Also: It has to be the spline dEdx, since one wants to get the sigma for the assumption(!)
908 // of such a particle, which by assumption then has this dEdx value
910 // For the eta correction, do NOT take the multiplicity corrected value of dEdx
911 Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
913 if (dEdxExpected < 1.)
916 Double_t tanTheta = GetTrackTanTheta(track);
917 Int_t binX = fhEtaSigmaPar1->GetXaxis()->FindFixBin(tanTheta);
918 Int_t binY = fhEtaSigmaPar1->GetYaxis()->FindFixBin(1. / dEdxExpected);
922 if (binX > fhEtaSigmaPar1->GetXaxis()->GetNbins())
923 binX = fhEtaSigmaPar1->GetXaxis()->GetNbins();
927 if (binY > fhEtaSigmaPar1->GetYaxis()->GetNbins())
928 binY = fhEtaSigmaPar1->GetYaxis()->GetNbins();
930 return fhEtaSigmaPar1->GetBinContent(binX, binY);
934 //_________________________________________________________________________
935 Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
938 // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track.
946 ETPCgainScenario gainScenario = kGainScenarioInvalid;
947 TSpline3* responseFunction = 0x0;
949 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
952 return GetSigmaPar1Fast(track, species, dEdx, responseFunction);
956 //_________________________________________________________________________
957 Bool_t AliTPCPIDResponse::SetEtaCorrMap(TH2D* hMap)
960 // Load map for TPC eta correction (a copy is stored and will be deleted automatically).
961 // If hMap is 0x0,the eta correction will be disabled and kFALSE is returned.
962 // If the map can be set, kTRUE is returned.
973 fhEtaCorr = (TH2D*)(hMap->Clone());
974 fhEtaCorr->SetDirectory(0);
980 //_________________________________________________________________________
981 Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0)
984 // Load map for TPC sigma map (a copy is stored and will be deleted automatically):
985 // Parameter 1 is stored as a 2D map (1/dEdx vs. tanTheta_local) and parameter 0 is
986 // a constant. If hSigmaPar1Map is 0x0, the old sigma parametrisation will be used
987 // (and sigmaPar0 is ignored!) and kFALSE is returned.
988 // If the map can be set, sigmaPar0 is also set and kTRUE will be returned.
991 delete fhEtaSigmaPar1;
993 if (!hSigmaPar1Map) {
994 fhEtaSigmaPar1 = 0x0;
1000 fhEtaSigmaPar1 = (TH2D*)(hSigmaPar1Map->Clone());
1001 fhEtaSigmaPar1->SetDirectory(0);
1002 fSigmaPar0 = sigmaPar0;
1008 //_________________________________________________________________________
1009 Double_t AliTPCPIDResponse::GetTrackTanTheta(const AliVTrack *track) const
1011 // Extract the tanTheta from the information available in the AliVTrack
1013 // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()).
1014 // However, this value is not available for AODs and, thus, not for AliVTrack.
1015 // Fortunately, the following formula allows to approximate the local tanTheta with the
1016 // global theta angle -> This is for by far most of the tracks the same, but gives at
1017 // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok.
1020 const AliExternalTrackParam* innerParam = track->GetInnerParam();
1021 Double_t tanTheta = 0;
1023 tanTheta = innerParam->GetTgl();
1025 tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
1027 // Constant in formula for B in kGauss (factor 0.1 to convert B from Tesla to kGauss),
1028 // pT in GeV/c (factor c*1e-9 to arrive at GeV/c) and curvature in 1/cm (factor 0.01 to get from m to cm)
1029 const Double_t constant = TMath::C()* 1e-9 * 0.1 * 0.01;
1030 const Double_t curvature = fMagField * constant / track->Pt(); // in 1./cm
1032 Double_t averageddzdr = 0.;
1035 for (Double_t r = 85; r < 245; r++) {
1036 Double_t sinPhiLocal = TMath::Abs(r*curvature*0.5);
1038 // Cut on |sin(phi)| as during reco
1039 if (TMath::Abs(sinPhiLocal) <= 0.95) {
1040 const Double_t phiLocal = TMath::ASin(sinPhiLocal);
1041 const Double_t tanPhiLocal = TMath::Tan(phiLocal);
1043 averageddzdr += tanTheta * TMath::Sqrt(1. + tanPhiLocal * tanPhiLocal);
1049 averageddzdr /= nParts;
1051 AliError("Problems during determination of dz/dr. Returning pure tanTheta as best estimate!");
1055 //printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\ntanThetaGlobalFromTheta/tanTheta/Averageddzdr: %f / %f / %f\n\n",
1056 // track->Pt(), constant, fMagField, 1./curvature, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, averageddzdr);
1058 return averageddzdr;
1062 // Alternatively (in average, the effect is found to be negligable!):
1063 // Take local tanTheta from TPC inner wall, if available (currently only for ESDs available)
1064 //const AliExternalTrackParam* innerParam = track->GetInnerParam();
1066 // return innerParam->GetTgl();
1069 return TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
1073 //_________________________________________________________________________
1074 Double_t AliTPCPIDResponse::GetMultiplicityCorrectionFast(const AliVTrack *track, Double_t dEdxExpected, Int_t multiplicity) const
1076 // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
1078 // Calculate the multiplicity correction factor for this track for the given multiplicity.
1079 // The parameter dEdxExpected should take into account the eta correction already!
1081 // Multiplicity depends on pure dEdx. Therefore, correction factor depends indirectly on eta
1082 // => Use eta corrected value for dEdexExpected.
1084 if (dEdxExpected <= 0 || multiplicity <= 0)
1087 const Double_t dEdxExpectedInv = 1. / dEdxExpected;
1088 Double_t relSlope = fCorrFuncMultiplicity->Eval(dEdxExpectedInv);
1090 const Double_t tanTheta = GetTrackTanTheta(track);
1091 relSlope += fCorrFuncMultiplicityTanTheta->Eval(tanTheta);
1093 return (1. + relSlope * multiplicity);
1097 //_________________________________________________________________________
1098 Double_t AliTPCPIDResponse::GetMultiplicityCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
1101 // Get multiplicity correction for the given track (w.r.t. the mulitplicity of the current event)
1104 //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse
1106 // No correction in case of multiplicity <= 0
1107 if (fCurrentEventMultiplicity <= 0)
1112 ETPCgainScenario gainScenario = kGainScenarioInvalid;
1113 TSpline3* responseFunction = 0x0;
1115 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
1118 //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid?
1119 // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course)
1120 Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
1122 return GetMultiplicityCorrectionFast(track, dEdxExpected, fCurrentEventMultiplicity);
1126 //_________________________________________________________________________
1127 Double_t AliTPCPIDResponse::GetMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
1130 // Get multiplicity corrected dEdx for the given track. For the correction, the expected dEdx of
1131 // the specified species will be used. If the species is set to AliPID::kUnknown, the
1132 // dEdx of the track is used instead.
1133 // WARNING: In the latter case, the correction might not be as good as if the
1134 // expected dEdx is used, which is the way the correction factor is designed
1136 // In any case, one has to decide carefully to which expected signal one wants to
1137 // compare the corrected value - to the corrected or uncorrected.
1138 // Anyhow, a safer way of looking e.g. at the n-sigma is to call
1139 // the corresponding function GetNumberOfSigmas!
1142 //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse
1146 ETPCgainScenario gainScenario = kGainScenarioInvalid;
1147 TSpline3* responseFunction = 0x0;
1149 // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case
1150 // it is not used anyway, so this causes no trouble.
1151 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
1155 // No correction in case of multiplicity <= 0
1156 if (fCurrentEventMultiplicity <= 0)
1159 Double_t multiplicityCorr = 0.;
1161 // TODO Normally, one should use the eta corrected values in BOTH of the following cases. Does it make sense to use the uncorrected dEdx values?
1162 if (species < AliPID::kUnknown) {
1163 // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course).
1164 // However, one needs the eta corrected value!
1165 Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
1166 multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdxSplines, fCurrentEventMultiplicity);
1169 // One needs the eta corrected value to determine the multiplicity correction factor!
1170 Double_t etaCorr = GetEtaCorrectionFast(track, dEdx);
1171 multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdx * etaCorr, fCurrentEventMultiplicity);
1174 if (multiplicityCorr <= 0)
1177 return dEdx / multiplicityCorr;
1181 //_________________________________________________________________________
1182 Double_t AliTPCPIDResponse::GetEtaAndMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species,
1183 ETPCdEdxSource dedxSource) const
1186 // Get multiplicity and eta corrected dEdx for the given track. For the correction,
1187 // the expected dEdx of the specified species will be used. If the species is set
1188 // to AliPID::kUnknown, the dEdx of the track is used instead.
1189 // WARNING: In the latter case, the correction might not be as good as if the
1190 // expected dEdx is used, which is the way the correction factor is designed
1192 // In any case, one has to decide carefully to which expected signal one wants to
1193 // compare the corrected value - to the corrected or uncorrected.
1194 // Anyhow, a safer way of looking e.g. at the n-sigma is to call
1195 // the corresponding function GetNumberOfSigmas!
1200 ETPCgainScenario gainScenario = kGainScenarioInvalid;
1201 TSpline3* responseFunction = 0x0;
1203 // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case
1204 // it is not used anyway, so this causes no trouble.
1205 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
1208 Double_t multiplicityCorr = 0.;
1209 Double_t etaCorr = 0.;
1211 if (species < AliPID::kUnknown) {
1212 // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course)
1213 Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
1214 etaCorr = GetEtaCorrectionFast(track, dEdxSplines);
1215 multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdxSplines * etaCorr, fCurrentEventMultiplicity);
1218 etaCorr = GetEtaCorrectionFast(track, dEdx);
1219 multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdx * etaCorr, fCurrentEventMultiplicity);
1222 if (multiplicityCorr <= 0 || etaCorr <= 0)
1225 return dEdx / multiplicityCorr / etaCorr;
1229 //_________________________________________________________________________
1230 Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrectionFast(Double_t dEdxExpected, Int_t multiplicity) const
1232 // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
1234 // Calculate the multiplicity sigma correction factor for the corresponding expected dEdx and for the given multiplicity.
1235 // The parameter dEdxExpected should take into account the eta correction already!
1237 // Multiplicity dependence of sigma depends on the real dEdx at zero multiplicity,
1238 // i.e. the eta (only) corrected dEdxExpected value has to be used
1239 // since all maps etc. have been created for ~zero multiplicity
1241 if (dEdxExpected <= 0 || multiplicity <= 0)
1244 Double_t relSigmaSlope = fCorrFuncSigmaMultiplicity->Eval(1. / dEdxExpected);
1246 return (1. + relSigmaSlope * multiplicity);
1250 //_________________________________________________________________________
1251 Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
1254 // Get multiplicity sigma correction for the given track (w.r.t. the mulitplicity of the current event)
1257 //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse
1259 // No correction in case of multiplicity <= 0
1260 if (fCurrentEventMultiplicity <= 0)
1265 ETPCgainScenario gainScenario = kGainScenarioInvalid;
1266 TSpline3* responseFunction = 0x0;
1268 if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
1271 //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid?
1272 // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course)
1273 Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
1275 return GetMultiplicitySigmaCorrectionFast(dEdxExpected, fCurrentEventMultiplicity);
1279 //_________________________________________________________________________
1280 void AliTPCPIDResponse::ResetMultiplicityCorrectionFunctions()
1282 // Default values: No correction, i.e. overall correction factor should be one
1284 // This function should always return 0.0, if multcorr disabled
1285 fCorrFuncMultiplicity->SetParameter(0, 0.);
1286 fCorrFuncMultiplicity->SetParameter(1, 0.);
1287 fCorrFuncMultiplicity->SetParameter(2, 0.);
1288 fCorrFuncMultiplicity->SetParameter(3, 0.);
1289 fCorrFuncMultiplicity->SetParameter(4, 0.);
1291 // This function should always return 0., if multCorr disabled
1292 fCorrFuncMultiplicityTanTheta->SetParameter(0, 0.);
1293 fCorrFuncMultiplicityTanTheta->SetParameter(1, 0.);
1294 fCorrFuncMultiplicityTanTheta->SetParameter(2, 0.);
1296 // This function should always return 0.0, if mutlcorr disabled
1297 fCorrFuncSigmaMultiplicity->SetParameter(0, 0.);
1298 fCorrFuncSigmaMultiplicity->SetParameter(1, 0.);
1299 fCorrFuncSigmaMultiplicity->SetParameter(2, 0.);
1300 fCorrFuncSigmaMultiplicity->SetParameter(3, 0.);
1304 //_________________________________________________________________________
1305 Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner,
1306 Double_t* trackPositionOuter,
1312 //calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses
1313 //for OROC chamber numbers add 36
1314 //returned angles are between (0,2pi)
1316 inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]);
1317 outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]);
1319 if (inphi<0) {inphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
1320 if (outphi<0) {outphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
1322 in = sectorNumber(inphi);
1323 out = sectorNumber(outphi);
1325 //for the C side (positive z) offset by 18
1326 if (trackPositionInner[2]>0.0)
1335 //_____________________________________________________________________________
1336 Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const
1338 //calculate sector number
1339 const Float_t width=TMath::TwoPi()/18.;
1340 return TMath::Floor(phi/width);
1344 //_____________________________________________________________________________
1345 void AliTPCPIDResponse::Print(Option_t* /*option*/) const
1348 fResponseFunctions.Print();
1352 //_____________________________________________________________________________
1353 AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const
1355 //status of the track: if it crosses any bad chambers return kChamberOff;
1356 //IROC:layer=1, OROC:layer=2
1357 if (layer<1 || layer>2) layer=1;
1362 Float_t innerRadius = (layer==1)?83.0:133.7;
1363 Float_t outerRadius = (layer==1)?133.5:247.7;
1365 /////////////////////////////////////////////////////////////////////////////
1366 //find out where track enters and leaves the layer.
1368 Double_t trackPositionInner[3];
1369 Double_t trackPositionOuter[3];
1371 //if there is no inner param this could mean we're using the AOD track,
1372 //we still can extrapolate from the vertex - so use those params.
1373 const AliExternalTrackParam* ip = track->GetInnerParam();
1374 if (ip) track=dynamic_cast<const AliVTrack*>(ip);
1376 Bool_t trackAtInner = track->GetXYZAt(innerRadius, fMagField, trackPositionInner);
1377 Bool_t trackAtOuter = track->GetXYZAt(outerRadius, fMagField, trackPositionOuter);
1381 //if we dont even enter inner radius we do nothing and return invalid
1386 return kChamberInvalid;
1391 //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position
1392 Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter);
1393 Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]);
1394 if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius)
1396 //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]);
1404 return kChamberInvalid;
1409 if (!sectorNumbersInOut(trackPositionInner,
1414 out)) return kChamberInvalid;
1416 /////////////////////////////////////////////////////////////////////////////
1417 //now we have the location of the track we can check
1418 //if it is in a good/bad chamber
1420 Bool_t sideA = kTRUE;
1422 if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE;
1427 if (TMath::Abs(in-out)>9)
1429 if (TMath::Max(in,out)==out)
1434 Float_t tmpphi=outphi;
1439 inphi-=TMath::TwoPi();
1443 if (TMath::Max(in,out)==in)
1448 Float_t tmpphi=outphi;
1454 Float_t trackLengthInBad = 0.;
1455 Float_t trackLengthInLowGain = 0.;
1456 Float_t trackLengthTotal = TMath::Abs(outphi-inphi);
1457 Float_t lengthFractionInBadSectors = 0.;
1459 const Float_t sectorWidth = TMath::TwoPi()/18.;
1461 for (Int_t i=in; i<=out; i++)
1464 if (i<0) j+=18; //correct for the negative values
1465 if (!sideA) j+=18; //move to the correct side
1467 Float_t deltaPhi = 0.;
1468 Float_t phiEdge=sectorWidth*i;
1469 if (inphi>phiEdge) {deltaPhi=phiEdge+sectorWidth-inphi;}
1470 else if ((outphi>=phiEdge) && (outphi<(phiEdge+sectorWidth))) {deltaPhi=outphi-phiEdge;}
1471 else {deltaPhi=sectorWidth;}
1473 Float_t v = fVoltageMap[(layer==1)?(j):(j+36)];
1474 if (v<=fBadIROCthreshhold)
1476 trackLengthInBad+=deltaPhi;
1477 lengthFractionInBadSectors=1.;
1479 if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold)
1481 trackLengthInLowGain+=deltaPhi;
1482 lengthFractionInBadSectors=1.;
1486 //for now low gain and bad (off) chambers are treated equally
1487 if (trackLengthTotal>0)
1488 lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal;
1490 //printf("### side: %s, pt: %.2f, pz: %.2f, in: %i, out: %i, phiIN: %.2f, phiOUT: %.2f, rIN: %.2f, rOUT: %.2f\n",(sideA)?"A":"C",track->Pt(),track->Pz(),in,out,inphi,outphi,innerRadius,outerRadius);
1492 if (lengthFractionInBadSectors>fMaxBadLengthFraction)
1494 //printf("%%%%%%%% %s kChamberLowGain\n",(layer==1)?"IROC":"OROC");
1495 return kChamberLowGain;
1498 return kChamberHighGain;
1502 //_____________________________________________________________________________
1503 Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const
1505 //return the radius of the outermost padrow containing a cluster in TPC
1507 const TBits* clusterMap=track->GetTPCClusterMapPtr();
1508 if (!clusterMap) return 0.;
1510 //from AliTPCParam, radius of first IROC row
1511 const Float_t rfirstIROC = 8.52249984741210938e+01;
1512 const Float_t padrowHeightIROC = 0.75;
1513 const Float_t rfirstOROC0 = 1.35100006103515625e+02;
1514 const Float_t padrowHeightOROC0 = 1.0;
1515 const Float_t rfirstOROC1 = 1.99350006103515625e+02;
1516 const Float_t padrowHeightOROC1 = 1.5;
1518 Int_t maxPadRow=160;
1519 while ((--maxPadRow)>0 && !clusterMap->TestBitNumber(maxPadRow)){}
1520 if (maxPadRow>126) return rfirstOROC1+(maxPadRow-126-1)*padrowHeightOROC1;
1521 if (maxPadRow>62) return rfirstOROC0+(maxPadRow-62-1)*padrowHeightOROC0;
1522 if (maxPadRow>0) return rfirstIROC+(maxPadRow-1)*padrowHeightIROC;
1527 //_____________________________________________________________________________
1528 Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const
1530 //calculate the coordinates of the apex of the track
1534 track->GetPxPyPz(p);
1535 Double_t r = 1./track->OneOverPt()/0.0299792458/magField; //signed - will determine the direction of b
1536 //printf("b: %.2f, x:%.2f, y:%.2f, pt: %.2f, px:%.2f, py%.2f, r: %.2f\n",magField, x[0],x[1],track->Pt(), p[0],p[1],r);
1537 //find orthogonal vector (Gram-Schmidt)
1538 Double_t alpha = (p[0]*x[0] + p[1]*x[1])/(p[0]*p[0] + p[1]*p[1]);
1540 b[0]=x[0]-alpha*p[0];
1541 b[1]=x[1]-alpha*p[1];
1543 Double_t norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
1544 if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
1551 //printf("center: x:%.2f, y:%.2f\n",b[0],b[1]);
1553 norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
1554 if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
1556 position[0]=b[0]+b[0]*TMath::Abs(r)/norm;
1557 position[1]=b[1]+b[1]*TMath::Abs(r)/norm;
1562 Double_t AliTPCPIDResponse::EvaldEdxSpline(Double_t bg,Int_t entry){
1564 // Evaluate the dEdx response for given entry
1566 TSpline * spline = (TSpline*)fSplineArray.At(entry);
1567 if (spline) return spline->Eval(bg);
1572 Bool_t AliTPCPIDResponse::RegisterSpline(const char * name, Int_t index){
1574 // register spline to be used for drawing comparisons
1576 TFile * fTPCBB = TFile::Open("$ALICE_ROOT/OADB/COMMON/PID/data/TPCPIDResponse.root");
1577 TObjArray *arrayTPCPID= (TObjArray*) fTPCBB->Get("TPCPIDResponse");
1578 if (fSplineArray.GetEntriesFast()<index) fSplineArray.Expand(index*2);
1581 spline = (TSpline3*)arrayTPCPID->FindObject(name);
1582 if (spline) fSplineArray.AddAt(spline->Clone(),index);