#ifndef THelix3D_hh #define THelix3D_hh class TTrajectory3D; #include "Angle.hh" #include "TTrajectory3D.hh" class TTrajectoryPoint; class THelix3D : public TTrajectory3D { public: // ****** Constructors and destructor THelix3D(); // Construct from particle momentum, position, // field, and charge. THelix3D(const TVector3 &, const TVector3 &, double q, double field); // Copy Constructor THelix3D(const THelix3D &right); // Construct on five helix parameters, // phi0 should be in the range [0,2pi) THelix3D(double cotTheta, double curvature, double z0, double d0, Angle phi0); virtual ~THelix3D(); // ****** init methods void Init(TTrajectoryPoint* point, double q, double field); // phi0 in [0,2pi) void Init(double cotTheta, double curvature, double z0, double d0, double phi0); // ****** modifiers void setCotTheta(double cotTheta) { _cotTheta=cotTheta; _isStale=true; } void setCurvature(double curvature) { _curvature=curvature; _isStale=true; } void setZ0(double z0) { _z0 = z0; _isStale=true; } void setD0(double d0) { _d0=d0; _isStale=true; } // phi0 [0,2pi] void setPhi0(Angle phi0) { _phi0=phi0; // Range automatically handled by Angle // while(_phi0<0) _phi0 += 2*M_PI; // while(_phi0>2*M_PI) _phi0 -= 2*M_PI; _isStale=true; } // Assignment operator const THelix3D & operator=(const THelix3D &right); // Relational operators. bool operator==(const THelix3D & right) const; bool operator!=(const THelix3D & right) const; // Get Position as a function of // (three-dimensional) path length virtual void GetPosition(TVector3& pos, double s = 0.0) const; // Get Direction as a function of // (three-dimensional) path length virtual void GetDirection(TVector3& dir, double s = 0.0) const; // Get the second derivative of the // helix vs (three-dimensional) path // length virtual TVector3 getSecondDerivative(double s = 0.0) const; // Get both position and direction // at once. virtual void GetPoint (TTrajectoryPoint* point, double s = 0.0) const; virtual void GetLocation(TTrajectory3D::Location & loc, double s = 0.0) const; // Get pathlength at fixed // rho=sqrt(x^2 + y^2) virtual double GetPathLengthAtRho(double rho) const; virtual double GetPathLengthAtZ (double z) const; /////////////////////////////////////////////////////////////////////////////// // KCDF: analytical computation of helix/plane intersection. // // What we really compute is the intersection of a line and // a circle (projected helix) in the x-y plane. // // >>>>>>>>>> W A R N I N G W A R N I N G W A R N I N G <<<<<<<<<< // > < // > We assume the plane to be parallel or perpendicular < // > to the z axis (i.e. the B-field), < // > since otherwise there is no analytical solution. (One would end up < // > with an equation of type cos(alpha) == alpha.) < // > Although we know this assumption doesn´t hold exactly, we think it < // > is a reasonable first approximation. < // > In cases when our assumption has to be dropped, one can use the < // > intersection point computed here as a *good* starting point of a < // > numerical method, or one uses another analytical method to compute < // > the intersection of the tangent line at the point and the plane. < // > We plan to use one of these approaches in the near future, but < // > this is NOT YET IMPLEMENTED! < // > For the time being, we invoke the old numerical < // > TTrajectory3D::newIntersectionWith in such circumstances. < // > < // >>>>>>>>>> W A R N I N G W A R N I N G W A R N I N G <<<<<<<<<< // // Kurt Rinnert, 08/31/1998 ////////////////////////////////////////////////////////////////////////////// // Location* newIntersectionWith(const HepPlane3D& plane) const; // Get certain parameters as a function // of two-dimensional R. Angle getPhiAtR(double r) const; double getZAtR(double r) const; double getL2DAtR(double r) const; double getCosAlphaAtR(double r) const; // Get signed 1/radius double getInverseRadius() const; // Get unsigned radius double getRadius() const; // Get the turning angle as a // function of path length SignedAngle getTurningAngle(double s) const; // Get the Curvature double getCurvature() const; // Get helicity, positive for a // counterclockwise helix double getHelicity() const; double getCotTheta() const; Angle getPhi0() const; double getD0() const; double getZ0() const; // Get sign of the z component of // angular momentum about origin double getSignLz() const; // Get sines and cosines of Phi0 // and Theta double getSinPhi0() const; double getCosPhi0() const; double getSinTheta() const; double getCosTheta() const; // Set Parameters. void setParameters(const double p[]) { setCotTheta(p[0]); setCurvature(p[1]); setZ0(p[2]); setD0(p[3]); setPhi0(p[4]); _isStale=true; } // Get the parameters as a vector. void getParameters(double p[]) const { p[0]=getCotTheta(); p[1]=getCurvature(); p[2]=getZ0(); p[3]=getD0(); p[4]=getPhi0(); } // Create a helix from a vector static THelix3D create(const double v[]) { return THelix3D(v[0],v[1],v[2],v[3],v[4]); } // Return size of the parameter space // (=5) static unsigned int getParameterSpaceSize() { return 5; } private: // This is the THelix3D: double _cotTheta; double _curvature; double _z0; double _d0; Angle _phi0; // This is the cache mutable bool _isStale; // ! mutable double _sinPhi0; // ! mutable double _cosPhi0; // ! mutable double _sinTheta; // ! mutable double _cosTheta; // ! mutable double _s; // ! mutable double _aa; // ! mutable double _ss; // ! mutable double _cc; // ! double _vParameters[5]; // ! mutable bool _centerIsValid; // ! needed by newIntersectionWith KR mutable double _m_x; // ! mutable double _m_y; // ! //needed whenever _sinPhi0, _cosPh0, // _sinTheta, or _cosTheta is used. void _refreshCache() const; // neede whenever _ss or _cc are used. void _cacheSinesAndCosines(double s) const; ClassDef(THelix3D,1) }; #endif