Helix.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
#pragma once
 
#include <tracking/trackingUtilities/geometry/PerigeeCircle.h>
#include <tracking/trackingUtilities/geometry/SZLine.h>
 
#include <tracking/trackingUtilities/geometry/HelixParameters.h>
#include <tracking/trackingUtilities/geometry/VectorUtil.h>
 
#include <Math/Vector3D.h>
#include <Math/Vector2D.h>
 
#include <iosfwd>
#include <cmath>
 
namespace Belle2 {
  namespace TrackingUtilities {

    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 ROOT::Math::XYVector& 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 ROOT::Math::XYZVector& point, bool firstPeriod = true) const;

      double arcLength2DToXY(const ROOT::Math::XYVector& point) const
      {
        return circleXY().arcLengthTo(point);
      }

      double arcLength2DToCylindricalR(double cylindricalR) const
      {
        return circleXY().arcLengthToCylindricalR(cylindricalR);
      }

      ROOT::Math::XYZVector closest(const ROOT::Math::XYZVector& point, bool firstPeriod = true) const
      {
        double arcLength2D = arcLength2DToClosest(point, firstPeriod);
        return atArcLength2D(arcLength2D);
      }

      ROOT::Math::XYZVector closestXY(const ROOT::Math::XYVector& pointXY) const
      {
        double arcLength2D = arcLength2DToXY(pointXY);
        return atArcLength2D(arcLength2D);
      }

      double distance(const ROOT::Math::XYZVector& point) const
      {
        return VectorUtil::Distance(point, closest(point));
      }

      double distanceXY(const ROOT::Math::XYVector& point) const
      {
        return m_circleXY.distance(point);
      }

      double passiveMoveBy(const ROOT::Math::XYZVector& by)
      {
        // First keep the necessary shift of the perpendicular travel distance to the new perigee
        // point.
        double byS = circleXY().arcLengthTo(VectorUtil::getXYVector(by));
        m_circleXY.passiveMoveBy(VectorUtil::getXYVector(by));
        ROOT::Math::XYVector bySZ(byS, by.z());
        m_szLine.passiveMoveBy(bySZ);
        return byS;
      }

      HelixJacobian passiveMoveByJacobian(const ROOT::Math::XYZVector& by) const;

      double shiftPeriod(int nPeriods)
      {
        double arcLength2D = nPeriods * fabs(perimeterXY());
        m_szLine.passiveMoveBy(ROOT::Math::XYVector(arcLength2D, 0.0));
        return arcLength2D;
      }

      ROOT::Math::XYZVector atArcLength2D(double s) const
      {
        const auto& tmp = circleXY().atArcLength(s);
        return ROOT::Math::XYZVector(tmp.X(), tmp.Y(), szLine().map(s));
      }

      ROOT::Math::XYZVector atZ(double z) const
      {
        const auto& tmp = xyAtZ(z);
        return ROOT::Math::XYZVector(tmp.X(), tmp.Y(), z);
      }

      ROOT::Math::XYVector xyAtZ(double z) const
      {
        return ROOT::Math::XYVector(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();
      }

      ROOT::Math::XYVector perigeeXY() const
      {
        return circleXY().perigee();
      }

      ROOT::Math::XYZVector perigee() const
      {
        const auto& tmp = perigeeXY();
        return ROOT::Math::XYZVector(tmp.X(), tmp.Y(), 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();
      }

      ROOT::Math::XYVector centerXY() const
      {
        return circleXY().center();
      }

      ROOT::Math::XYZVector tangential() const
      {
        const auto& tmp = phi0Vec();
        return VectorUtil::unit(ROOT::Math::XYZVector(tmp.X(), tmp.Y(), tanLambda()));
      }

      const ROOT::Math::XYVector& 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);
  }
}