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