Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4RepleteEofM Class Reference

#include <G4RepleteEofM.hh>

+ Inheritance diagram for G4RepleteEofM:

Public Member Functions

 G4RepleteEofM (G4Field *, G4int nvar=8)
 
 ~G4RepleteEofM ()
 
void SetChargeMomentumMass (G4ChargeState particleCharge, G4double MomentumXc, G4double mass)
 
void EvaluateRhsGivenB (const G4double y[], const G4double Field[], G4double dydx[]) const
 
void SetAnomaly (G4double a)
 
G4double GetAnomaly () const
 
void SetBField ()
 
void SetEField ()
 
void SetgradB ()
 
void SetSpin ()
 
- Public Member Functions inherited from G4EquationOfMotion
 G4EquationOfMotion (G4Field *Field)
 
virtual ~G4EquationOfMotion ()
 
virtual void EvaluateRhsGivenB (const G4double y[], const G4double B[3], G4double dydx[]) const =0
 
virtual void SetChargeMomentumMass (G4ChargeState particleCharge, G4double MomentumXc, G4double MassXc2)=0
 
void RightHandSide (const G4double y[], G4double dydx[]) const
 
void EvaluateRhsReturnB (const G4double y[], G4double dydx[], G4double Field[]) const
 
void GetFieldValue (const G4double Point[4], G4double Field[]) const
 
const G4FieldGetFieldObj () const
 
G4FieldGetFieldObj ()
 
void SetFieldObj (G4Field *pField)
 

Detailed Description

Definition at line 44 of file G4RepleteEofM.hh.

Constructor & Destructor Documentation

◆ G4RepleteEofM()

G4RepleteEofM::G4RepleteEofM ( G4Field field,
G4int  nvar = 8 
)

Definition at line 40 of file G4RepleteEofM.cc.

41 : G4EquationOfMotion( field ), fNvar(nvar)
42{
43 fGfield = field->IsGravityActive();
44}
G4bool IsGravityActive() const
Definition: G4Field.hh:101

◆ ~G4RepleteEofM()

G4RepleteEofM::~G4RepleteEofM ( )

Definition at line 46 of file G4RepleteEofM.cc.

47{
48}

Member Function Documentation

◆ EvaluateRhsGivenB()

void G4RepleteEofM::EvaluateRhsGivenB ( const G4double  y[],
const G4double  Field[],
G4double  dydx[] 
) const

Definition at line 77 of file G4RepleteEofM.cc.

