Fixes for bug #52499: Field polarities inconsistiency
[u/mrichter/AliRoot.git] / STEER / AliTOFpidESD.cxx
CommitLineData
b9d9108c 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//-----------------------------------------------------------------//
17// //
18// Implementation of the TOF PID class //
19// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //
20// //
21//-----------------------------------------------------------------//
22
23#include "TMath.h"
24#include "AliLog.h"
25
26#include "AliESDtrack.h"
27#include "AliESDEvent.h"
28
29#include "AliTOFpidESD.h"
30
31ClassImp(AliTOFpidESD)
32
33//_________________________________________________________________________
34AliTOFpidESD::AliTOFpidESD():
35 fSigma(0),
36 fRange(0),
37 fPmax(0), // zero at 0.5 GeV/c for pp
38 fTime0(0)
39{
40}
41//_________________________________________________________________________
42AliTOFpidESD::AliTOFpidESD(Double_t *param):
43 fSigma(param[0]),
44 fRange(param[1]),
45 fPmax(0), // zero at 0.5 GeV/c for pp
46 fTime0(0)
47{
48 //
49 // The main constructor
50 //
51 //
52
53 //fPmax=TMath::Exp(-0.5*3*3)/fSigma; // ~3 sigma at 0.5 GeV/c for PbPb
54}
55
56//_________________________________________________________________________
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 //
63 const Double_t km=0.5; // "reference" momentum (GeV/c)
64
65 Double_t ref2=km*km*km*km/(km*km + mass*mass);// "reference" (p*beta)^2
66 Double_t p2beta2=p*p*p*p/(p*p + mass*mass);
67
68 return fPmax*ref2/p2beta2;
69}
70
71//_________________________________________________________________________
72Int_t AliTOFpidESD::MakePID(AliESDEvent *event, Double_t timeZero)
73{
74 //
75 // This function calculates the "detector response" PID probabilities
76 // Just for a bare hint...
77
78 AliDebug(1,Form("TOF PID Parameters: Sigma (ps)= %f, Range= %f",fSigma,fRange));
79 AliDebug(1,"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ \n");
80
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
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;
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;
139 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
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;
145
146 // Check the mismatching
147 Double_t mass = AliPID::ParticleMass(j);
148 Double_t pm = GetMismatchProbability(t->GetP(),mass);
149 if (p[j]>pm) mismatch = kFALSE;
150
151 // Check for particles heavier than (AliPID::kSPECIES - 1)
152 if (tof < (time[j] + fRange*sig)) heavy=kFALSE;
153
154 }
155
156 if (mismatch)
157 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
158
159 t->SetTOFpid(p);
160
161 if (heavy) t->ResetStatus(AliESDtrack::kTOFpid);
162
163 }
164
165 delete[] tracks;
166
167 return 0;
168}
169
170//_________________________________________________________________________
171Bool_t AliTOFpidESD::ExpectedSignals(const AliESDtrack *t,
172 Double_t s[],Int_t n) const
173{
174 //
175 // Return the expected PID signals for the involved particle species
176 //
177
178 if (n > AliPID::kSPECIESN) return kFALSE;
179 if ( !t->IsOn(AliESDtrack::kTIME) ) return kFALSE;
180
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);
208 }
209 return kTRUE;
210}
211
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 //
220
221 Double_t time[AliPID::kSPECIESN];
222 if ( !ExpectedSignals(t,time,n) ) return kFALSE;
223
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}
231
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 //
240
241 if (Int_t(n) >= AliPID::kSPECIESN) return -1.;
242 if ( !t->IsOn(AliESDtrack::kTIME) ) return -1.;
243
244 Double_t time[AliPID::kSPECIESN];
245 t->GetIntegratedTimes(time);
246
247 return time[n];
248}
249
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 //
259
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);
271}
272
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}