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 //-----------------------------------------------------------------
28 #include <TObjArray.h>
34 #include "AliExternalTrackParam.h"
35 #include "AliVTrack.h"
36 #include "AliTPCPIDResponse.h"
37 #include "AliTPCdEdxInfo.h"
39 ClassImp(AliTPCPIDResponse);
41 const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]=
43 "", //default - no name
49 //_________________________________________________________________________
50 AliTPCPIDResponse::AliTPCPIDResponse():
61 fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
63 fLowGainIROCthreshold(-40),
64 fBadIROCthreshhold(-70),
65 fLowGainOROCthreshold(-40),
66 fBadOROCthreshhold(-40),
67 fMaxBadLengthFraction(0.5),
68 fCurrentResponseFunction(NULL),
71 fCurrentGainScenario(kGainScenarioInvalid),
75 // The default constructor
77 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;}
80 //_________________________________________________________________________
81 AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
92 fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
94 fLowGainIROCthreshold(-40),
95 fBadIROCthreshhold(-70),
96 fLowGainOROCthreshold(-40),
97 fBadOROCthreshhold(-40),
98 fMaxBadLengthFraction(0.5),
99 fCurrentResponseFunction(NULL),
102 fCurrentGainScenario(kGainScenarioInvalid),
106 // The main constructor
108 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];}
111 //_________________________________________________________________________
112 AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that):
122 fUseDatabase(that.fUseDatabase),
123 fResponseFunctions(that.fResponseFunctions),
124 fVoltageMap(that.fVoltageMap),
125 fLowGainIROCthreshold(that.fLowGainIROCthreshold),
126 fBadIROCthreshhold(that.fBadIROCthreshhold),
127 fLowGainOROCthreshold(that.fLowGainOROCthreshold),
128 fBadOROCthreshhold(that.fBadOROCthreshhold),
129 fMaxBadLengthFraction(that.fMaxBadLengthFraction),
130 fCurrentResponseFunction(NULL),
133 fCurrentGainScenario(kGainScenarioInvalid),
134 fMagField(that.fMagField)
137 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
140 //_________________________________________________________________________
141 AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
144 if (&that==this) return *this;
145 TNamed::operator=(that);
152 fUseDatabase=that.fUseDatabase;
153 fResponseFunctions=that.fResponseFunctions;
154 fVoltageMap=that.fVoltageMap;
155 fLowGainIROCthreshold=that.fLowGainIROCthreshold;
156 fBadIROCthreshhold=that.fBadIROCthreshhold;
157 fLowGainOROCthreshold=that.fLowGainOROCthreshold;
158 fBadOROCthreshhold=that.fBadOROCthreshhold;
159 fMaxBadLengthFraction=that.fMaxBadLengthFraction;
160 fMagField=that.fMagField;
161 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
165 //_________________________________________________________________________
166 Double_t AliTPCPIDResponse::Bethe(Double_t betaGamma) const {
168 // This is the Bethe-Bloch function normalised to 1 at the minimum
170 // Simulated and reconstructed Bethe-Bloch differs
171 // Simulated curve is the dNprim/dx
172 // Reconstructed is proportianal dNtot/dx
173 // Temporary fix for production - Simple linear correction function
174 // Future 2 Bethe Bloch formulas needed
176 // 2. for reconstructed PID
179 // const Float_t kmeanCorrection =0.1;
181 AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5);
185 //_________________________________________________________________________
186 void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1,
192 // Set the parameters of the ALEPH Bethe-Bloch formula
201 //_________________________________________________________________________
202 void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) {
204 // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint)
206 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;}
209 //_________________________________________________________________________
210 Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom,
211 AliPID::EParticleType n) const {
213 // Calculates the expected PID signal as the function of
214 // the information stored in the track, for the specified particle type
216 // At the moment, these signals are just the results of calling the
217 // Bethe-Bloch formula.
218 // This can be improved. By taking into account the number of
219 // assigned clusters and/or the track dip angle, for example.
222 Double_t mass=AliPID::ParticleMassZ(n);
223 if (!fUseDatabase) return Bethe(mom/mass);
225 const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
227 if (!responseFunction) return Bethe(mom/mass);
228 //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
229 // !!! Splines for light nuclei need to be normalised to this factor !!!
230 const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3);
231 return fMIP*responseFunction->Eval(mom/mass)*chargeFactor;
235 //_________________________________________________________________________
236 Double_t AliTPCPIDResponse::GetExpectedSigma(const Float_t mom,
238 AliPID::EParticleType n) const {
240 // Calculates the expected sigma of the PID signal as the function of
241 // the information stored in the track, for the specified particle type
245 return GetExpectedSignal(mom,n)*fRes0[0]*sqrt(1. + fResN2[0]/nPoints);
247 return GetExpectedSignal(mom,n)*fRes0[0];
250 ////////////////////////////////////////////////////NEW//////////////////////////////
252 //_________________________________________________________________________
253 void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario) {
255 // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint)
257 if ((Int_t)gainScenario>(Int_t)fgkNumberOfGainScenarios) return; //TODO: better handling!
258 fRes0[gainScenario]=res0;
259 fResN2[gainScenario]=resN2;
262 //_________________________________________________________________________
263 Double_t AliTPCPIDResponse::GetExpectedSignal(Double_t mom,
264 AliPID::EParticleType species,
265 const TSpline3* responseFunction) const
267 // Calculates the expected PID signal as the function of
268 // the information stored in the track, for the specified particle type
270 // At the moment, these signals are just the results of calling the
271 // Bethe-Bloch formula.
272 // This can be improved. By taking into account the number of
273 // assigned clusters and/or the track dip angle, for example.
277 Double_t mass=AliPID::ParticleMass(species);
279 if (!fUseDatabase) return Bethe(mom/mass);
281 if (!responseFunction) return Bethe(mom/mass);
282 return fMIP*responseFunction->Eval(mom/mass);
285 //_________________________________________________________________________
286 Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
287 AliPID::EParticleType species,
288 ETPCdEdxSource dedxSource)
290 // Calculates the expected PID signal as the function of
291 // the information stored in the track, for the specified particle type
293 // At the moment, these signals are just the results of calling the
294 // Bethe-Bloch formula.
295 // This can be improved. By taking into account the number of
296 // assigned clusters and/or the track dip angle, for example.
299 Double_t mom = track->GetTPCmomentum();
300 Double_t mass=AliPID::ParticleMass(species);
302 if (!fUseDatabase) return Bethe(mom/mass);
304 const TSpline3* responseFunction = GetResponseFunction(track,species,dedxSource);
305 if (!responseFunction) return Bethe(mom/mass);
306 return fMIP*responseFunction->Eval(mom/mass);
310 //_________________________________________________________________________
311 TSpline3* AliTPCPIDResponse::GetResponseFunction( AliPID::EParticleType type,
312 AliTPCPIDResponse::ETPCgainScenario gainScenario ) const
314 //get response function
315 return dynamic_cast<TSpline3*>(fResponseFunctions.At(ResponseFunctionIndex(type,gainScenario)));
318 //_________________________________________________________________________
319 TSpline3* AliTPCPIDResponse::GetResponseFunction( const AliVTrack* track,
320 AliPID::EParticleType species,
321 ETPCdEdxSource dedxSource )
323 //the splines are stored in an array, different scenarios
325 if (ResponseFunctiondEdxN(track,
328 return fCurrentResponseFunction;
333 //_________________________________________________________________________
334 void AliTPCPIDResponse::ResetSplines()
336 //reset the array with splines
337 for (Int_t i=0;i<fResponseFunctions.GetEntriesFast();i++)
339 fResponseFunctions.AddAt(NULL,i);
342 //_________________________________________________________________________
343 Int_t AliTPCPIDResponse::ResponseFunctionIndex( AliPID::EParticleType species,
344 ETPCgainScenario gainScenario ) const
346 //get the index in fResponseFunctions given type and scenario
347 return Int_t(species)+Int_t(gainScenario)*fgkNumberOfParticleSpecies;
350 //_________________________________________________________________________
351 void AliTPCPIDResponse::SetResponseFunction( TObject* o,
352 AliPID::EParticleType species,
353 ETPCgainScenario gainScenario )
355 fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario));
358 //_________________________________________________________________________
359 Double_t AliTPCPIDResponse::GetExpectedSigma( Double_t mom,
361 AliPID::EParticleType species,
362 ETPCgainScenario gainScenario,
363 const TSpline3* responseFunction) const
365 // Calculates the expected sigma of the PID signal as the function of
366 // the information stored in the track, for the specified particle type
370 if (!responseFunction) return 999;
372 return GetExpectedSignal(mom,species,responseFunction) *
373 fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints);
375 return GetExpectedSignal(mom,species,responseFunction)*fRes0[gainScenario];
378 //_________________________________________________________________________
379 Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
380 AliPID::EParticleType species,
381 ETPCdEdxSource dedxSource)
383 // Calculates the expected sigma of the PID signal as the function of
384 // the information stored in the track, for the specified particle type
388 if (!ResponseFunctiondEdxN(track,
390 dedxSource)) return 998; //TODO: better handling!
392 return GetExpectedSigma( track->GetTPCmomentum(),
395 fCurrentGainScenario,
396 fCurrentResponseFunction );
399 //_________________________________________________________________________
400 Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
401 AliPID::EParticleType species,
402 ETPCdEdxSource dedxSource)
404 //calculates the number of sigmas of the PID signal from the expected value
405 //for a gicven particle species in the presence of multiple gain scenarios
408 if (!ResponseFunctiondEdxN(track,
410 dedxSource)) return 997; //TODO: better handling!
412 Double_t mom = track->GetTPCmomentum();
413 Double_t bethe=GetExpectedSignal(mom,species,fCurrentResponseFunction);
414 Double_t sigma=GetExpectedSigma( mom,
417 fCurrentGainScenario,
418 fCurrentResponseFunction );
419 return (fCurrentdEdx-bethe)/sigma;
422 //_________________________________________________________________________
423 Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
424 AliPID::EParticleType species,
425 ETPCdEdxSource dedxSource )
427 // Calculates the right parameters for PID
428 // dEdx parametrization for the proper gain scenario, dEdx
429 // and NPoints used for dEdx
430 // based on the track geometry (which chambers it crosses) for the specified particle type
431 // and preferred source of dedx.
432 // returns true on success
434 Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used)
435 Char_t ncl[3]; //same
436 Char_t nrows[3]; //same
437 AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo();
439 if (!dEdxInfo && dedxSource!=kdEdxDefault) //in one case its ok if we dont have the info
441 AliError("AliTPCdEdxInfo not available");
442 InvalidateCurrentValues();
446 if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows);
448 //check if we cross a bad OROC in which case we reject
449 EChamberStatus trackStatus = TrackStatus(track,2);
450 if (trackStatus==kChamberOff || trackStatus==kChamberLowGain)
452 InvalidateCurrentValues();
460 fCurrentdEdx = track->GetTPCsignal();
461 fCurrentNPoints = track->GetTPCsignalN();
462 fCurrentGainScenario = kDefault;
467 fCurrentdEdx = signal[3];
468 fCurrentNPoints = ncl[2]+ncl[1];
469 fCurrentGainScenario = kOROChigh;
474 //if we cross a bad IROC we use OROC dedx, if we dont we use combined
475 EChamberStatus status = TrackStatus(track,1);
476 if (status!=kChamberHighGain)
478 fCurrentdEdx = signal[3];
479 fCurrentNPoints = ncl[2]+ncl[1];
480 fCurrentGainScenario = kOROChigh;
484 fCurrentdEdx = track->GetTPCsignal();
485 fCurrentNPoints = track->GetTPCsignalN();
486 fCurrentGainScenario = kALLhigh;
494 fCurrentGainScenario = kGainScenarioInvalid;
498 TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,fCurrentGainScenario));
499 fCurrentResponseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
503 //_________________________________________________________________________
504 Bool_t AliTPCPIDResponse::sectorNumbersInOut(const AliVTrack* track,
505 Double_t innerRadius,
506 Double_t outerRadius,
512 //calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses
513 //for OROC chamber numbers add 36
514 //returned angles are between (0,2pi)
516 Double_t trackPositionInner[3];
517 Double_t trackPositionOuter[3];
519 Bool_t trackAtInner=kTRUE;
520 Bool_t trackAtOuter=kTRUE;
521 const AliExternalTrackParam* ip = track->GetInnerParam();
523 //if there is no inner param this could mean we're using the AOD track,
524 //we still can extrapolate from the vertex - so use those params.
527 if (!track->GetXYZAt(innerRadius, fMagField, trackPositionInner))
530 if (!track->GetXYZAt(outerRadius, fMagField, trackPositionOuter))
535 //if we dont even enter inner radius we do nothing
545 //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position
546 Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter);
547 Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]);
548 if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius)
550 //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]);
562 inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]);
563 outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]);
565 if (inphi<0) {inphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
566 if (outphi<0) {outphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
568 in = sectorNumber(inphi);
569 out = sectorNumber(outphi);
571 //for the C side (positive z) offset by 18
572 if (trackPositionInner[2]>0.0)
580 //_____________________________________________________________________________
581 Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const
583 //calculate sector number
584 const Float_t width=TMath::TwoPi()/18.;
585 return TMath::Floor(phi/width);
588 //_____________________________________________________________________________
589 void AliTPCPIDResponse::Print(Option_t* /*option*/) const
592 fResponseFunctions.Print();
595 //_____________________________________________________________________________
596 AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const
598 //status of the track: if it crosses any bad chambers return kChamberOff;
599 //IROC:layer=1, OROC:layer=2
600 if (layer<1 || layer>2) layer=1;
605 Float_t innerRadius = (layer==1)?83.0:133.7;
606 Float_t outerRadius = (layer==1)?133.5:247.7;
607 if (!sectorNumbersInOut(track,
613 out)) return kChamberOff;
615 Bool_t sideA = kTRUE;
617 if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE;
622 if (TMath::Abs(in-out)>9)
624 if (TMath::Max(in,out)==out)
629 Float_t tmpphi=outphi;
634 inphi-=TMath::TwoPi();
638 if (TMath::Max(in,out)==in)
643 Float_t tmpphi=outphi;
649 Float_t trackLengthInBad = 0.;
650 Float_t trackLengthInLowGain = 0.;
651 Float_t trackLengthTotal = TMath::Abs(outphi-inphi);
653 const Float_t sectorWidth = TMath::TwoPi()/18.;
655 for (Int_t i=in; i<=out; i++)
658 if (i<0) j+=18; //correct for the negative values
659 if (!sideA) j+=18; //move to the correct side
661 Float_t deltaPhi = 0.;
662 Float_t phiEdge=sectorWidth*i;
663 if (inphi>phiEdge) {deltaPhi=phiEdge+sectorWidth-inphi;}
664 else if ((outphi>=phiEdge) && (outphi<(phiEdge+sectorWidth))) {deltaPhi=outphi-phiEdge;}
665 else {deltaPhi=sectorWidth;}
667 Float_t v = fVoltageMap[(layer==1)?(j):(j+36)];
668 if (v<=fBadIROCthreshhold) trackLengthInBad+=deltaPhi;
669 if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold) trackLengthInLowGain+=deltaPhi;
672 //for now low gain and bad (off) chambers are treated equally
673 Float_t lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal;
675 //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);
677 if (lengthFractionInBadSectors>fMaxBadLengthFraction)
679 //printf("%%%%%%%% %s kChamberLowGain\n",(layer==1)?"IROC":"OROC");
680 return kChamberLowGain;
683 return kChamberHighGain;
686 //_____________________________________________________________________________
687 Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const
689 //return the radius of the outermost padrow containing a cluster in TPC
691 const TBits* clusterMap=track->GetTPCClusterMapPtr();
692 if (!clusterMap) return 0.;
694 //from AliTPCParam, radius of first IROC row
695 const Float_t rfirstIROC = 8.52249984741210938e+01;
696 const Float_t padrowHeightIROC = 0.75;
697 const Float_t rfirstOROC0 = 1.35100006103515625e+02;
698 const Float_t padrowHeightOROC0 = 1.0;
699 const Float_t rfirstOROC1 = 1.99350006103515625e+02;
700 const Float_t padrowHeightOROC1 = 1.5;
703 while ((--maxPadRow)>0 && !clusterMap->TestBitNumber(maxPadRow)){}
704 if (maxPadRow>126) return rfirstOROC1+(maxPadRow-126-1)*padrowHeightOROC1;
705 if (maxPadRow>62) return rfirstOROC0+(maxPadRow-62-1)*padrowHeightOROC0;
706 if (maxPadRow>0) return rfirstIROC+(maxPadRow-1)*padrowHeightIROC;
710 //_____________________________________________________________________________
711 void AliTPCPIDResponse::InvalidateCurrentValues()
713 //make the "current" stored values from the last processed track invalid
714 fCurrentResponseFunction=NULL;
717 fCurrentGainScenario=kGainScenarioInvalid;
720 //_____________________________________________________________________________
721 Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const
723 //calculate the coordinates of the apex of the track
728 Double_t r = 1./track->OneOverPt()/0.0299792458/magField; //signed - will determine the direction of b
729 //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);
730 //find orthogonal vector (Gram-Schmidt)
731 Double_t alpha = (p[0]*x[0] + p[1]*x[1])/(p[0]*p[0] + p[1]*p[1]);
733 b[0]=x[0]-alpha*p[0];
734 b[1]=x[1]-alpha*p[1];
736 Double_t norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
737 if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
744 //printf("center: x:%.2f, y:%.2f\n",b[0],b[1]);
746 norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
747 if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
749 position[0]=b[0]+b[0]*TMath::Abs(r)/norm;
750 position[1]=b[1]+b[1]*TMath::Abs(r)/norm;