RungeKutta/src/CylindricalState.h

0001 #ifndef CylindricalState_H
0002 #define CylindricalState_H
0003
0004 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0005 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0006 #include "FWCore/Utilities/interface/Visibility.h"
0007 #include "RKSmallVector.h"
0008
0009 #include <iostream>
0010
0011 /**
0012 State for solving the equation of motion with radius (in cylindrical coordinates) as free variable.
0013 The dependent variables are
0014   phi     - azimuthal angle
0015   z       - z coordinate
0016   dphi/dr - derivative of phi versus r
0017   dz/dr   - derivative of z versus r
0018   q/p     - charge over momentum magnitude
0019
0020 The coordinate system is externally defined
0021 */
0022
0023 class dso_internal CylindricalState {
0024 public:
0025   typedef double Scalar;
0026   typedef RKSmallVector<Scalar, 5> Vector;
0027
0028   CylindricalState() {}
0029
0030   CylindricalState(const LocalPoint& pos, const LocalVector& mom, Scalar ch) {
0031     rho_ = pos.perp();
0032     Scalar cosphi = pos.x() / rho_;
0033     Scalar sinphi = pos.y() / rho_;
0034     Scalar p_rho = mom.x() * cosphi + mom.y() * sinphi;
0035     Scalar p_phi = -mom.x() * sinphi + mom.y() * cosphi;
0036
0037     par_(0) = pos.phi();
0038     par_(1) = pos.z();
0039     par_(2) = p_phi / (p_rho * rho_);
0040     par_(3) = mom.z() / p_rho;
0041     par_(4) = ch / mom.mag();
0042
0043     prSign_ = p_rho > 0 ? 1.0 : -1.0;
0044
0045     std::cout << "CylindricalState built from pos " << pos << " mom " << mom << " charge " << ch << std::endl;
0046     std::cout << "p_rho " << p_rho << " p_phi " << p_phi << " dphi_drho " << par_(2) << std::endl;
0047     std::cout << "Which results in                " << position() << " mom " << momentum() << " charge " << charge()
0048               << std::endl;
0049   }
0050
0051   CylindricalState(Scalar rho, const Vector& par, Scalar prSign) : par_(par), rho_(rho), prSign_(prSign) {}
0052
0053   const LocalPoint position() const { return LocalPoint(LocalPoint::Cylindrical(rho_, par_(0), par_(1))); }
0054
0055   const LocalVector momentum() const {
0056     Scalar cosphi = cos(par_(0));
0057     Scalar sinphi = sin(par_(0));
0058     Scalar Q = sqrt(1 + rho_ * rho_ * par_(2) * par_(2) + par_(3) * par_(3));
0059     Scalar P = std::abs(1. / par_(4));
0060     Scalar p_rho = prSign_ * P / Q;
0061     Scalar p_phi = rho_ * par_(2) * p_rho;
0062     Scalar p_z = par_(3) * p_rho;
0063     LocalVector result(p_rho * cosphi - p_phi * sinphi, p_rho * sinphi + p_phi * cosphi, p_z);
0064     return result;
0065   }
0066
0067   const Vector& parameters() const { return par_; }
0068
0069   Scalar charge() const { return par_(4) > 0 ? 1 : -1; }
0070
0071   Scalar rho() const { return rho_; }
0072
0073   double prSign() const { return prSign_; }
0074
0075 private:
0076   Vector par_;
0077   Scalar rho_;
0078   Scalar prSign_;  ///< sign of local p_r
0079 };
0080
0081 #endif