ROOT logo
/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//-------------------------------------------------------------------------
//                Implementation of the AliHelix class
//        Origin: Marian Ivanov, CERN, marian.ivanov@cern.ch
//-------------------------------------------------------------------------


#include "AliHelix.h"
#include "AliKalmanTrack.h"
#include "AliTracker.h"
#include "TMath.h"
ClassImp(AliHelix)


//_______________________________________________________________________
AliHelix::AliHelix()
{
  //
  // Default constructor
  //
  for (Int_t i =0;i<9;i++) fHelix[i]=0;
}

//_______________________________________________________________________
AliHelix::AliHelix(const AliHelix &t):TObject(t){
  //
  //
  for (Int_t i=0;i<9;i++) 
    fHelix[i]=t.fHelix[i];
}

AliHelix::AliHelix(const AliKalmanTrack &t)
{
  //
  // 
  Double_t alpha,x,cs,sn;
  t.GetExternalParameters(x,fHelix); 
  alpha=t.GetAlpha();
  //
  //circle parameters
  //PH Sometimes fP4 and fHelix[4] are very big and the calculation
  //PH of the Sqrt cannot be done. To be investigated...
  fHelix[4]=fHelix[4]/(-1000/0.299792458/AliTracker::GetBz());    // C
  cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);

  Double_t xc, yc, rc;
  rc  =  1/fHelix[4];
  xc  =  x-fHelix[2]*rc;
  Double_t dummy = 1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4];
  if (dummy<0) {
    AliError(Form("The argument of the Sqrt is %f => set to 0\n",dummy));
    dummy = 0;
  }
  yc  =  fHelix[0]+TMath::Sqrt(dummy)/fHelix[4];
  
  fHelix[6] = xc*cs - yc*sn;
  fHelix[7] = xc*sn + yc*cs;
  fHelix[8] =  TMath::Abs(rc);
  //
  //
  fHelix[5]=x*cs - fHelix[0]*sn;            // x0
  fHelix[0]=x*sn + fHelix[0]*cs;            // y0
  //fHelix[1]=                               // z0
  fHelix[2]=TMath::ATan2(-(fHelix[5]-fHelix[6]),fHelix[0]-fHelix[7]); // phi0
  if (fHelix[4]>0) fHelix[2]-=TMath::Pi();

  //fHelix[3]=                               // tgl
  //
  //
  fHelix[5]   = fHelix[6];
  fHelix[0]   = fHelix[7];
}


AliHelix::AliHelix(const AliExternalTrackParam &t)
{
  //
  // 
  Double_t alpha,x,cs,sn;
  const Double_t *param =t.GetParameter(); 
  for (Int_t i=0;i<5;i++) fHelix[i]=param[i]; 
  x = t.GetX();
  alpha=t.GetAlpha();
  //
  //circle parameters
  //PH Sometimes fP4 and fHelix[4] are very big and the calculation
  //PH of the Sqrt cannot be done. To be investigated...
  fHelix[4]=fHelix[4]/(-1000/0.299792458/AliTracker::GetBz());    // C
  cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);

  Double_t xc, yc, rc;
  rc  =  1/fHelix[4];
  xc  =  x-fHelix[2]*rc;
  Double_t dummy = 1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4];
  if (dummy<0) {
    AliError(Form("The argument of the Sqrt is %f => set to 0\n",dummy));
    dummy = 0;
  }
  yc  =  fHelix[0]+TMath::Sqrt(dummy)/fHelix[4];
  
  fHelix[6] = xc*cs - yc*sn;
  fHelix[7] = xc*sn + yc*cs;
  fHelix[8] =  TMath::Abs(rc);
  //
  //
  fHelix[5]=x*cs - fHelix[0]*sn;            // x0
  fHelix[0]=x*sn + fHelix[0]*cs;            // y0
  //fHelix[1]=                               // z0
  fHelix[2]=TMath::ASin(fHelix[2]) + alpha; // phi0
  //fHelix[3]=                               // tgl
  //
  //
  fHelix[5]   = fHelix[6];
  fHelix[0]   = fHelix[7];
}

AliHelix::AliHelix(Double_t x[3], Double_t p[3], Double_t charge, Double_t conversion)
{
  //
  //
  //
  Double_t pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
  if (TMath::Abs(conversion)<0.00000001) 
    conversion = -1000/0.299792458/AliTracker::GetBz();
  //
  //  
  fHelix[4] = charge/(conversion*pt); // C
  fHelix[3] = p[2]/pt;    // tgl
  //  
  Double_t xc, yc, rc;
  rc  =  1/fHelix[4];
  xc  =  x[0]  -rc*p[1]/pt;
  yc  =  x[1]  +rc*p[0]/pt; 
  //
  fHelix[5] = x[0];   // x0
  fHelix[0] = x[1];   // y0
  fHelix[1] = x[2];   // z0
  //
  fHelix[6] = xc;
  fHelix[7] = yc;
  fHelix[8] = TMath::Abs(rc);
  //
  fHelix[5]=xc; 
  fHelix[0]=yc; 
  //
  if (TMath::Abs(p[1])<TMath::Abs(p[0])){     
    fHelix[2]=TMath::ASin(p[1]/pt);
    if (charge*yc<charge*x[1])  fHelix[2] = TMath::Pi()-fHelix[2];
  }
  else{
    fHelix[2]=TMath::ACos(p[0]/pt);
    if (charge*xc>charge*x[0])  fHelix[2] = -fHelix[2];
  }

}

void  AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion, Double_t *xr)
{
  // return  momentum at given phase
  Double_t x[3],g[3],gg[3];
  Evaluate(phase,x,g,gg);
  //  if (TMath::Abs(conversion)<0.0001) conversion = -1000/0.299792458/AliTracker::GetBz();
  if (TMath::Abs(conversion)<0.0001) conversion = TMath::Abs(1./kB2C/AliTracker::GetBz());

  Double_t mt = TMath::Sqrt(g[0]*g[0]+g[1]*g[1]);
  p[0] = fHelix[8]*g[0]/(mt*conversion);
  p[1] = fHelix[8]*g[1]/(mt*conversion);
  p[2] = fHelix[8]*g[2]/(mt*conversion);
  if (xr){
    xr[0] = x[0]; xr[1] = x[1]; xr[2] = x[2];
  }
}

