Helix.h

/**************************************************************************
 * BASF2 (Belle Analysis Framework 2)                                     *
 * Copyright(C) 2014 - Belle II Collaboration                             *
 *                                                                        *
 * Author: The Belle II Collaboration                                     *
 * Contributors: Oliver Frost                                             *
 *                                                                        *
 * This software is provided "as is" without any warranty.                *
 **************************************************************************/
#pragma once
 
#include <tracking/trackFindingCDC/geometry/PerigeeCircle.h>
#include <tracking/trackFindingCDC/geometry/SZLine.h>
 
#include <tracking/trackFindingCDC/geometry/HelixParameters.h>
#include <tracking/trackFindingCDC/geometry/Vector3D.h>
#include <tracking/trackFindingCDC/geometry/Vector2D.h>
 
#include <iosfwd>
#include <cmath>
 
namespace Belle2 {
  namespace TrackFindingCDC {
 
    class Helix {
 
    public:
      Helix()
        : m_circleXY()
        , m_szLine()
      {
      }
 
      Helix(const PerigeeCircle& circleXY, const SZLine& szLine)
        : m_circleXY(circleXY)
        , m_szLine(szLine)
      {
      }
 
      explicit Helix(const HelixParameters& parameters)
        : m_circleXY(HelixUtil::getPerigeeParameters(parameters))
        , m_szLine(HelixUtil::getSZParameters(parameters))
      {
      }
 
      Helix(double curvature, double phi0, double impact, double tanLambda, double z0)
        : m_circleXY(curvature, phi0, impact)
        , m_szLine(tanLambda, z0)
      {
      }
 
      Helix(double curvature, const Vector2D& phi0Vec, double impact, double tanLambda, double z0)
        : m_circleXY(curvature, phi0Vec, impact)
        , m_szLine(tanLambda, z0)
      {
      }
 
      void invalidate()
      {
        m_circleXY.invalidate();
        m_szLine.invalidate();
      }
 
      bool isInvalid() const
      {
        return circleXY().isInvalid() and szLine().isInvalid();
      }
 
      void reverse()
      {
        m_circleXY.reverse();
        m_szLine.reverse();
      }
 
      Helix reversed() const
      {
        return Helix(circleXY().reversed(), szLine().reversed());
      }
 
    public:
      double arcLength2DToClosest(const Vector3D& point, bool firstPeriod = true) const;
 
      double arcLength2DToXY(const Vector2D& point) const
      {
        return circleXY().arcLengthTo(point);
      }
 
      double arcLength2DToCylindricalR(double cylindricalR) const
      {
        return circleXY().arcLengthToCylindricalR(cylindricalR);
      }
 
      Vector3D closest(const Vector3D& point, bool firstPeriod = true) const
      {
        double arcLength2D = arcLength2DToClosest(point, firstPeriod);
        return atArcLength2D(arcLength2D);
      }
 
      Vector3D closestXY(const Vector2D& pointXY) const
      {
        double arcLength2D = arcLength2DToXY(pointXY);
        return atArcLength2D(arcLength2D);
      }
 
      double distance(const Vector3D& point) const
      {
        return point.distance(closest(point));
      }
 
      double distanceXY(const Vector2D& point) const
      {
        return m_circleXY.distance(point);
      }
 
      double passiveMoveBy(const Vector3D& by)
      {
        // First keep the necessary shift of the perpendicular travel distance to the new perigee
        // point.
        double byS = circleXY().arcLengthTo(by.xy());
        m_circleXY.passiveMoveBy(by.xy());
        Vector2D bySZ(byS, by.z());
        m_szLine.passiveMoveBy(bySZ);
        return byS;
      }
 
      HelixJacobian passiveMoveByJacobian(const Vector3D& by) const;
 
      double shiftPeriod(int nPeriods)
      {
        double arcLength2D = nPeriods * fabs(perimeterXY());
        m_szLine.passiveMoveBy(Vector2D(arcLength2D, 0.0));
        return arcLength2D;
      }
 
      Vector3D atArcLength2D(double s) const
      {
        return Vector3D(circleXY().atArcLength(s), szLine().map(s));
      }
 
      Vector3D atZ(double z) const
      {
        return Vector3D(xyAtZ(z), z);
      }
 
      Vector2D xyAtZ(double z) const
      {
        return Vector2D(circleXY().atArcLength(szLine().inverseMap(z)));
      }
 
      double minimalCylindricalR() const
      {
        return circleXY().minimalCylindricalR();
      }
 
      double maximalCylindricalR() const
      {
        return circleXY().maximalCylindricalR();
      }
 
      double curvatureXY() const
      {
        return circleXY().curvature();
      }
 
      double impactXY() const
      {
        return circleXY().impact();
      }
 
      double omega() const
      {
        return circleXY().omega();
      }
 
      double d0() const
      {
        return circleXY().d0();
      }
 
      Vector2D perigeeXY() const
      {
        return circleXY().perigee();
      }
 
      Vector3D perigee() const
      {
        return Vector3D(perigeeXY(), z0());
      }
 
      double tanLambda() const
      {
        return m_szLine.tanLambda();
      }
 
      double cotTheta() const
      {
        return tanLambda();
      }
 
      double z0() const
      {
        return m_szLine.z0();
      }
 
      double zPeriod() const
      {
        return tanLambda() * arcLength2DPeriod();
      }
 
      double arcLength2DPeriod() const
      {
        return circleXY().arcLengthPeriod();
      }
 
      double perimeterXY() const
      {
        return circleXY().perimeter();
      }
 
      double radiusXY() const
      {
        return circleXY().radius();
      }
 
      Vector2D centerXY() const
      {
        return circleXY().center();
      }
 
      Vector3D tangential() const
      {
        return Vector3D(phi0Vec(), tanLambda()).unit();
      }
 
      const Vector2D& phi0Vec() const
      {
        return circleXY().phi0Vec();
      }
 
      double phi0() const
      {
        return circleXY().phi0();
      }
 
      HelixParameters helixParameters() const
      {
        HelixParameters result;
        using namespace NHelixParameterIndices;
        result(c_Curv) = curvatureXY();
        result(c_Phi0) = phi0();
        result(c_I) = impactXY();
        result(c_TanL) = tanLambda();
        result(c_Z0) = z0();
        return result;
      }
 
      const PerigeeCircle& circleXY() const
      {
        return m_circleXY;
      }
 
      const SZLine& szLine() const
      {
        return m_szLine;
      }
 
    private:
      PerigeeCircle m_circleXY;
 
      SZLine m_szLine;
 
    };
 
    std::ostream& operator<<(std::ostream& output, const Helix& helix);
  }
}