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