Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Protected Member Functions | List of all members
G4MagHelicalStepper Class Referenceabstract

#include <G4MagHelicalStepper.hh>

Inheritance diagram for G4MagHelicalStepper:
Inheritance graph
[legend]
Collaboration diagram for G4MagHelicalStepper:
Collaboration graph
[legend]

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
 
- Public Member Functions inherited from G4MagIntegratorStepper
 G4MagIntegratorStepper (G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, bool isFSAL=false)
 
virtual ~G4MagIntegratorStepper ()
 
virtual void ComputeRightHandSide (const G4double y[], G4double dydx[])
 
void NormaliseTangentVector (G4double vec[6])
 
void NormalisePolarizationVector (G4double vec[12])
 
void RightHandSide (const double y[], double dydx[])
 
G4int GetNumberOfVariables () const
 
G4int GetNumberOfStateVariables () const
 
virtual G4int IntegratorOrder () const =0
 
G4int IntegrationOrder ()
 
G4EquationOfMotion * GetEquationOfMotion ()
 
void SetEquationOfMotion (G4EquationOfMotion *newEquation)
 
unsigned long GetfNoRHSCalls ()
 
void ResetfNORHSCalls ()
 
bool IsFSAL ()
 

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
 
- Protected Member Functions inherited from G4MagIntegratorStepper
void SetIntegrationOrder (int order)
 
void SetFSAL (bool flag=true)
 

Detailed Description

Definition at line 58 of file G4MagHelicalStepper.hh.

Constructor & Destructor Documentation

G4MagHelicalStepper::G4MagHelicalStepper ( G4Mag_EqRhs *  EqRhs)

Definition at line 47 of file G4MagHelicalStepper.cc.

48  : G4MagIntegratorStepper(EqRhs, 6), // integrate over 6 variables only !!
49  // position & velocity
50  fPtrMagEqOfMot(EqRhs), fAngCurve(0.), frCurve(0.), frHelix(0.)
51 {
52 }
G4MagIntegratorStepper::G4MagIntegratorStepper
G4MagIntegratorStepper(G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, bool isFSAL=false)
Definition: G4MagIntegratorStepper.cc:36
G4MagHelicalStepper::~G4MagHelicalStepper ( )
virtual

Definition at line 54 of file G4MagHelicalStepper.cc.

55 {
56 }

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.

