PerigeeCircle.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/GeneralizedCircle.h>
#include <tracking/trackFindingCDC/geometry/Line2D.h>
 
#include <tracking/trackFindingCDC/geometry/PerigeeParameters.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 <cmath>
#include <iosfwd>
 
namespace Belle2 {
  namespace TrackFindingCDC {
    class Circle2D;
 
    class PerigeeCircle {
 
    public:
      PerigeeCircle();
 
      PerigeeCircle(double curvature, const Vector2D& phi0Vec, double impact);
 
      PerigeeCircle(double curvature, double phi0, double impact);
 
      explicit PerigeeCircle(const PerigeeParameters& perigeeParameters);
 
    private:
      PerigeeCircle(double curvature, double phi0, const Vector2D& phi0Vec, double impact);
 
    public:
      explicit PerigeeCircle(const Line2D& n012);
 
      explicit PerigeeCircle(const GeneralizedCircle& n0123);
 
      explicit PerigeeCircle(const Circle2D& circle);
 
      static PerigeeCircle fromN(double n0, double n1, double n2, double n3 = 0);
 
      static PerigeeCircle fromN(double n0, const Vector2D& n12, double n3 = 0);
 
      static PerigeeCircle
      fromCenterAndRadius(const Vector2D& center,
                          double absRadius,
                          ERotation orientation = ERotation::c_CounterClockwise);
 
    public:
      void invalidate();
 
      bool isInvalid() const;
 
      bool isValid() const
      {
        return not isInvalid();
      }
 
      void reverse();
 
      PerigeeCircle reversed() const;
 
    public:
      void conformalTransform();
 
      PerigeeCircle conformalTransformed() const;
 
      void passiveMoveBy(const Vector2D& by);
 
      PerigeeJacobian passiveMoveByJacobian(const Vector2D& by) const;
 
      void passiveMoveByJacobian(const Vector2D& by, PerigeeJacobian& jacobian) const;
 
      Vector2D atArcLength(double arcLength) const;
 
      double arcLengthTo(const Vector2D& point) const;
 
      double arcLengthBetween(const Vector2D& from, const Vector2D& to) const;
 
      double arcLengthToCylindricalR(double cylindricalR) const;
 
      std::pair<Vector2D, Vector2D> atCylindricalR(double cylindricalR) const;
 
      Vector2D atCylindricalRForwardOf(const Vector2D& startPoint, double cylindricalR) const;
 
      EForwardBackward isForwardOrBackwardOf(const Vector2D& from, const Vector2D& to) const
      {
        return tangential(from).isForwardOrBackwardOf(to - from);
      }
 
      EForwardBackward isForwardOrBackward(const Vector2D& to) const
      {
        return tangential().isForwardOrBackwardOf(to);
      }
 
      Vector2D
      chooseNextForwardOf(const Vector2D& start, const Vector2D& end1, const Vector2D& end2) const;
      Vector2D closest(const Vector2D& point) const;
 
      double distance(const Vector2D& point) const
      {
        return distance(fastDistance(point));
      }
 
      double distance(double fastDistance) const;
 
      double fastDistance(const Vector2D& point) const;
 
      double fastImpact() const
      {
        return fastDistance(impact());
      }
 
      double fastDistance(double distance) const
      {
        return distance * (1.0 + distance * curvature() / 2);
      }
 
      ERightLeft isRightOrLeft(const Vector2D& point) const
      {
        return static_cast<ERightLeft>(sign(fastDistance(point)));
      }
 
      bool isLine() const
      {
        return curvature() == 0.0;
      }
 
      bool isCircle() const
      {
        return curvature() != 0.0;
      }
 
      ERotation orientation() const
      {
        return static_cast<ERotation>(sign(curvature()));
      }
 
      Vector2D gradient(const Vector2D& point) const
      {
        return (point - perigee()) * curvature() - phi0Vec().orthogonal();
      }
 
      Vector2D normal(const Vector2D& point) const
      {
        return gradient(point).unit();
      }
 
      Vector2D tangential(const Vector2D& point) const
      {
        return normal(point).orthogonal();
      }
 