80{
81
82 // Components of y:
83 // 0-2 dr/ds,
84 // 3-5 dp/ds - momentum derivatives
85 // 9-11 dSpin/ds = (1/beta) dSpin/dt - spin derivatives
86 //
87 // The BMT equation, following J.D.Jackson, Classical
88 // Electrodynamics, Second Edition,
89 // dS/dt = (e/mc) S \cross
90 // [ (g/2-1 +1/\gamma) B
91 // -(g/2-1)\gamma/(\gamma+1) (\beta \cdot B)\beta
92 // -(g/2-\gamma/(\gamma+1) \beta \cross E ]
93 // where
94 // S = \vec{s}, where S^2 = 1
95 // B = \vec{B}
96 // \beta = \vec{\beta} = \beta \vec{u} with u^2 = 1
97 // E = \vec{E}
98 //
99 // Field[0,1,2] are the magnetic field components
100 // Field[3,4,5] are the electric field components
101 // Field[6,7,8] are the gravity field components
102 // The Field[] array may trivially be extended to 18 components
103 // Field[ 9] == dB_x/dx; Field[10] == dB_y/dx; Field[11] == dB_z/dx
104 // Field[12] == dB_x/dy; Field[13] == dB_y/dy; Field[14] == dB_z/dy
105 // Field[15] == dB_x/dz; Field[16] == dB_y/dz; Field[17] == dB_z/dz
106
107 G4double momentum_mag_square = y[3]*y[3] + y[4]*y[4] + y[5]*y[5];
108 G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square );
109
110 G4double Energy = std::sqrt(momentum_mag_square + mass*mass);
111 G4double inverse_velocity = Energy*inv_momentum_magnitude/c_light;
112
113 G4double cof1 = ElectroMagCof*inv_momentum_magnitude;
114 G4double cof2 = Energy/c_light;
115 G4double cof3 = inv_momentum_magnitude*mass;
116
117 dydx[0] = y[3]*inv_momentum_magnitude; // (d/ds)x = Vx/V
118 dydx[1] = y[4]*inv_momentum_magnitude; // (d/ds)y = Vy/V
119 dydx[2] = y[5]*inv_momentum_magnitude; // (d/ds)z = Vz/V
120
121 dydx[3] = 0.;
122 dydx[4] = 0.;
123 dydx[5] = 0.;
124
125 G4double field[18] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
126
127 field[0] = Field[0];
128 field[1] = Field[1];
129 field[2] = Field[2];
130
131 // Force due to B field - Field[0,1,2]
132
133 if (fBfield)
134 {
135 if (charge != 0.)
136 {
137 dydx[3] += cof1*(y[4]*field[2] - y[5]*field[1]);
138 dydx[4] += cof1*(y[5]*field[0] - y[3]*field[2]);
139 dydx[5] += cof1*(y[3]*field[1] - y[4]*field[0]);
140 }
141 }
142
143 // add force due to E field - Field[3,4,5]
144
145 if (!fBfield)
146 {
147 field[3] = Field[0];
148 field[4] = Field[1];
149 field[5] = Field[2];
150 }
151 else
152 {
153 field[3] = Field[3];
154 field[4] = Field[4];
155 field[5] = Field[5];
156 }
157
158 if (fEfield)
159 {
160 if (charge != 0.)
161 {
162 dydx[3] += cof1*cof2*field[3];
163 dydx[4] += cof1*cof2*field[4];
164 dydx[5] += cof1*cof2*field[5];
165 }
166 }
167
168 // add force due to gravity field - Field[6,7,8]
169
170 if (!fBfield && !fEfield)
171 {
172 field[6] = Field[0];
173 field[7] = Field[1];
174 field[8] = Field[2];
175 }
176 else
177 {
178 field[6] = Field[6];
179 field[7] = Field[7];
180 field[8] = Field[8];
181 }
182
183 if (fGfield)
184 {
185 if (mass > 0.)
186 {
187 dydx[3] += field[6]*cof2*cof3/c_light;
188 dydx[4] += field[7]*cof2*cof3/c_light;
189 dydx[5] += field[8]*cof2*cof3/c_light;
190 }
191 }
192
193 // add force
194
195 if (!fBfield && !fEfield && !fGfield)
196 {
197 field[9] = Field[0];
198 field[10] = Field[1];
199 field[11] = Field[2];
200 field[12] = Field[3];
201 field[13] = Field[4];
202 field[14] = Field[5];
203 field[15] = Field[6];
204 field[16] = Field[7];
205 field[17] = Field[8];
206 }
207 else
208 {
209 field[9] = Field[9];
210 field[10] = Field[10];
211 field[11] = Field[11];
212 field[12] = Field[12];
213 field[13] = Field[13];
214 field[14] = Field[14];
215 field[15] = Field[15];
216 field[16] = Field[16];
217 field[17] = Field[17];
218 }
219
220 if (fgradB)
221 {
222 if (magMoment != 0.)
223 {
224 dydx[3] += magMoment*(y[9]*field[ 9]+y[10]*field[10]+y[11]*field[11])
225 *inv_momentum_magnitude*Energy;
226 dydx[4] += magMoment*(y[9]*field[12]+y[10]*field[13]+y[11]*field[14])
227 *inv_momentum_magnitude*Energy;
228 dydx[5] += magMoment*(y[9]*field[15]+y[10]*field[16]+y[11]*field[17])
229 *inv_momentum_magnitude*Energy;
230 }
231 }
232
233 dydx[6] = 0.; // not used
234
235 // Lab Time of flight
236 //
237 dydx[7] = inverse_velocity;
238
239 if (fNvar == 12)
240 {
241 dydx[ 8] = 0.; //not used
242
243 dydx[ 9] = 0.;
244 dydx[10] = 0.;
245 dydx[11] = 0.;
246 }
247
248 if (fSpin)
249 {
250 G4ThreeVector BField(0.,0.,0.);
251 if (fBfield)
252 {
253 G4ThreeVector F(field[0],field[1],field[2]);
254 BField = F;
255 }
256
257 G4ThreeVector EField(0.,0.,0.);
258 if (fEfield)
259 {
260 G4ThreeVector F(field[3],field[4],field[5]);
261 EField = F;
262 }
263
264 EField /= c_light;
265
266 G4ThreeVector u(y[3], y[4], y[5]);
267 u *= inv_momentum_magnitude;
268
269 G4double udb = anomaly*beta*gamma/(1.+gamma) * (BField * u);
270 G4double ucb = (anomaly+1./gamma)/beta;
271 G4double uce = anomaly + 1./(gamma+1.);
272
273 G4ThreeVector Spin(y[9],y[10],y[11]);
274
275 G4double pcharge;
276 if (charge == 0.) pcharge = 1.;
277 else pcharge = charge;
278
279 G4ThreeVector dSpin(0.,0.,0);
280 if (Spin.mag2() != 0.)
281 {
282 if (fBfield)
283 {
284 dSpin =
285 pcharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.cross(u)) );
286 }
287 if (fEfield)
288 {
289 dSpin -= pcharge*omegac*( uce*(u*(Spin*EField) - EField*(Spin*u)) );
290 // from Jackson
291 // -uce*Spin.cross(u.cross(EField)) );
292 // but this form has one less operation
293 }
294 }
295
296 dydx[ 9] = dSpin.x();
297 dydx[10] = dSpin.y();
298 dydx[11] = dSpin.z();
299 }
300
301 return;
302}
double G4double
Definition: G4Types.hh:83

