Eff C++
[u/mrichter/AliRoot.git] / TPC / AliTPCpolyTrack.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 // this class is used for a simple
17 // track fitting using the polynom
18 // bending in both directions are allowed
19 //-----------------------------------------------------------------
20
21
22 #include "AliTPCpolyTrack.h"
23 #include "TMath.h"
24
25 ClassImp(AliTPCpolyTrack)
26
27
28 AliTPCpolyTrack::AliTPCpolyTrack()
29   :TObject(),
30    fA(0), 
31    fB(0), 
32    fC(0), 
33    fD(0), 
34    fE(0), 
35    fF(0), 
36    fMaxX(0), 
37    fMinX(0), 
38    fSumW(0),    
39    fSumX(0),    
40    fSumX2(0),   
41    fSumX3(0),    
42    fSumX4(0),   
43    fSumY(0),    
44    fSumYX(0),   
45    fSumYX2(0),  
46    fSumZ(0),     
47    fSumZX(0),   
48    fSumZX2(0), 
49    fNPoints(0)
50 {
51   Reset();
52 }
53
54 void   AliTPCpolyTrack::Reset()
55 {
56   //
57   // reset track
58   fSumX = fSumX2= fSumX3=fSumX4 = fSumY=fSumYX=fSumYX2=fSumZ=fSumZX=fSumZX2=fSumW =0;
59   fNPoints = 0;
60 }
61
62 void AliTPCpolyTrack::AddPoint(Double_t x, Double_t y, Double_t z,Double_t sy, Double_t sz)
63 {
64   //
65   //
66   if (fNPoints==0){
67     fMaxX = x;
68     fMinX = x;
69   }else{
70     if (x>fMaxX) fMaxX=x;
71     if (x<fMinX) fMinX=x;
72   }
73
74   Double_t x2 = x*x; 
75   Double_t w = 2./(sy+sz);
76   fSumW += w;
77   //
78   fSumX       += x*w;
79   fSumX2      += x2*w;
80   fSumX3      += x2*x*w;
81   fSumX4      += x2*x2*w;
82   //
83   fSumY       +=y*w;
84   fSumYX      +=y*x*w;
85   fSumYX2     +=y*x2*w;
86   //
87   fSumZ       +=z*w;
88   fSumZX      +=z*x*w;
89   fSumZX2     +=z*x2*w;
90   //
91   fX[fNPoints] = x;
92   fY[fNPoints] = y;
93   fZ[fNPoints] = z;  
94   fSY[fNPoints] = sy;
95   fSZ[fNPoints] = sz;  
96
97   fNPoints++;
98   
99 }
100
101 void  AliTPCpolyTrack::UpdateParameters()
102 {
103   //
104   //
105   //Update fit parameters
106   if (fNPoints>4){
107     Fit2(fSumY,fSumYX,fSumYX2,fSumX,fSumX2,fSumX3,fSumX4,fSumW,fA,fB,fC);
108     //    Fit2(fSumZ,fSumZX,fSumZX2,fSumX,fSumX2,fSumX3,fSumX4,fNPoints,fD,fE,fF);
109     Fit1(fSumZ,fSumZX,fSumX,fSumX2,fSumW,fD,fE,fF);
110   }
111   else
112     {
113       Fit1(fSumY,fSumYX,fSumX,fSumX2,fSumW,fA,fB,fC);
114       Fit1(fSumZ,fSumZX,fSumX,fSumX2,fSumW,fD,fE,fF);
115     }
116 }
117
118 void  AliTPCpolyTrack::UpdateParameters(Int_t ny, Int_t nz)
119 {
120   //
121   //
122   //Update fit parameters
123   if (ny>1){
124     Fit2(fSumY,fSumYX,fSumYX2,fSumX,fSumX2,fSumX3,fSumX4,fSumW,fA,fB,fC);
125   }
126   else{ 
127     Fit1(fSumY,fSumYX,fSumX,fSumX2,fSumW,fA,fB,fC);
128   }
129   if (nz>1){
130     Fit2(fSumZ,fSumZX,fSumZX2,fSumX,fSumX2,fSumX3,fSumX4,fNPoints,fD,fE,fF);
131   }
132   else{
133     Fit1(fSumZ,fSumZX,fSumX,fSumX2,fSumW,fD,fE,fF);
134   }
135 }
136
137
138
139 void  AliTPCpolyTrack::Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
140             Double_t fSumX,  Double_t fSumX2, Double_t fSumX3, 
141             Double_t fSumX4, Double_t fSumW,
142             Double_t &a, Double_t &b, Double_t &c)
143 {
144   //fit of second order
145   Double_t det = 
146     fSumW* (fSumX2*fSumX4-fSumX3*fSumX3) -
147     fSumX*      (fSumX*fSumX4-fSumX3*fSumX2)+
148     fSumX2*     (fSumX*fSumX3-fSumX2*fSumX2);
149     
150   if (TMath::Abs(det)> 0.000000000000001) {    
151     a = 
152       (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
153        fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
154        fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det; 
155     b=
156       (fSumW*(fSumYX*fSumX4-fSumX3*fSumYX2)-
157       fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
158       fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
159     c=
160       (fSumW*(fSumX2*fSumYX2-fSumYX*fSumX3)-
161        fSumX*(fSumX*fSumYX2-fSumYX*fSumX2)+
162        fSumY*(fSumX*fSumX3-fSumX2*fSumX2))/det;  
163   }
164 }
165
166 void  AliTPCpolyTrack::Fit1(Double_t fSumY, Double_t fSumYX, 
167               Double_t fSumX,  Double_t fSumX2, 
168               Double_t fSumW, Double_t &a, Double_t &b, Double_t &c)
169 {
170   //
171   //
172   //
173   Double_t det = fSumW*fSumX2-fSumX*fSumX;
174   if (TMath::Abs(det)> 0.000000000000001) { 
175     b = (fSumW*fSumYX-fSumX*fSumY)/det;
176     a = (fSumX2*fSumY-fSumX*fSumYX)/det;
177     c = 0;
178   }else{
179     a =fSumYX/fSumX;
180     b =0;
181     c =0;
182   }
183
184 }
185
186 void AliTPCpolyTrack::Refit(AliTPCpolyTrack &track, Double_t deltay, Double_t deltaz)
187 {
188   //
189   // refit with cut on distortion
190   //
191   track.Reset();
192   //first refit to temporary
193   AliTPCpolyTrack track0;
194   track0.Reset();
195   for (Int_t i=0;i<fNPoints;i++){
196     Double_t y,z;
197     GetFitPoint(fX[i],y,z);
198     if ( (TMath::Abs(y-fY[i])<deltay)&&(TMath::Abs(z-fZ[i])<deltaz)){
199       track0.AddPoint(fX[i],y,z);
200     }
201   }
202   if (track0.GetN()>2) 
203     track0.UpdateParameters();
204   else 
205     return;
206   //
207   for (Int_t i=0;i<fNPoints;i++){
208     Double_t y,z;
209     track0.GetFitPoint(fX[i],y,z);
210     if ( (TMath::Abs(y-fY[i])<deltay)&&(TMath::Abs(z-fZ[i])<deltaz)){
211       track.AddPoint(fX[i],y,z);
212     }
213   }
214   if (track.GetN()>2) 
215     track.UpdateParameters();
216
217 }
218
219 void AliTPCpolyTrack::Refit(AliTPCpolyTrack &track, Double_t deltay, Double_t deltaz, Int_t nfirst, Int_t ny, Int_t nz)
220 {
221   //
222   // refit with cut on distortion
223   //
224   track.Reset();
225   //first refit to temporary
226   AliTPCpolyTrack track0;
227   track0.Reset();
228   for (Int_t i=0;i<fNPoints;i++){
229     Double_t y,z;
230     GetFitPoint(fX[i],y,z);
231     if ( (TMath::Abs(y-fY[i])<deltay)&&(TMath::Abs(z-fZ[i])<deltaz)){
232       track0.AddPoint(fX[i],y,z);
233     }    
234   }
235   if (track0.GetN()>2){ 
236     if (track0.GetN()>nfirst)
237       track0.UpdateParameters(ny,nz);
238     else 
239       track0.UpdateParameters(1,1);
240   }
241   else 
242     return;
243   //
244   for (Int_t i=0;i<fNPoints;i++){
245     Double_t y,z;
246     track0.GetFitPoint(fX[i],y,z);
247     if ( (TMath::Abs(y-fY[i])<deltay)&&(TMath::Abs(z-fZ[i])<deltaz)){
248       track.AddPoint(fX[i],y,z);
249     }
250   }
251   if (track.GetN()>2) 
252     track.UpdateParameters(ny,nz);
253
254 }
255