G4MagHelicalStepper Class Reference

G4MagHelicalStepper is an abstract base class for integrator of particle's equation of motion, used in tracking in space dependent magnetic field, and for a set of steppers which use the helix as 'first order' solution. More...

#include <G4MagHelicalStepper.hh>

Inheritance diagram for G4MagHelicalStepper:

Public Member Functions
	G4MagHelicalStepper (G4Mag_EqRhs *EqRhs)
	~G4MagHelicalStepper () override=default
	G4MagHelicalStepper (const G4MagHelicalStepper &)=delete
G4MagHelicalStepper &	operator= (const G4MagHelicalStepper &)=delete
void	Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[]) override
virtual void	DumbStepper (const G4double y[], G4ThreeVector Bfld, G4double h, G4double yout[])=0
G4double	DistChord () const override
Public Member Functions inherited from G4MagIntegratorStepper
	G4MagIntegratorStepper (G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, G4bool isFSAL=false)
virtual	~G4MagIntegratorStepper ()=default
	G4MagIntegratorStepper (const G4MagIntegratorStepper &)=delete
G4MagIntegratorStepper &	operator= (const G4MagIntegratorStepper &)=delete
void	NormaliseTangentVector (G4double vec[6])
void	NormalisePolarizationVector (G4double vec[12])
void	RightHandSide (const G4double y[], G4double dydx[]) const
void	RightHandSide (const G4double y[], G4double dydx[], G4double field[]) const
G4int	GetNumberOfVariables () const
G4int	GetNumberOfStateVariables () const
virtual G4int	IntegratorOrder () const =0
virtual G4StepperType	StepperType () const
G4int	IntegrationOrder ()
G4EquationOfMotion *	GetEquationOfMotion ()
const G4EquationOfMotion *	GetEquationOfMotion () const
void	SetEquationOfMotion (G4EquationOfMotion *newEquation)
unsigned long	GetfNoRHSCalls ()
void	ResetfNORHSCalls ()
G4bool	IsFSAL () const
G4bool	isQSS () const
void	SetIsQSS (G4bool val)

Protected Member Functions
void	LinearStep (const G4double yIn[], G4double h, G4double yHelix[]) const
void	AdvanceHelix (const G4double yIn[], const G4ThreeVector &Bfld, G4double h, G4double yHelix[], G4double yHelix2[]=nullptr)
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
Protected Member Functions inherited from G4MagIntegratorStepper
void	SetIntegrationOrder (G4int order)
void	SetFSAL (G4bool flag=true)

Detailed Description

G4MagHelicalStepper is an abstract base class for integrator of particle's equation of motion, used in tracking in space dependent magnetic field, and for a set of steppers which use the helix as 'first order' solution.

Definition at line 57 of file G4MagHelicalStepper.hh.

Constructor & Destructor Documentation

G4MagHelicalStepper() [1/2]

G4MagHelicalStepper::G4MagHelicalStepper

(

G4Mag_EqRhs *

EqRhs

)

Constructor for G4MagHelicalStepper.

Parameters

[in]

EqRhs

Pointer to the provided equation of motion.

Definition at line 45 of file G4MagHelicalStepper.cc.

   : G4MagIntegratorStepper(EqRhs, 6), // integrate over 6 variables only !!
                                       // position & velocity
     fPtrMagEqOfMot(EqRhs)
{
}

Referenced by G4ExactHelixStepper::G4ExactHelixStepper(), G4HelixExplicitEuler::G4HelixExplicitEuler(), G4HelixHeum::G4HelixHeum(), G4HelixImplicitEuler::G4HelixImplicitEuler(), G4HelixMixedStepper::G4HelixMixedStepper(), G4HelixSimpleRunge::G4HelixSimpleRunge(), G4MagHelicalStepper(), and operator=().

~G4MagHelicalStepper()

G4MagHelicalStepper::~G4MagHelicalStepper ( )

overridedefault

Default Destructor.

G4MagHelicalStepper() [2/2]

G4MagHelicalStepper::G4MagHelicalStepper ( const G4MagHelicalStepper & )

delete

Copy constructor and assignment operator not allowed.

Member Function Documentation

AdvanceHelix()

void G4MagHelicalStepper::AdvanceHelix ( const G4double yIn[], const G4ThreeVector & Bfld, G4double h, G4double yHelix[], G4double yHelix2[] = nullptr )

protected

A first order Step along a helix inside the field.

Definition at line 53 of file G4MagHelicalStepper.cc.

