HF/D0_helix/StHelix.h

0001 /**
0002  * \class StHelix
0003  * \author Thomas Ullrich, Sep 26 1997
0004  *
0005  * Parametrization of a helix. Can also cope with straight tracks, i.e.
0006  * with zero curvature. This represents only the mathematical model of
0007  * a helix. See the SCL user guide for more.
0008  */
0009 /***************************************************************************
0010  *
0011  * $Id: StHelix.hh,v 1.14 2015/04/22 18:02:01 ullrich Exp $
0012  *
0013  * Author: Thomas Ullrich, Sep 1997
0014  ***************************************************************************
0015  *
0016  * Description: Parametrization of a helix
0017  *
0018  ***************************************************************************
0019  *
0020  * $Log: StHelix.hh,v $
0021  * Revision 1.14  2015/04/22 18:02:01  ullrich
0022  * Added two default argument to dca of two helices for HFT.
0023  *
0024  * Revision 1.13  2010/10/18 21:55:11  fisyak
0025  * Warn off for gcc4.5.1 64bits
0026  *
0027  * Revision 1.12  2006/06/29 15:53:22  ullrich
0028  * Added direction vector at given pathlength (written by Yuri).
0029  *
0030  * Revision 1.11  2005/10/13 22:23:27  genevb
0031  * NoSolution is public
0032  *
0033  * Revision 1.10  2005/07/06 18:49:56  fisyak
0034  * Replace StHelixD, StLorentzVectorD,StLorentzVectorF,StMatrixD,StMatrixF,StPhysicalHelixD,StThreeVectorD,StThreeVectorF by templated version
0035  *
0036  * Revision 1.9  2004/12/02 02:51:16  ullrich
0037  * Added option to pathLenghth() and distance() to search for
0038  * DCA only within one period. Default stays as it was.
0039  *
0040  * Revision 1.8  2003/10/30 20:06:46  perev
0041  * Check of quality added
0042  *
0043  * Revision 1.7  2002/06/21 17:49:25  genevb
0044  * Some minor speed improvements
0045  *
0046  * Revision 1.6  2002/04/24 02:41:55  ullrich
0047  * Restored old format.
0048  *
0049  **************************************************************************/
0050
0051 #ifndef ST_HELIX_HH
0052 #define ST_HELIX_HH
0053
0054 #include <math.h>
0055 #include <utility>
0056 #include <algorithm>
0057 #include <iostream>
0058 #include "TVector3.h"
0059 #include "SystemOfUnits.h"
0060 #if !defined(ST_NO_NAMESPACES)
0061 using std::pair;
0062 using std::swap;
0063 using std::max;
0064 #endif
0065 using namespace std;
0066
0067 class StHelix {
0068 public:
0069     /// curvature, dip angle, phase, origin, h
0070     StHelix(double c, double dip, double phase,
0071         const TVector3& o, int h=-1);
0072
0073     virtual ~StHelix();
0074     // StHelix(const StHelix&);         // use default
0075     // StHelix& operator=(const StHelix&);  // use default
0076
0077     double       dipAngle()   const;
0078     double       curvature()  const;    /// 1/R in xy-plane
0079     double       phase()      const;    /// aziumth in xy-plane measured from ring center
0080     double       xcenter()    const;    /// x-center of circle in xy-plane
0081     double       ycenter()    const;    /// y-center of circle in xy-plane
0082     int          h()          const;    /// -sign(q*B);
0083
0084     const TVector3& origin() const; /// starting point
0085
0086     void setParameters(double c, double dip, double phase, const TVector3& o, int h);
0087
0088     /// coordinates of helix at point s
0089     double       x(double s)  const;
0090     double       y(double s)  const;
0091     double       z(double s)  const;
0092
0093     TVector3  at(double s) const;
0094
0095     /// pointing vector of helix at point s
0096     double       cx(double s)  const;
0097     double       cy(double s)  const;
0098     double       cz(double s = 0)  const;
0099
0100     TVector3  cat(double s) const;
0101
0102     /// returns period length of helix
0103     double       period()       const;
0104
0105     /// path length at given r (cylindrical r)
0106     pair<double, double> pathLength(double r)   const;
0107
0108     /// path length at given r (cylindrical r, cylinder axis at x,y)
0109     pair<double, double> pathLength(double r, double x, double y);
0110
0111     /// path length at distance of closest approach to a given point
0112     double       pathLength(const TVector3& p, bool scanPeriods = true) const;
0113
0114     /// path length at intersection with plane
0115     double       pathLength(const TVector3& r,
0116                 const TVector3& n) const;
0117
0118     /// path length at distance of closest approach in the xy-plane to a given point
0119     double       pathLength(double x, double y) const;
0120
0121     /// path lengths at dca between two helices
0122     pair<double, double> pathLengths(const StHelix&,
0123                                      double minStepSize = 10*micrometer,
0124                                      double minRange = 10*centimeter) const;
0125
0126     /// minimal distance between point and helix
0127     double       distance(const TVector3& p, bool scanPeriods = true) const;
0128
0129     /// checks for valid parametrization
0130     bool         valid(double world = 1.e+5) const {return !bad(world);}
0131     int            bad(double world = 1.e+5) const;
0132
0133     /// move the origin along the helix to s which becomes then s=0
0134     virtual void moveOrigin(double s);
0135
0136     static const double NoSolution;
0137
0138 protected:
0139     StHelix();
0140
0141     void setCurvature(double);  /// performs also various checks
0142     void setPhase(double);
0143     void setDipAngle(double);
0144
0145     /// value of S where distance in x-y plane is minimal
0146     double fudgePathLength(const TVector3&) const;
0147
0148 protected:
0149     bool                   mSingularity;    // true for straight line case (B=0)
0150     TVector3               mOrigin;
0151     double                 mDipAngle;
0152     double                 mCurvature;
0153     double                 mPhase;
0154     int                    mH;          // -sign(q*B);
0155
0156     double                 mCosDipAngle;
0157     double                 mSinDipAngle;
0158     double                 mCosPhase;
0159     double                 mSinPhase;
0160
0161 #ifdef __ROOT__
0162   ClassDef(StHelix,1)
0163 #endif
0164 };
0165
0166 //
0167 //     Non-member functions
0168 //
0169 int operator== (const StHelix&, const StHelix&);
0170 int operator!= (const StHelix&, const StHelix&);
0171 ostream& operator<<(ostream&, const StHelix&);
0172
0173 //
0174 //     Inline functions
0175 //
0176 inline int StHelix::h() const {return mH;}
0177
0178 inline double StHelix::dipAngle() const {return mDipAngle;}
0179
0180 inline double StHelix::curvature() const {return mCurvature;}
0181
0182 inline double StHelix::phase() const {return mPhase;}
0183
0184 inline double StHelix::x(double s) const
0185 {
0186     if (mSingularity)
0187     return mOrigin.x() - s*mCosDipAngle*mSinPhase;
0188     else
0189     return mOrigin.x() + (cos(mPhase + s*mH*mCurvature*mCosDipAngle)-mCosPhase)/mCurvature;
0190 }
0191
0192 inline double StHelix::y(double s) const
0193 {
0194     if (mSingularity)
0195     return mOrigin.y() + s*mCosDipAngle*mCosPhase;
0196     else
0197     return mOrigin.y() + (sin(mPhase + s*mH*mCurvature*mCosDipAngle)-mSinPhase)/mCurvature;
0198 }
0199
0200 inline double StHelix::z(double s) const
0201 {
0202     return mOrigin.z() + s*mSinDipAngle;
0203 }
0204
0205 inline double StHelix::cx(double s)  const
0206 {
0207     if (mSingularity)
0208     return -mCosDipAngle*mSinPhase;
0209     else
0210     return -sin(mPhase + s*mH*mCurvature*mCosDipAngle)*mH*mCosDipAngle;
0211 }
0212
0213 inline double StHelix::cy(double s)  const
0214 {
0215     if (mSingularity)
0216     return mCosDipAngle*mCosPhase;
0217     else
0218     return cos(mPhase + s*mH*mCurvature*mCosDipAngle)*mH*mCosDipAngle;
0219 }
0220
0221 inline double StHelix::cz(double /* s */)  const
0222 {
0223     return mSinDipAngle;
0224 }
0225
0226 inline const TVector3& StHelix::origin() const {return mOrigin;}
0227
0228 inline TVector3 StHelix::at(double s) const
0229 {
0230     return TVector3(x(s), y(s), z(s));
0231 }
0232
0233 inline TVector3 StHelix::cat(double s) const
0234 {
0235     return TVector3(cx(s), cy(s), cz(s));
0236 }
0237
0238 inline double StHelix::pathLength(double X, double Y) const
0239 {
0240     return fudgePathLength(TVector3(X, Y, 0));
0241 }
0242 inline int StHelix::bad(double WorldSize) const
0243 {
0244
0245 //    int ierr;
0246     if (!::finite(mDipAngle    ))   return   11;
0247     if (!::finite(mCurvature   ))   return   12;
0248
0249 //    ierr = mOrigin.bad(WorldSize);
0250 //    if (ierr)                           return    3+ierr*100;
0251
0252     if (::fabs(mDipAngle)  >1.58)   return   21;
0253     double qwe = ::fabs(::fabs(mDipAngle)-M_PI/2);
0254     if (qwe < 1./WorldSize      )   return   31;
0255
0256     if (::fabs(mCurvature) > WorldSize) return   22;
0257     if (mCurvature < 0          )   return   32;
0258
0259     if (abs(mH) != 1            )       return   24;
0260
0261     return 0;
0262 }
0263
0264 #endif