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 //-----------------------------------------------------------------
28 #include "AliESDpid.h"
29 #include "AliESDEvent.h"
30 #include "AliESDtrack.h"
34 Int_t AliESDpid::MakePID(AliESDEvent *event, Bool_t TPConly, Float_t TimeZeroTOF) const {
36 // Calculate probabilities for all detectors, except if TPConly==kTRUE
39 // Option TPConly==kTRUE is used during reconstruction,
40 // because ITS tracking uses TPC pid
41 // HMPID and TRD pid are done in detector reconstructors
45 Float_t TimeZeroTOF = 0;
47 TimeZeroTOF = event->GetT0();
49 Int_t nTrk=event->GetNumberOfTracks();
50 for (Int_t iTrk=0; iTrk<nTrk; iTrk++) {
51 AliESDtrack *track=event->GetTrack(iTrk);
55 MakeTOFPID(track, TimeZeroTOF);
56 //MakeHMPIDPID(track);
63 //_________________________________________________________________________
64 void AliESDpid::MakeTPCPID(AliESDtrack *track) const
67 // TPC pid using bethe-bloch and gaussian response
69 if ((track->GetStatus()&AliESDtrack::kTPCin )==0)
70 if ((track->GetStatus()&AliESDtrack::kTPCout)==0) return;
72 Double_t mom = track->GetP();
73 const AliExternalTrackParam *in=track->GetInnerParam();
74 if (in) mom = in->GetP();
76 Double_t p[AliPID::kSPECIES];
77 Double_t dedx=track->GetTPCsignal();
78 Bool_t mismatch=kTRUE, heavy=kTRUE;
80 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
81 AliPID::EParticleType type=AliPID::EParticleType(j);
82 Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
83 Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
84 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
85 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
87 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
91 // Check for particles heavier than (AliPID::kSPECIES - 1)
92 if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
97 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
101 if (heavy) track->ResetStatus(AliESDtrack::kTPCpid);
104 //_________________________________________________________________________
105 void AliESDpid::MakeITSPID(AliESDtrack *track) const
109 // Two options, depending on fITSPIDmethod:
110 // 1) Truncated mean method
111 // 2) Likelihood, using charges measured in all 4 layers and
112 // Landau+gaus response functions
115 if ((track->GetStatus()&AliESDtrack::kITSin)==0 &&
116 (track->GetStatus()&AliESDtrack::kITSout)==0) return;
118 Double_t mom=track->GetP();
119 if (fITSPIDmethod == kITSTruncMean) {
120 Double_t dedx=track->GetITSsignal();
122 Bool_t mismatch=kTRUE, heavy=kTRUE;
123 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
124 Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
125 Double_t bethe=fITSResponse.Bethe(mom,mass);
126 Double_t sigma=fITSResponse.GetResolution(bethe);
127 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
128 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
130 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
134 // Check for particles heavier than (AliPID::kSPECIES - 1)
135 if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
140 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
144 if (heavy) track->ResetStatus(AliESDtrack::kITSpid);
146 else { // Likelihood method
147 Double_t condprobfun[AliPID::kSPECIES];
149 track->GetITSdEdxSamples(qclu);
150 fITSResponse.GetITSProbabilities(mom,qclu,condprobfun);
151 track->SetITSpid(condprobfun);
155 //_________________________________________________________________________
156 void AliESDpid::MakeTOFPID(AliESDtrack *track, Float_t TimeZeroTOF) const
159 // TOF PID using gaussian response
161 if ((track->GetStatus()&AliESDtrack::kTOFout)==0) return;
162 if ((track->GetStatus()&AliESDtrack::kTIME)==0) return;
164 Double_t time[AliPID::kSPECIESN];
165 track->GetIntegratedTimes(time);
167 Double_t sigma[AliPID::kSPECIES];
168 for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
169 sigma[iPart] = fTOFResponse.GetExpectedSigma(track->GetP(),time[iPart],AliPID::ParticleMass(iPart));
172 AliDebugGeneral("AliESDpid::MakeTOFPID",2,
173 Form("Expected TOF signals [ps]: %f %f %f %f %f",
174 time[AliPID::kElectron],
178 time[AliPID::kProton]));
180 AliDebugGeneral("AliESDpid::MakeTOFPID",2,
181 Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
182 sigma[AliPID::kElectron],
183 sigma[AliPID::kMuon],
184 sigma[AliPID::kPion],
185 sigma[AliPID::kKaon],
186 sigma[AliPID::kProton]
189 Double_t tof = track->GetTOFsignal() - TimeZeroTOF;
191 Double_t p[AliPID::kSPECIES];
192 Bool_t mismatch = kTRUE, heavy = kTRUE;
193 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
194 Double_t sig = sigma[j];
195 if (TMath::Abs(tof-time[j]) > fRange*sig) {
196 p[j] = TMath::Exp(-0.5*fRange*fRange)/sig;
198 p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
200 // Check the mismatching
201 Double_t mass = AliPID::ParticleMass(j);
202 Double_t pm = fTOFResponse.GetMismatchProbability(track->GetP(),mass);
203 if (p[j]>pm) mismatch = kFALSE;
205 // Check for particles heavier than (AliPID::kSPECIES - 1)
206 if (tof < (time[j] + fRange*sig)) heavy=kFALSE;
211 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
215 if (heavy) track->ResetStatus(AliESDtrack::kTOFpid);
217 //_________________________________________________________________________
218 void AliESDpid::MakeTRDPID(AliESDtrack *track) const
221 // Method to recalculate the TRD PID probabilities
223 if((track->GetStatus()&AliESDtrack::kTRDout)==0) return;
224 Double_t prob[AliPID::kSPECIES]; Float_t mom[6];
225 Double_t dedx[48]; // Allocate space for the maximum number of TRD slices
226 for(Int_t ilayer = 0; ilayer < 6; ilayer++){
227 mom[ilayer] = track->GetTRDmomentum(ilayer);
228 for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
229 dedx[ilayer*track->GetNumberOfTRDslices()+islice] = track->GetTRDslice(ilayer, islice);
232 fTRDResponse.GetResponse(track->GetNumberOfTRDslices(), dedx, mom, prob);
233 track->SetTRDpid(prob);
235 //_________________________________________________________________________
236 void AliESDpid::CombinePID(AliESDtrack *track) const
239 // Combine the information of various detectors
240 // to determine the Particle Identification
242 Int_t ns=AliPID::kSPECIES;
243 Double_t p[10]={1.,1.,1.,1.,1.,1.,1.,1.,1.,1.};
245 if (track->IsOn(AliESDtrack::kITSpid)) {
248 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
251 if (track->IsOn(AliESDtrack::kTPCpid)) {
254 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
257 if (track->IsOn(AliESDtrack::kTRDpid)) {
260 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
263 if (track->IsOn(AliESDtrack::kTOFpid)) {
266 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
269 if (track->IsOn(AliESDtrack::kHMPIDpid)) {
271 track->GetHMPIDpid(d);
272 for (Int_t j=0; j<ns; j++) p[j]*=d[j];