Geant4  10.02.p03
Public Member Functions | Private Attributes | List of all members
G4RepleteEofM Class Reference

#include <G4RepleteEofM.hh>

Inheritance diagram for G4RepleteEofM:
Collaboration 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
 
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 G4Field * GetFieldObj () const
 
void SetFieldObj (G4Field *pField)
 

Private Attributes

G4int fNvar
 
G4bool fBfield
 
G4bool fEfield
 
G4bool fGfield
 
G4bool fgradB
 
G4bool fSpin
 
G4double charge
 
G4double mass
 
G4double magMoment
 
G4double spin
 
G4double ElectroMagCof
 
G4double omegac
 
G4double anomaly
 
G4double beta
 
G4double gamma
 

Detailed Description

Definition at line 50 of file G4RepleteEofM.hh.

Constructor & Destructor Documentation

◆ G4RepleteEofM()

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

Definition at line 45 of file G4RepleteEofM.cc.

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.)
52 {
53  fGfield = field->IsGravityActive();
54 }
G4RepleteEofM::charge
G4double charge
Definition: G4RepleteEofM.hh:82
G4RepleteEofM::omegac
G4double omegac
Definition: G4RepleteEofM.hh:86
G4RepleteEofM::anomaly
G4double anomaly
Definition: G4RepleteEofM.hh:86
G4Field::IsGravityActive
G4bool IsGravityActive() const
Definition: G4Field.hh:98
G4RepleteEofM::gamma
G4double gamma
Definition: G4RepleteEofM.hh:87
G4RepleteEofM::magMoment
G4double magMoment
Definition: G4RepleteEofM.hh:82
G4EquationOfMotion::G4EquationOfMotion
G4EquationOfMotion(G4Field *Field)
G4RepleteEofM::ElectroMagCof
G4double ElectroMagCof
Definition: G4RepleteEofM.hh:84
Here is the call graph for this function:

◆ ~G4RepleteEofM()

G4RepleteEofM::~G4RepleteEofM ( )

Definition at line 56 of file G4RepleteEofM.cc.

57 {
58 }

Member Function Documentation

◆ EvaluateRhsGivenB()

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

Definition at line 87 of file G4RepleteEofM.cc.

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

◆ GetAnomaly()

G4double G4RepleteEofM::GetAnomaly ( ) const
inline

Definition at line 68 of file G4RepleteEofM.hh.

68 { return anomaly; }
G4RepleteEofM::anomaly
G4double anomaly
Definition: G4RepleteEofM.hh:86

◆ SetAnomaly()

void G4RepleteEofM::SetAnomaly ( G4double  a)
inline

Definition at line 67 of file G4RepleteEofM.hh.

67 { anomaly = a; }
G4RepleteEofM::anomaly
G4double anomaly
Definition: G4RepleteEofM.hh:86

◆ SetBField()

void G4RepleteEofM::SetBField ( )
inline

Definition at line 71 of file G4RepleteEofM.hh.

71 {fBfield = true;}

◆ SetChargeMomentumMass()

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

Implements G4EquationOfMotion.

Definition at line 61 of file G4RepleteEofM.cc.

