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