Belle II Software  release-05-01-25
SinEqLine.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2014 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Oliver Frost *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 #include <tracking/trackFindingCDC/numerics/SinEqLine.h>
11 
12 #include <framework/logging/Logger.h>
13 
14 #include <cmath>
15 
16 using namespace Belle2;
17 using namespace TrackFindingCDC;
18 
19 double SinEqLine::computeSmallestPositiveRoot(int maxIHalfPeriod) const
20 {
23  if (root > 0) return root;
24 
28  if (root > 0) return root;
29 
30  for (int iHalfPeriod = 1; iHalfPeriod <= maxIHalfPeriod; ++iHalfPeriod) {
31  root = computeRootLargerThanExtemumInHalfPeriod(iHalfPeriod);
32 
34  if (root > 0) return root;
35  }
36 
37  return NAN;
38 }
39 
41 {
42  if (hasLargeSlope()) {
43  return computeRootForLargeSlope();
44 
45  } else {
46  double lowerX = computeExtremumXInHalfPeriod(iHalfPeriod);
47  double upperX = computeExtremumXInHalfPeriod(iHalfPeriod + 1);
48  return computeRootInInterval(lowerX, upperX);
49  }
50 }
51 
53 {
54  if (not hasLargeSlope()) return NAN;
55 
56  double xLineIsPlusOne = (1.0 - getIntercept()) / getSlope();
57  double xLineIsMinusOne = (-1.0 - getIntercept()) / getSlope();
58 
59  if (xLineIsPlusOne > xLineIsMinusOne) {
60  return computeRootInInterval(xLineIsMinusOne, xLineIsPlusOne);
61 
62  } else if (xLineIsPlusOne < xLineIsMinusOne) {
63  return computeRootInInterval(xLineIsPlusOne, xLineIsMinusOne);
64  } else {
65  return NAN;
66  }
67 }
68 
69 double SinEqLine::computeRootInInterval(double lowerX, double upperX) const
70 {
71  if (not(lowerX < upperX)) return NAN;
72 
73  Vector2D lower(lowerX, map(lowerX));
74  Vector2D upper(upperX, map(upperX));
75 
78  if (isConverged(lower, upper)) {
79  B2INFO("Coverage early");
80  return getConvergedBound(lower, upper);
81  }
82 
84  if (not changesSign(lower, upper)) {
85  return NAN;
86  }
87 
88  Vector2D last(lower);
89  Vector2D current(upper);
90 
91  Vector2D next;
92  next.setX(secantX(last, current));
93  next.setY(map(next.x()));
94 
95  // Should always succeed since we checked everything before
96  // cppcheck-suppress unreadVariable
97  bool updatedBound = updateBounds(lower, upper, next);
98  if (not updatedBound) return NAN;
99 
100  while (not isConverged(lower, upper)) {
101  // swap accepted values
102  last.set(current);
103  current.set(next);
104 
105  next.setX(newtonX(current));
106  next.setY(map(next.x()));
107 
108  updatedBound = updateBounds(lower, upper, next);
109 
110  if (not updatedBound) {
111 
112  // fall back to sekant.
113  next.setX(secantX(last, current));
114  next.setY(map(next.x()));
115 
116  updatedBound = updateBounds(lower, upper, next);
117 
118  if (not updatedBound) {
119  // fallback to interval devision.
120  next.setX(middleX(lower, upper));
121  next.setY(map(next.x()));
122  updatedBound = updateBounds(lower, upper, next);
123 
124  if (not updatedBound) break;
125  }
126  }
127  }
128 
129  if (isConverged(lower, upper)) {
130  return getConvergedBound(lower, upper);
131 
132  } else {
133  return NAN;
134  }
135 }
136 
137 double SinEqLine::middleX(const Vector2D& lower, const Vector2D& upper)
138 {
139  return (lower.x() + upper.x()) / 2.0;
140 }
141 
142 double SinEqLine::secantX(const Vector2D& lower, const Vector2D& upper)
143 {
144  return (lower.x() * upper.y() - upper.x() * lower.y()) / (upper.y() - lower.y());
145 }
146 
147 double SinEqLine::newtonX(const Vector2D& pos) const
148 {
149  return pos.x() - pos.y() / gradient(pos.x());
150 }
151 
152 bool SinEqLine::updateBounds(Vector2D& lower, Vector2D& upper, const Vector2D& next)
153 {
155  if (not isBetween(lower, next, upper)) return false;
156 
157  EIncDec incDecInfo = getEIncDec(lower, upper);
158  if (incDecInfo == EIncDec::c_Increasing) {
159  if (next.y() > 0.0 and upper.y() > 0.0) {
160  upper.set(next);
161  return true;
162 
163  } else if (next.y() <= 0.0 and lower.y() <= 0.0) {
164  lower.set(next);
165  return true;
166 
167  } else {
168  return false;
169  }
170 
171  } else if (incDecInfo == EIncDec::c_Decreasing) {
172  if (next.y() >= 0 and lower.y() >= 0.0) {
173  lower.set(next);
174  return true;
175 
176  } else if (next.y() < 0 and upper.y() < 0) {
177  upper.set(next);
178  return true;
179 
180  } else {
181  return false;
182  }
183  } else {
184  return false;
185  }
186 }
187 
188 double SinEqLine::computeExtremumXInHalfPeriod(int iHalfPeriod) const
189 {
190  const double slope = getSlope();
191  double extremumInFirstHalfPeriod = acos(slope);
192 
193  double extremumInFirstPeriod =
194  isEven(iHalfPeriod) ? extremumInFirstHalfPeriod : 2 * M_PI - extremumInFirstHalfPeriod;
195 
196  int iPeriod = getIPeriodFromIHalfPeriod(iHalfPeriod);
197  return extremumInFirstPeriod + 2 * M_PI * iPeriod;
198 }
Belle2::TrackFindingCDC::SinEqLine::getSlope
double getSlope() const
Getter for the slope.
Definition: SinEqLine.h:170
Belle2::TrackFindingCDC::SinEqLine::changesSign
static bool changesSign(const Vector2D &lower, const Vector2D &upper)
Checks if the function changes sign in the intervall.
Definition: SinEqLine.h:139
Belle2::TrackFindingCDC::SinEqLine::isConverged
static bool isConverged(const Vector2D &lower, const Vector2D &upper)
Check if the intervall has shrunk close enough to the solution.
Definition: SinEqLine.h:114
Belle2::TrackFindingCDC::Vector2D::setY
void setY(const double y)
Setter for the y coordinate.
Definition: Vector2D.h:624
Belle2::TrackFindingCDC::SinEqLine::newtonX
double newtonX(const Vector2D &pos) const
Shrinking method of the newton algorithm return the next candidate root.
Definition: SinEqLine.cc:147
Belle2::TrackFindingCDC::Vector2D
A two dimensional vector which is equipped with functions for correct handeling of orientation relat...
Definition: Vector2D.h:37
Belle2::TrackFindingCDC::SinEqLine::getIPeriodFromIHalfPeriod
static int getIPeriodFromIHalfPeriod(int iHalfPeriod)
Helper function to translate the index of the half period to index of the containing period.
Definition: SinEqLine.h:161
Belle2::TrackFindingCDC::NIncDec::EIncDec
EIncDec
Enumeration to represent the distinct possibilities of the right left passage information.
Definition: EIncDec.h:35
Belle2::TrackFindingCDC::Vector2D::y
double y() const
Getter for the y coordinate.
Definition: Vector2D.h:619
Belle2::TrackFindingCDC::SinEqLine::updateBounds
static bool updateBounds(Vector2D &lower, Vector2D &upper, const Vector2D &next)
Replaces the lower or upper bound inplace if the next candidate position is valid and within the inte...
Definition: SinEqLine.cc:152
Belle2::TrackFindingCDC::SinEqLine::map
double map(const double x) const
Interpreting as the function f this method carries out the translation from x to y coordinates.
Definition: SinEqLine.h:70
Belle2::TrackFindingCDC::SinEqLine::isBetween
static bool isBetween(const Vector2D &lower, const Vector2D &next, const Vector2D &upper)
Check is next position is within the intervall given by lower and upper.
Definition: SinEqLine.h:110
Belle2::TrackFindingCDC::SinEqLine::gradient
double gradient(const double x) const
Interpreting as the function f this method calculates the gradient as need in Newtons algorithms.
Definition: SinEqLine.h:74
Belle2::TrackFindingCDC::SinEqLine::computeSmallestPositiveRoot
double computeSmallestPositiveRoot(int maxIHalfPeriod=5) const
Definition: SinEqLine.cc:19
Belle2::TrackFindingCDC::SinEqLine::getIntercept
double getIntercept() const
Getter for the intercept.
Definition: SinEqLine.h:174
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TrackFindingCDC::SinEqLine::middleX
static double middleX(const Vector2D &lower, const Vector2D &upper)
Simple fall back shrinking method using trivial devision of the intervall.
Definition: SinEqLine.cc:137
Belle2::TrackFindingCDC::SinEqLine::secantX
static double secantX(const Vector2D &lower, const Vector2D &upper)
Fall back shrinking method to the secant algorithm.
Definition: SinEqLine.cc:142
Belle2::TrackFindingCDC::SinEqLine::computeRootLargerThanExtemumInHalfPeriod
double computeRootLargerThanExtemumInHalfPeriod(int iHalfPeriod) const
Computes the solution that is addressed by the given half period index.
Definition: SinEqLine.cc:40
Belle2::TrackFindingCDC::SinEqLine::hasLargeSlope
bool hasLargeSlope() const
Indicates that the slope is so large such that the function has no local exterma.
Definition: SinEqLine.h:166
Belle2::TrackFindingCDC::SinEqLine::computeRootInInterval
double computeRootInInterval(double lowerX, double upperX) const
Computes the solution in between the given x values. The x values are generally choosen consecutive l...
Definition: SinEqLine.cc:69
Belle2::TrackFindingCDC::Vector2D::x
double x() const
Getter for the x coordinate.
Definition: Vector2D.h:609
Belle2::TrackFindingCDC::SinEqLine::getEIncDec
static EIncDec getEIncDec(const Vector2D &lower, const Vector2D &upper)
Determines if the function is increasing or decreasing in the intervall.
Definition: SinEqLine.h:143
Belle2::TrackFindingCDC::Vector2D::setX
void setX(const double x)
Setter for the x coordinate.
Definition: Vector2D.h:614
Belle2::TrackFindingCDC::SinEqLine::getConvergedBound
static double getConvergedBound(const Vector2D &lower, const Vector2D &upper)
Returns the better solution x from the bounds of the intervall.
Definition: SinEqLine.h:120
Belle2::TrackFindingCDC::Vector2D::set
void set(const double first, const double second)
Setter for both coordinate.
Definition: Vector2D.h:664
Belle2::TrackFindingCDC::SinEqLine::computeRootForLargeSlope
double computeRootForLargeSlope() const
Compute single solution in the case that fabs(slope) >= 1.
Definition: SinEqLine.cc:52
Belle2::TrackFindingCDC::SinEqLine::computeExtremumXInHalfPeriod
double computeExtremumXInHalfPeriod(int iHalfPeriod) const
Get the local extermum that is located in the half period indicated by the given index.
Definition: SinEqLine.cc:188