Data members for the material budget, temporary solution (M.Ivanov)
[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 //           Implementation of the ESD track class
17 //   ESD = Event Summary Data
18 //   This is the class to deal with during the phisics analysis of data
19 //      Origin: Iouri Belikov, CERN
20 //      e-mail: Jouri.Belikov@cern.ch
21 //-----------------------------------------------------------------
22
23 #include "TMath.h"
24
25 #include "AliESDtrack.h"
26 #include "AliKalmanTrack.h"
27 #include "AliLog.h"
28
29 ClassImp(AliESDtrack)
30
31 //_______________________________________________________________________
32 AliESDtrack::AliESDtrack() : 
33   TObject(),
34   fFlags(0),
35   fLabel(0),
36   fID(0),
37   fTrackLength(0),
38   fD(0),
39   fZ(0),
40   fStopVertex(0),
41   fRalpha(0),
42   fRx(0),
43   fCalpha(0),
44   fCx(0),
45   fCchi2(1e10),
46   fIalpha(0),
47   fIx(0),
48   fTalpha(0),
49   fTx(0),
50   fITSchi2(0),
51   fITSncls(0),
52   fITSsignal(0),
53   fITSLabel(0),
54   fITSFakeRatio(0),
55   fITStrack(0),
56   fTPCchi2(0),
57   fTPCncls(0),
58   fTPCClusterMap(159),//number of padrows
59   fTPCsignal(0),
60   fTPCLabel(0),
61   fTRDchi2(0),
62   fTRDncls(0),
63   fTRDncls0(0),
64   fTRDsignal(0),
65   fTRDLabel(0),
66   fTRDQuality(0),
67   fTRDBudget(0),
68   fTRDtrack(0),
69   fTOFchi2(0),
70   fTOFindex(0),
71   fTOFsignal(-1),
72   fPHOSsignal(-1),
73   fEMCALsignal(-1),
74   fRICHchi2(1e10),
75   fRICHncls(0),
76   fRICHindex(0),
77   fRICHsignal(-1),
78   fRICHtheta(0),
79   fRICHphi(0),
80   fRICHdx(0),
81   fRICHdy(0)
82 {
83   //
84   // The default ESD constructor 
85   //
86   for (Int_t i=0; i<AliPID::kSPECIES; i++) {
87     fTrackTime[i]=0.;
88     fR[i]=1.;
89     fITSr[i]=1.;
90     fTPCr[i]=1.;
91     fTRDr[i]=1.;
92     fTOFr[i]=1.;
93     fRICHr[i]=1.;
94   }
95   
96   for (Int_t i=0; i<AliPID::kSPECIESN; i++) {
97     fPHOSr[i]  = 1.;
98     fEMCALr[i] = 1.;
99   }
100
101  
102   fPHOSpos[0]=fPHOSpos[1]=fPHOSpos[2]=0.;
103   fEMCALpos[0]=fEMCALpos[1]=fEMCALpos[2]=0.;
104   Int_t i;
105   for (i=0; i<5; i++)  { 
106     fRp[i]=fCp[i]=fIp[i]=fTp[i]=0.;
107   }
108   for (i=0; i<15; i++) { 
109     fRc[i]=fCc[i]=fIc[i]=fTc[i]=0.;  
110   }
111   for (i=0; i<6; i++)  { fITSindex[i]=0; }
112   for (i=0; i<180; i++){ fTPCindex[i]=0; }
113   for (i=0; i<3;i++)   { fKinkIndexes[i]=0;}
114   for (i=0; i<3;i++)   { fV0Indexes[i]=-1;}
115   for (i=0; i<130; i++) { fTRDindex[i]=0; }
116   for (i=0;i<kNPlane;i++) {fTRDsignals[i]=0.; fTRDTimBin[i]=-1;}
117   for (Int_t i=0;i<4;i++) {fTPCPoints[i]=-1;}
118   for (Int_t i=0;i<3;i++) {fTOFLabel[i]=-1;}
119   for (Int_t i=0;i<10;i++) {fTOFInfo[i]=-1;}
120   fTPCLabel = 0;
121   fTRDLabel = 0;
122   fTRDQuality =0;
123   fTRDBudget =0;
124   fITSLabel = 0;
125   fITStrack = 0;
126   fTRDtrack = 0;  
127 }
128
129 //_______________________________________________________________________
130 AliESDtrack::AliESDtrack(const AliESDtrack& track):
131   TObject(track),
132   fFlags(track.fFlags),
133   fLabel(track.fLabel),
134   fID(track.fID),
135   fTrackLength(track.fTrackLength),
136   fD(track.fD),
137   fZ(track.fZ),
138   fStopVertex(track.fStopVertex),
139   fRalpha(track.fRalpha),
140   fRx(track.fRx),
141   fCalpha(track.fCalpha),
142   fCx(track.fCx),
143   fCchi2(track.fCchi2),
144   fIalpha(track.fIalpha),
145   fIx(track.fIx),
146   fTalpha(track.fTalpha),
147   fTx(track.fTx),
148   fITSchi2(track.fITSchi2),
149   fITSncls(track.fITSncls),
150   fITSsignal(track.fITSsignal),
151   fITSLabel(track.fITSLabel),
152   fITSFakeRatio(track.fITSFakeRatio),
153   fITStrack(0),    //coping separatelly - in user code
154   fTPCchi2(track.fTPCchi2),
155   fTPCncls(track.fTPCncls),
156   fTPCClusterMap(track.fTPCClusterMap),
157   fTPCsignal(track.fTPCsignal),
158   fTPCLabel(track.fTPCLabel),
159   fTRDchi2(track.fTRDchi2),
160   fTRDncls(track.fTRDncls),
161   fTRDncls0(track.fTRDncls0),
162   fTRDsignal(track.fTRDsignal),
163   fTRDLabel(track.fTRDLabel),
164   fTRDQuality(track.fTRDQuality),
165   fTRDBudget(track.fTRDBudget),
166   fTRDtrack(0),
167   fTOFchi2(track.fTOFchi2),
168   fTOFindex(track.fTOFindex),
169   fTOFsignal(track.fTOFsignal),
170   fPHOSsignal(track.fPHOSsignal),
171   fEMCALsignal(track.fEMCALsignal),
172   fRICHchi2(track.fRICHchi2),
173   fRICHncls(track.fRICHncls),
174   fRICHindex(track.fRICHindex),
175   fRICHsignal(track.fRICHsignal),
176   fRICHtheta(track.fRICHtheta),
177   fRICHphi(track.fRICHphi),
178   fRICHdx(track.fRICHdx),
179   fRICHdy(track.fRICHdy)
180 {
181   //
182   //copy constructor
183   //
184   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i] =track.fTrackTime[i];
185   for (Int_t i=0;i<AliPID::kSPECIES;i++)  fR[i] =track.fR[i];
186   //
187   for (Int_t i=0;i<5;i++) fRp[i] =track.fRp[i];
188   for (Int_t i=0;i<15;i++) fRc[i] =track.fRc[i];
189   //
190   for (Int_t i=0;i<5;i++) fCp[i] =track.fCp[i];
191   for (Int_t i=0;i<15;i++)  fCc[i] =track.fCc[i];
192   //
193   for (Int_t i=0;i<5;i++) fIp[i] =track.fIp[i];
194   for (Int_t i=0;i<15;i++)  fIc[i] =track.fIc[i];
195   //
196   for (Int_t i=0;i<5;i++) fTp[i] =track.fTp[i];
197   for (Int_t i=0;i<15;i++)  fTc[i] =track.fTc[i];
198   //
199   for (Int_t i=0;i<12;i++) fITSchi2MIP[i] =track.fITSchi2MIP[i];
200   for (Int_t i=0;i<6;i++) fITSindex[i]=track.fITSindex[i];    
201   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=track.fITSr[i]; 
202   //
203   for (Int_t i=0;i<180;i++) fTPCindex[i]=track.fTPCindex[i];  
204   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=track.fTPCr[i]; 
205   for (Int_t i=0;i<4;i++) {fTPCPoints[i]=track.fTPCPoints[i];}
206   for (Int_t i=0; i<3;i++)   { fKinkIndexes[i]=track.fKinkIndexes[i];}
207   for (Int_t i=0; i<3;i++)   { fV0Indexes[i]=track.fV0Indexes[i];}
208   //
209   for (Int_t i=0;i<130;i++) fTRDindex[i]=track.fTRDindex[i];   
210   for (Int_t i=0;i<kNPlane;i++) {
211       fTRDsignals[i]=track.fTRDsignals[i]; 
212       fTRDTimBin[i]=track.fTRDTimBin[i];
213   }
214   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i]=track.fTRDr[i]; 
215   //
216   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i]=track.fTOFr[i];
217   for (Int_t i=0;i<3;i++) fTOFLabel[i]=track.fTOFLabel[i];
218   for (Int_t i=0;i<10;i++) fTOFInfo[i]=track.fTOFInfo[i];
219   //
220   for (Int_t i=0;i<3;i++) fPHOSpos[i]=track.fPHOSpos[i]; 
221   for (Int_t i=0;i<AliPID::kSPECIESN;i++) fPHOSr[i]=track.fPHOSr[i]; 
222   //
223   for (Int_t i=0;i<3;i++) fEMCALpos[i]=track.fEMCALpos[i]; 
224   for (Int_t i=0;i<AliPID::kSPECIESN;i++) fEMCALr[i]=track.fEMCALr[i]; 
225   //
226   for (Int_t i=0;i<AliPID::kSPECIES;i++) fRICHr[i]=track.fRICHr[i];
227 }
228 //_______________________________________________________________________
229 AliESDtrack::~AliESDtrack(){ 
230   //
231   // This is destructor according Coding Conventrions 
232   //
233   //printf("Delete track\n");
234   delete fITStrack;
235   delete fTRDtrack;  
236 }
237
238 //_______________________________________________________________________
239 void AliESDtrack::MakeMiniESDtrack(){
240   // Resets everything except
241   // fFlags: Reconstruction status flags 
242   // fLabel: Track label
243   // fID:  Unique ID of the track
244   // fD: Impact parameter in XY-plane
245   // fZ: Impact parameter in Z 
246   // fR[AliPID::kSPECIES]: combined "detector response probability"
247   // Running track parameters
248   // fRalpha: track rotation angle
249   // fRx: X-coordinate of the track reference plane 
250   // fRp[5]: external track parameters  
251   // fRc[15]: external cov. matrix of the track parameters
252   
253   fTrackLength = 0;
254   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i] = 0;
255   fStopVertex = 0;
256
257   // Reset track parameters constrained to the primary vertex
258   fCalpha = 0;
259   fCx = 0;
260   for (Int_t i=0;i<5;i++) fCp[i] = 0;
261   for (Int_t i=0;i<15;i++)  fCc[i] = 0;
262   fCchi2 = 0;
263
264   // Reset track parameters at the inner wall of TPC
265   fIalpha = 0;
266   fIx = 0;
267   for (Int_t i=0;i<5;i++) fIp[i] = 0;
268   for (Int_t i=0;i<15;i++)  fIc[i] = 0;
269
270   // Reset track parameters at the inner wall of the TRD
271   fTalpha = 0;
272   fTx = 0;
273   for (Int_t i=0;i<5;i++) fTp[i] = 0;
274   for (Int_t i=0;i<15;i++)  fTc[i] = 0;
275
276   // Reset ITS track related information
277   fITSchi2 = 0;
278   for (Int_t i=0;i<12;i++) fITSchi2MIP[i] = 0;
279   fITSncls = 0;       
280   for (Int_t i=0;i<6;i++) fITSindex[i]= 0;    
281   fITSsignal = 0;     
282   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]= 0; 
283   fITSLabel = 0;       
284   fITSFakeRatio = 0;   
285   fITStrack =0;
286
287   // Reset TPC related track information
288   fTPCchi2 = 0;       
289   fTPCncls = 0;       
290   for (Int_t i=0;i<180;i++) fTPCindex[i] = 0;  
291   fTPCClusterMap = 0;  
292   fTPCsignal= 0;      
293   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0; 
294   fTPCLabel=0;       
295   for (Int_t i=0;i<4;i++) fTPCPoints[i] = 0;
296   for (Int_t i=0; i<3;i++)   fKinkIndexes[i] = 0;
297   for (Int_t i=0; i<3;i++)   fV0Indexes[i] = 0;
298
299   // Reset TRD related track information
300   fTRDchi2 = 0;        
301   fTRDncls = 0;       
302   fTRDncls0 = 0;       
303   for (Int_t i=0;i<130;i++) fTRDindex[i] = 0;   
304   fTRDsignal = 0;      
305   for (Int_t i=0;i<kNPlane;i++) {
306       fTRDsignals[i] = 0; 
307       fTRDTimBin[i]  = 0;
308   }
309   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0; 
310   fTRDLabel = 0;       
311   fTRDtrack = 0; 
312   fTRDQuality  = 0;
313   fTRDBudget  = 0;
314
315   // Reset TOF related track information
316   fTOFchi2 = 0;        
317   fTOFindex = 0;       
318   fTOFsignal = 0;      
319   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i] = 0;
320   for (Int_t i=0;i<3;i++) fTOFLabel[i] = 0;
321   for (Int_t i=0;i<10;i++) fTOFInfo[i] = 0;
322
323   // Reset PHOS related track information
324   for (Int_t i=0;i<3;i++) fPHOSpos[i] = 0; 
325   fPHOSsignal = 0; 
326   for (Int_t i=0;i<AliPID::kSPECIESN;i++) fPHOSr[i] = 0;
327  
328   // Reset EMCAL related track information
329   for (Int_t i=0;i<3;i++) fEMCALpos[i] = 0; 
330   fEMCALsignal = 0; 
331   for (Int_t i=0;i<AliPID::kSPECIESN;i++) fEMCALr[i] = 0;
332  
333   // Reset RICH related track information
334   fRICHchi2 = 0;     
335   fRICHncls = 0;     
336   fRICHindex = 0;     
337   fRICHsignal = 0;     
338   for (Int_t i=0;i<AliPID::kSPECIES;i++) fRICHr[i] = 0;
339   fRICHtheta = 0;     
340   fRICHphi = 0;      
341   fRICHdx = 0;     
342   fRICHdy = 0;      
343
344
345 //_______________________________________________________________________
346 Double_t AliESDtrack::GetMass() const {
347   // Returns the mass of the most probable particle type
348   Float_t max=0.;
349   Int_t k=-1;
350   for (Int_t i=0; i<AliPID::kSPECIES; i++) {
351     if (fR[i]>max) {k=i; max=fR[i];}
352   }
353   if (k==0) { // dE/dx "crossing points" in the TPC
354      Double_t p=GetP();
355      if ((p>0.38)&&(p<0.48))
356         if (fR[0]<fR[3]*10.) return AliPID::ParticleMass(AliPID::kKaon);
357      if ((p>0.75)&&(p<0.85))
358         if (fR[0]<fR[4]*10.) return AliPID::ParticleMass(AliPID::kProton);
359      return 0.00051;
360   }
361   if (k==1) return AliPID::ParticleMass(AliPID::kMuon); 
362   if (k==2||k==-1) return AliPID::ParticleMass(AliPID::kPion);
363   if (k==3) return AliPID::ParticleMass(AliPID::kKaon);
364   if (k==4) return AliPID::ParticleMass(AliPID::kProton);
365   AliWarning("Undefined mass !");
366   return AliPID::ParticleMass(AliPID::kPion);
367 }
368
369 //_______________________________________________________________________
370 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags) {
371   //
372   // This function updates track's running parameters 
373   //
374   Bool_t rc=kTRUE;
375
376   SetStatus(flags);
377   fLabel=t->GetLabel();
378
379   if (t->IsStartedTimeIntegral()) {
380     SetStatus(kTIME);
381     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
382     SetIntegratedLength(t->GetIntegratedLength());
383   }
384
385   fRalpha=t->GetAlpha();
386   t->GetExternalParameters(fRx,fRp);
387   t->GetExternalCovariance(fRc);
388
389   switch (flags) {
390     
391   case kITSin: case kITSout: case kITSrefit:
392     fITSncls=t->GetNumberOfClusters();
393     fITSchi2=t->GetChi2();
394     for (Int_t i=0;i<fITSncls;i++) fITSindex[i]=t->GetClusterIndex(i);
395     fITSsignal=t->GetPIDsignal();
396     fITSLabel = t->GetLabel();
397     fITSFakeRatio = t->GetFakeRatio();
398     break;
399     
400   case kTPCin: case kTPCrefit:
401     fTPCLabel = t->GetLabel();
402     fIalpha=fRalpha;
403     fIx=fRx;    
404     {
405       Int_t i;
406       for (i=0; i<5; i++) fIp[i]=fRp[i];
407       for (i=0; i<15;i++) fIc[i]=fRc[i];
408     }
409   case kTPCout:
410   
411     fTPCncls=t->GetNumberOfClusters();
412     fTPCchi2=t->GetChi2();
413     
414      {//prevrow must be declared in separate namespace, otherwise compiler cries:
415       //"jump to case label crosses initialization of `Int_t prevrow'"
416        Int_t prevrow = -1;
417        //       for (Int_t i=0;i<fTPCncls;i++) 
418        for (Int_t i=0;i<160;i++) 
419         {
420           fTPCindex[i]=t->GetClusterIndex(i);
421
422           // Piotr's Cluster Map for HBT  
423           // ### please change accordingly if cluster array is changing 
424           // to "New TPC Tracking" style (with gaps in array) 
425           Int_t idx = fTPCindex[i];
426           Int_t sect = (idx&0xff000000)>>24;
427           Int_t row = (idx&0x00ff0000)>>16;
428           if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
429
430           fTPCClusterMap.SetBitNumber(row,kTRUE);
431
432           //Fill the gap between previous row and this row with 0 bits
433           //In case  ###  pleas change it as well - just set bit 0 in case there 
434           //is no associated clusters for current "i"
435           if (prevrow < 0) 
436            {
437              prevrow = row;//if previous bit was not assigned yet == this is the first one
438            }
439           else
440            { //we don't know the order (inner to outer or reverse)
441              //just to be save in case it is going to change
442              Int_t n = 0, m = 0;
443              if (prevrow < row)
444               {
445                 n = prevrow;
446                 m = row;
447               }
448              else
449               {
450                 n = row;
451                 m = prevrow;
452               }
453
454              for (Int_t j = n+1; j < m; j++)
455               {
456                 fTPCClusterMap.SetBitNumber(j,kFALSE);
457               }
458              prevrow = row; 
459            }
460           // End Of Piotr's Cluster Map for HBT
461         }
462      }
463     fTPCsignal=t->GetPIDsignal();
464     {Double_t mass=t->GetMass();    // preliminary mass setting 
465     if (mass>0.5) fR[4]=1.;         //        used by
466     else if (mass<0.4) fR[2]=1.;    // the ITS reconstruction
467     else fR[3]=1.;}
468                      //
469     break;
470
471   case kTRDout: case kTRDin: case kTRDrefit:
472     fTRDLabel = t->GetLabel(); 
473     fTRDncls=t->GetNumberOfClusters();
474     fTRDchi2=t->GetChi2();
475     for (Int_t i=0;i<fTRDncls;i++) fTRDindex[i]=t->GetClusterIndex(i);
476     fTRDsignal=t->GetPIDsignal();
477     break;
478   case kTRDbackup:
479     t->GetExternalParameters(fTx,fTp);
480     t->GetExternalCovariance(fTc);
481     fTalpha = t->GetAlpha();
482     fTRDncls0 = t->GetNumberOfClusters(); 
483     break;
484   case kTOFin: 
485     break;
486   case kTOFout: 
487     break;
488   case kTRDStop:
489     break;
490   default: 
491     AliError("Wrong flag !");
492     return kFALSE;
493   }
494
495   return rc;
496 }
497
498 //_______________________________________________________________________
499 void 
500 AliESDtrack::SetConstrainedTrackParams(const AliKalmanTrack *t, Double_t chi2) {
501   //
502   // This function sets the constrained track parameters 
503   //
504   Int_t i;
505   Double_t x,buf[15];
506   fCalpha=t->GetAlpha();
507   t->GetExternalParameters(x,buf); fCx=x;
508   for (i=0; i<5; i++) fCp[i]=buf[i];
509   t->GetExternalCovariance(buf);
510   for (i=0; i<15; i++) fCc[i]=buf[i];
511   fCchi2=chi2;
512 }
513
514
515 //_______________________________________________________________________
516 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
517   //---------------------------------------------------------------------
518   // This function returns external representation of the track parameters
519   //---------------------------------------------------------------------
520   x=fRx;
521   for (Int_t i=0; i<5; i++) p[i]=fRp[i];
522 }
523
524 //_______________________________________________________________________
525 Bool_t AliESDtrack::
526 GetExternalParametersAt(Double_t x, Double_t b, Double_t p[5]) const {
527   //---------------------------------------------------------------------
528   // This function returns external track parameters extrapolated to
529   // the radial position "x" (cm) in the magnetic field "b" (kG)  
530   //---------------------------------------------------------------------
531   Double_t convconst=0.299792458*b/1000.;
532   Double_t dx=x-fRx;
533   Double_t f1=fRp[2], f2=f1 + dx*fRp[4]*convconst;
534
535   if (TMath::Abs(f2) >= 0.9999) return kFALSE;
536   
537   Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
538   p[0] = fRp[0] + dx*(f1+f2)/(r1+r2);
539   p[1] = fRp[1] + dx*(f1+f2)/(f1*r2 + f2*r1)*fRp[3];
540   p[2] = f2;
541   p[3] = fRp[3];
542   p[4] = fRp[4];
543
544   return kTRUE;
545 }
546
547 //_______________________________________________________________________
548 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
549   //---------------------------------------------------------------------
550   // This function returns external representation of the cov. matrix
551   //---------------------------------------------------------------------
552   for (Int_t i=0; i<15; i++) cov[i]=fRc[i];
553 }
554
555
556 //_______________________________________________________________________
557 void 
558 AliESDtrack::GetConstrainedExternalParameters(Double_t &x, Double_t p[5])const{
559   //---------------------------------------------------------------------
560   // This function returns the constrained external track parameters
561   //---------------------------------------------------------------------
562   x=fCx;
563   for (Int_t i=0; i<5; i++) p[i]=fCp[i];
564 }
565 //_______________________________________________________________________
566 void 
567 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
568   //---------------------------------------------------------------------
569   // This function returns the constrained external cov. matrix
570   //---------------------------------------------------------------------
571   for (Int_t i=0; i<15; i++) c[i]=fCc[i];
572 }
573
574
575 Double_t AliESDtrack::GetP() const {
576   //---------------------------------------------------------------------
577   // This function returns the track momentum
578   // Results for (nearly) straight tracks are meaningless !
579   //---------------------------------------------------------------------
580   if (TMath::Abs(fRp[4])<=0) return 0;
581   Double_t pt=1./TMath::Abs(fRp[4]);
582   return pt*TMath::Sqrt(1.+ fRp[3]*fRp[3]);
583 }
584
585 //_______________________________________________________________________
586 Double_t AliESDtrack::GetD(Double_t b, Double_t x, Double_t y) const {
587   //------------------------------------------------------------------
588   // This function calculates the transverse impact parameter
589   // with respect to a point with global coordinates (x,y)
590   // in the magnetic field "b" (kG)
591   //------------------------------------------------------------------
592   Double_t convconst=0.299792458*b/1000.;
593   Double_t rp4=fRp[4]*convconst;
594
595   Double_t xt=fRx, yt=fRp[0];
596
597   Double_t sn=TMath::Sin(fRalpha), cs=TMath::Cos(fRalpha);
598   Double_t a = x*cs + y*sn;
599   y = -x*sn + y*cs; x=a;
600   xt-=x; yt-=y;
601
602   sn=rp4*xt - fRp[2]; cs=rp4*yt + TMath::Sqrt(1.- fRp[2]*fRp[2]);
603   a=2*(xt*fRp[2] - yt*TMath::Sqrt(1.- fRp[2]*fRp[2]))-rp4*(xt*xt + yt*yt);
604   if (rp4<0) a=-a;
605   return a/(1 + TMath::Sqrt(sn*sn + cs*cs));
606 }
607
608 Bool_t Local2GlobalMomentum(Double_t p[3],Double_t alpha) {
609   //----------------------------------------------------------------
610   // This function performs local->global transformation of the
611   // track momentum.
612   // When called, the arguments are:
613   //    p[0] = 1/pt of the track;
614   //    p[1] = sine of local azim. angle of the track momentum;
615   //    p[2] = tangent of the track momentum dip angle;
616   //   alpha - rotation angle. 
617   // The result is returned as:
618   //    p[0] = px
619   //    p[1] = py
620   //    p[2] = pz
621   // Results for (nearly) straight tracks are meaningless !
622   //----------------------------------------------------------------
623   if (TMath::Abs(p[0])<=0)        return kFALSE;
624   if (TMath::Abs(p[1])> 0.999999) return kFALSE;
625
626   Double_t pt=1./TMath::Abs(p[0]);
627   Double_t cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
628   Double_t r=TMath::Sqrt(1 - p[1]*p[1]);
629   p[0]=pt*(r*cs - p[1]*sn); p[1]=pt*(p[1]*cs + r*sn); p[2]=pt*p[2];
630
631   return kTRUE;
632 }
633
634 Bool_t Local2GlobalPosition(Double_t r[3],Double_t alpha) {
635   //----------------------------------------------------------------
636   // This function performs local->global transformation of the
637   // track position.
638   // When called, the arguments are:
639   //    r[0] = local x
640   //    r[1] = local y
641   //    r[2] = local z
642   //   alpha - rotation angle. 
643   // The result is returned as:
644   //    r[0] = global x
645   //    r[1] = global y
646   //    r[2] = global z
647   //----------------------------------------------------------------
648   Double_t cs=TMath::Cos(alpha), sn=TMath::Sin(alpha), x=r[0];
649   r[0]=x*cs - r[1]*sn; r[1]=x*sn + r[1]*cs;
650
651   return kTRUE;
652 }
653
654 Bool_t AliESDtrack::GetConstrainedPxPyPz(Double_t *p) const {
655   //---------------------------------------------------------------------
656   // This function returns the constrained global track momentum components
657   // Results for (nearly) straight tracks are meaningless !
658   //---------------------------------------------------------------------
659   p[0]=fCp[4]; p[1]=fCp[2]; p[2]=fCp[3];
660   return Local2GlobalMomentum(p,fCalpha);
661 }  
662
663 Bool_t AliESDtrack::GetConstrainedXYZ(Double_t *r) const {
664   //---------------------------------------------------------------------
665   // This function returns the constrained global track position
666   //---------------------------------------------------------------------
667   r[0]=fCx; r[1]=fCp[0]; r[2]=fCp[1];
668   return Local2GlobalPosition(r,fCalpha);
669 }
670
671 Bool_t AliESDtrack::GetPxPyPz(Double_t *p) const {
672   //---------------------------------------------------------------------
673   // This function returns the global track momentum components
674   // Results for (nearly) straight tracks are meaningless !
675   //---------------------------------------------------------------------
676   p[0]=fRp[4]; p[1]=fRp[2]; p[2]=fRp[3];
677   return Local2GlobalMomentum(p,fRalpha);
678 }
679
680 Bool_t AliESDtrack::GetXYZ(Double_t *r) const {
681   //---------------------------------------------------------------------
682   // This function returns the global track position
683   //---------------------------------------------------------------------
684   r[0]=fRx; r[1]=fRp[0]; r[2]=fRp[1];
685   return Local2GlobalPosition(r,fRalpha);
686 }
687
688 void AliESDtrack::GetCovariance(Double_t cv[21]) const {
689   //---------------------------------------------------------------------
690   // This function returns the global covariance matrix of the track params
691   // 
692   // Cov(x,x) ... :   cv[0]
693   // Cov(y,x) ... :   cv[1]  cv[2]
694   // Cov(z,x) ... :   cv[3]  cv[4]  cv[5]
695   // Cov(px,x)... :   cv[6]  cv[7]  cv[8]  cv[9]
696   // Cov(py,x)... :   cv[10] cv[11] cv[12] cv[13] cv[14]
697   // Cov(pz,x)... :   cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
698   //
699   // Results for (nearly) straight tracks are meaningless !
700   //---------------------------------------------------------------------
701   if (TMath::Abs(fRp[4])<=0) {
702      for (Int_t i=0; i<21; i++) cv[i]=0.;
703      return;
704   }
705   if (TMath::Abs(fRp[2]) > 0.999999) {
706      for (Int_t i=0; i<21; i++) cv[i]=0.;
707      return;
708   }
709   Double_t pt=1./TMath::Abs(fRp[4]);
710   Double_t cs=TMath::Cos(fRalpha), sn=TMath::Sin(fRalpha);
711   Double_t r=TMath::Sqrt(1-fRp[2]*fRp[2]);
712
713   Double_t m00=-sn, m10=cs;
714   Double_t m23=-pt*(sn + fRp[2]*cs/r), m43=-pt*pt*(r*cs - fRp[2]*sn);
715   Double_t m24= pt*(cs - fRp[2]*sn/r), m44=-pt*pt*(r*sn + fRp[2]*cs);
716   Double_t m35=pt, m45=-pt*pt*fRp[3];
717
718   cv[0]=fRc[0]*m00*m00;
719   cv[1]=fRc[0]*m00*m10; 
720   cv[2]=fRc[0]*m10*m10;
721   cv[3]=fRc[1]*m00; 
722   cv[4]=fRc[1]*m10; 
723   cv[5]=fRc[2];
724   cv[6]=m00*(fRc[3]*m23+fRc[10]*m43); 
725   cv[7]=m10*(fRc[3]*m23+fRc[10]*m43); 
726   cv[8]=fRc[4]*m23+fRc[11]*m43; 
727   cv[9]=m23*(fRc[5]*m23+fRc[12]*m43)+m43*(fRc[12]*m23+fRc[14]*m43);
728   cv[10]=m00*(fRc[3]*m24+fRc[10]*m44); 
729   cv[11]=m10*(fRc[3]*m24+fRc[10]*m44); 
730   cv[12]=fRc[4]*m24+fRc[11]*m44; 
731   cv[13]=m23*(fRc[5]*m24+fRc[12]*m44)+m43*(fRc[12]*m24+fRc[14]*m44);
732   cv[14]=m24*(fRc[5]*m24+fRc[12]*m44)+m44*(fRc[12]*m24+fRc[14]*m44);
733   cv[15]=m00*(fRc[6]*m35+fRc[10]*m45); 
734   cv[16]=m10*(fRc[6]*m35+fRc[10]*m45); 
735   cv[17]=fRc[7]*m35+fRc[11]*m45; 
736   cv[18]=m23*(fRc[8]*m35+fRc[12]*m45)+m43*(fRc[13]*m35+fRc[14]*m45);
737   cv[19]=m24*(fRc[8]*m35+fRc[12]*m45)+m44*(fRc[13]*m35+fRc[14]*m45); 
738   cv[20]=m35*(fRc[9]*m35+fRc[13]*m45)+m45*(fRc[13]*m35+fRc[14]*m45);
739 }
740
741 Bool_t AliESDtrack::GetInnerPxPyPz(Double_t *p) const {
742   //---------------------------------------------------------------------
743   // This function returns the global track momentum components
744   // af the entrance of the TPC
745   //---------------------------------------------------------------------
746   p[0]=fIp[4]; p[1]=fIp[2]; p[2]=fIp[3];
747   return Local2GlobalMomentum(p,fIalpha);
748 }
749
750 Bool_t AliESDtrack::GetInnerXYZ(Double_t *r) const {
751   //---------------------------------------------------------------------
752   // This function returns the global track position
753   // af the entrance of the TPC
754   //---------------------------------------------------------------------
755   if (fIx==0) return kFALSE;
756   r[0]=fIx; r[1]=fIp[0]; r[2]=fIp[1];
757   return Local2GlobalPosition(r,fIalpha);
758 }
759
760 void AliESDtrack::GetInnerExternalParameters(Double_t &x, Double_t p[5]) const 
761 {
762   //skowron
763  //---------------------------------------------------------------------
764   // This function returns external representation of the track parameters at Inner Layer of TPC
765   //---------------------------------------------------------------------
766   x=fIx;
767   for (Int_t i=0; i<5; i++) p[i]=fIp[i];
768 }
769 void AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const
770 {
771  //skowron
772  //---------------------------------------------------------------------
773  // This function returns external representation of the cov. matrix at Inner Layer of TPC
774  //---------------------------------------------------------------------
775  for (Int_t i=0; i<15; i++) cov[i]=fIc[i];
776  
777 }
778
779 void  AliESDtrack::GetTRDExternalParameters(Double_t &x, Double_t&alpha, Double_t p[5], Double_t cov[15]) const
780 {
781   //
782   //this function returns TRD parameters
783   //
784   x=fTx;
785   alpha = fTalpha; 
786   for (Int_t i=0; i<5; i++) p[i]=fTp[i];
787   for (Int_t i=0; i<15; i++) cov[i]=fTc[i];
788 }
789
790 Bool_t AliESDtrack::GetPxPyPzAt(Double_t x, Double_t b, Double_t *p) const {
791   //---------------------------------------------------------------------
792   // This function returns the global track momentum extrapolated to
793   // the radial position "x" (cm) in the magnetic field "b" (kG)
794   //---------------------------------------------------------------------
795   Double_t convconst=0.299792458*b/1000.;
796   p[0]=fRp[4]; 
797   p[1]=fRp[2]+(x-fRx)*fRp[4]*convconst; 
798   p[2]=fRp[3];
799   return Local2GlobalMomentum(p,fRalpha);
800 }
801
802 Bool_t AliESDtrack::GetXYZAt(Double_t x, Double_t b, Double_t *r) const {
803   //---------------------------------------------------------------------
804   // This function returns the global track position extrapolated to
805   // the radial position "x" (cm) in the magnetic field "b" (kG)
806   //---------------------------------------------------------------------
807   Double_t convconst=0.299792458*b/1000.;
808   Double_t dx=x-fRx;
809   Double_t f1=fRp[2], f2=f1 + dx*fRp[4]*convconst;
810
811   if (TMath::Abs(f2) >= 0.9999) return kFALSE;
812   
813   Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
814   r[0] = x;
815   r[1] = fRp[0] + dx*(f1+f2)/(r1+r2);
816   r[2] = fRp[1] + dx*(f1+f2)/(f1*r2 + f2*r1)*fRp[3];
817   return Local2GlobalPosition(r,fRalpha);
818 }
819
820 //_______________________________________________________________________
821 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
822   // Returns the array with integrated times for each particle hypothesis
823   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
824 }
825
826 //_______________________________________________________________________
827 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
828   // Sets the array with integrated times for each particle hypotesis
829   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
830 }
831
832 //_______________________________________________________________________
833 void AliESDtrack::SetITSpid(const Double_t *p) {
834   // Sets values for the probability of each particle type (in ITS)
835   for (Int_t i=0; i<AliPID::kSPECIES; i++) fITSr[i]=p[i];
836   SetStatus(AliESDtrack::kITSpid);
837 }
838
839 void AliESDtrack::SetITSChi2MIP(const Float_t *chi2mip){
840   for (Int_t i=0; i<12; i++) fITSchi2MIP[i]=chi2mip[i];
841 }
842 //_______________________________________________________________________
843 void AliESDtrack::GetITSpid(Double_t *p) const {
844   // Gets the probability of each particle type (in ITS)
845   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
846 }
847
848 //_______________________________________________________________________
849 Int_t AliESDtrack::GetITSclusters(UInt_t *idx) const {
850   //---------------------------------------------------------------------
851   // This function returns indices of the assgined ITS clusters 
852   //---------------------------------------------------------------------
853   for (Int_t i=0; i<fITSncls; i++) idx[i]=fITSindex[i];
854   return fITSncls;
855 }
856
857 //_______________________________________________________________________
858 Int_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
859   //---------------------------------------------------------------------
860   // This function returns indices of the assgined ITS clusters 
861   //---------------------------------------------------------------------
862   if (idx!=0)
863     for (Int_t i=0; i<180; i++) idx[i]=fTPCindex[i];  // MI I prefer some constant
864   return fTPCncls;
865 }
866
867 Float_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
868   //
869   // GetDensity of the clusters on given region between row0 and row1
870   // Dead zone effect takin into acoount
871   //
872   Int_t good  = 0;
873   Int_t found = 0;
874   //  
875   for (Int_t i=row0;i<=row1;i++){     
876     Int_t index = fTPCindex[i];
877     if (index!=-1)  good++;             // track outside of dead zone
878     if (index>0)    found++;
879   }
880   Float_t density=0.5;
881   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
882   return density;
883 }
884
885 //_______________________________________________________________________
886 void AliESDtrack::SetTPCpid(const Double_t *p) {  
887   // Sets values for the probability of each particle type (in TPC)
888   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTPCr[i]=p[i];
889   SetStatus(AliESDtrack::kTPCpid);
890 }
891
892 //_______________________________________________________________________
893 void AliESDtrack::GetTPCpid(Double_t *p) const {
894   // Gets the probability of each particle type (in TPC)
895   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
896 }
897
898 //_______________________________________________________________________
899 Int_t AliESDtrack::GetTRDclusters(UInt_t *idx) const {
900   //---------------------------------------------------------------------
901   // This function returns indices of the assgined TRD clusters 
902   //---------------------------------------------------------------------
903   if (idx!=0)
904     for (Int_t i=0; i<130; i++) idx[i]=fTRDindex[i];  // MI I prefer some constant
905   return fTRDncls;
906 }
907
908 //_______________________________________________________________________
909 void AliESDtrack::SetTRDpid(const Double_t *p) {  
910   // Sets values for the probability of each particle type (in TRD)
911   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTRDr[i]=p[i];
912   SetStatus(AliESDtrack::kTRDpid);
913 }
914
915 //_______________________________________________________________________
916 void AliESDtrack::GetTRDpid(Double_t *p) const {
917   // Gets the probability of each particle type (in TRD)
918   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
919 }
920
921 //_______________________________________________________________________
922 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
923 {
924   // Sets the probability of particle type iSpecies to p (in TRD)
925   fTRDr[iSpecies] = p;
926 }
927
928 Float_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
929 {
930   // Returns the probability of particle type iSpecies (in TRD)
931   return fTRDr[iSpecies];
932 }
933
934 //_______________________________________________________________________
935 void AliESDtrack::SetTOFpid(const Double_t *p) {  
936   // Sets the probability of each particle type (in TOF)
937   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTOFr[i]=p[i];
938   SetStatus(AliESDtrack::kTOFpid);
939 }
940
941 //_______________________________________________________________________
942 void AliESDtrack::SetTOFLabel(const Int_t *p) {  
943   // Sets  (in TOF)
944   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
945 }
946
947 //_______________________________________________________________________
948 void AliESDtrack::GetTOFpid(Double_t *p) const {
949   // Gets probabilities of each particle type (in TOF)
950   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
951 }
952
953 //_______________________________________________________________________
954 void AliESDtrack::GetTOFLabel(Int_t *p) const {
955   // Gets (in TOF)
956   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
957 }
958
959 //_______________________________________________________________________
960 void AliESDtrack::GetTOFInfo(Float_t *info) const {
961   // Gets (in TOF)
962   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
963 }
964
965 //_______________________________________________________________________
966 void AliESDtrack::SetTOFInfo(Float_t*info) {
967   // Gets (in TOF)
968   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
969 }
970
971
972
973 //_______________________________________________________________________
974 void AliESDtrack::SetPHOSpid(const Double_t *p) {  
975   // Sets the probability of each particle type (in PHOS)
976   for (Int_t i=0; i<AliPID::kSPECIESN; i++) fPHOSr[i]=p[i];
977   SetStatus(AliESDtrack::kPHOSpid);
978 }
979
980 //_______________________________________________________________________
981 void AliESDtrack::GetPHOSpid(Double_t *p) const {
982   // Gets probabilities of each particle type (in PHOS)
983   for (Int_t i=0; i<AliPID::kSPECIESN; i++) p[i]=fPHOSr[i];
984 }
985
986 //_______________________________________________________________________
987 void AliESDtrack::SetEMCALpid(const Double_t *p) {  
988   // Sets the probability of each particle type (in EMCAL)
989   for (Int_t i=0; i<AliPID::kSPECIESN; i++) fEMCALr[i]=p[i];
990   SetStatus(AliESDtrack::kEMCALpid);
991 }
992
993 //_______________________________________________________________________
994 void AliESDtrack::GetEMCALpid(Double_t *p) const {
995   // Gets probabilities of each particle type (in EMCAL)
996   for (Int_t i=0; i<AliPID::kSPECIESN; i++) p[i]=fEMCALr[i];
997 }
998
999 //_______________________________________________________________________
1000 void AliESDtrack::SetRICHpid(const Double_t *p) {  
1001   // Sets the probability of each particle type (in RICH)
1002   for (Int_t i=0; i<AliPID::kSPECIES; i++) fRICHr[i]=p[i];
1003   SetStatus(AliESDtrack::kRICHpid);
1004 }
1005
1006 //_______________________________________________________________________
1007 void AliESDtrack::GetRICHpid(Double_t *p) const {
1008   // Gets probabilities of each particle type (in RICH)
1009   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fRICHr[i];
1010 }
1011
1012
1013
1014 //_______________________________________________________________________
1015 void AliESDtrack::SetESDpid(const Double_t *p) {  
1016   // Sets the probability of each particle type for the ESD track
1017   for (Int_t i=0; i<AliPID::kSPECIES; i++) fR[i]=p[i];
1018   SetStatus(AliESDtrack::kESDpid);
1019 }
1020
1021 //_______________________________________________________________________
1022 void AliESDtrack::GetESDpid(Double_t *p) const {
1023   // Gets probability of each particle type for the ESD track
1024   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
1025 }
1026
1027 //_______________________________________________________________________
1028 void AliESDtrack::Print(Option_t *) const {
1029   // Prints info on the track
1030   
1031   printf("ESD track info\n") ; 
1032   Double_t p[AliPID::kSPECIESN] ; 
1033   Int_t index = 0 ; 
1034   if( IsOn(kITSpid) ){
1035     printf("From ITS: ") ; 
1036     GetITSpid(p) ; 
1037     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1038       printf("%f, ", p[index]) ;
1039     printf("\n           signal = %f\n", GetITSsignal()) ;
1040   } 
1041   if( IsOn(kTPCpid) ){
1042     printf("From TPC: ") ; 
1043     GetTPCpid(p) ; 
1044     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1045       printf("%f, ", p[index]) ;
1046     printf("\n           signal = %f\n", GetTPCsignal()) ;
1047   }
1048   if( IsOn(kTRDpid) ){
1049     printf("From TRD: ") ; 
1050     GetTRDpid(p) ; 
1051     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1052       printf("%f, ", p[index]) ;
1053     printf("\n           signal = %f\n", GetTRDsignal()) ;
1054   }
1055   if( IsOn(kTOFpid) ){
1056     printf("From TOF: ") ; 
1057     GetTOFpid(p) ; 
1058     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1059       printf("%f, ", p[index]) ;
1060     printf("\n           signal = %f\n", GetTOFsignal()) ;
1061   }
1062   if( IsOn(kRICHpid) ){
1063     printf("From RICH: ") ; 
1064     GetRICHpid(p) ; 
1065     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1066       printf("%f, ", p[index]) ;
1067     printf("\n           signal = %f\n", GetRICHsignal()) ;
1068   }
1069   if( IsOn(kPHOSpid) ){
1070     printf("From PHOS: ") ; 
1071     GetPHOSpid(p) ; 
1072     for(index = 0 ; index < AliPID::kSPECIESN; index++) 
1073       printf("%f, ", p[index]) ;
1074     printf("\n           signal = %f\n", GetPHOSsignal()) ;
1075   }
1076   if( IsOn(kEMCALpid) ){
1077     printf("From EMCAL: ") ; 
1078     GetEMCALpid(p) ; 
1079     for(index = 0 ; index < AliPID::kSPECIESN; index++) 
1080       printf("%f, ", p[index]) ;
1081     printf("\n           signal = %f\n", GetEMCALsignal()) ;
1082   }
1083