Geant4 11.2.2
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)
 
 ~G4MagHelicalStepper () override
 
 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
 
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

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}
G4MagIntegratorStepper::G4MagIntegratorStepper
G4MagIntegratorStepper(G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, G4bool isFSAL=false)
Definition G4MagIntegratorStepper.cc:38

◆ ~G4MagHelicalStepper()

G4MagHelicalStepper::~G4MagHelicalStepper ( )
overridedefault

◆ G4MagHelicalStepper() [2/2]

G4MagHelicalStepper::G4MagHelicalStepper ( const G4MagHelicalStepper & )
delete

Member Function Documentation

◆ AdvanceHelix()

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

Definition at line 55 of file G4MagHelicalStepper.cc.

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

◆ DistChord()

G4double G4MagHelicalStepper::DistChord ( ) const
overridevirtual

Implements G4MagIntegratorStepper.

Definition at line 233 of file G4MagHelicalStepper.cc.

234{
235 // Check whether h/R > pi !!
236 // Method DistLine is good only for < pi
237
238 G4double Ang=GetAngCurve();
239 if(Ang<=pi)
240 {
241 return GetRadHelix()*(1-std::cos(0.5*Ang));
242 }
243
244 if(Ang<twopi)
245 {
246 return GetRadHelix()*(1+std::cos(0.5*(twopi-Ang)));
247 }
248 else // return Diameter of projected circle
249 {
250 return 2*GetRadHelix();
251 }
252}
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(), G4ExactHelixStepper::Stepper(), G4HelixExplicitEuler::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[] )
overridevirtual

Implements G4MagIntegratorStepper.

Definition at line 186 of file G4MagHelicalStepper.cc.

191{
192 const G4int nvar = 6;
193
194 // correction for Richardson Extrapolation.
195 // G4double correction = 1. / ( (1 << IntegratorOrder()) -1 );
196
197 G4double yTemp[8], yIn[8] ;
198 G4ThreeVector Bfld_initial, Bfld_midpoint;
199
200 // Saving yInput because yInput and yOut can be aliases for same array
201 //
202 for(G4int i=0; i<nvar; ++i)
203 {
204 yIn[i]=yInput[i];
205 }
206
207 G4double h = hstep * 0.5;
208
209 MagFieldEvaluate(yIn, Bfld_initial) ;
210
211 // Do two half steps
212 //
213 DumbStepper(yIn, Bfld_initial, h, yTemp);
214 MagFieldEvaluate(yTemp, Bfld_midpoint) ;
215 DumbStepper(yTemp, Bfld_midpoint, h, yOut);
216
217 // Do a full Step
218 //
219 h = hstep ;
220 DumbStepper(yIn, Bfld_initial, h, yTemp);
221
222 // Error estimation
223 //
224 for(G4int i=0; i<nvar; ++i)
225 {
226 yErr[i] = yOut[i] - yTemp[i] ;
227 }
228
229 return;
230}
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: