#include <G4MagHelicalStepper.hh>
Inheritance diagram for G4MagHelicalStepper:
Public Member Functions | |
G4MagHelicalStepper (G4Mag_EqRhs *EqRhs) | |
virtual | ~G4MagHelicalStepper () |
virtual void | Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[]) |
virtual void | DumbStepper (const G4double y[], G4ThreeVector Bfld, G4double h, G4double yout[])=0 |
G4double | DistChord () const |
Protected Member Functions | |
void | LinearStep (const G4double yIn[], G4double h, G4double yHelix[]) const |
void | AdvanceHelix (const G4double yIn[], G4ThreeVector Bfld, G4double h, G4double yHelix[], G4double yHelix2[]=0) |
void | MagFieldEvaluate (const G4double y[], G4ThreeVector &Bfield) |
G4double | GetInverseCurve (const G4double Momentum, const G4double Bmag) |
void | SetAngCurve (const G4double Ang) |
G4double | GetAngCurve () const |
void | SetCurve (const G4double Curve) |
G4double | GetCurve () const |
void | SetRadHelix (const G4double Rad) |
G4double | GetRadHelix () const |
Definition at line 53 of file G4MagHelicalStepper.hh.
G4MagHelicalStepper::G4MagHelicalStepper | ( | G4Mag_EqRhs * | EqRhs | ) |
Definition at line 47 of file G4MagHelicalStepper.cc.
00048 : G4MagIntegratorStepper(EqRhs, 6), // integrate over 6 variables only !! 00049 // position & velocity 00050 fPtrMagEqOfMot(EqRhs), fAngCurve(0.), frCurve(0.), frHelix(0.) 00051 { 00052 }
G4MagHelicalStepper::~G4MagHelicalStepper | ( | ) | [virtual] |
void G4MagHelicalStepper::AdvanceHelix | ( | const G4double | yIn[], | |
G4ThreeVector | Bfld, | |||
G4double | h, | |||
G4double | yHelix[], | |||
G4double | yHelix2[] = 0 | |||
) | [protected] |
Definition at line 59 of file G4MagHelicalStepper.cc.
References G4Mag_EqRhs::FCof(), GetInverseCurve(), LinearStep(), SetAngCurve(), SetCurve(), and SetRadHelix().
Referenced by G4HelixSimpleRunge::DumbStepper(), G4HelixMixedStepper::DumbStepper(), G4HelixImplicitEuler::DumbStepper(), G4HelixHeum::DumbStepper(), G4HelixExplicitEuler::DumbStepper(), and G4ExactHelixStepper::DumbStepper().
00064 { 00065 // const G4int nvar = 6; 00066 00067 // OLD const G4double approc_limit = 0.05; 00068 // OLD approc_limit = 0.05 gives max.error=x^5/5!=(0.05)^5/5!=2.6*e-9 00069 // NEW approc_limit = 0.005 gives max.error=x^5/5!=2.6*e-14 00070 00071 const G4double approc_limit = 0.005; 00072 G4ThreeVector Bnorm, B_x_P, vperp, vpar; 00073 00074 G4double B_d_P; 00075 G4double B_v_P; 00076 G4double Theta; 00077 G4double R_1; 00078 G4double R_Helix; 00079 G4double CosT2, SinT2, CosT, SinT; 00080 G4ThreeVector positionMove, endTangent; 00081 00082 G4double Bmag = Bfld.mag(); 00083 const G4double *pIn = yIn+3; 00084 G4ThreeVector initVelocity= G4ThreeVector( pIn[0], pIn[1], pIn[2]); 00085 G4double velocityVal = initVelocity.mag(); 00086 G4ThreeVector initTangent = (1.0/velocityVal) * initVelocity; 00087 00088 R_1=GetInverseCurve(velocityVal,Bmag); 00089 00090 // for too small magnetic fields there is no curvature 00091 // (include momentum here) FIXME 00092 00093 if( (std::fabs(R_1) < 1e-10)||(Bmag<1e-12) ) 00094 { 00095 LinearStep( yIn, h, yHelix ); 00096 00097 // Store and/or calculate parameters for chord distance 00098 00099 SetAngCurve(1.); 00100 SetCurve(h); 00101 SetRadHelix(0.); 00102 } 00103 else 00104 { 00105 Bnorm = (1.0/Bmag)*Bfld; 00106 00107 // calculate the direction of the force 00108 00109 B_x_P = Bnorm.cross(initTangent); 00110 00111 // parallel and perp vectors 00112 00113 B_d_P = Bnorm.dot(initTangent); // this is the fraction of P parallel to B 00114 00115 vpar = B_d_P * Bnorm; // the component parallel to B 00116 vperp= initTangent - vpar; // the component perpendicular to B 00117 00118 B_v_P = std::sqrt( 1 - B_d_P * B_d_P); // Fraction of P perp to B 00119 00120 // calculate the stepping angle 00121 00122 Theta = R_1 * h; // * B_v_P; 00123 00124 // Trigonometrix 00125 00126 if( std::fabs(Theta) > approc_limit ) 00127 { 00128 SinT = std::sin(Theta); 00129 CosT = std::cos(Theta); 00130 } 00131 else 00132 { 00133 G4double Theta2 = Theta*Theta; 00134 G4double Theta3 = Theta2 * Theta; 00135 G4double Theta4 = Theta2 * Theta2; 00136 SinT = Theta - 1.0/6.0 * Theta3; 00137 CosT = 1 - 0.5 * Theta2 + 1.0/24.0 * Theta4; 00138 } 00139 00140 // the actual "rotation" 00141 00142 G4double R = 1.0 / R_1; 00143 00144 positionMove = R * ( SinT * vperp + (1-CosT) * B_x_P) + h * vpar; 00145 endTangent = CosT * vperp + SinT * B_x_P + vpar; 00146 00147 // Store the resulting position and tangent 00148 00149 yHelix[0] = yIn[0] + positionMove.x(); 00150 yHelix[1] = yIn[1] + positionMove.y(); 00151 yHelix[2] = yIn[2] + positionMove.z(); 00152 yHelix[3] = velocityVal * endTangent.x(); 00153 yHelix[4] = velocityVal * endTangent.y(); 00154 yHelix[5] = velocityVal * endTangent.z(); 00155 00156 // Store 2*h step Helix if exist 00157 00158 if(yHelix2) 00159 { 00160 SinT2 = 2.0 * SinT * CosT; 00161 CosT2 = 1.0 - 2.0 * SinT * SinT; 00162 endTangent = (CosT2 * vperp + SinT2 * B_x_P + vpar); 00163 positionMove = R * ( SinT2 * vperp + (1-CosT2) * B_x_P) + h*2 * vpar; 00164 00165 yHelix2[0] = yIn[0] + positionMove.x(); 00166 yHelix2[1] = yIn[1] + positionMove.y(); 00167 yHelix2[2] = yIn[2] + positionMove.z(); 00168 yHelix2[3] = velocityVal * endTangent.x(); 00169 yHelix2[4] = velocityVal * endTangent.y(); 00170 yHelix2[5] = velocityVal * endTangent.z(); 00171 } 00172 00173 // Store and/or calculate parameters for chord distance 00174 00175 G4double ptan=velocityVal*B_v_P; 00176 00177 G4double particleCharge = fPtrMagEqOfMot->FCof() / (eplus*c_light); 00178 R_Helix =std::abs( ptan/(fUnitConstant * particleCharge*Bmag)); 00179 00180 SetAngCurve(std::abs(Theta)); 00181 SetCurve(std::abs(R)); 00182 SetRadHelix(R_Helix); 00183 } 00184 }
G4double G4MagHelicalStepper::DistChord | ( | ) | const [virtual] |
Implements G4MagIntegratorStepper.
Reimplemented in G4ExactHelixStepper, G4HelixExplicitEuler, and G4HelixMixedStepper.
Definition at line 238 of file G4MagHelicalStepper.cc.
References GetAngCurve(), GetRadHelix(), and G4INCL::Math::pi.
00239 { 00240 // Check whether h/R > pi !! 00241 // Method DistLine is good only for < pi 00242 00243 G4double Ang=GetAngCurve(); 00244 if(Ang<=pi) 00245 { 00246 return GetRadHelix()*(1-std::cos(0.5*Ang)); 00247 } 00248 else 00249 { 00250 if(Ang<twopi) 00251 { 00252 return GetRadHelix()*(1+std::cos(0.5*(twopi-Ang))); 00253 } 00254 else // return Diameter of projected circle 00255 { 00256 return 2*GetRadHelix(); 00257 } 00258 } 00259 }
virtual void G4MagHelicalStepper::DumbStepper | ( | const G4double | y[], | |
G4ThreeVector | Bfld, | |||
G4double | h, | |||
G4double | yout[] | |||
) | [pure virtual] |
Implemented in G4ExactHelixStepper, G4HelixExplicitEuler, G4HelixHeum, G4HelixImplicitEuler, G4HelixMixedStepper, and G4HelixSimpleRunge.
G4double G4MagHelicalStepper::GetAngCurve | ( | ) | const [inline, protected] |
Definition at line 77 of file G4MagHelicalStepper.icc.
Referenced by DistChord(), G4HelixMixedStepper::DistChord(), G4HelixExplicitEuler::DistChord(), and G4ExactHelixStepper::DistChord().
G4double G4MagHelicalStepper::GetCurve | ( | ) | const [inline, protected] |
G4double G4MagHelicalStepper::GetInverseCurve | ( | const G4double | Momentum, | |
const G4double | Bmag | |||
) | [inline, protected] |
Definition at line 61 of file G4MagHelicalStepper.icc.
References G4Mag_EqRhs::FCof().
Referenced by AdvanceHelix().
00063 { 00064 00065 G4double inv_momentum = 1.0 / Momentum ; 00066 G4double particleCharge = fPtrMagEqOfMot->FCof() / (CLHEP::eplus*CLHEP::c_light); 00067 G4double fCoefficient = -fUnitConstant * particleCharge*inv_momentum; 00068 00069 return fCoefficient*Bmag; 00070 }
G4double G4MagHelicalStepper::GetRadHelix | ( | ) | const [inline, protected] |
Definition at line 99 of file G4MagHelicalStepper.icc.
Referenced by DistChord(), G4HelixMixedStepper::DistChord(), G4HelixExplicitEuler::DistChord(), and G4ExactHelixStepper::DistChord().
void G4MagHelicalStepper::LinearStep | ( | const G4double | yIn[], | |
G4double | h, | |||
G4double | yHelix[] | |||
) | const [inline, protected] |
Definition at line 34 of file G4MagHelicalStepper.icc.
Referenced by AdvanceHelix().
00037 { 00038 G4double momentum_val = std::sqrt(yIn[3]*yIn[3] + yIn[4]*yIn[4] + yIn[5]*yIn[5]) ; 00039 G4double inv_momentum = 1.0 / momentum_val ; 00040 G4double yDir[3]; 00041 // G4double h_div_momentum = 1.0 / momentum_val ; 00042 00043 for( G4int i = 0; i < 3; i++ ) { 00044 yDir[i] = inv_momentum * yIn[i+3]; 00045 yLinear[i] = yIn[i] + h * yDir[i]; 00046 // yLinear[i] = yIn[i] + h_div_momentum * yIn[i+3]; 00047 yLinear[i+3] = yIn[i+3]; 00048 } 00049 }
void G4MagHelicalStepper::MagFieldEvaluate | ( | const G4double | y[], | |
G4ThreeVector & | Bfield | |||
) | [inline, protected] |
Definition at line 52 of file G4MagHelicalStepper.icc.
References G4MagIntegratorStepper::GetEquationOfMotion().
Referenced by G4HelixSimpleRunge::DumbStepper(), G4HelixImplicitEuler::DumbStepper(), and G4HelixHeum::DumbStepper().
00054 { 00055 G4double B[3]; 00056 GetEquationOfMotion()-> GetFieldValue(y, B); 00057 Bfield= G4ThreeVector( B[0], B[1], B[2] ); 00058 }
void G4MagHelicalStepper::SetAngCurve | ( | const G4double | Ang | ) | [inline, protected] |
void G4MagHelicalStepper::SetCurve | ( | const G4double | Curve | ) | [inline, protected] |
void G4MagHelicalStepper::SetRadHelix | ( | const G4double | Rad | ) | [inline, protected] |