ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MagIntegratorDriver.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MagIntegratorDriver.cc
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 // G4MagInt_Driver implementation
27 //
28 // V.Grichine, 07.10.1996 - Created
29 // W.Wander, 28.01.1998 - Added ability for low order integrators
30 // J.Apostolakis, 08.11.2001 - Respect minimum step in AccurateAdvance
31 // --------------------------------------------------------------------
32 
33 #include <iomanip>
34 
35 #include "globals.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4GeometryTolerance.hh"
38 #include "G4MagIntegratorDriver.hh"
39 #include "G4FieldTrack.hh"
40 
41 // ---------------------------------------------------------
42 
43 // Constructor
44 //
45 G4MagInt_Driver::G4MagInt_Driver( G4double hminimum,
46  G4MagIntegratorStepper* pStepper,
47  G4int numComponents,
48  G4int statisticsVerbose)
49  : fNoIntegrationVariables(numComponents),
50  fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )),
51  fStatisticsVerboseLevel(statisticsVerbose)
52 {
53  // In order to accomodate "Laboratory Time", which is [7], fMinNoVars=8
54  // is required. For proper time of flight and spin, fMinNoVars must be 12
55 
56  RenewStepperAndAdjust( pStepper );
57  fMinimumStep = hminimum;
58 
59  fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
60 #ifdef G4DEBUG_FIELD
61  fVerboseLevel=2;
62 #endif
63 
64  if( (fVerboseLevel > 0) || (fStatisticsVerboseLevel > 1) )
65  {
66  G4cout << "MagIntDriver version: Accur-Adv: "
67  << "invE_nS, QuickAdv-2sqrt with Statistics "
68 #ifdef G4FLD_STATS
69  << " enabled "
70 #else
71  << " disabled "
72 #endif
73  << G4endl;
74  }
75 }
76 
77 // ---------------------------------------------------------
78 
79 // Destructor
80 //
81 G4MagInt_Driver::~G4MagInt_Driver()
82 {
83  if( fStatisticsVerboseLevel > 1 )
84  {
85  PrintStatisticsReport();
86  }
87 }
88 
89 // ---------------------------------------------------------
90 
91 G4bool
92 G4MagInt_Driver::AccurateAdvance(G4FieldTrack& y_current,
93  G4double hstep,
94  G4double eps,
95  G4double hinitial )
96 {
97  // Runge-Kutta driver with adaptive stepsize control. Integrate starting
98  // values at y_current over hstep x2 with accuracy eps.
99  // On output ystart is replaced by values at the end of the integration
100  // interval. RightHandSide is the right-hand side of ODE system.
101  // The source is similar to odeint routine from NRC p.721-722 .
102 
103  G4int nstp, i, no_warnings = 0;
104  G4double x, hnext, hdid, h;
105 
106 #ifdef G4DEBUG_FIELD
107  static G4int dbg = 1;
108  static G4int nStpPr = 50; // For debug printing of long integrations
109  G4double ySubStepStart[G4FieldTrack::ncompSVEC];
110  G4FieldTrack yFldTrkStart(y_current);
111 #endif
112 
113  G4double y[G4FieldTrack::ncompSVEC], dydx[G4FieldTrack::ncompSVEC];
114  G4double ystart[G4FieldTrack::ncompSVEC], yEnd[G4FieldTrack::ncompSVEC];
115  G4double x1, x2;
116  G4bool succeeded = true, lastStepSucceeded;
117 
118  G4double startCurveLength;
119 
120  G4int noFullIntegr = 0, noSmallIntegr = 0;
121  static G4ThreadLocal G4int noGoodSteps = 0; // Bad = chord > curve-len
122  const G4int nvar = fNoVars;
123 
124  G4FieldTrack yStartFT(y_current);
125 
126  // Ensure that hstep > 0
127  //
128  if( hstep <= 0.0 )
129  {
130  if( hstep == 0.0 )
131  {
132  std::ostringstream message;
133  message << "Proposed step is zero; hstep = " << hstep << " !";
134  G4Exception("G4MagInt_Driver::AccurateAdvance()",
135  "GeomField1001", JustWarning, message);
136  return succeeded;
137  }
138  else
139  {
140  std::ostringstream message;
141  message << "Invalid run condition." << G4endl
142  << "Proposed step is negative; hstep = " << hstep << "." << G4endl
143  << "Requested step cannot be negative! Aborting event.";
144  G4Exception("G4MagInt_Driver::AccurateAdvance()",
145  "GeomField0003", EventMustBeAborted, message);
146  return false;
147  }
148  }
149 
150  y_current.DumpToArray( ystart );
151 
152  startCurveLength= y_current.GetCurveLength();
153  x1= startCurveLength;
154  x2= x1 + hstep;
155 
156  if ( (hinitial > 0.0) && (hinitial < hstep)
157  && (hinitial > perMillion * hstep) )
158  {
159  h = hinitial;
160  }
161  else // Initial Step size "h" defaults to the full interval
162  {
163  h = hstep;
164  }
165 
166  x = x1;
167 
168  for ( i=0; i<nvar; ++i) { y[i] = ystart[i]; }
169 
170  G4bool lastStep= false;
171  nstp = 1;
172 
173  do
174  {
175  G4ThreeVector StartPos( y[0], y[1], y[2] );
176 
177 #ifdef G4DEBUG_FIELD
178  G4double xSubStepStart= x;
179  for (i=0; i<nvar; ++i) { ySubStepStart[i] = y[i]; }
180  yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
181  yFldTrkStart.SetCurveLength(x);
182 #endif
183 
184  pIntStepper->RightHandSide( y, dydx );
185  ++fNoTotalSteps;
186 
187  // Perform the Integration
188  //
189  if( h > fMinimumStep )
190  {
191  OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ;
192  //--------------------------------------
193  lastStepSucceeded = (hdid == h);
194 #ifdef G4DEBUG_FIELD
195  if (dbg>2)
196  {
197  PrintStatus( ySubStepStart, xSubStepStart, y, x, h, nstp); // Only
198  }
199 #endif
200  }
201  else
202  {
203  G4FieldTrack yFldTrk( G4ThreeVector(0,0,0),
204  G4ThreeVector(0,0,0), 0., 0., 0., 0. );
205  G4double dchord_step, dyerr, dyerr_len; // What to do with these ?
206  yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
207  yFldTrk.SetCurveLength( x );
208 
209  QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
210  //-----------------------------------------------------
211 
212  yFldTrk.DumpToArray(y);
213 
214 #ifdef G4FLD_STATS
215  ++fNoSmallSteps;
216  if ( dyerr_len > fDyerr_max ) { fDyerr_max = dyerr_len; }
217  fDyerrPos_smTot += dyerr_len;
218  fSumH_sm += h; // Length total for 'small' steps
219  if (nstp<=1) { ++fNoInitialSmallSteps; }
220 #endif
221 #ifdef G4DEBUG_FIELD
222  if (dbg>1)
223  {
224  if(fNoSmallSteps<2) { PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
225  G4cout << "Another sub-min step, no " << fNoSmallSteps
226  << " of " << fNoTotalSteps << " this time " << nstp << G4endl;
227  PrintStatus( ySubStepStart, x1, y, x, h, nstp); // Only this
228  G4cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h
229  << " epsilon= " << eps << " hstep= " << hstep
230  << " h= " << h << " hmin= " << fMinimumStep << G4endl;
231  }
232 #endif
233  if( h == 0.0 )
234  {
235  G4Exception("G4MagInt_Driver::AccurateAdvance()",
236  "GeomField0003", FatalException,
237  "Integration Step became Zero!");
238  }
239  dyerr = dyerr_len / h;
240  hdid = h;
241  x += hdid;
242 
243  // Compute suggested new step
244  hnext = ComputeNewStepSize( dyerr/eps, h);
245 
246  // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
247  lastStepSucceeded = (dyerr<= eps);
248  }
249 
250  if (lastStepSucceeded) { ++noFullIntegr; }
251  else { ++noSmallIntegr; }
252 
253  G4ThreeVector EndPos( y[0], y[1], y[2] );
254 
255 #ifdef G4DEBUG_FIELD
256  if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
257  {
258  if( nstp==nStpPr ) { G4cout << "***** Many steps ****" << G4endl; }
259  G4cout << "MagIntDrv: " ;
260  G4cout << "hdid=" << std::setw(12) << hdid << " "
261  << "hnext=" << std::setw(12) << hnext << " "
262  << "hstep=" << std::setw(12) << hstep << " (requested) "
263  << G4endl;
264  PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
265  }
266 #endif
267 
268  // Check the endpoint
269  G4double endPointDist= (EndPos-StartPos).mag();
270  if ( endPointDist >= hdid*(1.+perMillion) )
271  {
272  ++fNoBadSteps;
273 
274  // Issue a warning only for gross differences -
275  // we understand how small difference occur.
276  if ( endPointDist >= hdid*(1.+perThousand) )
277  {
278 #ifdef G4DEBUG_FIELD
279  if (dbg)
280  {
281  WarnEndPointTooFar ( endPointDist, hdid, eps, dbg );
282  G4cerr << " Total steps: bad " << fNoBadSteps
283  << " good " << noGoodSteps << " current h= " << hdid
284  << G4endl;
285  PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
286  }
287 #endif
288  ++no_warnings;
289  }
290  }
291  else
292  {
293  ++noGoodSteps;
294  }
295 // #endif
296 
297  // Avoid numerous small last steps
298  if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
299  {
300  // No more integration -- the next step will not happen
301  lastStep = true;
302  }
303  else
304  {
305  // Check the proposed next stepsize
306  if(std::fabs(hnext) <= Hmin())
307  {
308 #ifdef G4DEBUG_FIELD
309  // If simply a very small interval is being integrated, do not warn
310  if( (x < x2 * (1-eps) ) && // The last step can be small: OK
311  (std::fabs(hstep) > Hmin()) ) // and if we are asked, it's OK
312  {
313  if(dbg>0)
314  {
315  WarnSmallStepSize( hnext, hstep, h, x-x1, nstp );
316  PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
317  }
318  ++no_warnings;
319  }
320 #endif
321  // Make sure that the next step is at least Hmin.
322  h = Hmin();
323  }
324  else
325  {
326  h = hnext;
327  }
328 
329  // Ensure that the next step does not overshoot
330  if ( x+h > x2 )
331  { // When stepsize overshoots, decrease it!
332  h = x2 - x ; // Must cope with difficult rounding-error
333  } // issues if hstep << x2
334 
335  if ( h == 0.0 )
336  {
337  // Cannot progress - accept this as last step - by default
338  lastStep = true;
339 #ifdef G4DEBUG_FIELD
340  if (dbg>2)
341  {
342  int prec= G4cout.precision(12);
343  G4cout << "Warning: G4MagIntegratorDriver::AccurateAdvance"
344  << G4endl
345  << " Integration step 'h' became "
346  << h << " due to roundoff. " << G4endl
347  << " Calculated as difference of x2= "<< x2 << " and x=" << x
348  << " Forcing termination of advance." << G4endl;
349  G4cout.precision(prec);
350  }
351 #endif
352  }
353  }
354  } while ( ((nstp++)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
355  // Loop checking, 07.10.2016, J. Apostolakis
356 
357  // Have we reached the end ?
358  // --> a better test might be x-x2 > an_epsilon
359 
360  succeeded = (x>=x2); // If it was a "forced" last step
361 
362  for (i=0; i<nvar; ++i) { yEnd[i] = y[i]; }
363 
364  // Put back the values.
365  y_current.LoadFromArray( yEnd, fNoIntegrationVariables );
366  y_current.SetCurveLength( x );
367 
368  if(nstp > fMaxNoSteps)
369  {
370  ++no_warnings;
371  succeeded = false;
372 #ifdef G4DEBUG_FIELD
373  if (dbg)
374  {
375  WarnTooManyStep( x1, x2, x ); // Issue WARNING
376  PrintStatus( yEnd, x1, y, x, hstep, -nstp);
377  }
378 #endif
379  }
380 
381 #ifdef G4DEBUG_FIELD
382  if( dbg && no_warnings )
383  {
384  G4cerr << "G4MagIntegratorDriver exit status: no-steps " << nstp << G4endl;
385  PrintStatus( yEnd, x1, y, x, hstep, nstp);
386  }
387 #endif
388 
389  return succeeded;
390 } // end of AccurateAdvance ...........................
391 
392 // ---------------------------------------------------------
393 
394 void
395 G4MagInt_Driver::WarnSmallStepSize( G4double hnext, G4double hstep,
396  G4double h, G4double xDone,
397  G4int nstp)
398 {
399  static G4ThreadLocal G4int noWarningsIssued = 0;
400  const G4int maxNoWarnings = 10; // Number of verbose warnings
401  std::ostringstream message;
402  if( (noWarningsIssued < maxNoWarnings) || fVerboseLevel > 10 )
403  {
404  message << "The stepsize for the next iteration, " << hnext
405  << ", is too small - in Step number " << nstp << "." << G4endl
406  << "The minimum for the driver is " << Hmin() << G4endl
407  << "Requested integr. length was " << hstep << " ." << G4endl
408  << "The size of this sub-step was " << h << " ." << G4endl
409  << "The integrations has already gone " << xDone;
410  }
411  else
412  {
413  message << "Too small 'next' step " << hnext
414  << ", step-no: " << nstp << G4endl
415  << ", this sub-step: " << h
416  << ", req_tot_len: " << hstep
417  << ", done: " << xDone << ", min: " << Hmin();
418  }
419  G4Exception("G4MagInt_Driver::WarnSmallStepSize()", "GeomField1001",
420  JustWarning, message);
421  ++noWarningsIssued;
422 }
423 
424 // ---------------------------------------------------------
425 
426 void
427 G4MagInt_Driver::WarnTooManyStep( G4double x1start,
428  G4double x2end,
429  G4double xCurrent )
430 {
431  std::ostringstream message;
432  message << "The number of steps used in the Integration driver"
433  << " (Runge-Kutta) is too many." << G4endl
434  << "Integration of the interval was not completed !" << G4endl
435  << "Only a " << (xCurrent-x1start)*100/(x2end-x1start)
436  << " % fraction of it was done.";
437  G4Exception("G4MagInt_Driver::WarnTooManyStep()", "GeomField1001",
438  JustWarning, message);
439 }
440 
441 // ---------------------------------------------------------
442 
443 void
444 G4MagInt_Driver::WarnEndPointTooFar (G4double endPointDist,
445  G4double h ,
446  G4double eps,
447  G4int dbg)
448 {
449  static G4ThreadLocal G4double maxRelError = 0.0;
450  G4bool isNewMax, prNewMax;
451 
452  isNewMax = endPointDist > (1.0 + maxRelError) * h;
453  prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
454  if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
455 
456  if( dbg && (h > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance())
457  && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+eps) ) ) )
458  {
459  static G4ThreadLocal G4int noWarnings = 0;
460  std::ostringstream message;
461  if( (noWarnings++ < 10) || (dbg>2) )
462  {
463  message << "The integration produced an end-point which " << G4endl
464  << "is further from the start-point than the curve length."
465  << G4endl;
466  }
467  message << " Distance of endpoints = " << endPointDist
468  << ", curve length = " << h << G4endl
469  << " Difference (curveLen-endpDist)= " << (h - endPointDist)
470  << ", relative = " << (h-endPointDist) / h
471  << ", epsilon = " << eps;
472  G4Exception("G4MagInt_Driver::WarnEndPointTooFar()", "GeomField1001",
473  JustWarning, message);
474  }
475 }
476 
477 // ---------------------------------------------------------
478 
479 void
480 G4MagInt_Driver::OneGoodStep( G4double y[], // InOut
481  const G4double dydx[],
482  G4double& x, // InOut
483  G4double htry,
484  G4double eps_rel_max,
485  G4double& hdid, // Out
486  G4double& hnext ) // Out
487 
488 // Driver for one Runge-Kutta Step with monitoring of local truncation error
489 // to ensure accuracy and adjust stepsize. Input are dependent variable
490 // array y[0,...,5] and its derivative dydx[0,...,5] at the
491 // starting value of the independent variable x . Also input are stepsize
492 // to be attempted htry, and the required accuracy eps. On output y and x
493 // are replaced by their new values, hdid is the stepsize that was actually
494 // accomplished, and hnext is the estimated next stepsize.
495 // This is similar to the function rkqs from the book:
496 // Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
497 // Edition, by William H. Press, Saul A. Teukolsky, William T.
498 // Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
499 // 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
500 
501 {
502  G4double errmax_sq;
503  G4double h, htemp, xnew ;
504 
505  G4double yerr[G4FieldTrack::ncompSVEC], ytemp[G4FieldTrack::ncompSVEC];
506 
507  h = htry ; // Set stepsize to the initial trial value
508 
509  G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
510 
511  G4double errpos_sq = 0.0; // square of displacement error
512  G4double errvel_sq = 0.0; // square of momentum vector difference
513  G4double errspin_sq = 0.0; // square of spin vector difference
514 
515  static G4ThreadLocal G4int tot_no_trials=0;
516  const G4int max_trials=100;
517 
518  G4ThreeVector Spin(y[9],y[10],y[11]);
519  G4double spin_mag2 = Spin.mag2();
520  G4bool hasSpin = (spin_mag2 > 0.0);
521 
522  for (G4int iter=0; iter<max_trials; ++iter)
523  {
524  ++tot_no_trials;
525  pIntStepper-> Stepper(y,dydx,h,ytemp,yerr);
526  // *******
527  G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep);
528  G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos);
529 
530  // Evaluate accuracy
531  //
532  errpos_sq = sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
533  errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
534 
535  // Accuracy for momentum
536  G4double magvel_sq= sqr(y[3]) + sqr(y[4]) + sqr(y[5]) ;
537  G4double sumerr_sq = sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) ;
538  if( magvel_sq > 0.0 )
539  {
540  errvel_sq = sumerr_sq / magvel_sq;
541  }
542  else
543  {
544  std::ostringstream message;
545  message << "Found case of zero momentum." << G4endl
546  << "- iteration= " << iter << "; h= " << h;
547  G4Exception("G4MagInt_Driver::OneGoodStep()",
548  "GeomField1001", JustWarning, message);
549  errvel_sq = sumerr_sq;
550  }
551  errvel_sq *= inv_eps_vel_sq;
552  errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
553 
554  if( hasSpin )
555  {
556  // Accuracy for spin
557  errspin_sq = ( sqr(yerr[9]) + sqr(yerr[10]) + sqr(yerr[11]) )
558  / spin_mag2; // ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
559  errspin_sq *= inv_eps_vel_sq;
560  errmax_sq = std::max( errmax_sq, errspin_sq );
561  }
562 
563  if ( errmax_sq <= 1.0 ) { break; } // Step succeeded.
564 
565  // Step failed; compute the size of retrial Step.
566  htemp = GetSafety() * h * std::pow( errmax_sq, 0.5*GetPshrnk() );
567 
568  if (htemp >= 0.1*h) { h = htemp; } // Truncation error too large,
569  else { h = 0.1*h; } // reduce stepsize, but no more
570  // than a factor of 10
571  xnew = x + h;
572  if(xnew == x)
573  {
574  std::ostringstream message;
575  message << "Stepsize underflow in Stepper !" << G4endl
576  << "- Step's start x=" << x << " and end x= " << xnew
577  << " are equal !! " << G4endl
578  << " Due to step-size= " << h
579  << ". Note that input step was " << htry;
580  G4Exception("G4MagInt_Driver::OneGoodStep()",
581  "GeomField1001", JustWarning, message);
582  break;
583  }
584  }
585 
586  // Compute size of next Step
587  if (errmax_sq > errcon*errcon)
588  {
589  hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow());
590  }
591  else
592  {
593  hnext = max_stepping_increase*h ; // No more than a factor of 5 increase
594  }
595  x += (hdid = h);
596 
597  for(G4int k=0; k<fNoIntegrationVariables; ++k) { y[k] = ytemp[k]; }
598 
599  return;
600 }
601 
602 //----------------------------------------------------------------------
603 
604 // QuickAdvance just tries one Step - it does not ensure accuracy
605 //
606 G4bool G4MagInt_Driver::QuickAdvance(G4FieldTrack& y_posvel, // INOUT
607  const G4double dydx[],
608  G4double hstep, // In
609  G4double& dchord_step,
610  G4double& dyerr_pos_sq,
611  G4double& dyerr_mom_rel_sq )
612 {
613  G4Exception("G4MagInt_Driver::QuickAdvance()", "GeomField0001",
614  FatalException, "Not yet implemented.");
615 
616  // Use the parameters of this method, to please compiler
617  //
618  dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0];
619  dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
620  return true;
621 }
622 
623 //----------------------------------------------------------------------
624 
625 G4bool G4MagInt_Driver::QuickAdvance(G4FieldTrack& y_posvel, // INOUT
626  const G4double dydx[],
627  G4double hstep, // In
628  G4double& dchord_step,
629  G4double& dyerr )
630 {
631  G4double dyerr_pos_sq, dyerr_mom_rel_sq;
632  G4double yerr_vec[G4FieldTrack::ncompSVEC],
633  yarrin[G4FieldTrack::ncompSVEC], yarrout[G4FieldTrack::ncompSVEC];
634  G4double s_start;
635  G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
636 
637  static G4ThreadLocal G4int no_call = 0;
638  ++no_call;
639 
640  // Move data into array
641  y_posvel.DumpToArray( yarrin ); // yarrin <== y_posvel
642  s_start = y_posvel.GetCurveLength();
643 
644  // Do an Integration Step
645  pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ;
646 
647  // Estimate curve-chord distance
648  dchord_step= pIntStepper-> DistChord();
649 
650  // Put back the values. yarrout ==> y_posvel
651  y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );
652  y_posvel.SetCurveLength( s_start + hstep );
653 
654 #ifdef G4DEBUG_FIELD
655  if(fVerboseLevel>2)
656  {
657  G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
658  PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep, 1);
659  }
660 #endif
661 
662  // A single measure of the error
663  // TO-DO : account for energy, spin, ... ?
664  vel_mag_sq = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
665  inv_vel_mag_sq = 1.0 / vel_mag_sq;
666  dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
667  dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
668  dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
669 
670  // Calculate also the change in the momentum squared also ???
671  // G4double veloc_square = y_posvel.GetVelocity().mag2();
672  // ...
673 
674 #ifdef RETURN_A_NEW_STEP_LENGTH
675  // The following step cannot be done here because "eps" is not known.
676  dyerr_len = std::sqrt( dyerr_len_sq );
677  dyerr_len_sq /= eps ;
678 
679  // Look at the velocity deviation ?
680  // sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
681 
682  // Set suggested new step
683  hstep = ComputeNewStepSize( dyerr_len, hstep);
684 #endif
685 
686  if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) )
687  {
688  dyerr = std::sqrt(dyerr_pos_sq);
689  }
690  else
691  {
692  // Scale it to the current step size - for now
693  dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
694  }
695 
696  return true;
697 }
698 
699 // --------------------------------------------------------------------------
700 
701 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT
702 G4bool G4MagInt_Driver::QuickAdvance(G4double yarrin[], // In
703  const G4double dydx[],
704  G4double hstep, // In
705  G4double yarrout[],
706  G4double& dchord_step,
707  G4double& dyerr ) // In length
708 {
709  G4Exception("G4MagInt_Driver::QuickAdvance()", "GeomField0001",
710  FatalException, "Not yet implemented.");
711  dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
712  yarrout[0]= yarrin[0];
713 }
714 #endif
715 
716 // --------------------------------------------------------------------------
717 
718 // This method computes new step sizes - but does not limit changes to
719 // within certain factors
720 //
721 G4double G4MagInt_Driver::
722 ComputeNewStepSize(G4double errMaxNorm, // max error (normalised)
723  G4double hstepCurrent) // current step size
724 {
725  G4double hnew;
726 
727  // Compute size of next Step for a failed step
728  if(errMaxNorm > 1.0 )
729  {
730  // Step failed; compute the size of retrial Step.
731  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
732  }
733  else if(errMaxNorm > 0.0 )
734  {
735  // Compute size of next Step for a successful step
736  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
737  }
738  else
739  {
740  // if error estimate is zero (possible) or negative (dubious)
741  hnew = max_stepping_increase * hstepCurrent;
742  }
743 
744  return hnew;
745 }
746 
747 // ---------------------------------------------------------------------------
748 
749 // This method computes new step sizes limiting changes within certain factors
750 //
751 // It shares its logic with AccurateAdvance.
752 // They are kept separate currently for optimisation.
753 //
754 G4double
755 G4MagInt_Driver::ComputeNewStepSize_WithinLimits(
756  G4double errMaxNorm, // max error (normalised)
757  G4double hstepCurrent) // current step size
758 {
759  G4double hnew;
760 
761  // Compute size of next Step for a failed step
762  if (errMaxNorm > 1.0 )
763  {
764  // Step failed; compute the size of retrial Step.
765  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
766 
767  if (hnew < max_stepping_decrease*hstepCurrent)
768  {
769  hnew = max_stepping_decrease*hstepCurrent ;
770  // reduce stepsize, but no more
771  // than this factor (value= 1/10)
772  }
773  }
774  else
775  {
776  // Compute size of next Step for a successful step
777  if (errMaxNorm > errcon)
778  { hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()); }
779  else // No more than a factor of 5 increase
780  { hnew = max_stepping_increase * hstepCurrent; }
781  }
782  return hnew;
783 }
784 
785 // ---------------------------------------------------------------------------
786 
787 void G4MagInt_Driver::PrintStatus( const G4double* StartArr,
788  G4double xstart,
789  const G4double* CurrentArr,
790  G4double xcurrent,
791  G4double requestStep,
792  G4int subStepNo )
793  // Potentially add as arguments:
794  // <dydx> - as Initial Force
795  // stepTaken(hdid) - last step taken
796  // nextStep (hnext) - proposal for size
797 {
798  G4FieldTrack StartFT(G4ThreeVector(0,0,0),
799  G4ThreeVector(0,0,0), 0., 0., 0., 0. );
800  G4FieldTrack CurrentFT (StartFT);
801 
802  StartFT.LoadFromArray( StartArr, fNoIntegrationVariables);
803  StartFT.SetCurveLength( xstart);
804  CurrentFT.LoadFromArray( CurrentArr, fNoIntegrationVariables);
805  CurrentFT.SetCurveLength( xcurrent );
806 
807  PrintStatus(StartFT, CurrentFT, requestStep, subStepNo );
808 }
809 
810 // ---------------------------------------------------------------------------
811 
812 void G4MagInt_Driver::PrintStatus(const G4FieldTrack& StartFT,
813  const G4FieldTrack& CurrentFT,
814  G4double requestStep,
815  G4int subStepNo)
816 {
817  G4int verboseLevel= fVerboseLevel;
818  const G4int noPrecision = 5;
819  G4int oldPrec= G4cout.precision(noPrecision);
820  // G4cout.setf(ios_base::fixed,ios_base::floatfield);
821 
822  const G4ThreeVector StartPosition= StartFT.GetPosition();
823  const G4ThreeVector StartUnitVelocity= StartFT.GetMomentumDir();
824  const G4ThreeVector CurrentPosition= CurrentFT.GetPosition();
825  const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir();
826 
827  G4double DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
828 
829  G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
830  G4double subStepSize = step_len;
831 
832  if( (subStepNo <= 1) || (verboseLevel > 3) )
833  {
834  subStepNo = - subStepNo; // To allow printing banner
835 
836  G4cout << std::setw( 6) << " " << std::setw( 25)
837  << " G4MagInt_Driver: Current Position and Direction" << " "
838  << G4endl;
839  G4cout << std::setw( 5) << "Step#" << " "
840  << std::setw( 7) << "s-curve" << " "
841  << std::setw( 9) << "X(mm)" << " "
842  << std::setw( 9) << "Y(mm)" << " "
843  << std::setw( 9) << "Z(mm)" << " "
844  << std::setw( 8) << " N_x " << " "
845  << std::setw( 8) << " N_y " << " "
846  << std::setw( 8) << " N_z " << " "
847  << std::setw( 8) << " N^2-1 " << " "
848  << std::setw(10) << " N(0).N " << " "
849  << std::setw( 7) << "KinEner " << " "
850  << std::setw(12) << "Track-l" << " " // Add the Sub-step ??
851  << std::setw(12) << "Step-len" << " "
852  << std::setw(12) << "Step-len" << " "
853  << std::setw( 9) << "ReqStep" << " "
854  << G4endl;
855  }
856 
857  if( (subStepNo <= 0) )
858  {
859  PrintStat_Aux( StartFT, requestStep, 0.,
860  0, 0.0, 1.0);
861  }
862 
863  if( verboseLevel <= 3 )
864  {
865  G4cout.precision(noPrecision);
866  PrintStat_Aux( CurrentFT, requestStep, step_len,
867  subStepNo, subStepSize, DotStartCurrentVeloc );
868  }
869 
870  G4cout.precision(oldPrec);
871 }
872 
873 // ---------------------------------------------------------------------------
874 
875 void G4MagInt_Driver::PrintStat_Aux(const G4FieldTrack& aFieldTrack,
876  G4double requestStep,
877  G4double step_len,
878  G4int subStepNo,
879  G4double subStepSize,
880  G4double dotVeloc_StartCurr)
881 {
882  const G4ThreeVector Position = aFieldTrack.GetPosition();
883  const G4ThreeVector UnitVelocity = aFieldTrack.GetMomentumDir();
884 
885  if( subStepNo >= 0)
886  {
887  G4cout << std::setw( 5) << subStepNo << " ";
888  }
889  else
890  {
891  G4cout << std::setw( 5) << "Start" << " ";
892  }
893  G4double curveLen= aFieldTrack.GetCurveLength();
894  G4cout << std::setw( 7) << curveLen;
895  G4cout << std::setw( 9) << Position.x() << " "
896  << std::setw( 9) << Position.y() << " "
897  << std::setw( 9) << Position.z() << " "
898  << std::setw( 8) << UnitVelocity.x() << " "
899  << std::setw( 8) << UnitVelocity.y() << " "
900  << std::setw( 8) << UnitVelocity.z() << " ";
901  G4int oldprec= G4cout.precision(3);
902  G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 << " ";
903  G4cout.precision(6);
904  G4cout << std::setw(10) << dotVeloc_StartCurr << " ";
905  G4cout.precision(oldprec);
906  G4cout << std::setw( 7) << aFieldTrack.GetKineticEnergy();
907  G4cout << std::setw(12) << step_len << " ";
908 
909  static G4ThreadLocal G4double oldCurveLength = 0.0;
910  static G4ThreadLocal G4double oldSubStepLength = 0.0;
911  static G4ThreadLocal G4int oldSubStepNo = -1;
912 
913  G4double subStep_len = 0.0;
914  if( curveLen > oldCurveLength )
915  {
916  subStep_len= curveLen - oldCurveLength;
917  }
918  else if (subStepNo == oldSubStepNo)
919  {
920  subStep_len= oldSubStepLength;
921  }
922  oldCurveLength= curveLen;
923  oldSubStepLength= subStep_len;
924 
925  G4cout << std::setw(12) << subStep_len << " ";
926  G4cout << std::setw(12) << subStepSize << " ";
927  if( requestStep != -1.0 )
928  {
929  G4cout << std::setw( 9) << requestStep << " ";
930  }
931  else
932  {
933  G4cout << std::setw( 9) << " InitialStep " << " ";
934  }
935  G4cout << G4endl;
936 }
937 
938 // ---------------------------------------------------------------------------
939 
940 void G4MagInt_Driver::PrintStatisticsReport()
941 {
942  G4int noPrecBig = 6;
943  G4int oldPrec = G4cout.precision(noPrecBig);
944 
945  G4cout << "G4MagInt_Driver Statistics of steps undertaken. " << G4endl;
946  G4cout << "G4MagInt_Driver: Number of Steps: "
947  << " Total= " << fNoTotalSteps
948  << " Bad= " << fNoBadSteps
949  << " Small= " << fNoSmallSteps
950  << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
951  << G4endl;
952  G4cout.precision(oldPrec);
953 }
954 
955 // ---------------------------------------------------------------------------
956 
957 void G4MagInt_Driver::SetSmallestFraction(G4double newFraction)
958 {
959  if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
960  {
961  fSmallestFraction= newFraction;
962  }
963  else
964  {
965  std::ostringstream message;
966  message << "Smallest Fraction not changed. " << G4endl
967  << " Proposed value was " << newFraction << G4endl
968  << " Value must be between 1.e-8 and 1.e-16";
969  G4Exception("G4MagInt_Driver::SetSmallestFraction()",
970  "GeomField1001", JustWarning, message);
971  }
972 }
973 
974 void G4MagInt_Driver::
975 GetDerivatives(const G4FieldTrack& y_curr, G4double* dydx) const
976 {
977  G4double ytemp[G4FieldTrack::ncompSVEC];
978  y_curr.DumpToArray(ytemp);
979  pIntStepper->RightHandSide(ytemp, dydx);
980  // Avoid virtual call for GetStepper
981  // Was: GetStepper()->ComputeRightHandSide(ytemp, dydx);
982 }
983 
984 void G4MagInt_Driver::GetDerivatives(const G4FieldTrack& track,
985  G4double dydx[],
986  G4double field[]) const
987 {
988  G4double ytemp[G4FieldTrack::ncompSVEC];
989  track.DumpToArray(ytemp);
990  pIntStepper->RightHandSide(ytemp, dydx, field);
991 }
992 
993 G4EquationOfMotion* G4MagInt_Driver::GetEquationOfMotion()
994 {
995  return pIntStepper->GetEquationOfMotion();
996 }
997 
998 void G4MagInt_Driver::SetEquationOfMotion(G4EquationOfMotion *equation)
999 {
1000  pIntStepper->SetEquationOfMotion(equation);
1001 }
1002 
1003 const G4MagIntegratorStepper* G4MagInt_Driver::GetStepper() const
1004 {
1005  return pIntStepper;
1006 }
1007 
1008 G4MagIntegratorStepper* G4MagInt_Driver::GetStepper()
1009 {
1010  return pIntStepper;
1011 }
1012 
1013 void G4MagInt_Driver::
1014 RenewStepperAndAdjust(G4MagIntegratorStepper* pItsStepper)
1015 {
1016  pIntStepper = pItsStepper;
1017  ReSetParameters();
1018 }