PWG2/FEMTOSCOPY/AliFemto/Infrastructure/AliFmPhysicalHelix.cxx

   1 /***************************************************************************
   2  *
   3  * $Id$
   4  *
   5  * Author: Brian Lasiuk, Sep 1997
   6  ***************************************************************************
   7  *
   8  * Description:
   9  * Parametrization of a physical helix.
  10  *
  11  ***************************************************************************
  12  *
  13  * $Log$
  14  * Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
  15  * First version on CVS
  16  *
  17  * Revision 1.7  2005/07/06 18:49:56  fisyak
  18  * Replace StHelixD, StLorentzVectorD,StLorentzVectorF,StMatrixD,StMatrixF,AliFmPhysicalHelixD,StThreeVectorD,StThreeVectorF by templated version
  19  *
  20  * Revision 1.6  2003/09/02 17:59:35  perev
  21  * gcc 3.2 updates + WarnOff
  22  *
  23  * Revision 1.5  2002/06/21 17:49:26  genevb
  24  * Some minor speed improvements
  25  *
  26  * Revision 1.4  2002/02/20 00:56:23  ullrich
  27  * Added methods to calculate signed DCA.
  28  *
  29  * Revision 1.3  1999/02/24 11:42:18  ullrich
  30  * Fixed bug in momentum().
  31  *
  32  * Revision 1.2  1999/02/12 01:01:04  wenaus
  33  * Fix bug in momentum calculation
  34  *
  35  * Revision 1.1  1999/01/30 03:59:04  fisyak
  36  * Root Version of StarClassLibrary
  37  *
  38  * Revision 1.1  1999/01/23 00:29:21  ullrich
  39  * Initial Revision
  40  *
  41  **************************************************************************/
  42 #include <math.h>
  43 #include "AliFmHelix.h"
  44 #include "AliFmPhysicalHelix.h"
  45 #include "PhysicalConstants.h"
  46 #include "SystemOfUnits.h"
  47 #ifdef __ROOT__
  48 ClassImpT(AliFmPhysicalHelix,double);
  49 #endif
  50 AliFmPhysicalHelix::AliFmPhysicalHelix(){}
  51
  52 AliFmPhysicalHelix::~AliFmPhysicalHelix() { /* nop */ }
  53
  54 AliFmPhysicalHelix::AliFmPhysicalHelix(const AliFmThreeVector<double>& p,
  55                                        const AliFmThreeVector<double>& o,
  56                                        double B, double q)
  57 {
  58   mH = (q*B <= 0) ? 1 : -1;
  59   if(p.y() == 0 && p.x() == 0)
  60     setPhase((M_PI/4)*(1-2.*mH));
  61   else
  62     setPhase(atan2(p.y(),p.x())-mH*M_PI/2);
  63   setDipAngle(atan2(p.z(),p.perp()));
  64   mOrigin = o;
  65
  66 #ifndef ST_NO_NAMESPACES
  67   {
  68     using namespace units;
  69 #endif
  70     setCurvature(fabs((c_light*nanosecond/meter*q*B/tesla)/
  71                       (abs(p.mag())/GeV*mCosDipAngle)/meter));
  72 #ifndef ST_NO_NAMESPACES
  73   }
  74 #endif
  75 }
  76
  77 AliFmPhysicalHelix::AliFmPhysicalHelix(double c, double d, double phase,
  78                                  const AliFmThreeVector<double>& o, int h)
  79     : AliFmHelix(c, d, phase, o, h) { /* nop */}
  80
  81
  82 AliFmThreeVector<double> AliFmPhysicalHelix::momentum(double B) const
  83 {
  84     if (mSingularity)
  85         return(AliFmThreeVector<double>(0,0,0));
  86     else {
  87 #ifndef ST_NO_NAMESPACES
  88         {
  89             using namespace units;
  90 #endif
  91             double pt = GeV*fabs(c_light*nanosecond/meter*B/tesla)/(fabs(mCurvature)*meter);
  92
  93             return (AliFmThreeVector<double>(pt*cos(mPhase+mH*M_PI/2),   // pos part pos field
  94                                           pt*sin(mPhase+mH*M_PI/2),
  95                                           pt*tan(mDipAngle)));
  96 #ifndef ST_NO_NAMESPACES
  97         }
  98 #endif
  99     }
 100 }
 101
 102 AliFmThreeVector<double> AliFmPhysicalHelix::momentumAt(double S, double B) const
 103 {
 104     // Obtain phase-shifted momentum from phase-shift of origin
 105     double xc = this->xcenter();
 106     double yc = this->ycenter();
 107     double rx = (y(S)-yc)/(mOrigin.y()-yc);
 108     double ry = (x(S)-xc)/(mOrigin.x()-xc);
 109     return (this->momentum(B)).pseudoProduct(rx,ry,1.0);
 110 }
 111
 112 int AliFmPhysicalHelix::charge(double B) const
 113 {
 114     return (B > 0 ? -mH : mH);
 115 }
 116
 117 double AliFmPhysicalHelix::geometricSignedDistance(double x, double y)
 118 {
 119     // Geometric signed distance
 120     double thePath = this->pathLength(x,y);
 121     AliFmThreeVector<double> DCA2dPosition = this->at(thePath);
 122     DCA2dPosition.setZ(0);
 123     AliFmThreeVector<double> position(x,y,0);
 124     AliFmThreeVector<double> DCAVec = (DCA2dPosition-position);
 125     AliFmThreeVector<double> momVec;
 126     // Deal with straight tracks
 127     if (this->mSingularity) {
 128         momVec = this->at(1)- this->at(0);
 129         momVec.setZ(0);
 130     }
 131     else {
 132         momVec = this->momentumAt(thePath,1./tesla); // Don't care about Bmag.  Helicity is what matters.
 133         momVec.setZ(0);
 134     }
 135
 136     double cross = DCAVec.x()*momVec.y() - DCAVec.y()*momVec.x();
 137     double theSign = (cross>=0) ? 1. : -1.;
 138     return theSign*DCAVec.perp();
 139 }
 140
 141 double AliFmPhysicalHelix::curvatureSignedDistance(double x, double y)
 142 {
 143     // Protect against mH = 0 or zero field
 144     if (this->mSingularity || abs(this->mH)<=0) {
 145         return (this->geometricSignedDistance(x,y));
 146     }
 147     else {
 148         return (this->geometricSignedDistance(x,y))/(this->mH);
 149     }
 150
 151 }
 152
 153 double AliFmPhysicalHelix::geometricSignedDistance(const AliFmThreeVector<double>& pos)
 154 {
 155     double sdca2d = this->geometricSignedDistance(pos.x(),pos.y());
 156     double theSign = (sdca2d>=0) ? 1. : -1.;
 157     return (this->distance(pos))*theSign;
 158 }
 159
 160 double AliFmPhysicalHelix::curvatureSignedDistance(const AliFmThreeVector<double>& pos)
 161 {
 162     double sdca2d = this->curvatureSignedDistance(pos.x(),pos.y());
 163     double theSign = (sdca2d>=0) ? 1. : -1.;
 164     return (this->distance(pos))*theSign;
 165 }