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