Adding RICH and updating PHOS part (Yu.Belikov)
[u/mrichter/AliRoot.git] / STEER / AliESDtrack.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 //-----------------------------------------------------------------
17 //           Implementation of the ESD track class
18 //   ESD = Event Summary Data
19 //   This is the class to deal with during the phisical analysis of data
20 //      Origin: Iouri Belikov, CERN
21 //      e-mail: Jouri.Belikov@cern.ch
22 //-----------------------------------------------------------------
23
24 #include "TMath.h"
25
26 #include "AliESDtrack.h"
27 #include "AliKalmanTrack.h"
28
29 ClassImp(AliESDtrack)
30
31 //_______________________________________________________________________
32 AliESDtrack::AliESDtrack() : 
33 fFlags(0),
34 fLabel(0),
35 fTrackLength(0),
36 fStopVertex(0),
37 fRalpha(0),
38 fRx(0),
39 fCalpha(0),
40 fCx(0),
41 fCchi2(1e10),
42 fIalpha(0),
43 fIx(0),
44 fOalpha(0),
45 fOx(0),
46 fITSchi2(0),
47 fITSncls(0),
48 fITSsignal(0),
49 fTPCchi2(0),
50 fTPCncls(0),
51 fTPCClusterMap(159),//number of padrows
52 fTPCsignal(0),
53 fTRDchi2(0),
54 fTRDncls(0),
55 fTRDsignal(0),
56 fTOFchi2(0),
57 fTOFindex(0),
58 fTOFsignal(-1),
59 fPHOSsignal(-1),
60 fRICHsignal(-1)
61 {
62   //
63   // The default ESD constructor 
64   //
65   for (Int_t i=0; i<kSPECIES; i++) {
66     fTrackTime[i]=0.;
67     fR[i]=1.;
68     fITSr[i]=1.;
69     fTPCr[i]=1.;
70     fTRDr[i]=1.;
71     fTOFr[i]=1.;
72     fPHOSr[i]=1.;
73     fRICHr[i]=1.;
74   }
75   fPHOSpos[0]=fPHOSpos[1]=fPHOSpos[2]=0.;
76   Int_t i;
77   for (i=0; i<5; i++)  { fRp[i]=0.; fCp[i]=0.; fIp[i]=0.; fOp[i]=0.;}
78   for (i=0; i<15; i++) { fRc[i]=0.; fCc[i]=0.; fIc[i]=0.; fOc[i]=0.;   }
79   for (i=0; i<6; i++)  { fITSindex[i]=0; }
80   for (i=0; i<180; i++){ fTPCindex[i]=0; }
81   for (i=0; i<90; i++) { fTRDindex[i]=0; }
82   fTPCLabel = 0;
83   fTRDLabel = 0;
84   fITSLabel = 0;
85
86 }
87
88 //_______________________________________________________________________
89 Double_t AliESDtrack::GetMass() const {
90   // Returns the mass of the most probable particle type
91   Float_t max=0.;
92   Int_t k=-1;
93   for (Int_t i=0; i<kSPECIES; i++) {
94     if (fR[i]>max) {k=i; max=fR[i];}
95   }
96   if (k==0) return 0.00051;
97   if (k==1) return 0.10566;
98   if (k==2||k==-1) return 0.13957;
99   if (k==3) return 0.49368;
100   if (k==4) return 0.93827;
101   Warning("GetMass()","Undefined mass !");
102   return 0.13957;
103 }
104
105 //_______________________________________________________________________
106 Bool_t AliESDtrack::UpdateTrackParams(AliKalmanTrack *t, ULong_t flags) {
107   //
108   // This function updates track's running parameters 
109   //
110   SetStatus(flags);
111   fLabel=t->GetLabel();
112
113   if (t->IsStartedTimeIntegral()) {
114     SetStatus(kTIME);
115     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
116     SetIntegratedLength(t->GetIntegratedLength());
117   }
118
119   fRalpha=t->GetAlpha();
120   t->GetExternalParameters(fRx,fRp);
121   t->GetExternalCovariance(fRc);
122
123   switch (flags) {
124     
125   case kITSin: case kITSout: case kITSrefit:
126     fITSncls=t->GetNumberOfClusters();
127     fITSchi2=t->GetChi2();
128     for (Int_t i=0;i<fITSncls;i++) fITSindex[i]=t->GetClusterIndex(i);
129     fITSsignal=t->GetPIDsignal();
130     fITSLabel = t->GetLabel();
131     break;
132     
133   case kTPCin: case kTPCrefit:
134     fTPCLabel = t->GetLabel();
135     fIalpha=fRalpha;
136     fIx=fRx;
137     {
138       Int_t i;
139       for (i=0; i<5; i++) fIp[i]=fRp[i];
140       for (i=0; i<15;i++) fIc[i]=fRc[i];
141     }
142   case kTPCout:
143   
144     fTPCncls=t->GetNumberOfClusters();
145     fTPCchi2=t->GetChi2();
146     
147      {//prevrow must be declared in separate namespace, otherwise compiler cries:
148       //"jump to case label crosses initialization of `Int_t prevrow'"
149        Int_t prevrow = -1;
150        //       for (Int_t i=0;i<fTPCncls;i++) 
151        for (Int_t i=0;i<160;i++) 
152         {
153           fTPCindex[i]=t->GetClusterIndex(i);
154
155           // Piotr's Cluster Map for HBT  
156           // ### please change accordingly if cluster array is changing 
157           // to "New TPC Tracking" style (with gaps in array) 
158           Int_t idx = fTPCindex[i];
159           Int_t sect = (idx&0xff000000)>>24;
160           Int_t row = (idx&0x00ff0000)>>16;
161           if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
162
163           fTPCClusterMap.SetBitNumber(row,kTRUE);
164
165           //Fill the gap between previous row and this row with 0 bits
166           //In case  ###  pleas change it as well - just set bit 0 in case there 
167           //is no associated clusters for current "i"
168           if (prevrow < 0) 
169            {
170              prevrow = row;//if previous bit was not assigned yet == this is the first one
171            }
172           else
173            { //we don't know the order (inner to outer or reverse)
174              //just to be save in case it is going to change
175              Int_t n = 0, m = 0;
176              if (prevrow < row)
177               {
178                 n = prevrow;
179                 m = row;
180               }
181              else
182               {
183                 n = row;
184                 m = prevrow;
185               }
186
187              for (Int_t j = n+1; j < m; j++)
188               {
189                 fTPCClusterMap.SetBitNumber(j,kFALSE);
190               }
191              prevrow = row; 
192            }
193           // End Of Piotr's Cluster Map for HBT
194         }
195      }
196     fTPCsignal=t->GetPIDsignal();
197     {Double_t mass=t->GetMass();    // preliminary mass setting 
198     if (mass>0.5) fR[4]=1.;         //        used by
199     else if (mass<0.4) fR[2]=1.;    // the ITS reconstruction
200     else fR[3]=1.;}
201                      //
202     break;
203
204   case kTRDout:
205     { //requested by the PHOS  ("temporary solution")
206       Double_t r=460.;
207       if (t->PropagateTo(r,30.,0.)) {  
208          fOalpha=t->GetAlpha();
209          t->GetExternalParameters(fOx,fOp);
210          t->GetExternalCovariance(fOc);
211       }
212     }
213   case kTRDin: case kTRDrefit:
214     fTRDLabel = t->GetLabel();
215
216     fTRDncls=t->GetNumberOfClusters();
217     fTRDchi2=t->GetChi2();
218     for (Int_t i=0;i<fTRDncls;i++) fTRDindex[i]=t->GetClusterIndex(i);
219     fTRDsignal=t->GetPIDsignal();
220     break;
221   case kTRDStop:
222     break;
223   default: 
224     Error("UpdateTrackParams()","Wrong flag !\n");
225     return kFALSE;
226   }
227
228   return kTRUE;
229 }
230
231 //_______________________________________________________________________
232 void 
233 AliESDtrack::SetConstrainedTrackParams(AliKalmanTrack *t, Double_t chi2) {
234   //
235   // This function sets the constrained track parameters 
236   //
237   fCalpha=t->GetAlpha();
238   t->GetExternalParameters(fCx,fCp);
239   t->GetExternalCovariance(fCc);
240   fCchi2=chi2;
241 }
242
243
244 //_______________________________________________________________________
245 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
246   //---------------------------------------------------------------------
247   // This function returns external representation of the track parameters
248   //---------------------------------------------------------------------
249   x=fRx;
250   for (Int_t i=0; i<5; i++) p[i]=fRp[i];
251 }
252 //_______________________________________________________________________
253 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
254   //---------------------------------------------------------------------
255   // This function returns external representation of the cov. matrix
256   //---------------------------------------------------------------------
257   for (Int_t i=0; i<15; i++) cov[i]=fRc[i];
258 }
259
260
261 //_______________________________________________________________________
262 void 
263 AliESDtrack::GetConstrainedExternalParameters(Double_t &x, Double_t p[5])const{
264   //---------------------------------------------------------------------
265   // This function returns the constrained external track parameters
266   //---------------------------------------------------------------------
267   x=fCx;
268   for (Int_t i=0; i<5; i++) p[i]=fCp[i];
269 }
270 //_______________________________________________________________________
271 void 
272 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
273   //---------------------------------------------------------------------
274   // This function returns the constrained external cov. matrix
275   //---------------------------------------------------------------------
276   for (Int_t i=0; i<15; i++) c[i]=fCc[i];
277 }
278
279
280 Double_t AliESDtrack::GetP() const {
281   //---------------------------------------------------------------------
282   // This function returns the track momentum
283   //---------------------------------------------------------------------
284   Double_t lam=TMath::ATan(fRp[3]);
285   Double_t pt=1./TMath::Abs(fRp[4]);
286   return pt/TMath::Cos(lam);
287 }
288
289 void AliESDtrack::GetConstrainedPxPyPz(Double_t *p) const {
290   //---------------------------------------------------------------------
291   // This function returns the constrained global track momentum components
292   //---------------------------------------------------------------------
293   Double_t phi=TMath::ASin(fCp[2]) + fCalpha;
294   Double_t pt=1./TMath::Abs(fCp[4]);
295   p[0]=pt*TMath::Cos(phi); p[1]=pt*TMath::Sin(phi); p[2]=pt*fCp[3]; 
296 }
297 void AliESDtrack::GetConstrainedXYZ(Double_t *xyz) const {
298   //---------------------------------------------------------------------
299   // This function returns the global track position
300   //---------------------------------------------------------------------
301   Double_t phi=TMath::ATan2(fCp[0],fCx) + fCalpha;
302   Double_t r=TMath::Sqrt(fCx*fCx + fCp[0]*fCp[0]);
303   xyz[0]=r*TMath::Cos(phi); xyz[1]=r*TMath::Sin(phi); xyz[2]=fCp[1]; 
304 }
305
306 void AliESDtrack::GetPxPyPz(Double_t *p) const {
307   //---------------------------------------------------------------------
308   // This function returns the global track momentum components
309   //---------------------------------------------------------------------
310   Double_t phi=TMath::ASin(fRp[2]) + fRalpha;
311   Double_t pt=1./TMath::Abs(fRp[4]);
312   p[0]=pt*TMath::Cos(phi); p[1]=pt*TMath::Sin(phi); p[2]=pt*fRp[3]; 
313 }
314 void AliESDtrack::GetXYZ(Double_t *xyz) const {
315   //---------------------------------------------------------------------
316   // This function returns the global track position
317   //---------------------------------------------------------------------
318   Double_t phi=TMath::ATan2(fRp[0],fRx) + fRalpha;
319   Double_t r=TMath::Sqrt(fRx*fRx + fRp[0]*fRp[0]);
320   xyz[0]=r*TMath::Cos(phi); xyz[1]=r*TMath::Sin(phi); xyz[2]=fRp[1]; 
321 }
322
323
324 void AliESDtrack::GetInnerPxPyPz(Double_t *p) const {
325   //---------------------------------------------------------------------
326   // This function returns the global track momentum components
327   // af the entrance of the TPC
328   //---------------------------------------------------------------------
329   if (fIx==0) {p[0]=p[1]=p[2]=0.; return;}
330   Double_t phi=TMath::ASin(fIp[2]) + fIalpha;
331   Double_t pt=1./TMath::Abs(fIp[4]);
332   p[0]=pt*TMath::Cos(phi); p[1]=pt*TMath::Sin(phi); p[2]=pt*fIp[3]; 
333 }
334
335 void AliESDtrack::GetInnerXYZ(Double_t *xyz) const {
336   //---------------------------------------------------------------------
337   // This function returns the global track position
338   // af the entrance of the TPC
339   //---------------------------------------------------------------------
340   if (fIx==0) {xyz[0]=xyz[1]=xyz[2]=0.; return;}
341   Double_t phi=TMath::ATan2(fIp[0],fIx) + fIalpha;
342   Double_t r=TMath::Sqrt(fIx*fIx + fIp[0]*fIp[0]);
343   xyz[0]=r*TMath::Cos(phi); xyz[1]=r*TMath::Sin(phi); xyz[2]=fIp[1]; 
344 }
345
346 void AliESDtrack::GetInnerExternalParameters(Double_t &x, Double_t p[5]) const 
347 {
348   //skowron
349  //---------------------------------------------------------------------
350   // This function returns external representation of the track parameters at Inner Layer of TPC
351   //---------------------------------------------------------------------
352   x=fIx;
353   for (Int_t i=0; i<5; i++) p[i]=fIp[i];
354 }
355 void AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const
356 {
357  //skowron
358  //---------------------------------------------------------------------
359  // This function returns external representation of the cov. matrix at Inner Layer of TPC
360  //---------------------------------------------------------------------
361  for (Int_t i=0; i<15; i++) cov[i]=fIc[i];
362  
363 }
364
365 void AliESDtrack::GetOuterPxPyPz(Double_t *p) const {
366   //---------------------------------------------------------------------
367   // This function returns the global track momentum components
368   // af the radius of the PHOS
369   //---------------------------------------------------------------------
370   if (fOx==0) {p[0]=p[1]=p[2]=0.; return;}
371   Double_t phi=TMath::ASin(fOp[2]) + fOalpha;
372   Double_t pt=1./TMath::Abs(fOp[4]);
373   p[0]=pt*TMath::Cos(phi); p[1]=pt*TMath::Sin(phi); p[2]=pt*fOp[3]; 
374 }
375
376 void AliESDtrack::GetOuterXYZ(Double_t *xyz) const {
377   //---------------------------------------------------------------------
378   // This function returns the global track position
379   // af the radius of the PHOS
380   //---------------------------------------------------------------------
381   if (fOx==0) {xyz[0]=xyz[1]=xyz[2]=0.; return;}
382   Double_t phi=TMath::ATan2(fOp[0],fOx) + fOalpha;
383   Double_t r=TMath::Sqrt(fOx*fOx + fOp[0]*fOp[0]);
384   xyz[0]=r*TMath::Cos(phi); xyz[1]=r*TMath::Sin(phi); xyz[2]=fOp[1]; 
385 }
386
387 //_______________________________________________________________________
388 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
389   // Returns the array with integrated times for each particle hypothesis
390   for (Int_t i=0; i<kSPECIES; i++) times[i]=fTrackTime[i];
391 }
392
393 //_______________________________________________________________________
394 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
395   // Sets the array with integrated times for each particle hypotesis
396   for (Int_t i=0; i<kSPECIES; i++) fTrackTime[i]=times[i];
397 }
398
399 //_______________________________________________________________________
400 void AliESDtrack::SetITSpid(const Double_t *p) {
401   // Sets values for the probability of each particle type (in ITS)
402   for (Int_t i=0; i<kSPECIES; i++) fITSr[i]=p[i];
403   SetStatus(AliESDtrack::kITSpid);
404 }
405
406 //_______________________________________________________________________
407 void AliESDtrack::GetITSpid(Double_t *p) const {
408   // Gets the probability of each particle type (in ITS)
409   for (Int_t i=0; i<kSPECIES; i++) p[i]=fITSr[i];
410 }
411
412 //_______________________________________________________________________
413 Int_t AliESDtrack::GetITSclusters(UInt_t *idx) const {
414   //---------------------------------------------------------------------
415   // This function returns indices of the assgined ITS clusters 
416   //---------------------------------------------------------------------
417   for (Int_t i=0; i<fITSncls; i++) idx[i]=fITSindex[i];
418   return fITSncls;
419 }
420
421 //_______________________________________________________________________
422 Int_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
423   //---------------------------------------------------------------------
424   // This function returns indices of the assgined ITS clusters 
425   //---------------------------------------------------------------------
426   if (idx!=0)
427     for (Int_t i=0; i<180; i++) idx[i]=fTPCindex[i];  // MI I prefer some constant
428   return fTPCncls;
429 }
430
431 //_______________________________________________________________________
432 void AliESDtrack::SetTPCpid(const Double_t *p) {  
433   // Sets values for the probability of each particle type (in TPC)
434   for (Int_t i=0; i<kSPECIES; i++) fTPCr[i]=p[i];
435   SetStatus(AliESDtrack::kTPCpid);
436 }
437
438 //_______________________________________________________________________
439 void AliESDtrack::GetTPCpid(Double_t *p) const {
440   // Gets the probability of each particle type (in TPC)
441   for (Int_t i=0; i<kSPECIES; i++) p[i]=fTPCr[i];
442 }
443
444 //_______________________________________________________________________
445 Int_t AliESDtrack::GetTRDclusters(UInt_t *idx) const {
446   //---------------------------------------------------------------------
447   // This function returns indices of the assgined TRD clusters 
448   //---------------------------------------------------------------------
449   if (idx!=0)
450     for (Int_t i=0; i<90; i++) idx[i]=fTRDindex[i];  // MI I prefer some constant
451   return fTRDncls;
452 }
453
454 //_______________________________________________________________________
455 void AliESDtrack::SetTRDpid(const Double_t *p) {  
456   // Sets values for the probability of each particle type (in TRD)
457   for (Int_t i=0; i<kSPECIES; i++) fTRDr[i]=p[i];
458   SetStatus(AliESDtrack::kTRDpid);
459 }
460
461 //_______________________________________________________________________
462 void AliESDtrack::GetTRDpid(Double_t *p) const {
463   // Gets the probability of each particle type (in TRD)
464   for (Int_t i=0; i<kSPECIES; i++) p[i]=fTRDr[i];
465 }
466
467 //_______________________________________________________________________
468 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
469 {
470   // Sets the probability of particle type iSpecies to p (in TRD)
471   fTRDr[iSpecies] = p;
472 }
473
474 Float_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
475 {
476   // Returns the probability of particle type iSpecies (in TRD)
477   return fTRDr[iSpecies];
478 }
479
480 //_______________________________________________________________________
481 void AliESDtrack::SetTOFpid(const Double_t *p) {  
482   // Sets the probability of each particle type (in TOF)
483   for (Int_t i=0; i<kSPECIES; i++) fTOFr[i]=p[i];
484   SetStatus(AliESDtrack::kTOFpid);
485 }
486
487 //_______________________________________________________________________
488 void AliESDtrack::GetTOFpid(Double_t *p) const {
489   // Gets probabilities of each particle type (in TOF)
490   for (Int_t i=0; i<kSPECIES; i++) p[i]=fTOFr[i];
491 }
492
493
494
495 //_______________________________________________________________________
496 void AliESDtrack::SetPHOSpid(const Double_t *p) {  
497   // Sets the probability of each particle type (in PHOS)
498   for (Int_t i=0; i<kSPECIES; i++) fPHOSr[i]=p[i];
499   SetStatus(AliESDtrack::kPHOSpid);
500 }
501
502 //_______________________________________________________________________
503 void AliESDtrack::GetPHOSpid(Double_t *p) const {
504   // Gets probabilities of each particle type (in PHOS)
505   for (Int_t i=0; i<kSPECIES; i++) p[i]=fPHOSr[i];
506 }
507
508
509 //_______________________________________________________________________
510 void AliESDtrack::SetRICHpid(const Double_t *p) {  
511   // Sets the probability of each particle type (in RICH)
512   for (Int_t i=0; i<kSPECIES; i++) fRICHr[i]=p[i];
513   SetStatus(AliESDtrack::kRICHpid);
514 }
515
516 //_______________________________________________________________________
517 void AliESDtrack::GetRICHpid(Double_t *p) const {
518   // Gets probabilities of each particle type (in RICH)
519   for (Int_t i=0; i<kSPECIES; i++) p[i]=fRICHr[i];
520 }
521
522
523
524 //_______________________________________________________________________
525 void AliESDtrack::SetESDpid(const Double_t *p) {  
526   // Sets the probability of each particle type for the ESD track
527   for (Int_t i=0; i<kSPECIES; i++) fR[i]=p[i];
528   SetStatus(AliESDtrack::kESDpid);
529 }
530
531 //_______________________________________________________________________
532 void AliESDtrack::GetESDpid(Double_t *p) const {
533   // Gets probability of each particle type for the ESD track
534   for (Int_t i=0; i<kSPECIES; i++) p[i]=fR[i];
535 }
536