#include "CbmSttGeomHelix.h"
using namespace std;
#define fieldStrength 20. // kiloGauss
#define conversionFactor 0.00029979 // c * 10e-12


void CbmSttGeomHelix::Copy(CbmSttGeomHelix const &other)
{
    fAlpha = other.fAlpha;
    fRadius = other.fRadius;
    fDistanceZ = other.fDistanceZ;
    fDistance = other.fDistance;
    fPhi = other.fPhi;
    fCharge = other.fCharge;
    fTolerance = other.fTolerance;
    fMaxIterations = other.fMaxIterations;
}

CbmSttGeomHelix::CbmSttGeomHelix(HepVector const babarVector)
{
    ConstructHelixFromBabarVector(babarVector);
}

void CbmSttGeomHelix::ConstructHelixFromBabarVector(HepVector const babarVector)
{
    // some sign difference might be left here due to different angle definitions
    // look at this later
    fAlpha = babarVector[4] / babarVector[2];
    fRadius = 1 / babarVector[2];
    fDistanceZ = babarVector[3] - ((babarVector[1] * babarVector[4]) / babarVector[2]);
    fDistance = babarVector[0];
    fPhi = babarVector[1];

    fCharge = 0.;
    fTolerance = 0.000001;
    fMaxIterations = 500;
}

void CbmSttGeomHelix::Destroy()
{
}

CbmSttGeomHelix::CbmSttGeomHelix()
{
    fAlpha = 0.;
    fRadius = 0.;
    fDistanceZ = 0.;
    fDistance = 0.;
    fPhi = 0.;
    fCharge = 0.;

    fTolerance = 0.000001;
    fMaxIterations = 500;
}    

CbmSttGeomHelix::CbmSttGeomHelix(CbmSttTrack *track)
{
    fDistance = track->GetParamLast()->GetX();
    fPhi = track->GetParamLast()->GetY();
    fDistanceZ = track->GetParamLast()->GetZ();
    fRadius = track->GetParamLast()->GetTx();
    fAlpha = track->GetParamLast()->GetTy();
    fCharge = track->GetParamLast()->GetQp();

    fTolerance = 0.000001;
    fMaxIterations = 500;
}

CbmSttGeomHelix::CbmSttGeomHelix(CbmMCTrack *track)
{
    fCharge = TDatabasePDG::Instance()->GetParticle(track->GetPdgCode())->Charge() / 3.;
    fTolerance = 0.000001;
    fMaxIterations = 500;

    ConstructHelixFromMomenta(track->GetMomentum().X(), track->GetMomentum().Y(), track->GetMomentum().Z(), 
			      track->GetStartVertex().X(), track->GetStartVertex().Y(), track->GetStartVertex().Z());   

    CbmSttGeomPoint
	vertex(track->GetStartVertex().X(), track->GetStartVertex().Y(), track->GetStartVertex().Z());

    TVector3
	momentum = GetMomentumAt(vertex);
}

CbmSttGeomHelix::CbmSttGeomHelix(Double_t px, Double_t py, Double_t pz, Double_t vx, Double_t vy, Double_t vz, Double_t charge)
{
    fCharge = charge;
    fTolerance = 0.000001;
    fMaxIterations = 500;

    ConstructHelixFromMomenta(px, py, pz, vx, vy, vz);
}

