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