STEER/AliHelix.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 /* $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"
  26 #include "AliExternalTrackParam.h"
  27 #include "TMath.h"
  28 ClassImp(AliHelix)
  29
  30
  31 //_______________________________________________________________________
  32 AliHelix::AliHelix()
  33 {
  34   //
  35   // Default constructor
  36   //
  37   for (Int_t i =0;i<9;i++) fHelix[i]=0;
  38 }
  39
  40 //_______________________________________________________________________
  41 AliHelix::AliHelix(const AliHelix &t):TObject(t){
  42   //
  43   //
  44   for (Int_t i=0;i<9;i++)
  45     fHelix[i]=t.fHelix[i];
  46 }
  47
  48 AliHelix::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
  73   fHelix[2]=TMath::ATan2(-(fHelix[5]-fHelix[6]),fHelix[0]-fHelix[7]); // phi0
  74   if (fHelix[4]>0) fHelix[2]-=TMath::Pi();
  75
  76   //fHelix[3]=                               // tgl
  77   //
  78   //
  79   fHelix[5]   = fHelix[6];
  80   fHelix[0]   = fHelix[7];
  81 }
  82
  83
  84 AliHelix::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];
 117 }
 118
 119 AliHelix::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
 159 void  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
 171 void   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   //
 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]);
 192   else{
 193     if (angle[0]>0) angle[0] = 0;
 194     if (angle[0]<0) angle[0] = TMath::Pi();
 195   }
 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]);
 199   else
 200     angle[1]=0;
 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]);
 204   else
 205     angle[2]=0;
 206
 207
 208
 209
 210 }
 211
 212
 213 void 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
 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
 233 Int_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
 270 Double_t  AliHelix::GetPhase(Double_t x, Double_t y )
 271
 272 {
 273   //
 274   //calculate helix param at given x,y  point
 275   //
 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();
 285
 286   Double_t t     = (phase2-fHelix[2]);
 287   t/=fHelix[4];
 288   return t;
 289 }
 290
 291 Int_t AliHelix::GetPhase(Double_t /*r0*/, Double_t * /*t[2]*/)
 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
 303 Double_t  AliHelix::GetPhaseZ(Double_t z0)
 304 {
 305   //
 306   //
 307   return (z0-fHelix[1])/fHelix[3];
 308 }
 309
 310
 311 Int_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]);
 325   if (d<0.000000000001) return 0;
 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];
 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;
 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
 387
 388
 389 Int_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
 430 Int_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--) {
 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;
 501  }
 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;
 510  return 0;
 511
 512 }
 513
 514
 515 Int_t  AliHelix::ParabolicDCA2(AliHelix&h,  //helixes
 516                                Double_t &t1, Double_t &t2,
 517                                Double_t &R, Double_t &dist,  Double_t err[3], Int_t iter)
 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   //
 529   Double_t dx2=err[0];
 530   Double_t dy2=err[1];
 531   Double_t dz2=err[2];
 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