AliRoot Core  edcc906 (edcc906)
AliESDpid.cxx
Go to the documentation of this file.
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 /* $Id: AliESDpid.cxx 64123 2013-09-05 15:09:53Z morsch $ */
17 
18 //-----------------------------------------------------------------
19 // Implementation of the combined PID class
20 // For the Event Summary Data Class
21 // produced by the reconstruction process
22 // and containing information on the particle identification
23 // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
24 //-----------------------------------------------------------------
25 
26 #include "TArrayI.h"
27 #include "TArrayF.h"
28 
29 #include "AliLog.h"
30 #include "AliPID.h"
31 #include "AliTOFHeader.h"
32 #include "AliESDpid.h"
33 #include "AliESDEvent.h"
34 #include "AliESDtrack.h"
35 #include "AliMCEvent.h"
36 #include "AliTOFPIDParams.h"
37 
38 #include <AliDetectorPID.h>
39 
40 ClassImp(AliESDpid)
41 
42 Bool_t AliESDpid::fgUseElectronExclusionBands = kFALSE;
43 Int_t AliESDpid::fgNSpeciesForTracking = AliPID::kSPECIESC;
44 
45 Int_t AliESDpid::MakePID(AliESDEvent *event, Bool_t TPConly, Float_t /*timeZeroTOF*/) const {
46  //
47  // Calculate probabilities for all detectors, except if TPConly==kTRUE
48  // and combine PID
49  //
50  // Option TPConly==kTRUE is used during reconstruction,
51  // because ITS tracking uses TPC pid
52  // HMPID and TRD pid are done in detector reconstructors
53  //
54 
55  /*
56  Float_t timeZeroTOF = 0;
57  if (subtractT0)
58  timeZeroTOF = event->GetT0();
59  */
60  Int_t nTrk=event->GetNumberOfTracks();
61  for (Int_t iTrk=0; iTrk<nTrk; iTrk++) {
62  AliESDtrack *track=event->GetTrack(iTrk);
63  MakeTPCPID(track);
64  if (!TPConly) {
65  MakeITSPID(track);
66  //MakeTOFPID(track, timeZeroTOF);
67  //MakeHMPIDPID(track);
68  //MakeTRDPID(track);
69  }
70  CombinePID(track);
71  }
72  return 0;
73 }
74 
75 //_________________________________________________________________________
77 {
78  //
79  // TPC pid using bethe-bloch and gaussian response
80  //
81  if ((track->GetStatus()&AliESDtrack::kTPCin )==0)
82  if ((track->GetStatus()&AliESDtrack::kTPCout)==0) return;
83 
84  Double_t mom = track->GetP();
85  const AliExternalTrackParam *in=track->GetInnerParam();
86  if (in) mom = in->GetP();
87 
88  Double_t p[AliPID::kSPECIES];
89  Double_t dedx=track->GetTPCsignal();
90  Bool_t mismatch=kTRUE, heavy=kTRUE;
91 
92  for (Int_t j=0; j<AliPID::kSPECIES; j++) {
94  Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
95  Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
96  if (TMath::Abs(dedx-bethe) > fRange*sigma) {
97  p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
98  } else {
99  p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
100  mismatch=kFALSE;
101  }
102 
103  // Check for particles heavier than (AliPID::kSPECIES - 1)
104  if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
105 
106  }
107 
108  if (mismatch)
109  for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
110 
111  track->SetTPCpid(p);
112 
113  if (heavy) track->ResetStatus(AliESDtrack::kTPCpid);
114 
115 }
116 //_________________________________________________________________________
118 {
119  //
120  // ITS PID
121  // Two options, depending on fITSPIDmethod:
122  // 1) Truncated mean method
123  // 2) Likelihood, using charges measured in all 4 layers and
124  // Landau+gaus response functions
125  //
126 
127  if ((track->GetStatus()&AliESDtrack::kITSin)==0 &&
128  (track->GetStatus()&AliESDtrack::kITSout)==0) return;
129 
130  Double_t mom=track->GetP();
131  if (fITSPIDmethod == kITSTruncMean) {
132  Double_t dedx=track->GetITSsignal();
133  Bool_t isSA=kTRUE;
134  Double_t momITS=mom;
135  ULong_t trStatus=track->GetStatus();
136  if(trStatus&AliESDtrack::kTPCin) isSA=kFALSE;
137  UChar_t clumap=track->GetITSClusterMap();
138  Int_t nPointsForPid=0;
139  for(Int_t i=2; i<6; i++){
140  if(clumap&(1<<i)) ++nPointsForPid;
141  }
142 
143  if(nPointsForPid<3) { // track not to be used for combined PID purposes
145  return;
146  }
147 
148  Double_t p[AliPID::kSPECIES];
149 
150  Bool_t mismatch=kTRUE, heavy=kTRUE;
151  for (Int_t j=0; j<AliPID::kSPECIES; j++) {
152  Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
153  Double_t bethe=fITSResponse.Bethe(momITS,mass);
154  Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
155  if (TMath::Abs(dedx-bethe) > fRange*sigma) {
156  p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
157  } else {
158  p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
159  mismatch=kFALSE;
160  }
161 
162  // Check for particles heavier than (AliPID::kSPECIES - 1)
163  if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
164 
165  }
166 
167  if (mismatch)
168  for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
169 
170  track->SetITSpid(p);
171 
172  if (heavy) track->ResetStatus(AliESDtrack::kITSpid);
173  }
174  else { // Likelihood method
175  Double_t condprobfun[AliPID::kSPECIES];
176  Double_t qclu[4];
177  track->GetITSdEdxSamples(qclu);
178  fITSResponse.GetITSProbabilities(mom,qclu,condprobfun);
179  track->SetITSpid(condprobfun);
180  }
181 
182 }
183 //_________________________________________________________________________
184 void AliESDpid::MakeTOFPID(AliESDtrack *track, Float_t /*timeZeroTOF*/) const
185 {
186  //
187  // TOF PID using gaussian response
188  //
189 
190  if ((track->GetStatus()&AliESDtrack::kTOFout)==0) return;
191  if ((track->GetStatus()&AliESDtrack::kTIME)==0) return;
192  if ((track->GetStatus()&AliESDtrack::kITSin)==0) return;
193 
194  Int_t ibin = fTOFResponse.GetMomBin(track->GetP());
195  Float_t timezero = fTOFResponse.GetT0bin(ibin);
196 
197  Double_t time[AliPID::kSPECIESC];
198  track->GetIntegratedTimes(time);
199 
200  Double_t sigma[AliPID::kSPECIES];
201  for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
202  sigma[iPart] = fTOFResponse.GetExpectedSigma(track->GetP(),time[iPart],AliPID::ParticleMass(iPart));
203  }
204 
205  AliDebugGeneral("AliESDpid::MakeTOFPID",2,
206  Form("Expected TOF signals [ps]: %f %f %f %f %f",
207  time[AliPID::kElectron],
208  time[AliPID::kMuon],
209  time[AliPID::kPion],
210  time[AliPID::kKaon],
211  time[AliPID::kProton]));
212 
213  AliDebugGeneral("AliESDpid::MakeTOFPID",2,
214  Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
215  sigma[AliPID::kElectron],
216  sigma[AliPID::kMuon],
217  sigma[AliPID::kPion],
218  sigma[AliPID::kKaon],
219  sigma[AliPID::kProton]
220  ));
221 
222  Double_t tof = track->GetTOFsignal() - timezero;
223 
224  Double_t p[AliPID::kSPECIES];
225 // Bool_t mismatch = kTRUE;
226  Bool_t heavy = kTRUE;
227  for (Int_t j=0; j<AliPID::kSPECIES; j++) {
228  Double_t sig = sigma[j];
229  if (TMath::Abs(tof-time[j]) > (fRange+2)*sig) {
230  p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
231  } else
232  p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
233 
234  // Check the mismatching
235 
236 // Double_t mass = AliPID::ParticleMass(j);
237 // Double_t pm = fTOFResponse.GetMismatchProbability(track->GetP(),mass);
238 // if (p[j]>pm) mismatch = kFALSE;
239 
240  // Check for particles heavier than (AliPID::kSPECIES - 1)
241  if (tof < (time[j] + fRange*sig)) heavy=kFALSE;
242 
243  }
244 
245  /*
246  if (mismatch)
247  for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
248  */
249 
250  track->SetTOFpid(p);
251 
252  if (heavy) track->ResetStatus(AliESDtrack::kTOFpid);
253 
254  // kTOFmismatch flas is not set because deprecated from 18/02/2013
255  // if (!CheckTOFMatching(track)) track->SetStatus(AliESDtrack::kTOFmismatch);
256  // else track->ResetStatus(AliESDtrack::kTOFmismatch);
257 }
258 //_________________________________________________________________________
260 {
261  //
262  // Method to recalculate the TRD PID probabilities
263  //
264  Double_t prob[AliPID::kSPECIES];
266  track->SetTRDpid(prob);
267 }
268 //_________________________________________________________________________
270 {
271  //
272  // Combine the information of various detectors
273  // to determine the Particle Identification
274  //
275  Double_t p[AliPID::kSPECIES]={1.};
276 
277  if (track->IsOn(AliESDtrack::kITSpid)) {
278  Double_t d[AliPID::kSPECIES];
279  track->GetITSpid(d);
280  for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]*=d[j];
281  }
282 
283  if (track->IsOn(AliESDtrack::kTPCpid)) {
284  Double_t d[AliPID::kSPECIES];
285  track->GetTPCpid(d);
286  for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]*=d[j];
287  }
288 
289  if (track->IsOn(AliESDtrack::kTRDpid)) {
290  Double_t d[AliPID::kSPECIES];
291  track->GetTRDpid(d);
292  for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]*=d[j];
293  }
294 
295  if (track->IsOn(AliESDtrack::kTOFpid)) {
296  Double_t d[AliPID::kSPECIES];
297  track->GetTOFpid(d);
298  for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]*=d[j];
299  }
300 
301  if (track->IsOn(AliESDtrack::kHMPIDpid)) {
302  Double_t d[AliPID::kSPECIES];
303  track->GetHMPIDpid(d);
304  for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]*=d[j];
305  }
306 
307  track->SetESDpid(p);
308 }
309 //_________________________________________________________________________
311  //
312  // Check pid matching of TOF with TPC as reference
313  //
314  Bool_t status = kFALSE;
315 
316  Double_t exptimes[AliPID::kSPECIESC];
317  track->GetIntegratedTimes(exptimes);
318 
319  Float_t p = track->P();
320 
321  Float_t dedx = track->GetTPCsignal();
322  Float_t time = track->GetTOFsignal() - fTOFResponse.GetStartTime(p);
323 
324  Double_t ptpc[3];
325  track->GetInnerPxPyPz(ptpc);
326  Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]);
327 
328  for(Int_t i=0;i < AliPID::kSPECIES;i++){
330 
331  Float_t resolutionTOF = fTOFResponse.GetExpectedSigma(p, exptimes[i], AliPID::ParticleMass(i));
332  if(TMath::Abs(exptimes[i] - time) < fRange * resolutionTOF){
333  Float_t dedxExp = fTPCResponse.GetExpectedSignal(momtpc,type);
334  Float_t resolutionTPC = fTPCResponse.GetExpectedSigma(momtpc,track->GetTPCsignalN(),type);
335 
336  if(TMath::Abs(dedx - dedxExp) < fRangeTOFMismatch * resolutionTPC){
337  status = kTRUE;
338  }
339  }
340  }
341 
342  // for nuclei
343  Float_t resolutionTOFpr = fTOFResponse.GetExpectedSigma(p, exptimes[4], AliPID::ParticleMass(4));
344  if(!status && (exptimes[4] + fRange*resolutionTOFpr < time)) status = kTRUE;
345 
346 
347  return status;
348 }
349 
350 //_________________________________________________________________________
351 Float_t AliESDpid::GetSignalDeltaTOFold(const AliVParticle *track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/) const
352 {
353  //
354  // TOF signal - expected
355  //
356  AliVTrack *vtrack=(AliVTrack*)track;
357 
358  const Double_t expTime = fTOFResponse.GetExpectedSignal(vtrack,type);
359  Double_t tofTime = 99999;
360  if (fTuneMConData && ((fTuneMConDataMask & kDetTOF) == kDetTOF) ) tofTime = (Double_t)this->GetTOFsignalTunedOnData((AliVTrack*)vtrack);
361  else tofTime=vtrack->GetTOFsignal();
362  tofTime = tofTime - fTOFResponse.GetStartTime(vtrack->P());
363  Double_t delta=-9999.;
364 
365  if (!ratio) delta=tofTime-expTime;
366  else if (expTime>1.e-20) delta=tofTime/expTime;
367 
368  return delta;
369 }
370 
371 //_________________________________________________________________________
373 {
374  //
375  // Number of sigma implementation for the TOF
376  //
377 
378  AliVTrack *vtrack=(AliVTrack*)track;
379  Double_t tofTime = 99999;
380  if (fTuneMConData && ((fTuneMConDataMask & kDetTOF) == kDetTOF) ) tofTime = (Double_t)this->GetTOFsignalTunedOnData((AliVTrack*)vtrack);
381  else tofTime=vtrack->GetTOFsignal();
382 
383  Double_t expTime = fTOFResponse.GetExpectedSignal(vtrack,type);
384  return (tofTime - fTOFResponse.GetStartTime(vtrack->P()) - expTime)/fTOFResponse.GetExpectedSigma(vtrack->P(),expTime,AliPID::ParticleMassZ(type));
385 }
386 
387 //_________________________________________________________________________
389 {
390  // assign mass for tracking
391  //
392 
393  // in principle AliPIDCombined could be used to also set priors
394  //AliPIDCombined pidProb;
395  //pidProb.SetDetectorMask(kDetTPC); // use only TPC information, couls also be changed
396  //pidProb.SetSelectedSpecies(AliPID::kSPECIESC); // number of species to use
397  //pidProb.SetDefaultTPCPriors(); // load default priors
398 
399  Double_t prob[AliPID::kSPECIESC]={0.};
400  // EDetPidStatus pidStatus=ComputePIDProbability(kTPC, esdtr, AliPID::kSPECIESC, prob);
402  // check if a valid signal was found, otherwise return pion mass
403  if (pidStatus==kDetNoSignal) { //kDetPidOk) {
405  return;
406  }
407 
408  // or with AliPIDCombined
409  // pidProb.ComputeProbabilities(esdtr, this, p);
410 
411  // find max probability
412  Float_t max=0.,min=1.e9;
413  Int_t pid=-1;
414  // for (Int_t i=0; i<AliPID::kSPECIESC; ++i) {
415  for (Int_t i=0; i<fgNSpeciesForTracking; ++i) { // ignore nuclei
416  if (prob[i]>max) {pid=i; max=prob[i];}
417  if (prob[i]<min) min=prob[i];
418  }
419 
420  //int pid = AliPID::kPion; // this should be substituted by real most probable TPC pid (e,mu -> pion) or poin if no PID possible
421 
422  //
423  if (pid>AliPID::kSPECIESC-1 || (min>=max)) pid = AliPID::kPion;
424  //
425  if (pid==0 && fgUseElectronExclusionBands) { // dE/dx "crossing points" in the TPC
426  Double_t p = esdtr->GetP();
427  if ((p>0.38)&&(p<0.48)) {
428  if (prob[0]<prob[3]*10.) pid = AliPID::kKaon;
429  }
430  else if ((p>0.75)&&(p<0.85)) {
431  if (prob[0]<prob[4]*10.) pid = AliPID::kProton;
432  }
433  }
434 
435  esdtr->SetPIDForTracking( pid );
436  //
437 }
438 
439 
440 //_________________________________________________________________________
442 {
443  // assign masses using for tracking
444  Int_t nTrk=event->GetNumberOfTracks();
445  for (Int_t iTrk=nTrk; iTrk--;) {
446  AliESDtrack *track = event->GetTrack(iTrk);
447  SetPIDForTracking(track);
448  }
449 }
450 
451 //_________________________________________________________________________
453 {
454  // disable electron tag in e/K and e/p crossing (a la Run1)
456  AliInfoClassF("Exclude electron tag in e/K and e/p crossings: %d",fgUseElectronExclusionBands);
457 }
458 
459 //_________________________________________________________________________
461 {
462  // set max number of species to consider for tracking
464  AliInfoClassF("First %d species will be considered for tracking PID",fgNSpeciesForTracking);
465 }
Double_t Bethe(Double_t p, Double_t mass, Bool_t isSA=kFALSE) const
AliTPCcalibPID * pid
Definition: CalibPID.C:69
Double_t GetExpectedSignal(const AliVTrack *track, AliPID::EParticleType type) const
const AliExternalTrackParam * GetInnerParam() const
Definition: AliESDtrack.h:132
virtual Float_t GetNumberOfSigmasTOFold(const AliVParticle *track, AliPID::EParticleType type) const
Definition: AliESDpid.cxx:372
ULong_t GetStatus() const
Definition: AliESDtrack.h:82
Float_t GetT0bin(Int_t ibin) const
Bool_t IsOn(Int_t mask) const
Definition: AliESDtrack.h:81
virtual Float_t GetSignalDeltaTOFold(const AliVParticle *track, AliPID::EParticleType type, Bool_t ratio=kFALSE) const
Definition: AliESDpid.cxx:351
void GetTOFpid(Double_t *p) const
void MakeTPCPID(AliESDtrack *track) const
Definition: AliESDpid.cxx:76
void SetPIDForTracking(Int_t pid)
Definition: AliESDtrack.h:100
void CombinePID(AliESDtrack *track) const
Definition: AliESDpid.cxx:269
virtual Double_t GetTOFsignal() const
Definition: AliVTrack.h:170
static Float_t ParticleMass(Int_t iType)
Definition: AliPID.h:56
virtual Double_t P() const =0
EParticleType
Definition: AliPID.h:27
void SetTRDpid(const Double_t *p)
Double_t GetResolution(Double_t bethe, Int_t nPtsForPid=4, Bool_t isSA=kFALSE, Double_t p=0., AliPID::EParticleType type=AliPID::kPion) const
Float_t p[]
Definition: kNNTest.C:133
void ResetStatus(ULong_t flags)
Definition: AliESDtrack.h:73
AliTPCfastTrack * track
void SetTPCpid(const Double_t *p)
Double_t GetExpectedSignal(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource=kdEdxDefault, Bool_t correctEta=kFALSE, Bool_t correctMultiplicity=kFALSE) const
EDetPidStatus GetComputeTRDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[], AliTRDPIDResponse::ETRDPIDMethod PIDmethod=AliTRDPIDResponse::kLQ1D) const
Float_t fRangeTOFMismatch
Definition: AliESDpid.h:67
void MakePIDForTracking(AliESDEvent *event) const
Definition: AliESDpid.cxx:441
static void SetUseElectronExclusionBands(Bool_t val)
Definition: AliESDpid.cxx:452
void SetTOFpid(const Double_t *p)
ITSPIDmethod fITSPIDmethod
void GetITSdEdxSamples(Double_t s[4]) const
virtual Float_t GetTOFsignalTunedOnData(const AliVTrack *t) const
Double_t GetExpectedSigma(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource=kdEdxDefault, Bool_t correctEta=kFALSE, Bool_t correctMultiplicity=kFALSE) const
EDetPidStatus ComputePIDProbability(EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
void SetESDpid(const Double_t *p)
static Int_t fgNSpeciesForTracking
Definition: AliESDpid.h:71
void GetTPCpid(Double_t *p) const
AliTOFPIDResponse fTOFResponse
Double_t GetTPCsignal() const
Definition: AliESDtrack.h:259
Bool_t GetInnerPxPyPz(Double_t *p) const
Definition: AliESDtrack.h:128
Definition: AliPID.h:19
Double_t GetExpectedSigma(Float_t mom, Float_t tof, Float_t massZ) const
void GetITSpid(Double_t *p) const
void GetITSProbabilities(Float_t mom, Double_t qclu[4], Double_t condprobfun[AliPID::kSPECIES], Bool_t isMC=kFALSE) const
UShort_t GetTPCsignalN() const
Definition: AliESDtrack.h:262
AliITSPIDResponse fITSResponse
static void SetNSpeciesForTracking(Int_t n)
Definition: AliESDpid.cxx:460
void SetPIDForTracking(AliESDtrack *track) const
Definition: AliESDpid.cxx:388
Int_t GetMomBin(Float_t p) const
void GetTRDpid(Double_t *p) const
#define AliInfoClassF(message,...)
Definition: AliLog.h:504
Float_t GetStartTime(Float_t mom) const
Double_t GetITSsignal() const
Definition: AliESDtrack.h:183
void MakeTOFPID(AliESDtrack *track, Float_t) const
Definition: AliESDpid.cxx:184
void GetHMPIDpid(Double_t *p) const
void GetIntegratedTimes(Double_t *times, Int_t nspec=AliPID::kSPECIES) const
Bool_t CheckTOFMatching(AliESDtrack *track) const
Definition: AliESDpid.cxx:310
UChar_t GetITSClusterMap() const
Definition: AliESDtrack.h:197
static Bool_t fgUseElectronExclusionBands
MC event handler.
Definition: AliESDpid.h:70
void SetITSpid(const Double_t *p)
AliTPCPIDResponse fTPCResponse
static Float_t ParticleMassZ(Int_t iType)
Definition: AliPID.h:60
#define AliDebugGeneral(scope, logLevel, message)
Definition: AliLog.h:327
void MakeITSPID(AliESDtrack *track) const
Definition: AliESDpid.cxx:117
Double_t GetTOFsignal() const
void MakeTRDPID(AliESDtrack *track) const
Definition: AliESDpid.cxx:259