Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4MagInt_Driver Class Reference

#include <G4MagIntegratorDriver.hh>

Public Member Functions

G4bool AccurateAdvance (G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
 
G4bool QuickAdvance (G4FieldTrack &y_val, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr)
 
G4bool QuickAdvance (G4FieldTrack &y_posvel, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr_pos_sq, G4double &dyerr_mom_rel_sq)
 
 G4MagInt_Driver (G4double hminimum, G4MagIntegratorStepper *pItsStepper, G4int numberOfComponents=6, G4int statisticsVerbosity=1)
 
 ~G4MagInt_Driver ()
 
G4double GetHmin () const
 
G4double Hmin () const
 
G4double GetSafety () const
 
G4double GetPshrnk () const
 
G4double GetPgrow () const
 
G4double GetErrcon () const
 
void GetDerivatives (const G4FieldTrack &y_curr, G4double dydx[])
 
void RenewStepperAndAdjust (G4MagIntegratorStepper *pItsStepper)
 
void ReSetParameters (G4double new_safety=0.9)
 
void SetSafety (G4double valS)
 
void SetPshrnk (G4double valPs)
 
void SetPgrow (G4double valPg)
 
void SetErrcon (G4double valEc)
 
G4double ComputeAndSetErrcon ()
 
const G4MagIntegratorStepperGetStepper () const
 
G4MagIntegratorStepperGetStepper ()
 
void OneGoodStep (G4double ystart[], const G4double dydx[], G4double &x, G4double htry, G4double eps, G4double &hdid, G4double &hnext)
 
G4double ComputeNewStepSize (G4double errMaxNorm, G4double hstepCurrent)
 
G4double ComputeNewStepSize_WithinLimits (G4double errMaxNorm, G4double hstepCurrent)
 
G4int GetMaxNoSteps () const
 
void SetMaxNoSteps (G4int val)
 
void SetHmin (G4double newval)
 
void SetVerboseLevel (G4int newLevel)
 
G4double GetVerboseLevel () const
 
G4double GetSmallestFraction () const
 
void SetSmallestFraction (G4double val)
 

Protected Member Functions

void WarnSmallStepSize (G4double hnext, G4double hstep, G4double h, G4double xDone, G4int noSteps)
 
void WarnTooManyStep (G4double x1start, G4double x2end, G4double xCurrent)
 
void WarnEndPointTooFar (G4double endPointDist, G4double hStepSize, G4double epsilonRelative, G4int debugFlag)
 
void PrintStatus (const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)
 
void PrintStatus (const G4FieldTrack &StartFT, const G4FieldTrack &CurrentFT, G4double requestStep, G4int subStepNo)
 
void PrintStat_Aux (const G4FieldTrack &aFieldTrack, G4double requestStep, G4double actualStep, G4int subStepNo, G4double subStepSize, G4double dotVelocities)
 
void PrintStatisticsReport ()
 

Detailed Description

Definition at line 48 of file G4MagIntegratorDriver.hh.

Constructor & Destructor Documentation

G4MagInt_Driver::G4MagInt_Driver ( G4double  hminimum,
G4MagIntegratorStepper pItsStepper,
G4int  numberOfComponents = 6,
G4int  statisticsVerbosity = 1 
)

Definition at line 69 of file G4MagIntegratorDriver.cc.

73  : fSmallestFraction( 1.0e-12 ),
74  fNoIntegrationVariables(numComponents),
75  fMinNoVars(12),
76  fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )),
77  fStatisticsVerboseLevel(statisticsVerbose),
78  fNoTotalSteps(0), fNoBadSteps(0), fNoSmallSteps(0),
79  fNoInitialSmallSteps(0),
80  fDyerr_max(0.0), fDyerr_mx2(0.0),
81  fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0),
82  fSumH_sm(0.0), fSumH_lg(0.0),
83  fVerboseLevel(0)
84 {
85  // In order to accomodate "Laboratory Time", which is [7], fMinNoVars=8
86  // is required. For proper time of flight and spin, fMinNoVars must be 12
87 
88  RenewStepperAndAdjust( pStepper );
89  fMinimumStep= hminimum;
90  fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
91 #ifdef G4DEBUG_FIELD
92  fVerboseLevel=2;
93 #endif
94 
95  if( (fVerboseLevel > 0) || (fStatisticsVerboseLevel > 1) )
96  {
97  G4cout << "MagIntDriver version: Accur-Adv: "
98  << "invE_nS, QuickAdv-2sqrt with Statistics "
99 #ifdef G4FLD_STATS
100  << " enabled "
101 #else
102  << " disabled "
103 #endif
104  << G4endl;
105  }
106 }
void RenewStepperAndAdjust(G4MagIntegratorStepper *pItsStepper)
virtual G4int IntegratorOrder() const =0
G4GLOB_DLL std::ostream G4cout
T max(const T t1, const T t2)
brief Return the largest of the two arguments
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

G4MagInt_Driver::~G4MagInt_Driver ( )

Definition at line 112 of file G4MagIntegratorDriver.cc.

113 {
114  if( fStatisticsVerboseLevel > 1 )
115  {
117  }
118 }

Here is the call graph for this function:

Member Function Documentation