void CbmSttGeomHelix::ConstructHelixFromMomenta(Double_t px, Double_t py, Double_t pz, Double_t vx, Double_t vy, Double_t vz)
{
    Double_t
        mcTrackPt = sqrt(px * px + py * py),
        mcTrackTheta = atan(mcTrackPt / pz);
    
    if (pz < 0.)
    {
        if (mcTrackPt < 0.)
            mcTrackTheta -= dPi;
        else
            mcTrackTheta += dPi;
    }
    
    fRadius = sqrt(px * px + py * py) / (sqrt(fCharge * fCharge) * conversionFactor * fieldStrength);
    
    Double_t
        phiStart = atan(py / px),
        xStart = vx,
        yStart = vy;

    if (px < 0.)
    {
        if (py < 0.)
            phiStart -= dPi;
        else
            phiStart += dPi;
    }
    
    Double_t
        xCircleCenter = xStart + (GetChargeSign() * fRadius * sin(phiStart)),
        yCircleCenter = yStart - (GetChargeSign() * fRadius * cos(phiStart));
    
    fDistance = sqrt(xCircleCenter * xCircleCenter + yCircleCenter * yCircleCenter) - fRadius;
    fPhi = atan(yCircleCenter / xCircleCenter);
    
    if (xCircleCenter < 0.)
    {
        if (yCircleCenter < 0.)
            fPhi -= dPi;
        else
            fPhi += dPi;
    }

    Double_t
        mcTmpXcenter = (fRadius + fDistance) * cos(fPhi),
        mcTmpYcenter = (fRadius + fDistance) * sin(fPhi),
        mcTmpPhi = atan((vy - mcTmpYcenter) / (vx - mcTmpXcenter));
    
    if (vx - mcTmpXcenter < 0.)
    {
        if (vy - mcTmpYcenter < 0.)
            mcTmpPhi -= dPi;
        else
            mcTmpPhi += dPi;
    }

    fAlpha = (-GetChargeSign() * fRadius) / tan(mcTrackTheta);
    fDistanceZ = vz - fAlpha * mcTmpPhi;
}

Double_t CbmSttGeomHelix::GetPt() const
{
    Double_t
	Pt = (sqrt(fCharge * fCharge) * conversionFactor * fieldStrength) * fRadius;

    return Pt;
}

Double_t CbmSttGeomHelix::GetP() const
{
    if (GetChargeSign() == 0)
    {
	cerr << "<I> CbmSttGeomHelix: No charge determined, cannot calculate theta of track" << endl;
	return 0.;
    }
    
    Double_t
	theta = GetTheta(),
	Pt = GetPt();

    Double_t
	P =  Pt / sin(theta);
    
    return P;
}

Double_t CbmSttGeomHelix::GetTheta() const
{
    if (GetChargeSign() == 0)
    {
	cerr << "<I> CbmSttGeomHelix: No charge determined, cannot calculate theta of track" << endl;
    }

    Double_t
	theta = atan(fRadius / fAlpha);
        
    if (fAlpha < 0.)
    {
        if (fRadius < 0.)
            theta -= dPi;
        else
            theta += dPi;
    }
   
    return -1 * GetChargeSign() * theta;
}

Double_t CbmSttGeomHelix::GetPhiAt(CbmSttGeomPoint thePoint) 
{
    // find closest point to thePoint on helix
    CbmSttGeomPoint
	closestPoint = ClosestPointTo(thePoint);
    
    Double_t
	pointX = closestPoint.GetX() - GetAxisX(),
	pointY = closestPoint.GetY() - GetAxisY();
    
    // calculate phi angle of particle at evaluation point
    Double_t
	phi = atan(pointY / pointX);
 
    if (pointX < 0.)
    {
	if (pointY < 0.)
	    phi -= dPi;
	else
	    phi += dPi;
    }

    // momentum is perpendicular to the line from the helix 
    // center to the evaluation point, exact direction 
    // depends on charge
    phi -= GetChargeSign() * (dPi / 2.);

    return phi;
}

Double_t CbmSttGeomHelix::GetChargeSign() const
{
    Double_t
	retval;

    if (fCharge > 0.)
	retval = 1.;
    else if (fCharge < 0.)
	retval = -1.;
    else 
	retval = 0.;

    return retval;
}

TVector3 CbmSttGeomHelix::GetMomentumAt(CbmSttGeomPoint thePoint) 
{
    cout << "thePoint: " << thePoint.GetX() << "   " << thePoint.GetY() << "   " << thePoint.GetZ() << endl;

    // calculate theta angle of particle  
    Double_t
	phi = GetPhiAt(thePoint),
        theta = GetTheta(),
	P = GetP(),
	Px = P * cos(phi) * sin(theta),
	Py = P * sin(phi) * sin(theta),
	Pz = P * cos(theta);

    // return momentum vector
    TVector3
	retval(Px, Py, Pz);

    return retval;
}


