Preprocessor for AD
[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 #include "AliTOFHeader.h"
32
33 #include "AliAODTrack.h"
34
35 ClassImp(AliAODTrack)
36
37 //______________________________________________________________________________
38 AliAODTrack::AliAODTrack() : 
39   AliVTrack(),
40   fRAtAbsorberEnd(0.),
41   fChi2perNDF(-999.),
42   fChi2MatchTrigger(0.),
43   fPID(0),
44   fFlags(0),
45   fLabel(-999),
46   fTOFLabel(),
47   fTrackLength(0),
48   fITSMuonClusterMap(0),
49   fMUONtrigHitsMapTrg(0),
50   fMUONtrigHitsMapTrk(0),
51   fFilterMap(0),
52   fTPCFitMap(),
53   fTPCClusterMap(),
54   fTPCSharedMap(),
55   fTPCnclsF(0),
56   fTPCNCrossedRows(0),
57   fID(-999),
58   fCharge(-99),
59   fType(kUndef),
60   fPIDForTracking(AliPID::kPion),
61   fCaloIndex(kEMCALNoMatch),
62   fCovMatrix(NULL),
63   fDetPid(NULL),
64   fDetectorPID(NULL),
65   fProdVertex(NULL),
66   fTrackPhiOnEMCal(-999),
67   fTrackEtaOnEMCal(-999),
68   fTrackPtOnEMCal(-999),
69   fIsMuonGlobalTrack(kFALSE),    // AU
70   fTPCsignalTuned(0),
71   fTOFsignalTuned(99999),
72   fMFTClusterPattern(0),         // AU
73   fAODEvent(NULL)
74 {
75   // default constructor
76
77   SetP();
78   SetPosition((Float_t*)NULL);
79   SetXYAtDCA(-999., -999.);
80   SetPxPyPzAtDCA(-999., -999., -999.);
81   for (Int_t i = 0; i < 3; i++) {fTOFLabel[i] = -1;}
82 }
83
84 //______________________________________________________________________________
85 AliAODTrack::AliAODTrack(Short_t id,
86                          Int_t label, 
87                          Double_t p[3],
88                          Bool_t cartesian,
89                          Double_t x[3],
90                          Bool_t isDCA,
91                          Double_t covMatrix[21],
92                          Short_t charge,
93                          UChar_t itsClusMap,
94                          AliAODVertex *prodVertex,
95                          Bool_t usedForVtxFit,
96                          Bool_t usedForPrimVtxFit,
97                          AODTrk_t ttype,
98                          UInt_t selectInfo,
99                          Float_t chi2perNDF) :
100   AliVTrack(),
101   fRAtAbsorberEnd(0.),
102   fChi2perNDF(chi2perNDF),
103   fChi2MatchTrigger(0.),
104   fPID(0),
105   fFlags(0),
106   fLabel(label),
107   fTOFLabel(),
108   fTrackLength(0),
109   fITSMuonClusterMap(0),
110   fMUONtrigHitsMapTrg(0),
111   fMUONtrigHitsMapTrk(0),
112   fFilterMap(selectInfo),
113   fTPCFitMap(),
114   fTPCClusterMap(),
115   fTPCSharedMap(),
116   fTPCnclsF(0),
117   fTPCNCrossedRows(0),
118   fID(id),
119   fCharge(charge),
120   fType(ttype),
121   fPIDForTracking(AliPID::kPion),
122   fCaloIndex(kEMCALNoMatch),
123   fCovMatrix(NULL),
124   fDetPid(NULL),
125   fDetectorPID(NULL),
126   fProdVertex(prodVertex),
127   fTrackPhiOnEMCal(-999),
128   fTrackEtaOnEMCal(-999),
129   fTrackPtOnEMCal(-999),
130   fIsMuonGlobalTrack(kFALSE),    // AU
131   fTPCsignalTuned(0),
132   fTOFsignalTuned(99999),
133   fMFTClusterPattern(0),         // AU
134   fAODEvent(NULL)
135 {
136   // constructor
137  
138   SetP(p, cartesian);
139   SetPosition(x, isDCA);
140   SetXYAtDCA(-999., -999.);
141   SetPxPyPzAtDCA(-999., -999., -999.);
142   SetUsedForVtxFit(usedForVtxFit);
143   SetUsedForPrimVtxFit(usedForPrimVtxFit);
144   if(covMatrix) SetCovMatrix(covMatrix);
145   SetITSClusterMap(itsClusMap);
146   for (Int_t i=0;i<3;i++) {fTOFLabel[i]=-1;}
147 }
148
149 //______________________________________________________________________________
150 AliAODTrack::AliAODTrack(Short_t id,
151                          Int_t label, 
152                          Float_t p[3],
153                          Bool_t cartesian,
154                          Float_t x[3],
155                          Bool_t isDCA,
156                          Float_t covMatrix[21],
157                          Short_t charge,
158                          UChar_t itsClusMap,
159                          AliAODVertex *prodVertex,
160                          Bool_t usedForVtxFit,
161                          Bool_t usedForPrimVtxFit,
162                          AODTrk_t ttype,
163                          UInt_t selectInfo,
164                          Float_t chi2perNDF ) :
165   AliVTrack(),
166   fRAtAbsorberEnd(0.),
167   fChi2perNDF(chi2perNDF),
168   fChi2MatchTrigger(0.),
169   fPID(0),
170   fFlags(0),
171   fLabel(label),
172   fTOFLabel(),
173   fTrackLength(0),
174   fITSMuonClusterMap(0),
175   fMUONtrigHitsMapTrg(0),
176   fMUONtrigHitsMapTrk(0),
177   fFilterMap(selectInfo),
178   fTPCFitMap(),
179   fTPCClusterMap(),
180   fTPCSharedMap(),
181   fTPCnclsF(0),
182   fTPCNCrossedRows(0),
183   fID(id),
184   fCharge(charge),
185   fType(ttype),
186   fPIDForTracking(AliPID::kPion),
187   fCaloIndex(kEMCALNoMatch),
188   fCovMatrix(NULL),
189   fDetPid(NULL),
190   fDetectorPID(NULL),
191   fProdVertex(prodVertex),
192   fTrackPhiOnEMCal(-999),
193   fTrackEtaOnEMCal(-999),
194   fTrackPtOnEMCal(-999),
195   fIsMuonGlobalTrack(kFALSE),    // AU
196   fTPCsignalTuned(0),
197   fTOFsignalTuned(99999),
198   fMFTClusterPattern(0),         // AU
199   fAODEvent(NULL)
200 {
201   // constructor
202  
203   SetP(p, cartesian);
204   SetPosition(x, isDCA);
205   SetXYAtDCA(-999., -999.);
206   SetPxPyPzAtDCA(-999., -999., -999.);
207   SetUsedForVtxFit(usedForVtxFit);
208   SetUsedForPrimVtxFit(usedForPrimVtxFit);
209   if(covMatrix) SetCovMatrix(covMatrix);
210   SetITSClusterMap(itsClusMap);
211   for (Int_t i=0;i<3;i++) {fTOFLabel[i]=-1;}
212 }
213
214 //______________________________________________________________________________
215 AliAODTrack::~AliAODTrack() 
216 {
217   // destructor
218   delete fCovMatrix;
219   delete fDetPid;
220   delete fDetectorPID;
221   if (fPID) {delete[] fPID; fPID = 0;}
222 }
223
224
225 //______________________________________________________________________________
226 AliAODTrack::AliAODTrack(const AliAODTrack& trk) :
227   AliVTrack(trk),
228   fRAtAbsorberEnd(trk.fRAtAbsorberEnd),
229   fChi2perNDF(trk.fChi2perNDF),
230   fChi2MatchTrigger(trk.fChi2MatchTrigger),
231   fPID(0),
232   fFlags(trk.fFlags),
233   fLabel(trk.fLabel),
234   fTOFLabel(),
235   fTrackLength(trk.fTrackLength),
236   fITSMuonClusterMap(trk.fITSMuonClusterMap),
237   fMUONtrigHitsMapTrg(trk.fMUONtrigHitsMapTrg),
238   fMUONtrigHitsMapTrk(trk.fMUONtrigHitsMapTrk),
239   fFilterMap(trk.fFilterMap),
240   fTPCFitMap(trk.fTPCFitMap),
241   fTPCClusterMap(trk.fTPCClusterMap),
242   fTPCSharedMap(trk.fTPCSharedMap),
243   fTPCnclsF(trk.fTPCnclsF),
244   fTPCNCrossedRows(trk.fTPCNCrossedRows),
245   fID(trk.fID),
246   fCharge(trk.fCharge),
247   fType(trk.fType),
248   fPIDForTracking(trk.fPIDForTracking),
249   fCaloIndex(trk.fCaloIndex),
250   fCovMatrix(NULL),
251   fDetPid(NULL),
252   fDetectorPID(NULL),
253   fProdVertex(trk.fProdVertex),
254   fTrackPhiOnEMCal(trk.fTrackPhiOnEMCal),
255   fTrackEtaOnEMCal(trk.fTrackEtaOnEMCal),
256   fTrackPtOnEMCal(trk.fTrackPtOnEMCal),
257   fIsMuonGlobalTrack(trk.fIsMuonGlobalTrack),    // AU
258   fTPCsignalTuned(trk.fTPCsignalTuned),
259   fTOFsignalTuned(trk.fTOFsignalTuned),
260   fMFTClusterPattern(trk.fMFTClusterPattern),    // AU
261   fAODEvent(trk.fAODEvent)
262 {
263   // Copy constructor
264
265   trk.GetP(fMomentum);
266   trk.GetPosition(fPosition);
267   SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
268   SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
269   SetUsedForVtxFit(trk.GetUsedForVtxFit());
270   SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
271   if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
272   if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
273   SetPID(trk.fPID);
274   if (trk.fDetectorPID) fDetectorPID = new AliDetectorPID(*trk.fDetectorPID);
275   for (Int_t i = 0; i < 3; i++) {fTOFLabel[i] = trk.fTOFLabel[i];}  
276 }
277
278 //______________________________________________________________________________
279 AliAODTrack& AliAODTrack::operator=(const AliAODTrack& trk)
280 {
281   // Assignment operator
282   if(this!=&trk) {
283
284     AliVTrack::operator=(trk);
285
286     trk.GetP(fMomentum);
287     trk.GetPosition(fPosition);
288     SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
289     SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
290     fRAtAbsorberEnd    = trk.fRAtAbsorberEnd;
291     fChi2perNDF        = trk.fChi2perNDF;
292     fChi2MatchTrigger  = trk.fChi2MatchTrigger;
293     SetPID( trk.fPID );
294     fFlags             = trk.fFlags;
295     fLabel             = trk.fLabel;    
296     fTrackLength       = trk.fTrackLength;
297     fITSMuonClusterMap = trk.fITSMuonClusterMap;
298     fMUONtrigHitsMapTrg = trk.fMUONtrigHitsMapTrg;
299     fMUONtrigHitsMapTrk = trk.fMUONtrigHitsMapTrk;
300     fFilterMap         = trk.fFilterMap;
301     fTPCFitMap         = trk.fTPCFitMap;
302     fTPCClusterMap     = trk.fTPCClusterMap;
303     fTPCSharedMap      = trk.fTPCSharedMap;
304     fTPCnclsF          = trk.fTPCnclsF;
305     fTPCNCrossedRows   = trk.fTPCNCrossedRows;
306     fID                = trk.fID;
307     fCharge            = trk.fCharge;
308     fType              = trk.fType;
309     fPIDForTracking    = trk.fPIDForTracking;
310     fCaloIndex         = trk.fCaloIndex;
311     fTrackPhiOnEMCal   = trk.fTrackPhiOnEMCal;
312     fTrackEtaOnEMCal   = trk.fTrackEtaOnEMCal;
313     fTrackPtOnEMCal    = trk.fTrackPtOnEMCal;
314     fIsMuonGlobalTrack = trk.fIsMuonGlobalTrack;     // AU
315     fTPCsignalTuned    = trk.fTPCsignalTuned;
316     fTOFsignalTuned    = trk.fTOFsignalTuned;
317     fMFTClusterPattern = trk.fMFTClusterPattern;     // AU
318     
319     delete fCovMatrix;
320     if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
321     else fCovMatrix=NULL;
322
323
324     fProdVertex        = trk.fProdVertex;
325     SetUsedForVtxFit(trk.GetUsedForVtxFit());
326     SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
327
328     //detector raw signals
329     delete fDetPid;
330     if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
331     else fDetPid=NULL;
332
333     //calibrated PID cache
334     delete fDetectorPID;
335     fDetectorPID=0x0;
336     if (trk.fDetectorPID) fDetectorPID = new AliDetectorPID(*trk.fDetectorPID);
337     for (Int_t i = 0; i < 3; i++) {fTOFLabel[i] = trk.fTOFLabel[i];}  
338   }
339
340   return *this;
341 }
342
343 //______________________________________________________________________________
344 Double_t AliAODTrack::M(AODTrkPID_t pid) const
345 {
346   // Returns the mass.
347   // Masses for nuclei don't exist in the PDG tables, therefore they were put by hand.
348
349   switch (pid) {
350
351   case kElectron :
352     return 0.000510999; //TDatabasePDG::Instance()->GetParticle(11/*::kElectron*/)->Mass();
353     break;
354
355   case kMuon :
356     return 0.1056584; //TDatabasePDG::Instance()->GetParticle(13/*::kMuonMinus*/)->Mass();
357     break;
358
359   case kPion :
360     return 0.13957; //TDatabasePDG::Instance()->GetParticle(211/*::kPiPlus*/)->Mass();
361     break;
362
363   case kKaon :
364     return 0.4937; //TDatabasePDG::Instance()->GetParticle(321/*::kKPlus*/)->Mass();
365     break;
366
367   case kProton :
368     return 0.9382720; //TDatabasePDG::Instance()->GetParticle(2212/*::kProton*/)->Mass();
369     break;
370
371   case kDeuteron :
372     return 1.8756; //TDatabasePDG::Instance()->GetParticle(1000010020)->Mass();
373     break;
374
375   case kTriton :
376     return 2.8089; //TDatabasePDG::Instance()->GetParticle(1000010030)->Mass();
377     break;
378
379   case kHelium3 :
380     return 2.8084; //TDatabasePDG::Instance()->GetParticle(1000020030)->Mass();
381     break;
382
383   case kAlpha :
384     return 3.7274; //TDatabasePDG::Instance()->GetParticle(1000020040)->Mass();
385     break;
386
387   case kUnknown :
388     return -999.;
389     break;
390
391   default :
392     return -999.;
393   }
394 }
395
396 //______________________________________________________________________________
397 Double_t AliAODTrack::E(AODTrkPID_t pid) const
398 {
399   // Returns the energy of the particle of a given pid.
400   
401   if (pid != kUnknown) { // particle was identified
402     Double_t m = M(pid);
403     return TMath::Sqrt(P()*P() + m*m);
404   } else { // pid unknown
405     return -999.;
406   }
407 }
408
409 //______________________________________________________________________________
410 Double_t AliAODTrack::Y(AODTrkPID_t pid) const
411 {
412   // Returns the rapidity of a particle of a given pid.
413   
414   if (pid != kUnknown) { // particle was identified
415     Double_t e = E(pid);
416     Double_t pz = Pz();
417     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
418       return 0.5*TMath::Log((e+pz)/(e-pz));
419     } else { // energy not known or equal to pz
420       return -999.;
421     }
422   } else { // pid unknown
423     return -999.;
424   }
425 }
426
427 //______________________________________________________________________________
428 Double_t AliAODTrack::Y(Double_t m) const
429 {
430   // Returns the rapidity of a particle of a given mass.
431   
432   if (m >= 0.) { // mass makes sense
433     Double_t e = E(m);
434     Double_t pz = Pz();
435     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
436       return 0.5*TMath::Log((e+pz)/(e-pz));
437     } else { // energy not known or equal to pz
438       return -999.;
439     }
440   } else { // pid unknown
441     return -999.;
442   }
443 }
444
445 void AliAODTrack::SetTOFLabel(const Int_t *p) {  
446   // Sets  (in TOF)
447   for (Int_t i = 0; i < 3; i++) fTOFLabel[i]=p[i];
448 }
449
450 //_______________________________________________________________________
451 void AliAODTrack::GetTOFLabel(Int_t *p) const {
452   // Gets (in TOF)
453   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
454 }
455
456 //______________________________________________________________________________
457 AliAODTrack::AODTrkPID_t AliAODTrack::GetMostProbablePID() const 
458 {
459   // Returns the most probable PID array element.
460   
461   Int_t nPID = 10;
462   AODTrkPID_t loc = kUnknown;
463   Bool_t allTheSame = kTRUE;
464   if (fPID) {
465     Double_t max = 0.;
466     for (Int_t iPID = 0; iPID < nPID; iPID++) {
467       if (fPID[iPID] >= max) {
468         if (fPID[iPID] > max) {
469           allTheSame = kFALSE;
470           max = fPID[iPID];
471           loc = (AODTrkPID_t)iPID;
472         } else {
473           allTheSame = kTRUE;
474         }
475       }
476     }
477   }
478   return allTheSame ? AODTrkPID_t(GetPIDForTracking()) : loc;
479 }
480
481 //______________________________________________________________________________
482 void AliAODTrack::ConvertAliPIDtoAODPID()
483 {
484   // Converts AliPID array.
485   // The numbering scheme is the same for electrons, muons, pions, kaons, and protons.
486   // Everything else has to be set to zero.
487   if (fPID) {
488     fPID[kDeuteron] = 0.;
489     fPID[kTriton]   = 0.;
490     fPID[kHelium3]  = 0.;
491     fPID[kAlpha]    = 0.;
492     fPID[kUnknown]  = 0.;
493   }
494   return;
495 }
496
497
498 //______________________________________________________________________________
499 template <typename T> void AliAODTrack::SetP(const T *p, const Bool_t cartesian) 
500 {
501   // Set the momentum
502
503   if (p) {
504     if (cartesian) {
505       Double_t pt2 = p[0]*p[0] + p[1]*p[1];
506       Double_t pp  = TMath::Sqrt(pt2 + p[2]*p[2]);
507       
508       fMomentum[0] = TMath::Sqrt(pt2); // pt
509       fMomentum[1] = (pt2 != 0.) ? TMath::Pi()+TMath::ATan2(-p[1], -p[0]) : -999; // phi
510       fMomentum[2] = (pp != 0.) ? TMath::ACos(p[2] / pp) : -999.; // theta
511     } else {
512       fMomentum[0] = p[0];  // pt
513       fMomentum[1] = p[1];  // phi
514       fMomentum[2] = p[2];  // theta
515     }
516   } else {
517     fMomentum[0] = -999.;
518     fMomentum[1] = -999.;
519     fMomentum[2] = -999.;
520   }
521 }
522
523 /*
524 //______________________________________________________________________________
525 template <typename T> void AliAODTrack::SetPosition(const T *x, const Bool_t dca) 
526 {
527   // set the position
528
529   if (x) {
530     if (!dca) {
531       ResetBit(kIsDCA);
532
533       fPosition[0] = x[0];
534       fPosition[1] = x[1];
535       fPosition[2] = x[2];
536     } else {
537       SetBit(kIsDCA);
538       // don't know any better yet
539       fPosition[0] = -999.;
540       fPosition[1] = -999.;
541       fPosition[2] = -999.;
542     }
543   } else {
544     ResetBit(kIsDCA);
545
546     fPosition[0] = -999.;
547     fPosition[1] = -999.;
548     fPosition[2] = -999.;
549   }
550 }
551 */
552 //______________________________________________________________________________
553 void AliAODTrack::SetDCA(Double_t d, Double_t z) 
554 {
555   // set the dca
556   fPosition[0] = d;
557   fPosition[1] = z;
558   fPosition[2] = 0.;
559   SetBit(kIsDCA);
560 }
561
562 //______________________________________________________________________________
563 void AliAODTrack::Print(Option_t* /* option */) const
564 {
565   // prints information about AliAODTrack
566
567   printf("Object name: %s   Track type: %s\n", GetName(), GetTitle()); 
568   printf("        px = %f\n", Px());
569   printf("        py = %f\n", Py());
570   printf("        pz = %f\n", Pz());
571   printf("        pt = %f\n", Pt());
572   printf("      1/pt = %f\n", OneOverPt());
573   printf("     theta = %f\n", Theta());
574   printf("       phi = %f\n", Phi());
575   printf("  chi2/NDF = %f\n", Chi2perNDF());
576   printf("    charge = %d\n", Charge());
577 }
578
579 //______________________________________________________________________________
580 void AliAODTrack::SetMatchTrigger(Int_t matchTrig)
581 {
582   // Set the MUON trigger information
583   switch(matchTrig){
584     case 0: // 0 track does not match trigger
585       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
586       break;
587     case 1: // 1 track match but does not pass pt cut
588       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x40000000;
589       break;
590     case 2: // 2 track match Low pt cut
591       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x80000000;
592       break;
593     case 3: // 3 track match High pt cut
594       fITSMuonClusterMap=fITSMuonClusterMap|0xc0000000;
595       break;
596     default:
597       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
598       AliWarning(Form("unknown case for matchTrig: %d\n",matchTrig));
599   }
600 }
601
602 //______________________________________________________________________________
603 Bool_t AliAODTrack::HitsMuonChamber(Int_t MuonChamber, Int_t cathode) const
604 {
605   // return kTRUE if the track fires the given tracking or trigger chamber.
606   // If the chamber is a trigger one:
607   // - if cathode = 0 or 1, the track matches the corresponding cathode
608   // - if cathode = -1, the track matches both cathodes
609   
610   if (MuonChamber < 0) return kFALSE;
611   
612   if (MuonChamber < 10) return TESTBIT(GetMUONClusterMap(), MuonChamber);
613   
614   if (MuonChamber < 14) {
615     
616     if (cathode < 0) return TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber) &&
617                             TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber+4);
618     
619     if (cathode < 2) return TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber+(1-cathode)*4);
620     
621   }
622   
623   return kFALSE;
624 }
625
626 //______________________________________________________________________________
627 Bool_t AliAODTrack::MatchTriggerDigits() const
628 {
629   // return kTRUE if the track matches a digit on both planes of at least 2 trigger chambers
630   
631   Int_t nMatchedChambers = 0;
632   for (Int_t ich=10; ich<14; ich++) if (HitsMuonChamber(ich)) nMatchedChambers++;
633   
634   return (nMatchedChambers >= 2);
635 }
636
637 //______________________________________________________________________________
638 Int_t AliAODTrack::GetMuonTrigDevSign() const
639 {
640   /// Return the sign of the  MTR deviation
641
642   Int_t signInfo = (Int_t)((fMUONtrigHitsMapTrg>>30)&0x3);
643   // Dummy value for old AODs which do not have the info
644   if ( signInfo == 0 ) return -999;
645   return signInfo - 2;
646 }
647
648 //______________________________________________________________________________
649 Bool_t AliAODTrack::PropagateToDCA(const AliVVertex *vtx, 
650     Double_t b, Double_t maxd, Double_t dz[2], Double_t covar[3])
651 {
652   // compute impact parameters to the vertex vtx and their covariance matrix
653   // b is the Bz, needed to propagate correctly the track to vertex 
654   // only the track parameters are update after the propagation (pos and mom),
655   // not the covariance matrix. This is OK for propagation over short distance
656   // inside the beam pipe.
657   // return kFALSE is something went wrong
658
659   // allowed only for tracks inside the beam pipe
660   Float_t xstart2 = fPosition[0]*fPosition[0]+fPosition[1]*fPosition[1];
661   if(xstart2 > 3.*3.) { // outside beampipe radius
662     AliError("This method can be used only for propagation inside the beam pipe");
663     return kFALSE; 
664   }
665
666   // convert to AliExternalTrackParam
667   AliExternalTrackParam etp; etp.CopyFromVTrack(this);  
668
669   // propagate
670   if(!etp.PropagateToDCA(vtx,b,maxd,dz,covar)) return kFALSE;
671
672   // update track position and momentum
673   Double_t mom[3];
674   etp.GetPxPyPz(mom);
675   SetP(mom,kTRUE);
676   etp.GetXYZ(mom);
677   SetPosition(mom,kFALSE);
678
679
680   return kTRUE;
681 }
682
683 //______________________________________________________________________________
684 Bool_t AliAODTrack::GetPxPyPz(Double_t p[3]) const 
685 {
686     //---------------------------------------------------------------------
687     // This function returns the global track momentum components
688     //---------------------------------------------------------------------
689   p[0]=Px(); p[1]=Py(); p[2]=Pz();
690   return kTRUE;
691 }
692
693
694 //_______________________________________________________________________
695 Float_t AliAODTrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1, Int_t bitType ) const
696 {
697   //
698   // TPC cluster information 
699   // type 0: get fraction of found/findable clusters with neighbourhood definition
700   //      1: findable clusters with neighbourhood definition
701   //      2: found clusters
702   // bitType:
703   //      0 - all cluster used
704   //      1 - clusters  used for the kalman update
705   // definition of findable clusters:
706   //            a cluster is defined as findable if there is another cluster
707   //           within +- nNeighbours pad rows. The idea is to overcome threshold
708   //           effects with a very simple algorithm.
709   //
710
711   
712   Int_t found=0;
713   Int_t findable=0;
714   Int_t last=-nNeighbours;
715   const TBits & clusterMap = (bitType%2==0) ? fTPCClusterMap : fTPCFitMap;
716   
717   Int_t upperBound=clusterMap.GetNbits();
718   if (upperBound>row1) upperBound=row1;
719   for (Int_t i=row0; i<upperBound; ++i){
720     //look to current row
721     if (clusterMap[i]) {
722       last=i;
723       ++found;
724       ++findable;
725       continue;
726     }
727     //look to nNeighbours before
728     if ((i-last)<=nNeighbours) {
729       ++findable;
730       continue;
731     }
732     //look to nNeighbours after
733     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
734       if (clusterMap[j]){
735         ++findable;
736         break;
737       }
738     }
739   }
740   if (type==2) return found;
741   if (type==1) return findable;
742   
743   if (type==0){
744     Float_t fraction=0;
745     if (findable>0) 
746       fraction=(Float_t)found/(Float_t)findable;
747     else 
748       fraction=0;
749     return fraction;
750   }  
751   return 0;  // undefined type - default value
752 }
753
754
755 //______________________________________________________________________________
756 Double_t  AliAODTrack::GetTRDslice(Int_t plane, Int_t slice) const {
757   //
758   // return TRD Pid information
759   //
760   if (!fDetPid) return -1;
761   Double32_t *trdSlices=fDetPid->GetTRDslices();
762   if (!trdSlices) return -1;
763   if ((plane<0) || (plane>=kTRDnPlanes)) {
764     return -1.;
765   }
766
767   Int_t ns=fDetPid->GetTRDnSlices();
768   if ((slice<-1) || (slice>=ns)) {
769     return -1.;
770   }
771
772   if(slice>=0) return trdSlices[plane*ns + slice];
773
774   // return average of the dEdx measurements
775   Double_t q=0.; Double32_t *s = &trdSlices[plane*ns];
776   for (Int_t i=0; i<ns; i++, s++) if((*s)>0.) q+=(*s);
777   return q/ns;
778 }
779
780 //______________________________________________________________________________
781 UChar_t AliAODTrack::GetTRDntrackletsPID() const{
782   //
783   // return number of tracklets calculated from the slices
784   //
785   if(!fDetPid) return -1;
786   return fDetPid->GetTRDntrackletsPID();
787 }
788
789 //______________________________________________________________________________
790 UChar_t AliAODTrack::GetTRDncls(Int_t layer) const {
791   // 
792   // return number of TRD clusters
793   //
794   if(!fDetPid || layer > 5) return -1;
795   if(layer < 0) return fDetPid->GetTRDncls();
796   else return fDetPid->GetTRDncls(layer);
797 }
798
799 //______________________________________________________________________________
800 Double_t AliAODTrack::GetTRDmomentum(Int_t plane, Double_t */*sp*/) const
801 {
802   //Returns momentum estimation
803   // in TRD layer "plane".
804
805   if (!fDetPid) return -1;
806   const Double_t *trdMomentum=fDetPid->GetTRDmomentum();
807
808   if (!trdMomentum) {
809     return -1.;
810   }
811   if ((plane<0) || (plane>=kTRDnPlanes)) {
812     return -1.;
813   }
814
815   return trdMomentum[plane];
816 }
817
818 //_______________________________________________________________________
819 Int_t AliAODTrack::GetTOFBunchCrossing(Double_t b, Bool_t) const 
820 {
821   // Returns the number of bunch crossings after trigger (assuming 25ns spacing)
822   const double kSpacing = 25e3; // min interbanch spacing
823   const double kShift = 0;
824   Int_t bcid = kTOFBCNA; // defualt one
825   if (!IsOn(kTOFout) || !IsOn(kESDpid)) return bcid; // no info
826   //
827   double tdif = GetTOFsignal();
828   if (IsOn(kTIME)) { // integrated time info is there
829     int pid = (int)GetMostProbablePID();
830     double ttimes[10]; 
831     GetIntegratedTimes(ttimes, pid>=AliPID::kSPECIES ? AliPID::kSPECIESC : AliPID::kSPECIES);
832     tdif -= ttimes[pid];
833   }
834   else { // assume integrated time info from TOF radius and momentum
835     const double kRTOF = 385.;
836     const double kCSpeed = 3.e-2; // cm/ps
837     double p = P();
838     if (p<0.001) p = 1.0;
839     double m = M();
840     double path =  kRTOF;     // mean TOF radius
841     if (TMath::Abs(b)>kAlmost0) {  // account for curvature
842       double curv = Pt()/(b*kB2C);
843       if (curv>kAlmost0) {
844         double tgl = Pz()/Pt();
845         path = 2./curv*TMath::ASin(kRTOF*curv/2.)*TMath::Sqrt(1.+tgl*tgl);
846       }
847     }
848     tdif -= path/kCSpeed*TMath::Sqrt(1.+m*m/(p*p));
849   }
850   bcid = TMath::Nint((tdif - kShift)/kSpacing);
851   return bcid;
852 }
853
854 void AliAODTrack::SetDetectorPID(const AliDetectorPID *pid)
855 {
856   //
857   // Set the detector PID
858   //
859   if (fDetectorPID) delete fDetectorPID;
860   fDetectorPID=pid;
861   
862 }
863
864 //_____________________________________________________________________________
865 Double_t AliAODTrack::GetHMPIDsignal() const
866 {
867   if(fAODEvent->GetHMPIDringForTrackID(fID)) return fAODEvent->GetHMPIDringForTrackID(fID)->GetHmpSignal();
868   else return -999.;
869 }
870
871 //_____________________________________________________________________________
872 Double_t AliAODTrack::GetHMPIDoccupancy() const
873 {
874   if(fAODEvent->GetHMPIDringForTrackID(fID)) return fAODEvent->GetHMPIDringForTrackID(fID)->GetHmpOccupancy();
875   else return -999.;
876 }
877
878 //_____________________________________________________________________________
879 Int_t AliAODTrack::GetHMPIDcluIdx() const
880 {
881   if(fAODEvent->GetHMPIDringForTrackID(fID)) return fAODEvent->GetHMPIDringForTrackID(fID)->GetHmpCluIdx();
882   else return -999;
883 }
884
885 //_____________________________________________________________________________
886 void AliAODTrack::GetHMPIDtrk(Float_t &x, Float_t &y, Float_t &th, Float_t &ph) const
887 {
888   x = -999; y = -999.; th = -999.; ph = -999.;
889
890   const AliAODHMPIDrings *ring=fAODEvent->GetHMPIDringForTrackID(fID);
891   if(ring){
892     x  = ring->GetHmpTrackX();
893     y  = ring->GetHmpTrackY();
894     th = ring->GetHmpTrackTheta();
895     ph = ring->GetHmpTrackPhi();
896   }
897 }
898
899 //_____________________________________________________________________________
900 void AliAODTrack::GetHMPIDmip(Float_t &x,Float_t &y,Int_t &q, Int_t &nph) const
901 {
902   x = -999; y = -999.; q = -999; nph = -999;
903   
904   const AliAODHMPIDrings *ring=fAODEvent->GetHMPIDringForTrackID(fID);
905   if(ring){
906     x   = ring->GetHmpMipX();
907     y   = ring->GetHmpMipY();
908     q   = (Int_t)ring->GetHmpMipCharge();
909     nph = (Int_t)ring->GetHmpNumOfPhotonClusters();
910   }
911 }
912
913 //_____________________________________________________________________________
914 Bool_t AliAODTrack::GetOuterHmpPxPyPz(Double_t *p) const 
915
916  if(fAODEvent->GetHMPIDringForTrackID(fID)) {fAODEvent->GetHMPIDringForTrackID(fID)->GetHmpMom(p); return kTRUE;}
917  
918  else return kFALSE;      
919 }
920 //_____________________________________________________________________________
921 Bool_t AliAODTrack::GetXYZAt(Double_t x, Double_t b, Double_t *r) const
922 {
923   //---------------------------------------------------------------------
924   // This function returns the global track position extrapolated to
925   // the radial position "x" (cm) in the magnetic field "b" (kG)
926   //---------------------------------------------------------------------
927
928   //conversion of track parameter representation is
929   //based on the implementation of AliExternalTrackParam::Set(...)
930   //maybe some of this code can be moved to AliVTrack to avoid code duplication
931   Double_t alpha=0.0;
932   Double_t radPos2 = fPosition[0]*fPosition[0]+fPosition[1]*fPosition[1];  
933   Double_t radMax  = 45.; // approximately ITS outer radius
934   if (radPos2 < radMax*radMax) { // inside the ITS     
935     alpha = fMomentum[1]; //TMath::ATan2(fMomentum[1],fMomentum[0]); // fMom is pt,phi,theta!
936   } else { // outside the ITS
937      Float_t phiPos = TMath::Pi()+TMath::ATan2(-fPosition[1], -fPosition[0]);
938      alpha = 
939      TMath::DegToRad()*(20*((((Int_t)(phiPos*TMath::RadToDeg()))/20))+10);
940   }
941   //
942   // Get the vertex of origin and the momentum
943   TVector3 ver(fPosition[0],fPosition[1],fPosition[2]);
944   TVector3 mom(Px(),Py(),Pz());
945   //
946   // Rotate to the local coordinate system
947   ver.RotateZ(-alpha);
948   mom.RotateZ(-alpha);
949
950   Double_t param0 = ver.Y();
951   Double_t param1 = ver.Z();
952   Double_t param2 = TMath::Sin(mom.Phi());
953   Double_t param3 = mom.Pz()/mom.Pt();
954   Double_t param4 = TMath::Sign(1/mom.Pt(),(Double_t)fCharge);
955
956   //calculate the propagated coordinates
957   //this is based on AliExternalTrackParam::GetXYZAt(Double_t x, Double_t b, Double_t *r)
958   Double_t dx=x-ver.X();
959   if(TMath::Abs(dx)<=kAlmost0) return GetXYZ(r);
960
961   Double_t f1=param2;
962   Double_t f2=f1 + dx*param4*b*kB2C;
963
964   if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
965   if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
966   
967   Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
968   r[0] = x;
969   r[1] = param0 + dx*(f1+f2)/(r1+r2);
970   r[2] = param1 + dx*(r2 + f2*(f1+f2)/(r1+r2))*param3;//Thanks to Andrea & Peter
971   return Local2GlobalPosition(r,alpha);
972 }
973
974 //_____________________________________________________________________________
975 Bool_t AliAODTrack::GetXYZatR(Double_t xr,Double_t bz, Double_t *xyz, Double_t* alpSect) const
976 {
977   // This method has 3 modes of behaviour
978   // 1) xyz[3] array is provided but alpSect pointer is 0: calculate the position of track intersection 
979   //    with circle of radius xr and fill it in xyz array
980   // 2) alpSect pointer is provided: find alpha of the sector where the track reaches local coordinate xr
981   //    Note that in this case xr is NOT the radius but the local coordinate.
982   //    If the xyz array is provided, it will be filled by track lab coordinates at local X in this sector
983   // 3) Neither alpSect nor xyz pointers are provided: just check if the track reaches radius xr
984   //
985   //
986   Double_t alpha=0.0;
987   Double_t radPos2 = fPosition[0]*fPosition[0]+fPosition[1]*fPosition[1];  
988   Double_t radMax  = 45.; // approximately ITS outer radius
989   if (radPos2 < radMax*radMax) { // inside the ITS     
990     alpha = fMomentum[1]; //TMath::ATan2(fMomentum[1],fMomentum[0]); // fMom is pt,phi,theta!
991   } else { // outside the ITS
992      Float_t phiPos = TMath::Pi()+TMath::ATan2(-fPosition[1], -fPosition[0]);
993      alpha = 
994      TMath::DegToRad()*(20*((((Int_t)(phiPos*TMath::RadToDeg()))/20))+10);
995   }
996   //  
997   // Get the vertex of origin and the momentum
998   TVector3 ver(fPosition[0],fPosition[1],fPosition[2]);
999   TVector3 mom(Px(),Py(),Pz());
1000   //
1001   // Rotate to the local coordinate system
1002   ver.RotateZ(-alpha);
1003   mom.RotateZ(-alpha);
1004   //
1005   Double_t fx = ver.X();
1006   Double_t fy = ver.Y();
1007   Double_t fz = ver.Z();
1008   Double_t sn = TMath::Sin(mom.Phi());
1009   Double_t tgl = mom.Pz()/mom.Pt();
1010   Double_t crv = TMath::Sign(1/mom.Pt(),(Double_t)fCharge)*bz*kB2C;
1011   //
1012   if ( (TMath::Abs(bz))<kAlmost0Field ) crv=0.;
1013   //
1014   // general circle parameterization:
1015   // x = (r0+tR)cos(phi0) - tR cos(t+phi0)
1016   // y = (r0+tR)sin(phi0) - tR sin(t+phi0)
1017   // where qb is the sign of the curvature, tR is the track's signed radius and r0 
1018   // is the DCA of helix to origin
1019   //
1020   double tR = 1./crv;            // track radius signed
1021   double cs = TMath::Sqrt((1-sn)*(1+sn));
1022   double x0 = fx - sn*tR;        // helix center coordinates
1023   double y0 = fy + cs*tR;
1024   double phi0 = TMath::ATan2(y0,x0);  // angle of PCA wrt to the origin
1025   if (tR<0) phi0 += TMath::Pi();
1026   if      (phi0 > TMath::Pi()) phi0 -= 2.*TMath::Pi();
1027   else if (phi0 <-TMath::Pi()) phi0 += 2.*TMath::Pi();
1028   double cs0 = TMath::Cos(phi0);
1029   double sn0 = TMath::Sin(phi0);
1030   double r0 = x0*cs0 + y0*sn0 - tR; // DCA to origin
1031   double r2R = 1.+r0/tR;
1032   //
1033   //
1034   if (r2R<kAlmost0) return kFALSE;  // helix is centered at the origin, no specific intersection with other concetric circle
1035   if (!xyz && !alpSect) return kTRUE;
1036   double xr2R = xr/tR;
1037   double r2Ri = 1./r2R;
1038   // the intersection cos(t) = [1 + (r0/tR+1)^2 - (r0/tR)^2]/[2(1+r0/tR)]
1039   double cosT = 0.5*(r2R + (1-xr2R*xr2R)*r2Ri);
1040   if ( TMath::Abs(cosT)>kAlmost1 ) {
1041     //    printf("Does not reach : %f %f\n",r0,tR);
1042     return kFALSE; // track does not reach the radius xr
1043   }
1044   //
1045   double t = TMath::ACos(cosT);
1046   if (tR<0) t = -t;
1047   // intersection point
1048   double xyzi[3];
1049   xyzi[0] = x0 - tR*TMath::Cos(t+phi0);
1050   xyzi[1] = y0 - tR*TMath::Sin(t+phi0);
1051   if (xyz) { // if postition is requested, then z is needed:
1052     double t0 = TMath::ATan2(cs,-sn) - phi0;
1053     double z0 = fz - t0*tR*tgl;    
1054     xyzi[2] = z0 + tR*t*tgl;
1055   }
1056   else xyzi[2] = 0;
1057   //
1058   Local2GlobalPosition(xyzi,alpha);
1059   //
1060   if (xyz) {
1061     xyz[0] = xyzi[0];
1062     xyz[1] = xyzi[1];
1063     xyz[2] = xyzi[2];
1064   }
1065   //
1066   if (alpSect) {
1067     double &alp = *alpSect;
1068     // determine the sector of crossing
1069     double phiPos = TMath::Pi()+TMath::ATan2(-xyzi[1],-xyzi[0]);
1070     int sect = ((Int_t)(phiPos*TMath::RadToDeg()))/20;
1071     alp = TMath::DegToRad()*(20*sect+10);
1072     double x2r,f1,f2,r1,r2,dx,dy2dx,yloc=0, ylocMax = xr*TMath::Tan(TMath::Pi()/18); // min max Y within sector at given X
1073     //
1074     while(1) {
1075       Double_t ca=TMath::Cos(alp-alpha), sa=TMath::Sin(alp-alpha);
1076       if ((cs*ca+sn*sa)<0) {
1077         AliDebug(1,Form("Rotation to target sector impossible: local cos(phi) would become %.2f",cs*ca+sn*sa));
1078         return kFALSE;
1079       }
1080       //
1081       f1 = sn*ca - cs*sa;
1082       if (TMath::Abs(f1) >= kAlmost1) {
1083         AliDebug(1,Form("Rotation to target sector impossible: local sin(phi) would become %.2f",f1));
1084         return kFALSE;
1085       }
1086       //
1087       double tmpX =  fx*ca + fy*sa;
1088       double tmpY = -fx*sa + fy*ca;
1089       //
1090       // estimate Y at X=xr
1091       dx=xr-tmpX;
1092       x2r = crv*dx;
1093       f2=f1 + x2r;
1094       if (TMath::Abs(f2) >= kAlmost1) {
1095         AliDebug(1,Form("Propagation in target sector failed ! %.10e",f2));
1096         return kFALSE;
1097       }
1098       r1 = TMath::Sqrt((1.-f1)*(1.+f1));
1099       r2 = TMath::Sqrt((1.-f2)*(1.+f2));
1100       dy2dx = (f1+f2)/(r1+r2);
1101       yloc = tmpY + dx*dy2dx;
1102       if      (yloc>ylocMax)  {alp += 2*TMath::Pi()/18; sect++;}
1103       else if (yloc<-ylocMax) {alp -= 2*TMath::Pi()/18; sect--;}
1104       else break;
1105       if      (alp >= TMath::Pi()) alp -= 2*TMath::Pi();
1106       else if (alp < -TMath::Pi()) alp += 2*TMath::Pi();
1107       //      if (sect>=18) sect = 0;
1108       //      if (sect<=0) sect = 17;
1109     }
1110     //
1111     // if alpha was requested, then recalculate the position at intersection in sector
1112     if (xyz) {
1113       xyz[0] = xr;
1114       xyz[1] = yloc;
1115       if (TMath::Abs(x2r)<0.05) xyz[2] = fz + dx*(r2 + f2*dy2dx)*tgl;
1116       else {
1117         // for small dx/R the linear apporximation of the arc by the segment is OK,
1118         // but at large dx/R the error is very large and leads to incorrect Z propagation
1119         // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
1120         // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
1121         // Similarly, the rotation angle in linear in dx only for dx<<R
1122         double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx);   // distance from old position to new one
1123         double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
1124         xyz[2] = fz + rot/crv*tgl;
1125       }
1126       Local2GlobalPosition(xyz,alp);
1127     }
1128   }
1129   return kTRUE;    
1130   //
1131 }
1132
1133 //_______________________________________________________
1134 void  AliAODTrack::GetITSdEdxSamples(Double_t s[4]) const
1135 {
1136   // get ITS dedx samples
1137   if (!fDetPid) for (int i=4;i--;) s[i]=0;
1138   else          for (int i=4;i--;) s[i] = fDetPid->GetITSdEdxSample(i);
1139 }
1140
1141 //_____________________________________________
1142 Double_t AliAODTrack::GetMassForTracking() const
1143 {
1144   int pid = fPIDForTracking;
1145   if (pid<AliPID::kPion) pid = AliPID::kPion;
1146   double m = AliPID::ParticleMass(fPIDForTracking);
1147   return (fPIDForTracking==AliPID::kHe3 || fPIDForTracking==AliPID::kAlpha) ? -m : m;
1148 }
1149 //_______________________________________________________
1150 const AliTOFHeader* AliAODTrack::GetTOFHeader() const {
1151   return fAODEvent->GetTOFHeader();
1152 }
1153 //_______________________________________________________
1154 Int_t AliAODTrack::GetNcls(Int_t idet) const
1155 {
1156   // Get number of clusters by subdetector index
1157   //
1158   Int_t ncls = 0;
1159   switch(idet){
1160   case 0:
1161     ncls = GetITSNcls();
1162     break;
1163   case 1:
1164     ncls = (Int_t)GetTPCNcls();
1165     break;
1166   case 2:
1167     ncls = (Int_t)GetTRDncls();
1168     break;
1169   case 3:
1170     break;
1171     /*if (fTOFindex != -1)
1172       ncls = 1;*/
1173     break;
1174   case 4: //PHOS
1175     break;
1176   case 5: //HMPID
1177     break;
1178     if ((GetHMPIDcluIdx() >= 0) && (GetHMPIDcluIdx() < 7000000)) {
1179       if ((GetHMPIDcluIdx()%1000000 != 9999) && (GetHMPIDcluIdx()%1000000 != 99999)) {
1180         ncls = 1;
1181         }
1182     }    
1183     break;
1184   default:
1185     break;
1186   }
1187   return ncls;
1188 }