{
  // const G4int    nvar = 6;
 
  // OLD  const G4double approc_limit = 0.05;
  // OLD  approc_limit = 0.05 gives max.error=x^5/5!=(0.05)^5/5!=2.6*e-9
  // NEW  approc_limit = 0.005 gives max.error=x^5/5!=2.6*e-14
 
  const G4double approc_limit = 0.005;
  G4ThreeVector  Bnorm, B_x_P, vperp, vpar;
 
  G4double B_d_P;
  G4double B_v_P;
  G4double Theta;
  G4double R_1;
  G4double R_Helix;
  G4double CosT2, SinT2, CosT, SinT;
  G4ThreeVector positionMove, endTangent;
 
  G4double Bmag = Bfld.mag();
  const G4double* pIn = yIn+3;
  G4ThreeVector initVelocity = G4ThreeVector( pIn[0], pIn[1], pIn[2]);
  G4double      velocityVal = initVelocity.mag();
  G4ThreeVector initTangent = (1.0/velocityVal) * initVelocity;
  
  R_1 = GetInverseCurve(velocityVal,Bmag);
 
  // for too small magnetic fields there is no curvature
  // (include momentum here) FIXME
 
  if( (std::fabs(R_1) < 1e-10)||(Bmag<1e-12) )
  {
    LinearStep( yIn, h, yHelix );
 
    // Store and/or calculate parameters for chord distance
 
    SetAngCurve(1.);     
    SetCurve(h);
    SetRadHelix(0.);
  }
  else
  {
    Bnorm = (1.0/Bmag)*Bfld;
 
    // calculate the direction of the force
    
    B_x_P = Bnorm.cross(initTangent);
 
    // parallel and perp vectors
 
    B_d_P = Bnorm.dot(initTangent); // this is the fraction of P parallel to B
 
    vpar = B_d_P * Bnorm;       // the component parallel      to B
    vperp= initTangent - vpar;  // the component perpendicular to B
    
    B_v_P  = std::sqrt( 1 - B_d_P * B_d_P); // Fraction of P perp to B
 
    // calculate  the stepping angle
  
    Theta   = R_1 * h; // * B_v_P;
 
    // Trigonometrix
      
    if( std::fabs(Theta) > approc_limit )
    {
       SinT     = std::sin(Theta);
       CosT     = std::cos(Theta);
    }
    else
    {
      G4double Theta2 = Theta*Theta;
      G4double Theta3 = Theta2 * Theta;
      G4double Theta4 = Theta2 * Theta2;
      SinT     = Theta - 1.0/6.0 * Theta3;
      CosT     = 1 - 0.5 * Theta2 + 1.0/24.0 * Theta4;
    }
 
    // the actual "rotation"
 
    G4double R = 1.0 / R_1;
 
    positionMove  = R * ( SinT * vperp + (1-CosT) * B_x_P) + h * vpar;
    endTangent    = CosT * vperp + SinT * B_x_P + vpar;
 
    // Store the resulting position and tangent
 
    yHelix[0]   = yIn[0] + positionMove.x(); 
    yHelix[1]   = yIn[1] + positionMove.y(); 
    yHelix[2]   = yIn[2] + positionMove.z();
    yHelix[3] = velocityVal * endTangent.x();
    yHelix[4] = velocityVal * endTangent.y();
    yHelix[5] = velocityVal * endTangent.z();
 
    // Store 2*h step Helix if exist
 
    if(yHelix2 != nullptr)
    {
      SinT2     = 2.0 * SinT * CosT;
      CosT2     = 1.0 - 2.0 * SinT * SinT;
      endTangent    = (CosT2 * vperp + SinT2 * B_x_P + vpar);
      positionMove  = R * ( SinT2 * vperp + (1-CosT2) * B_x_P) + h*2 * vpar;
      
      yHelix2[0]   = yIn[0] + positionMove.x(); 
      yHelix2[1]   = yIn[1] + positionMove.y(); 
      yHelix2[2]   = yIn[2] + positionMove.z(); 
      yHelix2[3] = velocityVal * endTangent.x();
      yHelix2[4] = velocityVal * endTangent.y();
      yHelix2[5] = velocityVal * endTangent.z();
    }
 
    // Store and/or calculate parameters for chord distance
 
    G4double ptan=velocityVal*B_v_P;
 
    G4double particleCharge = fPtrMagEqOfMot->FCof() / (eplus*c_light); 
    R_Helix =std::abs( ptan/(fUnitConstant  * particleCharge*Bmag));
       
    SetAngCurve(std::abs(Theta));
    SetCurve(std::abs(R));
    SetRadHelix(R_Helix);
  }
}

Referenced by G4ExactHelixStepper::DumbStepper(), G4HelixExplicitEuler::DumbStepper(), G4HelixHeum::DumbStepper(), G4HelixImplicitEuler::DumbStepper(), G4HelixMixedStepper::DumbStepper(), G4HelixSimpleRunge::DumbStepper(), G4ExactHelixStepper::Stepper(), G4HelixExplicitEuler::Stepper(), and G4HelixMixedStepper::Stepper().

DistChord()

G4double G4MagHelicalStepper::DistChord ( ) const

overridevirtual

Estimates the maximum distance of curved solution and chord.

Implements G4MagIntegratorStepper.

Definition at line 231 of file G4MagHelicalStepper.cc.

{
  // Check whether h/R >  pi  !!
  // Method DistLine is good only for <  pi
 
  G4double Ang=GetAngCurve();
  if(Ang<=pi)
  {
    return GetRadHelix()*(1-std::cos(0.5*Ang));
  }
  
  if(Ang<twopi)
  {
    return GetRadHelix()*(1+std::cos(0.5*(twopi-Ang)));
  }
 
  // return Diameter of projected circle
  return 2*GetRadHelix();
}

