include/Geant4/G4TMagFieldEquation.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4TMagFieldEquation
0027 //
0028 // Class description:
0029 //
0030 // Templated version of equation of motion of a particle in a pure magnetic field.
0031 // Enables use of inlined code for field, equation, stepper, driver,
0032 // avoiding all virtual calls.
0033 //
0034 // Adapted from G4Mag_UsualEqRhs.hh
0035 // --------------------------------------------------------------------
0036 // Created: Josh Xie  (Google Summer of Code 2014 )
0037 // Adapted from G4Mag_UsualEqRhs
0038 //
0039 // #include "G4ChargeState.hh"
0040 #include "G4Mag_UsualEqRhs.hh"
0041
0042 template
0043 <class T_Field>
0044 class G4TMagFieldEquation : public G4Mag_UsualEqRhs
0045 {
0046   public:
0047
0048     G4TMagFieldEquation(T_Field* f)
0049        : G4Mag_UsualEqRhs(f)
0050     {
0051             itsField = f;
0052     }
0053
0054     virtual ~G4TMagFieldEquation(){;}
0055
0056     inline void GetFieldValue(const G4double Point[4],
0057                               G4double Field[]) const
0058     {
0059       itsField->T_Field::GetFieldValue(Point, Field);
0060     }
0061
0062     inline void TEvaluateRhsGivenB( const G4double y[],
0063                                     const G4double B[3],
0064                                     G4double dydx[] ) const
0065     {
0066       G4double momentum_mag_square = y[3]*y[3] + y[4]*y[4] + y[5]*y[5];
0067       G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square );
0068       G4double cof = FCof()*inv_momentum_magnitude;
0069
0070       dydx[0] = y[3]*inv_momentum_magnitude;       //  (d/ds)x = Vx/V
0071       dydx[1] = y[4]*inv_momentum_magnitude;       //  (d/ds)y = Vy/V
0072       dydx[2] = y[5]*inv_momentum_magnitude;       //  (d/ds)z = Vz/V
0073
0074       dydx[3] = cof*(y[4]*B[2] - y[5]*B[1]) ;  // Ax = a*(Vy*Bz - Vz*By)
0075       dydx[4] = cof*(y[5]*B[0] - y[3]*B[2]) ;  // Ay = a*(Vz*Bx - Vx*Bz)
0076       dydx[5] = cof*(y[3]*B[1] - y[4]*B[0]) ;  // Az = a*(Vx*By - Vy*Bx)
0077
0078       return ;
0079     }
0080
0081     __attribute__((always_inline))
0082     void RightHandSide(const G4double y[], G4double dydx[] )
0083     //  const
0084     {
0085       G4double Field[G4maximum_number_of_field_components];
0086       G4double  PositionAndTime[4];
0087       PositionAndTime[0] = y[0];
0088       PositionAndTime[1] = y[1];
0089       PositionAndTime[2] = y[2];
0090       PositionAndTime[3] = y[7];
0091       GetFieldValue(PositionAndTime, Field) ;
0092       TEvaluateRhsGivenB(y, Field, dydx);
0093     }
0094
0095 private:
0096   enum { G4maximum_number_of_field_components = 24 };
0097
0098   // Dependent objects
0099   T_Field *itsField;
0100 };
0101