5ea54212fb38bab7c21fc010e479a393f334c7c6
[u/mrichter/AliRoot.git] / TRD / AliTRDtrack.cxx
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
28 ClassImp(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 //_____________________________________________________________________________
39 AliTRDtrack::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 //_____________________________________________________________________________
91 AliTRDtrack::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 //_____________________________________________________________________________
178 AliTRDtrack::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 //_____________________________________________________________________________
244 AliTRDtrack::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 //_____________________________________________________________________________
301 AliTRDtrack::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 //____________________________________________________________________________
385 AliTRDtrack::~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 //____________________________________________________________________________
408 Float_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 //_____________________________________________________________________________
423 Int_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 //_____________________________________________________________________________
446 void 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((int)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 //_____________________________________________________________________________
482 void 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 //_____________________________________________________________________________
558 void 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 //_____________________________________________________________________________
612 void 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 //_____________________________________________________________________________
633 Float_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 //_____________________________________________________________________________
651 Bool_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 //_____________________________________________________________________________
754 Bool_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 //_____________________________________________________________________________
843 Bool_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 //_____________________________________________________________________________
877 Int_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 //_____________________________________________________________________________
916 Bool_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 //_____________________________________________________________________________
946 Bool_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 //_____________________________________________________________________________
966 Double_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 //_____________________________________________________________________________
985 void 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 //_____________________________________________________________________________
999 Int_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 //_____________________________________________________________________________
1020 Int_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 //_____________________________________________________________________________
1073 Int_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 //_____________________________________________________________________________
1129 Int_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 //_____________________________________________________________________________
1141 void 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 //_____________________________________________________________________________
1155 void 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 }