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