Corrected protection
[u/mrichter/AliRoot.git] / TOF / AliTOFpidESD.cxx
CommitLineData
c630aafd 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
0e46b9ae 16//-----------------------------------------------------------------//
17// //
18// Implementation of the TOF PID class //
19// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //
20// //
21//-----------------------------------------------------------------//
22
23#include "TMath.h"
02cdc6ac 24#include "AliLog.h"
0e46b9ae 25
c630aafd 26#include "AliESDtrack.h"
af885e0f 27#include "AliESDEvent.h"
c630aafd 28
0e46b9ae 29#include "AliTOFpidESD.h"
30
c630aafd 31ClassImp(AliTOFpidESD)
32
0f4a7374 33//_________________________________________________________________________
7a8614f3 34AliTOFpidESD::AliTOFpidESD():
655e379f 35 fSigma(0),
7a8614f3 36 fRange(0),
cc79c586 37 fPmax(0), // zero at 0.5 GeV/c for pp
38 fTime0(0)
655e379f 39{
40}
41//_________________________________________________________________________
42AliTOFpidESD::AliTOFpidESD(Double_t *param):
7a8614f3 43 fSigma(param[0]),
44 fRange(param[1]),
cc79c586 45 fPmax(0), // zero at 0.5 GeV/c for pp
46 fTime0(0)
7a8614f3 47{
0f4a7374 48 //
49 // The main constructor
50 //
7a8614f3 51 //
52
53 //fPmax=TMath::Exp(-0.5*3*3)/fSigma; // ~3 sigma at 0.5 GeV/c for PbPb
54}
55
cc79c586 56//_________________________________________________________________________
7a8614f3 57Double_t
58AliTOFpidESD::GetMismatchProbability(Double_t p, Double_t mass) const {
59 //
60 // Returns the probability of mismatching
61 // assuming 1/(p*beta)^2 scaling
62 //
62bdc654 63 const Double_t km=0.5; // "reference" momentum (GeV/c)
7a8614f3 64
62bdc654 65 Double_t ref2=km*km*km*km/(km*km + mass*mass);// "reference" (p*beta)^2
7a8614f3 66 Double_t p2beta2=p*p*p*p/(p*p + mass*mass);
0f4a7374 67
7a8614f3 68 return fPmax*ref2/p2beta2;
0f4a7374 69}
70
0f4a7374 71//_________________________________________________________________________
af885e0f 72Int_t AliTOFpidESD::MakePID(AliESDEvent *event, Double_t timeZero)
bc9f08da 73{
74 //
75 // This function calculates the "detector response" PID probabilities
76 // Just for a bare hint...
77
02cdc6ac 78 AliDebug(1,Form("TOF PID Parameters: Sigma (ps)= %f, Range= %f",fSigma,fRange));
79 AliDebug(1,"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ \n");
7a8614f3 80
cc79c586 81 fTime0=timeZero;
82 return MakePID(event);
83}
84
85//_________________________________________________________________________
86Int_t AliTOFpidESD::MakePID(AliESDEvent *event)
87{
88 //
89 // This function calculates the "detector response" PID probabilities
90 // Just for a bare hint...
91
bc9f08da 92 Int_t ntrk=event->GetNumberOfTracks();
93 AliESDtrack **tracks=new AliESDtrack*[ntrk];
94
95 Int_t i;
96 for (i=0; i<ntrk; i++) {
97 AliESDtrack *t=event->GetTrack(i);
98 tracks[i]=t;
99 }
100
101 for (i=0; i<ntrk; i++) {
102 AliESDtrack *t=tracks[i];
103 if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
104 if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
cc79c586 105
106 Double_t time[AliPID::kSPECIES];
107 if (!ExpectedSignals(t,time,AliPID::kSPECIES)) continue;
108 Double_t sigma[AliPID::kSPECIES];
109 if (!ExpectedSigmas(t,sigma,AliPID::kSPECIES)) continue;
110
111 AliDebug(2,Form("Expected TOF signals [ps]: %f %f %f %f %f",
112 GetExpectedSignal(t,AliPID::kElectron),
113 GetExpectedSignal(t,AliPID::kMuon),
114 GetExpectedSignal(t,AliPID::kPion),
115 GetExpectedSignal(t,AliPID::kKaon),
116 GetExpectedSignal(t,AliPID::kProton)
117 ));
118
119 AliDebug(2,Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
120 GetExpectedSigma(t,AliPID::kElectron),
121 GetExpectedSigma(t,AliPID::kMuon),
122 GetExpectedSigma(t,AliPID::kPion),
123 GetExpectedSigma(t,AliPID::kKaon),
124 GetExpectedSigma(t,AliPID::kProton)
125 ));
126
127 AliDebug(2,Form("Expected TOF std deviations [number of expected sigmas]: %f %f %f %f %f",
128 GetNumberOfSigmas(t,AliPID::kElectron),
129 GetNumberOfSigmas(t,AliPID::kMuon),
130 GetNumberOfSigmas(t,AliPID::kPion),
131 GetNumberOfSigmas(t,AliPID::kKaon),
132 GetNumberOfSigmas(t,AliPID::kProton)
133 ));
134
135 Double_t tof = t->GetTOFsignal() - fTime0;
136
137 Double_t p[AliPID::kSPECIES];
138 Bool_t mismatch = kTRUE, heavy = kTRUE;
bc9f08da 139 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
cc79c586 140 Double_t sig = sigma[j];
141 if (TMath::Abs(tof-time[j]) > fRange*sig) {
142 p[j] = TMath::Exp(-0.5*fRange*fRange)/sig;
143 } else
144 p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
7a8614f3 145
146 // Check the mismatching
cc79c586 147 Double_t mass = AliPID::ParticleMass(j);
148 Double_t pm = GetMismatchProbability(t->GetP(),mass);
149 if (p[j]>pm) mismatch = kFALSE;
7a8614f3 150
151 // Check for particles heavier than (AliPID::kSPECIES - 1)
cc79c586 152 if (tof < (time[j] + fRange*sig)) heavy=kFALSE;
7a8614f3 153
bc9f08da 154 }
7a8614f3 155
156 if (mismatch)
157 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
158
bc9f08da 159 t->SetTOFpid(p);
7a8614f3 160
161 if (heavy) t->ResetStatus(AliESDtrack::kTOFpid);
162
bc9f08da 163 }
164
165 delete[] tracks;
166
167 return 0;
168}
169
170//_________________________________________________________________________
cc79c586 171Bool_t AliTOFpidESD::ExpectedSignals(const AliESDtrack *t,
172 Double_t s[],Int_t n) const
c630aafd 173{
174 //
cc79c586 175 // Return the expected PID signals for the involved particle species
176 //
c630aafd 177
cc79c586 178 if (n > AliPID::kSPECIESN) return kFALSE;
179 if ( !t->IsOn(AliESDtrack::kTIME) ) return kFALSE;
3f83f224 180
cc79c586 181 Double_t time[AliPID::kSPECIESN];
182 t->GetIntegratedTimes(time);
183 for (Int_t i=0; i<n; i++) s[i]=time[i];
184 return kTRUE;
185}
186
187//_________________________________________________________________________
188Bool_t AliTOFpidESD::ExpectedSigmas(const AliESDtrack *t,
189 Double_t s[],Int_t n) const
190{
191 //
192 // Return the expected sigma of PID signals for the involved
193 // particle species.
194 // This approximate (but reasonable) formula takes into account the
195 // relative momentum resolution.
196 //
197
198 Double_t time[AliPID::kSPECIESN];
199 if ( !ExpectedSignals(t,time,n) ) return kFALSE;
200
201 Double_t mom = t->GetP();
202 Double_t dpp = 0.01; //mean relative pt resolution;
203 if (mom>0.5) dpp = 0.01*mom;
204 for (Int_t i=0; i<n; i++) {
205 Double_t mass = AliPID::ParticleMass(i);
206 Double_t sigma = dpp*time[i]/(1.+ mom*mom/(mass*mass));
207 s[i] = TMath::Sqrt(sigma*sigma + fSigma*fSigma);
3f83f224 208 }
cc79c586 209 return kTRUE;
210}
3f83f224 211
cc79c586 212//_________________________________________________________________________
213Bool_t AliTOFpidESD::NumberOfSigmas(const AliESDtrack *t,
214 Double_t s[],Int_t n) const
215{
216 //
217 // Returns the deviation of the actual PID signal from the expected
218 // signal, in units of expected sigmas.
219 //
7a8614f3 220
cc79c586 221 Double_t time[AliPID::kSPECIESN];
222 if ( !ExpectedSignals(t,time,n) ) return kFALSE;
7a8614f3 223
cc79c586 224 if ( !ExpectedSigmas(t,s,n) ) return kFALSE;
225
226 Double_t tof = t->GetTOFsignal() - fTime0;
227 for (Int_t i=0; i<n; i++) s[i] = (time[i]-tof)/s[i];
228
229 return kTRUE;
230}
7a8614f3 231
cc79c586 232//_________________________________________________________________________
233 Double_t AliTOFpidESD::GetExpectedSignal(const AliESDtrack *t,
234 AliPID::EParticleType n) const
235{
236 //
237 // Return the expected PID signal for the specified particle type.
238 // If the operation is not possible, return a negative value.
239 //
7a8614f3 240
cc79c586 241 if (Int_t(n) >= AliPID::kSPECIESN) return -1.;
242 if ( !t->IsOn(AliESDtrack::kTIME) ) return -1.;
7a8614f3 243
cc79c586 244 Double_t time[AliPID::kSPECIESN];
245 t->GetIntegratedTimes(time);
7a8614f3 246
cc79c586 247 return time[n];
248}
7a8614f3 249
cc79c586 250//_________________________________________________________________________
251 Double_t AliTOFpidESD::GetExpectedSigma(const AliESDtrack *t,
252 AliPID::EParticleType n) const
253{
254 //
255 // Return the expected sigma of the PID signal for the specified
256 // particle type.
257 // If the operation is not possible, return a negative value.
258 //
c630aafd 259
cc79c586 260 Double_t time[AliPID::kSPECIESN];
261 if ( !ExpectedSignals(t,time,AliPID::kSPECIESN) ) return -1.;
262
263 Double_t mom = t->GetP();
264 Double_t dpp = 0.01; //mean relative pt resolution;
265 if (mom>0.5) dpp = 0.01*mom;
266
267 Double_t mass = AliPID::ParticleMass(n);
268 Double_t sigma = dpp*time[n]/(1.+ mom*mom/(mass*mass));
269
270 return TMath::Sqrt(sigma*sigma + fSigma*fSigma);
c630aafd 271}
0f4a7374 272
cc79c586 273//_________________________________________________________________________
274 Double_t AliTOFpidESD::GetNumberOfSigmas(const AliESDtrack *t,
275 AliPID::EParticleType n) const
276{
277 //
278 // Returns the deviation of the actual PID signal from the expected
279 // signal for the specified particle type, in units of expected
280 // sigmas.
281 // If the operation is not possible, return a negative value.
282 //
283
284 Double_t time=GetExpectedSignal(t,n);;
285 if (time < 0.) return -1.;
286
287 Double_t sigma=GetExpectedSigma(t,n);
288 if (sigma < 0.) return -1;
289
290 Double_t tof=t->GetTOFsignal() - fTime0;
291 return (time-tof)/sigma;
292}