Belle II Software  release-05-01-25
Helix.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/geometry/Helix.h>
11 
12 #include <tracking/trackFindingCDC/geometry/HelixParameters.h>
13 #include <tracking/trackFindingCDC/geometry/PerigeeCircle.h>
14 #include <tracking/trackFindingCDC/geometry/PerigeeParameters.h>
15 #include <tracking/trackFindingCDC/geometry/SZLine.h>
16 #include <tracking/trackFindingCDC/geometry/SZParameters.h>
17 
18 #include <tracking/trackFindingCDC/geometry/Vector3D.h>
19 
20 #include <boost/math/tools/minima.hpp>
21 
22 #include <utility>
23 #include <algorithm>
24 #include <limits>
25 #include <ostream>
26 #include <cstdint>
27 
28 using namespace Belle2;
29 using namespace TrackFindingCDC;
30 
31 double Helix::arcLength2DToClosest(const Vector3D& point, bool firstPeriod) const
32 {
33  // The point may happen to lie in the center of the helix.
34  double d0 = circleXY().distance(point.xy());
35  double denominator = 1 + curvatureXY() * d0;
36  if (denominator == 0) {
37  // When this happens we can optimise the z distance for the closest point
38  double arcLength2D = (point.z() - z0()) / tanLambda();
39  if (not std::isfinite(arcLength2D)) {
40  return 0.0;
41  } else {
42  return arcLength2D;
43  }
44  }
45 
46  // First approximation optimising the xy distance.
47  double arcLength2D = circleXY().arcLengthTo(point.xy());
48  // In case the helix is a flat circle with no extend into z this is the actual solution
49  if (tanLambda() == 0) {
50  return arcLength2D;
51  }
52 
53  double deltaZ = point.z() - szLine().map(arcLength2D);
54  if (not firstPeriod) {
55  if (fabs(deltaZ) > zPeriod() / 2) {
56  double newDeltaZ = std::remainder(deltaZ, zPeriod());
57  double nPeriodShift = (deltaZ - newDeltaZ) / zPeriod();
58  arcLength2D += nPeriodShift * perimeterXY();
59  deltaZ = newDeltaZ;
60  }
61  }
62 
63  using boost::math::tools::brent_find_minima;
64 
65  auto distance3D = [this, &point](const double & s) -> double {
66  Vector3D pos = atArcLength2D(s);
67  return pos.distance(point);
68  };
69 
70  double searchWidth = std::fmin(perimeterXY(), 2 * deltaZ / tanLambda());
71 
72  double lowerS = arcLength2D - searchWidth;
73  double upperS = arcLength2D + searchWidth;
74 
75  int bits = std::numeric_limits<double>::digits;
76  boost::uintmax_t nMaxIter = 100;
77 
78  std::pair<double, double> sBounds = brent_find_minima(distance3D, lowerS, upperS, bits, nMaxIter);
79 
80  // Stopped before iterations were exhausted?
81  // bool converged = nMaxIter > 0;
82 
83  double firstDistance = distance3D(sBounds.first);
84  double secondDistance = distance3D(sBounds.second);
85  if (firstDistance < secondDistance) {
86  return sBounds.first;
87  } else {
88  return sBounds.second;
89  }
90 }
91 
93 {
94  // Fills the upper left 3x3 corner.
95  PerigeeJacobian perigeeJacobian = circleXY().passiveMoveByJacobian(by.xy());
96  SZJacobian szJacobian = SZUtil::identity();
97  HelixJacobian jacobian = HelixUtil::stackBlocks(perigeeJacobian, szJacobian);
98 
99  double curv = curvatureXY();
100  double tanL = tanLambda();
101  double sArc = circleXY().arcLengthTo(by.xy());
102 
103  using namespace NHelixParameterIndices;
104  jacobian(c_Z0, c_Curv) = tanL * (jacobian(c_Phi0, c_Curv) - sArc) / curv;
105  jacobian(c_Z0, c_Phi0) = tanL * (jacobian(c_Phi0, c_Phi0) - 1.) / curv;
106  jacobian(c_Z0, c_I) = tanL * jacobian(c_Phi0, c_I) / curv;
107  jacobian(c_Z0, c_TanL) = sArc;
108 
109  return jacobian;
110 }
111 
112 std::ostream& TrackFindingCDC::operator<<(std::ostream& output, const Helix& helix)
113 {
114  return output << "Helix("
115  << "curv=" << helix.curvatureXY() << ","
116  << "phi0=" << helix.phi0() << ","
117  << "impact=" << helix.impactXY() << ","
118  << "tanL=" << helix.tanLambda() << ","
119  << "z0=" << helix.z0() << ")";
120 }
Belle2::TrackFindingCDC::Helix::curvatureXY
double curvatureXY() const
Getter for the signed curvature in the xy projection.
Definition: Helix.h:214
Belle2::TrackFindingCDC::PerigeeCircle::arcLengthTo
double arcLengthTo(const Vector2D &point) const
Calculates the arc length between the perigee and the given point.
Definition: PerigeeCircle.cc:238
Belle2::TrackFindingCDC::Helix
Extension of the generalized circle also caching the perigee coordinates.
Definition: Helix.h:38
Belle2::TrackFindingCDC::Helix::tanLambda
double tanLambda() const
Getter for the proportinality factor from arc length in xy space to z.
Definition: Helix.h:250
Belle2::operator<<
std::ostream & operator<<(std::ostream &output, const Helix &helix)
Output operator for debugging and the generation of unittest error messages.
Definition: Helix.cc:632
Belle2::TrackFindingCDC::Helix::arcLength2DToClosest
double arcLength2DToClosest(const Vector3D &point, bool firstPeriod=true) const
Calculates the perpendicular travel distance at which the helix has the closest approach to the given...
Definition: Helix.cc:31
Belle2::TrackFindingCDC::UncertainParametersUtil< SZUtil, ESZParameter >::identity
static CovarianceMatrix identity()
Returns an identity matrix.
Definition: UncertainParameters.icc.h:65
Belle2::TrackFindingCDC::Helix::szLine
const SZLine & szLine() const
Getter for the projection into xy space.
Definition: Helix.h:336
Belle2::TrackFindingCDC::Helix::passiveMoveByJacobian
HelixJacobian passiveMoveByJacobian(const Vector3D &by) const
Computes the Jacobi matrix for a move of the coordinate system by the given vector.
Definition: Helix.cc:92
Belle2::TrackFindingCDC::PerigeeCircle::passiveMoveByJacobian
PerigeeJacobian passiveMoveByJacobian(const Vector2D &by) const
Computes the Jacobi matrix for a move of the coordinate system by the given vector.
Definition: PerigeeCircle.cc:185
Belle2::TrackFindingCDC::Helix::perimeterXY
double perimeterXY() const
Getter for the perimeter of the circle in the xy projection.
Definition: Helix.h:281
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TrackFindingCDC::Helix::z0
double z0() const
Getter for z coordinate at the perigee point of the helix.
Definition: Helix.h:262
Belle2::TrackFindingCDC::Vector3D
A three dimensional vector.
Definition: Vector3D.h:34
Belle2::TrackFindingCDC::SZLine::map
double map(double s) const
Maps the two dimensional arc length s to z.
Definition: SZLine.h:126
Belle2::TrackFindingCDC::Helix::zPeriod
double zPeriod() const
Getter for the distance in z at which the two points on the helix coincide in the xy projection.
Definition: Helix.h:269
Belle2::TrackFindingCDC::Helix::circleXY
const PerigeeCircle & circleXY() const
Getter for the projection into xy space.
Definition: Helix.h:330
Belle2::TrackFindingCDC::Vector3D::xy
const Vector2D & xy() const
Getter for the xy projected vector ( reference ! )
Definition: Vector3D.h:500
Belle2::TrackFindingCDC::Helix::impactXY
double impactXY() const
Getter for the signed distance to the z axes at the perigee point.
Definition: Helix.h:220
Belle2::TrackFindingCDC::PlainMatrix
A matrix implementation to be used as an interface typ through out the track finder.
Definition: PlainMatrix.h:50
Belle2::TrackFindingCDC::Vector3D::z
double z() const
Getter for the z coordinate.
Definition: Vector3D.h:488
Belle2::TrackFindingCDC::PerigeeCircle::distance
double distance(const Vector2D &point) const
Getter for the proper signed distance of the point to the circle.
Definition: PerigeeCircle.h:220
Belle2::TrackFindingCDC::Helix::atArcLength2D
Vector3D atArcLength2D(double s) const
Calculates the point, which lies at the give perpendicular travel distance (counted from the perigee)
Definition: Helix.h:184
Belle2::TrackFindingCDC::HelixUtil::stackBlocks
static HelixUtil::CovarianceMatrix stackBlocks(const PerigeeUtil::CovarianceMatrix &perigeeCov, const SZUtil::CovarianceMatrix &szCov)
Combine covariance matrices from the xy space and the sz space in their respective blocks.
Definition: HelixParameters.cc:60
Belle2::TrackFindingCDC::Helix::d0
double d0() const
Getter for the signed distance to the z axes at the perigee point.
Definition: Helix.h:232
Belle2::TrackFindingCDC::Helix::phi0
double phi0() const
Getter for the azimuth angle of the direction of flight at the perigee.
Definition: Helix.h:311