G4bool G4MagInt_Driver::AccurateAdvance ( G4FieldTrack y_current,
G4double  hstep,
G4double  eps,
G4double  hinitial = 0.0 
)

Definition at line 127 of file G4MagIntegratorDriver.cc.

131 {
132  // Runge-Kutta driver with adaptive stepsize control. Integrate starting
133  // values at y_current over hstep x2 with accuracy eps.
134  // On output ystart is replaced by values at the end of the integration
135  // interval. RightHandSide is the right-hand side of ODE system.
136  // The source is similar to odeint routine from NRC p.721-722 .
137 
138  G4int nstp, i, no_warnings=0;
139  G4double x, hnext, hdid, h;
140 
141 #ifdef G4DEBUG_FIELD
142  static G4int dbg=1;
143  static G4int nStpPr=50; // For debug printing of long integrations
144  G4double ySubStepStart[G4FieldTrack::ncompSVEC];
145  G4FieldTrack yFldTrkStart(y_current);
146 #endif
147 
150  G4double x1, x2;
151  G4bool succeeded = true, lastStepSucceeded;
152 
153  G4double startCurveLength;
154 
155  G4int noFullIntegr=0, noSmallIntegr = 0 ;
156  static G4ThreadLocal G4int noGoodSteps =0 ; // Bad = chord > curve-len
157  const G4int nvar= fNoVars;
158 
159  G4FieldTrack yStartFT(y_current);
160 
161  // Ensure that hstep > 0
162  //
163  if( hstep <= 0.0 )
164  {
165  if(hstep==0.0)
166  {
167  std::ostringstream message;
168  message << "Proposed step is zero; hstep = " << hstep << " !";
169  G4Exception("G4MagInt_Driver::AccurateAdvance()",
170  "GeomField1001", JustWarning, message);
171  return succeeded;
172  }
173  else
174  {
175  std::ostringstream message;
176  message << "Invalid run condition." << G4endl
177  << "Proposed step is negative; hstep = " << hstep << "." << G4endl
178  << "Requested step cannot be negative! Aborting event.";
179  G4Exception("G4MagInt_Driver::AccurateAdvance()",
180  "GeomField0003", EventMustBeAborted, message);
181  return false;
182  }
183  }
184 
185  y_current.DumpToArray( ystart );
186 
187  startCurveLength= y_current.GetCurveLength();
188  x1= startCurveLength;
189  x2= x1 + hstep;
190 
191  if ( (hinitial > 0.0) && (hinitial < hstep)
192  && (hinitial > perMillion * hstep) )
193  {
194  h = hinitial;
195  }
196  else // Initial Step size "h" defaults to the full interval
197  {
198  h = hstep;
199  }
200 
201  x = x1;
202 
203  for (i=0;i<nvar;i++) { y[i] = ystart[i]; }
204 
205  G4bool lastStep= false;
206  nstp=1;
207 
208  do
209  {
210  G4ThreeVector StartPos( y[0], y[1], y[2] );
211 
212 #ifdef G4DEBUG_FIELD
213  G4double xSubStepStart= x;
214  for (i=0;i<nvar;i++) { ySubStepStart[i] = y[i]; }
215  yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
216  yFldTrkStart.SetCurveLength(x);
217 #endif
218 
219  // Old method - inline call to Equation of Motion
220  // pIntStepper->RightHandSide( y, dydx );
221  // New method allows to cache field, or state (eg momentum magnitude)
222  pIntStepper->ComputeRightHandSide( y, dydx );
223  fNoTotalSteps++;
224 
225  // Perform the Integration
226  //
227  if( h > fMinimumStep )
228  {
229  OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ;
230  //--------------------------------------
231  lastStepSucceeded= (hdid == h);
232 #ifdef G4DEBUG_FIELD
233  if (dbg>2) {
234  PrintStatus( ySubStepStart, xSubStepStart, y, x, h, nstp); // Only
235  }
236 #endif
237  }
238  else
239  {
240  G4FieldTrack yFldTrk( G4ThreeVector(0,0,0),
241  G4ThreeVector(0,0,0), 0., 0., 0., 0. );
242  G4double dchord_step, dyerr, dyerr_len; // What to do with these ?
243  yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
244  yFldTrk.SetCurveLength( x );
245 
246  QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
247  //-----------------------------------------------------
248 
249  yFldTrk.DumpToArray(y);
250 
251 #ifdef G4FLD_STATS
252  fNoSmallSteps++;
253  if ( dyerr_len > fDyerr_max) { fDyerr_max= dyerr_len; }
254  fDyerrPos_smTot += dyerr_len;
255  fSumH_sm += h; // Length total for 'small' steps
256  if (nstp<=1) { fNoInitialSmallSteps++; }
257 #endif
258 #ifdef G4DEBUG_FIELD
259  if (dbg>1)
260  {
261  if(fNoSmallSteps<2) { PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
262  G4cout << "Another sub-min step, no " << fNoSmallSteps
263  << " of " << fNoTotalSteps << " this time " << nstp << G4endl;
264  PrintStatus( ySubStepStart, x1, y, x, h, nstp); // Only this
265  G4cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h
266  << " epsilon= " << eps << " hstep= " << hstep
267  << " h= " << h << " hmin= " << fMinimumStep << G4endl;
268  }
269 #endif
270  if( h == 0.0 )
271  {
272  G4Exception("G4MagInt_Driver::AccurateAdvance()",
273  "GeomField0003", FatalException,
274  "Integration Step became Zero!");
275  }
276  dyerr = dyerr_len / h;
277  hdid= h;
278  x += hdid;
279 
280  // Compute suggested new step
281  hnext= ComputeNewStepSize( dyerr/eps, h);
282 
283  // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
284  lastStepSucceeded= (dyerr<= eps);
285  }
286 
287  if (lastStepSucceeded) { noFullIntegr++; }
288  else { noSmallIntegr++; }
289 
290  G4ThreeVector EndPos( y[0], y[1], y[2] );
291 
292 #ifdef G4DEBUG_FIELD
293  if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
294  {
295  if( nstp==nStpPr ) { G4cout << "***** Many steps ****" << G4endl; }
296  G4cout << "MagIntDrv: " ;
297  G4cout << "hdid=" << std::setw(12) << hdid << " "
298  << "hnext=" << std::setw(12) << hnext << " "
299  << "hstep=" << std::setw(12) << hstep << " (requested) "
300  << G4endl;
301  PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
302  }
303 #endif
304 
305  // Check the endpoint
306  G4double endPointDist= (EndPos-StartPos).mag();
307  if ( endPointDist >= hdid*(1.+perMillion) )
308  {
309  fNoBadSteps++;
310 
311  // Issue a warning only for gross differences -
312  // we understand how small difference occur.
313  if ( endPointDist >= hdid*(1.+perThousand) )
314  {
315 #ifdef G4DEBUG_FIELD
316  if (dbg)
317  {
318  WarnEndPointTooFar ( endPointDist, hdid, eps, dbg );
319  G4cerr << " Total steps: bad " << fNoBadSteps
320  << " good " << noGoodSteps << " current h= " << hdid
321  << G4endl;
322  PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
323  }
324 #endif
325  no_warnings++;
326  }
327  }
328  else
329  {
330  noGoodSteps ++;
331  }
332 // #endif
333 
334  // Avoid numerous small last steps
335  if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
336  {
337  // No more integration -- the next step will not happen
338  lastStep = true;
339  }
340  else
341  {
342  // Check the proposed next stepsize
343  if(std::fabs(hnext) <= Hmin())
344  {
345 #ifdef G4DEBUG_FIELD
346  // If simply a very small interval is being integrated, do not warn
347  if( (x < x2 * (1-eps) ) && // The last step can be small: OK
348  (std::fabs(hstep) > Hmin()) ) // and if we are asked, it's OK
349  {
350  if(dbg>0)
351  {
352  WarnSmallStepSize( hnext, hstep, h, x-x1, nstp );
353  PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
354  }
355  no_warnings++;
356  }
357 #endif
358  // Make sure that the next step is at least Hmin.
359  h = Hmin();
360  }
361  else
362  {
363  h = hnext;
364  }
365 
366  // Ensure that the next step does not overshoot
367  if ( x+h > x2 )
368  { // When stepsize overshoots, decrease it!
369  h = x2 - x ; // Must cope with difficult rounding-error
370  } // issues if hstep << x2
371 
372  if ( h == 0.0 )
373  {
374  // Cannot progress - accept this as last step - by default
375  lastStep = true;
376 #ifdef G4DEBUG_FIELD
377  if (dbg>2)
378  {
379  int prec= G4cout.precision(12);
380  G4cout << "Warning: G4MagIntegratorDriver::AccurateAdvance"
381  << G4endl
382  << " Integration step 'h' became "
383  << h << " due to roundoff. " << G4endl
384  << " Calculated as difference of x2= "<< x2 << " and x=" << x
385  << " Forcing termination of advance." << G4endl;
386  G4cout.precision(prec);
387  }
388 #endif
389  }
390  }
391  } while ( ((nstp++)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
392  // Loop checking, 07.10.2016, J. Apostolakis
393 
394  // Have we reached the end ?
395  // --> a better test might be x-x2 > an_epsilon
396 
397  succeeded= (x>=x2); // If it was a "forced" last step
398 
399  for (i=0;i<nvar;i++) { yEnd[i] = y[i]; }
400 
401  // Put back the values.
402  y_current.LoadFromArray( yEnd, fNoIntegrationVariables );
403  y_current.SetCurveLength( x );
404 
405  if(nstp > fMaxNoSteps)
406  {
407  no_warnings++;
408  succeeded = false;
409 #ifdef G4DEBUG_FIELD
410  if (dbg)
411  {
412  WarnTooManyStep( x1, x2, x ); // Issue WARNING
413  PrintStatus( yEnd, x1, y, x, hstep, -nstp);
414  }
415 #endif
416  }
417 
418 #ifdef G4DEBUG_FIELD
419  if( dbg && no_warnings )
420  {
421  G4cerr << "G4MagIntegratorDriver exit status: no-steps " << nstp <<G4endl;
422  PrintStatus( yEnd, x1, y, x, hstep, nstp);
423  }
424 #endif
425 
426  return succeeded;
427 } // end of AccurateAdvance ...........................
static constexpr double perMillion
Definition: G4SIunits.hh:334
virtual void ComputeRightHandSide(const G4double y[], G4double dydx[])
G4double GetCurveLength() const
CLHEP::Hep3Vector G4ThreeVector
void SetCurveLength(G4double nCurve_s)
void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent)
void WarnEndPointTooFar(G4double endPointDist, G4double hStepSize, G4double epsilonRelative, G4int debugFlag)
static const G4double eps
#define G4ThreadLocal
Definition: tls.hh:89
G4double Hmin() const
int G4int
Definition: G4Types.hh:78
void OneGoodStep(G4double ystart[], const G4double dydx[], G4double &x, G4double htry, G4double eps, G4double &hdid, G4double &hnext)
static const double prec
Definition: RanecuEngine.cc:58
G4GLOB_DLL std::ostream G4cout
void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)
bool G4bool
Definition: G4Types.hh:79
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4bool QuickAdvance(G4FieldTrack &y_val, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr)
void DumpToArray(G4double valArr[ncompSVEC]) const
void WarnSmallStepSize(G4double hnext, G4double hstep, G4double h, G4double xDone, G4int noSteps)
void LoadFromArray(const G4double valArr[ncompSVEC], G4int noVarsIntegrated)
G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent)
#define G4endl
Definition: G4ios.hh:61
static constexpr double perThousand
Definition: G4SIunits.hh:333
double G4double
Definition: G4Types.hh:76
G4GLOB_DLL std::ostream G4cerr

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4MagInt_Driver::ComputeAndSetErrcon ( )
inline
G4double G4MagInt_Driver::ComputeNewStepSize ( G4double  errMaxNorm,
G4double  hstepCurrent 
)