CbmSttGeomHelix::CbmSttGeomHelix(Double_t distance, Double_t phi, Double_t radius, Double_t zdistance, Double_t alpha, Double_t charge)
{
    fAlpha = alpha;
    fRadius = radius;
    fDistanceZ = zdistance;
    fDistance = distance;
    fPhi = phi;
    fCharge = charge;

    fTolerance = 0.000001;
    fMaxIterations = 500;
}    

CbmSttGeomHelix::~CbmSttGeomHelix()
{
    Destroy();
}    

CbmSttGeomHelix::CbmSttGeomHelix(CbmSttGeomHelix const &other)
{
    if (this != &other)
    {
	Destroy();
	Copy(other);
    }
}

void CbmSttGeomHelix::operator=(CbmSttGeomHelix const &other)
{
    Copy(other);
}

Double_t CbmSttGeomHelix::DistanceTo(CbmSttGeomPoint myPoint)
{
    CbmSttGeomPoint
	closestPoint = ClosestPointTo(myPoint);
    
    return closestPoint.DistanceTo(myPoint); 
}

CbmSttGeomPoint CbmSttGeomHelix::ClosestPointTo(CbmSttGeomPoint myPoint)
{
    SetupMinimization(myPoint);
    
    Double_t
	phi = GetPhiEstimate(myPoint);

    if (!Minimize(phi))
    {
	cout << "no dice" << endl;
    }

    CbmSttGeomPoint
	closestPoint = At(phi);
    
    return closestPoint; 
}

CbmSttGeomPoint CbmSttGeomHelix::At(Double_t phi) const
{
    CbmSttGeomPoint
	retval(fRadius * cos(phi) + GetAxisX(), fRadius * sin(phi) + GetAxisY(), fAlpha * phi + fDistanceZ);

    return retval;
}

void CbmSttGeomHelix::SetupMinimization(CbmSttGeomPoint myPoint)
{
    Double_t
	xtrans = (fDistance + fRadius) * cos(fPhi),
	ytrans = (fDistance + fRadius) * sin(fPhi),
	ztrans = fDistanceZ;
    
    // transform point into helix frame
    myPoint.Transform(-xtrans, -ytrans, -ztrans);
    
    // assign the private members
    m_a0 = fAlpha * -1 * myPoint.GetZ();
    m_a1 = fAlpha * fAlpha;
    m_a2 = fRadius * myPoint.GetX();
    m_a3 = fRadius * -1 * myPoint.GetY();
    m_a4 = 0.;
    m_a5 = 0.;
    m_a6 = 0.;
    m_a7 = 0.;
}

Double_t CbmSttGeomHelix::GetAxisX() const
{
    return (fDistance + fRadius) * cos(fPhi);
}

Double_t CbmSttGeomHelix::GetAxisY() const
{
    return (fDistance + fRadius) * sin(fPhi);
}

void CbmSttGeomHelix::SetupMinimization(CbmSttGeomLine myLine)
{
    Double_t
	xtrans = GetAxisX(),
	ytrans = GetAxisY(),
	ztrans = fDistanceZ;
    
    // transform line into helix frame
    myLine.Transform(-xtrans, -ytrans, -ztrans);

    // lamba(t) = l0 + t * l  .... already in HelixFrame
    TVector3
	zaxis(0., 0., -1.),
	l0(myLine.GetX1(), myLine.GetY1(), myLine.GetZ1()),
	l(myLine.GetX2() - myLine.GetX1(), 
	  myLine.GetY2() - myLine.GetY1(), 
	  myLine.GetZ2() - myLine.GetZ1());

    Double_t
	inverseL2 = 1. / l.Mag2(),
	Ltilde2 = l0.Dot(l),
	Lratio = inverseL2 * Ltilde2; 
    
    // assign the private members
    m_a0 = fAlpha * (-1 * l0.z() + Lratio * l.z());
    m_a1 = fAlpha * fAlpha * (1. - l.z() * l.z() * inverseL2);
    m_a2 = fRadius * (l0.x() - l.x() * Lratio - inverseL2 * l.y() * l.z() * fAlpha);
    m_a3 = fRadius * (-1 * l0.y() + l.y() * Lratio - inverseL2 * l.x() * l.z() * fAlpha);
    m_a4 = 0.5 * inverseL2 * fRadius * fRadius * (l.x() * l.x() - l.y() * l.y());
    m_a5 = -1 * inverseL2 * fRadius * fRadius * l.x() * l.y();
    m_a6 = fAlpha * fRadius * inverseL2 * l.z() * l.x();
    m_a7 = -1 * m_a6 / l.x() * l.y();
}

