RungeKutta/src/CurvilinearLorentzForce.icc

0001 #include "CurvilinearLorentzForce.h"
0002 #include "CurvilinearState.h"
0003 #include "RKLocalFieldProvider.h"
0004
0005 #include <exception>
0006
0007 class CurvilinearLorentzForceException : public std::exception {
0008 public:
0009   CurvilinearLorentzForceException() throw() {}
0010   ~CurvilinearLorentzForceException() throw() override {}
0011 };
0012
0013 template <typename T, int N>
0014 typename CurvilinearLorentzForce<T, N>::Vector CurvilinearLorentzForce<T, N>::operator()(Scalar z,
0015                                                                                          const Vector& state) const {
0016   // derivatives in case Z is the free parameter
0017
0018   RKLocalFieldProvider::Vector B = theField.inTesla(state[0], state[1], z);
0019   double k = 2.99792458e-3;  // conversion to [cm]
0020
0021   //    cout << "CurvilinearLorentzForce: field at (" << state[0] << ", "
0022   //         <<  state[1] << ", " << z << ") is " << B << endl;
0023
0024   /// this should be changed when the ref frame of the solution is not the global one!
0025   //CurvilinearState::Vector3D B( MagneticField::inTesla(GlobalPoint(start.position())));
0026
0027   double dxdz = state[2];
0028   double dydz = state[3];
0029   double lambda = state[4];
0030
0031   double R = sqrt(1.0 + dxdz * dxdz + dydz * dydz);
0032
0033   // cout << "CurvilinearLorentzForce: R " << R << " dxdz " << dxdz << " dydz " << dydz << endl;
0034   if (R > 10) {
0035     throw CurvilinearLorentzForceException();
0036   }
0037
0038   double d2x_dz2 = k * lambda * R * (dxdz * dydz * B.x() - (1.0 + dxdz * dxdz) * B.y() + dydz * B.z());
0039   double d2y_dz2 = k * lambda * R * ((1 + dydz * dydz) * B.x() - dxdz * dydz * B.y() - dxdz * B.z());
0040
0041   // the derivative of q/p is zero -- momentum conservation if no energy loss
0042   Vector res;
0043   res[0] = dxdz;
0044   res[1] = dydz;
0045   res[2] = d2x_dz2;
0046   res[3] = d2y_dz2;
0047   res[4] = 0;
0048
0049   // cout << "CurvilinearLorentzForce: derivatives are " << res;
0050
0051   return res;
0052 }