Definition at line 764 of file G4MagIntegratorDriver.cc.

767 {
768  G4double hnew;
769 
770  // Compute size of next Step for a failed step
771  if(errMaxNorm > 1.0 )
772  {
773  // Step failed; compute the size of retrial Step.
774  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
775  } else if(errMaxNorm > 0.0 ) {
776  // Compute size of next Step for a successful step
777  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
778  } else {
779  // if error estimate is zero (possible) or negative (dubious)
780  hnew = max_stepping_increase * hstepCurrent;
781  }
782 
783  return hnew;
784 }
G4double GetSafety() const
G4double GetPgrow() const
G4double GetPshrnk() const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4MagInt_Driver::ComputeNewStepSize_WithinLimits ( G4double  errMaxNorm,
G4double  hstepCurrent 
)

Definition at line 794 of file G4MagIntegratorDriver.cc.

797 {
798  G4double hnew;
799 
800  // Compute size of next Step for a failed step
801  if (errMaxNorm > 1.0 )
802  {
803  // Step failed; compute the size of retrial Step.
804  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
805 
806  if (hnew < max_stepping_decrease*hstepCurrent)
807  {
808  hnew = max_stepping_decrease*hstepCurrent ;
809  // reduce stepsize, but no more
810  // than this factor (value= 1/10)
811  }
812  }
813  else
814  {
815  // Compute size of next Step for a successful step
816  if (errMaxNorm > errcon)
817  { hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()); }
818  else // No more than a factor of 5 increase
819  { hnew = max_stepping_increase * hstepCurrent; }
820  }
821  return hnew;
822 }
G4double GetSafety() const
G4double GetPgrow() const
G4double GetPshrnk() const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4MagInt_Driver::GetDerivatives ( const G4FieldTrack y_curr,
G4double  dydx[] 
)
inline
G4double G4MagInt_Driver::GetErrcon ( ) const
inline
G4double G4MagInt_Driver::GetHmin ( ) const
inline
G4int G4MagInt_Driver::GetMaxNoSteps ( ) const
inline
G4double G4MagInt_Driver::GetPgrow ( ) const
inline