64 {
65  charge = particleCharge.GetCharge();
66  mass = particleMass;
67  magMoment = particleCharge.GetMagneticDipoleMoment();
68  spin = particleCharge.GetSpin();
69 
70  ElectroMagCof = eplus*charge*c_light;
71  omegac = (eplus/mass)*c_light;
72 
73  G4double muB = 0.5*eplus*hbar_Planck/(mass/c_squared);
74 
75  G4double g_BMT;
76  if ( spin != 0. ) g_BMT = (std::abs(magMoment)/muB)/spin;
77  else g_BMT = 2.;
78 
79  anomaly = (g_BMT - 2.)/2.;
80 
81  G4double E = std::sqrt(sqr(MomentumXc)+sqr(mass));
82  beta = MomentumXc/E;
83  gamma = E/mass;
84 }
python.hepunit.c_squared
float c_squared
Definition: hepunit.py:258
G4RepleteEofM::charge
G4double charge
Definition: G4RepleteEofM.hh:82
G4RepleteEofM::omegac
G4double omegac
Definition: G4RepleteEofM.hh:86
G4RepleteEofM::anomaly
G4double anomaly
Definition: G4RepleteEofM.hh:86
G4RepleteEofM::gamma
G4double gamma
Definition: G4RepleteEofM.hh:87
G4ChargeState::GetSpin
G4double GetSpin() const
Definition: G4ChargeState.hh:109
G4ChargeState::GetCharge
G4double GetCharge() const
Definition: G4ChargeState.hh:65
G4RepleteEofM::magMoment
G4double magMoment
Definition: G4RepleteEofM.hh:82
G4ChargeState::GetMagneticDipoleMoment
G4double GetMagneticDipoleMoment() const
Definition: G4ChargeState.hh:74
G4RepleteEofM::ElectroMagCof
G4double ElectroMagCof
Definition: G4RepleteEofM.hh:84
python.hepunit.hbar_Planck
float hbar_Planck
Definition: hepunit.py:264
sqr
T sqr(const T &x)
Definition: templates.hh:145
G4double
double G4double
Definition: G4Types.hh:76
eplus
static const double eplus
Definition: G4SIunits.hh:196
python.hepunit.c_light
float c_light
Definition: hepunit.py:257
Here is the call graph for this function:

◆ SetEField()

void G4RepleteEofM::SetEField ( )
inline

Definition at line 72 of file G4RepleteEofM.hh.

72 {fEfield = true;}

◆ SetgradB()

void G4RepleteEofM::SetgradB ( )
inline

Definition at line 73 of file G4RepleteEofM.hh.

73 {fgradB = true;}

◆ SetSpin()

void G4RepleteEofM::SetSpin ( )
inline

Definition at line 74 of file G4RepleteEofM.hh.

74 {fSpin = true;}

Member Data Documentation

◆ anomaly

G4double G4RepleteEofM::anomaly
private

Definition at line 86 of file G4RepleteEofM.hh.

◆ beta

G4double G4RepleteEofM::beta
private

Definition at line 87 of file G4RepleteEofM.hh.

◆ charge

G4double G4RepleteEofM::charge
private

Definition at line 82 of file G4RepleteEofM.hh.

◆ ElectroMagCof

G4double G4RepleteEofM::ElectroMagCof
private

Definition at line 84 of file G4RepleteEofM.hh.

◆ fBfield

G4bool G4RepleteEofM::fBfield
private

Definition at line 80 of file G4RepleteEofM.hh.

◆ fEfield

G4bool G4RepleteEofM::fEfield
private

Definition at line 80 of file G4RepleteEofM.hh.

◆ fGfield

G4bool G4RepleteEofM::fGfield
private

Definition at line 80 of file G4RepleteEofM.hh.

◆ fgradB

G4bool G4RepleteEofM::fgradB
private

Definition at line 80 of file G4RepleteEofM.hh.

◆ fNvar

G4int G4RepleteEofM::fNvar
private

Definition at line 78 of file G4RepleteEofM.hh.

◆ fSpin

G4bool G4RepleteEofM::fSpin
private

Definition at line 80 of file G4RepleteEofM.hh.

◆ gamma

G4double G4RepleteEofM::gamma
private

Definition at line 87 of file G4RepleteEofM.hh.

◆ magMoment

G4double G4RepleteEofM::magMoment
private

Definition at line 82 of file G4RepleteEofM.hh.

◆ mass

G4double G4RepleteEofM::mass
private

Definition at line 82 of file G4RepleteEofM.hh.

◆ omegac

G4double G4RepleteEofM::omegac
private

Definition at line 86 of file G4RepleteEofM.hh.

◆ spin

G4double G4RepleteEofM::spin
private

Definition at line 82 of file G4RepleteEofM.hh.

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