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