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/trackFindingCDC/geometry/Line2D.h>
#include <tracking/trackFindingCDC/geometry/Vector2D.h>
 
#include <tracking/trackFindingCDC/numerics/EForwardBackward.h>
#include <tracking/trackFindingCDC/numerics/ERightLeft.h>
#include <tracking/trackFindingCDC/numerics/ERotation.h>
#include <tracking/trackFindingCDC/numerics/ESign.h>
 
#include <utility>
#include <iosfwd>
#include <cmath>
 
namespace Belle2 {
  namespace TrackFindingCDC {
    class Circle2D;
 
    class GeneralizedCircle {
 
    public:
      GeneralizedCircle();
 
      GeneralizedCircle(double n0, double n1, double n2, double n3 = 0);
 
      GeneralizedCircle(double n0, const Vector2D& n12, double n3 = 0);
 
      explicit GeneralizedCircle(const Line2D& n012);
 
      explicit GeneralizedCircle(const Circle2D& circle);
 
      static GeneralizedCircle
      fromCenterAndRadius(const Vector2D& center,
                          double absRadius,
                          ERotation orientation = ERotation::c_CounterClockwise);
 
      static GeneralizedCircle
      fromPerigeeParameters(double curvature, const Vector2D& 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 Vector2D& n12)
      {
        m_n12.setXY(n12);
      }
 
      void setN3(const double n3)
      {
        m_n3 = n3;
      }
 
    public:
      void setCenterAndRadius(const Vector2D& center,
                              double absRadius,
                              ERotation orientation = ERotation::c_CounterClockwise);
 
      void setPerigeeParameters(double curvature, const Vector2D& tangential, double impact);
 
      void
      setPerigeeParameters(const double curvature, const double tangentialPhi, const double impact)
      {
        setPerigeeParameters(curvature, Vector2D::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 Vector2D& 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 Vector2D& 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 Vector2D& n12() const
      {
        return m_n12;
      }
 
      double n3() const
      {
        return m_n3;
      }
 
    public:
      bool isInvalid() const
      {
        return n0() == 0 and n12().isNull() and n3() == 0;
      }
 
      bool isValid() const
      {
        return not isInvalid();
      }
 
      double normalizationSquared() const
      {
        return n12().normSquared() - 4 * n0() * n3();
      }
 
    public:
      GeneralizedCircle reversed() const
      {
        return GeneralizedCircle(-n0(), -n12(), -n3());
      }
 
    public:
      GeneralizedCircle conformalTransformed() const
      {
        return GeneralizedCircle(-n3(), -n12(), -n0());
      }
 
    public:
      Vector2D gradient(const Vector2D& point) const
      {
        return point * (2.0 * n3()) + n12();
      }
 
      Vector2D normal(const Vector2D& point) const
      {
        return gradient(point).unit();
      }
 
      Vector2D tangential(const Vector2D& point) const
      {
        return normal(point).orthogonal();
      }
 
      Vector2D closest(const Vector2D& point) const;
 
      Vector2D perigee() const;
 
      Vector2D apogee() const;
 
      EForwardBackward isForwardOrBackwardOf(const Vector2D& from, const Vector2D& to) const
      {
        Vector2D difference = to - from;
        Vector2D tangentialAtFrom = tangential(from);
        return tangentialAtFrom.isForwardOrBackwardOf(difference);
      }
 
      Vector2D
      chooseNextForwardOf(const Vector2D& start, const Vector2D& end1, const Vector2D& end2) const;
 
      std::pair<Belle2::TrackFindingCDC::Vector2D, Belle2::TrackFindingCDC::Vector2D>
      atCylindricalR(double cylindricalR) const;
 
      Vector2D atCylindricalRForwardOf(const Vector2D& startPoint, double cylindricalR) const;
 
      double fastDistance(const Vector2D& point) const
      {
        return n0() + point.dot(n12()) + point.normSquared() * n3();
      }
 
      double fastImpact() const
      {
        return n0();
      }
 
      double distance(const Vector2D& 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();
      }
 
      Vector2D tangential() const
      {
        return tangential(Vector2D(0.0, 0.0)).unit();
      }
 
      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 Vector2D& point) const
      {
        return fabs(distance(point));
      }
 
      ERightLeft isRightOrLeft(const Vector2D& point) const
      {
        return static_cast<ERightLeft>(sign(fastDistance(point)));
      }
 
      bool isLeft(const Vector2D& rhs) const
      {
        return isRightOrLeft(rhs) == ERightLeft::c_Left;
      }
 
      bool isRight(const Vector2D& 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();
      }
 
      Vector2D center() const
      {
        return n12().divided(-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 Vector2D& from, const Vector2D& to) const;
 
      double arcLengthTo(const Vector2D& 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<Vector2D, Vector2D> intersections(const GeneralizedCircle& generalizedCircle) const;
 
      Vector2D atArcLength(double arcLength) const;
 
    private:
      // Order of this parameters make them easier to initialize
 
      double m_n3;
 
      Vector2D m_n12;
 
      double m_n0;
    };
 
    std::ostream& operator<<(std::ostream& output, const GeneralizedCircle& circle);
  }
}