| 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 // G4RepleteEofM implementation
27 //
28 // Created: P.Gumplinger, 08.04.2013
29 // -------------------------------------------
30
31 #include "G4RepleteEofM.hh"
32 #include "G4Field.hh"
33 #include "G4ThreeVector.hh"
34 #include "globals.hh"
35
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38
39
40 G4RepleteEofM::G4RepleteEofM( G4Field* field,
41 : G4EquationOfMotion( field ), fNvar
42 {
43 fGfield = field->IsGravityActive();
44 }
45
46 G4RepleteEofM::~G4RepleteEofM() = default;
47
48 void
49 G4RepleteEofM::SetChargeMomentumMass(G4ChargeS
50 G4double Momentu
51 G4double particl
52 {
53 charge = particleCharge.GetCharge();
54 mass = particleMass;
55 magMoment = particleCharge.GetMagneticDipol
56 spin = particleCharge.GetSpin();
57
58 ElectroMagCof = eplus*charge*c_light;
59 omegac = (eplus/mass)*c_light;
60
61 G4double muB = 0.5*eplus*hbar_Planck/(mass/
62
63 G4double g_BMT;
64 if ( spin != 0. )
65 {
66 g_BMT = (std::abs(magMoment)/muB)/spin;
67 }
68 else
69 {
70 g_BMT = 2.;
71 }
72
73 anomaly = (g_BMT - 2.)/2.;
74
75 G4double E = std::sqrt(sqr(MomentumXc)+sqr(
76 beta = MomentumXc/E;
77 gamma = E/mass;
78 }
79
80 void
81 G4RepleteEofM::EvaluateRhsGivenB( const G4doub
82 const G4doub
83 G4doub
84 {
85
86 // Components of y:
87 // 0-2 dr/ds,
88 // 3-5 dp/ds - momentum derivatives
89 // 9-11 dSpin/ds = (1/beta) dSpin/dt - s
90 //
91 // The BMT equation, following J.D.Jackson,
92 // Electrodynamics, Second Edition,
93 // dS/dt = (e/mc) S \cross
94 // [ (g/2-1 +1/\gamma) B
95 // -(g/2-1)\gamma/(\gamma+1)
96 // -(g/2-\gamma/(\gamma+1) \b
97 // where
98 // S = \vec{s}, where S^2 = 1
99 // B = \vec{B}
100 // \beta = \vec{\beta} = \beta \vec{u} with
101 // E = \vec{E}
102 //
103 // Field[0,1,2] are the magnetic field comp
104 // Field[3,4,5] are the electric field comp
105 // Field[6,7,8] are the gravity field comp
106 // The Field[] array may trivially be exten
107 // Field[ 9] == dB_x/dx; Field[10] == dB_y/
108 // Field[12] == dB_x/dy; Field[13] == dB_y/
109 // Field[15] == dB_x/dz; Field[16] == dB_y/
110
111 G4double momentum_mag_square = y[3]*y[3] +
112 G4double inv_momentum_magnitude = 1.0 / std
113
114 G4double Energy = std::sqrt(momentum_mag_sq
115 G4double inverse_velocity = Energy*inv_mome
116
117 G4double cof1 = ElectroMagCof*inv_momentum_
118 G4double cof2 = Energy/c_light;
119 G4double cof3 = inv_momentum_magnitude*mass
120
121 dydx[0] = y[3]*inv_momentum_magnitude;
122 dydx[1] = y[4]*inv_momentum_magnitude;
123 dydx[2] = y[5]*inv_momentum_magnitude;
124
125 dydx[3] = 0.;
126 dydx[4] = 0.;
127 dydx[5] = 0.;
128
129 G4double field[18] = {0.,0.,0.,0.,0.,0.,0.,
130
131 field[0] = Field[0];
132 field[1] = Field[1];
133 field[2] = Field[2];
134
135 // Force due to B field - Field[0,1,2]
136
137 if (fBfield)
138 {
139 if (charge != 0.)
140 {
141 dydx[3] += cof1*(y[4]*field[2] - y[5]
142 dydx[4] += cof1*(y[5]*field[0] - y[3]
143 dydx[5] += cof1*(y[3]*field[1] - y[4]
144 }
145 }
146
147 // add force due to E field - Field[3,4,5]
148
149 if (!fBfield)
150 {
151 field[3] = Field[0];
152 field[4] = Field[1];
153 field[5] = Field[2];
154 }
155 else
156 {
157 field[3] = Field[3];
158 field[4] = Field[4];
159 field[5] = Field[5];
160 }
161
162 if (fEfield)
163 {
164 if (charge != 0.)
165 {
166 dydx[3] += cof1*cof2*field[3];
167 dydx[4] += cof1*cof2*field[4];
168 dydx[5] += cof1*cof2*field[5];
169 }
170 }
171
172 // add force due to gravity field - Field[6
173
174 if (!fBfield && !fEfield)
175 {
176 field[6] = Field[0];
177 field[7] = Field[1];
178 field[8] = Field[2];
179 }
180 else
181 {
182 field[6] = Field[6];
183 field[7] = Field[7];
184 field[8] = Field[8];
185 }
186
187 if (fGfield)
188 {
189 if (mass > 0.)
190 {
191 dydx[3] += field[6]*cof2*cof3/c_light
192 dydx[4] += field[7]*cof2*cof3/c_light
193 dydx[5] += field[8]*cof2*cof3/c_light
194 }
195 }
196
197 // add force
198
199 if (!fBfield && !fEfield && !fGfield)
200 {
201 field[9] = Field[0];
202 field[10] = Field[1];
203 field[11] = Field[2];
204 field[12] = Field[3];
205 field[13] = Field[4];
206 field[14] = Field[5];
207 field[15] = Field[6];
208 field[16] = Field[7];
209 field[17] = Field[8];
210 }
211 else
212 {
213 field[9] = Field[9];
214 field[10] = Field[10];
215 field[11] = Field[11];
216 field[12] = Field[12];
217 field[13] = Field[13];
218 field[14] = Field[14];
219 field[15] = Field[15];
220 field[16] = Field[16];
221 field[17] = Field[17];
222 }
223
224 if (fgradB)
225 {
226 if (magMoment != 0.)
227 {
228 dydx[3] += magMoment*(y[9]*field[ 9]+
229
230 dydx[4] += magMoment*(y[9]*field[12]+
231
232 dydx[5] += magMoment*(y[9]*field[15]+
233
234 }
235 }
236
237 dydx[6] = 0.; // not used
238
239 // Lab Time of flight
240 //
241 dydx[7] = inverse_velocity;
242
243 if (fNvar == 12)
244 {
245 dydx[ 8] = 0.; //not used
246
247 dydx[ 9] = 0.;
248 dydx[10] = 0.;
249 dydx[11] = 0.;
250 }
251
252 if (fSpin)
253 {
254 G4ThreeVector BField(0.,0.,0.);
255 if (fBfield)
256 {
257 G4ThreeVector F(field[0],field[1],fie
258 BField = F;
259 }
260
261 G4ThreeVector EField(0.,0.,0.);
262 if (fEfield)
263 {
264 G4ThreeVector F(field[3],field[4],fie
265 EField = F;
266 }
267
268 EField /= c_light;
269
270 G4ThreeVector u(y[3], y[4], y[5]);
271 u *= inv_momentum_magnitude;
272
273 G4double udb = anomaly*beta*gamma/(1.+ga
274 G4double ucb = (anomaly+1./gamma)/beta;
275 G4double uce = anomaly + 1./(gamma+1.);
276
277 G4ThreeVector Spin(y[9],y[10],y[11]);
278
279 G4double pcharge;
280 if (charge == 0.)
281 {
282 pcharge = 1.;
283 }
284 else
285 {
286 pcharge = charge;
287 }
288
289 G4ThreeVector dSpin(0.,0.,0);
290 if (Spin.mag2() != 0.)
291 {
292 if (fBfield)
293 {
294 dSpin =
295 pcharge*omegac*( ucb*(Spin.cross(
296 }
297 if (fEfield)
298 {
299 dSpin -= pcharge*omegac*( uce*(u*(
300 // from Jackson
301 // -uce*Spin.cross(u.cross(EFiel
302 // but this form has one less op
303 }
304 }
305
306 dydx[ 9] = dSpin.x();
307 dydx[10] = dSpin.y();
308 dydx[11] = dSpin.z();
309 }
310
311 return;
312 }
313