DumbStepper()

virtual void G4MagHelicalStepper::DumbStepper ( const G4double y[], G4ThreeVector Bfld, G4double h, G4double yout[] )

pure virtual

Same as Stepper() function above, but should perform a 'dump' step without error calculation. To be implemented in concrete derived classes.

Implemented in G4ExactHelixStepper, G4HelixExplicitEuler, G4HelixHeum, G4HelixImplicitEuler, G4HelixMixedStepper, and G4HelixSimpleRunge.

Referenced by Stepper().

GetAngCurve()

G4double G4MagHelicalStepper::GetAngCurve ( ) const

inlineprotected

Referenced by G4ExactHelixStepper::DistChord(), G4HelixExplicitEuler::DistChord(), G4HelixMixedStepper::DistChord(), DistChord(), and G4HelixExplicitEuler::Stepper().

GetCurve()

G4double G4MagHelicalStepper::GetCurve ( ) const

inlineprotected

GetInverseCurve()

G4double G4MagHelicalStepper::GetInverseCurve ( const G4double Momentum, const G4double Bmag )

inlineprotected

Evaluates inverse of Curvature of Track.

Referenced by AdvanceHelix(), and G4HelixMixedStepper::Stepper().

GetRadHelix()

G4double G4MagHelicalStepper::GetRadHelix ( ) const

inlineprotected

Referenced by G4ExactHelixStepper::DistChord(), G4HelixExplicitEuler::DistChord(), G4HelixMixedStepper::DistChord(), and DistChord().

LinearStep()

void G4MagHelicalStepper::LinearStep ( const G4double yIn[], G4double h, G4double yHelix[] ) const

inlineprotected

Performs a linear Step in regions without magnetic field.

Referenced by AdvanceHelix().

MagFieldEvaluate()

void G4MagHelicalStepper::MagFieldEvaluate ( const G4double y[], G4ThreeVector & Bfield )

inlineprotected

Evaluates the field at a certain point.

Referenced by G4HelixHeum::DumbStepper(), G4HelixImplicitEuler::DumbStepper(), G4HelixSimpleRunge::DumbStepper(), G4ExactHelixStepper::Stepper(), G4HelixExplicitEuler::Stepper(), G4HelixMixedStepper::Stepper(), and Stepper().

operator=()

G4MagHelicalStepper & G4MagHelicalStepper::operator= ( const G4MagHelicalStepper & )

delete

SetAngCurve()

void G4MagHelicalStepper::SetAngCurve ( const G4double Ang )

inlineprotected

Modifiers and accessors for storing and using the parameters of a track: radius of curve, Stepping angle, Radius of projected helix.

Referenced by AdvanceHelix(), G4HelixExplicitEuler::Stepper(), and G4HelixMixedStepper::Stepper().

SetCurve()

void G4MagHelicalStepper::SetCurve ( const G4double Curve )

inlineprotected

Referenced by AdvanceHelix(), and G4HelixMixedStepper::Stepper().

SetRadHelix()

void G4MagHelicalStepper::SetRadHelix ( const G4double Rad )

inlineprotected

Referenced by AdvanceHelix().

Stepper()

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

overridevirtual

The stepper for the Runge Kutta integration. The stepsize is fixed, with the step size given by 'h'. Integrates ODE starting values y[0 to 6]. Outputs yout[] and its estimated error yerr[].

Parameters

[in]	y	Starting values array of integration variables.
[in]	dydx	Derivatives array.
[in]	h	The given step size.
[out]	yout	Integration output.
[out]	yerr	The estimated error.

Implements G4MagIntegratorStepper.

Definition at line 184 of file G4MagHelicalStepper.cc.

{  
   const G4int nvar = 6;
 
   // correction for Richardson Extrapolation.
   // G4double  correction = 1. / ( (1 << IntegratorOrder()) -1 );
   
   G4double      yTemp[8], yIn[8] ;
   G4ThreeVector Bfld_initial, Bfld_midpoint;
   
   //  Saving yInput because yInput and yOut can be aliases for same array
   //
   for(G4int i=0; i<nvar; ++i)
   {
     yIn[i]=yInput[i];
   }
 
   G4double h = hstep * 0.5; 
 
   MagFieldEvaluate(yIn, Bfld_initial) ;      
 
   // Do two half steps
   //
   DumbStepper(yIn, Bfld_initial, h, yTemp);
   MagFieldEvaluate(yTemp, Bfld_midpoint) ;     
   DumbStepper(yTemp, Bfld_midpoint, h, yOut); 
 
   // Do a full Step
   //
   h = hstep ;
   DumbStepper(yIn, Bfld_initial, h, yTemp);
 
   // Error estimation
   //
   for(G4int i=0; i<nvar; ++i)
   {
     yErr[i] = yOut[i] - yTemp[i] ;
   }
   
   return;
}

---

The documentation for this class was generated from the following files:

• G4MagHelicalStepper.hh
• G4MagHelicalStepper.cc