Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MagIntegratorDriver.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4MagIntegratorDriver.cc 101384 2016-11-16 11:03:44Z gcosmo $
28 //
29 //
30 //
31 // Implementation for class G4MagInt_Driver
32 // Tracking in space dependent magnetic field
33 //
34 // History of major changes:
35 // 8 Nov 01 J. Apostolakis: Respect minimum step in AccurateAdvance
36 // 27 Jul 99 J. Apostolakis: Ensured that AccurateAdvance does not loop
37 // due to very small eps & step size (precision)
38 // 28 Jan 98 W. Wander: Added ability for low order integrators
39 // 7 Oct 96 V. Grichine First version
40 // --------------------------------------------------------------------
41 
42 #include <iomanip>
43 
44 #include "globals.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4GeometryTolerance.hh"
47 #include "G4MagIntegratorDriver.hh"
48 #include "G4FieldTrack.hh"
49 
50 // Stepsize can increase by no more than 5.0
51 // and decrease by no more than 1/10. = 0.1
52 //
53 const G4double G4MagInt_Driver::max_stepping_increase = 5.0;
54 const G4double G4MagInt_Driver::max_stepping_decrease = 0.1;
55 
56 // The (default) maximum number of steps is Base
57 // divided by the order of Stepper
58 //
59 const G4int G4MagInt_Driver::fMaxStepBase = 250; // Was 5000
60 
61 #ifndef G4NO_FIELD_STATISTICS
62 #define G4FLD_STATS 1
63 #endif
64 
65 // ---------------------------------------------------------
66 
67 // Constructor
68 //
70  G4MagIntegratorStepper *pStepper,
71  G4int numComponents,
72  G4int statisticsVerbose)
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 }
107 
108 // ---------------------------------------------------------
109 
110 // Destructor
111 //
113 {
114  if( fStatisticsVerboseLevel > 1 )
115  {
117  }
118 }
119 
120 // To add much printing for debugging purposes, uncomment the following
121 // and set verbose level to 1 or higher value !
122 // #define G4DEBUG_FIELD 1
123 
124 // ---------------------------------------------------------
125 
126 G4bool
128  G4double hstep,
129  G4double eps,
130  G4double hinitial )
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 ...........................
428 
429 // ---------------------------------------------------------
430 
431 void
433  G4double h, G4double xDone,
434  G4int nstp)
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 }
460 
461 // ---------------------------------------------------------
462 
463 void
465  G4double x2end,
466  G4double xCurrent)
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 }
477 
478 // ---------------------------------------------------------
479 
480 void
482  G4double h ,
483  G4double eps,
484  G4int dbg)
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 }
513 
514 // ---------------------------------------------------------
515 
516 void
518  const G4double dydx[],
519  G4double& x, // InOut
520  G4double htry,
521  G4double eps_rel_max,
522  G4double& hdid, // Out
523  G4double& hnext ) // Out
524 
525 // Driver for one Runge-Kutta Step with monitoring of local truncation error
526 // to ensure accuracy and adjust stepsize. Input are dependent variable
527 // array y[0,...,5] and its derivative dydx[0,...,5] at the
528 // starting value of the independent variable x . Also input are stepsize
529 // to be attempted htry, and the required accuracy eps. On output y and x
530 // are replaced by their new values, hdid is the stepsize that was actually
531 // accomplished, and hnext is the estimated next stepsize.
532 // This is similar to the function rkqs from the book:
533 // Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
534 // Edition, by William H. Press, Saul A. Teukolsky, William T.
535 // Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
536 // 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
537 
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 .............................
639 
640 //----------------------------------------------------------------------
641 
642 // QuickAdvance just tries one Step - it does not ensure accuracy
643 //
645  G4FieldTrack& y_posvel, // INOUT
646  const G4double dydx[],
647  G4double hstep, // In
648  G4double& dchord_step,
649  G4double& dyerr_pos_sq,
650  G4double& dyerr_mom_rel_sq )
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 }
660 
661 //----------------------------------------------------------------------
662 
664  G4FieldTrack& y_posvel, // INOUT
665  const G4double dydx[],
666  G4double hstep, // In
667  G4double& dchord_step,
668  G4double& dyerr )
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 }
739 
740 // --------------------------------------------------------------------------
741 
742 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT
744  G4double yarrin[], // In
745  const G4double dydx[],
746  G4double hstep, // In
747  G4double yarrout[],
748  G4double& dchord_step,
749  G4double& dyerr ) // In length
750 {
751  G4Exception("G4MagInt_Driver::QuickAdvance()", "GeomField0001",
752  FatalException, "Not yet implemented.");
753  dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
754  yarrout[0]= yarrin[0];
755 }
756 #endif
757 
758 // --------------------------------------------------------------------------
759 
760 // This method computes new step sizes - but does not limit changes to
761 // within certain factors
762 //
763 G4double
765  G4double errMaxNorm, // max error (normalised)
766  G4double hstepCurrent) // current step size
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 }
785 
786 // ---------------------------------------------------------------------------
787 
788 // This method computes new step sizes limiting changes within certain factors
789 //
790 // It shares its logic with AccurateAdvance.
791 // They are kept separate currently for optimisation.
792 //
793 G4double
795  G4double errMaxNorm, // max error (normalised)
796  G4double hstepCurrent) // current step size
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 }
823 
824 // ---------------------------------------------------------------------------
825 
826 void G4MagInt_Driver::PrintStatus( const G4double* StartArr,
827  G4double xstart,
828  const G4double* CurrentArr,
829  G4double xcurrent,
830  G4double requestStep,
831  G4int subStepNo)
832  // Potentially add as arguments:
833  // <dydx> - as Initial Force
834  // stepTaken(hdid) - last step taken
835  // nextStep (hnext) - proposal for size
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 }
848 
849 // ---------------------------------------------------------------------------
850 
852  const G4FieldTrack& StartFT,
853  const G4FieldTrack& CurrentFT,
854  G4double requestStep,
855  G4int subStepNo)
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 }
926 
927 // ---------------------------------------------------------------------------
928 
930  const G4FieldTrack& aFieldTrack,
931  G4double requestStep,
932  G4double step_len,
933  G4int subStepNo,
934  G4double subStepSize,
935  G4double dotVeloc_StartCurr)
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 }
992 
993 // ---------------------------------------------------------------------------
994 
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 }
1040 
1041 // ---------------------------------------------------------------------------
1042 
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 }
static constexpr double perMillion
Definition: G4SIunits.hh:334
virtual void ComputeRightHandSide(const G4double y[], G4double dydx[])
G4double GetCurveLength() const
CLHEP::Hep3Vector G4ThreeVector
double x() const
G4double GetKineticEnergy() const
double dot(const Hep3Vector &) const
void SetCurveLength(G4double nCurve_s)
void RenewStepperAndAdjust(G4MagIntegratorStepper *pItsStepper)
void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent)
void WarnEndPointTooFar(G4double endPointDist, G4double hStepSize, G4double epsilonRelative, G4int debugFlag)
G4double GetSurfaceTolerance() const
const G4ThreeVector & GetMomentumDir() const
static const G4double eps
tuple x
Definition: test.py:50
G4MagInt_Driver(G4double hminimum, G4MagIntegratorStepper *pItsStepper, G4int numberOfComponents=6, G4int statisticsVerbosity=1)
#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)
G4double GetSafety() const
double z() const
static const double prec
Definition: RanecuEngine.cc:58
virtual G4int IntegratorOrder() const =0
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)
G4double ComputeNewStepSize_WithinLimits(G4double errMaxNorm, G4double hstepCurrent)
bool G4bool
Definition: G4Types.hh:79
G4double GetPgrow() const
void PrintStat_Aux(const G4FieldTrack &aFieldTrack, G4double requestStep, G4double actualStep, G4int subStepNo, G4double subStepSize, G4double dotVelocities)
G4double GetPshrnk() const
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)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void LoadFromArray(const G4double valArr[ncompSVEC], G4int noVarsIntegrated)
double y() const
void SetSmallestFraction(G4double val)
double mag2() const
G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent)
#define G4endl
Definition: G4ios.hh:61
static constexpr double perThousand
Definition: G4SIunits.hh:333
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
static G4GeometryTolerance * GetInstance()
G4GLOB_DLL std::ostream G4cerr