STEER/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 "AliTracker.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   //PH Sometimes fP4 and fHelix[4] are very big and the calculation
  58   //PH of the Sqrt cannot be done. To be investigated...
  59   fHelix[4]=fHelix[4]/(-1000/0.299792458/AliTracker::GetBz());    // C
  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;
  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];
  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
  80   fHelix[2]=TMath::ATan2(-(fHelix[5]-fHelix[6]),fHelix[0]-fHelix[7]); // phi0
  81   if (fHelix[4]>0) fHelix[2]-=TMath::Pi();
  82
  83   //fHelix[3]=                               // tgl
  84   //
  85   //
  86   fHelix[5]   = fHelix[6];
  87   fHelix[0]   = fHelix[7];
  88 }
  89
  90
  91 AliHelix::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];
  98   x = t.GetX();
  99   alpha=t.GetAlpha();
 100   //
 101   //circle parameters
 102   //PH Sometimes fP4 and fHelix[4] are very big and the calculation
 103   //PH of the Sqrt cannot be done. To be investigated...
 104   fHelix[4]=fHelix[4]/(-1000/0.299792458/AliTracker::GetBz());    // C
 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;
 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];
 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];
 131 }
 132
 133 AliHelix::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)
 140     conversion = -1000/0.299792458/AliTracker::GetBz();
 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
 173 void  AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion, Double_t *xr)
 174 {
 175   // return  momentum at given phase
 176   Double_t x[3],g[3],gg[3];
 177   Evaluate(phase,x,g,gg);
 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
 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);
 185   if (xr){
 186     xr[0] = x[0]; xr[1] = x[1]; xr[2] = x[2];
 187   }
 188 }
 189
 190 void   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   //
 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]);
 211   else{
 212     if (angle[0]>0) angle[0] = 0;
 213     if (angle[0]<0) angle[0] = TMath::Pi();
 214   }
 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]);
 218   else
 219     angle[1]=0;
 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]);
 223   else
 224     angle[2]=0;
 225
 226
 227
 228
 229 }
 230
 231
 232 void 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
 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
 252 Int_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
 289 Double_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
 309 Double_t  AliHelix::GetPhase(Double_t x, Double_t y )
 310
 311 {
 312   //
 313   //calculate helix param at given x,y  point
 314   //
 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();
 324
 325   Double_t t     = (phase2-fHelix[2]);
 326   t/=fHelix[4];
 327   return t;
 328 }
 329
 330 Int_t AliHelix::GetPhase(Double_t /*r0*/, Double_t * /*t[2]*/)
 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
 342 Double_t  AliHelix::GetPhaseZ(Double_t z0)
 343 {
 344   //
 345   //
 346   return (z0-fHelix[1])/fHelix[3];
 347 }
 348
 349
 350 Int_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]};
 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
 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]);
 369   if (d<0.000000000001) return 0;
 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];
 378     //    return 0;
 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];
 382     return 1;
 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];
 389     phase[1][1] = phase[0][1] = h.GetPhase(xx,yy);
 390     //
 391     Double_t xx2 = c2[0]*c1[2]/d+ fHelix[5];
 392     Double_t yy2 = c2[1]*c1[2]/d+ fHelix[0];
 393     phase[1][0] = phase[0][0] = GetPhase(xx2,yy2);
 394     ri[1] = ri[0] = xx*xx+yy*yy;
 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];
 403     phase[1][1] = phase[0][1] = GetPhase(xx,yy);
 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];
 407     phase[1][0] = phase[0][0] = h.GetPhase(xx2,yy2);
 408     ri[1] = ri[0] = xx*xx+yy*yy;
 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
 431
 432
 433 Int_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
 474 Int_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--) {
 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;
 535   //     if (div==0){
 536 //      div =1;
 537 //       }
 538 //       if (div>512) {
 539 //      break;
 540 //       }
 541     }
 542     dm=dd;
 543     t1+=dt1;
 544     t2+=dt2;
 545  }
 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;
 554  return 0;
 555
 556 }
 557
 558
 559 Int_t  AliHelix::ParabolicDCA2(AliHelix&h,  //helixes
 560                                Double_t &t1, Double_t &t2,
 561                                Double_t &R, Double_t &dist,  Double_t err[3], Int_t iter)
 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   //
 573   Double_t dx2=err[0];
 574   Double_t dy2=err[1];
 575   Double_t dz2=err[2];
 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