Here is the caller graph for this function:

G4double G4MagInt_Driver::GetPshrnk ( ) const
inline

Here is the caller graph for this function:

G4double G4MagInt_Driver::GetSafety ( ) const
inline

Here is the caller graph for this function:

G4double G4MagInt_Driver::GetSmallestFraction ( ) const
inline
const G4MagIntegratorStepper* G4MagInt_Driver::GetStepper ( ) const
inline

Here is the caller graph for this function:

G4MagIntegratorStepper* G4MagInt_Driver::GetStepper ( )
inline
G4double G4MagInt_Driver::GetVerboseLevel ( ) const
inline
G4double G4MagInt_Driver::Hmin ( ) const
inline

Here is the caller graph for this function:

void G4MagInt_Driver::OneGoodStep ( G4double  ystart[],
const G4double  dydx[],
G4double x,
G4double  htry,
G4double  eps,
G4double hdid,
G4double hnext 
)

Definition at line 517 of file G4MagIntegratorDriver.cc.

538 {
539  G4double errmax_sq;
540  G4double h, htemp, xnew ;
541 
543 
544  h = htry ; // Set stepsize to the initial trial value
545 
546  G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
547 
548  G4double errpos_sq=0.0; // square of displacement error
549  G4double errvel_sq=0.0; // square of momentum vector difference
550  G4double errspin_sq=0.0; // square of spin vector difference
551 
552  G4int iter;
553 
554  static G4ThreadLocal G4int tot_no_trials=0;
555  const G4int max_trials=100;
556 
557  G4ThreeVector Spin(y[9],y[10],y[11]);
558  G4double spin_mag2 =Spin.mag2() ;
559  G4bool hasSpin= (spin_mag2 > 0.0);
560 
561  for (iter=0; iter<max_trials ;iter++)
562  {
563  tot_no_trials++;
564  pIntStepper-> Stepper(y,dydx,h,ytemp,yerr);
565  // *******
566  G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep);
567  G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos);
568 
569  // Evaluate accuracy
570  //
571  errpos_sq = sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
572  errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
573 
574  // Accuracy for momentum
575  G4double magvel_sq= sqr(y[3]) + sqr(y[4]) + sqr(y[5]) ;
576  G4double sumerr_sq = sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) ;
577  if( magvel_sq > 0.0 ) {
578  errvel_sq = sumerr_sq / magvel_sq;
579  }else{
580  G4cerr << "** G4MagIntegrationDriver: found case of zero momentum."
581  << " iteration= " << iter << " h= " << h << G4endl;
582  errvel_sq = sumerr_sq;
583  }
584  errvel_sq *= inv_eps_vel_sq;
585  errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
586 
587  if( hasSpin )
588  {
589  // Accuracy for spin
590  errspin_sq = ( sqr(yerr[9]) + sqr(yerr[10]) + sqr(yerr[11]) )
591  / spin_mag2; // ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
592  errspin_sq *= inv_eps_vel_sq;
593  errmax_sq = std::max( errmax_sq, errspin_sq );
594  }
595 
596  if ( errmax_sq <= 1.0 ) { break; } // Step succeeded.
597 
598  // Step failed; compute the size of retrial Step.
599  htemp = GetSafety()*h* std::pow( errmax_sq, 0.5*GetPshrnk() );
600 
601  if (htemp >= 0.1*h) { h = htemp; } // Truncation error too large,
602  else { h = 0.1*h; } // reduce stepsize, but no more
603  // than a factor of 10
604  xnew = x + h;
605  if(xnew == x)
606  {
607  G4cerr << "G4MagIntegratorDriver::OneGoodStep:" << G4endl
608  << " Stepsize underflow in Stepper " << G4endl ;
609  G4cerr << " Step's start x=" << x << " and end x= " << xnew
610  << " are equal !! " << G4endl
611  <<" Due to step-size= " << h
612  << " . Note that input step was " << htry << G4endl;
613  break;
614  }
615  }
616 
617 #ifdef G4FLD_STATS
618  // Sum of squares of position error // and momentum dir (underestimated)
619  fSumH_lg += h;
620  fDyerrPos_lgTot += errpos_sq;
621  fDyerrVel_lgTot += errvel_sq * h * h;
622 #endif
623 
624  // Compute size of next Step
625  if (errmax_sq > errcon*errcon)
626  {
627  hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow());
628  }
629  else
630  {
631  hnext = max_stepping_increase*h ; // No more than a factor of 5 increase
632  }
633  x += (hdid = h);
634 
635  for(G4int k=0;k<fNoIntegrationVariables;k++) { y[k] = ytemp[k]; }
636 
637  return;
638 } // end of OneGoodStep .............................
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
G4double GetSafety() const
bool G4bool
Definition: G4Types.hh:79
G4double GetPgrow() const
G4double GetPshrnk() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
#define G4endl
Definition: G4ios.hh:61
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
G4GLOB_DLL std::ostream G4cerr

