ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Mag_SpinEqRhs.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Mag_SpinEqRhs.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 // G4Mag_SpinEqRhs implementation
27 //
28 // Created: J.Apostolakis, P.Gumplinger - 08.02.1999
29 // --------------------------------------------------------------------
30 
31 #include "G4Mag_SpinEqRhs.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4MagneticField.hh"
35 #include "G4ThreeVector.hh"
36 
38  : G4Mag_EqRhs( MagField )
39 {
40 }
41 
43 {
44 }
45 
46 void
48  G4double MomentumXc,
49  G4double particleMass)
50 {
51  G4Mag_EqRhs::SetChargeMomentumMass( particleCharge, MomentumXc, mass);
52 
53  charge = particleCharge.GetCharge();
54  mass = particleMass;
55  magMoment = particleCharge.GetMagneticDipoleMoment();
56  spin = particleCharge.GetSpin();
57 
59 
61 
62  G4double g_BMT;
63  if ( spin != 0. ) g_BMT = (std::abs(magMoment)/muB)/spin;
64  else g_BMT = 2.;
65 
66  anomaly = (g_BMT - 2.)/2.;
67 
68  G4double E = std::sqrt(sqr(MomentumXc)+sqr(mass));
69  beta = MomentumXc/E;
70  gamma = E/mass;
71 }
72 
73 void
75  const G4double B[3],
76  G4double dydx[] ) const
77 {
78  G4double momentum_mag_square = sqr(y[3]) + sqr(y[4]) + sqr(y[5]);
79  G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square );
80  G4double cof = FCof()*inv_momentum_magnitude;
81 
82  dydx[0] = y[3] * inv_momentum_magnitude; // (d/ds)x = Vx/V
83  dydx[1] = y[4] * inv_momentum_magnitude; // (d/ds)y = Vy/V
84  dydx[2] = y[5] * inv_momentum_magnitude; // (d/ds)z = Vz/V
85 
86  if (charge == 0.)
87  {
88  dydx[3] = 0.;
89  dydx[4] = 0.;
90  dydx[5] = 0.;
91  }
92  else
93  {
94  dydx[3] = cof*(y[4]*B[2] - y[5]*B[1]) ; // Ax = a*(Vy*Bz - Vz*By)
95  dydx[4] = cof*(y[5]*B[0] - y[3]*B[2]) ; // Ay = a*(Vz*Bx - Vx*Bz)
96  dydx[5] = cof*(y[3]*B[1] - y[4]*B[0]) ; // Az = a*(Vx*By - Vy*Bx)
97  }
98 
99  G4ThreeVector u(y[3], y[4], y[5]);
100  u *= inv_momentum_magnitude;
101 
102  G4ThreeVector BField(B[0],B[1],B[2]);
103 
104  G4double udb = anomaly*beta*gamma/(1.+gamma) * (BField * u);
105  G4double ucb = (anomaly+1./gamma)/beta;
106 
107  // Initialise the values of dydx that we do not update.
108  dydx[6] = dydx[7] = dydx[8] = 0.0;
109 
110  G4ThreeVector Spin(y[9],y[10],y[11]);
111 
112  G4double pcharge;
113  if (charge == 0.)
114  {
115  pcharge = 1.;
116  }
117  else
118  {
119  pcharge = charge;
120  }
121 
122  G4ThreeVector dSpin(0.,0.,0.);
123  if (Spin.mag2() != 0.)
124  {
125  dSpin = pcharge*omegac*(ucb*(Spin.cross(BField))-udb*(Spin.cross(u)));
126  }
127 
128  dydx[9] = dSpin.x();
129  dydx[10] = dSpin.y();
130  dydx[11] = dSpin.z();
131 
132  return;
133 }