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 **************************************************************************/
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 //-----------------------------------------------------------------
32 #include "AliTOFHeader.h"
33 #include "AliESDpid.h"
34 #include "AliESDEvent.h"
35 #include "AliESDtrack.h"
36 #include "AliMCEvent.h"
38 #include <AliDetectorPID.h>
42 Int_t AliESDpid::MakePID(AliESDEvent *event, Bool_t TPConly, Float_t timeZeroTOF) const {
44 // Calculate probabilities for all detectors, except if TPConly==kTRUE
47 // Option TPConly==kTRUE is used during reconstruction,
48 // because ITS tracking uses TPC pid
49 // HMPID and TRD pid are done in detector reconstructors
53 Float_t timeZeroTOF = 0;
55 timeZeroTOF = event->GetT0();
57 Int_t nTrk=event->GetNumberOfTracks();
58 for (Int_t iTrk=0; iTrk<nTrk; iTrk++) {
59 AliESDtrack *track=event->GetTrack(iTrk);
63 MakeTOFPID(track, timeZeroTOF);
64 //MakeHMPIDPID(track);
71 //_________________________________________________________________________
72 Float_t AliESDpid::GetTPCsignalTunedOnData(const AliVTrack *t) const {
73 AliESDtrack *track = (AliESDtrack *) t;
74 Float_t dedx = track->GetTPCsignalTunedOnData();
76 if(dedx > 0) return dedx;
78 Double_t mom = t->GetTPCmomentum();
79 dedx = t->GetTPCsignal();
80 track->SetTPCsignalTunedOnData(dedx);
81 if(dedx < 20) return dedx;
83 AliPID::EParticleType type = AliPID::kPion;
85 AliMCEventHandler* eventHandler=dynamic_cast<AliMCEventHandler*>(fEventHandler);
87 AliMCEvent* mcEvent = eventHandler->MCEvent();
90 AliMCParticle *MCpart = (AliMCParticle *) mcEvent->GetTrack(TMath::Abs(t->GetLabel()));
91 TParticle *part = MCpart->Particle();
93 Int_t iS = TMath::Abs(part->GetPdgCode());
95 if(iS==AliPID::ParticleCode(AliPID::kElectron)){
96 type = AliPID::kElectron;
98 else if(iS==AliPID::ParticleCode(AliPID::kMuon)){
101 else if(iS==AliPID::ParticleCode(AliPID::kPion)){
102 type = AliPID::kPion;
104 else if(iS==AliPID::ParticleCode(AliPID::kKaon)){
105 type = AliPID::kKaon;
107 else if(iS==AliPID::ParticleCode(AliPID::kProton)){
108 type = AliPID::kProton;
110 else if(iS==AliPID::ParticleCode(AliPID::kDeuteron)){ // d
111 type = AliPID::kDeuteron;
113 else if(iS==AliPID::ParticleCode(AliPID::kTriton)){ // t
114 type = AliPID::kTriton;
116 else if(iS==AliPID::ParticleCode(AliPID::kHe3)){ // 3He
119 else if(iS==AliPID::ParticleCode(AliPID::kAlpha)){ // 4He
120 type = AliPID::kAlpha;
126 Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
127 Double_t sigma=fTPCResponse.GetExpectedSigma(mom,t->GetTPCsignalN(),type);
128 dedx = gRandom->Gaus(bethe,sigma);
129 if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
134 track->SetTPCsignalTunedOnData(dedx);
137 //_________________________________________________________________________
138 void AliESDpid::MakeTPCPID(AliESDtrack *track) const
141 // TPC pid using bethe-bloch and gaussian response
143 if ((track->GetStatus()&AliESDtrack::kTPCin )==0)
144 if ((track->GetStatus()&AliESDtrack::kTPCout)==0) return;
146 Double_t mom = track->GetP();
147 const AliExternalTrackParam *in=track->GetInnerParam();
148 if (in) mom = in->GetP();
150 Double_t p[AliPID::kSPECIES];
151 Double_t dedx=track->GetTPCsignal();
152 Bool_t mismatch=kTRUE, heavy=kTRUE;
154 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
155 AliPID::EParticleType type=AliPID::EParticleType(j);
156 Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
157 Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
158 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
159 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
161 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
165 // Check for particles heavier than (AliPID::kSPECIES - 1)
166 if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
171 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
175 if (heavy) track->ResetStatus(AliESDtrack::kTPCpid);
178 //_________________________________________________________________________
179 void AliESDpid::MakeITSPID(AliESDtrack *track) const
183 // Two options, depending on fITSPIDmethod:
184 // 1) Truncated mean method
185 // 2) Likelihood, using charges measured in all 4 layers and
186 // Landau+gaus response functions
189 if ((track->GetStatus()&AliESDtrack::kITSin)==0 &&
190 (track->GetStatus()&AliESDtrack::kITSout)==0) return;
192 Double_t mom=track->GetP();
193 if (fITSPIDmethod == kITSTruncMean) {
194 Double_t dedx=track->GetITSsignal();
197 ULong_t trStatus=track->GetStatus();
198 if(trStatus&AliESDtrack::kTPCin) isSA=kFALSE;
199 UChar_t clumap=track->GetITSClusterMap();
200 Int_t nPointsForPid=0;
201 for(Int_t i=2; i<6; i++){
202 if(clumap&(1<<i)) ++nPointsForPid;
205 if(nPointsForPid<3) { // track not to be used for combined PID purposes
206 track->ResetStatus(AliESDtrack::kITSpid);
212 Bool_t mismatch=kTRUE, heavy=kTRUE;
213 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
214 Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
215 Double_t bethe=fITSResponse.Bethe(momITS,mass);
216 Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
217 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
218 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
220 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
224 // Check for particles heavier than (AliPID::kSPECIES - 1)
225 if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
230 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
234 if (heavy) track->ResetStatus(AliESDtrack::kITSpid);
236 else { // Likelihood method
237 Double_t condprobfun[AliPID::kSPECIES];
239 track->GetITSdEdxSamples(qclu);
240 fITSResponse.GetITSProbabilities(mom,qclu,condprobfun);
241 track->SetITSpid(condprobfun);
245 //_________________________________________________________________________
246 void AliESDpid::MakeTOFPID(AliESDtrack *track, Float_t /*timeZeroTOF*/) const
249 // TOF PID using gaussian response
252 if ((track->GetStatus()&AliESDtrack::kTOFout)==0) return;
253 if ((track->GetStatus()&AliESDtrack::kTIME)==0) return;
254 if ((track->GetStatus()&AliESDtrack::kITSin)==0) return;
256 Int_t ibin = fTOFResponse.GetMomBin(track->GetP());
257 Float_t timezero = fTOFResponse.GetT0bin(ibin);
259 Double_t time[AliPID::kSPECIES];
260 track->GetIntegratedTimes(time);
262 Double_t sigma[AliPID::kSPECIES];
263 for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
264 sigma[iPart] = fTOFResponse.GetExpectedSigma(track->GetP(),time[iPart],AliPID::ParticleMass(iPart));
267 AliDebugGeneral("AliESDpid::MakeTOFPID",2,
268 Form("Expected TOF signals [ps]: %f %f %f %f %f",
269 time[AliPID::kElectron],
273 time[AliPID::kProton]));
275 AliDebugGeneral("AliESDpid::MakeTOFPID",2,
276 Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
277 sigma[AliPID::kElectron],
278 sigma[AliPID::kMuon],
279 sigma[AliPID::kPion],
280 sigma[AliPID::kKaon],
281 sigma[AliPID::kProton]
284 Double_t tof = track->GetTOFsignal() - timezero;
286 Double_t p[AliPID::kSPECIES];
287 Bool_t mismatch = kTRUE, heavy = kTRUE;
288 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
289 Double_t sig = sigma[j];
290 if (TMath::Abs(tof-time[j]) > (fRange+2)*sig) {
291 p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
293 p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
295 // Check the mismatching
297 Double_t mass = AliPID::ParticleMass(j);
298 Double_t pm = fTOFResponse.GetMismatchProbability(track->GetP(),mass);
299 if (p[j]>pm) mismatch = kFALSE;
301 // Check for particles heavier than (AliPID::kSPECIES - 1)
302 if (tof < (time[j] + fRange*sig)) heavy=kFALSE;
308 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
313 if (heavy) track->ResetStatus(AliESDtrack::kTOFpid);
315 if (!CheckTOFMatching(track)) track->SetStatus(AliESDtrack::kTOFmismatch);
316 else track->ResetStatus(AliESDtrack::kTOFmismatch);
318 //_________________________________________________________________________
319 void AliESDpid::MakeTRDPID(AliESDtrack *track) const
322 // Method to recalculate the TRD PID probabilities
324 Double_t prob[AliPID::kSPECIES];
325 ComputeTRDProbability(track, AliPID::kSPECIES, prob);
326 track->SetTRDpid(prob);
328 //_________________________________________________________________________
329 void AliESDpid::CombinePID(AliESDtrack *track) const
332 // Combine the information of various detectors
333 // to determine the Particle Identification
335 Int_t ns=AliPID::kSPECIES;
336 Double_t p[10]={1.,1.,1.,1.,1.,1.,1.,1.,1.,1.};
338 if (track->IsOn(AliESDtrack::kITSpid)) {
341 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
344 if (track->IsOn(AliESDtrack::kTPCpid)) {
347 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
350 if (track->IsOn(AliESDtrack::kTRDpid)) {
353 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
356 if (track->IsOn(AliESDtrack::kTOFpid)) {
359 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
362 if (track->IsOn(AliESDtrack::kHMPIDpid)) {
364 track->GetHMPIDpid(d);
365 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
370 //_________________________________________________________________________
371 Bool_t AliESDpid::CheckTOFMatching(AliESDtrack *track) const{
373 // Check pid matching of TOF with TPC as reference
375 Bool_t status = kFALSE;
377 Double_t exptimes[5];
378 track->GetIntegratedTimes(exptimes);
380 Float_t p = track->P();
382 Float_t dedx = track->GetTPCsignal();
383 Float_t time = track->GetTOFsignal() - fTOFResponse.GetStartTime(p);
386 track->GetInnerPxPyPz(ptpc);
387 Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]);
389 for(Int_t i=0;i < 5;i++){
390 AliPID::EParticleType type=AliPID::EParticleType(i);
392 Float_t resolutionTOF = fTOFResponse.GetExpectedSigma(p, exptimes[i], AliPID::ParticleMass(i));
393 if(TMath::Abs(exptimes[i] - time) < fRange * resolutionTOF){
394 Float_t dedxExp = fTPCResponse.GetExpectedSignal(momtpc,type);
395 Float_t resolutionTPC = fTPCResponse.GetExpectedSigma(momtpc,track->GetTPCsignalN(),type);
397 if(TMath::Abs(dedx - dedxExp) < fRangeTOFMismatch * resolutionTPC){
404 Float_t resolutionTOFpr = fTOFResponse.GetExpectedSigma(p, exptimes[4], AliPID::ParticleMass(4));
405 if(!status && (exptimes[4] + fRange*resolutionTOFpr < time)) status = kTRUE;
411 Float_t AliESDpid::NumberOfSigmasTOF(const AliVParticle *track, AliPID::EParticleType type, const Float_t /*timeZeroTOF*/) const
414 // Number of sigma implementation for the TOF
417 AliVTrack *vtrack=(AliVTrack*)track;
418 // look for cached value first
419 if (vtrack->GetDetectorPID()){
420 return vtrack->GetDetectorPID()->GetNumberOfSigmas(kTOF, type);
423 Double_t expTime = fTOFResponse.GetExpectedSignal(vtrack,type);
424 return (vtrack->GetTOFsignal() - fTOFResponse.GetStartTime(vtrack->P()) - expTime)/fTOFResponse.GetExpectedSigma(vtrack->P(),expTime,AliPID::ParticleMassZ(type));