Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HelixMixedStepper Class Reference

#include <G4HelixMixedStepper.hh>

Inheritance diagram for G4HelixMixedStepper:
Collaboration diagram for G4HelixMixedStepper:

Public Member Functions

 G4HelixMixedStepper (G4Mag_EqRhs *EqRhs, G4int StepperNumber=-1, G4double Angle_threshold=-1.0)
 
 ~G4HelixMixedStepper ()
 
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
 
void DumbStepper (const G4double y[], G4ThreeVector Bfld, G4double h, G4double yout[])
 
G4double DistChord () const
 
void SetVerbose (G4int newvalue)
 
void PrintCalls ()
 
G4MagIntegratorStepperSetupStepper (G4Mag_EqRhs *EqRhs, G4int StepperName)
 
void SetAngleThreshold (G4double val)
 
G4double GetAngleThreshold ()
 
G4int IntegratorOrder () const
 
- Public Member Functions inherited from G4MagHelicalStepper
 G4MagHelicalStepper (G4Mag_EqRhs *EqRhs)
 
virtual ~G4MagHelicalStepper ()
 
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
 
G4int IntegrationOrder ()
 
G4EquationOfMotionGetEquationOfMotion ()
 
void SetEquationOfMotion (G4EquationOfMotion *newEquation)
 
unsigned long GetfNoRHSCalls ()
 
void ResetfNORHSCalls ()
 
bool IsFSAL ()
 

Additional Inherited Members

- Protected Member Functions inherited from G4MagHelicalStepper
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 66 of file G4HelixMixedStepper.hh.

Constructor & Destructor Documentation

G4HelixMixedStepper::G4HelixMixedStepper ( G4Mag_EqRhs EqRhs,
G4int  StepperNumber = -1,
G4double  Angle_threshold = -1.0 
)

Definition at line 65 of file G4HelixMixedStepper.cc.

67  : G4MagHelicalStepper(EqRhs), fNumCallsRK4(0), fNumCallsHelix(0)
68 {
69  SetVerbose(1);
70  if( angleThreshold < 0.0 ){
71  fAngle_threshold= 0.33*pi;
72  }else{
73  fAngle_threshold= angleThreshold;
74  }
75 
76  if(stepperNumber<0)
77  stepperNumber=4; // Default is RK4 (original)
78  // stepperNumber=745; // Default is DormandPrince745 (ie DoPri5)
79  // stepperNumber=8; // Default is CashKarp
80 
81  fStepperNumber = stepperNumber; // Store the choice
82  fRK4Stepper = SetupStepper(EqRhs, fStepperNumber);
83 }
G4MagIntegratorStepper * SetupStepper(G4Mag_EqRhs *EqRhs, G4int StepperName)
G4MagHelicalStepper(G4Mag_EqRhs *EqRhs)
void SetVerbose(G4int newvalue)
static constexpr double pi
Definition: G4SIunits.hh:75

Here is the call graph for this function:

G4HelixMixedStepper::~G4HelixMixedStepper ( )

Definition at line 85 of file G4HelixMixedStepper.cc.

86 {
87  delete(fRK4Stepper);
88  if (fVerbose>0){ PrintCalls();};
89 }

Here is the call graph for this function:

Member Function Documentation

G4double G4HelixMixedStepper::DistChord ( ) const
virtual

Implements G4MagIntegratorStepper.

Definition at line 160 of file G4HelixMixedStepper.cc.

161 {
162  // Implementation : must check whether h/R > 2 pi !!
163  // If( h/R < pi) use G4LineSection::DistLine
164  // Else DistChord=R_helix
165  //
166  G4double distChord;
167  G4double Ang_curve=GetAngCurve();
168 
169  if(Ang_curve<=pi){
170  distChord=GetRadHelix()*(1-std::cos(0.5*Ang_curve));
171  }
172  else
173  {
174  if(Ang_curve<twopi){
175  distChord=GetRadHelix()*(1+std::cos(0.5*(twopi-Ang_curve)));
176  }
177  else{
178  distChord=2.*GetRadHelix();
179  }
180  }
181 
182  return distChord;
183 }
G4double GetRadHelix() const
static constexpr double twopi
Definition: G4SIunits.hh:76
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
G4double GetAngCurve() const

Here is the call graph for this function:

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

Implements G4MagHelicalStepper.

Definition at line 152 of file G4HelixMixedStepper.cc.