Double_t CbmSttGeomHelix::DistanceTo(CbmSttGeomLine myLine)
{
    if (myLine.IsParallelToZaxis())
    {
	return sqrt((myLine.GetX1() - GetAxisX()) * (myLine.GetX1() - GetAxisX()) + 
		    (myLine.GetY1() - GetAxisY()) * (myLine.GetY1() - GetAxisY())) - fRadius;
    }

    CbmSttGeomPoint
	closestPoint = ClosestPointTo(myLine);
    
    return closestPoint.DistanceTo(myLine); 
}

CbmSttGeomPoint CbmSttGeomHelix::ClosestPointTo(CbmSttGeomLine myLine)
{
    SetupMinimization(myLine);
    
    Double_t
	phi = GetPhiEstimate(myLine);

    if (!myLine.IsParallelToZaxis())
    {
	if (!Minimize(phi))
	{
	    cout << "no dice" << endl;
	}
    }

    CbmSttGeomPoint
	closestPoint = At(phi);

    return closestPoint; 
}

Bool_t CbmSttGeomHelix::Minimize(Double_t &phi)
{
    Bool_t 
	solfound = false; 
     
    for (Int_t istep = 0; istep < fMaxIterations; istep++)
     {
         phi = phi - f(phi) / fder(phi);

         if (fabs(f(phi)) < fTolerance)
	 {
	     solfound = true;
	     break;
         }
     }
     
     return solfound;
}

Double_t CbmSttGeomHelix::GetPhiEstimate(CbmSttGeomPoint point)
{
    Double_t
	xtrans = GetAxisX(),
	ytrans = GetAxisY(),
	ztrans = fDistanceZ;
    
    // transform point into helix frame
    point.Transform(-xtrans, -ytrans, -ztrans);

    // intersect line with helix in x-y plane
    CbmSttGeomLine
	myLine(point.GetX(), point.GetY(), 0., 0., 0., 0.);

    CbmSttGeomCircle
	myCircle(0., 0., fRadius);

    CbmSttGeomPoint
	intersectionPoint1,
	intersectionPoint2;
    
    myCircle.IntersectWith(myLine, intersectionPoint1, intersectionPoint2);
    
    Double_t
	intersectX1 = intersectionPoint1.GetX(),
	intersectY1 = intersectionPoint1.GetY(),
	intersectX2 = intersectionPoint2.GetX(),
	intersectY2 = intersectionPoint2.GetY();
    
    // calculate phi angle
    Double_t
	intersectAngle1 = atan(intersectY1 / intersectX1),
	intersectAngle2 = atan(intersectY2 / intersectX2);
 
    if (intersectX1 < 0.)
      {
	if (intersectY1 < 0.)
	  intersectAngle1 -= dPi;
	else
	  intersectAngle1 += dPi;
      }

    if (intersectX2 < 0.)
      {
	if (intersectY2 < 0.)
	  intersectAngle2 -= dPi;
	else
	  intersectAngle2 += dPi;
      }

    // obtain closest winding
    Int_t
	revNumber1 = ((Int_t) (((point.GetZ() - fAlpha * intersectAngle1) / (2 * dPi * fAlpha)) + 0.5)),
	revNumber2 = ((Int_t) (((point.GetZ() - fAlpha * intersectAngle2) / (2 * dPi * fAlpha)) + 0.5));
	
    Double_t
	helixAngle1 = intersectAngle1 + 2 * dPi * revNumber1,
	helixAngle2 = intersectAngle2 + 2 * dPi * revNumber2;

    CbmSttGeomPoint
	closestPoint1 = At(helixAngle1),
	closestPoint2 = At(helixAngle2);

    Double_t
	retval;

    if (closestPoint1.DistanceTo(point) < closestPoint2.DistanceTo(point))
    {
	retval = helixAngle1;
    }
    else
    {
	retval = helixAngle2;
    }

    // TODO: check one turn above and one below

    return retval;
}

