Add a new variable, that is the radial position of the track at the end
[u/mrichter/AliRoot.git] / STEER / AliAODTrack.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //-------------------------------------------------------------------------
19 //     AOD track implementation of AliVTrack
20 //     Author: Markus Oldenburg, CERN
21 //     Markus.Oldenburg@cern.ch
22 //-------------------------------------------------------------------------
23
24 #include "AliLog.h"
25 #include "AliAODTrack.h"
26
27 ClassImp(AliAODTrack)
28
29 //______________________________________________________________________________
30 AliAODTrack::AliAODTrack() : 
31   AliVTrack(),
32   fRAtAbsorberEnd(0.),
33   fChi2perNDF(-999.),
34   fChi2MatchTrigger(0.),
35   fFlags(0),
36   fLabel(-999),
37   fITSMuonClusterMap(0),
38   fFilterMap(0),
39   fTPCClusterMap(),
40   fTPCSharedMap(),
41   fID(-999),
42   fCharge(-99),
43   fType(kUndef),
44   fCovMatrix(NULL),
45   fDetPid(NULL),
46   fProdVertex(NULL)
47 {
48   // default constructor
49
50   SetP();
51   SetPosition((Float_t*)NULL);
52   SetXYAtDCA(-999., -999.);
53   SetPxPyPzAtDCA(-999., -999., -999.);
54   SetPID((Float_t*)NULL);
55 }
56
57 //______________________________________________________________________________
58 AliAODTrack::AliAODTrack(Short_t id,
59                          Int_t label, 
60                          Double_t p[3],
61                          Bool_t cartesian,
62                          Double_t x[3],
63                          Bool_t isDCA,
64                          Double_t covMatrix[21],
65                          Short_t charge,
66                          UChar_t itsClusMap,
67                          Double_t pid[10],
68                          AliAODVertex *prodVertex,
69                          Bool_t usedForVtxFit,
70                          Bool_t usedForPrimVtxFit,
71                          AODTrk_t ttype,
72                          UInt_t selectInfo,
73                          Float_t chi2perNDF) :
74   AliVTrack(),
75   fRAtAbsorberEnd(0.),
76   fChi2perNDF(chi2perNDF),
77   fChi2MatchTrigger(0.),
78   fFlags(0),
79   fLabel(label),
80   fITSMuonClusterMap(0),
81   fFilterMap(selectInfo),
82   fTPCClusterMap(),
83   fTPCSharedMap(),
84   fID(id),
85   fCharge(charge),
86   fType(ttype),
87   fCovMatrix(NULL),
88   fDetPid(NULL),
89   fProdVertex(prodVertex)
90 {
91   // constructor
92  
93   SetP(p, cartesian);
94   SetPosition(x, isDCA);
95   SetXYAtDCA(-999., -999.);
96   SetPxPyPzAtDCA(-999., -999., -999.);
97   SetUsedForVtxFit(usedForVtxFit);
98   SetUsedForPrimVtxFit(usedForPrimVtxFit);
99   if(covMatrix) SetCovMatrix(covMatrix);
100   SetPID(pid);
101   SetITSClusterMap(itsClusMap);
102 }
103
104 //______________________________________________________________________________
105 AliAODTrack::AliAODTrack(Short_t id,
106                          Int_t label, 
107                          Float_t p[3],
108                          Bool_t cartesian,
109                          Float_t x[3],
110                          Bool_t isDCA,
111                          Float_t covMatrix[21],
112                          Short_t charge,
113                          UChar_t itsClusMap,
114                          Float_t pid[10],
115                          AliAODVertex *prodVertex,
116                          Bool_t usedForVtxFit,
117                          Bool_t usedForPrimVtxFit,
118                          AODTrk_t ttype,
119                          UInt_t selectInfo,
120                          Float_t chi2perNDF) :
121   AliVTrack(),
122   fRAtAbsorberEnd(0.),
123   fChi2perNDF(chi2perNDF),
124   fChi2MatchTrigger(0.),
125   fFlags(0),
126   fLabel(label),
127   fITSMuonClusterMap(0),
128   fFilterMap(selectInfo),
129   fTPCClusterMap(),
130   fTPCSharedMap(),
131   fID(id),
132   fCharge(charge),
133   fType(ttype),
134   fCovMatrix(NULL),
135   fDetPid(NULL),
136   fProdVertex(prodVertex)
137 {
138   // constructor
139  
140   SetP(p, cartesian);
141   SetPosition(x, isDCA);
142   SetXYAtDCA(-999., -999.);
143   SetPxPyPzAtDCA(-999., -999., -999.);
144   SetUsedForVtxFit(usedForVtxFit);
145   SetUsedForPrimVtxFit(usedForPrimVtxFit);
146   if(covMatrix) SetCovMatrix(covMatrix);
147   SetPID(pid);
148   SetITSClusterMap(itsClusMap);
149 }
150
151 //______________________________________________________________________________
152 AliAODTrack::~AliAODTrack() 
153 {
154   // destructor
155   delete fCovMatrix;
156   delete fDetPid;
157 }
158
159
160 //______________________________________________________________________________
161 AliAODTrack::AliAODTrack(const AliAODTrack& trk) :
162   AliVTrack(trk),
163   fRAtAbsorberEnd(trk.fRAtAbsorberEnd),
164   fChi2perNDF(trk.fChi2perNDF),
165   fChi2MatchTrigger(trk.fChi2MatchTrigger),
166   fFlags(trk.fFlags),
167   fLabel(trk.fLabel),
168   fITSMuonClusterMap(trk.fITSMuonClusterMap),
169   fFilterMap(trk.fFilterMap),
170   fTPCClusterMap(trk.fTPCClusterMap),
171   fTPCSharedMap(trk.fTPCSharedMap),
172   fID(trk.fID),
173   fCharge(trk.fCharge),
174   fType(trk.fType),
175   fCovMatrix(NULL),
176   fDetPid(NULL),
177   fProdVertex(trk.fProdVertex)
178 {
179   // Copy constructor
180
181   trk.GetP(fMomentum);
182   trk.GetPosition(fPosition);
183   SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
184   SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
185   SetUsedForVtxFit(trk.GetUsedForVtxFit());
186   SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
187   if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
188   if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
189   SetPID(trk.fPID);
190 }
191
192 //______________________________________________________________________________
193 AliAODTrack& AliAODTrack::operator=(const AliAODTrack& trk)
194 {
195   // Assignment operator
196   if(this!=&trk) {
197
198     AliVTrack::operator=(trk);
199
200     trk.GetP(fMomentum);
201     trk.GetPosition(fPosition);
202     trk.GetPID(fPID);
203
204     SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
205     SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
206     
207     fRAtAbsorberEnd = trk.fRAtAbsorberEnd;
208     
209     fChi2perNDF = trk.fChi2perNDF;
210     fChi2MatchTrigger = trk.fChi2MatchTrigger;
211
212     fFlags = trk.fFlags;
213     fLabel = trk.fLabel;    
214     
215     fITSMuonClusterMap = trk.fITSMuonClusterMap;
216     fFilterMap = trk.fFilterMap;
217
218     fID = trk.fID;
219
220     fCharge = trk.fCharge;
221     fType = trk.fType;
222
223     delete fCovMatrix;
224     if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
225     else fCovMatrix=NULL;
226     fProdVertex = trk.fProdVertex;
227
228     SetUsedForVtxFit(trk.GetUsedForVtxFit());
229     SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
230
231     delete fDetPid;
232     if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
233     else fDetPid=NULL;
234   }
235
236   return *this;
237 }
238
239 //______________________________________________________________________________
240 Double_t AliAODTrack::M(AODTrkPID_t pid) const
241 {
242   // Returns the mass.
243   // Masses for nuclei don't exist in the PDG tables, therefore they were put by hand.
244
245   switch (pid) {
246
247   case kElectron :
248     return 0.000510999; //TDatabasePDG::Instance()->GetParticle(11/*::kElectron*/)->Mass();
249     break;
250
251   case kMuon :
252     return 0.1056584; //TDatabasePDG::Instance()->GetParticle(13/*::kMuonMinus*/)->Mass();
253     break;
254
255   case kPion :
256     return 0.13957; //TDatabasePDG::Instance()->GetParticle(211/*::kPiPlus*/)->Mass();
257     break;
258
259   case kKaon :
260     return 0.4937; //TDatabasePDG::Instance()->GetParticle(321/*::kKPlus*/)->Mass();
261     break;
262
263   case kProton :
264     return 0.9382720; //TDatabasePDG::Instance()->GetParticle(2212/*::kProton*/)->Mass();
265     break;
266
267   case kDeuteron :
268     return 1.8756; //TDatabasePDG::Instance()->GetParticle(1000010020)->Mass();
269     break;
270
271   case kTriton :
272     return 2.8089; //TDatabasePDG::Instance()->GetParticle(1000010030)->Mass();
273     break;
274
275   case kHelium3 :
276     return 2.8084; //TDatabasePDG::Instance()->GetParticle(1000020030)->Mass();
277     break;
278
279   case kAlpha :
280     return 3.7274; //TDatabasePDG::Instance()->GetParticle(1000020040)->Mass();
281     break;
282
283   case kUnknown :
284     return -999.;
285     break;
286
287   default :
288     return -999.;
289   }
290 }
291
292 //______________________________________________________________________________
293 Double_t AliAODTrack::E(AODTrkPID_t pid) const
294 {
295   // Returns the energy of the particle of a given pid.
296   
297   if (pid != kUnknown) { // particle was identified
298     Double_t m = M(pid);
299     return TMath::Sqrt(P()*P() + m*m);
300   } else { // pid unknown
301     return -999.;
302   }
303 }
304
305 //______________________________________________________________________________
306 Double_t AliAODTrack::Y(AODTrkPID_t pid) const
307 {
308   // Returns the rapidity of a particle of a given pid.
309   
310   if (pid != kUnknown) { // particle was identified
311     Double_t e = E(pid);
312     Double_t pz = Pz();
313     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
314       return 0.5*TMath::Log((e+pz)/(e-pz));
315     } else { // energy not known or equal to pz
316       return -999.;
317     }
318   } else { // pid unknown
319     return -999.;
320   }
321 }
322
323 //______________________________________________________________________________
324 Double_t AliAODTrack::Y(Double_t m) const
325 {
326   // Returns the rapidity of a particle of a given mass.
327   
328   if (m >= 0.) { // mass makes sense
329     Double_t e = E(m);
330     Double_t pz = Pz();
331     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
332       return 0.5*TMath::Log((e+pz)/(e-pz));
333     } else { // energy not known or equal to pz
334       return -999.;
335     }
336   } else { // pid unknown
337     return -999.;
338   }
339 }
340
341 //______________________________________________________________________________
342 AliAODTrack::AODTrkPID_t AliAODTrack::GetMostProbablePID() const 
343 {
344   // Returns the most probable PID array element.
345   
346   Int_t nPID = 10;
347   if (fPID) {
348     AODTrkPID_t loc = kUnknown;
349     Double_t max = 0.;
350     Bool_t allTheSame = kTRUE;
351     
352     for (Int_t iPID = 0; iPID < nPID; iPID++) {
353       if (fPID[iPID] >= max) {
354         if (fPID[iPID] > max) {
355           allTheSame = kFALSE;
356           max = fPID[iPID];
357           loc = (AODTrkPID_t)iPID;
358         } else {
359           allTheSame = kTRUE;
360         }
361       }
362     }
363     
364     return allTheSame ? kUnknown : loc;
365   } else {
366     return kUnknown;
367   }
368 }
369
370 //______________________________________________________________________________
371 void AliAODTrack::ConvertAliPIDtoAODPID()
372 {
373   // Converts AliPID array.
374   // The numbering scheme is the same for electrons, muons, pions, kaons, and protons.
375   // Everything else has to be set to zero.
376
377   fPID[kDeuteron] = 0.;
378   fPID[kTriton]   = 0.;
379   fPID[kHelium3]  = 0.;
380   fPID[kAlpha]    = 0.;
381   fPID[kUnknown]  = 0.;
382   
383   return;
384 }
385
386
387 //______________________________________________________________________________
388 template <class T> void AliAODTrack::SetP(const T *p, const Bool_t cartesian) 
389 {
390   // Set the momentum
391
392   if (p) {
393     if (cartesian) {
394       Double_t pt2 = p[0]*p[0] + p[1]*p[1];
395       Double_t pp  = TMath::Sqrt(pt2 + p[2]*p[2]);
396       
397       fMomentum[0] = TMath::Sqrt(pt2); // pt
398       fMomentum[1] = (pt2 != 0.) ? TMath::Pi()+TMath::ATan2(-p[1], -p[0]) : -999; // phi
399       fMomentum[2] = (pp != 0.) ? TMath::ACos(p[2] / pp) : -999.; // theta
400     } else {
401       fMomentum[0] = p[0];  // pt
402       fMomentum[1] = p[1];  // phi
403       fMomentum[2] = p[2];  // theta
404     }
405   } else {
406     fMomentum[0] = -999.;
407     fMomentum[1] = -999.;
408     fMomentum[2] = -999.;
409   }
410 }
411
412 //______________________________________________________________________________
413 template <class T> void AliAODTrack::SetPosition(const T *x, const Bool_t dca) 
414 {
415   // set the position
416
417   if (x) {
418     if (!dca) {
419       ResetBit(kIsDCA);
420
421       fPosition[0] = x[0];
422       fPosition[1] = x[1];
423       fPosition[2] = x[2];
424     } else {
425       SetBit(kIsDCA);
426       // don't know any better yet
427       fPosition[0] = -999.;
428       fPosition[1] = -999.;
429       fPosition[2] = -999.;
430     }
431   } else {
432     ResetBit(kIsDCA);
433
434     fPosition[0] = -999.;
435     fPosition[1] = -999.;
436     fPosition[2] = -999.;
437   }
438 }
439
440 //______________________________________________________________________________
441 void AliAODTrack::SetDCA(Double_t d, Double_t z) 
442 {
443   // set the dca
444   fPosition[0] = d;
445   fPosition[1] = z;
446   fPosition[2] = 0.;
447   SetBit(kIsDCA);
448 }
449
450 //______________________________________________________________________________
451 void AliAODTrack::Print(Option_t* /* option */) const
452 {
453   // prints information about AliAODTrack
454
455   printf("Object name: %s   Track type: %s\n", GetName(), GetTitle()); 
456   printf("        px = %f\n", Px());
457   printf("        py = %f\n", Py());
458   printf("        pz = %f\n", Pz());
459   printf("        pt = %f\n", Pt());
460   printf("      1/pt = %f\n", OneOverPt());
461   printf("     theta = %f\n", Theta());
462   printf("       phi = %f\n", Phi());
463   printf("  chi2/NDF = %f\n", Chi2perNDF());
464   printf("    charge = %d\n", Charge());
465 }
466
467 //______________________________________________________________________________
468 void AliAODTrack::SetMatchTrigger(Int_t matchTrig)
469 {
470   // Set the MUON trigger information
471   switch(matchTrig){
472     case 0: // 0 track does not match trigger
473       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
474       break;
475     case 1: // 1 track match but does not pass pt cut
476       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x40000000;
477       break;
478     case 2: // 2 track match Low pt cut
479       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x80000000;
480       break;
481     case 3: // 3 track match High pt cut
482       fITSMuonClusterMap=fITSMuonClusterMap|0xc0000000;
483       break;
484     default:
485       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
486       AliWarning(Form("unknown case for matchTrig: %d\n",matchTrig));
487   }
488 }
489
490 //______________________________________________________________________________
491 Bool_t AliAODTrack::HitsMuonChamber(Int_t MuonChamber, Int_t cathode) const
492 {
493   // return kTRUE if the track fires the given tracking or trigger chamber.
494   // If the chamber is a trigger one:
495   // - if cathode = 0 or 1, the track matches the corresponding cathode
496   // - if cathode = -1, the track matches both cathodes
497   
498   if (MuonChamber < 0) return kFALSE;
499   
500   if (MuonChamber < 10) return TESTBIT(GetMUONClusterMap(), MuonChamber);
501   
502   if (MuonChamber < 14) {
503     
504     if (cathode < 0) return TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber) &&
505                             TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber+4);
506     
507     if (cathode < 2) return TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber+(1-cathode)*4);
508     
509   }
510   
511   return kFALSE;
512 }
513
514 //______________________________________________________________________________
515 Bool_t AliAODTrack::MatchTriggerDigits() const
516 {
517   // return kTRUE if the track matches a digit on both planes of at least 2 trigger chambers
518   
519   Int_t nMatchedChambers = 0;
520   for (Int_t ich=10; ich<14; ich++) if (HitsMuonChamber(ich)) nMatchedChambers++;
521   
522   return (nMatchedChambers >= 2);
523 }
524