GeneralizedCircle.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/Line2D.h>
#include <tracking/trackingUtilities/geometry/VectorUtil.h>
 
#include <tracking/trackingUtilities/numerics/EForwardBackward.h>
#include <tracking/trackingUtilities/numerics/ERightLeft.h>
#include <tracking/trackingUtilities/numerics/ERotation.h>
#include <tracking/trackingUtilities/numerics/ESign.h>
 
#include <Math/Vector2D.h>
 
#include <utility>
#include <iosfwd>
#include <cmath>
 
namespace Belle2 {
 
  namespace TrackingUtilities {
    class Circle2D;

    class GeneralizedCircle {
 
    public:
      GeneralizedCircle() = default;

      GeneralizedCircle(double n0, double n1, double n2, double n3 = 0);

      GeneralizedCircle(double n0, const ROOT::Math::XYVector& n12, double n3 = 0);

      explicit GeneralizedCircle(const Line2D& n012);

      explicit GeneralizedCircle(const Circle2D& circle);

      static GeneralizedCircle
      fromCenterAndRadius(const ROOT::Math::XYVector& center,
                          double absRadius,
                          ERotation orientation = ERotation::c_CounterClockwise);

      static GeneralizedCircle
      fromPerigeeParameters(double curvature, const ROOT::Math::XYVector& tangential, double impact);

      static GeneralizedCircle
      fromPerigeeParameters(double curvature, double tangentialPhi, double impact);
 
    protected:
      void setN0(const double n0)
      {
        m_n0 = n0;
      }

      void setN1(const double n1)
      {
        m_n12.SetX(n1);
      }

      void setN2(const double n2)
      {
        m_n12.SetY(n2);
      }

      void setN12(const double n1, const double n2)
      {
        m_n12.SetXY(n1, n2);
      }

      void setN12(const ROOT::Math::XYVector& n12)
      {
        m_n12.SetXY(n12.X(), n12.Y());
      }

      void setN3(const double n3)
      {
        m_n3 = n3;
      }
 
    public:
      void setCenterAndRadius(const ROOT::Math::XYVector& center,
                              double absRadius,
                              ERotation orientation = ERotation::c_CounterClockwise);

      void setPerigeeParameters(double curvature, const ROOT::Math::XYVector& tangential, double impact);

      void
      setPerigeeParameters(const double curvature, const double tangentialPhi, const double impact)
      {
        setPerigeeParameters(curvature, VectorUtil::Phi(tangentialPhi), impact);
      }

      void setN(const double n0, const double n1, const double n2, const double n3 = 0.0)
      {
        setN0(n0);
        setN12(n1, n2);
        setN3(n3);
        normalize();
      }

      void setN(const double n0, const ROOT::Math::XYVector& n12, const double n3 = 0.0)
      {
        setN0(n0);
        setN12(n12);
        setN3(n3);
        normalize();
      }

      void setN(const Line2D& n012)
      {
        setN(n012.n0(), n012.n12());
      }

      void setN(const GeneralizedCircle& n0123)
      {
        setN(n0123.n0(), n0123.n12(), n0123.n3());
      }

      void invalidate()
      {
        setN(0.0, 0.0, 0.0, 0.0);
      }

      void reverse()
      {
        scaleN(-1);
      }

      void conformalTransform()
      {
        std::swap(m_n0, m_n3);
        reverse(); // Correct orientation
      }

      void passiveMoveBy(const ROOT::Math::XYVector& by)
      {
        setN(fastDistance(by), gradient(by), n3());
      }
 
    protected:
      void normalize()
      {
        double normalization_squared = normalizationSquared();
        if (normalization_squared > 0) scaleN(1.0 / std::sqrt(normalization_squared));
      }
 
    private:
      void scaleN(const double factor)
      {
        m_n0 *= factor;
        m_n12 *= factor;
        m_n3 *= factor;
      }
 
    public:
      double n0() const
      {
        return m_n0;
      }

      double n1() const
      {
        return m_n12.x();
      }

      double n2() const
      {
        return m_n12.y();
      }

      const ROOT::Math::XYVector& n12() const
      {
        return m_n12;
      }

      double n3() const
      {
        return m_n3;
      }
 
    public:
      bool isInvalid() const
      {
        return n0() == 0 and VectorUtil::isNull(n12()) and n3() == 0;
      }

      bool isValid() const
      {
        return not isInvalid();
      }

      double normalizationSquared() const
      {
        return n12().Mag2() - 4 * n0() * n3();
      }
 
    public:
      GeneralizedCircle reversed() const
      {
        return GeneralizedCircle(-n0(), -n12(), -n3());
      }
 
    public:
      GeneralizedCircle conformalTransformed() const
      {
        return GeneralizedCircle(-n3(), -n12(), -n0());
      }
 
    public:
      ROOT::Math::XYVector gradient(const ROOT::Math::XYVector& point) const
      {
        return point * (2.0 * n3()) + n12();
      }

      ROOT::Math::XYVector normal(const ROOT::Math::XYVector& point) const
      {
        return VectorUtil::unit(gradient(point));
      }

      ROOT::Math::XYVector tangential(const ROOT::Math::XYVector& point) const
      {
        return VectorUtil::Orthogonal(normal(point));
      }

      ROOT::Math::XYVector closest(const ROOT::Math::XYVector& point) const;

      ROOT::Math::XYVector perigee() const;

      ROOT::Math::XYVector apogee() const;

      EForwardBackward isForwardOrBackwardOf(const ROOT::Math::XYVector& from, const ROOT::Math::XYVector& to) const
      {
        ROOT::Math::XYVector difference = to - from;
        ROOT::Math::XYVector tangentialAtFrom = tangential(from);
        return VectorUtil::isForwardOrBackwardOf(tangentialAtFrom, difference);
      }

      ROOT::Math::XYVector
      chooseNextForwardOf(const ROOT::Math::XYVector& start, const ROOT::Math::XYVector& end1, const ROOT::Math::XYVector& end2) const;

      std::pair<ROOT::Math::XYVector, ROOT::Math::XYVector>
      atCylindricalR(double cylindricalR) const;

      ROOT::Math::XYVector atCylindricalRForwardOf(const ROOT::Math::XYVector& startPoint, double cylindricalR) const;

      double fastDistance(const ROOT::Math::XYVector& point) const
      {
        return n0() + point.Dot(n12()) + point.Mag2() * n3();
      }

      double fastImpact() const
      {
        return n0();
      }

      double distance(const ROOT::Math::XYVector& point) const;

      double distance(double fastDistance) const;

      double fastDistance(const double distance) const
      {
        return (distance * n3() + 1.0) * distance;
      }

      double impact() const
      {
        return distance(fastImpact());
      }

      double d0() const
      {
        return -impact();
      }

      ROOT::Math::XYVector tangential() const
      {
        return VectorUtil::unit(tangential(ROOT::Math::XYVector(0.0, 0.0)));
      }

      double tangentialPhi() const
      {
        return tangential().Phi();
      }

      double minimalCylindricalR() const
      {
        return std::fabs(impact());
      }

      double maximalCylindricalR() const
      {
        return std::fabs(impact() + 2 * radius());
      }

      double absDistance(const ROOT::Math::XYVector& point) const
      {
        return fabs(distance(point));
      }

      ERightLeft isRightOrLeft(const ROOT::Math::XYVector& point) const
      {
        return static_cast<ERightLeft>(sign(fastDistance(point)));
      }

      bool isLeft(const ROOT::Math::XYVector& rhs) const
      {
        return isRightOrLeft(rhs) == ERightLeft::c_Left;
      }

      bool isRight(const ROOT::Math::XYVector& rhs) const
      {
        return isRightOrLeft(rhs) == ERightLeft::c_Right;
      }

      bool isLine() const
      {
        return n3() == 0.0;
      }

      bool isCircle() const
      {
        return n3() != 0.0;
      }

      double radius() const
      {
        return 1 / curvature();
      }

      double absRadius() const
      {
        return fabs(radius());
      }

      double curvature() const
      {
        return 2 * n3();
      }

      double omega() const
      {
        return -curvature();
      }

      ROOT::Math::XYVector center() const
      {
        ROOT::Math::XYVector tmp = n12();
        return tmp / (-2. * n3());
      }

      double perimeter() const
      {
        return 2 * M_PI * radius();
      }

      ERotation orientation() const
      {
        return static_cast<ERotation>(sign(n3()));
      }

      double arcLengthPeriod() const
      {
        return std::fabs(perimeter());
      }

      double arcLengthBetween(const ROOT::Math::XYVector& from, const ROOT::Math::XYVector& to) const;

      double arcLengthTo(const ROOT::Math::XYVector& to) const;

      double arcLengthToCylindricalR(double cylindricalR) const;

      double arcLengthFactor(const double directDistance) const
      {
        return arcLengthFactor(directDistance, curvature());
      }

      static double arcLengthFactor(double directDistance, double curvature);

      std::pair<ROOT::Math::XYVector, ROOT::Math::XYVector> intersections(const GeneralizedCircle& generalizedCircle) const;

      ROOT::Math::XYVector atArcLength(double arcLength) const;
 
    private:
      // Order of this parameters make them easier to initialize

      double m_n3 = 0.0;

      ROOT::Math::XYVector m_n12;

      double m_n0 = 0.0;
    };

    std::ostream& operator<<(std::ostream& output, const GeneralizedCircle& circle);
  }
}