RungeKutta/src/CylindricalLorentzForce.icc

0001 #include "CylindricalLorentzForce.h"
0002 #include "CylindricalState.h"
0003 #include "RKLocalFieldProvider.h"
0004 #include "CylindricalState.h"
0005
0006 #include <exception>
0007
0008 class CylindricalLorentzForceException : public std::exception {
0009 public:
0010   CylindricalLorentzForceException() throw() {}
0011   ~CylindricalLorentzForceException() throw() override {}
0012 };
0013
0014 template <typename T, int N>
0015 typename CylindricalLorentzForce<T, N>::Vector CylindricalLorentzForce<T, N>::operator()(Scalar r,
0016                                                                                          const Vector& state) const {
0017   //    std::cout << "CylindricalLorentzForce called with r " << r << " state " << state << std::endl;
0018
0019   // derivatives in case Radius is the free parameter
0020   CylindricalState cstate(r, state, 1.);
0021   LocalPoint pos(cstate.position());
0022   RKLocalFieldProvider::Vector B = theField.inTesla(pos.x(), pos.y(), pos.z());
0023   double k = 2.99792458e-3;  // conversion to [cm]
0024
0025   double dphi_dr = state[2];
0026   double dz_dr = state[3];
0027   double lambda = state[4];
0028
0029   double Q = sqrt(1.0 + dphi_dr * dphi_dr * r * r + dz_dr * dz_dr);
0030
0031   //  std::cout << "CylindricalLorentzForce: Q " << Q << " dphi_dr " << dphi_dr << " dz_dr " << dz_dr << std::endl;
0032   if (Q > 1000) {
0033     throw CylindricalLorentzForceException();
0034   }
0035
0036   double sinphi = sin(state[0]);
0037   double cosphi = cos(state[0]);
0038   double B_r = B.x() * cosphi + B.y() * sinphi;
0039   double B_phi = -B.x() * sinphi + B.y() * cosphi;
0040   double dphi_dr2 = dphi_dr * dphi_dr;
0041
0042   //  std::cout << "B_r " << B_r << " B_phi " << B_phi << " B.z() " << B.z() << std::endl;
0043
0044   double d2phi_dr2 = -2 / r * dphi_dr - r * dphi_dr * dphi_dr2 +
0045                      k * lambda * Q / r * (dz_dr * B_r + r * dphi_dr * dz_dr * B_phi - (1 + r * r * dphi_dr2) * B.z());
0046
0047   //  std::cout << "-2/r * dphi_dr - r*dphi_dr*dphi_dr2 " << -2/r * dphi_dr - r*dphi_dr*dphi_dr2
0048   //       << "  -k*lambda*Q/r * (1+r*r*dphi_dr2)*B.z() " <<  -k*lambda*Q/r * (1+r*r*dphi_dr2)*B.z() << std::endl;
0049
0050   double d2z_dr2 = -r * dphi_dr2 * dz_dr +
0051                    k * lambda * Q * (-r * dphi_dr * B_r + (1 + dz_dr * dz_dr) * B_phi - r * dphi_dr * dz_dr * B.z());
0052
0053   // the derivative of q/p is zero -- momentum conservation if no energy loss
0054   Vector res;
0055   res[0] = dphi_dr;
0056   res[1] = dz_dr;
0057   res[2] = d2phi_dr2;
0058   res[3] = d2z_dr2;
0059   res[4] = 0;
0060
0061   //  std::cout << "CylindricalLorentzForce: derivatives are " << res << std::endl;
0062
0063   return res;
0064 }