Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Protected Member Functions | List of all members
G4MagHelicalStepper Class Referenceabstract

#include <G4MagHelicalStepper.hh>

+ Inheritance diagram for G4MagHelicalStepper:

Public Member Functions

 G4MagHelicalStepper (G4Mag_EqRhs *EqRhs)
 
virtual ~G4MagHelicalStepper ()
 
 G4MagHelicalStepper (const G4MagHelicalStepper &)=delete
 
G4MagHelicalStepper & operator= (const G4MagHelicalStepper &)=delete
 
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, G4bool isFSAL=false)
 
virtual ~G4MagIntegratorStepper ()=default
 
 G4MagIntegratorStepper (const G4MagIntegratorStepper &)=delete
 
G4MagIntegratorStepper & operator= (const G4MagIntegratorStepper &)=delete
 
virtual void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])=0
 
virtual G4double DistChord () const =0
 
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
 
G4int IntegrationOrder ()
 
G4EquationOfMotion * GetEquationOfMotion ()
 
const G4EquationOfMotion * GetEquationOfMotion () const
 
void SetEquationOfMotion (G4EquationOfMotion *newEquation)
 
unsigned long GetfNoRHSCalls ()
 
void ResetfNORHSCalls ()
 
G4bool IsFSAL () 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
 
- Protected Member Functions inherited from G4MagIntegratorStepper
void SetIntegrationOrder (G4int order)
 
void SetFSAL (G4bool flag=true)
 

Detailed Description

Definition at line 50 of file G4MagHelicalStepper.hh.

Constructor & Destructor Documentation

◆ G4MagHelicalStepper() [1/2]

G4MagHelicalStepper::G4MagHelicalStepper ( G4Mag_EqRhs *  EqRhs)

Definition at line 45 of file G4MagHelicalStepper.cc.

46 : G4MagIntegratorStepper(EqRhs, 6), // integrate over 6 variables only !!
47 // position & velocity
48 fPtrMagEqOfMot(EqRhs)
49{
50}

◆ ~G4MagHelicalStepper()

G4MagHelicalStepper::~G4MagHelicalStepper ( )
virtual

Definition at line 52 of file G4MagHelicalStepper.cc.

53{
54}

◆ G4MagHelicalStepper() [2/2]

G4MagHelicalStepper::G4MagHelicalStepper ( const G4MagHelicalStepper &  )
delete

Member Function Documentation

◆ AdvanceHelix()

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

Definition at line 57 of file G4MagHelicalStepper.cc.

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

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

◆ DistChord()

G4double G4MagHelicalStepper::DistChord ( ) const
virtual

Implements G4MagIntegratorStepper.

Definition at line 235 of file G4MagHelicalStepper.cc.

236{
237 // Check whether h/R > pi !!
238 // Method DistLine is good only for < pi
239
240 G4double Ang=GetAngCurve();
241 if(Ang<=pi)
242 {
243 return GetRadHelix()*(1-std::cos(0.5*Ang));
244 }
245 else
246 {
247 if(Ang<twopi)
248 {
249 return GetRadHelix()*(1+std::cos(0.5*(twopi-Ang)));
250 }
251 else // return Diameter of projected circle
252 {
253 return 2*GetRadHelix();
254 }
255 }
256}
G4MagHelicalStepper::GetRadHelix
G4double GetRadHelix() const
G4MagHelicalStepper::GetAngCurve
G4double GetAngCurve() const

◆ DumbStepper()

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

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

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

Referenced by AdvanceHelix().

◆ MagFieldEvaluate()

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

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

◆ operator=()

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

◆ SetAngCurve()

void G4MagHelicalStepper::SetAngCurve ( const G4double  Ang)
inlineprotected

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[] 
)
virtual

Implements G4MagIntegratorStepper.

Reimplemented in G4ExactHelixStepper, and G4HelixMixedStepper.

Definition at line 188 of file G4MagHelicalStepper.cc.

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