Here is the call graph for this function:

Here is the caller graph for this function:

void G4MagInt_Driver::PrintStat_Aux ( const G4FieldTrack aFieldTrack,
G4double  requestStep,
G4double  actualStep,
G4int  subStepNo,
G4double  subStepSize,
G4double  dotVelocities 
)
protected

Definition at line 929 of file G4MagIntegratorDriver.cc.

936 {
937  const G4ThreeVector Position= aFieldTrack.GetPosition();
938  const G4ThreeVector UnitVelocity= aFieldTrack.GetMomentumDir();
939 
940  if( subStepNo >= 0)
941  {
942  G4cout << std::setw( 5) << subStepNo << " ";
943  }
944  else
945  {
946  G4cout << std::setw( 5) << "Start" << " ";
947  }
948  G4double curveLen= aFieldTrack.GetCurveLength();
949  G4cout << std::setw( 7) << curveLen;
950  G4cout << std::setw( 9) << Position.x() << " "
951  << std::setw( 9) << Position.y() << " "
952  << std::setw( 9) << Position.z() << " "
953  << std::setw( 8) << UnitVelocity.x() << " "
954  << std::setw( 8) << UnitVelocity.y() << " "
955  << std::setw( 8) << UnitVelocity.z() << " ";
956  G4int oldprec= G4cout.precision(3);
957  G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 << " ";
958  G4cout.precision(6);
959  G4cout << std::setw(10) << dotVeloc_StartCurr << " ";
960  G4cout.precision(oldprec);
961  G4cout << std::setw( 7) << aFieldTrack.GetKineticEnergy();
962  G4cout << std::setw(12) << step_len << " ";
963 
964  static G4ThreadLocal G4double oldCurveLength= 0.0;
965  static G4ThreadLocal G4double oldSubStepLength= 0.0;
966  static G4ThreadLocal G4int oldSubStepNo= -1;
967 
968  G4double subStep_len=0.0;
969  if( curveLen > oldCurveLength )
970  {
971  subStep_len= curveLen - oldCurveLength;
972  }
973  else if (subStepNo == oldSubStepNo)
974  {
975  subStep_len= oldSubStepLength;
976  }
977  oldCurveLength= curveLen;
978  oldSubStepLength= subStep_len;
979 
980  G4cout << std::setw(12) << subStep_len << " ";
981  G4cout << std::setw(12) << subStepSize << " ";
982  if( requestStep != -1.0 )
983  {
984  G4cout << std::setw( 9) << requestStep << " ";
985  }
986  else
987  {
988  G4cout << std::setw( 9) << " InitialStep " << " ";
989  }
990  G4cout << G4endl;
991 }
G4double GetCurveLength() const
double x() const
G4double GetKineticEnergy() const
const G4ThreeVector & GetMomentumDir() const
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
double z() const
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
double y() const
double mag2() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4MagInt_Driver::PrintStatisticsReport ( )
protected