void   AliHelix::GetAngle(Double_t t1, AliHelix &h, Double_t t2, Double_t angle[3])
{
  //
  //
  //
  Double_t x1[3],g1[3],gg1[3];
  Double_t x2[3],g2[3],gg2[3];
  Evaluate(t1,x1,g1,gg1);
  h.Evaluate(t2,x2,g2,gg2);

  //
  Double_t norm1r = g1[0]*g1[0]+g1[1]*g1[1];
  Double_t norm1  = TMath::Sqrt(norm1r+g1[2]*g1[2]);
  norm1r         = TMath::Sqrt(norm1r);
  //
  Double_t norm2r = g2[0]*g2[0]+g2[1]*g2[1];
  Double_t norm2  = TMath::Sqrt(norm2r+g2[2]*g2[2]);
  norm2r         = TMath::Sqrt(norm2r);
  //
  angle[0]  = (g1[0]*g2[0]+g1[1]*g2[1])/(norm1r*norm2r);   // angle in phi projection
  if (TMath::Abs(angle[0])<1.) angle[0] = TMath::ACos(angle[0]);
  else{ 
    if (angle[0]>0) angle[0] = 0;
    if (angle[0]<0) angle[0] = TMath::Pi();
  }
  //
  angle[1]  = ((norm1r*norm2r)+g1[2]*g2[2])/(norm1*norm2); // angle in rz  projection
  if (TMath::Abs(angle[1])<1.) angle[1] = TMath::ACos(angle[1]);
  else 
    angle[1]=0;

  angle[2]  = (g1[0]*g2[0]+g1[1]*g2[1]+g1[2]*g2[2])/(norm1*norm2); //3D angle
  if (TMath::Abs(angle[2])<1.) angle[2] = TMath::ACos(angle[2]);
  else 
    angle[2]=0;

  
  

}


void AliHelix::Evaluate(Double_t t,
                     Double_t r[3],  //radius vector
                     Double_t g[3],  //first defivatives
                     Double_t gg[3]) //second derivatives
{
  //--------------------------------------------------------------------
  // Calculate position of a point on a track and some derivatives at given phase
  //--------------------------------------------------------------------
  Double_t phase=fHelix[4]*t+fHelix[2];
  Double_t sn=TMath::Sin(phase), cs=TMath::Cos(phase);

  r[0] = fHelix[5] + sn/fHelix[4];
  r[1] = fHelix[0] - cs/fHelix[4];  
  r[2] = fHelix[1] + fHelix[3]*t;

  g[0] = cs; g[1]=sn; g[2]=fHelix[3];
  
  gg[0]=-fHelix[4]*sn; gg[1]=fHelix[4]*cs; gg[2]=0.;
}

Int_t     AliHelix::GetClosestPhases(AliHelix &h, Double_t phase[2][2])
{
  //
  // get phases to minimize distances
  //
  Double_t xyz0[3];
  Double_t xyz1[3];

  for (Int_t i=0;i<2;i++){  
    Evaluate(phase[i][0]  ,xyz0);	
    h.Evaluate(phase[i][1],xyz1);
    Double_t mindist = TMath::Sqrt((xyz0[0]-xyz1[0])*(xyz0[0]-xyz1[0])+
				   (xyz0[1]-xyz1[1])*(xyz0[1]-xyz1[1])+
				   (xyz0[2]-xyz1[2])*(xyz0[2]-xyz1[2]));  
    Double_t tbest[2]={phase[i][0],phase[i][1]};
    for (Int_t i0=-1;i0<=1;i0++){
      Double_t t0 = ((phase[i][0]*fHelix[4])+i0*2.*TMath::Pi())/fHelix[4];
      Evaluate(t0,xyz0);
      for (Int_t i1=-1;i1<=1;i1++){
	Double_t t1 = ((phase[i][1]*h.fHelix[4])+i1*2.*TMath::Pi())/h.fHelix[4];    
	h.Evaluate(t1,xyz1);    
	Double_t dist = TMath::Sqrt((xyz0[0]-xyz1[0])*(xyz0[0]-xyz1[0])+
				    (xyz0[1]-xyz1[1])*(xyz0[1]-xyz1[1])+
				    (xyz0[2]-xyz1[2])*(xyz0[2]-xyz1[2])); 
	if (dist<=mindist){
	  tbest[0] = t0;
	  tbest[1] = t1;
	  mindist=dist;
	}
      }
    }
    phase[i][0] = tbest[0];
    phase[i][1] = tbest[1];
  } 
  return 1;
}

Double_t  AliHelix::GetPointAngle(AliHelix &h, Double_t phase[2], const Float_t * vertex)
{
  //
  // get point angle bettwen two helixes
  // 
  Double_t r0[3],p0[4];
  Double_t r1[3],p1[4];
  GetMomentum(phase[0],p0,1,r0);
  h.GetMomentum(phase[1],p1,1,r1);
  //
  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]};
  //intersection point - relative to the prim vertex
  Double_t p[3] = { p0[0]+p1[0], p0[1]+p1[1],p0[2]+p1[2]};
  // derivation vector
  Double_t normr = TMath::Sqrt(r[0]*r[0]+r[1]*r[1]+r[2]*r[2]);
  Double_t normp = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);   
  Double_t pointAngle = (r[0]*p[0]+r[1]*p[1]+r[2]*p[2])/(normr*normp);
  return pointAngle;
}

Double_t  AliHelix::GetPhase(Double_t x, Double_t y )
			
