SinEqLine.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/numerics/Modulo.h>
 
#include <tracking/trackingUtilities/geometry/Line2D.h>
 
#include <tracking/trackingUtilities/numerics/EIncDec.h>
 
#include <Math/Vector2D.h>
 
#include <cmath>
 
namespace Belle2 {
 
  namespace TrackingUtilities {

    class SinEqLine  {
 
    public:
      SinEqLine() :
        m_slope(0.0),
        m_intercept(0.0)
      {}
 

      explicit SinEqLine(const Line2D& line2D) :
        m_slope(line2D.slope()),
        m_intercept(line2D.intercept())
      {}

      SinEqLine(const double slope, const double intercept) :
        m_slope(slope),
        m_intercept(intercept)
      {}
 

      double map(const double x) const
      { return sin(x) - getSlope() * x - getIntercept(); }

      double gradient(const double x) const
      { return cos(x) - getSlope(); }

      static int getIHalfPeriod(const double x)
      { return floor(x / M_PI); }
 
      /* Computes the positive solution that has the smallest value of x.
      The additional parameter serves as a criterion to abort the search if the solutions is further away than the specified half period.
      */
      double computeSmallestPositiveRoot(int maxIHalfPeriod = 5) const;
 

      double computeRootLargerThanExtemumInHalfPeriod(int iHalfPeriod) const;

      double computeRootForLargeSlope() const;

      double computeRootInInterval(double lowerX, double upperX) const;
 
    private:
      double newtonX(const ROOT::Math::XYVector& pos) const;

      static double secantX(const ROOT::Math::XYVector& lower, const ROOT::Math::XYVector& upper);

      static double middleX(const ROOT::Math::XYVector& lower, const ROOT::Math::XYVector& upper);

      static bool updateBounds(ROOT::Math::XYVector& lower, ROOT::Math::XYVector& upper, const ROOT::Math::XYVector& next);

      static bool isBetween(const ROOT::Math::XYVector& lower, const ROOT::Math::XYVector& next, const ROOT::Math::XYVector& upper)
      { return lower.x() < next.x() and next.x() < upper.x(); }

      static bool isConverged(const ROOT::Math::XYVector& lower, const ROOT::Math::XYVector& upper)
      {
        return fabs(lower.y()) < 10e-7 or fabs(upper.y()) < 10e-7;
      }

      static double getConvergedBound(const ROOT::Math::XYVector& lower, const ROOT::Math::XYVector& upper)
      {
        if (not std::isfinite(lower.y()) or not std::isfinite(upper.y())) {
          return NAN;
        }
 
        if (fabs(lower.y()) <= fabs(upper.y())) {
          return lower.x();
        }
 
        if (fabs(lower.y()) > fabs(upper.y())) {
          return upper.x();
        }
 
        return NAN;
      }
 
    public:
      static bool changesSign(const ROOT::Math::XYVector& lower, const ROOT::Math::XYVector& upper)
      { return (lower.y() > 0 and upper.y() < 0) or (lower.y() < 0 and upper.y() > 0); }

      static EIncDec getEIncDec(const ROOT::Math::XYVector& lower, const ROOT::Math::XYVector& upper)
      {
        if (lower.y() < upper.y()) {
          return EIncDec::c_Increasing;
        } else if (lower.y() > upper.y()) {
          return EIncDec::c_Decreasing;
        } else if (lower.y() == upper.y()) {
          return EIncDec::c_Constant;
        } else {
          return EIncDec::c_Invalid;
        }
      }
 
    public:
      double computeExtremumXInHalfPeriod(int iHalfPeriod) const;

      static int getIPeriodFromIHalfPeriod(int iHalfPeriod)
      { return isEven(iHalfPeriod) ? iHalfPeriod / 2 : (iHalfPeriod - 1) / 2; }
 
    public:
      bool hasLargeSlope() const
      { return fabs(getSlope()) >= 1; }

      double getSlope() const
      { return m_slope; }

      double getIntercept() const
      { return m_intercept; }
 
    private:
      double m_slope;

      double m_intercept;
 
 
    };
 
  }
}