algorithms/reco/Helix.h

0001 // FIXME SPDX-License-Identifier: TBD
0002 // Copyright (C) 1997 - 2025 The STAR Collaboration, Xin Dong, Rongrong Ma
0003
0004 #pragma once
0005
0006 #include <edm4eic/ReconstructedParticleCollection.h>
0007 #include <edm4eic/TrackParametersCollection.h>
0008 #include <edm4eic/unit_system.h>
0009 #include <edm4hep/Vector3f.h>
0010 #include <stdlib.h>
0011 #include <cmath>
0012 #include <iterator>
0013 #include <utility>
0014
0015 namespace eicrecon {
0016
0017 class Helix {
0018   bool mSingularity; // true for straight line case (B=0)
0019   edm4hep::Vector3f mOrigin;
0020   double mDipAngle;
0021   double mCurvature;
0022   double mPhase;
0023   int mH; // -sign(q*B);
0024
0025   double mCosDipAngle;
0026   double mSinDipAngle;
0027   double mCosPhase;
0028   double mSinPhase;
0029
0030 public:
0031   /// curvature, dip angle, phase, origin, h
0032   Helix(const double c, const double dip, const double phase, const edm4hep::Vector3f& o,
0033         const int h = -1);
0034
0035   /// momentum, origin, b_field, charge
0036   Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q);
0037
0038   /// ReconstructParticle, b field
0039   Helix(const edm4eic::ReconstructedParticle& p, const double b_field);
0040
0041   ~Helix() = default;
0042
0043   double dipAngle() const;
0044   double curvature() const; /// 1/R in xy-plane
0045   double phase() const;     /// aziumth in xy-plane measured from ring center
0046   double xcenter() const;   /// x-center of circle in xy-plane
0047   double ycenter() const;   /// y-center of circle in xy-plane
0048   int h() const;            /// -sign(q*B);
0049
0050   const edm4hep::Vector3f& origin() const; /// starting point
0051
0052   void setParameters(double c, double dip, double phase, const edm4hep::Vector3f& o, int h);
0053
0054   void setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B,
0055                      const int q);
0056
0057   /// edm4eic::TrackParameters, b field
0058   void setParameters(const edm4eic::TrackParameters& trk, const double b_field);
0059
0060   /// coordinates of helix at point s
0061   double x(double s) const;
0062   double y(double s) const;
0063   double z(double s) const;
0064
0065   edm4hep::Vector3f at(double s) const;
0066
0067   /// pointing vector of helix at point s
0068   double cx(double s) const;
0069   double cy(double s) const;
0070   double cz(double s = 0) const;
0071
0072   edm4hep::Vector3f cat(double s) const;
0073
0074   /// returns period length of helix
0075   double period() const;
0076
0077   /// path length at given r (cylindrical r)
0078   std::pair<double, double> pathLength(double r) const;
0079
0080   /// path length at given r (cylindrical r, cylinder axis at x,y)
0081   std::pair<double, double> pathLength(double r, double x, double y);
0082
0083   /// path length at distance of closest approach to a given point
0084   double pathLength(const edm4hep::Vector3f& p, bool scanPeriods = true) const;
0085
0086   /// path length at intersection with plane
0087   double pathLength(const edm4hep::Vector3f& r, const edm4hep::Vector3f& n) const;
0088
0089   /// path length at distance of closest approach in the xy-plane to a given point
0090   double pathLength(double x, double y) const;
0091
0092   /// path lengths at dca between two helices
0093   std::pair<double, double> pathLengths(const Helix&, double minStepSize = 10 * edm4eic::unit::um,
0094                                         double minRange = 10 * edm4eic::unit::cm) const;
0095
0096   /// minimal distance between point and helix
0097   double distance(const edm4hep::Vector3f& p, bool scanPeriods = true) const;
0098
0099   /// checks for valid parametrization
0100   bool valid(double world = 1.e+5) const { return !bad(world); }
0101   int bad(double world = 1.e+5) const;
0102
0103   /// move the origin along the helix to s which becomes then s=0
0104   void moveOrigin(double s);
0105
0106   static const double NoSolution;
0107
0108   void setCurvature(double); /// performs also various checks
0109   void setPhase(double);
0110   void setDipAngle(double);
0111
0112   /// value of S where distance in x-y plane is minimal
0113   double fudgePathLength(const edm4hep::Vector3f&) const;
0114
0115   // Requires:  signed Magnetic Field
0116   edm4hep::Vector3f momentum(double) const;           // returns the momentum at origin
0117   edm4hep::Vector3f momentumAt(double, double) const; // returns momentum at S
0118   int charge(double) const;                           // returns charge of particle
0119   // 2d DCA to x,y point signed relative to curvature
0120   double curvatureSignedDistance(double x, double y);
0121   // 2d DCA to x,y point signed relative to rotation
0122   double geometricSignedDistance(double x, double y);
0123   // 3d DCA to 3d point signed relative to curvature
0124   double curvatureSignedDistance(const edm4hep::Vector3f&);
0125   // 3d DCA to 3d point signed relative to rotation
0126   double geometricSignedDistance(const edm4hep::Vector3f&);
0127
0128   //
0129   void Print() const;
0130
0131 }; // end class Helix
0132
0133 //
0134 //     Inline functions
0135 //
0136 inline int Helix::h() const { return mH; }
0137
0138 inline double Helix::dipAngle() const { return mDipAngle; }
0139
0140 inline double Helix::curvature() const { return mCurvature; }
0141
0142 inline double Helix::phase() const { return mPhase; }
0143
0144 inline double Helix::x(double s) const {
0145   if (mSingularity)
0146     return mOrigin.x - s * mCosDipAngle * mSinPhase;
0147   else
0148     return mOrigin.x + (cos(mPhase + s * mH * mCurvature * mCosDipAngle) - mCosPhase) / mCurvature;
0149 }
0150
0151 inline double Helix::y(double s) const {
0152   if (mSingularity)
0153     return mOrigin.y + s * mCosDipAngle * mCosPhase;
0154   else
0155     return mOrigin.y + (sin(mPhase + s * mH * mCurvature * mCosDipAngle) - mSinPhase) / mCurvature;
0156 }
0157
0158 inline double Helix::z(double s) const { return mOrigin.z + s * mSinDipAngle; }
0159
0160 inline double Helix::cx(double s) const {
0161   if (mSingularity)
0162     return -mCosDipAngle * mSinPhase;
0163   else
0164     return -sin(mPhase + s * mH * mCurvature * mCosDipAngle) * mH * mCosDipAngle;
0165 }
0166
0167 inline double Helix::cy(double s) const {
0168   if (mSingularity)
0169     return mCosDipAngle * mCosPhase;
0170   else
0171     return cos(mPhase + s * mH * mCurvature * mCosDipAngle) * mH * mCosDipAngle;
0172 }
0173
0174 inline double Helix::cz(double /* s */) const { return mSinDipAngle; }
0175
0176 inline const edm4hep::Vector3f& Helix::origin() const { return mOrigin; }
0177
0178 inline edm4hep::Vector3f Helix::at(double s) const { return edm4hep::Vector3f(x(s), y(s), z(s)); }
0179
0180 inline edm4hep::Vector3f Helix::cat(double s) const {
0181   return edm4hep::Vector3f(cx(s), cy(s), cz(s));
0182 }
0183
0184 inline double Helix::pathLength(double X, double Y) const {
0185   return fudgePathLength(edm4hep::Vector3f(X, Y, 0));
0186 }
0187 inline int Helix::bad(double WorldSize) const {
0188
0189   //    int ierr;
0190   if (!std::isfinite(mDipAngle))
0191     return 11;
0192   if (!std::isfinite(mCurvature))
0193     return 12;
0194
0195   //    ierr = mOrigin.bad(WorldSize);
0196   //    if (ierr)                           return    3+ierr*100;
0197
0198   if (std::abs(mDipAngle) > 1.58)
0199     return 21;
0200   double qwe = std::abs(std::abs(mDipAngle) - M_PI / 2);
0201   if (qwe < 1. / WorldSize)
0202     return 31;
0203
0204   if (std::abs(mCurvature) > WorldSize)
0205     return 22;
0206   if (mCurvature < 0)
0207     return 32;
0208
0209   if (std::abs(mH) != 1)
0210     return 24;
0211
0212   return 0;
0213 }
0214
0215 } // namespace eicrecon