]>
Commit | Line | Data |
---|---|---|
98937d93 | 1 | #include "TMatrixDSym.h" |
2 | #include "TMatrixD.h" | |
3 | #include "AliTrackFitterRieman.h" | |
24d4520d | 4 | #include "AliLog.h" |
98937d93 | 5 | |
6 | ClassImp(AliTrackFitterRieman) | |
7 | ||
8 | AliTrackFitterRieman::AliTrackFitterRieman(): | |
9 | AliTrackFitter() | |
10 | { | |
11 | // | |
12 | // default constructor | |
13 | // | |
14 | fAlpha = 0.; | |
15 | for (Int_t i=0;i<9;i++) fSumXY[i] = 0; | |
46ae650f | 16 | fSumYY = 0; |
98937d93 | 17 | for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; |
46ae650f | 18 | fSumZZ = 0; |
19 | fNUsed = 0; | |
98937d93 | 20 | fConv = kFALSE; |
21 | } | |
22 | ||
23 | ||
24 | AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner): | |
25 | AliTrackFitter(array,owner) | |
26 | { | |
27 | // | |
28 | // Constructor | |
29 | // | |
30 | fAlpha = 0.; | |
31 | for (Int_t i=0;i<9;i++) fSumXY[i] = 0; | |
46ae650f | 32 | fSumYY = 0; |
98937d93 | 33 | for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; |
46ae650f | 34 | fSumZZ = 0; |
35 | fNUsed = 0; | |
98937d93 | 36 | fConv = kFALSE; |
37 | } | |
38 | ||
39 | AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman): | |
40 | AliTrackFitter(rieman) | |
41 | { | |
42 | // | |
43 | // copy constructor | |
44 | // | |
45 | fAlpha = rieman.fAlpha; | |
46 | for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i]; | |
46ae650f | 47 | fSumYY = rieman.fSumYY; |
98937d93 | 48 | for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i]; |
46ae650f | 49 | fSumZZ = rieman.fSumZZ; |
50 | fNUsed = rieman.fNUsed; | |
98937d93 | 51 | fConv = rieman.fConv; |
52 | } | |
53 | ||
54 | //_____________________________________________________________________________ | |
55 | AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman) | |
56 | { | |
57 | // assignment operator | |
58 | // | |
59 | if(this==&rieman) return *this; | |
60 | ((AliTrackFitter *)this)->operator=(rieman); | |
61 | ||
62 | fAlpha = rieman.fAlpha; | |
63 | for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i]; | |
46ae650f | 64 | fSumYY = rieman.fSumYY; |
98937d93 | 65 | for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i]; |
46ae650f | 66 | fSumZZ = rieman.fSumZZ; |
67 | fNUsed = rieman.fNUsed; | |
98937d93 | 68 | fConv = rieman.fConv; |
69 | ||
70 | return *this; | |
71 | } | |
72 | ||
73 | AliTrackFitterRieman::~AliTrackFitterRieman() | |
74 | { | |
75 | // destructor | |
76 | // | |
77 | } | |
78 | ||
79 | void AliTrackFitterRieman::Reset() | |
80 | { | |
81 | // Reset the track parameters and | |
82 | // rieman sums | |
83 | AliTrackFitter::Reset(); | |
84 | fAlpha = 0.; | |
85 | for (Int_t i=0;i<9;i++) fSumXY[i] = 0; | |
46ae650f | 86 | fSumYY = 0; |
98937d93 | 87 | for (Int_t i=0;i<9;i++) fSumXZ[i] = 0; |
46ae650f | 88 | fSumZZ = 0; |
89 | fNUsed = 0; | |
98937d93 | 90 | fConv =kFALSE; |
91 | } | |
92 | ||
cc345ce3 | 93 | Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit, |
98937d93 | 94 | AliAlignObj::ELayerID layerRangeMin, |
95 | AliAlignObj::ELayerID layerRangeMax) | |
96 | { | |
97 | // Fit the track points. The method takes as an input | |
cc345ce3 | 98 | // the set of id's (volids) of the volumes in which |
99 | // one wants to calculate the residuals. | |
100 | // The following parameters are used to define the | |
98937d93 | 101 | // range of volumes to be used in the fitting |
102 | // As a result two AliTrackPointArray's obects are filled. | |
103 | // The first one contains the space points with | |
cc345ce3 | 104 | // volume id's from volids list. The second array of points represents |
105 | // the track extrapolations corresponding to the space points | |
98937d93 | 106 | // in the first array. The two arrays can be used to find |
cc345ce3 | 107 | // the residuals in the volids and consequently construct a |
98937d93 | 108 | // chi2 function to be minimized during the alignment |
109 | // procedures. For the moment the track extrapolation is taken | |
cc345ce3 | 110 | // at the space-point reference plane. The reference plane is |
111 | // found using the covariance matrix of the point | |
112 | // (assuming sigma(x)=0 at the reference coordinate system. | |
98937d93 | 113 | |
46ae650f | 114 | Reset(); |
98937d93 | 115 | |
116 | Int_t npoints = fPoints->GetNPoints(); | |
117 | if (npoints < 3) return kFALSE; | |
118 | ||
cc345ce3 | 119 | Bool_t isAlphaCalc = kFALSE; |
46ae650f | 120 | AliTrackPoint p,plocal; |
cc345ce3 | 121 | // fPoints->GetPoint(p,0); |
122 | // fAlpha = TMath::ATan2(p.GetY(),p.GetX()); | |
98937d93 | 123 | |
124 | Int_t npVolId = 0; | |
46ae650f | 125 | fNUsed = 0; |
98937d93 | 126 | Int_t *pindex = new Int_t[npoints]; |
cc345ce3 | 127 | fX = new Float_t[npoints]; |
128 | fY = new Float_t[npoints]; | |
129 | fZ = new Float_t[npoints]; | |
130 | fSy = new Float_t[npoints]; | |
131 | fSz = new Float_t[npoints]; | |
98937d93 | 132 | for (Int_t ipoint = 0; ipoint < npoints; ipoint++) |
133 | { | |
134 | fPoints->GetPoint(p,ipoint); | |
135 | UShort_t iVolId = p.GetVolumeID(); | |
cc345ce3 | 136 | if (FindVolId(volIds,iVolId)) { |
98937d93 | 137 | pindex[npVolId] = ipoint; |
138 | npVolId++; | |
139 | } | |
cc345ce3 | 140 | if (volIdsFit != 0x0) { |
141 | if (!FindVolId(volIdsFit,iVolId)) continue; | |
46ae650f | 142 | } |
143 | else { | |
144 | if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) || | |
145 | iVolId > AliAlignObj::LayerToVolUID(layerRangeMax, | |
c041444f | 146 | AliAlignObj::LayerSize(layerRangeMax))) continue; |
46ae650f | 147 | } |
cc345ce3 | 148 | if (!isAlphaCalc) { |
149 | fAlpha = p.GetAngle(); | |
150 | isAlphaCalc = kTRUE; | |
151 | } | |
46ae650f | 152 | plocal = p.Rotate(fAlpha); |
153 | AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(), | |
154 | TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5])); | |
155 | fNUsed++; | |
98937d93 | 156 | } |
157 | ||
24d4520d | 158 | if (npVolId == 0 || fNUsed < 3) { |
98937d93 | 159 | delete [] pindex; |
cc345ce3 | 160 | delete [] fX; |
161 | delete [] fY; | |
162 | delete [] fZ; | |
163 | delete [] fSy; | |
164 | delete [] fSz; | |
98937d93 | 165 | return kFALSE; |
166 | } | |
167 | ||
168 | Update(); | |
169 | ||
cc345ce3 | 170 | delete [] fX; |
171 | delete [] fY; | |
172 | delete [] fZ; | |
173 | delete [] fSy; | |
174 | delete [] fSz; | |
175 | ||
98937d93 | 176 | if (!fConv) { |
177 | delete [] pindex; | |
178 | return kFALSE; | |
179 | } | |
180 | ||
181 | if ((fParams[0] == 0) || | |
182 | ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) { | |
183 | delete [] pindex; | |
184 | return kFALSE; | |
185 | } | |
186 | ||
cc345ce3 | 187 | |
188 | if (fNUsed < fMinNPoints) { | |
189 | delete [] pindex; | |
190 | return kFALSE; | |
191 | } | |
192 | ||
46ae650f | 193 | fPVolId = new AliTrackPointArray(npVolId); |
194 | fPTrack = new AliTrackPointArray(npVolId); | |
98937d93 | 195 | AliTrackPoint p2; |
196 | for (Int_t ipoint = 0; ipoint < npVolId; ipoint++) | |
197 | { | |
198 | Int_t index = pindex[ipoint]; | |
199 | fPoints->GetPoint(p,index); | |
200 | if (GetPCA(p,p2)) { | |
cc345ce3 | 201 | Float_t xyz[3],xyz2[3]; |
202 | p.GetXYZ(xyz); p2.GetXYZ(xyz2); | |
203 | // printf("residuals %f %d %d %f %f %f %f %f %f\n",fChi2,fNUsed,fConv,xyz[0],xyz[1],xyz[2],xyz2[0]-xyz[0],xyz2[1]-xyz[1],xyz2[2]-xyz[2]); | |
46ae650f | 204 | fPVolId->AddPoint(ipoint,&p); |
205 | fPTrack->AddPoint(ipoint,&p2); | |
98937d93 | 206 | } |
207 | } | |
208 | ||
209 | delete [] pindex; | |
210 | ||
46ae650f | 211 | // debug info |
212 | // Float_t chi2 = 0, chi22 = 0; | |
213 | // for (Int_t ipoint = 0; ipoint < npoints; ipoint++) | |
214 | // { | |
215 | // fPoints->GetPoint(p,ipoint); | |
216 | // UShort_t iVolId = p.GetVolumeID(); | |
217 | // if (volIdFit != 0) { | |
218 | // if (iVolId != volIdFit) continue; | |
219 | // } | |
220 | // else { | |
221 | // if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) || | |
c041444f | 222 | // iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,AliAlignObj::LayerSize(layerRangeMax))) continue; |
46ae650f | 223 | // } |
224 | // plocal = p.Rotate(fAlpha); | |
225 | // Float_t delta = (fParams[0]*(plocal.GetX()*plocal.GetX()+plocal.GetY()*plocal.GetY())+ | |
226 | // 2.*plocal.GetX()*fParams[1]+ | |
227 | // fParams[2]- | |
228 | // 2.*plocal.GetY())/ | |
229 | // (2.*TMath::Sqrt(plocal.GetCov()[3])); | |
230 | // // Float_t delta2 = (fParams[3]+ | |
231 | // // plocal.GetX()*fParams[4]+ | |
232 | // // plocal.GetX()*plocal.GetX()*fParams[5]- | |
233 | // // plocal.GetZ())/ | |
234 | // // (TMath::Sqrt(plocal.GetCov()[5])); | |
235 | // Double_t r = TMath::Sqrt(plocal.GetX()*plocal.GetX()+plocal.GetY()*plocal.GetY()); | |
236 | // Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
237 | // Float_t delta2 = (fParams[3]+ | |
238 | // r*fParams[4]+r*r*r*fParams[4]*Rm1*Rm1/24- | |
239 | // plocal.GetZ())/ | |
240 | // (TMath::Sqrt(plocal.GetCov()[5])); | |
241 | // chi2 += delta*delta; | |
242 | // chi22 += delta2*delta2; | |
243 | // // printf("pulls %d %d %f %f\n",ipoint,iVolId,delta,delta2); | |
244 | ||
245 | // } | |
246 | // printf("My chi2 = %f + %f = %f\n",chi2,chi22,chi2+chi22); | |
247 | ||
98937d93 | 248 | return kTRUE; |
249 | } | |
250 | ||
251 | void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz) | |
252 | { | |
253 | // | |
254 | // Rieman update | |
255 | // | |
256 | //------------------------------------------------------ | |
257 | // XY direction | |
258 | // | |
259 | // (x-x0)^2+(y-y0)^2-R^2=0 ===> | |
260 | // | |
261 | // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> | |
262 | // | |
46ae650f | 263 | // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2 |
98937d93 | 264 | // |
265 | // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation | |
266 | // | |
267 | // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 | |
268 | // | |
269 | // final linear equation : a + u*b +t*c - v =0; | |
270 | // | |
271 | // Minimization : | |
272 | // | |
273 | // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min; | |
274 | // | |
275 | // where sigmai is the error of maesurement (a + ui*b +ti*c - vi) | |
276 | // | |
277 | // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t | |
278 | // | |
cc345ce3 | 279 | fX[fNUsed] = x; fY[fNUsed]=y; fZ[fNUsed]=z; fSy[fNUsed]=sy; fSz[fNUsed]=sz; |
98937d93 | 280 | // |
281 | // XY part | |
282 | // | |
283 | Double_t t = x*x+y*y; | |
284 | if (t<2) return; | |
285 | t = 1./t; | |
286 | Double_t u = 2.*x*t; | |
287 | Double_t v = 2.*y*t; | |
288 | Double_t error = 2.*sy*t; | |
289 | error *=error; | |
290 | Double_t weight = 1./error; | |
291 | fSumXY[0] +=weight; | |
292 | fSumXY[1] +=u*weight; fSumXY[2]+=v*weight; fSumXY[3]+=t*weight; | |
293 | fSumXY[4] +=u*u*weight; fSumXY[5]+=t*t*weight; | |
294 | fSumXY[6] +=u*v*weight; fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight; | |
46ae650f | 295 | fSumYY += v*v*weight; |
98937d93 | 296 | // |
297 | // XZ part | |
298 | // | |
cc345ce3 | 299 | if (1) { |
46ae650f | 300 | weight = 1./(sz*sz); |
cc345ce3 | 301 | // fSumXZ[0] +=weight; |
302 | // fSumXZ[1] +=weight*x; fSumXZ[2] +=weight*x*x; fSumXZ[3] +=weight*x*x*x; fSumXZ[4] += weight*x*x*x*x; | |
303 | // fSumXZ[5] +=weight*z; fSumXZ[6] +=weight*x*z; fSumXZ[7] +=weight*x*x*z; | |
46ae650f | 304 | fSumZZ += z*z*weight; |
305 | } | |
306 | else { | |
307 | weight = 1./(sz*sz); | |
308 | fSumXZ[0] +=weight; | |
309 | Double_t r = TMath::Sqrt(x*x+y*y); | |
310 | fSumXZ[1] +=weight*r; fSumXZ[2] +=weight*r*r; fSumXZ[3] +=weight*z; fSumXZ[4] += weight*r*z; | |
311 | // Now the auxulary sums | |
312 | fSumXZ[5] +=weight*r*r*r/24; fSumXZ[6] +=weight*r*r*r*r/12; fSumXZ[7] +=weight*r*r*r*z/24; | |
313 | fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24); | |
314 | fSumZZ += z*z*weight; | |
315 | } | |
98937d93 | 316 | } |
317 | ||
318 | void AliTrackFitterRieman::Update(){ | |
319 | // | |
320 | // Rieman update | |
321 | // | |
322 | // | |
323 | for (Int_t i=0;i<6;i++)fParams[i]=0; | |
46ae650f | 324 | fChi2 = 0; |
325 | fNdf = 0; | |
98937d93 | 326 | Int_t conv=0; |
327 | // | |
328 | // XY part | |
329 | // | |
330 | TMatrixDSym smatrix(3); | |
331 | TMatrixD sums(1,3); | |
332 | // | |
333 | // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2; | |
334 | // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut; | |
335 | // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt; | |
336 | ||
337 | smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5]; | |
338 | smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7]; | |
339 | sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8]; | |
340 | smatrix.Invert(); | |
341 | if (smatrix.IsValid()){ | |
342 | for (Int_t i=0;i<3;i++) | |
343 | for (Int_t j=0;j<=i;j++){ | |
344 | (*fCov)(i,j)=smatrix(i,j); | |
345 | } | |
346 | TMatrixD res = sums*smatrix; | |
347 | fParams[0] = res(0,0); | |
348 | fParams[1] = res(0,1); | |
349 | fParams[2] = res(0,2); | |
46ae650f | 350 | TMatrixD tmp = res*sums.T(); |
351 | fChi2 += fSumYY - tmp(0,0); | |
352 | fNdf += fNUsed - 3; | |
98937d93 | 353 | conv++; |
354 | } | |
355 | // | |
356 | // XZ part | |
357 | // | |
cc345ce3 | 358 | if (1) { |
359 | Double_t x0 = -fParams[1]/fParams[0]; | |
360 | Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
361 | ||
362 | for (Int_t i=0;i<fNUsed;i++){ | |
363 | Double_t phi = TMath::ASin((fX[i]-x0)*Rm1); | |
364 | Double_t phi0 = TMath::ASin((0.-x0)*Rm1); | |
365 | Double_t weight = 1/fSz[i]; | |
366 | weight *=weight; | |
367 | Double_t dphi = (phi-phi0)/Rm1; | |
368 | fSumXZ[0] +=weight; | |
369 | fSumXZ[1] +=weight*dphi; | |
370 | fSumXZ[2] +=weight*dphi*dphi; | |
371 | fSumXZ[3] +=weight*fZ[i]; | |
372 | fSumXZ[4] +=weight*dphi*fZ[i]; | |
373 | } | |
374 | ||
375 | TMatrixDSym smatrixz(2); | |
376 | smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2]; | |
46ae650f | 377 | smatrixz.Invert(); |
cc345ce3 | 378 | TMatrixD sumsxz(1,2); |
46ae650f | 379 | if (smatrixz.IsValid()){ |
cc345ce3 | 380 | sumsxz(0,0) = fSumXZ[3]; |
381 | sumsxz(0,1) = fSumXZ[4]; | |
46ae650f | 382 | TMatrixD res = sumsxz*smatrixz; |
383 | fParams[3] = res(0,0); | |
384 | fParams[4] = res(0,1); | |
cc345ce3 | 385 | fParams[5] = 0; |
386 | for (Int_t i=0;i<2;i++) | |
46ae650f | 387 | for (Int_t j=0;j<=i;j++){ |
cc345ce3 | 388 | (*fCov)(i+3,j+3)=smatrixz(i,j); |
46ae650f | 389 | } |
390 | TMatrixD tmp = res*sumsxz.T(); | |
391 | fChi2 += fSumZZ - tmp(0,0); | |
cc345ce3 | 392 | fNdf += fNUsed - 2; |
46ae650f | 393 | conv++; |
394 | } | |
395 | } | |
396 | else { | |
397 | Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
398 | fSumXZ[1] += fSumXZ[5]*Rm1*Rm1; | |
399 | fSumXZ[2] += fSumXZ[6]*Rm1*Rm1 + fSumXZ[8]*Rm1*Rm1*Rm1*Rm1; | |
400 | fSumXZ[4] += fSumXZ[7]*Rm1*Rm1; | |
401 | ||
402 | TMatrixDSym smatrixz(2); | |
403 | smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2]; | |
404 | smatrixz.Invert(); | |
405 | TMatrixD sumsxz(1,2); | |
406 | if (smatrixz.IsValid()){ | |
407 | sumsxz(0,0) = fSumXZ[3]; | |
408 | sumsxz(0,1) = fSumXZ[4]; | |
409 | TMatrixD res = sumsxz*smatrixz; | |
410 | fParams[3] = res(0,0); | |
411 | fParams[4] = res(0,1); | |
412 | fParams[5] = 0; | |
413 | for (Int_t i=0;i<2;i++) | |
414 | for (Int_t j=0;j<=i;j++){ | |
415 | (*fCov)(i+3,j+3)=smatrixz(i,j); | |
416 | } | |
417 | TMatrixD tmp = res*sumsxz.T(); | |
418 | fChi2 += fSumZZ - tmp(0,0); | |
419 | fNdf += fNUsed - 2; | |
420 | conv++; | |
98937d93 | 421 | } |
98937d93 | 422 | } |
423 | ||
424 | // (x-x0)^2+(y-y0)^2-R^2=0 ===> | |
425 | // | |
426 | // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==> | |
427 | // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2 | |
428 | // | |
429 | // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation | |
430 | // | |
431 | // next substition a = 1/y0 b = -x0/y0 c = -D0/y0 | |
432 | // final linear equation : a + u*b +t*c - v =0; | |
433 | // | |
434 | // fParam[0] = 1/y0 | |
435 | // fParam[1] = -x0/y0 | |
436 | // fParam[2] = - (R^2 - x0^2 - y0^2)/y0 | |
437 | if (conv>1) fConv =kTRUE; | |
438 | else | |
439 | fConv=kFALSE; | |
440 | } | |
441 | ||
442 | Double_t AliTrackFitterRieman::GetYat(Double_t x){ | |
443 | if (!fConv) return 0.; | |
444 | Double_t res2 = (x*fParams[0]+fParams[1]); | |
445 | res2*=res2; | |
446 | res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2; | |
447 | if (res2<0) return 0; | |
448 | res2 = TMath::Sqrt(res2); | |
449 | res2 = (1-res2)/fParams[0]; | |
450 | return res2; | |
451 | } | |
452 | ||
cc345ce3 | 453 | Double_t AliTrackFitterRieman::GetDYat(Double_t x) const { |
98937d93 | 454 | if (!fConv) return 0.; |
455 | Double_t x0 = -fParams[1]/fParams[0]; | |
456 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0; | |
457 | Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
458 | if ( 1./(Rm1*Rm1)-(x-x0)*(x-x0)<=0) return 0; | |
459 | Double_t res = (x-x0)/TMath::Sqrt(1./(Rm1*Rm1)-(x-x0)*(x-x0)); | |
460 | if (fParams[0]<0) res*=-1.; | |
461 | return res; | |
462 | } | |
463 | ||
464 | ||
465 | ||
cc345ce3 | 466 | Double_t AliTrackFitterRieman::GetZat(Double_t x) const { |
98937d93 | 467 | if (!fConv) return 0.; |
468 | Double_t x0 = -fParams[1]/fParams[0]; | |
469 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; | |
470 | Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
471 | Double_t phi = TMath::ASin((x-x0)*Rm1); | |
472 | Double_t phi0 = TMath::ASin((0.-x0)*Rm1); | |
473 | Double_t dphi = (phi-phi0); | |
474 | Double_t res = fParams[3]+fParams[4]*dphi/Rm1; | |
475 | return res; | |
476 | } | |
477 | ||
cc345ce3 | 478 | Double_t AliTrackFitterRieman::GetDZat(Double_t x) const { |
98937d93 | 479 | if (!fConv) return 0.; |
480 | Double_t x0 = -fParams[1]/fParams[0]; | |
481 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; | |
482 | Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
483 | Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*Rm1)); | |
484 | return res; | |
485 | } | |
486 | ||
487 | ||
488 | Double_t AliTrackFitterRieman::GetC(){ | |
489 | return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
490 | } | |
491 | ||
492 | Bool_t AliTrackFitterRieman::GetXYZat(Double_t r, Float_t *xyz){ | |
493 | if (!fConv) return kFALSE; | |
494 | Double_t res2 = (r*fParams[0]+fParams[1]); | |
495 | res2*=res2; | |
496 | res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2; | |
497 | if (res2<0) return kFALSE; | |
498 | res2 = TMath::Sqrt(res2); | |
499 | res2 = (1-res2)/fParams[0]; | |
500 | ||
501 | if (!fConv) return kFALSE; | |
502 | Double_t x0 = -fParams[1]/fParams[0]; | |
503 | if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0; | |
504 | Double_t Rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); | |
505 | Double_t phi = TMath::ASin((r-x0)*Rm1); | |
506 | Double_t phi0 = TMath::ASin((0.-x0)*Rm1); | |
507 | Double_t dphi = (phi-phi0); | |
508 | Double_t res = fParams[3]+fParams[4]*dphi/Rm1; | |
509 | ||
510 | Double_t sin = TMath::Sin(fAlpha); | |
511 | Double_t cos = TMath::Cos(fAlpha); | |
512 | xyz[0] = r*cos - res2*sin; | |
513 | xyz[1] = res2*cos + r*sin; | |
514 | xyz[2] = res; | |
515 | ||
516 | return kTRUE; | |
517 | } | |
518 | ||
519 | Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const | |
520 | { | |
521 | // Get the closest to a given spacepoint track trajectory point | |
522 | // Look for details in the description of the Fit() method | |
523 | ||
524 | if (!fConv) return kFALSE; | |
525 | ||
526 | // First X and Y coordinates | |
527 | Double_t sin = TMath::Sin(fAlpha); | |
528 | Double_t cos = TMath::Cos(fAlpha); | |
529 | // fParam[0] = 1/y0 | |
530 | // fParam[1] = -x0/y0 | |
531 | // fParam[2] = - (R^2 - x0^2 - y0^2)/y0 | |
532 | if (fParams[0] == 0) return kFALSE; | |
cc345ce3 | 533 | // Track parameters in the global coordinate system |
98937d93 | 534 | Double_t x0 = -fParams[1]/fParams[0]*cos - 1./fParams[0]*sin; |
535 | Double_t y0 = 1./fParams[0]*cos - fParams[1]/fParams[0]*sin; | |
536 | if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE; | |
537 | Double_t R = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/ | |
538 | fParams[0]; | |
539 | ||
cc345ce3 | 540 | // Define space-point refence plane |
541 | Double_t alphap = p.GetAngle(); | |
542 | Double_t sinp = TMath::Sin(alphap); | |
543 | Double_t cosp = TMath::Cos(alphap); | |
544 | Double_t x = p.GetX()*cosp + p.GetY()*sinp; | |
545 | Double_t y = p.GetY()*cosp - p.GetX()*sinp; | |
546 | Double_t x0p= x0*cosp + y0*sinp; | |
547 | Double_t y0p= y0*cosp - x0*sinp; | |
24d4520d | 548 | if ((R*R - (x-x0p)*(x-x0p))<0) { |
549 | AliWarning(Form("Track extrapolation failed ! (Track radius = %f, track circle x = %f, space-point x = %f, reference plane angle = %f\n",R,x0p,x,alphap)); | |
550 | return kFALSE; | |
551 | } | |
cc345ce3 | 552 | Double_t temp = TMath::Sqrt(R*R - (x-x0p)*(x-x0p)); |
553 | Double_t y1 = y0p + temp; | |
554 | Double_t y2 = y0p - temp; | |
555 | Double_t yprime = y1; | |
556 | if(TMath::Abs(y2-y) < TMath::Abs(y1-y)) yprime = y2; | |
557 | ||
558 | // Back to the global coordinate system | |
559 | Double_t xsecond = x*cosp - yprime*sinp; | |
560 | Double_t ysecond = yprime*cosp + x*sinp; | |
561 | ||
562 | // Now Z coordinate and track angles | |
563 | Double_t x2 = xsecond*cos + ysecond*sin; | |
564 | Double_t zsecond = GetZat(x2); | |
565 | Double_t dydx = GetDYat(x2); | |
566 | Double_t dzdx = GetDZat(x2); | |
567 | ||
568 | // Fill the cov matrix of the track extrapolation point | |
569 | Double_t cov[6] = {0,0,0,0,0,0}; | |
570 | Double_t sigmax = 100*100.; | |
571 | cov[0] = sigmax; cov[1] = sigmax*dydx; cov[2] = sigmax*dzdx; | |
572 | cov[3] = sigmax*dydx*dydx; cov[4] = sigmax*dydx*dzdx; | |
573 | cov[5] = sigmax*dzdx*dzdx; | |
574 | ||
575 | Float_t newcov[6]; | |
576 | newcov[0] = cov[0]*cos*cos- | |
577 | 2*cov[1]*sin*cos+ | |
578 | cov[3]*sin*sin; | |
579 | newcov[1] = cov[1]*(cos*cos-sin*sin)- | |
580 | (cov[3]-cov[0])*sin*cos; | |
581 | newcov[2] = cov[2]*cos- | |
582 | cov[4]*sin; | |
583 | newcov[3] = cov[0]*sin*sin+ | |
584 | 2*cov[1]*sin*cos+ | |
585 | cov[3]*cos*cos; | |
586 | newcov[4] = cov[4]*cos+ | |
587 | cov[2]*sin; | |
588 | newcov[5] = cov[5]; | |
589 | ||
590 | p2.SetXYZ(xsecond,ysecond,zsecond,newcov); | |
98937d93 | 591 | |
592 | return kTRUE; | |
593 | } |