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