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