SinEqLine.cc

/**************************************************************************
 * 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.                  *
 **************************************************************************/
#include <tracking/trackingUtilities/numerics/SinEqLine.h>
 
#include <framework/logging/Logger.h>
 
#include <cmath>
 
using namespace Belle2;
using namespace TrackingUtilities;
 
double SinEqLine::computeSmallestPositiveRoot(int maxIHalfPeriod) const
{
  double root = computeRootLargerThanExtemumInHalfPeriod(-1);
  if (root > 0) return root;

  root = computeRootLargerThanExtemumInHalfPeriod(0);
  if (root > 0) return root;
 
  for (int iHalfPeriod = 1; iHalfPeriod <= maxIHalfPeriod; ++iHalfPeriod) {
    root = computeRootLargerThanExtemumInHalfPeriod(iHalfPeriod);

    if (root > 0) return root;
  }
 
  return NAN;
}
 
double SinEqLine::computeRootLargerThanExtemumInHalfPeriod(int iHalfPeriod) const
{
  if (hasLargeSlope()) {
    return computeRootForLargeSlope();
 
  } else {
    double lowerX = computeExtremumXInHalfPeriod(iHalfPeriod);
    double upperX = computeExtremumXInHalfPeriod(iHalfPeriod + 1);
    return computeRootInInterval(lowerX, upperX);
  }
}
 
double SinEqLine::computeRootForLargeSlope() const
{
  if (not hasLargeSlope()) return NAN;
 
  double xLineIsPlusOne = (1.0 - getIntercept()) / getSlope();
  double xLineIsMinusOne = (-1.0 - getIntercept()) / getSlope();
 
  if (xLineIsPlusOne > xLineIsMinusOne) {
    return computeRootInInterval(xLineIsMinusOne, xLineIsPlusOne);
 
  } else if (xLineIsPlusOne < xLineIsMinusOne) {
    return computeRootInInterval(xLineIsPlusOne, xLineIsMinusOne);
  } else {
    return NAN;
  }
}
 
double SinEqLine::computeRootInInterval(double lowerX, double upperX) const
{
  if (not(lowerX < upperX)) return NAN;
 
  ROOT::Math::XYVector lower(lowerX, map(lowerX));
  ROOT::Math::XYVector upper(upperX, map(upperX));

  if (isConverged(lower, upper)) {
    B2INFO("Coverage early");
    return getConvergedBound(lower, upper);
  }

  if (not changesSign(lower, upper)) {
    return NAN;
  }
 
  ROOT::Math::XYVector last(lower);
  ROOT::Math::XYVector current(upper);
 
  ROOT::Math::XYVector next;
  next.SetX(secantX(last, current));
  next.SetY(map(next.x()));
 
  // Should  always succeed since we checked everything before
  bool updatedBound = updateBounds(lower, upper, next);
  if (not updatedBound) return NAN;
 
  while (not isConverged(lower, upper)) {
    // swap accepted values
    last.SetXY(current.X(), current.Y());
    current.SetXY(next.X(), next.Y());
 
    next.SetX(newtonX(current));
    next.SetY(map(next.x()));
 
    updatedBound = updateBounds(lower, upper, next);
 
    if (not updatedBound) {
 
      // fall back to secant.
      next.SetX(secantX(last, current));
      next.SetY(map(next.x()));
 
      updatedBound = updateBounds(lower, upper, next);
 
      if (not updatedBound) {
        // fallback to interval division.
        next.SetX(middleX(lower, upper));
        next.SetY(map(next.x()));
        updatedBound = updateBounds(lower, upper, next);
 
        if (not updatedBound) break;
      }
    }
  }
 
  if (isConverged(lower, upper)) {
    return getConvergedBound(lower, upper);
 
  } else {
    return NAN;
  }
}
 
double SinEqLine::middleX(const ROOT::Math::XYVector& lower, const ROOT::Math::XYVector& upper)
{
  return (lower.x() + upper.x()) / 2.0;
}
 
double SinEqLine::secantX(const ROOT::Math::XYVector& lower, const ROOT::Math::XYVector& upper)
{
  return (lower.x() * upper.y() - upper.x() * lower.y()) / (upper.y() - lower.y());
}
 
double SinEqLine::newtonX(const ROOT::Math::XYVector& pos) const
{
  return pos.x() - pos.y() / gradient(pos.x());
}
 
bool SinEqLine::updateBounds(ROOT::Math::XYVector& lower, ROOT::Math::XYVector& upper, const ROOT::Math::XYVector& next)
{
  if (not isBetween(lower, next, upper)) return false;
 
  EIncDec incDecInfo = getEIncDec(lower, upper);
  if (incDecInfo == EIncDec::c_Increasing) {
    if (next.y() > 0.0 and upper.y() > 0.0) {
      upper.SetXY(next.X(), next.Y());
      return true;
 
    } else if (next.y() <= 0.0 and lower.y() <= 0.0) {
      lower.SetXY(next.X(), next.Y());
      return true;
 
    } else {
      return false;
    }
 
  } else if (incDecInfo == EIncDec::c_Decreasing) {
    if (next.y() >= 0 and lower.y() >= 0.0) {
      lower.SetXY(next.X(), next.Y());
      return true;
 
    } else if (next.y() < 0 and upper.y() < 0) {
      upper.SetXY(next.X(), next.Y());
      return true;
 
    } else {
      return false;
    }
  } else {
    return false;
  }
}
 
double SinEqLine::computeExtremumXInHalfPeriod(int iHalfPeriod) const
{
  const double slope = getSlope();
  double extremumInFirstHalfPeriod = acos(slope);
 
  double extremumInFirstPeriod =
    isEven(iHalfPeriod) ? extremumInFirstHalfPeriod : 2 * M_PI - extremumInFirstHalfPeriod;
 
  int iPeriod = getIPeriodFromIHalfPeriod(iHalfPeriod);
  return extremumInFirstPeriod + 2 * M_PI * iPeriod;
}