PWG2/FEMTOSCOPY/AliFemto/AliFmPhysicalHelix.cxx

   1 ///////////////////////////////////////////////////////////////////////////
   2 //                                                                       //
   3 // AliFmHelix: a helper helix class                                      //
   4 // Includes all the operations and specifications of the helix. Can be   //
   5 // used to determine path lengths, distance of closest approach etc.     //
   6 //                                                                       //
   7 ///////////////////////////////////////////////////////////////////////////
   8 #include <math.h>
   9 #include "AliFmHelix.h"
  10 #include "AliFmPhysicalHelix.h"
  11 #include "PhysicalConstants.h"
  12 #include "SystemOfUnits.h"
  13 #ifdef __ROOT__
  14 ClassImpT(AliFmPhysicalHelix,double);
  15 #endif
  16 AliFmPhysicalHelix::AliFmPhysicalHelix(){}
  17
  18 AliFmPhysicalHelix::~AliFmPhysicalHelix() { /* nop */ }
  19
  20 AliFmPhysicalHelix::AliFmPhysicalHelix(const AliFmThreeVector<double>& p,
  21                                        const AliFmThreeVector<double>& o,
  22                                        double B, double q)
  23 {
  24   // Constructor from given parameters
  25   fH = (q*B <= 0) ? 1 : -1;
  26   if(p.y() == 0 && p.x() == 0)
  27     SetPhase((M_PI/4)*(1-2.*fH));
  28   else
  29     SetPhase(atan2(p.y(),p.x())-fH*M_PI/2);
  30   SetDipAngle(atan2(p.z(),p.perp()));
  31   fOrigin = o;
  32
  33 #ifndef ST_NO_NAMESPACES
  34   {
  35     using namespace units;
  36 #endif
  37     SetCurvature(fabs((c_light*nanosecond/meter*q*B/tesla)/
  38                       (abs(p.mag())/GeV*fCosDipAngle)/meter));
  39 #ifndef ST_NO_NAMESPACES
  40   }
  41 #endif
  42 }
  43
  44 AliFmPhysicalHelix::AliFmPhysicalHelix(double c, double d, double phase,
  45                                        const AliFmThreeVector<double>& o, int h)
  46   : AliFmHelix(c, d, phase, o, h) { /* nop */}
  47
  48
  49 AliFmThreeVector<double> AliFmPhysicalHelix::Momentum(double B) const
  50 {
  51   // momentum for given magnetic field
  52   if (fSingularity)
  53     return(AliFmThreeVector<double>(0,0,0));
  54   else {
  55 #ifndef ST_NO_NAMESPACES
  56     {
  57             using namespace units;
  58 #endif
  59             double pt = GeV*fabs(c_light*nanosecond/meter*B/tesla)/(fabs(fCurvature)*meter);
  60
  61             return (AliFmThreeVector<double>(pt*cos(fPhase+fH*M_PI/2),   // pos part pos field
  62                                           pt*sin(fPhase+fH*M_PI/2),
  63                                           pt*tan(fDipAngle)));
  64 #ifndef ST_NO_NAMESPACES
  65         }
  66 #endif
  67     }
  68 }
  69
  70 AliFmThreeVector<double> AliFmPhysicalHelix::MomentumAt(double S, double B) const
  71 {
  72     // Obtain phase-shifted momentum from phase-shift of origin
  73     double xc = this->XCenter();
  74     double yc = this->YCenter();
  75     double rx = (Y(S)-yc)/(fOrigin.y()-yc);
  76     double ry = (X(S)-xc)/(fOrigin.x()-xc);
  77     return (this->Momentum(B)).pseudoProduct(rx,ry,1.0);
  78 }
  79
  80 int AliFmPhysicalHelix::Charge(double B) const
  81 {
  82   // charge
  83     return (B > 0 ? -fH : fH);
  84 }
  85
  86 double AliFmPhysicalHelix::GeometricSignedDistance(double x, double y)
  87 {
  88     // Geometric signed distance
  89     double thePath = this->PathLength(x,y);
  90     AliFmThreeVector<double> tDCA2dPosition = this->At(thePath);
  91     tDCA2dPosition.setZ(0);
  92     AliFmThreeVector<double> position(x,y,0);
  93     AliFmThreeVector<double> tDCAVec = (tDCA2dPosition-position);
  94     AliFmThreeVector<double> momVec;
  95     // Deal with straight tracks
  96     if (this->fSingularity) {
  97         momVec = this->At(1)- this->At(0);
  98         momVec.setZ(0);
  99     }
 100     else {
 101         momVec = this->MomentumAt(thePath,1./tesla); // Don't care about Bmag.  Helicity is what matters.
 102         momVec.setZ(0);
 103     }
 104
 105     double cross = tDCAVec.x()*momVec.y() - tDCAVec.y()*momVec.x();
 106     double theSign = (cross>=0) ? 1. : -1.;
 107     return theSign*tDCAVec.perp();
 108 }
 109
 110 double AliFmPhysicalHelix::CurvatureSignedDistance(double x, double y)
 111 {
 112     // Protect against fH = 0 or zero field
 113     if (this->fSingularity || abs(this->fH)<=0) {
 114         return (this->GeometricSignedDistance(x,y));
 115     }
 116     else {
 117         return (this->GeometricSignedDistance(x,y))/(this->fH);
 118     }
 119
 120 }
 121
 122 double AliFmPhysicalHelix::GeometricSignedDistance(const AliFmThreeVector<double>& pos)
 123 {
 124   // Geometric distance
 125     double sdca2d = this->GeometricSignedDistance(pos.x(),pos.y());
 126     double theSign = (sdca2d>=0) ? 1. : -1.;
 127     return (this->Distance(pos))*theSign;
 128 }
 129
 130 double AliFmPhysicalHelix::CurvatureSignedDistance(const AliFmThreeVector<double>& pos)
 131 {
 132   // Distance with the sign dependent on curvature sigm
 133     double sdca2d = this->CurvatureSignedDistance(pos.x(),pos.y());
 134     double theSign = (sdca2d>=0) ? 1. : -1.;
 135     return (this->Distance(pos))*theSign;
 136 }