156 {
157  AdvanceHelix(yIn, Bfld, h, yOut);
158 }
void AdvanceHelix(const G4double yIn[], G4ThreeVector Bfld, G4double h, G4double yHelix[], G4double yHelix2[]=0)

Here is the call graph for this function:

G4double G4HelixMixedStepper::GetAngleThreshold ( )
inline

Definition at line 101 of file G4HelixMixedStepper.hh.

101 { return fAngle_threshold; }
G4int G4HelixMixedStepper::IntegratorOrder ( ) const
inlinevirtual

Implements G4MagIntegratorStepper.

Definition at line 103 of file G4HelixMixedStepper.hh.

103 { return 4; }
void G4HelixMixedStepper::PrintCalls ( )

Definition at line 186 of file G4HelixMixedStepper.cc.

187 {
188  G4cout << "In HelixMixedStepper::Number of calls to smallStepStepper = "
189  << fNumCallsRK4
190  << " and Number of calls to Helix = " << fNumCallsHelix << G4endl;
191 }
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61

Here is the caller graph for this function:

void G4HelixMixedStepper::SetAngleThreshold ( G4double  val)
inline

Definition at line 100 of file G4HelixMixedStepper.hh.

100 { fAngle_threshold= val;}
G4MagIntegratorStepper * G4HelixMixedStepper::SetupStepper ( G4Mag_EqRhs EqRhs,
G4int  StepperName 
)

Definition at line 194 of file G4HelixMixedStepper.cc.

195 {
196  G4MagIntegratorStepper* pStepper;
197  if (fVerbose>0) G4cout << " G4HelixMixedStepper: ";
198  switch ( StepperNumber )
199  {
200  // Robust, classic method
201  case 4:
202  pStepper = new G4ClassicalRK4( pE );
203  if (fVerbose>0) G4cout << "G4ClassicalRK4";
204  break;
205 
206  // Steppers with embedded estimation of error
207  case 8:
208  pStepper = new G4CashKarpRKF45( pE );
209  if (fVerbose>0) G4cout << "G4CashKarpRKF45";
210  break;
211  case 13:
212  pStepper = new G4NystromRK4( pE );
213  if (fVerbose>0) G4cout << "G4NystromRK4";
214  break;
215 
216  // Lowest order RK Stepper - experimental
217  case 1:
218  pStepper = new G4ImplicitEuler( pE );
219  if (fVerbose>0) G4cout << "G4ImplicitEuler";
220  break;
221 
222  // Lower order RK Steppers - ok overall, good for uneven fields
223  case 2:
224  pStepper = new G4SimpleRunge( pE );
225  if (fVerbose>0) G4cout << "G4SimpleRunge";
226  break;
227  case 3:
228  pStepper = new G4SimpleHeum( pE );
229  if (fVerbose>0) G4cout << "G4SimpleHeum";
230  break;
231  case 23:
232  pStepper = new G4BogackiShampine23( pE );
233  if (fVerbose>0) G4cout << "G4BogackiShampine23";
234  break;
235 
236  // Higher order RK Steppers
237  // for smoother fields and high accuracy requirements
238  case 45:
239  pStepper = new G4BogackiShampine45( pE );
240  if (fVerbose>0) G4cout << "G4BogackiShampine45";
241  break;
242  case 145:
243  pStepper = new G4TsitourasRK45( pE );
244  if (fVerbose>0) G4cout << "G4TsitourasRK45";
245  break;
246  case 745:
247  pStepper = new G4DormandPrince745( pE );
248  if (fVerbose>0) G4cout << "G4DormandPrince745";
249  break;
250 
251  // Helical Steppers
252  case 6:
253  pStepper = new G4HelixImplicitEuler( pE );
254  if (fVerbose>0) G4cout << "G4HelixImplicitEuler";
255  break;
256  case 7:
257  pStepper = new G4HelixSimpleRunge( pE );
258  if (fVerbose>0) G4cout << "G4HelixSimpleRunge";
259  break;
260  case 5:
261  pStepper = new G4HelixExplicitEuler( pE );
262  if (fVerbose>0) G4cout << "G4HelixExplicitEuler";
263  break; // Since Helix Explicit is used for long steps,
264  // this is useful only to measure overhead.
265  // Exact Helix - likely good only for cases of
266  // i) uniform field (potentially over small distances)
267  // ii) segmented uniform field (maybe)
268  case 9:
269  pStepper = new G4ExactHelixStepper( pE );
270  if (fVerbose>0) G4cout << "G4ExactHelixStepper";
271  break;
272  case 10:
273  pStepper = new G4RKG3_Stepper( pE );
274  if (fVerbose>0) G4cout << "G4RKG3_Stepper";
275  break;
276 
277  // Low Order Steppers - not good except for very weak fields
278  case 11:
279  pStepper = new G4ExplicitEuler( pE );
280  if (fVerbose>0) G4cout << "G4ExplicitEuler";
281  break;
282  case 12:
283  pStepper = new G4ImplicitEuler( pE );
284  if (fVerbose>0) G4cout << "G4ImplicitEuler";
285  break;
286 
287  case 0:
288  case -1:
289  default:
290  pStepper = new G4ClassicalRK4( pE );
291  if (fVerbose>0) G4cout << "G4ClassicalRK4 (Default)";
292  break;
293  }
294  if(fVerbose>0)
295  G4cout << " chosen as stepper for small steps in G4HelixMixedStepper."
296  << G4endl;
297 
298  return pStepper;
299 }
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61

