]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliESDpid.cxx
Compilation warning fixed.
[u/mrichter/AliRoot.git] / STEER / AliESDpid.cxx
CommitLineData
8c6a71ab 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
4f679a16 16/* $Id$ */
17
8c6a71ab 18//-----------------------------------------------------------------
19// Implementation of the combined PID class
4f679a16 20// For the Event Summary Data Class
21// produced by the reconstruction process
22// and containing information on the particle identification
8c6a71ab 23// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
24//-----------------------------------------------------------------
25
10d100d4 26#include "AliLog.h"
27#include "AliPID.h"
8c6a71ab 28#include "AliESDpid.h"
af885e0f 29#include "AliESDEvent.h"
8c6a71ab 30#include "AliESDtrack.h"
31
32ClassImp(AliESDpid)
33
10d100d4 34Int_t AliESDpid::MakePID(AliESDEvent *event, Bool_t TPConly, Float_t TimeZeroTOF) const {
35 //
36 // Calculate probabilities for all detectors, except if TPConly==kTRUE
37 // and combine PID
38 //
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
42 //
43
44 /*
45 Float_t TimeZeroTOF = 0;
46 if (subtractT0)
47 TimeZeroTOF = event->GetT0();
48 */
49 Int_t nTrk=event->GetNumberOfTracks();
50 for (Int_t iTrk=0; iTrk<nTrk; iTrk++) {
51 AliESDtrack *track=event->GetTrack(iTrk);
52 MakeTPCPID(track);
53 if (!TPConly) {
54 MakeITSPID(track);
55 MakeTOFPID(track, TimeZeroTOF);
56 //MakeHMPIDPID(track);
ce5d6908 57 //MakeTRDPID(track);
10d100d4 58 }
59 CombinePID(track);
60 }
61 return 0;
62}
8c6a71ab 63//_________________________________________________________________________
10d100d4 64void AliESDpid::MakeTPCPID(AliESDtrack *track) const
8c6a71ab 65{
4f679a16 66 //
10d100d4 67 // TPC pid using bethe-bloch and gaussian response
4f679a16 68 //
10d100d4 69 if ((track->GetStatus()&AliESDtrack::kTPCin )==0)
70 if ((track->GetStatus()&AliESDtrack::kTPCout)==0) return;
8c6a71ab 71
10d100d4 72 Double_t mom = track->GetP();
73 const AliExternalTrackParam *in=track->GetInnerParam();
74 if (in) mom = in->GetP();
8c6a71ab 75
10d100d4 76 Double_t p[AliPID::kSPECIES];
77 Double_t dedx=track->GetTPCsignal();
78 Bool_t mismatch=kTRUE, heavy=kTRUE;
8c6a71ab 79
10d100d4 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;
86 } else {
87 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
88 mismatch=kFALSE;
89 }
8c6a71ab 90
10d100d4 91 // Check for particles heavier than (AliPID::kSPECIES - 1)
92 if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
c630aafd 93
4a78b8c5 94 }
95
10d100d4 96 if (mismatch)
97 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
98
99 track->SetTPCpid(p);
100
101 if (heavy) track->ResetStatus(AliESDtrack::kTPCpid);
102
103}
104//_________________________________________________________________________
105void AliESDpid::MakeITSPID(AliESDtrack *track) const
106{
107 //
108 // ITS PID
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
113 //
114
115 if ((track->GetStatus()&AliESDtrack::kITSin)==0 &&
116 (track->GetStatus()&AliESDtrack::kITSout)==0) return;
117
118 Double_t mom=track->GetP();
119 if (fITSPIDmethod == kITSTruncMean) {
120 Double_t dedx=track->GetITSsignal();
15e979c9 121 Bool_t isSA=kTRUE;
122 Double_t momITS=mom;
123 ULong_t trStatus=track->GetStatus();
124 if(trStatus&AliESDtrack::kTPCin){
125 isSA=kFALSE;
126 const AliExternalTrackParam *in = track->GetInnerParam();
127 momITS=in->P();
128 }
129 UChar_t clumap=track->GetITSClusterMap();
130 Int_t nPointsForPid=0;
131 for(Int_t i=2; i<6; i++){
132 if(clumap&(1<<i)) ++nPointsForPid;
133 }
134
135 if(nPointsForPid<3) { // track not to be used for combined PID purposes
136 track->ResetStatus(AliESDtrack::kITSpid);
137 return;
138 }
139
10d100d4 140 Double_t p[10];
15e979c9 141
10d100d4 142 Bool_t mismatch=kTRUE, heavy=kTRUE;
143 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
144 Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
15e979c9 145 Double_t bethe=fITSResponse.Bethe(momITS,mass);
146 Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
10d100d4 147 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
148 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
149 } else {
150 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
151 mismatch=kFALSE;
152 }
153
154 // Check for particles heavier than (AliPID::kSPECIES - 1)
155 if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
156
4a78b8c5 157 }
c630aafd 158
10d100d4 159 if (mismatch)
160 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
161
162 track->SetITSpid(p);
163
164 if (heavy) track->ResetStatus(AliESDtrack::kITSpid);
165 }
166 else { // Likelihood method
167 Double_t condprobfun[AliPID::kSPECIES];
168 Double_t qclu[4];
169 track->GetITSdEdxSamples(qclu);
170 fITSResponse.GetITSProbabilities(mom,qclu,condprobfun);
171 track->SetITSpid(condprobfun);
8c6a71ab 172 }
7a8614f3 173
10d100d4 174}
175//_________________________________________________________________________
176void AliESDpid::MakeTOFPID(AliESDtrack *track, Float_t TimeZeroTOF) const
177{
178 //
179 // TOF PID using gaussian response
180 //
181 if ((track->GetStatus()&AliESDtrack::kTOFout)==0) return;
182 if ((track->GetStatus()&AliESDtrack::kTIME)==0) return;
183
184 Double_t time[AliPID::kSPECIESN];
185 track->GetIntegratedTimes(time);
186
187 Double_t sigma[AliPID::kSPECIES];
188 for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
189 sigma[iPart] = fTOFResponse.GetExpectedSigma(track->GetP(),time[iPart],AliPID::ParticleMass(iPart));
190 }
191
192 AliDebugGeneral("AliESDpid::MakeTOFPID",2,
193 Form("Expected TOF signals [ps]: %f %f %f %f %f",
194 time[AliPID::kElectron],
195 time[AliPID::kMuon],
196 time[AliPID::kPion],
197 time[AliPID::kKaon],
198 time[AliPID::kProton]));
199
200 AliDebugGeneral("AliESDpid::MakeTOFPID",2,
201 Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
202 sigma[AliPID::kElectron],
203 sigma[AliPID::kMuon],
204 sigma[AliPID::kPion],
205 sigma[AliPID::kKaon],
206 sigma[AliPID::kProton]
207 ));
208
209 Double_t tof = track->GetTOFsignal() - TimeZeroTOF;
210
211 Double_t p[AliPID::kSPECIES];
212 Bool_t mismatch = kTRUE, heavy = kTRUE;
213 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
214 Double_t sig = sigma[j];
215 if (TMath::Abs(tof-time[j]) > fRange*sig) {
216 p[j] = TMath::Exp(-0.5*fRange*fRange)/sig;
217 } else
218 p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
219
220 // Check the mismatching
221 Double_t mass = AliPID::ParticleMass(j);
222 Double_t pm = fTOFResponse.GetMismatchProbability(track->GetP(),mass);
223 if (p[j]>pm) mismatch = kFALSE;
224
225 // Check for particles heavier than (AliPID::kSPECIES - 1)
226 if (tof < (time[j] + fRange*sig)) heavy=kFALSE;
227
228 }
229
230 if (mismatch)
231 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
232
233 track->SetTOFpid(p);
234
235 if (heavy) track->ResetStatus(AliESDtrack::kTOFpid);
236}
237//_________________________________________________________________________
b439f460 238void AliESDpid::MakeTRDPID(AliESDtrack *track) const
239{
240 //
241 // Method to recalculate the TRD PID probabilities
242 //
243 if((track->GetStatus()&AliESDtrack::kTRDout)==0) return;
244 Double_t prob[AliPID::kSPECIES]; Float_t mom[6];
245 Double_t dedx[48]; // Allocate space for the maximum number of TRD slices
246 for(Int_t ilayer = 0; ilayer < 6; ilayer++){
247 mom[ilayer] = track->GetTRDmomentum(ilayer);
248 for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
249 dedx[ilayer*track->GetNumberOfTRDslices()+islice] = track->GetTRDslice(ilayer, islice);
250 }
251 }
252 fTRDResponse.GetResponse(track->GetNumberOfTRDslices(), dedx, mom, prob);
253 track->SetTRDpid(prob);
254}
255//_________________________________________________________________________
10d100d4 256void AliESDpid::CombinePID(AliESDtrack *track) const
257{
258 //
259 // Combine the information of various detectors
260 // to determine the Particle Identification
261 //
262 Int_t ns=AliPID::kSPECIES;
263 Double_t p[10]={1.,1.,1.,1.,1.,1.,1.,1.,1.,1.};
264
265 if (track->IsOn(AliESDtrack::kITSpid)) {
266 Double_t d[10];
267 track->GetITSpid(d);
268 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
269 }
270
271 if (track->IsOn(AliESDtrack::kTPCpid)) {
272 Double_t d[10];
273 track->GetTPCpid(d);
274 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
275 }
276
277 if (track->IsOn(AliESDtrack::kTRDpid)) {
278 Double_t d[10];
279 track->GetTRDpid(d);
280 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
281 }
282
283 if (track->IsOn(AliESDtrack::kTOFpid)) {
284 Double_t d[10];
285 track->GetTOFpid(d);
286 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
287 }
288
289 if (track->IsOn(AliESDtrack::kHMPIDpid)) {
290 Double_t d[10];
291 track->GetHMPIDpid(d);
292 for (Int_t j=0; j<ns; j++) p[j]*=d[j];
293 }
294
295 track->SetESDpid(p);
8c6a71ab 296}