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