| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4RepleteEofM implementation <<
27 // 26 //
28 // Created: P.Gumplinger, 08.04.2013 << 27 // $Id: G4RepleteEofM.cc $
>> 28 //
>> 29 //
>> 30 // This is the standard right-hand side for equation of motion.
>> 31 //
>> 32 // 08.04.2013 Peter Gumplinger
>> 33 //
29 // ------------------------------------------- 34 // -------------------------------------------------------------------
30 35
31 #include "G4RepleteEofM.hh" 36 #include "G4RepleteEofM.hh"
32 #include "G4Field.hh" 37 #include "G4Field.hh"
33 #include "G4ThreeVector.hh" 38 #include "G4ThreeVector.hh"
34 #include "globals.hh" 39 #include "globals.hh"
35 40
36 #include "G4PhysicalConstants.hh" 41 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh"
38 43
39 44
40 G4RepleteEofM::G4RepleteEofM( G4Field* field, 45 G4RepleteEofM::G4RepleteEofM( G4Field* field, G4int nvar )
41 : G4EquationOfMotion( field ), fNvar << 46 : G4EquationOfMotion( field ), fNvar(nvar),
>> 47 fBfield(false), fEfield(false), fGfield(false),
>> 48 fgradB(false), fSpin(false),
>> 49 charge(0.), mass(0.), magMoment(0.), spin(0.),
>> 50 ElectroMagCof(0.), omegac(0.), anomaly(0.),
>> 51 beta(0.), gamma(0.)
42 { 52 {
43 fGfield = field->IsGravityActive(); 53 fGfield = field->IsGravityActive();
44 } 54 }
45 55
46 G4RepleteEofM::~G4RepleteEofM() = default; << 56 G4RepleteEofM::~G4RepleteEofM()
>> 57 {
>> 58 }
47 59
48 void 60 void
49 G4RepleteEofM::SetChargeMomentumMass(G4ChargeS 61 G4RepleteEofM::SetChargeMomentumMass(G4ChargeState particleCharge, // e+ units
50 G4double Momentu 62 G4double MomentumXc,
51 G4double particl 63 G4double particleMass)
52 { 64 {
53 charge = particleCharge.GetCharge(); 65 charge = particleCharge.GetCharge();
54 mass = particleMass; 66 mass = particleMass;
55 magMoment = particleCharge.GetMagneticDipol 67 magMoment = particleCharge.GetMagneticDipoleMoment();
56 spin = particleCharge.GetSpin(); 68 spin = particleCharge.GetSpin();
57 69
58 ElectroMagCof = eplus*charge*c_light; 70 ElectroMagCof = eplus*charge*c_light;
59 omegac = (eplus/mass)*c_light; 71 omegac = (eplus/mass)*c_light;
60 72
61 G4double muB = 0.5*eplus*hbar_Planck/(mass/ 73 G4double muB = 0.5*eplus*hbar_Planck/(mass/c_squared);
62 74
63 G4double g_BMT; 75 G4double g_BMT;
64 if ( spin != 0. ) << 76 if ( spin != 0. ) g_BMT = (std::abs(magMoment)/muB)/spin;
65 { << 77 else g_BMT = 2.;
66 g_BMT = (std::abs(magMoment)/muB)/spin; <<
67 } <<
68 else <<
69 { <<
70 g_BMT = 2.; <<
71 } <<
72 78
73 anomaly = (g_BMT - 2.)/2.; 79 anomaly = (g_BMT - 2.)/2.;
74 80
75 G4double E = std::sqrt(sqr(MomentumXc)+sqr( 81 G4double E = std::sqrt(sqr(MomentumXc)+sqr(mass));
76 beta = MomentumXc/E; 82 beta = MomentumXc/E;
77 gamma = E/mass; 83 gamma = E/mass;
78 } 84 }
79 85
80 void 86 void
81 G4RepleteEofM::EvaluateRhsGivenB( const G4doub << 87 G4RepleteEofM::EvaluateRhsGivenB(const G4double y[],
82 const G4doub << 88 const G4double Field[],
83 G4doub << 89 G4double dydx[] ) const
84 { 90 {
85 91
86 // Components of y: 92 // Components of y:
87 // 0-2 dr/ds, 93 // 0-2 dr/ds,
88 // 3-5 dp/ds - momentum derivatives 94 // 3-5 dp/ds - momentum derivatives
89 // 9-11 dSpin/ds = (1/beta) dSpin/dt - s 95 // 9-11 dSpin/ds = (1/beta) dSpin/dt - spin derivatives
90 // 96 //
91 // The BMT equation, following J.D.Jackson, 97 // The BMT equation, following J.D.Jackson, Classical
92 // Electrodynamics, Second Edition, 98 // Electrodynamics, Second Edition,
93 // dS/dt = (e/mc) S \cross 99 // dS/dt = (e/mc) S \cross
94 // [ (g/2-1 +1/\gamma) B 100 // [ (g/2-1 +1/\gamma) B
95 // -(g/2-1)\gamma/(\gamma+1) 101 // -(g/2-1)\gamma/(\gamma+1) (\beta \cdot B)\beta
96 // -(g/2-\gamma/(\gamma+1) \b 102 // -(g/2-\gamma/(\gamma+1) \beta \cross E ]
97 // where 103 // where
98 // S = \vec{s}, where S^2 = 1 104 // S = \vec{s}, where S^2 = 1
99 // B = \vec{B} 105 // B = \vec{B}
100 // \beta = \vec{\beta} = \beta \vec{u} with 106 // \beta = \vec{\beta} = \beta \vec{u} with u^2 = 1
101 // E = \vec{E} 107 // E = \vec{E}
102 // 108 //
103 // Field[0,1,2] are the magnetic field comp 109 // Field[0,1,2] are the magnetic field components
104 // Field[3,4,5] are the electric field comp 110 // Field[3,4,5] are the electric field components
105 // Field[6,7,8] are the gravity field comp 111 // Field[6,7,8] are the gravity field components
106 // The Field[] array may trivially be exten 112 // The Field[] array may trivially be extended to 18 components
107 // Field[ 9] == dB_x/dx; Field[10] == dB_y/ 113 // Field[ 9] == dB_x/dx; Field[10] == dB_y/dx; Field[11] == dB_z/dx
108 // Field[12] == dB_x/dy; Field[13] == dB_y/ 114 // Field[12] == dB_x/dy; Field[13] == dB_y/dy; Field[14] == dB_z/dy
109 // Field[15] == dB_x/dz; Field[16] == dB_y/ 115 // Field[15] == dB_x/dz; Field[16] == dB_y/dz; Field[17] == dB_z/dz
110 116
111 G4double momentum_mag_square = y[3]*y[3] + 117 G4double momentum_mag_square = y[3]*y[3] + y[4]*y[4] + y[5]*y[5];
112 G4double inv_momentum_magnitude = 1.0 / std 118 G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square );
113 119
114 G4double Energy = std::sqrt(momentum_mag_sq 120 G4double Energy = std::sqrt(momentum_mag_square + mass*mass);
115 G4double inverse_velocity = Energy*inv_mome 121 G4double inverse_velocity = Energy*inv_momentum_magnitude/c_light;
116 122
117 G4double cof1 = ElectroMagCof*inv_momentum_ 123 G4double cof1 = ElectroMagCof*inv_momentum_magnitude;
118 G4double cof2 = Energy/c_light; 124 G4double cof2 = Energy/c_light;
119 G4double cof3 = inv_momentum_magnitude*mass 125 G4double cof3 = inv_momentum_magnitude*mass;
120 126
121 dydx[0] = y[3]*inv_momentum_magnitude; 127 dydx[0] = y[3]*inv_momentum_magnitude; // (d/ds)x = Vx/V
122 dydx[1] = y[4]*inv_momentum_magnitude; 128 dydx[1] = y[4]*inv_momentum_magnitude; // (d/ds)y = Vy/V
123 dydx[2] = y[5]*inv_momentum_magnitude; 129 dydx[2] = y[5]*inv_momentum_magnitude; // (d/ds)z = Vz/V
124 130
125 dydx[3] = 0.; 131 dydx[3] = 0.;
126 dydx[4] = 0.; 132 dydx[4] = 0.;
127 dydx[5] = 0.; 133 dydx[5] = 0.;
128 134
129 G4double field[18] = {0.,0.,0.,0.,0.,0.,0., 135 G4double field[18] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
130 136
131 field[0] = Field[0]; 137 field[0] = Field[0];
132 field[1] = Field[1]; 138 field[1] = Field[1];
133 field[2] = Field[2]; 139 field[2] = Field[2];
134 140
135 // Force due to B field - Field[0,1,2] 141 // Force due to B field - Field[0,1,2]
136 142
137 if (fBfield) << 143 if (fBfield) {
138 { << 144 if (charge != 0.) {
139 if (charge != 0.) <<
140 { <<
141 dydx[3] += cof1*(y[4]*field[2] - y[5] 145 dydx[3] += cof1*(y[4]*field[2] - y[5]*field[1]);
142 dydx[4] += cof1*(y[5]*field[0] - y[3] 146 dydx[4] += cof1*(y[5]*field[0] - y[3]*field[2]);
143 dydx[5] += cof1*(y[3]*field[1] - y[4] 147 dydx[5] += cof1*(y[3]*field[1] - y[4]*field[0]);
144 } 148 }
145 } 149 }
146 150
147 // add force due to E field - Field[3,4,5] 151 // add force due to E field - Field[3,4,5]
148 152
149 if (!fBfield) << 153 if (!fBfield) {
150 { <<
151 field[3] = Field[0]; 154 field[3] = Field[0];
152 field[4] = Field[1]; 155 field[4] = Field[1];
153 field[5] = Field[2]; 156 field[5] = Field[2];
154 } << 157 } else {
155 else <<
156 { <<
157 field[3] = Field[3]; 158 field[3] = Field[3];
158 field[4] = Field[4]; 159 field[4] = Field[4];
159 field[5] = Field[5]; 160 field[5] = Field[5];
160 } 161 }
161 162
162 if (fEfield) << 163 if (fEfield) {
163 { << 164 if (charge != 0.) {
164 if (charge != 0.) <<
165 { <<
166 dydx[3] += cof1*cof2*field[3]; 165 dydx[3] += cof1*cof2*field[3];
167 dydx[4] += cof1*cof2*field[4]; 166 dydx[4] += cof1*cof2*field[4];
168 dydx[5] += cof1*cof2*field[5]; 167 dydx[5] += cof1*cof2*field[5];
169 } 168 }
170 } 169 }
171 170
172 // add force due to gravity field - Field[6 171 // add force due to gravity field - Field[6,7,8]
173 172
174 if (!fBfield && !fEfield) << 173 if (!fBfield && !fEfield) {
175 { <<
176 field[6] = Field[0]; 174 field[6] = Field[0];
177 field[7] = Field[1]; 175 field[7] = Field[1];
178 field[8] = Field[2]; 176 field[8] = Field[2];
179 } << 177 } else {
180 else <<
181 { <<
182 field[6] = Field[6]; 178 field[6] = Field[6];
183 field[7] = Field[7]; 179 field[7] = Field[7];
184 field[8] = Field[8]; 180 field[8] = Field[8];
185 } 181 }
186 182
187 if (fGfield) << 183 if (fGfield) {
188 { << 184 if (mass > 0.) {
189 if (mass > 0.) <<
190 { <<
191 dydx[3] += field[6]*cof2*cof3/c_light 185 dydx[3] += field[6]*cof2*cof3/c_light;
192 dydx[4] += field[7]*cof2*cof3/c_light 186 dydx[4] += field[7]*cof2*cof3/c_light;
193 dydx[5] += field[8]*cof2*cof3/c_light 187 dydx[5] += field[8]*cof2*cof3/c_light;
194 } 188 }
195 } 189 }
196 190
197 // add force << 191 // add force due to ∇(µ⋅B) == (µ⋅∇)B when (∇xB) = 0
198 192
199 if (!fBfield && !fEfield && !fGfield) << 193 if (!fBfield && !fEfield && !fGfield) {
200 { <<
201 field[9] = Field[0]; 194 field[9] = Field[0];
202 field[10] = Field[1]; 195 field[10] = Field[1];
203 field[11] = Field[2]; 196 field[11] = Field[2];
204 field[12] = Field[3]; 197 field[12] = Field[3];
205 field[13] = Field[4]; 198 field[13] = Field[4];
206 field[14] = Field[5]; 199 field[14] = Field[5];
207 field[15] = Field[6]; 200 field[15] = Field[6];
208 field[16] = Field[7]; 201 field[16] = Field[7];
209 field[17] = Field[8]; 202 field[17] = Field[8];
210 } << 203 } else {
211 else <<
212 { <<
213 field[9] = Field[9]; 204 field[9] = Field[9];
214 field[10] = Field[10]; 205 field[10] = Field[10];
215 field[11] = Field[11]; 206 field[11] = Field[11];
216 field[12] = Field[12]; 207 field[12] = Field[12];
217 field[13] = Field[13]; 208 field[13] = Field[13];
218 field[14] = Field[14]; 209 field[14] = Field[14];
219 field[15] = Field[15]; 210 field[15] = Field[15];
220 field[16] = Field[16]; 211 field[16] = Field[16];
221 field[17] = Field[17]; 212 field[17] = Field[17];
222 } 213 }
223 214
224 if (fgradB) << 215 if (fgradB) {
225 { << 216 if (magMoment != 0.) {
226 if (magMoment != 0.) << 217
227 { << 218 // field[ 9] == dB_x/dx; field[10] == dB_y/dx; field[11] == dB_z/dx
>> 219 // field[12] == dB_x/dy; field[13] == dB_y/dy; field[14] == dB_z/dy
>> 220 // field[15] == dB_x/dz; field[16] == dB_y/dz; field[17] == dB_z/dz
>> 221
>> 222 // G4cout << "y[9]: " << y[9] << " y[10]: " << y[10] << " y[11]: " << y[11] << G4endl;
>> 223 // G4cout << "field[9]: " << field[9] << " field[10]: " << field[10] << " field[11]: " << field[11] << G4endl;
>> 224 // G4cout << "field[12]: " << field[12] << " field[13]: " << field[13] << " field[14]: " << field[14] << G4endl;
>> 225 // G4cout << "field[15]: " << field[15] << " field[16]: " << field[16] << " field[17]: " << field[17] << G4endl;
>> 226 // G4cout << "inv_momentum_magnitdue: " << inv_momentum_magnitude << " Energy: " << Energy << G4endl;
>> 227
228 dydx[3] += magMoment*(y[9]*field[ 9]+ 228 dydx[3] += magMoment*(y[9]*field[ 9]+y[10]*field[10]+y[11]*field[11])
229 229 *inv_momentum_magnitude*Energy;
230 dydx[4] += magMoment*(y[9]*field[12]+ 230 dydx[4] += magMoment*(y[9]*field[12]+y[10]*field[13]+y[11]*field[14])
231 231 *inv_momentum_magnitude*Energy;
232 dydx[5] += magMoment*(y[9]*field[15]+ 232 dydx[5] += magMoment*(y[9]*field[15]+y[10]*field[16]+y[11]*field[17])
233 233 *inv_momentum_magnitude*Energy;
>> 234
>> 235 // G4cout << "dydx[3,4,5] " << dydx[3] << " " << dydx[4] << " " << dydx[5] << G4endl;
234 } 236 }
235 } 237 }
236 238
237 dydx[6] = 0.; // not used << 239 dydx[6] = 0.; //not used
238 240
239 // Lab Time of flight 241 // Lab Time of flight
240 // <<
241 dydx[7] = inverse_velocity; 242 dydx[7] = inverse_velocity;
242 243
243 if (fNvar == 12) << 244 if (fNvar == 12) {
244 { <<
245 dydx[ 8] = 0.; //not used 245 dydx[ 8] = 0.; //not used
246 246
247 dydx[ 9] = 0.; 247 dydx[ 9] = 0.;
248 dydx[10] = 0.; 248 dydx[10] = 0.;
249 dydx[11] = 0.; 249 dydx[11] = 0.;
250 } 250 }
251 251
252 if (fSpin) << 252 if (fSpin) {
253 { << 253 // G4cout << "y[9,10,11] " << y[9] << " " << y[10] << " " << y[11] << G4endl;
254 G4ThreeVector BField(0.,0.,0.); 254 G4ThreeVector BField(0.,0.,0.);
255 if (fBfield) << 255 if (fBfield) {
256 { <<
257 G4ThreeVector F(field[0],field[1],fie 256 G4ThreeVector F(field[0],field[1],field[2]);
258 BField = F; 257 BField = F;
259 } 258 }
260 259
261 G4ThreeVector EField(0.,0.,0.); 260 G4ThreeVector EField(0.,0.,0.);
262 if (fEfield) << 261 if (fEfield) {
263 { <<
264 G4ThreeVector F(field[3],field[4],fie 262 G4ThreeVector F(field[3],field[4],field[5]);
265 EField = F; 263 EField = F;
266 } 264 }
267 265
268 EField /= c_light; 266 EField /= c_light;
269 267
270 G4ThreeVector u(y[3], y[4], y[5]); 268 G4ThreeVector u(y[3], y[4], y[5]);
271 u *= inv_momentum_magnitude; 269 u *= inv_momentum_magnitude;
272 270
273 G4double udb = anomaly*beta*gamma/(1.+ga 271 G4double udb = anomaly*beta*gamma/(1.+gamma) * (BField * u);
274 G4double ucb = (anomaly+1./gamma)/beta; 272 G4double ucb = (anomaly+1./gamma)/beta;
275 G4double uce = anomaly + 1./(gamma+1.); 273 G4double uce = anomaly + 1./(gamma+1.);
276 274
277 G4ThreeVector Spin(y[9],y[10],y[11]); 275 G4ThreeVector Spin(y[9],y[10],y[11]);
278 276
279 G4double pcharge; 277 G4double pcharge;
280 if (charge == 0.) << 278 if (charge == 0.) pcharge = 1.;
281 { << 279 else pcharge = charge;
282 pcharge = 1.; <<
283 } <<
284 else <<
285 { <<
286 pcharge = charge; <<
287 } <<
288 280
289 G4ThreeVector dSpin(0.,0.,0); 281 G4ThreeVector dSpin(0.,0.,0);
290 if (Spin.mag2() != 0.) << 282 if (Spin.mag2() != 0.) {
291 { << 283 if (fBfield) {
292 if (fBfield) <<
293 { <<
294 dSpin = 284 dSpin =
295 pcharge*omegac*( ucb*(Spin.cross( 285 pcharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.cross(u)) );
296 } 286 }
297 if (fEfield) << 287 if (fEfield) {
298 { << 288 dSpin -=
299 dSpin -= pcharge*omegac*( uce*(u*( << 289 // from Jackson
300 // from Jackson << 290 // -uce*Spin.cross(u.cross(EField)) );
301 // -uce*Spin.cross(u.cross(EFiel << 291 // but this form has one less operation
302 // but this form has one less op << 292 pcharge*omegac*( uce*(u*(Spin*EField) - EField*(Spin*u)) );
303 } 293 }
304 } 294 }
305 295
306 dydx[ 9] = dSpin.x(); 296 dydx[ 9] = dSpin.x();
307 dydx[10] = dSpin.y(); 297 dydx[10] = dSpin.y();
308 dydx[11] = dSpin.z(); 298 dydx[11] = dSpin.z();
>> 299
>> 300 // G4cout << "dydx[9,10,11] " << dydx[9] << " " << dydx[10] << " " << dydx[11] << G4endl;
309 } 301 }
310 302
311 return; << 303 return ;
312 } 304 }
313 305