G4MagHelicalStepper Class Reference

#include <G4MagHelicalStepper.hh>

Inheritance diagram for G4MagHelicalStepper:

G4MagIntegratorStepper G4ExactHelixStepper G4HelixExplicitEuler G4HelixHeum G4HelixImplicitEuler G4HelixMixedStepper G4HelixSimpleRunge

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

Detailed Description

Definition at line 53 of file G4MagHelicalStepper.hh.

Constructor & Destructor Documentation

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]

Definition at line 54 of file G4MagHelicalStepper.cc.

00055 {
00056 }

Member Function Documentation

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().

00078 {
00079 return fAngCurve;  
00080 
00081 }

G4double G4MagHelicalStepper::GetCurve (  )  const [inline, protected]

Definition at line 88 of file G4MagHelicalStepper.icc.

00089 {
00090 return frCurve;  
00091 
00092 }

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().

00100 {
00101 return frHelix;  
00102 
00103 }

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]

Definition at line 71 of file G4MagHelicalStepper.icc.

Referenced by AdvanceHelix().

00072 {                                                
00073 fAngCurve=Ang; 
00074     
00075 }

void G4MagHelicalStepper::SetCurve ( const G4double  Curve  )  [inline, protected]

Definition at line 83 of file G4MagHelicalStepper.icc.

Referenced by AdvanceHelix().

00084 {
00085  frCurve=Curve;  
00086 }

void G4MagHelicalStepper::SetRadHelix ( const G4double  Rad  )  [inline, protected]

Definition at line 94 of file G4MagHelicalStepper.icc.

Referenced by AdvanceHelix().

00095 {
00096  frHelix=Rad;  
00097 }

virtual void G4MagHelicalStepper::Stepper ( const G4double  y[],
const G4double  dydx[],
G4double  h,
G4double  yout[],
G4double  yerr[] 
) [virtual]

Implements G4MagIntegratorStepper.

Reimplemented in G4ExactHelixStepper, and G4HelixMixedStepper.

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:27 2013 for Geant4 by  doxygen 1.4.7