{
  //
  //calculate helix param at given x,y  point
  //
  //Double_t phase2 = TMath::ATan2((y-fHelix[0]), (x-fHelix[5]))- TMath::Pi()/2.;
  Double_t phase2 = TMath::ATan2(-(x-fHelix[5]),(y-fHelix[0]));
  Int_t sign = (fHelix[4]>0)? 1:-1;
  if (sign>0) phase2 = phase2-TMath::Pi();
  //
  Float_t delta = TMath::Nint((phase2-fHelix[2])/(2.*TMath::Pi()));
  phase2-= 2*TMath::Pi()*delta;
  if ( (phase2-fHelix[2])>TMath::Pi()) phase2 -=2.*TMath::Pi();
  if ( (phase2-fHelix[2])<-TMath::Pi()) phase2+=2.*TMath::Pi();

  Double_t t     = (phase2-fHelix[2]);
  t/=fHelix[4];
  return t;
}

Int_t AliHelix::GetPhase(Double_t /*r0*/, Double_t * /*t[2]*/) 
{
  //
  //calculate helix param at given r  point - return nearest point ()
  //
  // not implemented yet
  

  return 0;
}


Double_t  AliHelix::GetPhaseZ(Double_t z0)
{
  //
  //
  return (z0-fHelix[1])/fHelix[3];
}


Int_t    AliHelix::GetRPHIintersections(AliHelix &h, Double_t phase[2][2], Double_t ri[2], Double_t cut)
{
  //--------------------------------------------------------------------
  // This function returns  phase vectors with intesection between helix (0, 1 or 2)
  // in x-y plane projection  
  //--------------------------------------------------------------------
  //    
  //  Double_t * c1 = &fHelix[6];
  //Double_t * c2 = &(h.fHelix[6]);
  //  Double_t  c1[3] = {fHelix[5],fHelix[0],fHelix[8]};

  // PH initiaziation in case of return
  phase[0][0]=phase[0][1]=phase[1][0]=phase[1][1]=0;
  ri[0]=ri[1]=1000000;

  Double_t  c1[3] = {0,0,fHelix[8]};
  Double_t  c2[3] = {h.fHelix[5]-fHelix[5],h.fHelix[0]-fHelix[0],h.fHelix[8]};

  Double_t d  = TMath::Sqrt(c2[0]*c2[0]+c2[1]*c2[1]); 
  if (d<0.000000000001) return 0;
  //
  Double_t x0[2];
  Double_t y0[2];
  //  
  if ( d>=(c1[2]+c2[2])){
    if (d>=(c1[2]+c2[2]+cut)) return 0;
    x0[0] = (d+c1[2]-c2[2])*c2[0]/(2*d)+ fHelix[5];
    y0[0] = (d+c1[2]-c2[2])*c2[1]/(2*d)+ fHelix[0];
    //    return 0;
    phase[1][0] = phase[0][0] = GetPhase(x0[0],y0[0]);
    phase[1][1] = phase[0][1] = h.GetPhase(x0[0],y0[0]);
    ri[1] = ri[0] = x0[0]*x0[0]+y0[0]*y0[0];
    return 1;
  }
  if ( (d+c2[2])<c1[2]){
    if ( (d+c2[2])+cut<c1[2]) return 0;
    //
    Double_t xx = c2[0]+ c2[0]*c2[2]/d+ fHelix[5];
    Double_t yy = c2[1]+ c2[1]*c2[2]/d+ fHelix[0]; 
    phase[1][1] = phase[0][1] = h.GetPhase(xx,yy);
    //
    Double_t xx2 = c2[0]*c1[2]/d+ fHelix[5];
    Double_t yy2 = c2[1]*c1[2]/d+ fHelix[0]; 
    phase[1][0] = phase[0][0] = GetPhase(xx2,yy2);
    ri[1] = ri[0] = xx*xx+yy*yy;
    return 1;
  }

  if ( (d+c1[2])<c2[2]){
    if ( (d+c1[2])+cut<c2[2]) return 0;
    //
    Double_t xx = -c2[0]*c1[2]/d+ fHelix[5];
    Double_t yy = -c2[1]*c1[2]/d+ fHelix[0]; 
    phase[1][1] = phase[0][1] = GetPhase(xx,yy);
    //
    Double_t xx2 = c2[0]- c2[0]*c2[2]/d+ fHelix[5];
    Double_t yy2 = c2[1]- c2[1]*c2[2]/d+ fHelix[0]; 
    phase[1][0] = phase[0][0] = h.GetPhase(xx2,yy2);
    ri[1] = ri[0] = xx*xx+yy*yy;
    return 1;
  }

  Double_t d1 = (d*d+c1[2]*c1[2]-c2[2]*c2[2])/(2.*d);
  Double_t v1 = c1[2]*c1[2]-d1*d1;
  if (v1<0) return 0;
  v1 = TMath::Sqrt(v1);
  //
  x0[0] = (c2[0]*d1+c2[1]*v1)/d + fHelix[5];
  y0[0] = (c2[1]*d1-c2[0]*v1)/d + fHelix[0];            
  //
  x0[1] = (c2[0]*d1-c2[1]*v1)/d + fHelix[5];
  y0[1] = (c2[1]*d1+c2[0]*v1)/d + fHelix[0];      
  //
  for (Int_t i=0;i<2;i++){
    phase[i][0] = GetPhase(x0[i],y0[i]);
    phase[i][1] = h.GetPhase(x0[i],y0[i]);
    ri[i] = x0[i]*x0[i]+y0[i]*y0[i];    
  }      
  return 2;
} 



