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