Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | List of all members
G4RepleteEofM Class Reference

#include <G4RepleteEofM.hh>

Inheritance diagram for G4RepleteEofM:
Inheritance graph
[legend]
Collaboration diagram for G4RepleteEofM:
Collaboration graph
[legend]

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)
 

Detailed Description

Definition at line 50 of file G4RepleteEofM.hh.

Constructor & Destructor Documentation

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 }
G4EquationOfMotion::G4EquationOfMotion
G4EquationOfMotion(G4Field *Field)
G4Field::IsGravityActive
G4bool IsGravityActive() const
Definition: G4Field.hh:98

Here is the call graph for this function:

G4RepleteEofM::~G4RepleteEofM ( )

Definition at line 56 of file G4RepleteEofM.cc.

57 {
58 }

Member Function Documentation

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 }
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:

G4double G4RepleteEofM::GetAnomaly ( ) const
inline

Definition at line 68 of file G4RepleteEofM.hh.

68 { return anomaly; }
void G4RepleteEofM::SetAnomaly ( G4double  a)
inline

Definition at line 67 of file G4RepleteEofM.hh.

67 { anomaly = a; }
a
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void G4RepleteEofM::SetBField ( )
inline

Definition at line 71 of file G4RepleteEofM.hh.

71 {fBfield = true;}
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 }
G4ChargeState::GetCharge
G4double GetCharge() const
Definition: G4ChargeState.hh:65
eplus
static constexpr double eplus
Definition: G4SIunits.hh:199
G4ChargeState::GetSpin
G4double GetSpin() const
Definition: G4ChargeState.hh:109
G4ChargeState::GetMagneticDipoleMoment
G4double GetMagneticDipoleMoment() const
Definition: G4ChargeState.hh:74
sqr
T sqr(const T &x)
Definition: templates.hh:145
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:

void G4RepleteEofM::SetEField ( )
inline

Definition at line 72 of file G4RepleteEofM.hh.

72 {fEfield = true;}
void G4RepleteEofM::SetgradB ( )
inline

Definition at line 73 of file G4RepleteEofM.hh.

73 {fgradB = true;}
void G4RepleteEofM::SetSpin ( )
inline

Definition at line 74 of file G4RepleteEofM.hh.

74 {fSpin = true;}
The documentation for this class was generated from the following files: