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