Here is the caller graph for this function:

void G4HelixMixedStepper::SetVerbose ( G4int  newvalue)
inline

Definition at line 94 of file G4HelixMixedStepper.hh.

94 {fVerbose=newvalue;}

Here is the caller graph for this function:

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

Reimplemented from G4MagHelicalStepper.

Definition at line 91 of file G4HelixMixedStepper.cc.

96 {
97  //Estimation of the Stepping Angle
98  G4ThreeVector Bfld;
99  MagFieldEvaluate(yInput, Bfld);
100 
101  G4double Bmag = Bfld.mag();
102  const G4double *pIn = yInput+3;
103  G4ThreeVector initVelocity= G4ThreeVector( pIn[0], pIn[1], pIn[2]);
104  G4double velocityVal = initVelocity.mag();
105  G4double R_1;
106  G4double Ang_curve;
107 
108  R_1=std::abs(GetInverseCurve(velocityVal,Bmag));
109  Ang_curve=R_1*Step;
110  SetAngCurve(Ang_curve);
111  SetCurve(std::abs(1/R_1));
112 
113  if(Ang_curve< fAngle_threshold){
114  fNumCallsRK4++;
115  fRK4Stepper->Stepper(yInput,dydx,Step,yOut,yErr);
116  }
117  else
118  {
119  fNumCallsHelix++;
120  const G4int nvar = 6 ;
121  const G4int nvarMax = 8 ;
122  G4int i;
123  G4double yTemp[nvarMax], yIn[nvarMax], yTemp2[nvarMax];
124  G4ThreeVector Bfld_midpoint;
125 
126  // Saving yInput because yInput and yOut can be aliases for same array
127  for(i=0;i<nvar;i++) yIn[i]=yInput[i];
128 
129  G4double halfS = Step * 0.5;
130  // 1. Do first half step and full step
131  AdvanceHelix(yIn, Bfld, halfS, yTemp, yTemp2); // yTemp2 for s=2*h (halfS)
132  //**********
133 
134  MagFieldEvaluate(yTemp, Bfld_midpoint) ;
135 
136  // 2. Do second half step - with revised field
137  // NOTE: Could avoid this call if 'Bfld_midpoint == Bfld'
138  // or diff 'almost' zero
139  AdvanceHelix(yTemp, Bfld_midpoint, halfS, yOut);
140  // Not requesting y at s=2*h (halfS)
141  //**********
142 
143  // 3. Estimate the integration error
144  // should be (nearly) zero if Bfield= constant
145  for(i=0;i<nvar;i++) {
146  yErr[i] = yOut[i] - yTemp2[i] ;
147  }
148  }
149 }
void AdvanceHelix(const G4double yIn[], G4ThreeVector Bfld, G4double h, G4double yHelix[], G4double yHelix2[]=0)
CLHEP::Hep3Vector G4ThreeVector
virtual void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])=0
void MagFieldEvaluate(const G4double y[], G4ThreeVector &Bfield)
int G4int
Definition: G4Types.hh:78
G4double GetInverseCurve(const G4double Momentum, const G4double Bmag)
void SetAngCurve(const G4double Ang)
void SetCurve(const G4double Curve)
Definition: Step.hh:41
double G4double
Definition: G4Types.hh:76
double mag() const

Here is the call graph for this function:


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