void CbmSttGeomHelix::Transform(Double_t shiftX, Double_t shiftY, Double_t shiftZ)
{
    Double_t
	xcenter = (fDistance + fRadius) * cos(fPhi),
	ycenter = (fDistance + fRadius) * sin(fPhi);

    xcenter += shiftX;
    ycenter += shiftY;

    fDistance = sqrt(xcenter * xcenter + ycenter * ycenter) - fRadius;
    fPhi = atan(ycenter / xcenter);

    if (xcenter < 0.)
    {
	if (ycenter < 0.)
	    fPhi -= dPi;
	else
	    fPhi += dPi;
    }

    fDistanceZ += shiftZ;
}

Double_t CbmSttGeomHelix::GetPhiEstimate(CbmSttGeomHelix helix)
{
  Double_t
    xtrans = GetAxisX(),
    ytrans = GetAxisY(),
    ztrans = fDistanceZ;
  
  // transform other helix into this helix frame
  helix.Transform(-xtrans, -ytrans, -ztrans);
  
  // find the intersection of the two helices in 2D
  CbmSttGeomCircle
    myOtherCircle(helix.GetAxisX(), helix.GetAxisY(), helix.GetRadius()),
    myCircle(0., 0., fRadius);
  
  CbmSttGeomPoint
    intersection1,
    intersection2;

  Int_t
      intersections = myOtherCircle.IntersectWith(myCircle, intersection1, intersection2);

  Double_t
      intersectAngle,
      intersectAngle1 = atan(intersection1.GetY() / intersection1.GetX()),
      intersectAngle2 = atan(intersection2.GetY() / intersection2.GetX());
  
  if (intersection1.GetX() < 0.)
    {
      if (intersection1.GetY() < 0.)
	intersectAngle1 -= dPi;
      else
	intersectAngle1 += dPi;
    }

  if (intersection2.GetX() < 0.)
    {
      if (intersection2.GetY() < 0.)
	intersectAngle2 -= dPi;
      else
	intersectAngle2 += dPi;
    }

  // no intersections
  if (intersections == 0)
  {
      // the centers of the helices in the x-y plane
      CbmSttGeomPoint
	  thisCenter(0., 0., 0.),
	  otherCenter(helix.GetAxisX(), helix.GetAxisY(), 0.);

      // form the line between the two helix centers in the x-y plane
      CbmSttGeomLine
	  CenterCenterLine(thisCenter, otherCenter);
      
      // get intersections between the other helix and this line
      myOtherCircle.IntersectWith(CenterCenterLine, intersection1, intersection2);
      
      // select intersection which lies between the centers
      if (intersection1.IsBetween(thisCenter, otherCenter))
	  intersectAngle = intersectAngle1;
      else
	  intersectAngle = intersectAngle2;
  }
  else if (intersections == 1)
  {
      // in this case both answers must be the same, we choose the first
      intersectAngle = intersectAngle1;
  }
  else
  {
      // construct the two points on this helix 
      CbmSttGeomPoint
	  helixPoint1 = At(intersectAngle1),
	  helixPoint2 = At(intersectAngle2);

      // and calculate the distance to the other helix
      Double_t
	  distance1 = helix.DistanceTo(helixPoint1), 
	  distance2 = helix.DistanceTo(helixPoint2); 

      if (distance1 < distance2)
 	  intersectAngle = intersectAngle1;
      else
	  intersectAngle = intersectAngle2;
  }

  // since helices are parallel we can't obtain the winding, so we'll have to take the first one
  return intersectAngle;
}

void CbmSttGeomHelix::Draw(Double_t size, Int_t color) const
{
    Int_t
	steps = 300;

    for (Int_t teller = 0; teller < steps; teller++)
    {
	CbmSttGeomPoint
	    tmpPoint = At(teller * (6 * dPi / steps) - 3 * dPi);

	tmpPoint.Draw(size, color);
    }
}

void CbmSttGeomHelix::DrawFromVertex(CbmSttGeomPoint vertex, Double_t size, Int_t color) const
{
    Int_t
	steps = 300;

    for (Int_t teller = 0; teller < steps; teller++)
    {
	CbmSttGeomPoint
	    tmpPoint = At(teller * (6 * dPi / steps) - 3 * dPi);

	if (((GetTheta() < dPi / 2.) && (tmpPoint.GetZ() > vertex.GetZ())) ||
	    ((GetTheta() > dPi / 2.) && (tmpPoint.GetZ() < vertex.GetZ()))) 
	{
	    tmpPoint.Draw(size, color);
	}
    }
}

Double_t CbmSttGeomHelix::GetPhiEstimate(CbmSttGeomLine line)
{   
  Double_t
    xtrans = GetAxisX(),
    ytrans = GetAxisY(),
    ztrans = fDistanceZ;
  
  // transform line into helix frame
  line.Transform(-xtrans, -ytrans, -ztrans);
  
  // find the intersection of the line and track cylinder in 2D
  CbmSttGeomPoint
    startPoint(line.GetX1(), line.GetY1(), 0.),
    endPoint(line.GetX2(), line.GetY2(), 0.);
  
  CbmSttGeomLine
    myLine(startPoint, endPoint);
  
  CbmSttGeomCircle
    myCircle(0., 0., fRadius);
  
  CbmSttGeomPoint
    intersectionPoint1,
    intersectionPoint2;

  myCircle.IntersectWith(myLine, intersectionPoint1, intersectionPoint2);

  // get z-coordinate
  Double_t
      interSectDistance1 = startPoint.DistanceTo(intersectionPoint1),
      interSectDistance2 = startPoint.DistanceTo(intersectionPoint2),
      totalDistance = startPoint.DistanceTo(endPoint),
      scaleFactor1 = interSectDistance1 / totalDistance,
      scaleFactor2 = interSectDistance2 / totalDistance,
      intersectX1 = line.GetX1() + scaleFactor1 * (line.GetX2() - line.GetX1()),
      intersectY1 = line.GetY1() + scaleFactor1 * (line.GetY2() - line.GetY1()),
      intersectZ1 = line.GetZ1() + scaleFactor1 * (line.GetZ2() - line.GetZ1()),
      intersectX2 = line.GetX1() + scaleFactor2 * (line.GetX2() - line.GetX1()),
      intersectY2 = line.GetY1() + scaleFactor2 * (line.GetY2() - line.GetY1()),
      intersectZ2 = line.GetZ1() + scaleFactor2 * (line.GetZ2() - line.GetZ1());
  
  // calculate phi angle
  Double_t
      intersectAngle1 = atan(intersectionPoint1.GetY() / intersectionPoint1.GetX()),
      intersectAngle2 = atan(intersectionPoint2.GetY() / intersectionPoint2.GetX()); 
  
  if (intersectX1 < 0.)
    {
      if (intersectY1 < 0.)
	intersectAngle1 -= dPi;
      else
	intersectAngle1 += dPi;
    }
  if (intersectX2 < 0.)
    {
      if (intersectY2 < 0.)
	intersectAngle2 -= dPi;
      else
	intersectAngle2 += dPi;
    }
  
  // obtain closest winding
  Int_t
      revNumber1 = ((Int_t) (((intersectZ1 - fAlpha * intersectAngle1) / (2 * dPi * fAlpha)) + 0.5)),
      revNumber2 = ((Int_t) (((intersectZ2 - fAlpha * intersectAngle2) / (2 * dPi * fAlpha)) + 0.5));

  Double_t
      helixAngle1 = intersectAngle1 + 2 * dPi * revNumber1,
      helixAngle2 = intersectAngle2 + 2 * dPi * revNumber2;

  CbmSttGeomPoint
      closestPoint1 = At(helixAngle1),
      closestPoint2 = At(helixAngle2);

  Double_t
      retval;
  
  if (closestPoint1.DistanceTo(line) < closestPoint2.DistanceTo(line))
  {
      retval = helixAngle1;
  }
  else
  {
      retval = helixAngle2;
  }

  // TODO: check one turn above and one below
  
  return retval;
}

