]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/STEER/AliHelix.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Compatibility with ROOT trunk
[u/mrichter/AliRoot.git] / STEER / STEER / AliHelix.cxx
CommitLineData
7f572c00 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 AliHelix class
20// Origin: Marian Ivanov, CERN, marian.ivanov@cern.ch
21//-------------------------------------------------------------------------
22
23
24#include "AliHelix.h"
25#include "AliKalmanTrack.h"
6c94f330 26#include "AliTracker.h"
7f572c00 27#include "TMath.h"
28ClassImp(AliHelix)
29
30
31//_______________________________________________________________________
32AliHelix::AliHelix()
33{
34 //
35 // Default constructor
36 //
37 for (Int_t i =0;i<9;i++) fHelix[i]=0;
38}
39
40//_______________________________________________________________________
176aff27 41AliHelix::AliHelix(const AliHelix &t):TObject(t){
7f572c00 42 //
43 //
44 for (Int_t i=0;i<9;i++)
45 fHelix[i]=t.fHelix[i];
46}
47
48AliHelix::AliHelix(const AliKalmanTrack &t)
49{
50 //
51 //
52 Double_t alpha,x,cs,sn;
53 t.GetExternalParameters(x,fHelix);
54 alpha=t.GetAlpha();
55 //
56 //circle parameters
022ba35d 57 //PH Sometimes fP4 and fHelix[4] are very big and the calculation
58 //PH of the Sqrt cannot be done. To be investigated...
190fc854 59 fHelix[4]=fHelix[4]/(-1000/0.299792458/AliTracker::GetBz()); // C
7f572c00 60 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
61
62 Double_t xc, yc, rc;
63 rc = 1/fHelix[4];
64 xc = x-fHelix[2]*rc;
022ba35d 65 Double_t dummy = 1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4];
66 if (dummy<0) {
67 AliError(Form("The argument of the Sqrt is %f => set to 0\n",dummy));
68 dummy = 0;
69 }
70 yc = fHelix[0]+TMath::Sqrt(dummy)/fHelix[4];
7f572c00 71
72 fHelix[6] = xc*cs - yc*sn;
73 fHelix[7] = xc*sn + yc*cs;
74 fHelix[8] = TMath::Abs(rc);
75 //
76 //
77 fHelix[5]=x*cs - fHelix[0]*sn; // x0
78 fHelix[0]=x*sn + fHelix[0]*cs; // y0
79 //fHelix[1]= // z0
9be2fe3a 80 fHelix[2]=TMath::ATan2(-(fHelix[5]-fHelix[6]),fHelix[0]-fHelix[7]); // phi0
81 if (fHelix[4]>0) fHelix[2]-=TMath::Pi();
82
7f572c00 83 //fHelix[3]= // tgl
84 //
85 //
86 fHelix[5] = fHelix[6];
87 fHelix[0] = fHelix[7];
7f572c00 88}
89
51ad6848 90
91AliHelix::AliHelix(const AliExternalTrackParam &t)
92{
93 //
94 //
95 Double_t alpha,x,cs,sn;
96 const Double_t *param =t.GetParameter();
97 for (Int_t i=0;i<5;i++) fHelix[i]=param[i];
c9ec41e8 98 x = t.GetX();
99 alpha=t.GetAlpha();
51ad6848 100 //
101 //circle parameters
022ba35d 102 //PH Sometimes fP4 and fHelix[4] are very big and the calculation
103 //PH of the Sqrt cannot be done. To be investigated...
190fc854 104 fHelix[4]=fHelix[4]/(-1000/0.299792458/AliTracker::GetBz()); // C
51ad6848 105 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
106
107 Double_t xc, yc, rc;
108 rc = 1/fHelix[4];
109 xc = x-fHelix[2]*rc;
022ba35d 110 Double_t dummy = 1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4];
111 if (dummy<0) {
112 AliError(Form("The argument of the Sqrt is %f => set to 0\n",dummy));
113 dummy = 0;
114 }
115 yc = fHelix[0]+TMath::Sqrt(dummy)/fHelix[4];
51ad6848 116
117 fHelix[6] = xc*cs - yc*sn;
118 fHelix[7] = xc*sn + yc*cs;
119 fHelix[8] = TMath::Abs(rc);
120 //
121 //
122 fHelix[5]=x*cs - fHelix[0]*sn; // x0
123 fHelix[0]=x*sn + fHelix[0]*cs; // y0
124 //fHelix[1]= // z0
125 fHelix[2]=TMath::ASin(fHelix[2]) + alpha; // phi0
126 //fHelix[3]= // tgl
127 //
128 //
129 fHelix[5] = fHelix[6];
130 fHelix[0] = fHelix[7];
51ad6848 131}
132
7f572c00 133AliHelix::AliHelix(Double_t x[3], Double_t p[3], Double_t charge, Double_t conversion)
134{
135 //
136 //
137 //
138 Double_t pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
139 if (TMath::Abs(conversion)<0.00000001)
190fc854 140 conversion = -1000/0.299792458/AliTracker::GetBz();
7f572c00 141 //
142 //
143 fHelix[4] = charge/(conversion*pt); // C
144 fHelix[3] = p[2]/pt; // tgl
145 //
146 Double_t xc, yc, rc;
147 rc = 1/fHelix[4];
148 xc = x[0] -rc*p[1]/pt;
149 yc = x[1] +rc*p[0]/pt;
150 //
151 fHelix[5] = x[0]; // x0
152 fHelix[0] = x[1]; // y0
153 fHelix[1] = x[2]; // z0
154 //
155 fHelix[6] = xc;
156 fHelix[7] = yc;
157 fHelix[8] = TMath::Abs(rc);
158 //
159 fHelix[5]=xc;
160 fHelix[0]=yc;
161 //
162 if (TMath::Abs(p[1])<TMath::Abs(p[0])){
163 fHelix[2]=TMath::ASin(p[1]/pt);
164 if (charge*yc<charge*x[1]) fHelix[2] = TMath::Pi()-fHelix[2];
165 }
166 else{
167 fHelix[2]=TMath::ACos(p[0]/pt);
168 if (charge*xc>charge*x[0]) fHelix[2] = -fHelix[2];
169 }
170
171}
172
81e97e0d 173void AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion, Double_t *xr)
7f572c00 174{
175 // return momentum at given phase
176 Double_t x[3],g[3],gg[3];
177 Evaluate(phase,x,g,gg);
549102ce 178 // if (TMath::Abs(conversion)<0.0001) conversion = -1000/0.299792458/AliTracker::GetBz();
179 if (TMath::Abs(conversion)<0.0001) conversion = TMath::Abs(1./kB2C/AliTracker::GetBz());
180
7f572c00 181 Double_t mt = TMath::Sqrt(g[0]*g[0]+g[1]*g[1]);
182 p[0] = fHelix[8]*g[0]/(mt*conversion);
183 p[1] = fHelix[8]*g[1]/(mt*conversion);
184 p[2] = fHelix[8]*g[2]/(mt*conversion);
81e97e0d 185 if (xr){
186 xr[0] = x[0]; xr[1] = x[1]; xr[2] = x[2];
187 }
7f572c00 188}
189
190void AliHelix::GetAngle(Double_t t1, AliHelix &h, Double_t t2, Double_t angle[3])
191{
192 //
193 //
194 //
195 Double_t x1[3],g1[3],gg1[3];
196 Double_t x2[3],g2[3],gg2[3];
197 Evaluate(t1,x1,g1,gg1);
198 h.Evaluate(t2,x2,g2,gg2);
199
200 //
201 Double_t norm1r = g1[0]*g1[0]+g1[1]*g1[1];
202 Double_t norm1 = TMath::Sqrt(norm1r+g1[2]*g1[2]);
203 norm1r = TMath::Sqrt(norm1r);
204 //
205 Double_t norm2r = g2[0]*g2[0]+g2[1]*g2[1];
206 Double_t norm2 = TMath::Sqrt(norm2r+g2[2]*g2[2]);
207 norm2r = TMath::Sqrt(norm2r);
208 //
51ad6848 209 angle[0] = (g1[0]*g2[0]+g1[1]*g2[1])/(norm1r*norm2r); // angle in phi projection
210 if (TMath::Abs(angle[0])<1.) angle[0] = TMath::ACos(angle[0]);
9be2fe3a 211 else{
212 if (angle[0]>0) angle[0] = 0;
213 if (angle[0]<0) angle[0] = TMath::Pi();
214 }
51ad6848 215 //
216 angle[1] = ((norm1r*norm2r)+g1[2]*g2[2])/(norm1*norm2); // angle in rz projection
217 if (TMath::Abs(angle[1])<1.) angle[1] = TMath::ACos(angle[1]);
9be2fe3a 218 else
219 angle[1]=0;
51ad6848 220
221 angle[2] = (g1[0]*g2[0]+g1[1]*g2[1]+g1[2]*g2[2])/(norm1*norm2); //3D angle
222 if (TMath::Abs(angle[2])<1.) angle[2] = TMath::ACos(angle[2]);
9be2fe3a 223 else
224 angle[2]=0;
7f572c00 225
51ad6848 226
7f572c00 227
228
229}
230
231
232void AliHelix::Evaluate(Double_t t,
233 Double_t r[3], //radius vector
234 Double_t g[3], //first defivatives
235 Double_t gg[3]) //second derivatives
236{
237 //--------------------------------------------------------------------
238 // Calculate position of a point on a track and some derivatives at given phase
239 //--------------------------------------------------------------------
240 Double_t phase=fHelix[4]*t+fHelix[2];
241 Double_t sn=TMath::Sin(phase), cs=TMath::Cos(phase);
242
7f572c00 243 r[0] = fHelix[5] + sn/fHelix[4];
244 r[1] = fHelix[0] - cs/fHelix[4];
245 r[2] = fHelix[1] + fHelix[3]*t;
246
247 g[0] = cs; g[1]=sn; g[2]=fHelix[3];
248
249 gg[0]=-fHelix[4]*sn; gg[1]=fHelix[4]*cs; gg[2]=0.;
250}
251
9be2fe3a 252Int_t AliHelix::GetClosestPhases(AliHelix &h, Double_t phase[2][2])
253{
254 //
255 // get phases to minimize distances
256 //
257 Double_t xyz0[3];
258 Double_t xyz1[3];
259
260 for (Int_t i=0;i<2;i++){
261 Evaluate(phase[i][0] ,xyz0);
262 h.Evaluate(phase[i][1],xyz1);
263 Double_t mindist = TMath::Sqrt((xyz0[0]-xyz1[0])*(xyz0[0]-xyz1[0])+
264 (xyz0[1]-xyz1[1])*(xyz0[1]-xyz1[1])+
265 (xyz0[2]-xyz1[2])*(xyz0[2]-xyz1[2]));
266 Double_t tbest[2]={phase[i][0],phase[i][1]};
267 for (Int_t i0=-1;i0<=1;i0++){
268 Double_t t0 = ((phase[i][0]*fHelix[4])+i0*2.*TMath::Pi())/fHelix[4];
269 Evaluate(t0,xyz0);
270 for (Int_t i1=-1;i1<=1;i1++){
271 Double_t t1 = ((phase[i][1]*h.fHelix[4])+i1*2.*TMath::Pi())/h.fHelix[4];
272 h.Evaluate(t1,xyz1);
273 Double_t dist = TMath::Sqrt((xyz0[0]-xyz1[0])*(xyz0[0]-xyz1[0])+
274 (xyz0[1]-xyz1[1])*(xyz0[1]-xyz1[1])+
275 (xyz0[2]-xyz1[2])*(xyz0[2]-xyz1[2]));
276 if (dist<=mindist){
277 tbest[0] = t0;
278 tbest[1] = t1;
279 mindist=dist;
280 }
281 }
282 }
283 phase[i][0] = tbest[0];
284 phase[i][1] = tbest[1];
285 }
286 return 1;
287}
288
81e97e0d 289Double_t AliHelix::GetPointAngle(AliHelix &h, Double_t phase[2], const Float_t * vertex)
290{
291 //
292 // get point angle bettwen two helixes
293 //
294 Double_t r0[3],p0[4];
295 Double_t r1[3],p1[4];
296 GetMomentum(phase[0],p0,1,r0);
297 h.GetMomentum(phase[1],p1,1,r1);
298 //
299 Double_t r[3] = {(r0[0]+r1[0])*0.5-vertex[0],(r0[1]+r1[1])*0.5-vertex[1],(r0[2]+r1[2])*0.5-vertex[2]};
300 //intersection point - relative to the prim vertex
301 Double_t p[3] = { p0[0]+p1[0], p0[1]+p1[1],p0[2]+p1[2]};
302 // derivation vector
303 Double_t normr = TMath::Sqrt(r[0]*r[0]+r[1]*r[1]+r[2]*r[2]);
304 Double_t normp = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
305 Double_t pointAngle = (r[0]*p[0]+r[1]*p[1]+r[2]*p[2])/(normr*normp);
306 return pointAngle;
307}
308
7f572c00 309Double_t AliHelix::GetPhase(Double_t x, Double_t y )
310
311{
312 //
313 //calculate helix param at given x,y point
314 //
9be2fe3a 315 //Double_t phase2 = TMath::ATan2((y-fHelix[0]), (x-fHelix[5]))- TMath::Pi()/2.;
316 Double_t phase2 = TMath::ATan2(-(x-fHelix[5]),(y-fHelix[0]));
317 Int_t sign = (fHelix[4]>0)? 1:-1;
318 if (sign>0) phase2 = phase2-TMath::Pi();
319 //
320 Float_t delta = TMath::Nint((phase2-fHelix[2])/(2.*TMath::Pi()));
321 phase2-= 2*TMath::Pi()*delta;
322 if ( (phase2-fHelix[2])>TMath::Pi()) phase2 -=2.*TMath::Pi();
323 if ( (phase2-fHelix[2])<-TMath::Pi()) phase2+=2.*TMath::Pi();
7f572c00 324
9be2fe3a 325 Double_t t = (phase2-fHelix[2]);
326 t/=fHelix[4];
7f572c00 327 return t;
328}
329
176aff27 330Int_t AliHelix::GetPhase(Double_t /*r0*/, Double_t * /*t[2]*/)
7f572c00 331{
332 //
333 //calculate helix param at given r point - return nearest point ()
334 //
335 // not implemented yet
336
337
338 return 0;
339}
340
341
342Double_t AliHelix::GetPhaseZ(Double_t z0)
343{
344 //
345 //
346 return (z0-fHelix[1])/fHelix[3];
347}
348
349
350Int_t AliHelix::GetRPHIintersections(AliHelix &h, Double_t phase[2][2], Double_t ri[2], Double_t cut)
351{
352 //--------------------------------------------------------------------
353 // This function returns phase vectors with intesection between helix (0, 1 or 2)
354 // in x-y plane projection
355 //--------------------------------------------------------------------
356 //
357 // Double_t * c1 = &fHelix[6];
358 //Double_t * c2 = &(h.fHelix[6]);
359 // Double_t c1[3] = {fHelix[5],fHelix[0],fHelix[8]};
9e8f4343 360
361 // PH initiaziation in case of return
362 phase[0][0]=phase[0][1]=phase[1][0]=phase[1][1]=0;
363 ri[0]=ri[1]=1000000;
364
7f572c00 365 Double_t c1[3] = {0,0,fHelix[8]};
366 Double_t c2[3] = {h.fHelix[5]-fHelix[5],h.fHelix[0]-fHelix[0],h.fHelix[8]};
367
368 Double_t d = TMath::Sqrt(c2[0]*c2[0]+c2[1]*c2[1]);
51ad6848 369 if (d<0.000000000001) return 0;
7f572c00 370 //
371 Double_t x0[2];
372 Double_t y0[2];
373 //
374 if ( d>=(c1[2]+c2[2])){
375 if (d>=(c1[2]+c2[2]+cut)) return 0;
376 x0[0] = (d+c1[2]-c2[2])*c2[0]/(2*d)+ fHelix[5];
377 y0[0] = (d+c1[2]-c2[2])*c2[1]/(2*d)+ fHelix[0];
51ad6848 378 // return 0;
9e8f4343 379 phase[1][0] = phase[0][0] = GetPhase(x0[0],y0[0]);
380 phase[1][1] = phase[0][1] = h.GetPhase(x0[0],y0[0]);
381 ri[1] = ri[0] = x0[0]*x0[0]+y0[0]*y0[0];
51ad6848 382 return 1;
7f572c00 383 }
384 if ( (d+c2[2])<c1[2]){
385 if ( (d+c2[2])+cut<c1[2]) return 0;
386 //
387 Double_t xx = c2[0]+ c2[0]*c2[2]/d+ fHelix[5];
388 Double_t yy = c2[1]+ c2[1]*c2[2]/d+ fHelix[0];
9e8f4343 389 phase[1][1] = phase[0][1] = h.GetPhase(xx,yy);
7f572c00 390 //
391 Double_t xx2 = c2[0]*c1[2]/d+ fHelix[5];
392 Double_t yy2 = c2[1]*c1[2]/d+ fHelix[0];
9e8f4343 393 phase[1][0] = phase[0][0] = GetPhase(xx2,yy2);
394 ri[1] = ri[0] = xx*xx+yy*yy;
7f572c00 395 return 1;
396 }
397
398 if ( (d+c1[2])<c2[2]){
399 if ( (d+c1[2])+cut<c2[2]) return 0;
400 //
401 Double_t xx = -c2[0]*c1[2]/d+ fHelix[5];
402 Double_t yy = -c2[1]*c1[2]/d+ fHelix[0];
9e8f4343 403 phase[1][1] = phase[0][1] = GetPhase(xx,yy);
7f572c00 404 //
405 Double_t xx2 = c2[0]- c2[0]*c2[2]/d+ fHelix[5];
406 Double_t yy2 = c2[1]- c2[1]*c2[2]/d+ fHelix[0];
9e8f4343 407 phase[1][0] = phase[0][0] = h.GetPhase(xx2,yy2);
408 ri[1] = ri[0] = xx*xx+yy*yy;
7f572c00 409 return 1;
410 }
411
412 Double_t d1 = (d*d+c1[2]*c1[2]-c2[2]*c2[2])/(2.*d);
413 Double_t v1 = c1[2]*c1[2]-d1*d1;
414 if (v1<0) return 0;
415 v1 = TMath::Sqrt(v1);
416 //
417 x0[0] = (c2[0]*d1+c2[1]*v1)/d + fHelix[5];
418 y0[0] = (c2[1]*d1-c2[0]*v1)/d + fHelix[0];
419 //
420 x0[1] = (c2[0]*d1-c2[1]*v1)/d + fHelix[5];
421 y0[1] = (c2[1]*d1+c2[0]*v1)/d + fHelix[0];
422 //
423 for (Int_t i=0;i<2;i++){
424 phase[i][0] = GetPhase(x0[i],y0[i]);
425 phase[i][1] = h.GetPhase(x0[i],y0[i]);
426 ri[i] = x0[i]*x0[i]+y0[i]*y0[i];
427 }
428 return 2;
429}
430
7f572c00 431
432
433Int_t AliHelix::LinearDCA(AliHelix &h, Double_t &t1, Double_t &t2,
434 Double_t &R, Double_t &dist)
435{
436 //
437 //
438 // find intersection using linear approximation
439 Double_t r1[3],g1[3],gg1[3];
440 Double_t r2[3],g2[3],gg2[3];
441 //
442 Evaluate(t1,r1,g1,gg1);
443 h.Evaluate(t2,r2,g2,gg2);
444 //
445 Double_t g1_2 = g1[0]*g1[0] +g1[1]*g1[1] +g1[2]*g1[2];
446 Double_t g2_2 = g2[0]*g2[0] +g2[1]*g2[1] +g2[2]*g2[2];
447 Double_t g1x2 = g1[0]*g2[0] +g1[1]*g2[1] +g1[2]*g2[2];
448 Double_t det = g1_2*g2_2 - g1x2*g1x2;
449 //
450 if (TMath::Abs(det)>0){
451 //
452 Double_t r1g1 = r1[0]*g1[0] +r1[1]*g1[1] +r1[2]*g1[2];
453 Double_t r2g1 = r2[0]*g1[0] +r2[1]*g1[1] +r2[2]*g1[2];
454 Double_t r1g2 = r1[0]*g2[0] +r1[1]*g2[1] +r1[2]*g2[2];
455 Double_t r2g2 = r2[0]*g2[0] +r2[1]*g2[1] +r2[2]*g2[2];
456 //
457 Double_t dt = - ( g2_2*(r1g1-r2g1) - g1x2*(r1g2-r2g2)) / det;
458 Double_t dp = - ( g1_2*(r2g2-r1g2) - g1x2*(r2g1-r1g1)) / det;
459 //
460 t1+=dt;
461 t2+=dp;
462 Evaluate(t1,r1);
463 h.Evaluate(t2,r2);
464 //
465 dist = (r1[0]-r2[0])*(r1[0]-r2[0])+
466 (r1[1]-r2[1])*(r1[1]-r2[1])+
467 (r1[2]-r2[2])*(r1[2]-r2[2]);
468 R = ((r1[0]+r2[0])*(r1[0]+r2[0])+(r1[1]+r2[1])*(r1[1]+r2[1]))/4.;
469 }
470 return 0;
471}
472
473
7f572c00 474Int_t AliHelix::ParabolicDCA(AliHelix&h, //helixes
475 Double_t &t1, Double_t &t2,
476 Double_t &R, Double_t &dist, Int_t iter)
477{
478 //
479 //
480 // find intersection using linear fit
481 Double_t r1[3],g1[3],gg1[3];
482 Double_t r2[3],g2[3],gg2[3];
483 //
484 Evaluate(t1,r1,g1,gg1);
485 h.Evaluate(t2,r2,g2,gg2);
486
487 //
488 Double_t dx2=1.;
489 Double_t dy2=1.;
490 Double_t dz2=1.;
491 //
492 Double_t dx=r2[0]-r1[0], dy=r2[1]-r1[1], dz=r2[2]-r1[2];
493 Double_t dm=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
494 //
495
496 iter++;
497 while (iter--) {
9be2fe3a 498 Double_t gt1=-(dx*g1[0]/dx2 + dy*g1[1]/dy2 + dz*g1[2]/dz2);
499 Double_t gt2=+(dx*g2[0]/dx2 + dy*g2[1]/dy2 + dz*g2[2]/dz2);
500
501 Double_t h11=(g1[0]*g1[0] - dx*gg1[0])/dx2 +
502 (g1[1]*g1[1] - dy*gg1[1])/dy2 +
503 (g1[2]*g1[2] - dz*gg1[2])/dz2;
504 Double_t h22=(g2[0]*g2[0] + dx*gg2[0])/dx2 +
505 (g2[1]*g2[1] + dy*gg2[1])/dy2 +
506 (g2[2]*g2[2] + dz*gg2[2])/dz2;
507 Double_t h12=-(g1[0]*g2[0]/dx2 + g1[1]*g2[1]/dy2 + g1[2]*g2[2]/dz2);
508
509 Double_t det=h11*h22-h12*h12;
510
511 Double_t dt1,dt2;
512 if (TMath::Abs(det)<1.e-33) {
513 //(quasi)singular Hessian
514 dt1=-gt1; dt2=-gt2;
515 } else {
516 dt1=-(gt1*h22 - gt2*h12)/det;
517 dt2=-(h11*gt2 - h12*gt1)/det;
518 }
519
520 if ((dt1*gt1+dt2*gt2)>0) {dt1=-dt1; dt2=-dt2;}
521
522 //if (TMath::Abs(dt1)/(TMath::Abs(t1)+1.e-3) < 1.e-4)
523 // if (TMath::Abs(dt2)/(TMath::Abs(t2)+1.e-3) < 1.e-4) {
524 // break;
525 // }
526
527 Double_t dd=dm;
528 for (Int_t div=1 ; div<512 ; div*=2) {
529 Evaluate(t1+dt1,r1,g1,gg1);
530 h.Evaluate(t2+dt2,r2,g2,gg2);
531 dx=r2[0]-r1[0]; dy=r2[1]-r1[1]; dz=r2[2]-r1[2];
532 dd=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
533 if (dd<dm) break;
534 dt1*=0.5; dt2*=0.5;
cd507f9c 535 // if (div==0){
536// div =1;
537// }
538// if (div>512) {
539// break;
540// }
9be2fe3a 541 }
542 dm=dd;
543 t1+=dt1;
544 t2+=dt2;
7f572c00 545 }
7f572c00 546 Evaluate(t1,r1,g1,gg1);
547 h.Evaluate(t2,r2,g2,gg2);
548 //
549 dist = (r1[0]-r2[0])*(r1[0]-r2[0])+
550 (r1[1]-r2[1])*(r1[1]-r2[1])+
551 (r1[2]-r2[2])*(r1[2]-r2[2]);
552
553 R = ((r1[0]+r2[0])*(r1[0]+r2[0])+(r1[1]+r2[1])*(r1[1]+r2[1]))/4;
9be2fe3a 554 return 0;
7f572c00 555
556}
7f572c00 557
558
9be2fe3a 559Int_t AliHelix::ParabolicDCA2(AliHelix&h, //helixes
7f572c00 560 Double_t &t1, Double_t &t2,
9be2fe3a 561 Double_t &R, Double_t &dist, Double_t err[3], Int_t iter)
7f572c00 562{
563 //
564 //
565 // find intersection using linear fit
566 Double_t r1[3],g1[3],gg1[3];
567 Double_t r2[3],g2[3],gg2[3];
568 //
569 Evaluate(t1,r1,g1,gg1);
570 h.Evaluate(t2,r2,g2,gg2);
571
572 //
9be2fe3a 573 Double_t dx2=err[0];
574 Double_t dy2=err[1];
575 Double_t dz2=err[2];
7f572c00 576 //
577 Double_t dx=r2[0]-r1[0], dy=r2[1]-r1[1], dz=r2[2]-r1[2];
578 Double_t dm=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
579 //
580
581 iter++;
582 while (iter--) {
583 Double_t gt1=-(dx*g1[0]/dx2 + dy*g1[1]/dy2 + dz*g1[2]/dz2);
584 Double_t gt2=+(dx*g2[0]/dx2 + dy*g2[1]/dy2 + dz*g2[2]/dz2);
585
586 Double_t h11=(g1[0]*g1[0] - dx*gg1[0])/dx2 +
587 (g1[1]*g1[1] - dy*gg1[1])/dy2 +
588 (g1[2]*g1[2] - dz*gg1[2])/dz2;
589 Double_t h22=(g2[0]*g2[0] + dx*gg2[0])/dx2 +
590 (g2[1]*g2[1] + dy*gg2[1])/dy2 +
591 (g2[2]*g2[2] + dz*gg2[2])/dz2;
592 Double_t h12=-(g1[0]*g2[0]/dx2 + g1[1]*g2[1]/dy2 + g1[2]*g2[2]/dz2);
593
594 Double_t det=h11*h22-h12*h12;
595
596 Double_t dt1,dt2;
597 if (TMath::Abs(det)<1.e-33) {
598 //(quasi)singular Hessian
599 dt1=-gt1; dt2=-gt2;
600 } else {
601 dt1=-(gt1*h22 - gt2*h12)/det;
602 dt2=-(h11*gt2 - h12*gt1)/det;
603 }
604
605 if ((dt1*gt1+dt2*gt2)>0) {dt1=-dt1; dt2=-dt2;}
606
607 //if (TMath::Abs(dt1)/(TMath::Abs(t1)+1.e-3) < 1.e-4)
608 // if (TMath::Abs(dt2)/(TMath::Abs(t2)+1.e-3) < 1.e-4) {
609 // break;
610 // }
611
612 Double_t dd=dm;
613 for (Int_t div=1 ; div<512 ; div*=2) {
614 Evaluate(t1+dt1,r1,g1,gg1);
615 h.Evaluate(t2+dt2,r2,g2,gg2);
616 dx=r2[0]-r1[0]; dy=r2[1]-r1[1]; dz=r2[2]-r1[2];
617 dd=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
618 if (dd<dm) break;
619 dt1*=0.5; dt2*=0.5;
620 if (div==0){
621 div =1;
622 }
623 if (div>512) {
624 break;
625 }
626 }
627 dm=dd;
628 t1+=dt1;
629 t2+=dt2;
630 }
631 Evaluate(t1,r1,g1,gg1);
632 h.Evaluate(t2,r2,g2,gg2);
633 //
634 dist = (r1[0]-r2[0])*(r1[0]-r2[0])+
635 (r1[1]-r2[1])*(r1[1]-r2[1])+
636 (r1[2]-r2[2])*(r1[2]-r2[2]);
637
638 R = ((r1[0]+r2[0])*(r1[0]+r2[0])+(r1[1]+r2[1])*(r1[1]+r2[1]))/4;
639 return 0;
640
641}
642