ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ConstRK4.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ConstRK4.hh
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4ConstRK4
27 //
28 // class description:
29 //
30 // G4ConstRK4 performs the integration of one step with error calculation
31 // in constant magnetic field. The integration method is the same as in
32 // ClassicalRK4. The field value is assumed constant for the step.
33 // This field evaluation is called only once per step.
34 // G4ConstRK4 can be used only for magnetic fields.
35 
36 // Created: J.Apostolakis, T.Nikitina - 18.09.2008
37 // -------------------------------------------------------------------
38 #ifndef G4CONSTRK4_HH
39 #define G4CONSTRK4_HH
40 
41 #include "G4MagErrorStepper.hh"
42 #include "G4EquationOfMotion.hh"
43 #include "G4Mag_EqRhs.hh"
44 
46 {
47  public: // with description
48 
49  G4ConstRK4(G4Mag_EqRhs* EquationMotion, G4int numberOfStateVariables=8);
50  ~G4ConstRK4();
51 
52  G4ConstRK4(const G4ConstRK4&) = delete;
53  G4ConstRK4& operator=(const G4ConstRK4&) = delete;
54  // Copy constructor and assignment operator not allowed
55 
56  void Stepper( const G4double y[],
57  const G4double dydx[],
58  G4double h,
59  G4double yout[],
60  G4double yerr[] );
61  void DumbStepper( const G4double yIn[],
62  const G4double dydx[],
63  G4double h,
64  G4double yOut[] ) ;
65  G4double DistChord() const;
66 
67  inline void RightHandSideConst(const G4double y[],
68  G4double dydx[] ) const;
69 
70  inline void GetConstField(const G4double y[], G4double Field[]);
71 
72  public: // without description
73 
74  G4int IntegratorOrder() const { return 4; }
75 
76  private:
77 
79  // Data stored in order to find the chord
80  G4double *dydxm, *dydxt, *yt; // scratch space - not state
82  G4Mag_EqRhs* fEq = nullptr;
83  G4double Field[3];
84 };
85 
86 // Inline methods
87 
89  G4double dydx[] ) const
90 {
91 
92  G4double momentum_mag_square = y[3]*y[3] + y[4]*y[4] + y[5]*y[5];
93  G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square );
94 
95  G4double cof = fEq->FCof()*inv_momentum_magnitude;
96 
97  dydx[0] = y[3]*inv_momentum_magnitude; // (d/ds)x = Vx/V
98  dydx[1] = y[4]*inv_momentum_magnitude; // (d/ds)y = Vy/V
99  dydx[2] = y[5]*inv_momentum_magnitude; // (d/ds)z = Vz/V
100 
101  dydx[3] = cof*(y[4]*Field[2] - y[5]*Field[1]) ; // Ax = a*(Vy*Bz - Vz*By)
102  dydx[4] = cof*(y[5]*Field[0] - y[3]*Field[2]) ; // Ay = a*(Vz*Bx - Vx*Bz)
103  dydx[5] = cof*(y[3]*Field[1] - y[4]*Field[0]) ; // Az = a*(Vx*By - Vy*Bx)
104 }
105 
107 {
108  G4double PositionAndTime[4];
109 
110  PositionAndTime[0] = y[0];
111  PositionAndTime[1] = y[1];
112  PositionAndTime[2] = y[2];
113  // Global Time
114  PositionAndTime[3] = y[7];
115  fEq -> GetFieldValue(PositionAndTime, B);
116 }
117 
118 #endif