Double_t CbmSttGeomHelix::f(Double_t phi)
{
    return (m_a0 + m_a1 * phi + m_a2 * sin(phi) + 
	    m_a3 * cos(phi) + m_a4 * sin(2. * phi) +
	    m_a5 * cos(2. * phi) + m_a6 * phi * sin(phi) +
	    m_a7 * phi * cos(phi)); 
}

Double_t CbmSttGeomHelix::fder(Double_t phi)
{
    return (m_a1 + (m_a2 + m_a7) * cos(phi) + (-1 * m_a3 + m_a6) * sin(phi) +
	    2. * m_a4 * cos(2. * phi) - 2. * m_a5 * sin(2. * phi) + 
	    m_a6 * phi * cos(phi) - m_a7 * phi * sin(phi));
}

Double_t CbmSttGeomHelix::DistanceTo(CbmSttGeomHelix myHelix)
{
/*
  double s1, s2;
  
  //
  //  First step: get dca in the xy-plane as start value
  //

  double dx = myHelix.xcenter() - xcenter();
  double dy = h.ycenter() - ycenter();
  double dd = ::sqrt(dx*dx + dy*dy);
  double r1 = 1/curvature();
  double r2 = 1/h.curvature();
  
  double cosAlpha = (r1*r1 + dd*dd - r2*r2)/(2*r1*dd);
  
  double s;
  double x, y;

  if (fabs(cosAlpha) < 1) 
    {           
      // two solutions
      double 
	sinAlpha = sin(acos(cosAlpha));

      x = xcenter() + r1*(cosAlpha*dx - sinAlpha*dy)/dd;
      y = ycenter() + r1*(sinAlpha*dx + cosAlpha*dy)/dd;
      s = pathLength(x, y);

      x = xcenter() + r1*(cosAlpha*dx + sinAlpha*dy)/dd;
      y = ycenter() + r1*(cosAlpha*dy - sinAlpha*dx)/dd;

      double 
	a = pathLength(x, y);
      
      if (h.distance(at(a)) < h.distance(at(s))) 
	s = a;
    }
  else 
    {                              // no intersection (or exactly one)
      x = xcenter() + r1 * dx / dd;
      y = ycenter() + r1 * dy / dd;
      s = pathLength(x, y);
    }
  
  //
  //   Second step: scan in decreasing intervals around seed 's'
  // 
  const double MinStepSize = 10*micrometer;
  const double MinRange    = 10*centimeter;    
  double dmin              = h.distance(at(s));
  double range             = max(2*dmin, MinRange);
  double ds                = range/10;
  double slast=-999999, ss, d;
  s1 = s - range/2.;
  s2 = s + range/2.;
  
  while (ds > MinStepSize) 
    {
      for (ss=s1; ss<s2+ds; ss+=ds) 
	{
	  d = h.distance(at(ss));
	  if (d < dmin) 
	    {
	      dmin = d;
	      s = ss;
	    }
	  slast = ss;
	}

      //
      //  In the rare cases where the minimum is at the
      //  the border of the current range we shift the range
      //  and start all over, i.e we do not decrease 'ds'.
      //  Else we decrease the search intervall around the
      //  current minimum and redo the scan in smaller steps.
      //

      if (s == s1) 
	{
	  d = 0.8*(s2-s1);
	  s1 -= d;
	  s2 -= d;
	}
      else if (s == slast) 
	{
	  d = 0.8*(s2-s1);
	  s1 += d;
	  s2 += d;
	}
      else 
	{           
	  s1 = s-ds;
	  s2 = s+ds;
	  ds /= 10;
	}
    }
  return pair<double, double>(s, h.pathLength(at(s)));
*/
    return 0.;
}

				 
ClassImp(CbmSttGeomHelix)