HF/D0_helix/StPhysicalHelix.cxx

0001 /***************************************************************************
0002  *
0003  * $Id: StPhysicalHelix.cc,v 1.8 2007/08/20 23:24:54 perev Exp $
0004  *
0005  * Author: Brian Lasiuk, Sep 1997
0006  ***************************************************************************
0007  *
0008  * Description:
0009  * Parametrization of a physical helix.
0010  *
0011  ***************************************************************************
0012  *
0013  * $Log: StPhysicalHelix.cc,v $
0014  * Revision 1.8  2007/08/20 23:24:54  perev
0015  * BugFix #1027
0016  *
0017  * Revision 1.7  2005/07/06 18:49:56  fisyak
0018  * Replace StHelixD, StLorentzVectorD,StLorentzVectorF,StMatrixD,StMatrixF,StPhysicalHelixD,StThreeVectorD,StThreeVectorF by templated version
0019  *
0020  * Revision 1.6  2003/09/02 17:59:35  perev
0021  * gcc 3.2 updates + WarnOff
0022  *
0023  * Revision 1.5  2002/06/21 17:49:26  genevb
0024  * Some minor speed improvements
0025  *
0026  * Revision 1.4  2002/02/20 00:56:23  ullrich
0027  * Added methods to calculate signed DCA.
0028  *
0029  * Revision 1.3  1999/02/24 11:42:18  ullrich
0030  * Fixed bug in momentum().
0031  *
0032  * Revision 1.2  1999/02/12 01:01:04  wenaus
0033  * Fix bug in momentum calculation
0034  *
0035  * Revision 1.1  1999/01/30 03:59:04  fisyak
0036  * Root Version of StarClassLibrary
0037  *
0038  * Revision 1.1  1999/01/23 00:29:21  ullrich
0039  * Initial Revision
0040  *
0041  **************************************************************************/
0042 #include <math.h>
0043 #include "StHelix.h"
0044 #include "StPhysicalHelix.h"
0045 #include "PhysicalConstants.h"
0046 #include "SystemOfUnits.h"
0047 #ifdef __ROOT__
0048 ClassImpT(StPhysicalHelix,double);
0049 #endif
0050 StPhysicalHelix::StPhysicalHelix(){}
0051
0052 StPhysicalHelix::~StPhysicalHelix() { /* nop */ }
0053
0054 StPhysicalHelix::StPhysicalHelix(const TVector3& p,
0055                  const TVector3& o,
0056                  double B, double q)
0057 {
0058     mH = (q*B <= 0) ? 1 : -1;
0059     if(p.y() == 0 && p.x() == 0)
0060     setPhase((M_PI/4)*(1-2.*mH));
0061     else
0062     setPhase(atan2(p.y(),p.x())-mH*M_PI/2);
0063     setDipAngle(atan2(p.z(),p.Perp()));
0064     mOrigin = o;
0065
0066 #ifndef ST_NO_NAMESPACES
0067     {
0068     using namespace units;
0069 #endif
0070     setCurvature(fabs((c_light*nanosecond/meter*q*B/tesla)/
0071               (p.Mag()/GeV*mCosDipAngle)/meter));
0072 #ifndef ST_NO_NAMESPACES
0073     }
0074 #endif
0075 }
0076
0077 StPhysicalHelix::StPhysicalHelix(double c, double d, double phase,
0078                  const TVector3& o, int h)
0079     : StHelix(c, d, phase, o, h) { /* nop */}
0080
0081
0082 TVector3 StPhysicalHelix::momentum(double B) const
0083 {
0084     if (mSingularity)
0085     return(TVector3(0,0,0));
0086     else {
0087 #ifndef ST_NO_NAMESPACES
0088     {
0089         using namespace units;
0090 #endif
0091         double pt = GeV*fabs(c_light*nanosecond/meter*B/tesla)/(fabs(mCurvature)*meter);
0092
0093         return (TVector3(pt*cos(mPhase+mH*M_PI/2),   // pos part pos field
0094                       pt*sin(mPhase+mH*M_PI/2),
0095                       pt*tan(mDipAngle)));
0096 #ifndef ST_NO_NAMESPACES
0097     }
0098 #endif
0099     }
0100 }
0101
0102 TVector3 StPhysicalHelix::momentumAt(double S, double B) const
0103 {
0104     // Obtain phase-shifted momentum from phase-shift of origin
0105     StPhysicalHelix tmp(*this);
0106     tmp.moveOrigin(S);
0107     return tmp.momentum(B);
0108 }
0109
0110 int StPhysicalHelix::charge(double B) const
0111 {
0112     return (B > 0 ? -mH : mH);
0113 }
0114
0115 double StPhysicalHelix::geometricSignedDistance(double x, double y)
0116 {
0117     // Geometric signed distance
0118     double thePath = this->pathLength(x,y);
0119     TVector3 DCA2dPosition = this->at(thePath);
0120     DCA2dPosition.SetZ(0);
0121     TVector3 position(x,y,0);
0122     TVector3 DCAVec = (DCA2dPosition-position);
0123     TVector3 momVec;
0124     // Deal with straight tracks
0125     if (this->mSingularity) {
0126     momVec = this->at(1)- this->at(0);
0127     momVec.SetZ(0);
0128     }
0129     else {
0130     momVec = this->momentumAt(thePath,1./tesla); // Don't care about Bmag.  Helicity is what matters.
0131     momVec.SetZ(0);
0132     }
0133
0134     double cross = DCAVec.x()*momVec.y() - DCAVec.y()*momVec.x();
0135     double theSign = (cross>=0) ? 1. : -1.;
0136     return theSign*DCAVec.Perp();
0137 }
0138
0139 double StPhysicalHelix::curvatureSignedDistance(double x, double y)
0140 {
0141     // Protect against mH = 0 or zero field
0142     if (this->mSingularity || abs(this->mH)<=0) {
0143     return (this->geometricSignedDistance(x,y));
0144     }
0145     else {
0146     return (this->geometricSignedDistance(x,y))/(this->mH);
0147     }
0148
0149 }
0150
0151 double StPhysicalHelix::geometricSignedDistance(const TVector3& pos)
0152 {
0153     double sdca2d = this->geometricSignedDistance(pos.x(),pos.y());
0154     double theSign = (sdca2d>=0) ? 1. : -1.;
0155     return (this->distance(pos))*theSign;
0156 }
0157
0158 double StPhysicalHelix::curvatureSignedDistance(const TVector3& pos)
0159 {
0160     double sdca2d = this->curvatureSignedDistance(pos.x(),pos.y());
0161     double theSign = (sdca2d>=0) ? 1. : -1.;
0162     return (this->distance(pos))*theSign;
0163 }