64 {
65  // const G4int nvar = 6;
66 
67  // OLD const G4double approc_limit = 0.05;
68  // OLD approc_limit = 0.05 gives max.error=x^5/5!=(0.05)^5/5!=2.6*e-9
69  // NEW approc_limit = 0.005 gives max.error=x^5/5!=2.6*e-14
70 
71  const G4double approc_limit = 0.005;
72  G4ThreeVector Bnorm, B_x_P, vperp, vpar;
73 
74  G4double B_d_P;
75  G4double B_v_P;
76  G4double Theta;
77  G4double R_1;
78  G4double R_Helix;
79  G4double CosT2, SinT2, CosT, SinT;
80  G4ThreeVector positionMove, endTangent;
81 
82  G4double Bmag = Bfld.mag();
83  const G4double *pIn = yIn+3;
84  G4ThreeVector initVelocity= G4ThreeVector( pIn[0], pIn[1], pIn[2]);
85  G4double velocityVal = initVelocity.mag();
86  G4ThreeVector initTangent = (1.0/velocityVal) * initVelocity;
87 
88  R_1=GetInverseCurve(velocityVal,Bmag);
89 
90  // for too small magnetic fields there is no curvature
91  // (include momentum here) FIXME
92 
93  if( (std::fabs(R_1) < 1e-10)||(Bmag<1e-12) )
94  {
95  LinearStep( yIn, h, yHelix );
96 
97  // Store and/or calculate parameters for chord distance
98 
99  SetAngCurve(1.);
100  SetCurve(h);
101  SetRadHelix(0.);
102  }
103  else
104  {
105  Bnorm = (1.0/Bmag)*Bfld;
106 
107  // calculate the direction of the force
108 
109  B_x_P = Bnorm.cross(initTangent);
110 
111  // parallel and perp vectors
112 
113  B_d_P = Bnorm.dot(initTangent); // this is the fraction of P parallel to B
114 
115  vpar = B_d_P * Bnorm; // the component parallel to B
116  vperp= initTangent - vpar; // the component perpendicular to B
117 
118  B_v_P = std::sqrt( 1 - B_d_P * B_d_P); // Fraction of P perp to B
119 
120  // calculate the stepping angle
121 
122  Theta = R_1 * h; // * B_v_P;
123 
124  // Trigonometrix
125 
126  if( std::fabs(Theta) > approc_limit )
127  {
128  SinT = std::sin(Theta);
129  CosT = std::cos(Theta);
130  }
131  else
132  {
133  G4double Theta2 = Theta*Theta;
134  G4double Theta3 = Theta2 * Theta;
135  G4double Theta4 = Theta2 * Theta2;
136  SinT = Theta - 1.0/6.0 * Theta3;
137  CosT = 1 - 0.5 * Theta2 + 1.0/24.0 * Theta4;
138  }
139 
140  // the actual "rotation"
141 
142  G4double R = 1.0 / R_1;
143 
144  positionMove = R * ( SinT * vperp + (1-CosT) * B_x_P) + h * vpar;
145  endTangent = CosT * vperp + SinT * B_x_P + vpar;
146 
147  // Store the resulting position and tangent
148 
149  yHelix[0] = yIn[0] + positionMove.x();
150  yHelix[1] = yIn[1] + positionMove.y();
151  yHelix[2] = yIn[2] + positionMove.z();
152  yHelix[3] = velocityVal * endTangent.x();
153  yHelix[4] = velocityVal * endTangent.y();
154  yHelix[5] = velocityVal * endTangent.z();
155 
156  // Store 2*h step Helix if exist
157 
158  if(yHelix2)
159  {
160  SinT2 = 2.0 * SinT * CosT;
161  CosT2 = 1.0 - 2.0 * SinT * SinT;
162  endTangent = (CosT2 * vperp + SinT2 * B_x_P + vpar);
163  positionMove = R * ( SinT2 * vperp + (1-CosT2) * B_x_P) + h*2 * vpar;
164 
165  yHelix2[0] = yIn[0] + positionMove.x();
166  yHelix2[1] = yIn[1] + positionMove.y();
167  yHelix2[2] = yIn[2] + positionMove.z();
168  yHelix2[3] = velocityVal * endTangent.x();
169  yHelix2[4] = velocityVal * endTangent.y();
170  yHelix2[5] = velocityVal * endTangent.z();
171  }
172 
173  // Store and/or calculate parameters for chord distance
174 
175  G4double ptan=velocityVal*B_v_P;
176 
177  G4double particleCharge = fPtrMagEqOfMot->FCof() / (eplus*c_light);
178  R_Helix =std::abs( ptan/(fUnitConstant * particleCharge*Bmag));
179 
180  SetAngCurve(std::abs(Theta));
181  SetCurve(std::abs(R));
182  SetRadHelix(R_Helix);
183  }
184 }
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
G4MagHelicalStepper::LinearStep
void LinearStep(const G4double yIn[], G4double h, G4double yHelix[]) const
CLHEP::Hep3Vector::z
double z() const
G4MagHelicalStepper::GetInverseCurve
G4double GetInverseCurve(const G4double Momentum, const G4double Bmag)
G4Mag_EqRhs::FCof
G4double FCof() const
Definition: G4Mag_EqRhs.hh:84
eplus
static constexpr double eplus
Definition: G4SIunits.hh:199
G4MagHelicalStepper::SetRadHelix
void SetRadHelix(const G4double Rad)
G4MagHelicalStepper::SetAngCurve
void SetAngCurve(const G4double Ang)
G4MagHelicalStepper::SetCurve
void SetCurve(const G4double Curve)
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
G4double
double G4double
Definition: G4Types.hh:76
CLHEP::Hep3Vector::mag
double mag() const
python.hepunit.c_light
float c_light
Definition: hepunit.py:257

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4MagHelicalStepper::DistChord ( ) const
virtual

