MagneticField/src/CMSTMagFieldEquation.h

0001 //
0002 // copied from G4TMagFieldEquation
0003 //
0004 // Class description:
0005 //
0006 // Templated version of equation of motion of a particle in a pure magnetic field.
0007 // Enables use of inlined code for field, equation, stepper, driver,
0008 // avoiding all virtual calls.
0009 //
0010 // Adapted from G4Mag_UsualEqRhs.hh
0011 // --------------------------------------------------------------------
0012 // Created: Josh Xie  (Google Summer of Code 2014 )
0013 // Adapted from G4Mag_UsualEqRhs
0014 //
0015 // Modified by Vincenzo innocente 7/2022 to adapt to specific CMS field
0016 //
0017 #ifndef CMSTMAGFIELDEQUATION_H
0018 #define CMSTMAGFIELDEQUATION_H
0019
0020 #include "G4Mag_UsualEqRhs.hh"
0021 #include <cassert>
0022
0023 template <class T_Field>
0024 class CMSTMagFieldEquation final : public G4Mag_UsualEqRhs {
0025 public:
0026   CMSTMagFieldEquation(T_Field* f) : G4Mag_UsualEqRhs(f) {
0027     assert(f);
0028     itsField = f;
0029   }
0030
0031   ~CMSTMagFieldEquation() override { ; }
0032
0033   inline void GetFieldValueCMS(const G4double Point[], G4double Field[]) const {
0034     itsField->GetFieldValue(Point, Field);
0035   }
0036
0037   inline void TEvaluateRhsGivenB(const G4double y[],
0038                                  G4double inv_momentum_magnitude,
0039                                  const G4double B[3],
0040                                  G4double dydx[]) const {
0041     G4double cof = FCof() * inv_momentum_magnitude;
0042
0043     dydx[0] = y[3] * inv_momentum_magnitude;  //  (d/ds)x = Vx/V
0044     dydx[1] = y[4] * inv_momentum_magnitude;  //  (d/ds)y = Vy/V
0045     dydx[2] = y[5] * inv_momentum_magnitude;  //  (d/ds)z = Vz/V
0046
0047     dydx[3] = cof * (y[4] * B[2] - y[5] * B[1]);  // Ax = a*(Vy*Bz - Vz*By)
0048     dydx[4] = cof * (y[5] * B[0] - y[3] * B[2]);  // Ay = a*(Vz*Bx - Vx*Bz)
0049     dydx[5] = cof * (y[3] * B[1] - y[4] * B[0]);  // Az = a*(Vx*By - Vy*Bx)
0050
0051     return;
0052   }
0053
0054   __attribute__((always_inline)) void TRightHandSide(const G4double y[],
0055                                                      G4double inv_momentum_magnitude,
0056                                                      G4double dydx[]) const {
0057     // CMS field  is three dimentional
0058     G4double Field[3];
0059     GetFieldValueCMS(y, Field);
0060     TEvaluateRhsGivenB(y, inv_momentum_magnitude, Field, dydx);
0061   }
0062
0063 private:
0064   // Dependent objects
0065   T_Field* itsField;
0066 };
0067
0068 #endif