PWG2/FEMTOSCOPY/AliFemto/AliFmHelix.h

   1 ///////////////////////////////////////////////////////////////////////////
   2 //                                                                       //
   3 // AliFmHelix: a helper helix class                                      //
   4 //                                                                       //
   5 ///////////////////////////////////////////////////////////////////////////
   6
   7 #ifndef ALIFMHELIX_H
   8 #define ALIFMHELIX_H
   9
  10 #include <math.h>
  11 #include <utility>
  12 #include <algorithm>
  13 #include "AliFmThreeVector.h"
  14 #if !defined(ST_NO_NAMESPACES)
  15 using std::pair;
  16 using std::swap;
  17 using std::max;
  18 #endif
  19
  20 #ifdef WIN32
  21 #include "gcc2vs.h"
  22 #endif
  23
  24 class AliFmHelix {
  25 public:
  26     /// curvature, dip angle, phase, origin, h
  27     AliFmHelix(double c, double dip, double phase,
  28                const AliFmThreeVector<double>& o, int h=-1);
  29
  30     virtual ~AliFmHelix();
  31     // AliFmHelix(const AliFmHelix&);                   // use default
  32     // AliFmHelix& operator=(const AliFmHelix&);        // use default
  33
  34     double       DipAngle()   const;
  35     double       Curvature()  const;    /// 1/R in xy-plane
  36     double       Phase()      const;    /// aziumth in xy-plane measured from ring center
  37     double       XCenter()    const;    /// x-center of circle in xy-plane
  38     double       YCenter()    const;    /// y-center of circle in xy-plane
  39     int          H()          const;    /// -sign(q*B);
  40
  41     const AliFmThreeVector<double>& Origin() const;     /// starting point
  42
  43     void SetParameters(double c, double dip, double phase, const AliFmThreeVector<double>& o, int h);
  44
  45     double       X(double s)  const;
  46     double       Y(double s)  const;
  47     double       Z(double s)  const;
  48
  49     AliFmThreeVector<double>  At(double s) const;
  50
  51     /// returns period length of helix
  52     double       Period()       const;
  53
  54     /// path length at given r (cylindrical r)
  55     pair<double, double> PathLength(double r)   const;
  56
  57     /// path length at given r (cylindrical r, cylinder axis at x,y)
  58     pair<double, double> PathLength(double r, double x, double y, bool scanPeriods = true);
  59
  60     /// path length at distance of closest approach to a given point
  61     double       PathLength(const AliFmThreeVector<double>& p, bool scanPeriods = true) const;
  62
  63     /// path length at intersection with plane
  64     double       PathLength(const AliFmThreeVector<double>& r,
  65                             const AliFmThreeVector<double>& n) const;
  66
  67     /// path length at distance of closest approach in the xy-plane to a given point
  68     double       PathLength(double x, double y) const;
  69
  70     /// path lengths at dca between two helices
  71     pair<double, double> PathLengths(const AliFmHelix& h, bool scanPeriods = true) const;
  72
  73     /// minimal distance between point and helix
  74     double       Distance(const AliFmThreeVector<double>& p, bool scanPeriods = true) const;
  75
  76     /// checks for valid parametrization
  77     bool         Valid(double world = 1.e+5) const {return !Bad(world);}
  78     int            Bad(double world = 1.e+5) const;
  79
  80     /// move the origin along the helix to s which becomes then s=0
  81     virtual void MoveOrigin(double s);
  82
  83     static const double fgkNoSolution;
  84
  85 protected:
  86     AliFmHelix();
  87
  88     void SetCurvature(double d);        /// performs also various checks
  89     void SetPhase(double d);
  90     void SetDipAngle(double d);
  91
  92     /// value of S where distance in x-y plane is minimal
  93     double FudgePathLength(const AliFmThreeVector<double>& v) const;
  94
  95 protected:
  96     bool                   fSingularity;        // true for straight line case (B=0)
  97     AliFmThreeVector<double>  fOrigin;          // Helix origin
  98     double                 fDipAngle;           // Helix dip angle
  99     double                 fCurvature;          // curvature
 100     double                 fPhase;              // phase
 101     int                    fH;                  // -sign(q*B);
 102
 103     double                 fCosDipAngle;        // cosine of the dip angle
 104     double                 fSinDipAngle;        // sine of the dip angle
 105     double                 fCosPhase;           // cosine of the phase
 106     double                 fSinPhase;           // sine of the phase
 107 #ifdef __ROOT__
 108   ClassDef(AliFmHelix,1)
 109 #endif
 110 };
 111
 112 //
 113 //     Non-member functions
 114 //
 115 int operator== (const AliFmHelix&, const AliFmHelix&);
 116 int operator!= (const AliFmHelix&, const AliFmHelix&);
 117 ostream& operator<<(ostream&, const AliFmHelix&);
 118
 119 //
 120 //     Inline functions
 121 //
 122 inline int AliFmHelix::H() const {return fH;}
 123
 124 inline double AliFmHelix::DipAngle() const {return fDipAngle;}
 125
 126 inline double AliFmHelix::Curvature() const {return fCurvature;}
 127
 128 inline double AliFmHelix::Phase() const {return fPhase;}
 129
 130 inline double AliFmHelix::X(double s) const
 131 {
 132     if (fSingularity)
 133         return fOrigin.x() - s*fCosDipAngle*fSinPhase;
 134     else
 135         return fOrigin.x() + (cos(fPhase + s*fH*fCurvature*fCosDipAngle)-fCosPhase)/fCurvature;
 136 }
 137
 138 inline double AliFmHelix::Y(double s) const
 139 {
 140     if (fSingularity)
 141         return fOrigin.y() + s*fCosDipAngle*fCosPhase;
 142     else
 143         return fOrigin.y() + (sin(fPhase + s*fH*fCurvature*fCosDipAngle)-fSinPhase)/fCurvature;
 144 }
 145
 146 inline double AliFmHelix::Z(double s) const
 147 {
 148     return fOrigin.z() + s*fSinDipAngle;
 149 }
 150
 151 inline const AliFmThreeVector<double>& AliFmHelix::Origin() const {return fOrigin;}
 152
 153 inline AliFmThreeVector<double> AliFmHelix::At(double s) const
 154 {
 155     return AliFmThreeVector<double>(X(s), Y(s), Z(s));
 156 }
 157
 158 inline double AliFmHelix::PathLength(double x, double y) const
 159 {
 160     return FudgePathLength(AliFmThreeVector<double>(x, y, 0));
 161 }
 162 inline int AliFmHelix::Bad(double WorldSize) const
 163 {
 164
 165     int ierr;
 166     if (!::finite(fDipAngle    ))       return   11;
 167     if (!::finite(fCurvature   ))       return   12;
 168
 169     ierr = fOrigin.bad(WorldSize);
 170     if (ierr)                           return    3+ierr*100;
 171
 172     if (::fabs(fDipAngle)  >1.58)       return   21;
 173     double qwe = ::fabs(::fabs(fDipAngle)-M_PI/2);
 174     if (qwe < 1./WorldSize      )       return   31;
 175
 176     if (::fabs(fCurvature) > WorldSize) return   22;
 177     if (fCurvature < 0          )       return   32;
 178
 179     if (abs(fH) != 1            )       return   24;
 180
 181     return 0;
 182 }
 183
 184 #endif