Implements G4MagIntegratorStepper.

Definition at line 238 of file G4MagHelicalStepper.cc.

239 {
240  // Check whether h/R > pi !!
241  // Method DistLine is good only for < pi
242 
243  G4double Ang=GetAngCurve();
244  if(Ang<=pi)
245  {
246  return GetRadHelix()*(1-std::cos(0.5*Ang));
247  }
248  else
249  {
250  if(Ang<twopi)
251  {
252  return GetRadHelix()*(1+std::cos(0.5*(twopi-Ang)));
253  }
254  else // return Diameter of projected circle
255  {
256  return 2*GetRadHelix();
257  }
258  }
259 }
G4MagHelicalStepper::GetRadHelix
G4double GetRadHelix() const
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4double
double G4double
Definition: G4Types.hh:76
G4MagHelicalStepper::GetAngCurve
G4double GetAngCurve() const

Here is the call graph for this function:

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

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

Here is the caller graph for this function:

G4double G4MagHelicalStepper::GetAngCurve ( ) const
inlineprotected

Here is the caller graph for this function:

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

Here is the caller graph for this function:

G4double G4MagHelicalStepper::GetRadHelix ( ) const
inlineprotected

Here is the caller graph for this function:

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

Here is the caller graph for this function:

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

Here is the caller graph for this function:

void G4MagHelicalStepper::SetAngCurve ( const G4double  Ang)
inlineprotected

Here is the caller graph for this function:

void G4MagHelicalStepper::SetCurve ( const G4double  Curve)
inlineprotected

Here is the caller graph for this function:

void G4MagHelicalStepper::SetRadHelix ( const G4double  Rad)
inlineprotected

Here is the caller graph for this function:

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

Implements G4MagIntegratorStepper.

Reimplemented in G4HelixMixedStepper, and G4ExactHelixStepper.

Definition at line 192 of file G4MagHelicalStepper.cc.

197 {
198  const G4int nvar = 6;
199 
200  G4int i;
201 
202  // correction for Richardson Extrapolation.
203  // G4double correction = 1. / ( (1 << IntegratorOrder()) -1 );
204 
205  G4double yTemp[8], yIn[8] ;
206  G4ThreeVector Bfld_initial, Bfld_midpoint;
207 
208  // Saving yInput because yInput and yOut can be aliases for same array
209 
210  for(i=0;i<nvar;i++) { yIn[i]=yInput[i]; }
211 
212  G4double h = hstep * 0.5;
213 
214  MagFieldEvaluate(yIn, Bfld_initial) ;
215 
216  // Do two half steps
217 
218  DumbStepper(yIn, Bfld_initial, h, yTemp);
219  MagFieldEvaluate(yTemp, Bfld_midpoint) ;
220  DumbStepper(yTemp, Bfld_midpoint, h, yOut);
221 
222  // Do a full Step
223 
224  h = hstep ;
225  DumbStepper(yIn, Bfld_initial, h, yTemp);
226 
227  // Error estimation
228 
229  for(i=0;i<nvar;i++)
230  {
231  yErr[i] = yOut[i] - yTemp[i] ;
232  }
233 
234  return;
235 }
G4MagHelicalStepper::DumbStepper
virtual void DumbStepper(const G4double y[], G4ThreeVector Bfld, G4double h, G4double yout[])=0
G4MagHelicalStepper::MagFieldEvaluate
void MagFieldEvaluate(const G4double y[], G4ThreeVector &Bfield)
G4int
int G4int
Definition: G4Types.hh:78
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

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