Definition at line 995 of file G4MagIntegratorDriver.cc.

996 {
997  G4int noPrecBig= 6;
998  G4int oldPrec= G4cout.precision(noPrecBig);
999 
1000  G4cout << "G4MagInt_Driver Statistics of steps undertaken. " << G4endl;
1001  G4cout << "G4MagInt_Driver: Number of Steps: "
1002  << " Total= " << fNoTotalSteps
1003  << " Bad= " << fNoBadSteps
1004  << " Small= " << fNoSmallSteps
1005  << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
1006  << G4endl;
1007 
1008 #ifdef G4FLD_STATS
1009  G4cout << "MID dyerr: "
1010  << " maximum= " << fDyerr_max
1011  << " Sum small= " << fDyerrPos_smTot
1012  << " std::sqrt(Sum large^2): pos= " << std::sqrt(fDyerrPos_lgTot)
1013  << " vel= " << std::sqrt( fDyerrVel_lgTot )
1014  << " Total h-distance: small= " << fSumH_sm
1015  << " large= " << fSumH_lg
1016  << G4endl;
1017 
1018 #if 0
1019  G4int noPrecSmall=4;
1020  // Single line precis of statistics ... optional
1021  G4cout.precision(noPrecSmall);
1022  G4cout << "MIDnums: " << fMinimumStep
1023  << " " << fNoTotalSteps
1024  << " " << fNoSmallSteps
1025  << " " << fNoSmallSteps-fNoInitialSmallSteps
1026  << " " << fNoBadSteps
1027  << " " << fDyerr_max
1028  << " " << fDyerr_mx2
1029  << " " << fDyerrPos_smTot
1030  << " " << fSumH_sm
1031  << " " << fDyerrPos_lgTot
1032  << " " << fDyerrVel_lgTot
1033  << " " << fSumH_lg
1034  << G4endl;
1035 #endif
1036 #endif
1037 
1038  G4cout.precision(oldPrec);
1039 }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61

Here is the caller graph for this function:

void G4MagInt_Driver::PrintStatus ( const G4double StartArr,
G4double  xstart,
const G4double CurrentArr,
G4double  xcurrent,
G4double  requestStep,
G4int  subStepNo 
)
protected

Definition at line 826 of file G4MagIntegratorDriver.cc.

836 {
837  G4FieldTrack StartFT(G4ThreeVector(0,0,0),
838  G4ThreeVector(0,0,0), 0., 0., 0., 0. );
839  G4FieldTrack CurrentFT (StartFT);
840 
841  StartFT.LoadFromArray( StartArr, fNoIntegrationVariables);
842  StartFT.SetCurveLength( xstart);
843  CurrentFT.LoadFromArray( CurrentArr, fNoIntegrationVariables);
844  CurrentFT.SetCurveLength( xcurrent );
845 
846  PrintStatus(StartFT, CurrentFT, requestStep, subStepNo );
847 }
CLHEP::Hep3Vector G4ThreeVector
void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)

Here is the call graph for this function:

Here is the caller graph for this function:

void G4MagInt_Driver::PrintStatus ( const G4FieldTrack StartFT,
const G4FieldTrack CurrentFT,
G4double  requestStep,
G4int  subStepNo 
)
protected

Definition at line 851 of file G4MagIntegratorDriver.cc.