Int_t   AliHelix::LinearDCA(AliHelix &h, Double_t &t1, Double_t &t2, 
		      Double_t &R, Double_t &dist)
{
  //
  //
  // find intersection using linear approximation
  Double_t r1[3],g1[3],gg1[3];
  Double_t r2[3],g2[3],gg2[3];
  //
  Evaluate(t1,r1,g1,gg1);
  h.Evaluate(t2,r2,g2,gg2);
  // 
  Double_t g1_2 = g1[0]*g1[0] +g1[1]*g1[1] +g1[2]*g1[2];
  Double_t g2_2 = g2[0]*g2[0] +g2[1]*g2[1] +g2[2]*g2[2];
  Double_t g1x2 = g1[0]*g2[0] +g1[1]*g2[1] +g1[2]*g2[2];  
  Double_t det  = g1_2*g2_2   - g1x2*g1x2;
  //  
  if (TMath::Abs(det)>0){
    //
    Double_t r1g1 = r1[0]*g1[0] +r1[1]*g1[1] +r1[2]*g1[2];    
    Double_t r2g1 = r2[0]*g1[0] +r2[1]*g1[1] +r2[2]*g1[2];      
    Double_t r1g2 = r1[0]*g2[0] +r1[1]*g2[1] +r1[2]*g2[2];
    Double_t r2g2 = r2[0]*g2[0] +r2[1]*g2[1] +r2[2]*g2[2];
    //    
    Double_t dt    = - ( g2_2*(r1g1-r2g1) - g1x2*(r1g2-r2g2)) / det;      
    Double_t dp    = - ( g1_2*(r2g2-r1g2) - g1x2*(r2g1-r1g1)) / det;
    //
    t1+=dt;
    t2+=dp;
    Evaluate(t1,r1);
    h.Evaluate(t2,r2);
    //
    dist = (r1[0]-r2[0])*(r1[0]-r2[0])+
				  (r1[1]-r2[1])*(r1[1]-r2[1])+
				  (r1[2]-r2[2])*(r1[2]-r2[2]);    
    R = ((r1[0]+r2[0])*(r1[0]+r2[0])+(r1[1]+r2[1])*(r1[1]+r2[1]))/4.;
  }     
  return 0;
}


Int_t  AliHelix::ParabolicDCA(AliHelix&h,  //helixes
			       Double_t &t1, Double_t &t2, 
			       Double_t &R, Double_t &dist, Int_t iter)
{
  //
  //
  // find intersection using linear fit
  Double_t r1[3],g1[3],gg1[3];
  Double_t r2[3],g2[3],gg2[3];
  //
  Evaluate(t1,r1,g1,gg1);
  h.Evaluate(t2,r2,g2,gg2);

  //
  Double_t dx2=1.;
  Double_t dy2=1.;
  Double_t dz2=1.;
  //
  Double_t dx=r2[0]-r1[0], dy=r2[1]-r1[1], dz=r2[2]-r1[2];
  Double_t dm=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
  //

 iter++;
 while (iter--) {
    Double_t gt1=-(dx*g1[0]/dx2 + dy*g1[1]/dy2 + dz*g1[2]/dz2);
    Double_t gt2=+(dx*g2[0]/dx2 + dy*g2[1]/dy2 + dz*g2[2]/dz2);
    
    Double_t h11=(g1[0]*g1[0] - dx*gg1[0])/dx2 + 
      (g1[1]*g1[1] - dy*gg1[1])/dy2 +
      (g1[2]*g1[2] - dz*gg1[2])/dz2;
    Double_t h22=(g2[0]*g2[0] + dx*gg2[0])/dx2 + 
      (g2[1]*g2[1] + dy*gg2[1])/dy2 +
      (g2[2]*g2[2] + dz*gg2[2])/dz2;
    Double_t h12=-(g1[0]*g2[0]/dx2 + g1[1]*g2[1]/dy2 + g1[2]*g2[2]/dz2);
    
    Double_t det=h11*h22-h12*h12;
    
    Double_t dt1,dt2;
    if (TMath::Abs(det)<1.e-33) {
      //(quasi)singular Hessian
      dt1=-gt1; dt2=-gt2;
    } else {
      dt1=-(gt1*h22 - gt2*h12)/det; 
      dt2=-(h11*gt2 - h12*gt1)/det;
    }
    
    if ((dt1*gt1+dt2*gt2)>0) {dt1=-dt1; dt2=-dt2;}
    
    //if (TMath::Abs(dt1)/(TMath::Abs(t1)+1.e-3) < 1.e-4)
    //  if (TMath::Abs(dt2)/(TMath::Abs(t2)+1.e-3) < 1.e-4) {
    //	break;
    //  }
    
    Double_t dd=dm;
    for (Int_t div=1 ; div<512 ; div*=2) {
      Evaluate(t1+dt1,r1,g1,gg1);
      h.Evaluate(t2+dt2,r2,g2,gg2);
      dx=r2[0]-r1[0]; dy=r2[1]-r1[1]; dz=r2[2]-r1[2];
      dd=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
      if (dd<dm) break;
      dt1*=0.5; dt2*=0.5;
  //     if (div==0){
// 	div =1;
//       }
//       if (div>512) {	  
// 	break;
//       }   
    }
    dm=dd;
    t1+=dt1;
    t2+=dt2;
 }
 Evaluate(t1,r1,g1,gg1);
 h.Evaluate(t2,r2,g2,gg2);
 //
 dist = (r1[0]-r2[0])*(r1[0]-r2[0])+
   (r1[1]-r2[1])*(r1[1]-r2[1])+
   (r1[2]-r2[2])*(r1[2]-r2[2]);    
 
 R = ((r1[0]+r2[0])*(r1[0]+r2[0])+(r1[1]+r2[1])*(r1[1]+r2[1]))/4;
 return 0;
 
}


