00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include "G4Mag_SpinEqRhs.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4MagneticField.hh"
00043 #include "G4ThreeVector.hh"
00044
00045 G4Mag_SpinEqRhs::G4Mag_SpinEqRhs( G4MagneticField* MagField )
00046 : G4Mag_EqRhs( MagField ), omegac(0.), anomaly(0.0011659208),
00047 pcharge(0.), E(0.), gamma(0.), beta(0.)
00048 {
00049 }
00050
00051 G4Mag_SpinEqRhs::~G4Mag_SpinEqRhs()
00052 {
00053 }
00054
00055 void
00056 G4Mag_SpinEqRhs::SetChargeMomentumMass(G4double particleCharge,
00057 G4double MomentumXc,
00058 G4double particleMass)
00059 {
00060
00061 G4Mag_EqRhs::SetChargeMomentumMass(particleCharge, MomentumXc, particleMass);
00062
00063 omegac = (eplus/particleMass)*c_light;
00064
00065 pcharge = particleCharge;
00066
00067 E = std::sqrt(sqr(MomentumXc)+sqr(particleMass));
00068 beta = MomentumXc/E;
00069 gamma = E/particleMass;
00070
00071 G4double neutronAnomaly = -2.9156797;
00072 if (pcharge==0.) SetAnomaly(neutronAnomaly);
00073 }
00074
00075 void
00076 G4Mag_SpinEqRhs::EvaluateRhsGivenB( const G4double y[],
00077 const G4double B[3],
00078 G4double dydx[] ) const
00079 {
00080 G4double momentum_mag_square = sqr(y[3]) + sqr(y[4]) + sqr(y[5]);
00081 G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square );
00082 G4double cof = FCof()*inv_momentum_magnitude;
00083
00084 dydx[0] = y[3] * inv_momentum_magnitude;
00085 dydx[1] = y[4] * inv_momentum_magnitude;
00086 dydx[2] = y[5] * inv_momentum_magnitude;
00087
00088 if (pcharge == 0.) {
00089 dydx[3] = 0.;
00090 dydx[4] = 0.;
00091 dydx[5] = 0.;
00092 } else {
00093 dydx[3] = cof*(y[4]*B[2] - y[5]*B[1]) ;
00094 dydx[4] = cof*(y[5]*B[0] - y[3]*B[2]) ;
00095 dydx[5] = cof*(y[3]*B[1] - y[4]*B[0]) ;
00096 }
00097
00098 G4ThreeVector u(y[3], y[4], y[5]);
00099 u *= inv_momentum_magnitude;
00100
00101 G4ThreeVector BField(B[0],B[1],B[2]);
00102
00103 G4double udb = anomaly*beta*gamma/(1.+gamma) * (BField * u);
00104 G4double ucb = (anomaly+1./gamma)/beta;
00105
00106
00107 dydx[6] = dydx[7] = dydx[8] = 0.0;
00108
00109 G4ThreeVector Spin(y[9],y[10],y[11]);
00110
00111 G4ThreeVector dSpin;
00112
00113 if (pcharge == 0.) {
00114
00115 dSpin = omegac*(ucb*(Spin.cross(BField))-udb*(Spin.cross(u)));
00116 } else {
00117 dSpin = pcharge*omegac*(ucb*(Spin.cross(BField))-udb*(Spin.cross(u)));
00118 }
00119
00120 dydx[ 9] = dSpin.x();
00121 dydx[10] = dSpin.y();
00122 dydx[11] = dSpin.z();
00123
00124 return ;
00125 }