Helix.h

Go to the documentation of this file.

// FIXME SPDX-License-Identifier: TBD
// Copyright (C) 1997 - 2025 The STAR Collaboration, Xin Dong, Rongrong Ma
 
#pragma once
 
#include <edm4eic/ReconstructedParticleCollection.h>
#include <edm4eic/TrackParametersCollection.h>
#include <edm4eic/unit_system.h>
#include <edm4hep/Vector3f.h>
#include <stdlib.h>
#include <cmath>
#include <iterator>
#include <utility>
 
namespace eicrecon {
 
class Helix {
  bool mSingularity; // true for straight line case (B=0)
  edm4hep::Vector3f mOrigin;
  double mDipAngle;
  double mCurvature;
  double mPhase;
  int mH; // -sign(q*B);
 
  double mCosDipAngle;
  double mSinDipAngle;
  double mCosPhase;
  double mSinPhase;
 
public:
  Helix(const double c, const double dip, const double phase, const edm4hep::Vector3f& o,
        const int h = -1);
 
  Helix(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B, const int q);
 
  Helix(const edm4eic::ReconstructedParticle& p, const double b_field);
 
  ~Helix() = default;
 
  double dipAngle() const;
  double curvature() const; 
  double phase() const;     
  double xcenter() const;   
  double ycenter() const;   
  int h() const;            
 
  const edm4hep::Vector3f& origin() const; 
 
  void setParameters(double c, double dip, double phase, const edm4hep::Vector3f& o, int h);
 
  void setParameters(const edm4hep::Vector3f& p, const edm4hep::Vector3f& o, const double B,
                     const int q);
 
  void setParameters(const edm4eic::TrackParameters& trk, const double b_field);
 
  double x(double s) const;
  double y(double s) const;
  double z(double s) const;
 
  edm4hep::Vector3f at(double s) const;
 
  double cx(double s) const;
  double cy(double s) const;
  double cz(double s = 0) const;
 
  edm4hep::Vector3f cat(double s) const;
 
  double period() const;
 
  std::pair<double, double> pathLength(double r) const;
 
  std::pair<double, double> pathLength(double r, double x, double y);
 
  double pathLength(const edm4hep::Vector3f& p, bool scanPeriods = true) const;
 
  double pathLength(const edm4hep::Vector3f& r, const edm4hep::Vector3f& n) const;
 
  double pathLength(double x, double y) const;
 
  std::pair<double, double> pathLengths(const Helix&, double minStepSize = 10 * edm4eic::unit::um,
                                        double minRange = 10 * edm4eic::unit::cm) const;
 
  double distance(const edm4hep::Vector3f& p, bool scanPeriods = true) const;
 
  bool valid(double world = 1.e+5) const { return !bad(world); }
  int bad(double world = 1.e+5) const;
 
  void moveOrigin(double s);
 
  static const double NoSolution;
 
  void setCurvature(double); 
  void setPhase(double);
  void setDipAngle(double);
 
  double fudgePathLength(const edm4hep::Vector3f&) const;
 
  // Requires:  signed Magnetic Field
  edm4hep::Vector3f momentum(double) const;           // returns the momentum at origin
  edm4hep::Vector3f momentumAt(double, double) const; // returns momentum at S
  int charge(double) const;                           // returns charge of particle
  // 2d DCA to x,y point signed relative to curvature
  double curvatureSignedDistance(double x, double y);
  // 2d DCA to x,y point signed relative to rotation
  double geometricSignedDistance(double x, double y);
  // 3d DCA to 3d point signed relative to curvature
  double curvatureSignedDistance(const edm4hep::Vector3f&);
  // 3d DCA to 3d point signed relative to rotation
  double geometricSignedDistance(const edm4hep::Vector3f&);
 
  //
  void Print() const;
 
}; // end class Helix
 
//
//     Inline functions
//
inline int Helix::h() const { return mH; }
 
inline double Helix::dipAngle() const { return mDipAngle; }
 
inline double Helix::curvature() const { return mCurvature; }
 
inline double Helix::phase() const { return mPhase; }
 
inline double Helix::x(double s) const {
  if (mSingularity)
    return mOrigin.x - s * mCosDipAngle * mSinPhase;
  else
    return mOrigin.x + (cos(mPhase + s * mH * mCurvature * mCosDipAngle) - mCosPhase) / mCurvature;
}
 
inline double Helix::y(double s) const {
  if (mSingularity)
    return mOrigin.y + s * mCosDipAngle * mCosPhase;
  else
    return mOrigin.y + (sin(mPhase + s * mH * mCurvature * mCosDipAngle) - mSinPhase) / mCurvature;
}
 
inline double Helix::z(double s) const { return mOrigin.z + s * mSinDipAngle; }
 
inline double Helix::cx(double s) const {
  if (mSingularity)
    return -mCosDipAngle * mSinPhase;
  else
    return -sin(mPhase + s * mH * mCurvature * mCosDipAngle) * mH * mCosDipAngle;
}
 
inline double Helix::cy(double s) const {
  if (mSingularity)
    return mCosDipAngle * mCosPhase;
  else
    return cos(mPhase + s * mH * mCurvature * mCosDipAngle) * mH * mCosDipAngle;
}
 
inline double Helix::cz(double /* s */) const { return mSinDipAngle; }
 
inline const edm4hep::Vector3f& Helix::origin() const { return mOrigin; }
 
inline edm4hep::Vector3f Helix::at(double s) const { return edm4hep::Vector3f(x(s), y(s), z(s)); }
 
inline edm4hep::Vector3f Helix::cat(double s) const {
  return edm4hep::Vector3f(cx(s), cy(s), cz(s));
}
 
inline double Helix::pathLength(double X, double Y) const {
  return fudgePathLength(edm4hep::Vector3f(X, Y, 0));
}
inline int Helix::bad(double WorldSize) const {
 
  //    int ierr;
  if (!std::isfinite(mDipAngle))
    return 11;
  if (!std::isfinite(mCurvature))
    return 12;
 
  //    ierr = mOrigin.bad(WorldSize);
  //    if (ierr)                           return    3+ierr*100;
 
  if (std::abs(mDipAngle) > 1.58)
    return 21;
  double qwe = std::abs(std::abs(mDipAngle) - M_PI / 2);
  if (qwe < 1. / WorldSize)
    return 31;
 
  if (std::abs(mCurvature) > WorldSize)
    return 22;
  if (mCurvature < 0)
    return 32;
 
  if (std::abs(mH) != 1)
    return 24;
 
  return 0;
}
 
} // namespace eicrecon