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