RungeKutta/src/CurvilinearState.h

0001 #ifndef CurvilinearState_H
0002 #define CurvilinearState_H
0003
0004 #include "FWCore/Utilities/interface/Visibility.h"
0005 #include "DataFormats/GeometryVector/interface/Basic2DVector.h"
0006 #include "DataFormats/GeometryVector/interface/Basic3DVector.h"
0007 #include "VectorDoublet.h"
0008
0009 /**
0010 State for solving the equation of motion with Z as free variable.
0011 The dependent variables are
0012   x      - x coordinate
0013   y      - y coordinate
0014   dx/dz  - derivative of x versus z
0015   dy/dz  - derivative of y versus z
0016   q/p    - charge over momentum magnitude
0017
0018 The coordinate system is externally defined
0019 */
0020
0021 class dso_internal CurvilinearState {
0022 public:
0023   typedef double Scalar;
0024   typedef Basic2DVector<Scalar> Vector2D;
0025   typedef Basic3DVector<Scalar> Vector3D;
0026   typedef VectorDoublet<Vector2D, Vector3D> Vector;
0027
0028   CurvilinearState() {}
0029
0030   CurvilinearState(const Vector& v, Scalar z, Scalar pzsign) : par_(v), z_(z), pzSign_(pzsign) {}
0031
0032   CurvilinearState(const Vector3D& pos, const Vector3D& p, Scalar ch)
0033       : par_(Vector2D(pos.x(), pos.y()), Vector3D(p.x() / p.z(), p.y() / p.z(), ch / p.mag())),
0034         z_(pos.z()),
0035         pzSign_(p.z() > 0. ? 1. : -1.) {}
0036
0037   const Vector3D position() const { return Vector3D(par_.first().x(), par_.first().y(), z_); }
0038
0039   const Vector3D momentum() const {
0040     Scalar p = 1. / fabs(par_.second().z());
0041     if (p > 1.e9)
0042       p = 1.e9;
0043     Scalar dxdz = par_.second().x();
0044     Scalar dydz = par_.second().y();
0045     Scalar dz = pzSign_ / sqrt(1. + dxdz * dxdz + dydz * dydz);
0046     Scalar dx = dz * dxdz;
0047     Scalar dy = dz * dydz;
0048     return Vector3D(dx * p, dy * p, dz * p);
0049   }
0050
0051   const Vector& parameters() const { return par_; }
0052
0053   Scalar charge() const { return par_.second().z() > 0 ? 1 : -1; }
0054
0055   Scalar z() const { return z_; }
0056
0057   double pzSign() const { return pzSign_; }
0058
0059 private:
0060   Vector par_;
0061   Scalar z_;
0062   Scalar pzSign_;  ///< sign of local pz
0063 };
0064
0065 inline CurvilinearState operator+(const CurvilinearState& a, const CurvilinearState& b) {
0066   return CurvilinearState(a.parameters() + b.parameters(), a.z() + b.z(), a.pzSign());
0067 }
0068
0069 inline CurvilinearState operator-(const CurvilinearState& a, const CurvilinearState& b) {
0070   return CurvilinearState(a.parameters() - b.parameters(), a.z() - b.z(), a.pzSign());
0071 }
0072
0073 inline CurvilinearState operator*(const CurvilinearState& v, const CurvilinearState::Scalar& s) {
0074   return CurvilinearState(v.parameters() * s, v.z() * s, v.pzSign());
0075 }
0076 inline CurvilinearState operator*(const CurvilinearState::Scalar& s, const CurvilinearState& v) {
0077   return CurvilinearState(v.parameters() * s, v.z() * s, v.pzSign());
0078 }
0079
0080 inline CurvilinearState operator/(const CurvilinearState& v, const CurvilinearState::Scalar& s) {
0081   return CurvilinearState(v.parameters() / s, v.z() / s, v.pzSign());
0082 }
0083
0084 #endif