]>
Commit | Line | Data |
---|---|---|
46d29e70 | 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 | * * | |
43bfc8af | 7 | * Permission to use , copy, modify and distribute this software and its * |
46d29e70 | 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 | ||
0fa7dfa7 | 16 | /* $Id$ */ |
46d29e70 | 17 | |
a2cb5b3d | 18 | #include <Riostream.h> |
fe33d239 | 19 | |
53c17fbf | 20 | #include <TMath.h> |
e1e6896f | 21 | #include <TVector2.h> |
46d29e70 | 22 | |
6c94f330 | 23 | #include "AliTracker.h" |
53c17fbf | 24 | #include "AliESDtrack.h" |
fe33d239 | 25 | |
46d29e70 | 26 | #include "AliTRDgeometry.h" |
27 | #include "AliTRDcluster.h" | |
28 | #include "AliTRDtrack.h" | |
53c17fbf | 29 | #include "AliTRDtracklet.h" |
b3a5a838 | 30 | |
43bfc8af | 31 | // A. Bercuci - used for PID calculations |
32 | #include "AliTRDcalibDB.h" | |
720a0a16 | 33 | #include "Cal/AliTRDCalPID.h" |
43bfc8af | 34 | |
46d29e70 | 35 | ClassImp(AliTRDtrack) |
36 | ||
53c17fbf | 37 | /////////////////////////////////////////////////////////////////////////////// |
38 | // // | |
39 | // Represents a reconstructed TRD track // | |
40 | // Local TRD Kalman track // | |
41 | // // | |
42 | /////////////////////////////////////////////////////////////////////////////// | |
7ad19338 | 43 | |
fe33d239 | 44 | //_____________________________________________________________________________ |
45 | AliTRDtrack::AliTRDtrack() | |
46 | :AliKalmanTrack() | |
47 | ,fSeedLab(-1) | |
48 | ,fdEdx(0) | |
49 | ,fDE(0) | |
79143cbf | 50 | ,fClusterOwner(kFALSE) |
fe33d239 | 51 | ,fStopped(kFALSE) |
52 | ,fLhElectron(0) | |
53 | ,fNWrong(0) | |
54 | ,fNRotate(0) | |
55 | ,fNCross(0) | |
56 | ,fNExpected(0) | |
57 | ,fNLast(0) | |
58 | ,fNExpectedLast(0) | |
59 | ,fNdedx(0) | |
60 | ,fChi2Last(1e10) | |
61 | ,fBackupTrack(0x0) | |
f146da82 | 62 | { |
fe33d239 | 63 | // |
64 | // AliTRDtrack default constructor | |
65 | // | |
66 | ||
67 | for (Int_t i = 0; i < kNplane; i++) { | |
68 | for (Int_t j = 0; j < kNslice; j++) { | |
69 | fdEdxPlane[i][j] = 0.0; | |
6d45eaef | 70 | } |
f146da82 | 71 | fTimBinPlane[i] = -1; |
af26ce80 | 72 | fMom[i] = -1.; |
73 | fSnp[i] = 0.; | |
74 | fTgl[i] = 0.; | |
75 | } | |
fe33d239 | 76 | |
77 | for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) { | |
78 | fIndex[i] = 0; | |
6c94f330 | 79 | fIndexBackup[i] = 0; |
fe33d239 | 80 | fdQdl[i] = 0; |
af26ce80 | 81 | fClusters[i] = 0x0; |
82 | } | |
fe33d239 | 83 | |
84 | for (Int_t i = 0; i < 3; i++) { | |
85 | fBudget[i] = 0; | |
f146da82 | 86 | } |
fe33d239 | 87 | |
f146da82 | 88 | } |
6d45eaef | 89 | |
46d29e70 | 90 | //_____________________________________________________________________________ |
43bfc8af | 91 | AliTRDtrack::AliTRDtrack(AliTRDcluster *c, Int_t index |
fe33d239 | 92 | , const Double_t p[5], const Double_t cov[15] |
93 | , Double_t x, Double_t alpha) | |
94 | :AliKalmanTrack() | |
95 | ,fSeedLab(-1) | |
96 | ,fdEdx(0) | |
97 | ,fDE(0) | |
79143cbf | 98 | ,fClusterOwner(kFALSE) |
fe33d239 | 99 | ,fStopped(kFALSE) |
100 | ,fLhElectron(0) | |
101 | ,fNWrong(0) | |
102 | ,fNRotate(0) | |
103 | ,fNCross(0) | |
104 | ,fNExpected(0) | |
105 | ,fNLast(0) | |
106 | ,fNExpectedLast(0) | |
107 | ,fNdedx(0) | |
108 | ,fChi2Last(1e10) | |
43bfc8af | 109 | ,fBackupTrack(0x0) |
15ed8ba1 | 110 | { |
fe33d239 | 111 | // |
112 | // The main AliTRDtrack constructor. | |
113 | // | |
114 | ||
af26ce80 | 115 | Double_t cnv = 1.0 / (GetBz() * kB2C); |
fe33d239 | 116 | |
117 | Double_t pp[5] = { p[0] | |
118 | , p[1] | |
119 | , x*p[4] - p[2] | |
120 | , p[3] | |
121 | , p[4]*cnv }; | |
6c94f330 | 122 | |
af26ce80 | 123 | Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[ 5]; |
124 | Double_t c32 = x*cov[13] - cov[ 8]; | |
125 | Double_t c20 = x*cov[10] - cov[ 3]; | |
126 | Double_t c21 = x*cov[11] - cov[ 4]; | |
127 | Double_t c42 = x*cov[14] - cov[12]; | |
6c94f330 | 128 | |
af26ce80 | 129 | Double_t cc[15] = { cov[ 0] |
130 | , cov[ 1], cov[ 2] | |
fe33d239 | 131 | , c20, c21, c22 |
af26ce80 | 132 | , cov[ 6], cov[ 7], c32, cov[ 9] |
fe33d239 | 133 | , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv }; |
6c94f330 | 134 | |
135 | Set(x,alpha,pp,cc); | |
5443e65e | 136 | SetNumberOfClusters(1); |
af26ce80 | 137 | fIndex[0] = index; |
138 | fClusters[0] = c; | |
5443e65e | 139 | |
fe33d239 | 140 | for (Int_t i = 0; i < kNplane; i++) { |
141 | for (Int_t j = 0; j < kNslice; j++) { | |
142 | fdEdxPlane[i][j] = 0.0; | |
6d45eaef | 143 | } |
144 | fTimBinPlane[i] = -1; | |
af26ce80 | 145 | fMom[i] = -1.; |
146 | fSnp[i] = 0.; | |
147 | fTgl[i] = 0.; | |
eab5961e | 148 | } |
a2b90f83 | 149 | |
5443e65e | 150 | Double_t q = TMath::Abs(c->GetQ()); |
fe33d239 | 151 | Double_t s = GetSnp(); |
152 | Double_t t = GetTgl(); | |
153 | if (s*s < 1) { | |
154 | q *= TMath::Sqrt((1-s*s)/(1+t*t)); | |
155 | } | |
53c17fbf | 156 | |
af26ce80 | 157 | fdQdl[0] = q; |
fe33d239 | 158 | for (UInt_t i = 1; i < kMAXCLUSTERSPERTRACK; i++) { |
159 | fdQdl[i] = 0; | |
160 | fIndex[i] = 0; | |
161 | fIndexBackup[i] = 0; | |
af26ce80 | 162 | fClusters[i] = 0x0; |
15ed8ba1 | 163 | } |
fe33d239 | 164 | |
165 | for (Int_t i = 0; i < 3;i++) { | |
166 | fBudget[i] = 0; | |
167 | } | |
168 | ||
46d29e70 | 169 | } |
170 | ||
171 | //_____________________________________________________________________________ | |
43bfc8af | 172 | AliTRDtrack::AliTRDtrack(const AliTRDtrack &t/*, const Bool_t owner*/) |
fe33d239 | 173 | :AliKalmanTrack(t) |
174 | ,fSeedLab(t.GetSeedLabel()) | |
175 | ,fdEdx(t.fdEdx) | |
176 | ,fDE(t.fDE) | |
79143cbf | 177 | ,fClusterOwner(kTRUE) |
fe33d239 | 178 | ,fStopped(t.fStopped) |
179 | ,fLhElectron(0) | |
180 | ,fNWrong(t.fNWrong) | |
181 | ,fNRotate(t.fNRotate) | |
182 | ,fNCross(t.fNCross) | |
183 | ,fNExpected(t.fNExpected) | |
184 | ,fNLast(t.fNLast) | |
185 | ,fNExpectedLast(t.fNExpectedLast) | |
186 | ,fNdedx(t.fNdedx) | |
187 | ,fChi2Last(t.fChi2Last) | |
43bfc8af | 188 | ,fBackupTrack(0x0) |
53c17fbf | 189 | { |
46d29e70 | 190 | // |
191 | // Copy constructor. | |
192 | // | |
fe33d239 | 193 | |
194 | for (Int_t i = 0; i < kNplane; i++) { | |
195 | for (Int_t j = 0; j < kNslice; j++) { | |
6d45eaef | 196 | fdEdxPlane[i][j] = t.fdEdxPlane[i][j]; |
197 | } | |
198 | fTimBinPlane[i] = t.fTimBinPlane[i]; | |
199 | fTracklets[i] = t.fTracklets[i]; | |
af26ce80 | 200 | fMom[i] = t.fMom[i]; |
201 | fSnp[i] = t.fSnp[i]; | |
202 | fTgl[i] = t.fTgl[i]; | |
eab5961e | 203 | } |
46d29e70 | 204 | |
fe33d239 | 205 | Int_t n = t.GetNumberOfClusters(); |
6c94f330 | 206 | SetNumberOfClusters(n); |
fe33d239 | 207 | |
208 | for (Int_t i = 0; i < n; i++) { | |
209 | fIndex[i] = t.fIndex[i]; | |
210 | fIndexBackup[i] = t.fIndex[i]; | |
211 | fdQdl[i] = t.fdQdl[i]; | |
af26ce80 | 212 | if (fClusterOwner && t.fClusters[i]) { |
213 | fClusters[i] = new AliTRDcluster(*(t.fClusters[i])); | |
214 | } | |
215 | else { | |
216 | fClusters[i] = t.fClusters[i]; | |
217 | } | |
218 | } | |
fe33d239 | 219 | |
220 | for (UInt_t i = n; i < kMAXCLUSTERSPERTRACK; i++) { | |
221 | fdQdl[i] = 0; | |
222 | fIndex[i] = 0; | |
223 | fIndexBackup[i] = 0; | |
af26ce80 | 224 | fClusters[i] = 0x0; |
03b0452e | 225 | } |
15ed8ba1 | 226 | |
fe33d239 | 227 | for (Int_t i = 0; i < 3;i++) { |
228 | fBudget[i] = t.fBudget[i]; | |
15ed8ba1 | 229 | } |
fe33d239 | 230 | |
5443e65e | 231 | } |
232 | ||
233 | //_____________________________________________________________________________ | |
fe33d239 | 234 | AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t /*alpha*/) |
235 | :AliKalmanTrack(t) | |
236 | ,fSeedLab(-1) | |
237 | ,fdEdx(t.GetPIDsignal()) | |
238 | ,fDE(0) | |
79143cbf | 239 | ,fClusterOwner(kFALSE) |
fe33d239 | 240 | ,fStopped(kFALSE) |
241 | ,fLhElectron(0.0) | |
242 | ,fNWrong(0) | |
243 | ,fNRotate(0) | |
244 | ,fNCross(0) | |
245 | ,fNExpected(0) | |
246 | ,fNLast(0) | |
247 | ,fNExpectedLast(0) | |
248 | ,fNdedx(0) | |
249 | ,fChi2Last(0.0) | |
250 | ,fBackupTrack(0x0) | |
53c17fbf | 251 | { |
5443e65e | 252 | // |
fe33d239 | 253 | // Constructor from AliTPCtrack or AliITStrack |
5443e65e | 254 | // |
255 | ||
6c94f330 | 256 | SetLabel(t.GetLabel()); |
fe33d239 | 257 | SetChi2(0.0); |
6c94f330 | 258 | SetMass(t.GetMass()); |
5443e65e | 259 | SetNumberOfClusters(0); |
260 | ||
fe33d239 | 261 | for (Int_t i = 0; i < kNplane; i++) { |
262 | for (Int_t j = 0; j < kNslice; j++) { | |
6d45eaef | 263 | fdEdxPlane[i][j] = 0.0; |
264 | } | |
eab5961e | 265 | fTimBinPlane[i] = -1; |
af26ce80 | 266 | fMom[i] = -1.; |
267 | fSnp[i] = 0.; | |
268 | fTgl[i] = 0.; | |
eab5961e | 269 | } |
5443e65e | 270 | |
fe33d239 | 271 | for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) { |
272 | fdQdl[i] = 0; | |
273 | fIndex[i] = 0; | |
274 | fIndexBackup[i] = 0; | |
af26ce80 | 275 | fClusters[i] = 0x0; |
79e94bf8 | 276 | } |
4f1c04d3 | 277 | |
fe33d239 | 278 | for (Int_t i = 0; i < 3; i++) { |
279 | fBudget[i] = 0; | |
280 | } | |
281 | ||
79e94bf8 | 282 | } |
53c17fbf | 283 | |
79e94bf8 | 284 | //_____________________________________________________________________________ |
fe33d239 | 285 | AliTRDtrack::AliTRDtrack(const AliESDtrack &t) |
286 | :AliKalmanTrack() | |
287 | ,fSeedLab(-1) | |
288 | ,fdEdx(t.GetTRDsignal()) | |
289 | ,fDE(0) | |
79143cbf | 290 | ,fClusterOwner(kFALSE) |
fe33d239 | 291 | ,fStopped(kFALSE) |
292 | ,fLhElectron(0) | |
293 | ,fNWrong(0) | |
294 | ,fNRotate(0) | |
295 | ,fNCross(0) | |
296 | ,fNExpected(0) | |
297 | ,fNLast(0) | |
298 | ,fNExpectedLast(0) | |
299 | ,fNdedx(0) | |
300 | ,fChi2Last(1e10) | |
301 | ,fBackupTrack(0x0) | |
53c17fbf | 302 | { |
79e94bf8 | 303 | // |
304 | // Constructor from AliESDtrack | |
305 | // | |
fe33d239 | 306 | |
79e94bf8 | 307 | SetLabel(t.GetLabel()); |
fe33d239 | 308 | SetChi2(0.0); |
79e94bf8 | 309 | SetMass(t.GetMass()); |
1e9bb598 | 310 | SetNumberOfClusters(t.GetTRDclusters(fIndex)); |
15ed8ba1 | 311 | |
46e2d86c | 312 | Int_t ncl = t.GetTRDclusters(fIndexBackup); |
fe33d239 | 313 | for (UInt_t i = ncl; i < kMAXCLUSTERSPERTRACK; i++) { |
314 | fIndexBackup[i] = 0; | |
315 | fIndex[i] = 0; | |
15ed8ba1 | 316 | } |
fe33d239 | 317 | |
318 | for (Int_t i = 0; i < kNplane; i++) { | |
319 | for (Int_t j = 0; j < kNslice; j++) { | |
6d45eaef | 320 | fdEdxPlane[i][j] = t.GetTRDsignals(i,j); |
321 | } | |
eab5961e | 322 | fTimBinPlane[i] = t.GetTRDTimBin(i); |
af26ce80 | 323 | fMom[i] = -1.; |
324 | fSnp[i] = 0.; | |
325 | fTgl[i] = 0.; | |
eab5961e | 326 | } |
79e94bf8 | 327 | |
fe33d239 | 328 | const AliExternalTrackParam *par = &t; |
329 | if (t.GetStatus() & AliESDtrack::kTRDbackup) { | |
330 | par = t.GetOuterParam(); | |
331 | if (!par) { | |
332 | AliError("No backup info!"); | |
333 | par = &t; | |
334 | } | |
c5a8e3df | 335 | } |
6c94f330 | 336 | Set(par->GetX(),par->GetAlpha(),par->GetParameter(),par->GetCovariance()); |
79e94bf8 | 337 | |
fe33d239 | 338 | for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) { |
43bfc8af | 339 | fdQdl[i] = 0; |
af26ce80 | 340 | fClusters[i] = 0x0; |
341 | } | |
fe33d239 | 342 | |
343 | for (Int_t i = 0; i < 3; i++) { | |
344 | fBudget[i] = 0; | |
345 | } | |
5443e65e | 346 | |
fe33d239 | 347 | if ((t.GetStatus() & AliESDtrack::kTIME) == 0) { |
348 | return; | |
349 | } | |
c630aafd | 350 | |
c630aafd | 351 | StartTimeIntegral(); |
fe33d239 | 352 | Double_t times[10]; |
353 | t.GetIntegratedTimes(times); | |
354 | SetIntegratedTimes(times); | |
c630aafd | 355 | SetIntegratedLength(t.GetIntegratedLength()); |
356 | ||
16d9fbba | 357 | } |
358 | ||
53c17fbf | 359 | //____________________________________________________________________________ |
360 | AliTRDtrack::~AliTRDtrack() | |
3fad3d32 | 361 | { |
362 | // | |
53c17fbf | 363 | // Destructor |
364 | // | |
3fad3d32 | 365 | |
fe33d239 | 366 | if (fBackupTrack) { |
367 | delete fBackupTrack; | |
368 | } | |
43bfc8af | 369 | fBackupTrack = 0x0; |
af26ce80 | 370 | if (fClusterOwner) { |
371 | UInt_t icluster = 0; | |
372 | while ((icluster < kMAXCLUSTERSPERTRACK) && fClusters[icluster]) { | |
373 | delete fClusters[icluster]; | |
374 | fClusters[icluster] = 0x0; | |
375 | icluster++; | |
376 | } | |
43bfc8af | 377 | } |
3fad3d32 | 378 | |
53c17fbf | 379 | } |
380 | ||
53c17fbf | 381 | //____________________________________________________________________________ |
382 | Float_t AliTRDtrack::StatusForTOF() | |
7ad19338 | 383 | { |
53c17fbf | 384 | // |
385 | // Defines the status of the TOF extrapolation | |
386 | // | |
03b0452e | 387 | |
fe33d239 | 388 | // Definition of res ???? |
389 | Float_t res = (0.2 + 0.8 * (fN / (fNExpected + 5.0))) | |
390 | * (0.4 + 0.6 * fTracklets[5].GetN() / 20.0); | |
391 | res *= (0.25 + 0.8 * 40.0 / (40.0 + fBudget[2])); | |
03b0452e | 392 | return res; |
393 | ||
fe33d239 | 394 | // This part of the function is never reached ???? |
395 | // What defines these parameters ???? | |
af26ce80 | 396 | //Int_t status = 0; |
397 | //if (GetNumberOfClusters() < 20) return 0; | |
398 | //if ((fN > 110) && | |
399 | // (fChi2/(Float_t(fN)) < 3)) return 3; // Gold | |
400 | //if ((fNLast > 30) && | |
401 | // (fChi2Last/(Float_t(fNLast)) < 3)) return 3; // Gold | |
402 | //if ((fNLast > 20) && | |
403 | // (fChi2Last/(Float_t(fNLast)) < 2)) return 3; // Gold | |
404 | //if ((fNLast/(fNExpectedLast+3.0) > 0.8) && | |
405 | // (fChi2Last/Float_t(fNLast) < 5) && | |
406 | // (fNLast > 20)) return 2; // Silver | |
407 | //if ((fNLast > 5) && | |
408 | // (((fNLast+1.0)/(fNExpectedLast+1.0)) > 0.8) && | |
409 | // (fChi2Last/(fNLast-5.0) < 6)) return 1; | |
410 | // | |
411 | //return status; | |
53c17fbf | 412 | |
4f1c04d3 | 413 | } |
46d29e70 | 414 | |
415 | //_____________________________________________________________________________ | |
53c17fbf | 416 | Int_t AliTRDtrack::Compare(const TObject *o) const |
417 | { | |
418 | // | |
419 | // Compares tracks according to their Y2 or curvature | |
420 | // | |
46d29e70 | 421 | |
d49e463b | 422 | AliTRDtrack *t = (AliTRDtrack *) o; |
46d29e70 | 423 | |
d49e463b | 424 | Double_t co = TMath::Abs(t->GetC()); |
425 | Double_t c = TMath::Abs(GetC()); | |
46d29e70 | 426 | |
d49e463b | 427 | if (c > co) { |
428 | return 1; | |
429 | } | |
430 | else if (c < co) { | |
431 | return -1; | |
432 | } | |
46d29e70 | 433 | return 0; |
53c17fbf | 434 | |
46d29e70 | 435 | } |
436 | ||
a819a5f7 | 437 | //_____________________________________________________________________________ |
43bfc8af | 438 | void AliTRDtrack::CookdEdx(Double_t low, Double_t up) |
15ed8ba1 | 439 | { |
440 | // | |
d49e463b | 441 | // Calculates the truncated dE/dx within the "low" and "up" cuts. |
15ed8ba1 | 442 | // |
443 | ||
d49e463b | 444 | // Array to sort the dEdx values according to amplitude |
af26ce80 | 445 | Float_t sorted[kMAXCLUSTERSPERTRACK]; |
446 | fdEdx = 0.; | |
43bfc8af | 447 | |
15ed8ba1 | 448 | // Require at least 10 clusters for a dedx measurement |
c6011b06 | 449 | if (fNdedx < 10) return; |
15ed8ba1 | 450 | |
d49e463b | 451 | // Can fdQdl be negative ???? |
af26ce80 | 452 | for (Int_t i = 0; i < fNdedx; i++) { |
453 | sorted[i] = TMath::Abs(fdQdl[i]); | |
454 | } | |
15ed8ba1 | 455 | // Sort the dedx values by amplitude |
c6011b06 | 456 | Int_t *index = new Int_t[fNdedx]; |
457 | TMath::Sort(fNdedx, sorted, index, kFALSE); | |
458 | ||
43bfc8af | 459 | // Sum up the truncated charge between lower and upper bounds |
c6011b06 | 460 | Int_t nl = Int_t(low * fNdedx); |
461 | Int_t nu = Int_t( up * fNdedx); | |
af26ce80 | 462 | for (Int_t i = nl; i <= nu; i++) { |
463 | fdEdx += sorted[index[i]]; | |
464 | } | |
465 | fdEdx /= (nu - nl + 1.0); | |
466 | ||
467 | delete[] index; | |
43bfc8af | 468 | |
43bfc8af | 469 | } |
470 | ||
471 | //_____________________________________________________________________________ | |
44dbae42 | 472 | void AliTRDtrack::CookdEdxTimBin(const Int_t/* tid*/) |
43bfc8af | 473 | { |
474 | // | |
475 | // Set fdEdxPlane and fTimBinPlane and also get the | |
476 | // Time bin for Max. Cluster | |
af26ce80 | 477 | // |
478 | // Authors: | |
43bfc8af | 479 | // Prashant Shukla (shukla@physi.uni-heidelberg.de) |
480 | // Alexandru Bercuci (A.Bercuci@gsi.de) | |
af26ce80 | 481 | // |
43bfc8af | 482 | |
af26ce80 | 483 | Double_t maxcharge[AliESDtrack::kNPlane]; // max charge in chamber |
484 | // Number of clusters attached to track per chamber and slice | |
485 | Int_t nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice]; | |
486 | // Number of time bins in chamber | |
487 | Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBins(); | |
488 | Int_t plane; // Plane of current cluster | |
489 | Int_t tb; // Time bin of current cluster | |
490 | Int_t slice; // Current slice | |
491 | AliTRDcluster *cluster = 0x0; // Pointer to current cluster | |
492 | ||
493 | // Reset class and local counters/variables | |
43bfc8af | 494 | for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) { |
af26ce80 | 495 | fTimBinPlane[iPlane] = -1; |
496 | maxcharge[iPlane] = 0.; | |
497 | for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) { | |
43bfc8af | 498 | fdEdxPlane[iPlane][iSlice] = 0.; |
499 | nCluster[iPlane][iSlice] = 0; | |
500 | } | |
15ed8ba1 | 501 | } |
c6011b06 | 502 | |
af26ce80 | 503 | // Start looping over clusters attached to this track |
504 | for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) { | |
505 | ||
43bfc8af | 506 | cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]); |
507 | if(!cluster) continue; | |
508 | ||
af26ce80 | 509 | // Read info from current cluster |
510 | plane = AliTRDgeometry::GetPlane(cluster->GetDetector()); | |
43bfc8af | 511 | if (plane < 0 || plane >= AliESDtrack::kNPlane) { |
512 | AliError(Form("Wrong plane %d", plane)); | |
513 | continue; | |
514 | } | |
515 | ||
af26ce80 | 516 | tb = cluster->GetLocalTimeBin(); |
517 | if ((tb == 0) || (tb >= ntb)) { | |
518 | AliWarning(Form("time bin 0 or > %d in cluster %d", ntb, iClus)); | |
519 | AliInfo(Form("dQ/dl %f", fdQdl[iClus])); | |
520 | continue; | |
521 | } | |
43bfc8af | 522 | |
af26ce80 | 523 | slice = tb * AliESDtrack::kNSlice / ntb; |
524 | ||
525 | fdEdxPlane[plane][slice] += fdQdl[iClus]; | |
526 | if (fdQdl[iClus] > maxcharge[plane]) { | |
527 | maxcharge[plane] = fdQdl[iClus]; | |
43bfc8af | 528 | fTimBinPlane[plane] = tb; |
529 | } | |
43bfc8af | 530 | |
af26ce80 | 531 | nCluster[plane][slice]++; |
43bfc8af | 532 | |
af26ce80 | 533 | } // End of loop over cluster |
43bfc8af | 534 | |
535 | // Normalize fdEdxPlane to number of clusters and set track segments | |
536 | for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) { | |
537 | for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) { | |
af26ce80 | 538 | if (nCluster[iPlane][iSlice]) { |
539 | fdEdxPlane[iPlane][iSlice] /= nCluster[iPlane][iSlice]; | |
540 | } | |
43bfc8af | 541 | } |
af26ce80 | 542 | } |
43bfc8af | 543 | |
af26ce80 | 544 | } |
43bfc8af | 545 | |
44dbae42 | 546 | //_____________________________________________________________________________ |
547 | void AliTRDtrack::CookdEdxNN(Float_t *dedx){ | |
548 | ||
549 | // This function calcuates the input for the neural networks | |
550 | // which are used for the PID. | |
551 | ||
552 | Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();//number of time bins in chamber | |
553 | Int_t plane; // plane of current cluster | |
554 | Int_t tb; // time bin of current cluster | |
555 | Int_t slice; // curent slice | |
556 | AliTRDcluster *cluster = 0x0; // pointer to current cluster | |
557 | const Int_t lMLPscale = 16000; // scaling of the MLP input to be smaller than 1 | |
558 | ||
559 | // Reset class and local contors/variables | |
560 | for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++){ | |
561 | for (Int_t iSlice = 0; iSlice < kNMLPslice; iSlice++) { | |
562 | *(dedx + (kNMLPslice * iPlane) + iSlice) = 0.; | |
563 | } | |
564 | } | |
565 | ||
566 | // start looping over clusters attached to this track | |
567 | for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) { | |
568 | cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]); | |
569 | if(!cluster) continue; | |
570 | ||
571 | // Read info from current cluster | |
572 | plane = AliTRDgeometry::GetPlane(cluster->GetDetector()); | |
573 | if (plane < 0 || plane >= AliESDtrack::kNPlane) { | |
574 | AliError(Form("Wrong plane %d", plane)); | |
575 | continue; | |
576 | } | |
577 | ||
578 | tb = cluster->GetLocalTimeBin(); | |
579 | if(tb == 0 || tb >= ntb){ | |
580 | AliWarning(Form("time bin 0 or > %d in cluster %d", ntb, iClus)); | |
581 | AliInfo(Form("dQ/dl %f", fdQdl[iClus])); | |
582 | continue; | |
583 | } | |
584 | ||
585 | slice = tb * kNMLPslice / ntb; | |
586 | ||
587 | *(dedx+(kNMLPslice * plane) + slice) += fdQdl[iClus]/lMLPscale; | |
588 | } // End of loop over cluster | |
589 | ||
590 | } | |
591 | ||
592 | ||
43bfc8af | 593 | //_____________________________________________________________________________ |
594 | void AliTRDtrack::SetTrackSegmentDirMom(const Int_t plane) | |
595 | { | |
af26ce80 | 596 | // |
597 | // Set the momenta for a track segment in a given plane | |
598 | // | |
599 | ||
600 | if ((plane < 0) || | |
601 | (plane>= kNplane)) { | |
602 | AliError(Form("Trying to access out of range plane (%d)", plane)); | |
603 | return; | |
604 | } | |
43bfc8af | 605 | |
af26ce80 | 606 | fSnp[plane] = GetSnp(); |
607 | fTgl[plane] = GetTgl(); | |
608 | Double_t p[3]; | |
609 | GetPxPyPz(p); | |
610 | fMom[plane] = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); | |
611 | ||
43bfc8af | 612 | } |
613 | ||
614 | //_____________________________________________________________________________ | |
615 | Float_t AliTRDtrack::GetTrackLengthPlane(Int_t plane) const | |
616 | { | |
af26ce80 | 617 | // |
618 | // Calculate the track length of a track segment in a given plane | |
619 | // | |
43bfc8af | 620 | |
af26ce80 | 621 | if ((plane < 0) || (plane >= kNplane)) return 0.; |
43bfc8af | 622 | |
af26ce80 | 623 | return (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick()) |
624 | / TMath::Sqrt((1.0 - fSnp[plane]*fSnp[plane]) | |
625 | / (1.0 + fTgl[plane]*fTgl[plane])); | |
626 | ||
627 | } | |
43bfc8af | 628 | |
629 | //_____________________________________________________________________________ | |
44dbae42 | 630 | Bool_t AliTRDtrack::CookPID(Int_t &pidQuality) |
43bfc8af | 631 | { |
632 | // | |
633 | // This function calculates the PID probabilities based on TRD signals | |
634 | // | |
635 | // The method produces probabilities based on the charge | |
636 | // and the position of the maximum time bin in each layer. | |
637 | // The dE/dx information can be used as global charge or 2 to 3 | |
720a0a16 | 638 | // slices. Check AliTRDCalPID and AliTRDCalPIDRefMaker for the actual |
43bfc8af | 639 | // implementation. |
640 | // | |
641 | // Author | |
642 | // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007 | |
af26ce80 | 643 | // |
43bfc8af | 644 | |
af26ce80 | 645 | AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); |
646 | if (!calibration) { | |
647 | AliError("No access to calibration data"); | |
44dbae42 | 648 | return kFALSE; |
af26ce80 | 649 | } |
43bfc8af | 650 | |
44dbae42 | 651 | // Retrieve the CDB container class with the probability distributions |
652 | const AliTRDCalPID *pd = calibration->GetPIDObject(fPIDmethod == kNN ? 0 : 1); | |
653 | if (!pd) { | |
654 | AliError("No access to AliTRDCalPID"); | |
655 | return kFALSE; | |
656 | } | |
43bfc8af | 657 | |
44dbae42 | 658 | // Calculate the input for the NN if fPIDmethod is kNN |
659 | Float_t ldEdxNN[AliTRDCalPID::kNPlane * kNMLPslice], *dedx = 0x0; | |
660 | if(fPIDmethod == kNN) { | |
661 | CookdEdxNN(&ldEdxNN[0]); | |
662 | } | |
43bfc8af | 663 | |
44dbae42 | 664 | // Sets the a priori probabilities |
665 | for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) { | |
666 | fPID[ispec] = 1./AliPID::kSPECIES; | |
667 | } | |
43bfc8af | 668 | |
43bfc8af | 669 | |
44dbae42 | 670 | if(AliPID::kSPECIES != 5){ |
671 | AliError("Probabilities array defined only for 5 values. Please modify !!"); | |
672 | return kFALSE; | |
673 | } | |
43bfc8af | 674 | |
43bfc8af | 675 | |
44dbae42 | 676 | pidQuality = 0; |
677 | Float_t length; // track segment length in chamber | |
43bfc8af | 678 | |
44dbae42 | 679 | // Skip tracks which have no TRD signal at all |
680 | if (fdEdx == 0.) return kFALSE; | |
681 | ||
682 | for (Int_t iPlane = 0; iPlane < AliTRDgeometry::kNplan; iPlane++) { | |
683 | ||
684 | length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/TMath::Sqrt((1. - fSnp[iPlane]*fSnp[iPlane]) / (1. + fTgl[iPlane]*fTgl[iPlane])); | |
685 | ||
686 | // check data | |
687 | if((fdEdxPlane[iPlane][0] + fdEdxPlane[iPlane][1] + fdEdxPlane[iPlane][2]) <= 0. | |
688 | || fTimBinPlane[iPlane] == -1.) continue; | |
689 | ||
690 | // this track segment has fulfilled all requierments | |
691 | pidQuality++; | |
692 | ||
693 | // Get the probabilities for the different particle species | |
694 | for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) { | |
695 | switch(fPIDmethod){ | |
696 | case kLQ: | |
697 | dedx = fdEdxPlane[iPlane]; | |
698 | break; | |
699 | case kNN: | |
700 | dedx = &ldEdxNN[iPlane*kNMLPslice]; | |
701 | break; | |
702 | } | |
703 | fPID[iSpecies] *= pd->GetProbability(iSpecies, fMom[iPlane], dedx, length, iPlane); | |
704 | } | |
705 | } | |
706 | if(pidQuality == 0) return kTRUE; | |
43bfc8af | 707 | |
44dbae42 | 708 | // normalize probabilities |
709 | Double_t probTotal = 0.; | |
710 | for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += fPID[iSpecies]; | |
43bfc8af | 711 | |
44dbae42 | 712 | if(probTotal <= 0.) { |
713 | AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference histograms."); | |
714 | return kFALSE; | |
715 | } | |
716 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) fPID[iSpecies] /= probTotal; | |
43bfc8af | 717 | |
44dbae42 | 718 | return kTRUE; |
43bfc8af | 719 | } |
720 | ||
44dbae42 | 721 | |
46d29e70 | 722 | //_____________________________________________________________________________ |
826fe01c | 723 | Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho) |
46d29e70 | 724 | { |
fe33d239 | 725 | // |
826fe01c | 726 | // Propagates this track to a reference plane defined by "xk" [cm] |
727 | // correcting for the mean crossed material. | |
fe33d239 | 728 | // |
826fe01c | 729 | // "xx0" - thickness/rad.length [units of the radiation length] |
730 | // "xrho" - thickness*density [g/cm^2] | |
731 | // | |
fe33d239 | 732 | |
733 | if (xk == GetX()) { | |
734 | return kTRUE; | |
735 | } | |
15ed8ba1 | 736 | |
fe33d239 | 737 | Double_t oldX = GetX(); |
738 | Double_t oldY = GetY(); | |
739 | Double_t oldZ = GetZ(); | |
15ed8ba1 | 740 | |
fe33d239 | 741 | Double_t bz = GetBz(); |
46d29e70 | 742 | |
fe33d239 | 743 | if (!AliExternalTrackParam::PropagateTo(xk,bz)) { |
744 | return kFALSE; | |
745 | } | |
9c9d2487 | 746 | |
fe33d239 | 747 | Double_t x = GetX(); |
748 | Double_t y = GetY(); | |
749 | Double_t z = GetZ(); | |
826fe01c | 750 | |
fe33d239 | 751 | if (oldX < xk) { |
826fe01c | 752 | xrho = -xrho; |
fe33d239 | 753 | if (IsStartedTimeIntegral()) { |
8234f064 | 754 | Double_t l2 = TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY) + (z-oldZ)*(z-oldZ)); |
fe33d239 | 755 | Double_t crv = GetC(); |
756 | if (TMath::Abs(l2*crv) > 0.0001) { | |
757 | // Make correction for curvature if neccesary | |
758 | l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY)); | |
759 | l2 = 2.0 * TMath::ASin(l2 * crv) / crv; | |
760 | l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ)); | |
761 | } | |
762 | AddTimeStep(l2); | |
6c94f330 | 763 | } |
46d29e70 | 764 | } |
15ed8ba1 | 765 | |
826fe01c | 766 | if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0,xrho,GetMass())) { |
fe33d239 | 767 | return kFALSE; |
768 | } | |
15ed8ba1 | 769 | |
fe33d239 | 770 | { |
771 | ||
772 | // Energy losses************************ | |
6c23ffed | 773 | Double_t p2 = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt()); |
fe33d239 | 774 | Double_t beta2 = p2 / (p2 + GetMass()*GetMass()); |
d248a39d | 775 | if (beta2<1.0e-10 || (5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0) { |
fe33d239 | 776 | return kFALSE; |
777 | } | |
778 | ||
779 | Double_t dE = 0.153e-3 / beta2 | |
780 | * (log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2) | |
826fe01c | 781 | * xrho; |
782 | fBudget[0] += xrho; | |
fe33d239 | 783 | |
784 | /* | |
785 | // Suspicious part - think about it ? | |
786 | Double_t kinE = TMath::Sqrt(p2); | |
787 | if (dE > 0.8*kinE) dE = 0.8 * kinE; // | |
788 | if (dE < 0) dE = 0.0; // Not valid region for Bethe bloch | |
789 | */ | |
790 | ||
791 | fDE += dE; | |
792 | ||
793 | /* | |
794 | // Suspicious ! I.B. | |
795 | Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE)); // Energy loss fluctuation | |
796 | Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2); | |
797 | fCcc += sigmac2; | |
798 | fCee += fX*fX * sigmac2; | |
799 | */ | |
6c94f330 | 800 | |
0d5b5c27 | 801 | } |
802 | ||
6c94f330 | 803 | return kTRUE; |
6d45eaef | 804 | |
46d29e70 | 805 | } |
806 | ||
46d29e70 | 807 | //_____________________________________________________________________________ |
fe33d239 | 808 | Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index |
809 | , Double_t h01) | |
46d29e70 | 810 | { |
fe33d239 | 811 | // |
6c94f330 | 812 | // Assignes found cluster to the track and updates track information |
fe33d239 | 813 | // |
46d29e70 | 814 | |
b8dc2353 | 815 | Bool_t fNoTilt = kTRUE; |
fe33d239 | 816 | if (TMath::Abs(h01) > 0.003) { |
817 | fNoTilt = kFALSE; | |
818 | } | |
819 | ||
820 | // Add angular effect to the error contribution - MI | |
6c94f330 | 821 | Float_t tangent2 = GetSnp()*GetSnp(); |
fe33d239 | 822 | if (tangent2 < 0.90000) { |
823 | tangent2 = tangent2 / (1.0 - tangent2); | |
b8dc2353 | 824 | } |
fe33d239 | 825 | //Float_t errang = tangent2 * 0.04; |
b8dc2353 | 826 | |
fe33d239 | 827 | Double_t p[2] = {c->GetY(), c->GetZ() }; |
828 | //Double_t cov[3] = {c->GetSigmaY2()+errang, 0.0, c->GetSigmaZ2()*100.0 }; | |
829 | Double_t sy2 = c->GetSigmaY2() * 4.0; | |
830 | Double_t sz2 = c->GetSigmaZ2() * 4.0; | |
af26ce80 | 831 | Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 }; |
46e2d86c | 832 | |
fe33d239 | 833 | if (!AliExternalTrackParam::Update(p,cov)) { |
834 | return kFALSE; | |
835 | } | |
46d29e70 | 836 | |
af26ce80 | 837 | Int_t n = GetNumberOfClusters(); |
fe33d239 | 838 | fIndex[n] = index; |
b8dc2353 | 839 | SetNumberOfClusters(n+1); |
af26ce80 | 840 | |
b8dc2353 | 841 | SetChi2(GetChi2()+chisq); |
5443e65e | 842 | |
af26ce80 | 843 | return kTRUE; |
fe33d239 | 844 | |
845 | } | |
53c17fbf | 846 | |
46e2d86c | 847 | //_____________________________________________________________________________ |
af26ce80 | 848 | Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index |
44dbae42 | 849 | , Double_t h01, Int_t /*plane*/, Int_t tid) |
fe33d239 | 850 | { |
851 | // | |
46e2d86c | 852 | // Assignes found cluster to the track and updates track information |
fe33d239 | 853 | // |
854 | ||
46e2d86c | 855 | Bool_t fNoTilt = kTRUE; |
fe33d239 | 856 | if (TMath::Abs(h01) > 0.003) { |
857 | fNoTilt = kFALSE; | |
858 | } | |
859 | ||
860 | // Add angular effect to the error contribution and make correction - MI | |
6c94f330 | 861 | Double_t tangent2 = GetSnp()*GetSnp(); |
862 | if (tangent2 < 0.90000){ | |
fe33d239 | 863 | tangent2 = tangent2 / (1.0-tangent2); |
15ed8ba1 | 864 | } |
6c94f330 | 865 | Double_t tangent = TMath::Sqrt(tangent2); |
fe33d239 | 866 | if (GetSnp() < 0) { |
867 | tangent *= -1; | |
868 | } | |
869 | ||
870 | // | |
871 | // Is the following still needed ???? | |
872 | // | |
873 | ||
6c94f330 | 874 | // Double_t correction = 0*plane; |
875 | /* | |
876 | Double_t errang = tangent2*0.04; // | |
877 | Double_t errsys =0.025*0.025*20; //systematic error part | |
15ed8ba1 | 878 | |
6c94f330 | 879 | Float_t extend =1; |
880 | if (c->GetNPads()==4) extend=2; | |
881 | */ | |
882 | //if (c->GetNPads()==5) extend=3; | |
883 | //if (c->GetNPads()==6) extend=3; | |
884 | //if (c->GetQ()<15) return 1; | |
15ed8ba1 | 885 | |
c5a8e3df | 886 | /* |
6c94f330 | 887 | if (corrector!=0){ |
3c625a9b | 888 | //if (0){ |
889 | correction = corrector->GetCorrection(plane,c->GetLocalTimeBin(),tangent); | |
890 | if (TMath::Abs(correction)>0){ | |
891 | //if we have info | |
892 | errang = corrector->GetSigma(plane,c->GetLocalTimeBin(),tangent); | |
893 | errang *= errang; | |
894 | errang += tangent2*0.04; | |
895 | } | |
896 | } | |
c5a8e3df | 897 | */ |
6c94f330 | 898 | // |
899 | //Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.); | |
900 | /* | |
901 | { | |
902 | Double_t dy=c->GetY() - GetY(), dz=c->GetZ() - GetZ(); | |
903 | printf("%e %e %e %e\n",dy,dz,padlength/2,h01); | |
c84a5e9e | 904 | } |
6c94f330 | 905 | */ |
15ed8ba1 | 906 | |
fe33d239 | 907 | Double_t p[2] = { c->GetY(), c->GetZ() }; |
908 | //Double_t cov[3]={ (c->GetSigmaY2()+errang+errsys)*extend, 0.0, c->GetSigmaZ2()*10000.0 }; | |
909 | Double_t sy2 = c->GetSigmaY2() * 4.0; | |
910 | Double_t sz2 = c->GetSigmaZ2() * 4.0; | |
911 | Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 }; | |
46e2d86c | 912 | |
fe33d239 | 913 | if (!AliExternalTrackParam::Update(p,cov)) { |
914 | return kFALSE; | |
915 | } | |
916 | ||
af26ce80 | 917 | // Register cluster to track |
918 | Int_t n = GetNumberOfClusters(); | |
919 | fIndex[n] = index; | |
920 | fClusters[n] = c; | |
46e2d86c | 921 | SetNumberOfClusters(n+1); |
fe33d239 | 922 | SetChi2(GetChi2() + chisq); |
923 | ||
924 | return kTRUE; | |
46e2d86c | 925 | |
53c17fbf | 926 | } |
7ad19338 | 927 | |
fe33d239 | 928 | // //_____________________________________________________________________________ |
929 | // Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet) | |
930 | // { | |
931 | // // | |
932 | // // Assignes found tracklet to the track and updates track information | |
933 | // // | |
934 | // // Can this be removed ???? | |
935 | // // | |
936 | // | |
937 | // Double_t r00=(tracklet.GetTrackletSigma2()), r01=0., r11= 10000.; | |
938 | // r00+=fCyy; r01+=fCzy; r11+=fCzz; | |
939 | // // | |
940 | // Double_t det=r00*r11 - r01*r01; | |
941 | // Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det; | |
942 | // // | |
943 | ||
944 | // Double_t dy=tracklet.GetY() - fY, dz=tracklet.GetZ() - fZ; | |
7ad19338 | 945 | |
15ed8ba1 | 946 | |
fe33d239 | 947 | // Double_t s00 = tracklet.GetTrackletSigma2(); // error pad |
948 | // Double_t s11 = 100000; // error pad-row | |
949 | // Float_t h01 = tracklet.GetTilt(); | |
950 | // // | |
951 | // // r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00; | |
952 | // r00 = fCyy + fCzz*h01*h01+s00; | |
953 | // // r01 = fCzy + fCzz*h01; | |
954 | // r01 = fCzy ; | |
955 | // r11 = fCzz + s11; | |
956 | // det = r00*r11 - r01*r01; | |
957 | // // inverse matrix | |
958 | // tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det; | |
959 | ||
960 | // Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11; | |
961 | // Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11; | |
962 | // Double_t k20=fCey*r00+fCez*r01, k21=fCey*r01+fCez*r11; | |
963 | // Double_t k30=fCty*r00+fCtz*r01, k31=fCty*r01+fCtz*r11; | |
964 | // Double_t k40=fCcy*r00+fCcz*r01, k41=fCcy*r01+fCcz*r11; | |
15ed8ba1 | 965 | |
fe33d239 | 966 | // // K matrix |
967 | // // k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01); | |
968 | // // k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01); | |
969 | // // k20=fCey*r00+fCez*(r01+h01*r00),k21=fCey*r01+fCez*(r11+h01*r01); | |
970 | // // k30=fCty*r00+fCtz*(r01+h01*r00),k31=fCty*r01+fCtz*(r11+h01*r01); | |
971 | // // k40=fCcy*r00+fCcz*(r01+h01*r00),k41=fCcy*r01+fCcz*(r11+h01*r01); | |
972 | // // | |
973 | // //Update measurement | |
974 | // Double_t cur=fC + k40*dy + k41*dz, eta=fE + k20*dy + k21*dz; | |
975 | // // cur=fC + k40*dy + k41*dz; eta=fE + k20*dy + k21*dz; | |
976 | // if (TMath::Abs(cur*fX-eta) >= 0.90000) { | |
977 | // //Int_t n=GetNumberOfClusters(); | |
978 | // // if (n>4) cerr<<n<<" AliTRDtrack warning: Filtering failed !\n"; | |
979 | // return 0; | |
980 | // } | |
981 | // // k01+=h01*k00; | |
982 | // // k11+=h01*k10; | |
983 | // // k21+=h01*k20; | |
984 | // // k31+=h01*k30; | |
985 | // // k41+=h01*k40; | |
986 | ||
987 | ||
988 | // fY += k00*dy + k01*dz; | |
989 | // fZ += k10*dy + k11*dz; | |
990 | // fE = eta; | |
991 | // fT += k30*dy + k31*dz; | |
992 | // fC = cur; | |
7ad19338 | 993 | |
6c94f330 | 994 | |
fe33d239 | 995 | // //Update covariance |
996 | // // | |
997 | // // | |
998 | // Double_t oldyy = fCyy, oldzz = fCzz; //, oldee=fCee, oldcc =fCcc; | |
999 | // Double_t oldzy = fCzy, oldey = fCey, oldty=fCty, oldcy =fCcy; | |
1000 | // Double_t oldez = fCez, oldtz = fCtz, oldcz=fCcz; | |
1001 | // //Double_t oldte = fCte, oldce = fCce; | |
1002 | // //Double_t oldct = fCct; | |
1003 | ||
1004 | // fCyy-=k00*oldyy+k01*oldzy; | |
1005 | // fCzy-=k10*oldyy+k11*oldzy; | |
1006 | // fCey-=k20*oldyy+k21*oldzy; | |
1007 | // fCty-=k30*oldyy+k31*oldzy; | |
1008 | // fCcy-=k40*oldyy+k41*oldzy; | |
1009 | // // | |
1010 | // fCzz-=k10*oldzy+k11*oldzz; | |
1011 | // fCez-=k20*oldzy+k21*oldzz; | |
1012 | // fCtz-=k30*oldzy+k31*oldzz; | |
1013 | // fCcz-=k40*oldzy+k41*oldzz; | |
1014 | // // | |
1015 | // fCee-=k20*oldey+k21*oldez; | |
1016 | // fCte-=k30*oldey+k31*oldez; | |
1017 | // fCce-=k40*oldey+k41*oldez; | |
1018 | // // | |
1019 | // fCtt-=k30*oldty+k31*oldtz; | |
1020 | // fCct-=k40*oldty+k41*oldtz; | |
1021 | // // | |
6c23ffed | 1022 | |
fe33d239 | 1023 | // fCcc-=k40*oldcy+k41*oldcz; |
1024 | // // | |
15ed8ba1 | 1025 | |
fe33d239 | 1026 | // //Int_t n=GetNumberOfClusters(); |
1027 | // //fIndex[n]=index; | |
1028 | // //SetNumberOfClusters(n+1); | |
6c94f330 | 1029 | |
fe33d239 | 1030 | // //SetChi2(GetChi2()+chisq); |
1031 | // // cerr<<"in update: fIndex["<<fN<<"] = "<<index<<endl; | |
7ad19338 | 1032 | |
fe33d239 | 1033 | // return 1; |
7ad19338 | 1034 | |
fe33d239 | 1035 | // } |
c84a5e9e | 1036 | |
46d29e70 | 1037 | //_____________________________________________________________________________ |
6c94f330 | 1038 | Bool_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute) |
46d29e70 | 1039 | { |
fe33d239 | 1040 | // |
6c94f330 | 1041 | // Rotates track parameters in R*phi plane |
1042 | // if absolute rotation alpha is in global system | |
1043 | // otherwise alpha rotation is relative to the current rotation angle | |
fe33d239 | 1044 | // |
1045 | ||
3fad3d32 | 1046 | if (absolute) { |
6c94f330 | 1047 | alpha -= GetAlpha(); |
3fad3d32 | 1048 | } |
1049 | else{ | |
1050 | fNRotate++; | |
1051 | } | |
46d29e70 | 1052 | |
6c94f330 | 1053 | return AliExternalTrackParam::Rotate(GetAlpha()+alpha); |
fe33d239 | 1054 | |
53c17fbf | 1055 | } |
7ad19338 | 1056 | |
46d29e70 | 1057 | //_____________________________________________________________________________ |
fd621f36 | 1058 | Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const |
46d29e70 | 1059 | { |
53c17fbf | 1060 | // |
1061 | // Returns the track chi2 | |
1062 | // | |
1063 | ||
fe33d239 | 1064 | Double_t p[2] = { c->GetY(), c->GetZ() }; |
1065 | Double_t sy2 = c->GetSigmaY2() * 4.0; | |
1066 | Double_t sz2 = c->GetSigmaZ2() * 4.0; | |
1067 | Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 }; | |
6c94f330 | 1068 | |
1069 | return AliExternalTrackParam::GetPredictedChi2(p,cov); | |
15ed8ba1 | 1070 | |
fe33d239 | 1071 | // |
1072 | // Can the following be removed ???? | |
1073 | // | |
6c94f330 | 1074 | /* |
1075 | Bool_t fNoTilt = kTRUE; | |
1076 | if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE; | |
15ed8ba1 | 1077 | |
6c94f330 | 1078 | return (c->GetY() - GetY())*(c->GetY() - GetY())/c->GetSigmaY2(); |
1079 | */ | |
15ed8ba1 | 1080 | |
6c94f330 | 1081 | /* |
1082 | Double_t chi2, dy, r00, r01, r11; | |
b8dc2353 | 1083 | |
6c94f330 | 1084 | if(fNoTilt) { |
1085 | dy=c->GetY() - fY; | |
1086 | r00=c->GetSigmaY2(); | |
1087 | chi2 = (dy*dy)/r00; | |
46d29e70 | 1088 | } |
b8dc2353 | 1089 | else { |
6c94f330 | 1090 | Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12); |
1091 | // | |
1092 | r00=c->GetSigmaY2(); r01=0.; r11=c->GetSigmaZ2(); | |
1093 | r00+=fCyy; r01+=fCzy; r11+=fCzz; | |
b8dc2353 | 1094 | |
6c94f330 | 1095 | Double_t det=r00*r11 - r01*r01; |
b8dc2353 | 1096 | if (TMath::Abs(det) < 1.e-10) { |
6c94f330 | 1097 | Int_t n=GetNumberOfClusters(); |
1098 | if (n>4) cerr<<n<<" AliTRDtrack warning: Singular matrix !\n"; | |
b8dc2353 | 1099 | return 1e10; |
1100 | } | |
6c94f330 | 1101 | Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01; |
1102 | Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ; | |
4f1c04d3 | 1103 | Double_t tiltdz = dz; |
6c94f330 | 1104 | if (TMath::Abs(tiltdz)>padlength/2.) { |
1105 | tiltdz = TMath::Sign(padlength/2,dz); | |
4f1c04d3 | 1106 | } |
6c94f330 | 1107 | // dy=dy+h01*dz; |
1108 | dy=dy+h01*tiltdz; | |
a819a5f7 | 1109 | |
6c94f330 | 1110 | chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det; |
b8dc2353 | 1111 | } |
53c17fbf | 1112 | |
b8dc2353 | 1113 | return chi2; |
6c94f330 | 1114 | */ |
fe33d239 | 1115 | |
1116 | } | |
7ad19338 | 1117 | |
53c17fbf | 1118 | //_____________________________________________________________________________ |
16d9fbba | 1119 | void AliTRDtrack::MakeBackupTrack() |
1120 | { | |
1121 | // | |
53c17fbf | 1122 | // Creates a backup track |
16d9fbba | 1123 | // |
53c17fbf | 1124 | |
fe33d239 | 1125 | if (fBackupTrack) { |
1126 | delete fBackupTrack; | |
1127 | } | |
16d9fbba | 1128 | fBackupTrack = new AliTRDtrack(*this); |
4f1c04d3 | 1129 | |
1130 | } | |
1131 | ||
53c17fbf | 1132 | //_____________________________________________________________________________ |
1133 | Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z) | |
1134 | { | |
4f1c04d3 | 1135 | // |
fe33d239 | 1136 | // Find a prolongation at given x |
1137 | // Return 0 if it does not exist | |
1138 | // | |
1139 | ||
1140 | Double_t bz = GetBz(); | |
15ed8ba1 | 1141 | |
fe33d239 | 1142 | if (!AliExternalTrackParam::GetYAt(xk,bz,y)) { |
1143 | return 0; | |
1144 | } | |
1145 | if (!AliExternalTrackParam::GetZAt(xk,bz,z)) { | |
1146 | return 0; | |
1147 | } | |
4f1c04d3 | 1148 | |
6c94f330 | 1149 | return 1; |
fe33d239 | 1150 | |
16d9fbba | 1151 | } |
3fad3d32 | 1152 | |
53c17fbf | 1153 | //_____________________________________________________________________________ |
fe33d239 | 1154 | Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step) |
3fad3d32 | 1155 | { |
1156 | // | |
6c94f330 | 1157 | // Propagate track to given x position |
af26ce80 | 1158 | // Works inside of the 20 degree segmentation |
1159 | // (local cooordinate frame for TRD , TPC, TOF) | |
3fad3d32 | 1160 | // |
fe33d239 | 1161 | // Material budget from geo manager |
6c94f330 | 1162 | // |
fe33d239 | 1163 | |
1164 | Double_t xyz0[3]; | |
1165 | Double_t xyz1[3]; | |
1166 | Double_t y; | |
1167 | Double_t z; | |
1168 | ||
1169 | const Double_t kAlphac = TMath::Pi()/9.0; | |
15ed8ba1 | 1170 | const Double_t kTalphac = TMath::Tan(kAlphac*0.5); |
fe33d239 | 1171 | |
1172 | // Critical alpha - cross sector indication | |
af26ce80 | 1173 | Double_t dir = (GetX() > xr) ? -1.0 : 1.0; |
fe33d239 | 1174 | |
1175 | // Direction +- | |
1176 | for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) { | |
1177 | ||
6c94f330 | 1178 | GetXYZ(xyz0); |
3fad3d32 | 1179 | GetProlongation(x,y,z); |
fe33d239 | 1180 | xyz1[0] = x * TMath::Cos(GetAlpha()) + y * TMath::Sin(GetAlpha()); |
1181 | xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha()); | |
3fad3d32 | 1182 | xyz1[2] = z; |
1183 | Double_t param[7]; | |
826fe01c | 1184 | AliTracker::MeanMaterialBudget(xyz0,xyz1,param); |
fe33d239 | 1185 | |
1186 | if ((param[0] > 0) && | |
1187 | (param[1] > 0)) { | |
826fe01c | 1188 | PropagateTo(x,param[1],param[0]*param[4]); |
fe33d239 | 1189 | } |
1190 | ||
1191 | if (GetY() > GetX()*kTalphac) { | |
53c17fbf | 1192 | Rotate(-kAlphac); |
3fad3d32 | 1193 | } |
fe33d239 | 1194 | if (GetY() < -GetX()*kTalphac) { |
1195 | Rotate( kAlphac); | |
3fad3d32 | 1196 | } |
fe33d239 | 1197 | |
3fad3d32 | 1198 | } |
fe33d239 | 1199 | |
3fad3d32 | 1200 | PropagateTo(xr); |
53c17fbf | 1201 | |
3fad3d32 | 1202 | return 0; |
3fad3d32 | 1203 | |
53c17fbf | 1204 | } |
3fad3d32 | 1205 | |
53c17fbf | 1206 | //_____________________________________________________________________________ |
6c94f330 | 1207 | Int_t AliTRDtrack::PropagateToR(Double_t r,Double_t step) |
3fad3d32 | 1208 | { |
1209 | // | |
fe33d239 | 1210 | // Propagate track to the radial position |
1211 | // Rotation always connected to the last track position | |
3fad3d32 | 1212 | // |
fe33d239 | 1213 | |
1214 | Double_t xyz0[3]; | |
1215 | Double_t xyz1[3]; | |
1216 | Double_t y; | |
1217 | Double_t z; | |
1218 | ||
6c94f330 | 1219 | Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY()); |
fe33d239 | 1220 | // Direction +- |
af26ce80 | 1221 | Double_t dir = (radius > r) ? -1.0 : 1.0; |
fe33d239 | 1222 | |
1223 | for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) { | |
1224 | ||
6c94f330 | 1225 | GetXYZ(xyz0); |
3fad3d32 | 1226 | Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]); |
1227 | Rotate(alpha,kTRUE); | |
6c94f330 | 1228 | GetXYZ(xyz0); |
3fad3d32 | 1229 | GetProlongation(x,y,z); |
fe33d239 | 1230 | xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha); |
1231 | xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha); | |
3fad3d32 | 1232 | xyz1[2] = z; |
1233 | Double_t param[7]; | |
826fe01c | 1234 | AliTracker::MeanMaterialBudget(xyz0,xyz1,param); |
fe33d239 | 1235 | if (param[1] <= 0) { |
1236 | param[1] = 100000000; | |
1237 | } | |
826fe01c | 1238 | PropagateTo(x,param[1],param[0]*param[4]); |
fe33d239 | 1239 | |
3fad3d32 | 1240 | } |
fe33d239 | 1241 | |
6c94f330 | 1242 | GetXYZ(xyz0); |
3fad3d32 | 1243 | Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]); |
1244 | Rotate(alpha,kTRUE); | |
6c94f330 | 1245 | GetXYZ(xyz0); |
3fad3d32 | 1246 | GetProlongation(r,y,z); |
fe33d239 | 1247 | xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha); |
1248 | xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha); | |
3fad3d32 | 1249 | xyz1[2] = z; |
1250 | Double_t param[7]; | |
826fe01c | 1251 | AliTracker::MeanMaterialBudget(xyz0,xyz1,param); |
fe33d239 | 1252 | |
1253 | if (param[1] <= 0) { | |
1254 | param[1] = 100000000; | |
1255 | } | |
826fe01c | 1256 | PropagateTo(r,param[1],param[0]*param[4]); |
53c17fbf | 1257 | |
3fad3d32 | 1258 | return 0; |
53c17fbf | 1259 | |
1260 | } | |
1261 | ||
1262 | //_____________________________________________________________________________ | |
4e28c495 | 1263 | Int_t AliTRDtrack::GetSector() const |
53c17fbf | 1264 | { |
1265 | // | |
1266 | // Return the current sector | |
1267 | // | |
1268 | ||
fe33d239 | 1269 | return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha()) |
af26ce80 | 1270 | % AliTRDgeometry::kNsect; |
53c17fbf | 1271 | |
3fad3d32 | 1272 | } |
1273 | ||
53c17fbf | 1274 | //_____________________________________________________________________________ |
4e28c495 | 1275 | void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i) |
53c17fbf | 1276 | { |
1277 | // | |
1278 | // The sampled energy loss | |
1279 | // | |
1280 | ||
1281 | Double_t s = GetSnp(); | |
1282 | Double_t t = GetTgl(); | |
fe33d239 | 1283 | q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t)); |
53c17fbf | 1284 | fdQdl[i] = q; |
1285 | ||
1286 | } | |
1287 | ||
fe33d239 | 1288 | //_____________________________________________________________________________ |
4e28c495 | 1289 | void AliTRDtrack::SetSampledEdx(Float_t q) |
53c17fbf | 1290 | { |
1291 | // | |
1292 | // The sampled energy loss | |
1293 | // | |
1294 | ||
1295 | Double_t s = GetSnp(); | |
1296 | Double_t t = GetTgl(); | |
fe33d239 | 1297 | q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t)); |
53c17fbf | 1298 | fdQdl[fNdedx] = q; |
1299 | fNdedx++; | |
1300 | ||
1301 | } | |
1302 | ||
fe33d239 | 1303 | //_____________________________________________________________________________ |
1304 | Double_t AliTRDtrack::GetBz() const | |
1305 | { | |
15ed8ba1 | 1306 | // |
fe33d239 | 1307 | // Returns Bz component of the magnetic field (kG) |
15ed8ba1 | 1308 | // |
53c17fbf | 1309 | |
fe33d239 | 1310 | if (AliTracker::UniformField()) { |
1311 | return AliTracker::GetBz(); | |
1312 | } | |
1313 | Double_t r[3]; | |
1314 | GetXYZ(r); | |
1315 | ||
1316 | return AliTracker::GetBz(r); | |
53c17fbf | 1317 | |
fe33d239 | 1318 | } |