Int_t  AliHelix::ParabolicDCA2(AliHelix&h,  //helixes
			       Double_t &t1, Double_t &t2, 
			       Double_t &R, Double_t &dist,  Double_t err[3], Int_t iter)
{
  //
  //
  // find intersection using linear fit
  Double_t r1[3],g1[3],gg1[3];
  Double_t r2[3],g2[3],gg2[3];
  //
  Evaluate(t1,r1,g1,gg1);
  h.Evaluate(t2,r2,g2,gg2);

  //
  Double_t dx2=err[0];
  Double_t dy2=err[1];
  Double_t dz2=err[2];
  //
  Double_t dx=r2[0]-r1[0], dy=r2[1]-r1[1], dz=r2[2]-r1[2];
  Double_t dm=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
  //

 iter++;
 while (iter--) {
    Double_t gt1=-(dx*g1[0]/dx2 + dy*g1[1]/dy2 + dz*g1[2]/dz2);
    Double_t gt2=+(dx*g2[0]/dx2 + dy*g2[1]/dy2 + dz*g2[2]/dz2);
    
    Double_t h11=(g1[0]*g1[0] - dx*gg1[0])/dx2 + 
      (g1[1]*g1[1] - dy*gg1[1])/dy2 +
      (g1[2]*g1[2] - dz*gg1[2])/dz2;
    Double_t h22=(g2[0]*g2[0] + dx*gg2[0])/dx2 + 
      (g2[1]*g2[1] + dy*gg2[1])/dy2 +
      (g2[2]*g2[2] + dz*gg2[2])/dz2;
    Double_t h12=-(g1[0]*g2[0]/dx2 + g1[1]*g2[1]/dy2 + g1[2]*g2[2]/dz2);
    
    Double_t det=h11*h22-h12*h12;
    
    Double_t dt1,dt2;
    if (TMath::Abs(det)<1.e-33) {
      //(quasi)singular Hessian
      dt1=-gt1; dt2=-gt2;
    } else {
      dt1=-(gt1*h22 - gt2*h12)/det; 
      dt2=-(h11*gt2 - h12*gt1)/det;
    }
    
    if ((dt1*gt1+dt2*gt2)>0) {dt1=-dt1; dt2=-dt2;}
    
    //if (TMath::Abs(dt1)/(TMath::Abs(t1)+1.e-3) < 1.e-4)
    //  if (TMath::Abs(dt2)/(TMath::Abs(t2)+1.e-3) < 1.e-4) {
    //	break;
    //  }
    
    Double_t dd=dm;
    for (Int_t div=1 ; div<512 ; div*=2) {
      Evaluate(t1+dt1,r1,g1,gg1);
      h.Evaluate(t2+dt2,r2,g2,gg2);
      dx=r2[0]-r1[0]; dy=r2[1]-r1[1]; dz=r2[2]-r1[2];
      dd=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
      if (dd<dm) break;
      dt1*=0.5; dt2*=0.5;
      if (div==0){
	div =1;
      }
      if (div>512) {	  
	break;
      }   
    }
    dm=dd;
    t1+=dt1;
    t2+=dt2;
 }
 Evaluate(t1,r1,g1,gg1);
 h.Evaluate(t2,r2,g2,gg2);
 //
 dist = (r1[0]-r2[0])*(r1[0]-r2[0])+
   (r1[1]-r2[1])*(r1[1]-r2[1])+
   (r1[2]-r2[2])*(r1[2]-r2[2]);    
 
 R = ((r1[0]+r2[0])*(r1[0]+r2[0])+(r1[1]+r2[1])*(r1[1]+r2[1]))/4;
 return 0;
 
}

 AliHelix.cxx:1
 AliHelix.cxx:2
 AliHelix.cxx:3
 AliHelix.cxx:4
 AliHelix.cxx:5
 AliHelix.cxx:6
 AliHelix.cxx:7
 AliHelix.cxx:8
 AliHelix.cxx:9
 AliHelix.cxx:10
 AliHelix.cxx:11
 AliHelix.cxx:12
 AliHelix.cxx:13
 AliHelix.cxx:14
 AliHelix.cxx:15
 AliHelix.cxx:16
 AliHelix.cxx:17
 AliHelix.cxx:18
 AliHelix.cxx:19
 AliHelix.cxx:20
 AliHelix.cxx:21
 AliHelix.cxx:22
 AliHelix.cxx:23
 AliHelix.cxx:24
 AliHelix.cxx:25
 AliHelix.cxx:26
 AliHelix.cxx:27
 AliHelix.cxx:28
 AliHelix.cxx:29
 AliHelix.cxx:30
 AliHelix.cxx:31
 AliHelix.cxx:32
 AliHelix.cxx:33
 AliHelix.cxx:34
 AliHelix.cxx:35
 AliHelix.cxx:36
 AliHelix.cxx:37
 AliHelix.cxx:38
 AliHelix.cxx:39
 AliHelix.cxx:40
 AliHelix.cxx:41
 AliHelix.cxx:42
 AliHelix.cxx:43
 AliHelix.cxx:44
 AliHelix.cxx:45
 AliHelix.cxx:46
 AliHelix.cxx:47
 AliHelix.cxx:48
 AliHelix.cxx:49
 AliHelix.cxx:50
 AliHelix.cxx:51
 AliHelix.cxx:52
 AliHelix.cxx:53
 AliHelix.cxx:54
 AliHelix.cxx:55
 AliHelix.cxx:56
 AliHelix.cxx:57
 AliHelix.cxx:58
 AliHelix.cxx:59
 AliHelix.cxx:60
 AliHelix.cxx:61
 AliHelix.cxx:62
 AliHelix.cxx:63
 AliHelix.cxx:64
 AliHelix.cxx:65
 AliHelix.cxx:66
 AliHelix.cxx:67
 AliHelix.cxx:68
 AliHelix.cxx:69
 AliHelix.cxx:70
 AliHelix.cxx:71
 AliHelix.cxx:72
 AliHelix.cxx:73
 AliHelix.cxx:74
 AliHelix.cxx:75
 AliHelix.cxx:76
 AliHelix.cxx:77
 AliHelix.cxx:78
 AliHelix.cxx:79
 AliHelix.cxx:80
 AliHelix.cxx:81
 AliHelix.cxx:82
 AliHelix.cxx:83
 AliHelix.cxx:84
 AliHelix.cxx:85
 AliHelix.cxx:86
 AliHelix.cxx:87
 AliHelix.cxx:88
 AliHelix.cxx:89
 AliHelix.cxx:90
 AliHelix.cxx:91
 AliHelix.cxx:92
 AliHelix.cxx:93
 AliHelix.cxx:94
 AliHelix.cxx:95
 AliHelix.cxx:96
 AliHelix.cxx:97
 AliHelix.cxx:98
 AliHelix.cxx:99
 AliHelix.cxx:100
 AliHelix.cxx:101
 AliHelix.cxx:102
 AliHelix.cxx:103
 AliHelix.cxx:104
 AliHelix.cxx:105
 AliHelix.cxx:106
 AliHelix.cxx:107
 AliHelix.cxx:108
 AliHelix.cxx:109
 AliHelix.cxx:110
 AliHelix.cxx:111
 AliHelix.cxx:112
 AliHelix.cxx:113
 AliHelix.cxx:114
 AliHelix.cxx:115
 AliHelix.cxx:116
 AliHelix.cxx:117
 AliHelix.cxx:118
 AliHelix.cxx:119
 AliHelix.cxx:120
 AliHelix.cxx:121
 AliHelix.cxx:122
 AliHelix.cxx:123
 AliHelix.cxx:124
 AliHelix.cxx:125
 AliHelix.cxx:126
 AliHelix.cxx:127
 AliHelix.cxx:128
 AliHelix.cxx:129
 AliHelix.cxx:130
 AliHelix.cxx:131
 AliHelix.cxx:132
 AliHelix.cxx:133
 AliHelix.cxx:134
 AliHelix.cxx:135
 AliHelix.cxx:136
 AliHelix.cxx:137
 AliHelix.cxx:138
 AliHelix.cxx:139
 AliHelix.cxx:140
 AliHelix.cxx:141
 AliHelix.cxx:142
 AliHelix.cxx:143
 AliHelix.cxx:144
 AliHelix.cxx:145
 AliHelix.cxx:146
 AliHelix.cxx:147
 AliHelix.cxx:148
 AliHelix.cxx:149
 AliHelix.cxx:150
 AliHelix.cxx:151
 AliHelix.cxx:152
 AliHelix.cxx:153
 AliHelix.cxx:154
 AliHelix.cxx:155
 AliHelix.cxx:156
 AliHelix.cxx:157
 AliHelix.cxx:158
 AliHelix.cxx:159
 AliHelix.cxx:160
 AliHelix.cxx:161
 AliHelix.cxx:162
 AliHelix.cxx:163
 AliHelix.cxx:164
 AliHelix.cxx:165
 AliHelix.cxx:166
 AliHelix.cxx:167
 AliHelix.cxx:168
 AliHelix.cxx:169
 AliHelix.cxx:170
 AliHelix.cxx:171
 AliHelix.cxx:172
 AliHelix.cxx:173
 AliHelix.cxx:174
 AliHelix.cxx:175
 AliHelix.cxx:176
 AliHelix.cxx:177
 AliHelix.cxx:178
 AliHelix.cxx:179
 AliHelix.cxx:180
 AliHelix.cxx:181
 AliHelix.cxx:182
 AliHelix.cxx:183
 AliHelix.cxx:184
 AliHelix.cxx:185
 AliHelix.cxx:186
 AliHelix.cxx:187
 AliHelix.cxx:188
 AliHelix.cxx:189
 AliHelix.cxx:190
 AliHelix.cxx:191
 AliHelix.cxx:192
 AliHelix.cxx:193
 AliHelix.cxx:194
 AliHelix.cxx:195
 AliHelix.cxx:196
 AliHelix.cxx:197
 AliHelix.cxx:198
 AliHelix.cxx:199
 AliHelix.cxx:200
 AliHelix.cxx:201
 AliHelix.cxx:202
 AliHelix.cxx:203
 AliHelix.cxx:204
 AliHelix.cxx:205
 AliHelix.cxx:206
 AliHelix.cxx:207
 AliHelix.cxx:208
 AliHelix.cxx:209
 AliHelix.cxx:210
 AliHelix.cxx:211
 AliHelix.cxx:212
 AliHelix.cxx:213
 AliHelix.cxx:214
 AliHelix.cxx:215
 AliHelix.cxx:216
 AliHelix.cxx:217
 AliHelix.cxx:218
 AliHelix.cxx:219
 AliHelix.cxx:220
 AliHelix.cxx:221
 AliHelix.cxx:222
 AliHelix.cxx:223
 AliHelix.cxx:224
 AliHelix.cxx:225
 AliHelix.cxx:226
 AliHelix.cxx:227
 AliHelix.cxx:228
 AliHelix.cxx:229
 AliHelix.cxx:230
 AliHelix.cxx:231
 AliHelix.cxx:232
 AliHelix.cxx:233
 AliHelix.cxx:234
 AliHelix.cxx:235
 AliHelix.cxx:236
 AliHelix.cxx:237
 AliHelix.cxx:238
 AliHelix.cxx:239
 AliHelix.cxx:240
 AliHelix.cxx:241
 AliHelix.cxx:242
 AliHelix.cxx:243
 AliHelix.cxx:244
 AliHelix.cxx:245
 AliHelix.cxx:246
 AliHelix.cxx:247
 AliHelix.cxx:248
 AliHelix.cxx:249
 AliHelix.cxx:250
 AliHelix.cxx:251
 AliHelix.cxx:252
 AliHelix.cxx:253
 AliHelix.cxx:254
 AliHelix.cxx:255
 AliHelix.cxx:256
 AliHelix.cxx:257
 AliHelix.cxx:258
 AliHelix.cxx:259
 AliHelix.cxx:260
 AliHelix.cxx:261
 AliHelix.cxx:262
 AliHelix.cxx:263
 AliHelix.cxx:264
 AliHelix.cxx:265
 AliHelix.cxx:266
 AliHelix.cxx:267
 AliHelix.cxx:268
 AliHelix.cxx:269
 AliHelix.cxx:270
 AliHelix.cxx:271
 AliHelix.cxx:272
 AliHelix.cxx:273
 AliHelix.cxx:274
 AliHelix.cxx:275
 AliHelix.cxx:276
 AliHelix.cxx:277
 AliHelix.cxx:278
 AliHelix.cxx:279
 AliHelix.cxx:280
 AliHelix.cxx:281
 AliHelix.cxx:282
 AliHelix.cxx:283
 AliHelix.cxx:284
 AliHelix.cxx:285
 AliHelix.cxx:286
 AliHelix.cxx:287
 AliHelix.cxx:288
 AliHelix.cxx:289
 AliHelix.cxx:290
 AliHelix.cxx:291
 AliHelix.cxx:292
 AliHelix.cxx:293
 AliHelix.cxx:294
 AliHelix.cxx:295
 AliHelix.cxx:296
 AliHelix.cxx:297
 AliHelix.cxx:298
 AliHelix.cxx:299
 AliHelix.cxx:300
 AliHelix.cxx:301
 AliHelix.cxx:302
 AliHelix.cxx:303
 AliHelix.cxx:304
 AliHelix.cxx:305
 AliHelix.cxx:306
 AliHelix.cxx:307
 AliHelix.cxx:308
 AliHelix.cxx:309
 AliHelix.cxx:310
 AliHelix.cxx:311
 AliHelix.cxx:312
 AliHelix.cxx:313
 AliHelix.cxx:314
 AliHelix.cxx:315
 AliHelix.cxx:316
 AliHelix.cxx:317
 AliHelix.cxx:318
 AliHelix.cxx:319
 AliHelix.cxx:320
 AliHelix.cxx:321
 AliHelix.cxx:322
 AliHelix.cxx:323
 AliHelix.cxx:324
 AliHelix.cxx:325
 AliHelix.cxx:326
 AliHelix.cxx:327
 AliHelix.cxx:328
 AliHelix.cxx:329
 AliHelix.cxx:330
 AliHelix.cxx:331
 AliHelix.cxx:332
 AliHelix.cxx:333
 AliHelix.cxx:334
 AliHelix.cxx:335
 AliHelix.cxx:336
 AliHelix.cxx:337
 AliHelix.cxx:338
 AliHelix.cxx:339
 AliHelix.cxx:340
 AliHelix.cxx:341
 AliHelix.cxx:342
 AliHelix.cxx:343
 AliHelix.cxx:344
 AliHelix.cxx:345
 AliHelix.cxx:346
 AliHelix.cxx:347
 AliHelix.cxx:348
 AliHelix.cxx:349
 AliHelix.cxx:350
 AliHelix.cxx:351
 AliHelix.cxx:352
 AliHelix.cxx:353
 AliHelix.cxx:354
 AliHelix.cxx:355
 AliHelix.cxx:356
 AliHelix.cxx:357
 AliHelix.cxx:358
 AliHelix.cxx:359
 AliHelix.cxx:360
 AliHelix.cxx:361
 AliHelix.cxx:362
 AliHelix.cxx:363
 AliHelix.cxx:364
 AliHelix.cxx:365
 AliHelix.cxx:366
 AliHelix.cxx:367
 AliHelix.cxx:368
 AliHelix.cxx:369
 AliHelix.cxx:370
 AliHelix.cxx:371
 AliHelix.cxx:372
 AliHelix.cxx:373
 AliHelix.cxx:374
 AliHelix.cxx:375
 AliHelix.cxx:376
 AliHelix.cxx:377
 AliHelix.cxx:378
 AliHelix.cxx:379
 AliHelix.cxx:380
 AliHelix.cxx:381
 AliHelix.cxx:382
 AliHelix.cxx:383
 AliHelix.cxx:384
 AliHelix.cxx:385
 AliHelix.cxx:386
 AliHelix.cxx:387
 AliHelix.cxx:388
 AliHelix.cxx:389
 AliHelix.cxx:390
 AliHelix.cxx:391
 AliHelix.cxx:392
 AliHelix.cxx:393
 AliHelix.cxx:394
 AliHelix.cxx:395
 AliHelix.cxx:396
 AliHelix.cxx:397
 AliHelix.cxx:398
 AliHelix.cxx:399
 AliHelix.cxx:400
 AliHelix.cxx:401
 AliHelix.cxx:402
 AliHelix.cxx:403
 AliHelix.cxx:404
 AliHelix.cxx:405
 AliHelix.cxx:406
 AliHelix.cxx:407
 AliHelix.cxx:408
 AliHelix.cxx:409
 AliHelix.cxx:410
 AliHelix.cxx:411
 AliHelix.cxx:412
 AliHelix.cxx:413
 AliHelix.cxx:414
 AliHelix.cxx:415
 AliHelix.cxx:416
 AliHelix.cxx:417
 AliHelix.cxx:418
 AliHelix.cxx:419
 AliHelix.cxx:420
 AliHelix.cxx:421
 AliHelix.cxx:422
 AliHelix.cxx:423
 AliHelix.cxx:424
 AliHelix.cxx:425
 AliHelix.cxx:426
 AliHelix.cxx:427
 AliHelix.cxx:428
 AliHelix.cxx:429
 AliHelix.cxx:430
 AliHelix.cxx:431
 AliHelix.cxx:432
 AliHelix.cxx:433
 AliHelix.cxx:434
 AliHelix.cxx:435
 AliHelix.cxx:436
 AliHelix.cxx:437
 AliHelix.cxx:438
 AliHelix.cxx:439
 AliHelix.cxx:440
 AliHelix.cxx:441
 AliHelix.cxx:442
 AliHelix.cxx:443
 AliHelix.cxx:444
 AliHelix.cxx:445
 AliHelix.cxx:446
 AliHelix.cxx:447
 AliHelix.cxx:448
 AliHelix.cxx:449
 AliHelix.cxx:450
 AliHelix.cxx:451
 AliHelix.cxx:452
 AliHelix.cxx:453
 AliHelix.cxx:454
 AliHelix.cxx:455
 AliHelix.cxx:456
 AliHelix.cxx:457
 AliHelix.cxx:458
 AliHelix.cxx:459
 AliHelix.cxx:460
 AliHelix.cxx:461
 AliHelix.cxx:462
 AliHelix.cxx:463
 AliHelix.cxx:464
 AliHelix.cxx:465
 AliHelix.cxx:466
 AliHelix.cxx:467
 AliHelix.cxx:468
 AliHelix.cxx:469
 AliHelix.cxx:470
 AliHelix.cxx:471
 AliHelix.cxx:472
 AliHelix.cxx:473
 AliHelix.cxx:474
 AliHelix.cxx:475
 AliHelix.cxx:476
 AliHelix.cxx:477
 AliHelix.cxx:478
 AliHelix.cxx:479
 AliHelix.cxx:480
 AliHelix.cxx:481
 AliHelix.cxx:482
 AliHelix.cxx:483
 AliHelix.cxx:484
 AliHelix.cxx:485
 AliHelix.cxx:486
 AliHelix.cxx:487
 AliHelix.cxx:488
 AliHelix.cxx:489
 AliHelix.cxx:490
 AliHelix.cxx:491
 AliHelix.cxx:492
 AliHelix.cxx:493
 AliHelix.cxx:494
 AliHelix.cxx:495
 AliHelix.cxx:496
 AliHelix.cxx:497
 AliHelix.cxx:498
 AliHelix.cxx:499
 AliHelix.cxx:500
 AliHelix.cxx:501
 AliHelix.cxx:502
 AliHelix.cxx:503
 AliHelix.cxx:504
 AliHelix.cxx:505
 AliHelix.cxx:506
 AliHelix.cxx:507
 AliHelix.cxx:508
 AliHelix.cxx:509
 AliHelix.cxx:510
 AliHelix.cxx:511
 AliHelix.cxx:512
 AliHelix.cxx:513
 AliHelix.cxx:514
 AliHelix.cxx:515
 AliHelix.cxx:516
 AliHelix.cxx:517
 AliHelix.cxx:518
 AliHelix.cxx:519
 AliHelix.cxx:520
 AliHelix.cxx:521
 AliHelix.cxx:522
 AliHelix.cxx:523
 AliHelix.cxx:524
 AliHelix.cxx:525
 AliHelix.cxx:526
 AliHelix.cxx:527
 AliHelix.cxx:528
 AliHelix.cxx:529
 AliHelix.cxx:530
 AliHelix.cxx:531
 AliHelix.cxx:532
 AliHelix.cxx:533
 AliHelix.cxx:534
 AliHelix.cxx:535
 AliHelix.cxx:536
 AliHelix.cxx:537
 AliHelix.cxx:538
 AliHelix.cxx:539
 AliHelix.cxx:540
 AliHelix.cxx:541
 AliHelix.cxx:542
 AliHelix.cxx:543
 AliHelix.cxx:544
 AliHelix.cxx:545
 AliHelix.cxx:546
 AliHelix.cxx:547
 AliHelix.cxx:548
 AliHelix.cxx:549
 AliHelix.cxx:550
 AliHelix.cxx:551
 AliHelix.cxx:552
 AliHelix.cxx:553
 AliHelix.cxx:554
 AliHelix.cxx:555
 AliHelix.cxx:556
 AliHelix.cxx:557
 AliHelix.cxx:558
 AliHelix.cxx:559
 AliHelix.cxx:560
 AliHelix.cxx:561
 AliHelix.cxx:562
 AliHelix.cxx:563
 AliHelix.cxx:564
 AliHelix.cxx:565
 AliHelix.cxx:566
 AliHelix.cxx:567
 AliHelix.cxx:568
 AliHelix.cxx:569
 AliHelix.cxx:570
 AliHelix.cxx:571
 AliHelix.cxx:572
 AliHelix.cxx:573
 AliHelix.cxx:574
 AliHelix.cxx:575
 AliHelix.cxx:576
 AliHelix.cxx:577
 AliHelix.cxx:578
 AliHelix.cxx:579
 AliHelix.cxx:580
 AliHelix.cxx:581
 AliHelix.cxx:582
 AliHelix.cxx:583
 AliHelix.cxx:584
 AliHelix.cxx:585
 AliHelix.cxx:586
 AliHelix.cxx:587
 AliHelix.cxx:588
 AliHelix.cxx:589
 AliHelix.cxx:590
 AliHelix.cxx:591
 AliHelix.cxx:592
 AliHelix.cxx:593
 AliHelix.cxx:594
 AliHelix.cxx:595
 AliHelix.cxx:596
 AliHelix.cxx:597
 AliHelix.cxx:598
 AliHelix.cxx:599
 AliHelix.cxx:600
 AliHelix.cxx:601
 AliHelix.cxx:602
 AliHelix.cxx:603
 AliHelix.cxx:604
 AliHelix.cxx:605
 AliHelix.cxx:606
 AliHelix.cxx:607
 AliHelix.cxx:608
 AliHelix.cxx:609
 AliHelix.cxx:610
 AliHelix.cxx:611
 AliHelix.cxx:612
 AliHelix.cxx:613
 AliHelix.cxx:614
 AliHelix.cxx:615
 AliHelix.cxx:616
 AliHelix.cxx:617
 AliHelix.cxx:618
 AliHelix.cxx:619
 AliHelix.cxx:620
 AliHelix.cxx:621
 AliHelix.cxx:622
 AliHelix.cxx:623
 AliHelix.cxx:624
 AliHelix.cxx:625
 AliHelix.cxx:626
 AliHelix.cxx:627
 AliHelix.cxx:628
 AliHelix.cxx:629
 AliHelix.cxx:630
 AliHelix.cxx:631
 AliHelix.cxx:632
 AliHelix.cxx:633
 AliHelix.cxx:634
 AliHelix.cxx:635
 AliHelix.cxx:636
 AliHelix.cxx:637
 AliHelix.cxx:638
 AliHelix.cxx:639
 AliHelix.cxx:640
 AliHelix.cxx:641
 AliHelix.cxx:642