856 {
857  G4int verboseLevel= fVerboseLevel;
858  static G4ThreadLocal G4int noPrecision= 5;
859  G4int oldPrec= G4cout.precision(noPrecision);
860  // G4cout.setf(ios_base::fixed,ios_base::floatfield);
861 
862  const G4ThreeVector StartPosition= StartFT.GetPosition();
863  const G4ThreeVector StartUnitVelocity= StartFT.GetMomentumDir();
864  const G4ThreeVector CurrentPosition= CurrentFT.GetPosition();
865  const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir();
866 
867  G4double DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
868 
869  G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
870  G4double subStepSize = step_len;
871 
872  if( (subStepNo <= 1) || (verboseLevel > 3) )
873  {
874  subStepNo = - subStepNo; // To allow printing banner
875 
876  G4cout << std::setw( 6) << " " << std::setw( 25)
877  << " G4MagInt_Driver: Current Position and Direction" << " "
878  << G4endl;
879  G4cout << std::setw( 5) << "Step#" << " "
880  << std::setw( 7) << "s-curve" << " "
881  << std::setw( 9) << "X(mm)" << " "
882  << std::setw( 9) << "Y(mm)" << " "
883  << std::setw( 9) << "Z(mm)" << " "
884  << std::setw( 8) << " N_x " << " "
885  << std::setw( 8) << " N_y " << " "
886  << std::setw( 8) << " N_z " << " "
887  << std::setw( 8) << " N^2-1 " << " "
888  << std::setw(10) << " N(0).N " << " "
889  << std::setw( 7) << "KinEner " << " "
890  << std::setw(12) << "Track-l" << " " // Add the Sub-step ??
891  << std::setw(12) << "Step-len" << " "
892  << std::setw(12) << "Step-len" << " "
893  << std::setw( 9) << "ReqStep" << " "
894  << G4endl;
895  }
896 
897  if( (subStepNo <= 0) )
898  {
899  PrintStat_Aux( StartFT, requestStep, 0.,
900  0, 0.0, 1.0);
901  //*************
902  }
903 
904  if( verboseLevel <= 3 )
905  {
906  G4cout.precision(noPrecision);
907  PrintStat_Aux( CurrentFT, requestStep, step_len,
908  subStepNo, subStepSize, DotStartCurrentVeloc );
909  //*************
910  }
911 
912  else // if( verboseLevel > 3 )
913  {
914  // Multi-line output
915 
916  // G4cout << "Current Position is " << CurrentPosition << G4endl
917  // << " and UnitVelocity is " << CurrentUnitVelocity << G4endl;
918  // G4cout << "Step taken was " << step_len
919  // << " out of PhysicalStep= " << requestStep << G4endl;
920  // G4cout << "Final safety is: " << safety << G4endl;
921  // G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
922  // << G4endl << G4endl;
923  }
924  G4cout.precision(oldPrec);
925 }
G4double GetCurveLength() const
double dot(const Hep3Vector &) const
const G4ThreeVector & GetMomentumDir() const
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
void PrintStat_Aux(const G4FieldTrack &aFieldTrack, G4double requestStep, G4double actualStep, G4int subStepNo, G4double subStepSize, G4double dotVelocities)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4bool G4MagInt_Driver::QuickAdvance ( G4FieldTrack y_val,
const G4double  dydx[],
G4double  hstep,
G4double dchord_step,
G4double dyerr 
)

Definition at line 663 of file G4MagIntegratorDriver.cc.

669 {
670  G4double dyerr_pos_sq, dyerr_mom_rel_sq;
673  G4double s_start;
674  G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
675 
676  static G4ThreadLocal G4int no_call=0;
677  no_call ++;
678 
679  // Move data into array
680  y_posvel.DumpToArray( yarrin ); // yarrin <== y_posvel
681  s_start = y_posvel.GetCurveLength();
682 
683  // Do an Integration Step
684  pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ;
685  // *******
686 
687  // Estimate curve-chord distance
688  dchord_step= pIntStepper-> DistChord();
689  // *********
690 
691  // Put back the values. yarrout ==> y_posvel
692  y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );
693  y_posvel.SetCurveLength( s_start + hstep );
694 
695 #ifdef G4DEBUG_FIELD
696  if(fVerboseLevel>2)
697  {
698  G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
699  PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep, 1);
700  }
701 #endif
702 
703  // A single measure of the error
704  // TO-DO : account for energy, spin, ... ?
705  vel_mag_sq = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
706  inv_vel_mag_sq = 1.0 / vel_mag_sq;
707  dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
708  dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
709  dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
710 
711  // Calculate also the change in the momentum squared also ???
712  // G4double veloc_square = y_posvel.GetVelocity().mag2();
713  // ...
714 
715 #ifdef RETURN_A_NEW_STEP_LENGTH
716  // The following step cannot be done here because "eps" is not known.
717  dyerr_len = std::sqrt( dyerr_len_sq );
718  dyerr_len_sq /= eps ;
719 
720  // Look at the velocity deviation ?
721  // sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
722 
723  // Set suggested new step
724  hstep= ComputeNewStepSize( dyerr_len, hstep);
725 #endif
726 
727  if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) )
728  {
729  dyerr = std::sqrt(dyerr_pos_sq);
730  }
731  else
732  {
733  // Scale it to the current step size - for now
734  dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
735  }
736 
737  return true;
738 }
static const G4double eps
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)
G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent)
#define G4endl
Definition: G4ios.hh:61
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4bool G4MagInt_Driver::QuickAdvance ( G4FieldTrack y_posvel,
const G4double  dydx[],
G4double  hstep,
G4double dchord_step,
G4double dyerr_pos_sq,
G4double dyerr_mom_rel_sq 
)

Definition at line 644 of file G4MagIntegratorDriver.cc.

651 {
652  G4Exception("G4MagInt_Driver::QuickAdvance()", "GeomField0001",
653  FatalException, "Not yet implemented.");
654 
655  // Use the parameters of this method, to please compiler
656  dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0];
657  dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
658  return true;
659 }
G4ThreeVector GetPosition() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double mag2() const

