Geant4  10.01.p01
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 81175 2014-05-22 07:39:10Z 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 //
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;
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 ?
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  // Have we reached the end ?
393  // --> a better test might be x-x2 > an_epsilon
394 
395  succeeded= (x>=x2); // If it was a "forced" last step
396 
397  for (i=0;i<nvar;i++) { yEnd[i] = y[i]; }
398 
399  // Put back the values.
400  y_current.LoadFromArray( yEnd, fNoIntegrationVariables );
401  y_current.SetCurveLength( x );
402 
403  if(nstp > fMaxNoSteps)
404  {
405  no_warnings++;
406  succeeded = false;
407 #ifdef G4DEBUG_FIELD
408  if (dbg)
409  {
410  WarnTooManyStep( x1, x2, x ); // Issue WARNING
411  PrintStatus( yEnd, x1, y, x, hstep, -nstp);
412  }
413 #endif
414  }
415 
416 #ifdef G4DEBUG_FIELD
417  if( dbg && no_warnings )
418  {
419  G4cerr << "G4MagIntegratorDriver exit status: no-steps " << nstp <<G4endl;
420  PrintStatus( yEnd, x1, y, x, hstep, nstp);
421  }
422 #endif
423 
424  return succeeded;
425 } // end of AccurateAdvance ...........................
426 
427 // ---------------------------------------------------------
428 
429 void
431  G4double h, G4double xDone,
432  G4int nstp)
433 {
434  static G4ThreadLocal G4int noWarningsIssued =0;
435  const G4int maxNoWarnings = 10; // Number of verbose warnings
436  std::ostringstream message;
437  if( (noWarningsIssued < maxNoWarnings) || fVerboseLevel > 10 )
438  {
439  message << "The stepsize for the next iteration, " << hnext
440  << ", is too small - in Step number " << nstp << "." << G4endl
441  << "The minimum for the driver is " << Hmin() << G4endl
442  << "Requested integr. length was " << hstep << " ." << G4endl
443  << "The size of this sub-step was " << h << " ." << G4endl
444  << "The integrations has already gone " << xDone;
445  }
446  else
447  {
448  message << "Too small 'next' step " << hnext
449  << ", step-no: " << nstp << G4endl
450  << ", this sub-step: " << h
451  << ", req_tot_len: " << hstep
452  << ", done: " << xDone << ", min: " << Hmin();
453  }
454  G4Exception("G4MagInt_Driver::WarnSmallStepSize()", "GeomField1001",
455  JustWarning, message);
456  noWarningsIssued++;
457 }
458 
459 // ---------------------------------------------------------
460 
461 void
463  G4double x2end,
464  G4double xCurrent)
465 {
466  std::ostringstream message;
467  message << "The number of steps used in the Integration driver"
468  << " (Runge-Kutta) is too many." << G4endl
469  << "Integration of the interval was not completed !" << G4endl
470  << "Only a " << (xCurrent-x1start)*100/(x2end-x1start)
471  << " % fraction of it was done.";
472  G4Exception("G4MagInt_Driver::WarnTooManyStep()", "GeomField1001",
473  JustWarning, message);
474 }
475 
476 // ---------------------------------------------------------
477 
478 void
480  G4double h ,
481  G4double eps,
482  G4int dbg)
483 {
484  static G4ThreadLocal G4double maxRelError=0.0;
485  G4bool isNewMax, prNewMax;
486 
487  isNewMax = endPointDist > (1.0 + maxRelError) * h;
488  prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
489  if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
490 
492  && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+eps) ) ) )
493  {
494  static G4ThreadLocal G4int noWarnings = 0;
495  std::ostringstream message;
496  if( (noWarnings ++ < 10) || (dbg>2) )
497  {
498  message << "The integration produced an end-point which " << G4endl
499  << "is further from the start-point than the curve length."
500  << G4endl;
501  }
502  message << " Distance of endpoints = " << endPointDist
503  << ", curve length = " << h << G4endl
504  << " Difference (curveLen-endpDist)= " << (h - endPointDist)
505  << ", relative = " << (h-endPointDist) / h
506  << ", epsilon = " << eps;
507  G4Exception("G4MagInt_Driver::WarnEndPointTooFar()", "GeomField1001",
508  JustWarning, message);
509  }
510 }
511 
512 // ---------------------------------------------------------
513 
514 void
516  const G4double dydx[],
517  G4double& x, // InOut
518  G4double htry,
519  G4double eps_rel_max,
520  G4double& hdid, // Out
521  G4double& hnext ) // Out
522 
523 // Driver for one Runge-Kutta Step with monitoring of local truncation error
524 // to ensure accuracy and adjust stepsize. Input are dependent variable
525 // array y[0,...,5] and its derivative dydx[0,...,5] at the
526 // starting value of the independent variable x . Also input are stepsize
527 // to be attempted htry, and the required accuracy eps. On output y and x
528 // are replaced by their new values, hdid is the stepsize that was actually
529 // accomplished, and hnext is the estimated next stepsize.
530 // This is similar to the function rkqs from the book:
531 // Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
532 // Edition, by William H. Press, Saul A. Teukolsky, William T.
533 // Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
534 // 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
535 
536 {
537  G4double errmax_sq;
538  G4double h, htemp, xnew ;
539 
541 
542  h = htry ; // Set stepsize to the initial trial value
543 
544  G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
545 
546  G4double errpos_sq=0.0; // square of displacement error
547  G4double errvel_sq=0.0; // square of momentum vector difference
548  G4double errspin_sq=0.0; // square of spin vector difference
549 
550  G4int iter;
551 
552  static G4ThreadLocal G4int tot_no_trials=0;
553  const G4int max_trials=100;
554 
555  G4ThreeVector Spin(y[9],y[10],y[11]);
556  G4double spin_mag2 =Spin.mag2() ;
557  G4bool hasSpin= (spin_mag2 > 0.0);
558 
559  for (iter=0; iter<max_trials ;iter++)
560  {
561  tot_no_trials++;
562  pIntStepper-> Stepper(y,dydx,h,ytemp,yerr);
563  // *******
564  G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep);
565  G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos);
566 
567  // Evaluate accuracy
568  //
569  errpos_sq = sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
570  errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
571 
572  // Accuracy for momentum
573  G4double magvel_sq= sqr(y[3]) + sqr(y[4]) + sqr(y[5]) ;
574  G4double sumerr_sq = sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) ;
575  if( magvel_sq > 0.0 ) {
576  errvel_sq = sumerr_sq / magvel_sq;
577  }else{
578  G4cerr << "** G4MagIntegrationDriver: found case of zero momentum."
579  << " iteration= " << iter << " h= " << h << G4endl;
580  errvel_sq = sumerr_sq;
581  }
582  errvel_sq *= inv_eps_vel_sq;
583  errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
584 
585  if( hasSpin )
586  {
587  // Accuracy for spin
588  errspin_sq = ( sqr(yerr[9]) + sqr(yerr[10]) + sqr(yerr[11]) )
589  / spin_mag2; // ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
590  errspin_sq *= inv_eps_vel_sq;
591  errmax_sq = std::max( errmax_sq, errspin_sq );
592  }
593 
594  if ( errmax_sq <= 1.0 ) { break; } // Step succeeded.
595 
596  // Step failed; compute the size of retrial Step.
597  htemp = GetSafety()*h* std::pow( errmax_sq, 0.5*GetPshrnk() );
598 
599  if (htemp >= 0.1*h) { h = htemp; } // Truncation error too large,
600  else { h = 0.1*h; } // reduce stepsize, but no more
601  // than a factor of 10
602  xnew = x + h;
603  if(xnew == x)
604  {
605  G4cerr << "G4MagIntegratorDriver::OneGoodStep:" << G4endl
606  << " Stepsize underflow in Stepper " << G4endl ;
607  G4cerr << " Step's start x=" << x << " and end x= " << xnew
608  << " are equal !! " << G4endl
609  <<" Due to step-size= " << h
610  << " . Note that input step was " << htry << G4endl;
611  break;
612  }
613  }
614 
615 #ifdef G4FLD_STATS
616  // Sum of squares of position error // and momentum dir (underestimated)
617  fSumH_lg += h;
618  fDyerrPos_lgTot += errpos_sq;
619  fDyerrVel_lgTot += errvel_sq * h * h;
620 #endif
621 
622  // Compute size of next Step
623  if (errmax_sq > errcon*errcon)
624  {
625  hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow());
626  }
627  else
628  {
629  hnext = max_stepping_increase*h ; // No more than a factor of 5 increase
630  }
631  x += (hdid = h);
632 
633  for(G4int k=0;k<fNoIntegrationVariables;k++) { y[k] = ytemp[k]; }
634 
635  return;
636 } // end of OneGoodStep .............................
637 
638 //----------------------------------------------------------------------
639 
640 // QuickAdvance just tries one Step - it does not ensure accuracy
641 //
643  G4FieldTrack& y_posvel, // INOUT
644  const G4double dydx[],
645  G4double hstep, // In
646  G4double& dchord_step,
647  G4double& dyerr_pos_sq,
648  G4double& dyerr_mom_rel_sq )
649 {
650  G4Exception("G4MagInt_Driver::QuickAdvance()", "GeomField0001",
651  FatalException, "Not yet implemented.");
652 
653  // Use the parameters of this method, to please compiler
654  dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0];
655  dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
656  return true;
657 }
658 
659 //----------------------------------------------------------------------
660 
662  G4FieldTrack& y_posvel, // INOUT
663  const G4double dydx[],
664  G4double hstep, // In
665  G4double& dchord_step,
666  G4double& dyerr )
667 {
668  G4double dyerr_pos_sq, dyerr_mom_rel_sq;
671  G4double s_start;
672  G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
673 
674  static G4ThreadLocal G4int no_call=0;
675  no_call ++;
676 
677  // Move data into array
678  y_posvel.DumpToArray( yarrin ); // yarrin <== y_posvel
679  s_start = y_posvel.GetCurveLength();
680 
681  // Do an Integration Step
682  pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ;
683  // *******
684 
685  // Estimate curve-chord distance
686  dchord_step= pIntStepper-> DistChord();
687  // *********
688 
689  // Put back the values. yarrout ==> y_posvel
690  y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );
691  y_posvel.SetCurveLength( s_start + hstep );
692 
693 #ifdef G4DEBUG_FIELD
694  if(fVerboseLevel>2)
695  {
696  G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
697  PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep, 1);
698  }
699 #endif
700 
701  // A single measure of the error
702  // TO-DO : account for energy, spin, ... ?
703  vel_mag_sq = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
704  inv_vel_mag_sq = 1.0 / vel_mag_sq;
705  dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
706  dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
707  dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
708 
709  // Calculate also the change in the momentum squared also ???
710  // G4double veloc_square = y_posvel.GetVelocity().mag2();
711  // ...
712 
713 #ifdef RETURN_A_NEW_STEP_LENGTH
714  // The following step cannot be done here because "eps" is not known.
715  dyerr_len = std::sqrt( dyerr_len_sq );
716  dyerr_len_sq /= eps ;
717 
718  // Look at the velocity deviation ?
719  // sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
720 
721  // Set suggested new step
722  hstep= ComputeNewStepSize( dyerr_len, hstep);
723 #endif
724 
725  if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) )
726  {
727  dyerr = std::sqrt(dyerr_pos_sq);
728  }
729  else
730  {
731  // Scale it to the current step size - for now
732  dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
733  }
734 
735  return true;
736 }
737 
738 // --------------------------------------------------------------------------
739 
740 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT
742  G4double yarrin[], // In
743  const G4double dydx[],
744  G4double hstep, // In
745  G4double yarrout[],
746  G4double& dchord_step,
747  G4double& dyerr ) // In length
748 {
749  G4Exception("G4MagInt_Driver::QuickAdvance()", "GeomField0001",
750  FatalException, "Not yet implemented.");
751  dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
752  yarrout[0]= yarrin[0];
753 }
754 #endif
755 
756 // --------------------------------------------------------------------------
757 
758 // This method computes new step sizes - but does not limit changes to
759 // within certain factors
760 //
761 G4double
763  G4double errMaxNorm, // max error (normalised)
764  G4double hstepCurrent) // current step size
765 {
766  G4double hnew;
767 
768  // Compute size of next Step for a failed step
769  if(errMaxNorm > 1.0 )
770  {
771  // Step failed; compute the size of retrial Step.
772  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
773  } else if(errMaxNorm > 0.0 ) {
774  // Compute size of next Step for a successful step
775  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
776  } else {
777  // if error estimate is zero (possible) or negative (dubious)
778  hnew = max_stepping_increase * hstepCurrent;
779  }
780 
781  return hnew;
782 }
783 
784 // ---------------------------------------------------------------------------
785 
786 // This method computes new step sizes limiting changes within certain factors
787 //
788 // It shares its logic with AccurateAdvance.
789 // They are kept separate currently for optimisation.
790 //
791 G4double
793  G4double errMaxNorm, // max error (normalised)
794  G4double hstepCurrent) // current step size
795 {
796  G4double hnew;
797 
798  // Compute size of next Step for a failed step
799  if (errMaxNorm > 1.0 )
800  {
801  // Step failed; compute the size of retrial Step.
802  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
803 
804  if (hnew < max_stepping_decrease*hstepCurrent)
805  {
806  hnew = max_stepping_decrease*hstepCurrent ;
807  // reduce stepsize, but no more
808  // than this factor (value= 1/10)
809  }
810  }
811  else
812  {
813  // Compute size of next Step for a successful step
814  if (errMaxNorm > errcon)
815  { hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()); }
816  else // No more than a factor of 5 increase
817  { hnew = max_stepping_increase * hstepCurrent; }
818  }
819  return hnew;
820 }
821 
822 // ---------------------------------------------------------------------------
823 
824 void G4MagInt_Driver::PrintStatus( const G4double* StartArr,
825  G4double xstart,
826  const G4double* CurrentArr,
827  G4double xcurrent,
828  G4double requestStep,
829  G4int subStepNo)
830  // Potentially add as arguments:
831  // <dydx> - as Initial Force
832  // stepTaken(hdid) - last step taken
833  // nextStep (hnext) - proposal for size
834 {
835  G4FieldTrack StartFT(G4ThreeVector(0,0,0),
836  G4ThreeVector(0,0,0), 0., 0., 0., 0. );
837  G4FieldTrack CurrentFT (StartFT);
838 
839  StartFT.LoadFromArray( StartArr, fNoIntegrationVariables);
840  StartFT.SetCurveLength( xstart);
841  CurrentFT.LoadFromArray( CurrentArr, fNoIntegrationVariables);
842  CurrentFT.SetCurveLength( xcurrent );
843 
844  PrintStatus(StartFT, CurrentFT, requestStep, subStepNo );
845 }
846 
847 // ---------------------------------------------------------------------------
848 
850  const G4FieldTrack& StartFT,
851  const G4FieldTrack& CurrentFT,
852  G4double requestStep,
853  G4int subStepNo)
854 {
855  G4int verboseLevel= fVerboseLevel;
856  static G4ThreadLocal G4int noPrecision= 5;
857  G4int oldPrec= G4cout.precision(noPrecision);
858  // G4cout.setf(ios_base::fixed,ios_base::floatfield);
859 
860  const G4ThreeVector StartPosition= StartFT.GetPosition();
861  const G4ThreeVector StartUnitVelocity= StartFT.GetMomentumDir();
862  const G4ThreeVector CurrentPosition= CurrentFT.GetPosition();
863  const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir();
864 
865  G4double DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
866 
867  G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
868  G4double subStepSize = step_len;
869 
870  if( (subStepNo <= 1) || (verboseLevel > 3) )
871  {
872  subStepNo = - subStepNo; // To allow printing banner
873 
874  G4cout << std::setw( 6) << " " << std::setw( 25)
875  << " G4MagInt_Driver: Current Position and Direction" << " "
876  << G4endl;
877  G4cout << std::setw( 5) << "Step#" << " "
878  << std::setw( 7) << "s-curve" << " "
879  << std::setw( 9) << "X(mm)" << " "
880  << std::setw( 9) << "Y(mm)" << " "
881  << std::setw( 9) << "Z(mm)" << " "
882  << std::setw( 8) << " N_x " << " "
883  << std::setw( 8) << " N_y " << " "
884  << std::setw( 8) << " N_z " << " "
885  << std::setw( 8) << " N^2-1 " << " "
886  << std::setw(10) << " N(0).N " << " "
887  << std::setw( 7) << "KinEner " << " "
888  << std::setw(12) << "Track-l" << " " // Add the Sub-step ??
889  << std::setw(12) << "Step-len" << " "
890  << std::setw(12) << "Step-len" << " "
891  << std::setw( 9) << "ReqStep" << " "
892  << G4endl;
893  }
894 
895  if( (subStepNo <= 0) )
896  {
897  PrintStat_Aux( StartFT, requestStep, 0.,
898  0, 0.0, 1.0);
899  //*************
900  }
901 
902  if( verboseLevel <= 3 )
903  {
904  G4cout.precision(noPrecision);
905  PrintStat_Aux( CurrentFT, requestStep, step_len,
906  subStepNo, subStepSize, DotStartCurrentVeloc );
907  //*************
908  }
909 
910  else // if( verboseLevel > 3 )
911  {
912  // Multi-line output
913 
914  // G4cout << "Current Position is " << CurrentPosition << G4endl
915  // << " and UnitVelocity is " << CurrentUnitVelocity << G4endl;
916  // G4cout << "Step taken was " << step_len
917  // << " out of PhysicalStep= " << requestStep << G4endl;
918  // G4cout << "Final safety is: " << safety << G4endl;
919  // G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
920  // << G4endl << G4endl;
921  }
922  G4cout.precision(oldPrec);
923 }
924 
925 // ---------------------------------------------------------------------------
926 
928  const G4FieldTrack& aFieldTrack,
929  G4double requestStep,
930  G4double step_len,
931  G4int subStepNo,
932  G4double subStepSize,
933  G4double dotVeloc_StartCurr)
934 {
935  const G4ThreeVector Position= aFieldTrack.GetPosition();
936  const G4ThreeVector UnitVelocity= aFieldTrack.GetMomentumDir();
937 
938  if( subStepNo >= 0)
939  {
940  G4cout << std::setw( 5) << subStepNo << " ";
941  }
942  else
943  {
944  G4cout << std::setw( 5) << "Start" << " ";
945  }
946  G4double curveLen= aFieldTrack.GetCurveLength();
947  G4cout << std::setw( 7) << curveLen;
948  G4cout << std::setw( 9) << Position.x() << " "
949  << std::setw( 9) << Position.y() << " "
950  << std::setw( 9) << Position.z() << " "
951  << std::setw( 8) << UnitVelocity.x() << " "
952  << std::setw( 8) << UnitVelocity.y() << " "
953  << std::setw( 8) << UnitVelocity.z() << " ";
954  G4int oldprec= G4cout.precision(3);
955  G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 << " ";
956  G4cout.precision(6);
957  G4cout << std::setw(10) << dotVeloc_StartCurr << " ";
958  G4cout.precision(oldprec);
959  G4cout << std::setw( 7) << aFieldTrack.GetKineticEnergy();
960  G4cout << std::setw(12) << step_len << " ";
961 
962  static G4ThreadLocal G4double oldCurveLength= 0.0;
963  static G4ThreadLocal G4double oldSubStepLength= 0.0;
964  static G4ThreadLocal G4int oldSubStepNo= -1;
965 
966  G4double subStep_len=0.0;
967  if( curveLen > oldCurveLength )
968  {
969  subStep_len= curveLen - oldCurveLength;
970  }
971  else if (subStepNo == oldSubStepNo)
972  {
973  subStep_len= oldSubStepLength;
974  }
975  oldCurveLength= curveLen;
976  oldSubStepLength= subStep_len;
977 
978  G4cout << std::setw(12) << subStep_len << " ";
979  G4cout << std::setw(12) << subStepSize << " ";
980  if( requestStep != -1.0 )
981  {
982  G4cout << std::setw( 9) << requestStep << " ";
983  }
984  else
985  {
986  G4cout << std::setw( 9) << " InitialStep " << " ";
987  }
988  G4cout << G4endl;
989 }
990 
991 // ---------------------------------------------------------------------------
992 
994 {
995  G4int noPrecBig= 6;
996  G4int oldPrec= G4cout.precision(noPrecBig);
997 
998  G4cout << "G4MagInt_Driver Statistics of steps undertaken. " << G4endl;
999  G4cout << "G4MagInt_Driver: Number of Steps: "
1000  << " Total= " << fNoTotalSteps
1001  << " Bad= " << fNoBadSteps
1002  << " Small= " << fNoSmallSteps
1003  << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
1004  << G4endl;
1005 
1006 #ifdef G4FLD_STATS
1007  G4cout << "MID dyerr: "
1008  << " maximum= " << fDyerr_max
1009  << " Sum small= " << fDyerrPos_smTot
1010  << " std::sqrt(Sum large^2): pos= " << std::sqrt(fDyerrPos_lgTot)
1011  << " vel= " << std::sqrt( fDyerrVel_lgTot )
1012  << " Total h-distance: small= " << fSumH_sm
1013  << " large= " << fSumH_lg
1014  << G4endl;
1015 
1016 #if 0
1017  G4int noPrecSmall=4;
1018  // Single line precis of statistics ... optional
1019  G4cout.precision(noPrecSmall);
1020  G4cout << "MIDnums: " << fMinimumStep
1021  << " " << fNoTotalSteps
1022  << " " << fNoSmallSteps
1024  << " " << fNoBadSteps
1025  << " " << fDyerr_max
1026  << " " << fDyerr_mx2
1027  << " " << fDyerrPos_smTot
1028  << " " << fSumH_sm
1029  << " " << fDyerrPos_lgTot
1030  << " " << fDyerrVel_lgTot
1031  << " " << fSumH_lg
1032  << G4endl;
1033 #endif
1034 #endif
1035 
1036  G4cout.precision(oldPrec);
1037 }
1038 
1039 // ---------------------------------------------------------------------------
1040 
1042 {
1043  if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
1044  {
1045  fSmallestFraction= newFraction;
1046  }
1047  else
1048  {
1049  G4cerr << "Warning: SmallestFraction not changed. " << G4endl
1050  << " Proposed value was " << newFraction << G4endl
1051  << " Value must be between 1.e-8 and 1.e-16" << G4endl;
1052  }
1053 }
virtual void ComputeRightHandSide(const G4double y[], G4double dydx[])
G4double GetCurveLength() const
CLHEP::Hep3Vector G4ThreeVector
G4MagIntegratorStepper * pIntStepper
G4double GetKineticEnergy() 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
static const double perThousand
Definition: G4SIunits.hh:297
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
static const G4double max_stepping_increase
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
static const double perMillion
Definition: G4SIunits.hh:298
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)
void SetSmallestFraction(G4double val)
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
static const G4int fMaxStepBase
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
static G4GeometryTolerance * GetInstance()
G4GLOB_DLL std::ostream G4cerr
const G4int fNoIntegrationVariables
static const G4double max_stepping_decrease