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