◆ GetAnomaly()

G4double G4RepleteEofM::GetAnomaly ( ) const
inline

Definition at line 62 of file G4RepleteEofM.hh.

62{ return anomaly; }

◆ SetAnomaly()

void G4RepleteEofM::SetAnomaly ( G4double  a)
inline

Definition at line 61 of file G4RepleteEofM.hh.

61{ anomaly = a; }

◆ SetBField()

void G4RepleteEofM::SetBField ( )
inline

Definition at line 65 of file G4RepleteEofM.hh.

65{ fBfield = true; }

◆ SetChargeMomentumMass()

void G4RepleteEofM::SetChargeMomentumMass ( G4ChargeState  particleCharge,
G4double  MomentumXc,
G4double  mass 
)
virtual

Implements G4EquationOfMotion.

Definition at line 51 of file G4RepleteEofM.cc.

54{
55 charge = particleCharge.GetCharge();
56 mass = particleMass;
57 magMoment = particleCharge.GetMagneticDipoleMoment();
58 spin = particleCharge.GetSpin();
59
60 ElectroMagCof = eplus*charge*c_light;
61 omegac = (eplus/mass)*c_light;
62
63 G4double muB = 0.5*eplus*hbar_Planck/(mass/c_squared);
64
65 G4double g_BMT;
66 if ( spin != 0. ) g_BMT = (std::abs(magMoment)/muB)/spin;
67 else g_BMT = 2.;
68
69 anomaly = (g_BMT - 2.)/2.;
70
71 G4double E = std::sqrt(sqr(MomentumXc)+sqr(mass));
72 beta = MomentumXc/E;
73 gamma = E/mass;
74}
G4double GetCharge() const
G4double GetMagneticDipoleMoment() const
G4double GetSpin() const
T sqr(const T &x)
Definition: templates.hh:128

◆ SetEField()

void G4RepleteEofM::SetEField ( )
inline

Definition at line 66 of file G4RepleteEofM.hh.

66{ fEfield = true; }

◆ SetgradB()

void G4RepleteEofM::SetgradB ( )
inline

Definition at line 67 of file G4RepleteEofM.hh.

67{ fgradB = true; }

◆ SetSpin()

void G4RepleteEofM::SetSpin ( )
inline

Definition at line 68 of file G4RepleteEofM.hh.

68{ fSpin = true; }

The documentation for this class was generated from the following files: