coverity fix
[u/mrichter/AliRoot.git] / TPC / AliTPCtrack.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 /* $Id$ */
17
18 //-----------------------------------------------------------------
19 //           Implementation of the TPC track class
20 //        This class is used by the AliTPCtracker class
21 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
22 //-----------------------------------------------------------------
23
24 #include <Riostream.h>
25
26 #include "AliTPCtrack.h"
27 #include "AliCluster.h"
28 #include "AliTracker.h"
29 #include "AliESDtrack.h"
30 #include "AliESDVertex.h"
31 #include "TTreeStream.h"
32 #include  "AliTPCRecoParam.h"
33 #include  "AliTPCReconstructor.h"
34 ClassImp(AliTPCtrack)
35
36 //_________________________________________________________________________
37 AliTPCtrack::AliTPCtrack(): 
38   AliKalmanTrack(),
39   fdEdx(0),
40   fSdEdx(1e10),
41   fNFoundable(0),
42   fBConstrain(kFALSE),
43   fLastPoint(-1),
44   fFirstPoint(-1),
45   fRemoval(0),
46   fTrackType(0),
47   fLab2(-1),
48   fNShared(0),
49   fReference()
50 {
51   //-------------------------------------------------
52   // default constructor
53   //-------------------------------------------------
54   for (Int_t i=0; i<kMaxRow;i++) fIndex[i]=-2;
55   for (Int_t i=0; i<4;i++) fPoints[i]=0.;
56   for (Int_t i=0; i<12;i++) fKinkPoint[i]=0.;
57   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
58   for (Int_t i=0; i<3;i++) fV0Indexes[i]=0;
59 }
60
61 //_________________________________________________________________________
62
63
64
65 AliTPCtrack::AliTPCtrack(Double_t x, Double_t alpha, const Double_t p[5],
66                          const Double_t cov[15], Int_t index) :
67   AliKalmanTrack(),
68   fdEdx(0),
69   fSdEdx(1e10),
70   fNFoundable(0),
71   fBConstrain(kFALSE),
72   fLastPoint(0),
73   fFirstPoint(0),
74   fRemoval(0),
75   fTrackType(0),
76   fLab2(0),
77   fNShared(0),
78   fReference()
79 {
80   //-----------------------------------------------------------------
81   // This is the main track constructor.
82   //-----------------------------------------------------------------
83   Double_t cnv=1./(AliTracker::GetBz()*kB2C);
84
85   Double_t pp[5]={
86     p[0],
87     p[1],
88     x*p[4] - p[2],
89     p[3],
90     p[4]*cnv
91   };
92
93   Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[5];
94   Double_t c32 = x*cov[13] - cov[8];
95   Double_t c20 = x*cov[10] - cov[3], 
96            c21 = x*cov[11] - cov[4], c42 = x*cov[14] - cov[12];
97
98   Double_t cc[15]={
99     cov[0 ],
100     cov[1 ],     cov[2 ],
101     c20,         c21,         c22,
102     cov[6 ],     cov[7 ],     c32,     cov[9 ],
103     cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv
104   };
105
106   Double_t mostProbablePt=AliExternalTrackParam::GetMostProbablePt();
107   Double_t p0=TMath::Sign(1/mostProbablePt,pp[4]);
108   Double_t w0=cc[14]/(cc[14] + p0*p0), w1=p0*p0/(cc[14] + p0*p0);
109   pp[4] = w0*p0 + w1*pp[4]; 
110   cc[10]*=w1; cc[11]*=w1; cc[12]*=w1; cc[13]*=w1; cc[14]*=w1;
111
112   Set(x,alpha,pp,cc);
113
114   SetNumberOfClusters(1);
115   
116   fIndex[0]=index;
117   for (Int_t i=1; i<kMaxRow;i++) fIndex[i]=-2;
118   for (Int_t i=0; i<4;i++) fPoints[i]=0.;
119   for (Int_t i=0; i<12;i++) fKinkPoint[i]=0.;
120   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
121   for (Int_t i=0; i<3;i++) fV0Indexes[i]=0;
122 }
123
124 //_____________________________________________________________________________
125 AliTPCtrack::AliTPCtrack(const AliESDtrack& t, TTreeSRedirector *pcstream) :
126   AliKalmanTrack(),
127   fdEdx(t.GetTPCsignal()),
128   fSdEdx(1e10),
129   fNFoundable(0),
130   fBConstrain(kFALSE),
131   fLastPoint(0),
132   fFirstPoint(0),
133   fRemoval(0),
134   fTrackType(0),
135   fLab2(0),
136   fNShared(0),
137   fReference()
138 {
139   //-----------------------------------------------------------------
140   // Conversion AliESDtrack -> AliTPCtrack.
141   //-----------------------------------------------------------------
142   const Double_t kmaxC[4]={10,10,0.1,0.1};  // cuts on the rms /fP0,fP1,fP2,fP3
143   SetNumberOfClusters(t.GetTPCclusters(fIndex));
144   SetLabel(t.GetLabel());
145   SetMass(t.GetMass());
146   for (Int_t i=0; i<4;i++) fPoints[i]=0.;
147   for (Int_t i=0; i<12;i++) fKinkPoint[i]=0.;
148   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
149   for (Int_t i=0; i<3;i++) fV0Indexes[i]=0;
150   //
151   // choose parameters to start
152   //
153   const AliTPCRecoParam * recoParam = AliTPCReconstructor::GetRecoParam();
154   Int_t reject=0;
155   AliExternalTrackParam param(t);
156
157   const AliExternalTrackParam  *tpcout=(t.GetFriendTrack())? ((AliESDfriendTrack*)(t.GetFriendTrack()))->GetTPCOut():0;
158   const AliExternalTrackParam  *tpcin = t.GetInnerParam();
159   const AliExternalTrackParam  *tpc=(tpcout)?tpcout:tpcin;
160   if (!tpc) tpc=&param;
161   Bool_t isOK=recoParam->GetUseOuterDetectors();
162   if (param.GetCovariance()[0]>kmaxC[0]*kmaxC[0]) isOK=kFALSE;
163   if (param.GetCovariance()[2]>kmaxC[1]*kmaxC[1]) isOK=kFALSE;
164   if (param.GetCovariance()[5]>kmaxC[2]*kmaxC[2]) isOK=kFALSE;
165   if (param.GetCovariance()[9]>kmaxC[3]*kmaxC[3]) isOK=kFALSE;
166   param.Rotate(tpc->GetAlpha());
167   Double_t oldX=param.GetX(),  oldY=param.GetY(),  oldZ=param.GetZ();
168   if (!isOK ){
169     param=*tpc;
170     isOK=kTRUE;
171     reject=1;
172   }
173   isOK=AliTracker::PropagateTrackToBxByBz(&param,tpc->GetX(),t.GetMass(),2.,kFALSE);
174   if (param.GetCovariance()[0]>kmaxC[0]*kmaxC[0]) isOK=kFALSE;
175   if (param.GetCovariance()[2]>kmaxC[1]*kmaxC[1]) isOK=kFALSE;
176   if (param.GetCovariance()[5]>kmaxC[2]*kmaxC[2]) isOK=kFALSE;
177   if (param.GetCovariance()[9]>kmaxC[3]*kmaxC[3]) isOK=kFALSE;
178   if (!isOK){
179     param=*tpc;
180     isOK=kTRUE;
181     reject=2;
182   }
183   if (reject>0){
184     param.ResetCovariance(4.);  // reset covariance if start from backup param
185   }
186   //
187   //
188   if (pcstream){
189     AliExternalTrackParam dummy;
190     AliExternalTrackParam *ptpc=(AliExternalTrackParam *)tpc;
191     //    if (!ptpc) ptpc=&dummy;
192     AliESDtrack *esd= (AliESDtrack *)&t;
193     (*pcstream)<<"trackP"<<
194       "reject="<<reject<<   // flag - rejection of current esd track parameters
195       "esd.="<<esd<<        // original esd track
196       "tr.="<<&param<<      // starting track parameters
197       "out.="<<ptpc<<       // backup tpc parameters
198       "\n";
199   }
200
201   Set(param.GetX(),param.GetAlpha(),param.GetParameter(),param.GetCovariance());
202
203   if ((t.GetStatus()&AliESDtrack::kTIME) == 0) return;
204   StartTimeIntegral();
205   Double_t times[10]; t.GetIntegratedTimes(times); SetIntegratedTimes(times);
206   SetIntegratedLength(t.GetIntegratedLength());
207
208   if (GetX()>oldX) {
209      Double_t dX=GetX()-oldX, dY=GetY()-oldY, dZ=GetZ()-oldZ;
210      Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ);
211      AddTimeStep(d);
212   }
213 }
214
215 //_____________________________________________________________________________
216 AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) :
217   AliKalmanTrack(t),
218   fdEdx(t.fdEdx),
219   fSdEdx(t.fSdEdx),
220   fNFoundable(t.fNFoundable),
221   fBConstrain(t.fBConstrain),
222   fLastPoint(t.fLastPoint),
223   fFirstPoint(t.fFirstPoint),
224   fRemoval(t.fRemoval),
225   fTrackType(t.fTrackType),
226   fLab2(t.fLab2),
227   fNShared(t.fNShared),
228   fReference(t.fReference)
229
230 {
231   //-----------------------------------------------------------------
232   // This is a track copy constructor.
233   //-----------------------------------------------------------------
234   Set(t.GetX(),t.GetAlpha(),t.GetParameter(),t.GetCovariance());
235
236   for (Int_t i=0; i<kMaxRow; i++) fIndex[i]=t.fIndex[i];
237   for (Int_t i=0; i<4;i++) fPoints[i]=t.fPoints[i];
238   for (Int_t i=0; i<12;i++) fKinkPoint[i]=t.fKinkPoint[i];
239   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=t.fKinkIndexes[i];
240   for (Int_t i=0; i<3;i++) fV0Indexes[i]=t.fV0Indexes[i];
241 }
242
243 AliTPCtrack& AliTPCtrack::operator=(const AliTPCtrack& o){
244   if(this!=&o){
245     AliKalmanTrack::operator=(o);
246     fdEdx = o.fdEdx;
247     for(Int_t i = 0;i<kMaxRow;++i)fIndex[i] = o.fIndex[i];
248     for(Int_t i = 0;i<4;++i)fPoints[i] = o.fPoints[i];
249     fSdEdx = o.fSdEdx;
250     fNFoundable = o.fNFoundable;
251     fBConstrain = o.fBConstrain;
252     fLastPoint  = o.fLastPoint;
253     fFirstPoint = o.fFirstPoint;
254     fTrackType  = o.fTrackType;
255     fLab2       = o.fLab2;
256     fNShared    = o.fNShared;
257     fReference  = o.fReference;
258     for(Int_t i = 0;i<12;++i) fKinkPoint[i] = o.fKinkPoint[i];
259
260     for(Int_t i = 0;i<3;++i){
261       fKinkIndexes[i] = o.fKinkIndexes[i];
262       fV0Indexes[i] = o.fV0Indexes[i];
263     }
264   }
265   return *this;
266
267 }
268
269
270 //_____________________________________________________________________________
271 Int_t AliTPCtrack::Compare(const TObject *o) const {
272   //-----------------------------------------------------------------
273   // This function compares tracks according to the their curvature
274   //-----------------------------------------------------------------
275   AliTPCtrack *t=(AliTPCtrack*)o;
276   //Double_t co=t->OneOverPt();
277   //Double_t c = OneOverPt();
278   Double_t co=t->GetSigmaY2()*t->GetSigmaZ2();
279   Double_t c =GetSigmaY2()*GetSigmaZ2();
280   if (c>co) return 1;
281   else if (c<co) return -1;
282   return 0;
283 }
284
285 Double_t AliTPCtrack::GetPredictedChi2(const AliCluster *c) const {
286   //-----------------------------------------------------------------
287   // This function calculates a predicted chi2 increment.
288   //-----------------------------------------------------------------
289   Double_t p[2]={c->GetY(), c->GetZ()};
290   Double_t cov[3]={c->GetSigmaY2(), 0., c->GetSigmaZ2()};
291   return AliExternalTrackParam::GetPredictedChi2(p,cov);
292 }
293
294 //_____________________________________________________________________________
295 Bool_t AliTPCtrack::PropagateTo(Double_t xk,Double_t rho,Double_t x0) {
296   //-----------------------------------------------------------------
297   //  This function propagates a track to a reference plane x=xk.
298   //  rho - density of the crossed matrial (g/cm^3)
299   //  x0  - radiation length of the crossed material (g/cm^2) 
300   //-----------------------------------------------------------------
301   //
302   Double_t bz=GetBz();
303   Double_t zat=0;
304   if (!GetZAt(xk, bz,zat)) return kFALSE;
305   if (TMath::Abs(zat)>250.){
306     // Don't propagate track outside of the fiducial volume - material budget not proper one
307     //
308     //AliWarning("Propagate outside of fiducial volume");
309     return kFALSE;
310   }
311
312   Double_t oldX=GetX(), oldY=GetY(), oldZ=GetZ();
313   //if (!AliExternalTrackParam::PropagateTo(xk,bz)) return kFALSE;
314   Double_t b[3]; GetBxByBz(b);
315   if (!AliExternalTrackParam::PropagateToBxByBz(xk,b)) return kFALSE;
316
317   Double_t d = TMath::Sqrt((GetX()-oldX)*(GetX()-oldX) + 
318                            (GetY()-oldY)*(GetY()-oldY) + 
319                            (GetZ()-oldZ)*(GetZ()-oldZ));
320   if (IsStartedTimeIntegral() && GetX()>oldX) AddTimeStep(d);
321
322   if (oldX < xk) d = -d;
323   if (!AliExternalTrackParam::CorrectForMeanMaterial(d*rho/x0,d*rho,GetMass(),
324       kFALSE,AliExternalTrackParam::BetheBlochGas)) return kFALSE;
325
326   return kTRUE;
327 }
328
329 //_____________________________________________________________________________
330 Bool_t 
331 AliTPCtrack::PropagateToVertex(const AliESDVertex *v,Double_t rho,Double_t x0) 
332 {
333   //-----------------------------------------------------------------
334   // This function propagates tracks to the vertex
335   // rho - density of the crossed matrial (g/cm3)
336   // x0  - radiation length of the crossed material (g/cm2) 
337   //-----------------------------------------------------------------
338   Double_t oldX=GetX(), oldY=GetY(), oldZ=GetZ();
339
340   //Double_t bz=GetBz();
341   //if (!PropagateToDCA(v,bz,kVeryBig)) return kFALSE;
342   Double_t b[3]; GetBxByBz(b);
343   if (!PropagateToDCABxByBz(v,b,kVeryBig)) return kFALSE;
344
345   Double_t d = TMath::Sqrt((GetX()-oldX)*(GetX()-oldX) + 
346                            (GetY()-oldY)*(GetY()-oldY) + 
347                            (GetZ()-oldZ)*(GetZ()-oldZ));
348
349   if (oldX < GetX()) d = -d;
350   if (!AliExternalTrackParam::CorrectForMeanMaterial(d*rho/x0,d*rho,GetMass(),
351       kFALSE,AliExternalTrackParam::BetheBlochGas)) return kFALSE;
352
353   return kTRUE;
354 }
355
356 //_____________________________________________________________________________
357 Bool_t AliTPCtrack::Update(const AliCluster *c, Double_t chisq, Int_t index) {
358   //-----------------------------------------------------------------
359   // This function associates a cluster with this track.
360   //-----------------------------------------------------------------
361   Double_t p[2]={c->GetY(), c->GetZ()};
362   Double_t cov[3]={c->GetSigmaY2(), 0., c->GetSigmaZ2()};
363
364   if (!AliExternalTrackParam::Update(p,cov)) return kFALSE;
365
366   AliTracker::FillResiduals(this,p,cov,c->GetVolumeId());
367
368   Int_t n=GetNumberOfClusters();
369   fIndex[n]=index;
370   SetNumberOfClusters(n+1);
371   SetChi2(GetChi2()+chisq);
372
373   return kTRUE;
374 }
375
376 ////////////////////////////////////////////////////////////////////////
377 // MI ADDITION
378
379 Float_t AliTPCtrack::Density(Int_t row0, Int_t row1)
380 {
381   //
382   // calculate cluster density
383   Int_t good  = 0;
384   Int_t found = 0;
385   //if (row0<fFirstPoint) row0 = fFirstPoint;
386   if (row1>fLastPoint) row1 = fLastPoint;
387
388   
389   for (Int_t i=row0;i<=row1;i++){ 
390     //    Int_t index = fClusterIndex[i];
391     Int_t index = fIndex[i];
392     if (index!=-1)  good++;
393     if (index>0)    found++;
394   }
395   Float_t density=0;
396   if (good>0) density = Float_t(found)/Float_t(good);
397   return density;
398 }
399
400
401 Float_t AliTPCtrack::Density2(Int_t row0, Int_t row1)
402 {
403   //
404   // calculate cluster density
405   Int_t good  = 0;
406   Int_t found = 0;
407   //  
408   for (Int_t i=row0;i<=row1;i++){     
409     Int_t index = fIndex[i];
410     if (index!=-1)  good++;
411     if (index>0)    found++;
412   }
413   Float_t density=0;
414   if (good>0) density = Float_t(found)/Float_t(good);
415   return density;
416 }
417
418 void  AliTPCtrack::UpdatePoints()
419 {
420   //--------------------------------------------------
421   //calculates first ,amx dens and last points
422   //--------------------------------------------------
423   Float_t density[160];
424   for (Int_t i=0;i<160;i++) density[i]=-1.;
425   fPoints[0]= 160;
426   fPoints[1] = -1;
427   //
428   Int_t ngood=0;
429   Int_t undeff=0;
430   Int_t nall =0;
431   Int_t range=20;
432   for (Int_t i=0;i<160;i++){
433     Int_t last = i-range;
434     if (nall<range) nall++;
435     if (last>=0){
436       if (fIndex[last]>0&& (fIndex[last]&0x8000)==0) ngood--;
437       if (fIndex[last]==-1) undeff--;
438     }
439     if (fIndex[i]>0&& (fIndex[i]&0x8000)==0)   ngood++;
440     if (fIndex[i]==-1) undeff++;
441     if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
442   }
443   Float_t maxdens=0;
444   Int_t indexmax =0;
445   for (Int_t i=0;i<160;i++){
446     if (density[i]<0) continue;
447     if (density[i]>maxdens){
448       maxdens=density[i];
449       indexmax=i;
450     }
451   }
452   //
453   //max dens point
454   fPoints[3] = maxdens;
455   fPoints[1] = indexmax;
456   //
457   // last point
458   for (Int_t i=indexmax;i<160;i++){
459     if (density[i]<0) continue;
460     if (density[i]<maxdens/2.) {
461       break;
462     }
463     fPoints[2]=i;
464   }
465   //
466   // first point
467   for (Int_t i=indexmax;i>0;i--){
468     if (density[i]<0) continue;
469     if (density[i]<maxdens/2.) {
470       break;
471     }
472     fPoints[0]=i;
473   }
474   //
475 }
476