| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4Mag_SpinEqRhs implementation
27 //
28 // Created: J.Apostolakis, P.Gumplinger - 08.0
29 // -------------------------------------------
30
31 #include "G4Mag_SpinEqRhs.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4MagneticField.hh"
35 #include "G4ThreeVector.hh"
36
37 G4Mag_SpinEqRhs::G4Mag_SpinEqRhs( G4MagneticFi
38 : G4Mag_EqRhs( MagField )
39 {
40 }
41
42 G4Mag_SpinEqRhs::~G4Mag_SpinEqRhs() = default;
43
44 void
45 G4Mag_SpinEqRhs::SetChargeMomentumMass(G4Charg
46 G4doubl
47 G4doubl
48 {
49 G4Mag_EqRhs::SetChargeMomentumMass( particl
50
51 charge = particleCharge.GetCharge();
52 mass = particleMass;
53 magMoment = particleCharge.GetMagneticDipol
54 spin = particleCharge.GetSpin();
55
56 omegac = (eplus/mass)*c_light;
57
58 G4double muB = 0.5*eplus*hbar_Planck/(mass/
59
60 G4double g_BMT;
61 if ( spin != 0. )
62 {
63 g_BMT = (std::abs(magMoment)/muB)/spin;
64 }
65 else
66 {
67 g_BMT = 2.;
68 }
69
70 anomaly = (g_BMT - 2.)/2.;
71
72 G4double E = std::sqrt(sqr(MomentumXc)+sqr(
73 beta = MomentumXc/E;
74 gamma = E/mass;
75 }
76
77 void
78 G4Mag_SpinEqRhs::EvaluateRhsGivenB( const G4do
79 const G4do
80 G4do
81 {
82 G4double momentum_mag_square = sqr(y[3]) +
83 G4double inv_momentum_magnitude = 1.0 / std
84 G4double cof = FCof()*inv_momentum_magnitud
85
86 dydx[0] = y[3] * inv_momentum_magnitude;
87 dydx[1] = y[4] * inv_momentum_magnitude;
88 dydx[2] = y[5] * inv_momentum_magnitude;
89
90 if (charge == 0.)
91 {
92 dydx[3] = 0.;
93 dydx[4] = 0.;
94 dydx[5] = 0.;
95 }
96 else
97 {
98 dydx[3] = cof*(y[4]*B[2] - y[5]*B[1]) ;
99 dydx[4] = cof*(y[5]*B[0] - y[3]*B[2]) ;
100 dydx[5] = cof*(y[3]*B[1] - y[4]*B[0]) ;
101 }
102
103 G4ThreeVector u(y[3], y[4], y[5]);
104 u *= inv_momentum_magnitude;
105
106 G4ThreeVector BField(B[0],B[1],B[2]);
107
108 G4double udb = anomaly*beta*gamma/(1.+gamma
109 G4double ucb = (anomaly+1./gamma)/beta;
110
111 // Initialise the values of dydx that we do
112 dydx[6] = dydx[7] = dydx[8] = 0.0;
113
114 G4ThreeVector Spin(y[9],y[10],y[11]);
115
116 G4double pcharge;
117 if (charge == 0.)
118 {
119 pcharge = 1.;
120 }
121 else
122 {
123 pcharge = charge;
124 }
125
126 G4ThreeVector dSpin(0.,0.,0.);
127 if (Spin.mag2() != 0.)
128 {
129 dSpin = pcharge*omegac*(ucb*(Spin.cross(B
130 }
131
132 dydx[9] = dSpin.x();
133 dydx[10] = dSpin.y();
134 dydx[11] = dSpin.z();
135
136 return;
137 }
138