      const Vector2D& tangential() const
      {
        return phi0Vec();
      }
 
      Vector2D perigee() const
      {
        return phi0Vec().orthogonal() * impact();
      }
 
      Vector2D center() const
      {
        return phi0Vec().orthogonal() * (impact() + radius());
      }
 
      Vector2D apogee() const
      {
        return phi0Vec().orthogonal() * (impact() + 2 * radius());
      }
 
      double minimalCylindricalR() const
      {
        return fabs(impact());
      }
 
      double maximalCylindricalR() const
      {
        return fabs(impact() + 2 * radius());
      }
 
      double arcLengthPeriod() const
      {
        return std::fabs(perimeter());
      }
 
      double perimeter() const
      {
        return 2 * M_PI * radius();
      }
 
      double radius() const
      {
        return 1 / curvature();
      }
 
      double absRadius() const
      {
        return fabs(radius());
      }
 
      void setCenterAndRadius(const Vector2D& center,
                              double absRadius,
                              ERotation orientation = ERotation::c_CounterClockwise);
 
      double n0() const
      {
        return impact() * (impact() * curvature() / 2.0 + 1.0);
      }
 
      Vector2D n12() const
      {
        return -phi0Vec().orthogonal() * (1 + curvature() * impact());
      }
 
      double n1() const
      {
        return n12().x();
      }
      //{ return phi0Vec().y() * (1 + curvature() * impact()); }
 
      double n2() const
      {
        return n12().y();
      }
      //{ return -(phi0Vec().x()) * (1 + curvature() * impact()); }
 
      double n3() const
      {
        return curvature() / 2.0;
      }
 
      void setN(double n0, double n1, double n2, double n3 = 0.0)
      {
        setN(n0, Vector2D(n1, n2), n3);
      }
 
      void setN(double n0, const Vector2D& n12, double n3 = 0.0);
 
      void setN(const Line2D& n012)
      {
        setN(n012.n0(), n012.n12());
      }
 
      void setN(const GeneralizedCircle& n0123)
      {
        setN(n0123.n0(), n0123.n12(), n0123.n3());
      }
 
      double omega() const
      {
        return -curvature();
      }
 
      double d0() const
      {
        return -impact();
      }
 
      double curvature() const
      {
        return m_curvature;
      }
 
      double phi0() const
      {
        return m_phi0;
      }
 
      const Vector2D& phi0Vec() const
      {
        return m_phi0Vec;
      }
 
      double impact() const
      {
        return m_impact;
      }
 
      PerigeeParameters perigeeParameters() const
      {
        using namespace NPerigeeParameterIndices;
        PerigeeParameters result;
        result(c_Curv) = curvature();
        result(c_Phi0) = phi0();
        result(c_I) = impact();
        return result;
      }
 
      void setCurvature(double curvature)
      {
        m_curvature = curvature;
      }
 
      void setPhi0(double phi0)
      {
        m_phi0 = phi0;
        m_phi0Vec = Vector2D::Phi(phi0);
      }
 
      void setPhi0(const Vector2D& phi0Vec)
      {
        m_phi0 = phi0Vec.phi();
        m_phi0Vec = phi0Vec.unit();
      }
 
      void setImpact(double impact)
      {
        m_impact = impact;
      }
 
      void setPerigeeParameters(double curvature, const Vector2D& phi0Vec, double impact)
      {
        m_impact = impact;
        m_phi0 = phi0Vec.phi();
        m_phi0Vec = phi0Vec.unit();
        m_curvature = curvature;
      }
 
      void setPerigeeParameters(double curvature, double phi0, double impact)
      {
        m_impact = impact;
        m_phi0 = phi0;
        m_phi0Vec = Vector2D::Phi(phi0);
        m_curvature = curvature;
      }
 
    private:
      double arcLengthAtDeltaLength(double delta, double dr) const;
 
      double arcLengthAtSecantLength(double secantLength) const;
 
    private:
      double m_curvature;
 
      double m_phi0;
 
      Vector2D m_phi0Vec;
 
      double m_impact;
 
    };
 
    std::ostream& operator<<(std::ostream& output, const PerigeeCircle& circle);
  }
}