Here is the call graph for this function:

void G4MagInt_Driver::RenewStepperAndAdjust ( G4MagIntegratorStepper pItsStepper)
inline

Here is the caller graph for this function:

void G4MagInt_Driver::ReSetParameters ( G4double  new_safety = 0.9)
inline
void G4MagInt_Driver::SetErrcon ( G4double  valEc)
inline
void G4MagInt_Driver::SetHmin ( G4double  newval)
inline
void G4MagInt_Driver::SetMaxNoSteps ( G4int  val)
inline
void G4MagInt_Driver::SetPgrow ( G4double  valPg)
inline
void G4MagInt_Driver::SetPshrnk ( G4double  valPs)
inline
void G4MagInt_Driver::SetSafety ( G4double  valS)
inline
void G4MagInt_Driver::SetSmallestFraction ( G4double  val)

Definition at line 1043 of file G4MagIntegratorDriver.cc.

1044 {
1045  if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
1046  {
1047  fSmallestFraction= newFraction;
1048  }
1049  else
1050  {
1051  G4cerr << "Warning: SmallestFraction not changed. " << G4endl
1052  << " Proposed value was " << newFraction << G4endl
1053  << " Value must be between 1.e-8 and 1.e-16" << G4endl;
1054  }
1055 }
#define G4endl
Definition: G4ios.hh:61
G4GLOB_DLL std::ostream G4cerr
void G4MagInt_Driver::SetVerboseLevel ( G4int  newLevel)
inline

Here is the caller graph for this function:

void G4MagInt_Driver::WarnEndPointTooFar ( G4double  endPointDist,
G4double  hStepSize,
G4double  epsilonRelative,
G4int  debugFlag 
)
protected

Definition at line 481 of file G4MagIntegratorDriver.cc.

485 {
486  static G4ThreadLocal G4double maxRelError=0.0;
487  G4bool isNewMax, prNewMax;
488 
489  isNewMax = endPointDist > (1.0 + maxRelError) * h;
490  prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
491  if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
492 
494  && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+eps) ) ) )
495  {
496  static G4ThreadLocal G4int noWarnings = 0;
497  std::ostringstream message;
498  if( (noWarnings ++ < 10) || (dbg>2) )
499  {
500  message << "The integration produced an end-point which " << G4endl
501  << "is further from the start-point than the curve length."
502  << G4endl;
503  }
504  message << " Distance of endpoints = " << endPointDist
505  << ", curve length = " << h << G4endl
506  << " Difference (curveLen-endpDist)= " << (h - endPointDist)
507  << ", relative = " << (h-endPointDist) / h
508  << ", epsilon = " << eps;
509  G4Exception("G4MagInt_Driver::WarnEndPointTooFar()", "GeomField1001",
510  JustWarning, message);
511  }
512 }
G4double GetSurfaceTolerance() const
static const G4double eps
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
bool G4bool
Definition: G4Types.hh:79
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
static G4GeometryTolerance * GetInstance()

Here is the call graph for this function:

Here is the caller graph for this function:

void G4MagInt_Driver::WarnSmallStepSize ( G4double  hnext,
G4double  hstep,
G4double  h,
G4double  xDone,
G4int  noSteps 
)
protected

Definition at line 432 of file G4MagIntegratorDriver.cc.

435 {
436  static G4ThreadLocal G4int noWarningsIssued =0;
437  const G4int maxNoWarnings = 10; // Number of verbose warnings
438  std::ostringstream message;
439  if( (noWarningsIssued < maxNoWarnings) || fVerboseLevel > 10 )
440  {
441  message << "The stepsize for the next iteration, " << hnext
442  << ", is too small - in Step number " << nstp << "." << G4endl
443  << "The minimum for the driver is " << Hmin() << G4endl
444  << "Requested integr. length was " << hstep << " ." << G4endl
445  << "The size of this sub-step was " << h << " ." << G4endl
446  << "The integrations has already gone " << xDone;
447  }
448  else
449  {
450  message << "Too small 'next' step " << hnext
451  << ", step-no: " << nstp << G4endl
452  << ", this sub-step: " << h
453  << ", req_tot_len: " << hstep
454  << ", done: " << xDone << ", min: " << Hmin();
455  }
456  G4Exception("G4MagInt_Driver::WarnSmallStepSize()", "GeomField1001",
457  JustWarning, message);
458  noWarningsIssued++;
459 }
#define G4ThreadLocal
Definition: tls.hh:89
G4double Hmin() const
int G4int
Definition: G4Types.hh:78
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

void G4MagInt_Driver::WarnTooManyStep ( G4double  x1start,
G4double  x2end,
G4double  xCurrent 
)
protected

Definition at line 464 of file G4MagIntegratorDriver.cc.

467 {
468  std::ostringstream message;
469  message << "The number of steps used in the Integration driver"
470  << " (Runge-Kutta) is too many." << G4endl
471  << "Integration of the interval was not completed !" << G4endl
472  << "Only a " << (xCurrent-x1start)*100/(x2end-x1start)
473  << " % fraction of it was done.";
474  G4Exception("G4MagInt_Driver::WarnTooManyStep()", "GeomField1001",
475  JustWarning, message);
476 }
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:


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