adding TPCFitMap to the AODtrack
[u/mrichter/AliRoot.git] / STEER / AOD / AliAODTrack.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2007, 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 //-------------------------------------------------------------------------
19 //     AOD track implementation of AliVTrack
20 //     Author: Markus Oldenburg, CERN
21 //     Markus.Oldenburg@cern.ch
22 //-------------------------------------------------------------------------
23
24 #include "AliLog.h"
25 #include "AliExternalTrackParam.h"
26 #include "AliVVertex.h"
27 #include "AliAODTrack.h"
28
29 ClassImp(AliAODTrack)
30
31 //______________________________________________________________________________
32 AliAODTrack::AliAODTrack() : 
33   AliVTrack(),
34   fRAtAbsorberEnd(0.),
35   fChi2perNDF(-999.),
36   fChi2MatchTrigger(0.),
37   fFlags(0),
38   fLabel(-999),
39   fITSMuonClusterMap(0),
40   fFilterMap(0),
41   fTPCFitMap(),
42   fTPCClusterMap(),
43   fTPCSharedMap(),
44   fTPCnclsF(0),
45   fID(-999),
46   fCharge(-99),
47   fType(kUndef),
48   fCaloIndex(kEMCALNoMatch),
49   fCovMatrix(NULL),
50   fDetPid(NULL),
51   fProdVertex(NULL)
52 {
53   // default constructor
54
55   SetP();
56   SetPosition((Float_t*)NULL);
57   SetXYAtDCA(-999., -999.);
58   SetPxPyPzAtDCA(-999., -999., -999.);
59   SetPID((Float_t*)NULL);
60 }
61
62 //______________________________________________________________________________
63 AliAODTrack::AliAODTrack(Short_t id,
64                          Int_t label, 
65                          Double_t p[3],
66                          Bool_t cartesian,
67                          Double_t x[3],
68                          Bool_t isDCA,
69                          Double_t covMatrix[21],
70                          Short_t charge,
71                          UChar_t itsClusMap,
72                          Double_t pid[10],
73                          AliAODVertex *prodVertex,
74                          Bool_t usedForVtxFit,
75                          Bool_t usedForPrimVtxFit,
76                          AODTrk_t ttype,
77                          UInt_t selectInfo,
78                          Float_t chi2perNDF) :
79   AliVTrack(),
80   fRAtAbsorberEnd(0.),
81   fChi2perNDF(chi2perNDF),
82   fChi2MatchTrigger(0.),
83   fFlags(0),
84   fLabel(label),
85   fITSMuonClusterMap(0),
86   fFilterMap(selectInfo),
87   fTPCFitMap(),
88   fTPCClusterMap(),
89   fTPCSharedMap(),
90   fTPCnclsF(0),
91   fID(id),
92   fCharge(charge),
93   fType(ttype),
94   fCaloIndex(kEMCALNoMatch),
95   fCovMatrix(NULL),
96   fDetPid(NULL),
97   fProdVertex(prodVertex)
98 {
99   // constructor
100  
101   SetP(p, cartesian);
102   SetPosition(x, isDCA);
103   SetXYAtDCA(-999., -999.);
104   SetPxPyPzAtDCA(-999., -999., -999.);
105   SetUsedForVtxFit(usedForVtxFit);
106   SetUsedForPrimVtxFit(usedForPrimVtxFit);
107   if(covMatrix) SetCovMatrix(covMatrix);
108   SetPID(pid);
109   SetITSClusterMap(itsClusMap);
110 }
111
112 //______________________________________________________________________________
113 AliAODTrack::AliAODTrack(Short_t id,
114                          Int_t label, 
115                          Float_t p[3],
116                          Bool_t cartesian,
117                          Float_t x[3],
118                          Bool_t isDCA,
119                          Float_t covMatrix[21],
120                          Short_t charge,
121                          UChar_t itsClusMap,
122                          Float_t pid[10],
123                          AliAODVertex *prodVertex,
124                          Bool_t usedForVtxFit,
125                          Bool_t usedForPrimVtxFit,
126                          AODTrk_t ttype,
127                          UInt_t selectInfo,
128                          Float_t chi2perNDF) :
129   AliVTrack(),
130   fRAtAbsorberEnd(0.),
131   fChi2perNDF(chi2perNDF),
132   fChi2MatchTrigger(0.),
133   fFlags(0),
134   fLabel(label),
135   fITSMuonClusterMap(0),
136   fFilterMap(selectInfo),
137   fTPCFitMap(),
138   fTPCClusterMap(),
139   fTPCSharedMap(),
140   fTPCnclsF(0),
141   fID(id),
142   fCharge(charge),
143   fType(ttype),
144   fCaloIndex(kEMCALNoMatch),
145   fCovMatrix(NULL),
146   fDetPid(NULL),
147   fProdVertex(prodVertex)
148 {
149   // constructor
150  
151   SetP(p, cartesian);
152   SetPosition(x, isDCA);
153   SetXYAtDCA(-999., -999.);
154   SetPxPyPzAtDCA(-999., -999., -999.);
155   SetUsedForVtxFit(usedForVtxFit);
156   SetUsedForPrimVtxFit(usedForPrimVtxFit);
157   if(covMatrix) SetCovMatrix(covMatrix);
158   SetPID(pid);
159   SetITSClusterMap(itsClusMap);
160 }
161
162 //______________________________________________________________________________
163 AliAODTrack::~AliAODTrack() 
164 {
165   // destructor
166   delete fCovMatrix;
167   delete fDetPid;
168 }
169
170
171 //______________________________________________________________________________
172 AliAODTrack::AliAODTrack(const AliAODTrack& trk) :
173   AliVTrack(trk),
174   fRAtAbsorberEnd(trk.fRAtAbsorberEnd),
175   fChi2perNDF(trk.fChi2perNDF),
176   fChi2MatchTrigger(trk.fChi2MatchTrigger),
177   fFlags(trk.fFlags),
178   fLabel(trk.fLabel),
179   fITSMuonClusterMap(trk.fITSMuonClusterMap),
180   fFilterMap(trk.fFilterMap),
181   fTPCFitMap(trk.fTPCFitMap),
182   fTPCClusterMap(trk.fTPCClusterMap),
183   fTPCSharedMap(trk.fTPCSharedMap),
184   fTPCnclsF(trk.fTPCnclsF),
185   fID(trk.fID),
186   fCharge(trk.fCharge),
187   fType(trk.fType),
188   fCaloIndex(trk.fCaloIndex),
189   fCovMatrix(NULL),
190   fDetPid(NULL),
191   fProdVertex(trk.fProdVertex)
192 {
193   // Copy constructor
194
195   trk.GetP(fMomentum);
196   trk.GetPosition(fPosition);
197   SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
198   SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
199   SetUsedForVtxFit(trk.GetUsedForVtxFit());
200   SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
201   if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
202   if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
203   SetPID(trk.fPID);
204 }
205
206 //______________________________________________________________________________
207 AliAODTrack& AliAODTrack::operator=(const AliAODTrack& trk)
208 {
209   // Assignment operator
210   if(this!=&trk) {
211
212     AliVTrack::operator=(trk);
213
214     trk.GetP(fMomentum);
215     trk.GetPosition(fPosition);
216     SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
217     SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
218     fRAtAbsorberEnd    = trk.fRAtAbsorberEnd;
219     fChi2perNDF        = trk.fChi2perNDF;
220     fChi2MatchTrigger  = trk.fChi2MatchTrigger;
221     trk.GetPID(fPID);
222     fFlags             = trk.fFlags;
223     fLabel             = trk.fLabel;    
224     fITSMuonClusterMap = trk.fITSMuonClusterMap;
225     fFilterMap         = trk.fFilterMap;
226     fTPCFitMap         = trk.fTPCFitMap;
227     fTPCClusterMap     = trk.fTPCClusterMap;
228     fTPCSharedMap      = trk.fTPCSharedMap;
229     fTPCnclsF          = trk.fTPCnclsF;
230     fID                = trk.fID;
231     fCharge            = trk.fCharge;
232     fType              = trk.fType;
233     fCaloIndex         = trk.fCaloIndex;
234
235     delete fCovMatrix;
236     if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
237     else fCovMatrix=NULL;
238
239
240     fProdVertex        = trk.fProdVertex;
241     SetUsedForVtxFit(trk.GetUsedForVtxFit());
242     SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
243
244     delete fDetPid;
245     if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
246     else fDetPid=NULL;
247
248
249   }
250
251   return *this;
252 }
253
254 //______________________________________________________________________________
255 Double_t AliAODTrack::M(AODTrkPID_t pid) const
256 {
257   // Returns the mass.
258   // Masses for nuclei don't exist in the PDG tables, therefore they were put by hand.
259
260   switch (pid) {
261
262   case kElectron :
263     return 0.000510999; //TDatabasePDG::Instance()->GetParticle(11/*::kElectron*/)->Mass();
264     break;
265
266   case kMuon :
267     return 0.1056584; //TDatabasePDG::Instance()->GetParticle(13/*::kMuonMinus*/)->Mass();
268     break;
269
270   case kPion :
271     return 0.13957; //TDatabasePDG::Instance()->GetParticle(211/*::kPiPlus*/)->Mass();
272     break;
273
274   case kKaon :
275     return 0.4937; //TDatabasePDG::Instance()->GetParticle(321/*::kKPlus*/)->Mass();
276     break;
277
278   case kProton :
279     return 0.9382720; //TDatabasePDG::Instance()->GetParticle(2212/*::kProton*/)->Mass();
280     break;
281
282   case kDeuteron :
283     return 1.8756; //TDatabasePDG::Instance()->GetParticle(1000010020)->Mass();
284     break;
285
286   case kTriton :
287     return 2.8089; //TDatabasePDG::Instance()->GetParticle(1000010030)->Mass();
288     break;
289
290   case kHelium3 :
291     return 2.8084; //TDatabasePDG::Instance()->GetParticle(1000020030)->Mass();
292     break;
293
294   case kAlpha :
295     return 3.7274; //TDatabasePDG::Instance()->GetParticle(1000020040)->Mass();
296     break;
297
298   case kUnknown :
299     return -999.;
300     break;
301
302   default :
303     return -999.;
304   }
305 }
306
307 //______________________________________________________________________________
308 Double_t AliAODTrack::E(AODTrkPID_t pid) const
309 {
310   // Returns the energy of the particle of a given pid.
311   
312   if (pid != kUnknown) { // particle was identified
313     Double_t m = M(pid);
314     return TMath::Sqrt(P()*P() + m*m);
315   } else { // pid unknown
316     return -999.;
317   }
318 }
319
320 //______________________________________________________________________________
321 Double_t AliAODTrack::Y(AODTrkPID_t pid) const
322 {
323   // Returns the rapidity of a particle of a given pid.
324   
325   if (pid != kUnknown) { // particle was identified
326     Double_t e = E(pid);
327     Double_t pz = Pz();
328     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
329       return 0.5*TMath::Log((e+pz)/(e-pz));
330     } else { // energy not known or equal to pz
331       return -999.;
332     }
333   } else { // pid unknown
334     return -999.;
335   }
336 }
337
338 //______________________________________________________________________________
339 Double_t AliAODTrack::Y(Double_t m) const
340 {
341   // Returns the rapidity of a particle of a given mass.
342   
343   if (m >= 0.) { // mass makes sense
344     Double_t e = E(m);
345     Double_t pz = Pz();
346     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
347       return 0.5*TMath::Log((e+pz)/(e-pz));
348     } else { // energy not known or equal to pz
349       return -999.;
350     }
351   } else { // pid unknown
352     return -999.;
353   }
354 }
355
356 //______________________________________________________________________________
357 AliAODTrack::AODTrkPID_t AliAODTrack::GetMostProbablePID() const 
358 {
359   // Returns the most probable PID array element.
360   
361   Int_t nPID = 10;
362   AODTrkPID_t loc = kUnknown;
363   Double_t max = 0.;
364   Bool_t allTheSame = kTRUE;
365   
366   for (Int_t iPID = 0; iPID < nPID; iPID++) {
367     if (fPID[iPID] >= max) {
368       if (fPID[iPID] > max) {
369         allTheSame = kFALSE;
370         max = fPID[iPID];
371         loc = (AODTrkPID_t)iPID;
372       } else {
373         allTheSame = kTRUE;
374       }
375     }
376   }
377   return allTheSame ? kUnknown : loc;
378 }
379
380 //______________________________________________________________________________
381 void AliAODTrack::ConvertAliPIDtoAODPID()
382 {
383   // Converts AliPID array.
384   // The numbering scheme is the same for electrons, muons, pions, kaons, and protons.
385   // Everything else has to be set to zero.
386
387   fPID[kDeuteron] = 0.;
388   fPID[kTriton]   = 0.;
389   fPID[kHelium3]  = 0.;
390   fPID[kAlpha]    = 0.;
391   fPID[kUnknown]  = 0.;
392   
393   return;
394 }
395
396
397 //______________________________________________________________________________
398 template <class T> void AliAODTrack::SetP(const T *p, const Bool_t cartesian) 
399 {
400   // Set the momentum
401
402   if (p) {
403     if (cartesian) {
404       Double_t pt2 = p[0]*p[0] + p[1]*p[1];
405       Double_t pp  = TMath::Sqrt(pt2 + p[2]*p[2]);
406       
407       fMomentum[0] = TMath::Sqrt(pt2); // pt
408       fMomentum[1] = (pt2 != 0.) ? TMath::Pi()+TMath::ATan2(-p[1], -p[0]) : -999; // phi
409       fMomentum[2] = (pp != 0.) ? TMath::ACos(p[2] / pp) : -999.; // theta
410     } else {
411       fMomentum[0] = p[0];  // pt
412       fMomentum[1] = p[1];  // phi
413       fMomentum[2] = p[2];  // theta
414     }
415   } else {
416     fMomentum[0] = -999.;
417     fMomentum[1] = -999.;
418     fMomentum[2] = -999.;
419   }
420 }
421
422 //______________________________________________________________________________
423 template <class T> void AliAODTrack::SetPosition(const T *x, const Bool_t dca) 
424 {
425   // set the position
426
427   if (x) {
428     if (!dca) {
429       ResetBit(kIsDCA);
430
431       fPosition[0] = x[0];
432       fPosition[1] = x[1];
433       fPosition[2] = x[2];
434     } else {
435       SetBit(kIsDCA);
436       // don't know any better yet
437       fPosition[0] = -999.;
438       fPosition[1] = -999.;
439       fPosition[2] = -999.;
440     }
441   } else {
442     ResetBit(kIsDCA);
443
444     fPosition[0] = -999.;
445     fPosition[1] = -999.;
446     fPosition[2] = -999.;
447   }
448 }
449
450 //______________________________________________________________________________
451 void AliAODTrack::SetDCA(Double_t d, Double_t z) 
452 {
453   // set the dca
454   fPosition[0] = d;
455   fPosition[1] = z;
456   fPosition[2] = 0.;
457   SetBit(kIsDCA);
458 }
459
460 //______________________________________________________________________________
461 void AliAODTrack::Print(Option_t* /* option */) const
462 {
463   // prints information about AliAODTrack
464
465   printf("Object name: %s   Track type: %s\n", GetName(), GetTitle()); 
466   printf("        px = %f\n", Px());
467   printf("        py = %f\n", Py());
468   printf("        pz = %f\n", Pz());
469   printf("        pt = %f\n", Pt());
470   printf("      1/pt = %f\n", OneOverPt());
471   printf("     theta = %f\n", Theta());
472   printf("       phi = %f\n", Phi());
473   printf("  chi2/NDF = %f\n", Chi2perNDF());
474   printf("    charge = %d\n", Charge());
475 }
476
477 //______________________________________________________________________________
478 void AliAODTrack::SetMatchTrigger(Int_t matchTrig)
479 {
480   // Set the MUON trigger information
481   switch(matchTrig){
482     case 0: // 0 track does not match trigger
483       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
484       break;
485     case 1: // 1 track match but does not pass pt cut
486       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x40000000;
487       break;
488     case 2: // 2 track match Low pt cut
489       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x80000000;
490       break;
491     case 3: // 3 track match High pt cut
492       fITSMuonClusterMap=fITSMuonClusterMap|0xc0000000;
493       break;
494     default:
495       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
496       AliWarning(Form("unknown case for matchTrig: %d\n",matchTrig));
497   }
498 }
499
500 //______________________________________________________________________________
501 Bool_t AliAODTrack::HitsMuonChamber(Int_t MuonChamber, Int_t cathode) const
502 {
503   // return kTRUE if the track fires the given tracking or trigger chamber.
504   // If the chamber is a trigger one:
505   // - if cathode = 0 or 1, the track matches the corresponding cathode
506   // - if cathode = -1, the track matches both cathodes
507   
508   if (MuonChamber < 0) return kFALSE;
509   
510   if (MuonChamber < 10) return TESTBIT(GetMUONClusterMap(), MuonChamber);
511   
512   if (MuonChamber < 14) {
513     
514     if (cathode < 0) return TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber) &&
515                             TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber+4);
516     
517     if (cathode < 2) return TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber+(1-cathode)*4);
518     
519   }
520   
521   return kFALSE;
522 }
523
524 //______________________________________________________________________________
525 Bool_t AliAODTrack::MatchTriggerDigits() const
526 {
527   // return kTRUE if the track matches a digit on both planes of at least 2 trigger chambers
528   
529   Int_t nMatchedChambers = 0;
530   for (Int_t ich=10; ich<14; ich++) if (HitsMuonChamber(ich)) nMatchedChambers++;
531   
532   return (nMatchedChambers >= 2);
533 }
534
535 //______________________________________________________________________________
536 Bool_t AliAODTrack::PropagateToDCA(const AliVVertex *vtx, 
537     Double_t b, Double_t maxd, Double_t dz[2], Double_t covar[3])
538 {
539   // compute impact parameters to the vertex vtx and their covariance matrix
540   // b is the Bz, needed to propagate correctly the track to vertex 
541   // only the track parameters are update after the propagation (pos and mom),
542   // not the covariance matrix. This is OK for propagation over short distance
543   // inside the beam pipe.
544   // return kFALSE is something went wrong
545
546   // convert to AliExternalTrackParam
547   AliExternalTrackParam etp; etp.CopyFromVTrack(this);  
548
549   Float_t xstart = etp.GetX();
550   if(xstart>3.) {
551     AliError("This method can be used only for propagation inside the beam pipe");
552     return kFALSE; 
553   }
554
555   if(!etp.PropagateToDCA(vtx,b,maxd,dz,covar)) return kFALSE;
556
557   // update track position and momentum
558   Double_t mom[3];
559   etp.GetPxPyPz(mom);
560   SetP(mom,kTRUE);
561   etp.GetXYZ(mom);
562   SetPosition(mom,kFALSE);
563
564
565   return kTRUE;
566 }
567
568 //______________________________________________________________________________
569 Bool_t AliAODTrack::GetPxPyPz(Double_t p[3]) const 
570 {
571     //---------------------------------------------------------------------
572     // This function returns the global track momentum components
573     //---------------------------------------------------------------------
574   p[0]=Px(); p[1]=Py(); p[2]=Pz();
575   return kTRUE;
576 }
577
578 //______________________________________________________________________________
579 Float_t AliAODTrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
580 {
581   //
582   // TPC cluster information
583   // type 0: get fraction of found/findable clusters with neighbourhood definition
584   //      1: findable clusters with neighbourhood definition
585   //      2: found clusters
586   //
587   // definition of findable clusters:
588   //            a cluster is defined as findable if there is another cluster
589   //           within +- nNeighbours pad rows. The idea is to overcome threshold
590   //           effects with a very simple algorithm.
591   //
592   
593   if (type==2) return fTPCClusterMap.CountBits();
594   
595   Int_t found=0;
596   Int_t findable=0;
597   Int_t last=-nNeighbours;
598   
599   for (Int_t i=row0; i<row1; ++i){
600     //look to current row
601     if (fTPCClusterMap[i]) {
602       last=i;
603       ++found;
604       ++findable;
605       continue;
606     }
607     //look to nNeighbours before
608     if ((i-last)<=nNeighbours) {
609       ++findable;
610       continue;
611     }
612     //look to nNeighbours after
613     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
614       if (fTPCClusterMap[j]){
615         ++findable;
616         break;
617       }
618     }
619   }
620   if (type==1) return findable;
621   
622   if (type==0){
623     Float_t fraction=0;
624     if (findable>0)
625       fraction=(Float_t)found/(Float_t)findable;
626     else
627       fraction=0;
628     return fraction;
629   }
630   return 0;  // undefined type - default value
631 }
632
633 //______________________________________________________________________________
634 Double_t  AliAODTrack::GetTRDslice(Int_t plane, Int_t slice) const {
635   //
636   // return TRD Pid information
637   //
638   if (!fDetPid) return -1;
639   Double32_t *trdSlices=fDetPid->GetTRDsignal();
640   if (!trdSlices) return -1;
641   if ((plane<0) || (plane>=kTRDnPlanes)) {
642     return -1.;
643   }
644
645   Int_t ns=fDetPid->GetTRDnSlices()/kTRDnPlanes;
646   if ((slice<-1) || (slice>=ns)) {
647     return -1.;
648   }
649
650   if(slice>=0) return trdSlices[plane*ns + slice];
651
652   // return average of the dEdx measurements
653   Double_t q=0.; Double32_t *s = &trdSlices[plane*ns];
654   for (Int_t i=0; i<ns; i++, s++) if((*s)>0.) q+=(*s);
655   return q/ns;
656 }
657
658 //______________________________________________________________________________
659 UChar_t AliAODTrack::GetTRDntrackletsPID() const{
660   //
661   // return number of tracklets calculated from the slices
662   //
663   if(!fDetPid) return -1;
664   return fDetPid->GetTRDntrackletsPID();
665 }
666
667 //______________________________________________________________________________
668 UChar_t AliAODTrack::GetTRDncls(Int_t layer) const {
669   // 
670   // return number of TRD clusters
671   //
672   if(!fDetPid || layer > 5) return -1;
673   if(layer < 0) return fDetPid->GetTRDncls();
674   else return fDetPid->GetTRDncls(layer);
675 }
676
677 //______________________________________________________________________________
678 Double_t AliAODTrack::GetTRDmomentum(Int_t plane, Double_t */*sp*/) const
679 {
680   //Returns momentum estimation
681   // in TRD layer "plane".
682
683   if (!fDetPid) return -1;
684   const Float_t *trdMomentum=fDetPid->GetTRDmomentum();
685
686   if (!trdMomentum) {
687     return -1.;
688   }
689   if ((plane<0) || (plane>=kTRDnPlanes)) {
690     return -1.;
691   }
692
693   return trdMomentum[plane];
694 }
695
696 //_______________________________________________________________________
697 Int_t AliAODTrack::GetTOFBunchCrossing(Double_t b) const 
698 {
699   // Returns the number of bunch crossings after trigger (assuming 25ns spacing)
700   const double kSpacing = 25e3; // min interbanch spacing
701   const double kShift = 0;
702   Int_t bcid = kTOFBCNA; // defualt one
703   if (!IsOn(kTOFout) || !IsOn(kESDpid)) return bcid; // no info
704   //
705   double tdif = GetTOFsignal();
706   if (IsOn(kTIME)) { // integrated time info is there
707     int pid = (int)GetMostProbablePID();
708     double ttimes[10]; 
709     GetIntegratedTimes(ttimes);
710     tdif -= ttimes[pid];
711   }
712   else { // assume integrated time info from TOF radius and momentum
713     const double kRTOF = 385.;
714     const double kCSpeed = 3.e-2; // cm/ps
715     double p = P();
716     if (p<0.001) p = 1.0;
717     double m = M();
718     double path =  kRTOF;     // mean TOF radius
719     if (TMath::Abs(b)>kAlmost0) {  // account for curvature
720       double curv = Pt()/(b*kB2C);
721       if (curv>kAlmost0) {
722         double tgl = Pz()/Pt();
723         path = 2./curv*TMath::ASin(kRTOF*curv/2.)*TMath::Sqrt(1.+tgl*tgl);
724       }
725     }
726     tdif -= path/kCSpeed*TMath::Sqrt(1.+m*m/(p*p));
727   }
728   bcid = TMath::Nint((tdif - kShift)/kSpacing);
729   return bcid;
730 }