MpdHelix.h 6.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219
  1. /***************************************************************************
  2. * Armen Kechechyan, September 2009
  3. ***************************************************************************
  4. * Description:
  5. * Parametrization of a helix
  6. **************************************************************************/
  7. #ifndef MPD_HELIX_H
  8. #define MPD_HELIX_H
  9. #include "TVector3.h"
  10. #include <math.h>
  11. #include <utility>
  12. #include <algorithm>
  13. #if !defined(NO_NAMESPACES)
  14. using std::pair;
  15. using std::swap;
  16. using std::max;
  17. #endif
  18. class MpdHelix {
  19. public:
  20. /// curvature, dip angle, phase, origin, h
  21. MpdHelix(double c, double dip, double phase,
  22. const TVector3 o, int h=-1);
  23. MpdHelix(TVector3 mom, TVector3 o, Double_t charge, Double_t Bz=0.5);
  24. virtual ~MpdHelix();
  25. // MpdHelix(const MpdHelix&); // use default
  26. // MpdHelix& operator=(const MpdHelix&); // use default
  27. double dipAngle() const;
  28. double curvature() const; /// 1/R in xy-plane
  29. double phase() const; /// aziumth in xy-plane measured from ring center
  30. double xcenter() const; /// x-center of circle in xy-plane
  31. double ycenter() const; /// y-center of circle in xy-plane
  32. int h() const; /// -sign(q*B);
  33. const TVector3 origin() const; /// starting point
  34. void setParameters(double c, double dip, double phase, const TVector3 o, int h);
  35. /// coordinates of helix at point s
  36. double x(double s) const;
  37. double y(double s) const;
  38. double z(double s) const;
  39. TVector3 at(double s) const;
  40. /// pointing vector of helix at point s
  41. double cx(double s) const;
  42. double cy(double s) const;
  43. double cz(double s) const;
  44. TVector3 cat(double s) const;
  45. /// returns period length of helix
  46. double period() const;
  47. /// path length at given r (cylindrical r)
  48. pair<double, double> pathLength(double r) const;
  49. /// path length at given r (cylindrical r, cylinder axis at x,y)
  50. pair<double, double> pathLength(double r, double x, double y);
  51. /// path length at distance of closest approach to a given point
  52. double pathLength(const TVector3 p, bool scanPeriods = true) const;
  53. /// path length at intersection with plane
  54. double pathLength(const TVector3 r, const TVector3 n) const;
  55. /// path length at distance of closest approach in the xy-plane to a given point
  56. double pathLength(double x, double y) const;
  57. /// path lengths at dca between two helices
  58. pair<double, double> pathLengths(const MpdHelix&) const;
  59. /// minimal distance between point and helix
  60. double distance(const TVector3 p, bool scanPeriods = true) const;
  61. /// checks for valid parametrization
  62. bool valid(double world = 1.e+5) const {return !bad(world);}
  63. int bad(double world = 1.e+5) const;
  64. /// move the origin along the helix to s which becomes then s=0
  65. virtual void moveOrigin(double s);
  66. static const double NoSolution;
  67. protected:
  68. MpdHelix();
  69. void setCurvature(double); /// performs also various checks
  70. void setPhase(double);
  71. void setDipAngle(double);
  72. /// value of S where distance in x-y plane is minimal
  73. double fudgePathLength(const TVector3) const;
  74. protected:
  75. bool mSingularity; // true for straight line case (B=0)
  76. TVector3 mOrigin;
  77. double mDipAngle;
  78. double mCurvature;
  79. double mPhase;
  80. int mH; // -sign(q*B);
  81. double mCosDipAngle;
  82. double mSinDipAngle;
  83. double mCosPhase;
  84. double mSinPhase;
  85. #ifdef __ROOT__
  86. ClassDef(MpdHelix,1)
  87. #endif
  88. };
  89. //
  90. // Non-member functions
  91. //
  92. int operator== (const MpdHelix&, const MpdHelix&);
  93. int operator!= (const MpdHelix&, const MpdHelix&);
  94. //ostream& operator<<(ostream&, const MpdHelix&);
  95. //
  96. // Inline functions
  97. //
  98. inline int MpdHelix::h() const {return mH;}
  99. inline double MpdHelix::dipAngle() const {return mDipAngle;}
  100. inline double MpdHelix::curvature() const {return mCurvature;}
  101. inline double MpdHelix::phase() const {return mPhase;}
  102. inline double MpdHelix::x(double s) const
  103. {
  104. if (mSingularity)
  105. return mOrigin.X() - s*mCosDipAngle*mSinPhase;
  106. else
  107. return mOrigin.X() + (cos(mPhase + s*mH*mCurvature*mCosDipAngle)-mCosPhase)/mCurvature;
  108. }
  109. inline double MpdHelix::y(double s) const
  110. {
  111. if (mSingularity)
  112. return mOrigin.Y() + s*mCosDipAngle*mCosPhase;
  113. else
  114. return mOrigin.Y() + (sin(mPhase + s*mH*mCurvature*mCosDipAngle)-mSinPhase)/mCurvature;
  115. }
  116. inline double MpdHelix::z(double s) const
  117. {
  118. return mOrigin.Z() + s*mSinDipAngle;
  119. }
  120. inline double MpdHelix::cx(double s) const
  121. {
  122. if (mSingularity)
  123. return -mCosDipAngle*mSinPhase;
  124. else
  125. return -sin(mPhase + s*mH*mCurvature*mCosDipAngle)*mH*mCosDipAngle;
  126. }
  127. inline double MpdHelix::cy(double s) const
  128. {
  129. if (mSingularity)
  130. return mCosDipAngle*mCosPhase;
  131. else
  132. return cos(mPhase + s*mH*mCurvature*mCosDipAngle)*mH*mCosDipAngle;
  133. }
  134. inline double MpdHelix::cz(double s) const
  135. {
  136. return mSinDipAngle;
  137. }
  138. inline const TVector3 MpdHelix::origin() const {return mOrigin;}
  139. inline TVector3 MpdHelix::at(double s) const
  140. {
  141. return TVector3(x(s), y(s), z(s));
  142. }
  143. inline TVector3 MpdHelix::cat(double s) const
  144. {
  145. return TVector3(cx(s), cy(s), cz(s));
  146. }
  147. inline double MpdHelix::pathLength(double xx, double yy) const
  148. {
  149. return fudgePathLength(TVector3(xx, yy, 0));
  150. }
  151. inline int MpdHelix::bad(double WorldSize) const
  152. {
  153. int ierr;
  154. if (!::finite(mDipAngle )) return 11;
  155. if (!::finite(mCurvature )) return 12;
  156. // ierr = mOrigin.bad(WorldSize);
  157. ierr = 1;
  158. if (ierr) return 3+ierr*100;
  159. if (::fabs(mDipAngle) >1.58) return 21;
  160. double qwe = ::fabs(::fabs(mDipAngle)-M_PI/2);
  161. if (qwe < 1./WorldSize ) return 31;
  162. if (::fabs(mCurvature) > WorldSize) return 22;
  163. if (mCurvature < 0 ) return 32;
  164. if (abs(mH) != 1 ) return 24;
  165. return 0;
  166. }
  167. #endif