RungeKutta/src/CartesianLorentzForce.h

0001 #ifndef CartesianLorentzForce_H
0002 #define CartesianLorentzForce_H
0003
0004 #include "FWCore/Utilities/interface/Visibility.h"
0005 #include "RKDerivative.h"
0006 #include "RKLocalFieldProvider.h"
0007
0008 /// Derivative calculation for the 6D cartesian case.
0009
0010 class dso_internal CartesianLorentzForce final : public RKDerivative<double, 6> {
0011 public:
0012   typedef RKDerivative<double, 6> Base;
0013   typedef Base::Scalar Scalar;
0014   typedef Base::Vector Vector;
0015
0016   CartesianLorentzForce(const RKLocalFieldProvider& field, float ch) : theField(field), theCharge(ch) {}
0017
0018   Vector operator()(Scalar z, const Vector& state) const override;
0019
0020 private:
0021   const RKLocalFieldProvider& theField;
0022   float theCharge;
0023 };
0024
0025 #include "CartesianStateAdaptor.h"
0026 inline CartesianLorentzForce::Vector CartesianLorentzForce::operator()(Scalar z, const Vector& state) const {
0027   // derivatives in case S is the free parameter
0028   CartesianStateAdaptor start(state);
0029   auto bfield = theField.inTesla(RKLocalFieldProvider::LocalPoint(start.position()));
0030   constexpr float k = 2.99792458e-3;  // conversion to [cm]
0031
0032   /// Derivative d(pos)/ds is simply normalized momentum
0033   auto dpos = start.momentum().unit();
0034
0035   /// Lorentz force in absence of electric field
0036   auto dmom = (k * theCharge) * dpos.cross(bfield);
0037
0038   return CartesianStateAdaptor::rkstate(dpos, dmom);
0039 }
0040
0041 #endif