Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4OldMagIntDriver.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // G4OldMagIntDriver -- same behaviour as old G4MagInt_Driver
 27 //
 28 // V.Grichine, 07.10.1996 - Created
 29 // J.Apostolakis, 08.11.2001 - Respect minimum step in AccurateAdvance
 30 // --------------------------------------------------------------------
 31 
 32 #include <iomanip>
 33 
 34 #include "globals.hh"
 35 #include "G4SystemOfUnits.hh"
 36 #include "G4GeometryTolerance.hh"
 37 #include "G4OldMagIntDriver.hh"
 38 #include "G4FieldTrack.hh"
 39 
 40 #ifdef   G4DEBUG_FIELD
 41 #include "G4DriverReporter.hh"
 42 #endif
 43 
 44 // ---------------------------------------------------------
 45 
 46 //  Constructor
 47 //
 48 G4OldMagIntDriver::G4OldMagIntDriver( G4double                hminimum, 
 49                                   G4MagIntegratorStepper* pStepper,
 50                                   G4int                   numComponents,
 51                                   G4int                   statisticsVerbose)
 52   : fNoIntegrationVariables(numComponents), 
 53     fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )),
 54     fStatisticsVerboseLevel(statisticsVerbose)
 55 {  
 56   // In order to accomodate "Laboratory Time", which is [7], fMinNoVars=8
 57   // is required. For proper time of flight and spin,  fMinNoVars must be 12
 58 
 59   RenewStepperAndAdjust( pStepper );
 60   fMinimumStep = hminimum;
 61 
 62   fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
 63 #ifdef G4DEBUG_FIELD
 64   fVerboseLevel=2;
 65 #endif
 66 
 67   if( (fVerboseLevel > 0) || (fStatisticsVerboseLevel > 1) )
 68   {
 69     G4cout << "MagIntDriver version: Accur-Adv: "
 70            << "invE_nS, QuickAdv-2sqrt with Statistics "
 71 #ifdef G4FLD_STATS
 72            << " enabled "
 73 #else
 74            << " disabled "
 75 #endif
 76            << G4endl;
 77   }
 78 }
 79 
 80 // ---------------------------------------------------------
 81 
 82 //  Destructor
 83 //
 84 G4OldMagIntDriver::~G4OldMagIntDriver()
 85 { 
 86   if( fStatisticsVerboseLevel > 1 )
 87   {
 88     PrintStatisticsReport();
 89   }
 90 }
 91 
 92 // ---------------------------------------------------------
 93 
 94 G4bool
 95 G4OldMagIntDriver::AccurateAdvance(G4FieldTrack& y_current,
 96                                  G4double      hstep,
 97                                  G4double      eps,
 98                                  G4double      hinitial )
 99 {
100   // Runge-Kutta driver with adaptive stepsize control. Integrate starting
101   // values at y_current over hstep x2 with accuracy eps. 
102   // On output ystart is replaced by values at the end of the integration 
103   // interval. RightHandSide is the right-hand side of ODE system. 
104   // The source is similar to odeint routine from NRC p.721-722 .
105 
106   G4int nstp, i;
107   G4double x, hnext, hdid, h;
108 
109 #ifdef G4DEBUG_FIELD
110   G4int no_warnings = 0;
111   static G4int dbg = 1;
112   G4double ySubStepStart[G4FieldTrack::ncompSVEC];
113   G4FieldTrack  yFldTrkStart(y_current);
114 #endif
115 
116   G4double y[G4FieldTrack::ncompSVEC], dydx[G4FieldTrack::ncompSVEC];
117   G4double ystart[G4FieldTrack::ncompSVEC], yEnd[G4FieldTrack::ncompSVEC]; 
118   G4double  x1, x2;
119   G4bool succeeded = true;
120 
121   G4double startCurveLength;
122 
123   const G4int nvar = fNoVars;
124 
125   G4FieldTrack yStartFT(y_current);
126 
127   //  Ensure that hstep > 0
128   //
129   if( hstep <= 0.0 )
130   { 
131     if( hstep == 0.0 )
132     {
133       std::ostringstream message;
134       message << "Proposed step is zero; hstep = " << hstep << " !";
135       G4Exception("G4OldMagIntDriver::AccurateAdvance()", 
136                   "GeomField1001", JustWarning, message);
137       return succeeded; 
138     }
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("G4OldMagIntDriver::AccurateAdvance()", 
145                 "GeomField0003", EventMustBeAborted, message);
146     return false;
147   }
148 
149   y_current.DumpToArray( ystart );
150 
151   startCurveLength= y_current.GetCurveLength();
152   x1= startCurveLength; 
153   x2= x1 + hstep;
154 
155   if ( (hinitial > 0.0) && (hinitial < hstep)
156     && (hinitial > perMillion * hstep) )
157   {
158      h = hinitial;
159   }
160   else  //  Initial Step size "h" defaults to the full interval
161   {
162      h = hstep;
163   }
164 
165   x = x1;
166 
167   for ( i=0; i<nvar; ++i)  { y[i] = ystart[i]; }
168 
169 #ifdef G4DEBUG_FIELD  
170   // G4cout << "IDriver: hstep = " << hstep << " hinitial= " << hinitial << " h = " << h << G4endl;
171   G4cout << "IDriver::AccurAdv called. " 
172          << " Input: hstep = " << hstep << " hinitial= " << hinitial << " , current: h = " << h << G4endl;        
173 #endif
174   
175   G4bool lastStep= false;
176   nstp = 1;
177 
178   do
179   {
180     G4ThreeVector StartPos( y[0], y[1], y[2] );
181 
182 #ifdef G4DEBUG_FIELD
183     G4double xSubStepStart= x; 
184     for (i=0; i<nvar; ++i)  { ySubStepStart[i] = y[i]; }
185     yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
186     yFldTrkStart.SetCurveLength(x);
187     if(dbg) // Debug
188        G4cout << "----- Iteration = " << nstp-1 << G4endl;
189 #endif
190 
191     pIntStepper->RightHandSide( y, dydx );
192     ++fNoTotalSteps;
193 
194     // Perform the Integration
195     //      
196     if( h > fMinimumStep )
197     { 
198       OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ;
199       //--------------------------------------
200 
201 #ifdef G4DEBUG_FIELD
202       if (dbg) // (dbg>2)
203       {
204         G4cout << "IntegrationDriver -- after OneGoodStep / requesting step = " << h << G4endl;
205         // PrintStatus( ySubStepStart, xSubStepStart, y, x, h,  nstp); // Only        
206         G4DriverReporter::PrintStatus( ySubStepStart, xSubStepStart, y, x, h,  nstp, nvar);
207       }
208 #endif
209     }
210     else
211     {
212       G4FieldTrack yFldTrk( G4ThreeVector(0,0,0), 
213                             G4ThreeVector(0,0,0), 0., 0., 0., 0. );
214       G4double dchord_step, dyerr, dyerr_len;   // What to do with these ?
215       yFldTrk.LoadFromArray(y, fNoIntegrationVariables); 
216       yFldTrk.SetCurveLength( x );
217 
218       QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
219       //-----------------------------------------------------
220 
221       yFldTrk.DumpToArray(y);    
222 
223 #ifdef G4FLD_STATS
224       ++fNoSmallSteps; 
225       if ( dyerr_len > fDyerr_max )  { fDyerr_max = dyerr_len; }
226       fDyerrPos_smTot += dyerr_len;
227       fSumH_sm += h;  // Length total for 'small' steps
228       if (nstp<=1)  { ++fNoInitialSmallSteps; }
229 #endif
230 #ifdef G4DEBUG_FIELD
231       if (dbg>1)
232       {
233         if(fNoSmallSteps<2) { PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
234         G4cout << "Another sub-min step, no " << fNoSmallSteps 
235                << " of " << fNoTotalSteps << " this time " << nstp << G4endl; 
236         PrintStatus( ySubStepStart, x1, y, x, h,  nstp);   // Only this
237         G4cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h 
238                << " epsilon= " << eps << " hstep= " << hstep 
239                << " h= " << h << " hmin= " << fMinimumStep << G4endl;
240       }
241 #endif        
242       if( h == 0.0 )
243       { 
244         G4Exception("G4OldMagIntDriver::AccurateAdvance()",
245                     "GeomField0003", FatalException,
246                     "Integration Step became Zero!"); 
247       }
248       dyerr = dyerr_len / h;
249       hdid = h;
250       x += hdid;
251 
252       // Compute suggested new step
253       hnext = ComputeNewStepSize( dyerr/eps, h);
254 
255       // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
256     }
257 
258     G4ThreeVector EndPos( y[0], y[1], y[2] );
259 
260 #if (G4DEBUG_FIELD>1)
261     static G4int nStpPr = 50;   // For debug printing of long integrations
262     if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
263     {
264       if( nstp==nStpPr )  { G4cout << "***** Many steps ****" << G4endl; }
265       G4cout << "MagIntDrv: " ; 
266       G4cout << "hdid="  << std::setw(12) << hdid  << " "
267              << "hnext=" << std::setw(12) << hnext << " " 
268              << "hstep=" << std::setw(12) << hstep << " (requested) " 
269              << G4endl;
270       PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp); 
271     }
272 #endif
273 
274     // Check the endpoint
275     G4double endPointDist= (EndPos-StartPos).mag(); 
276     if ( endPointDist >= hdid*(1.+perMillion) )
277     {
278       ++fNoBadSteps;
279 
280       // Issue a warning only for gross differences -
281       // we understand how small difference occur.
282       if ( endPointDist >= hdid*(1.+perThousand) )
283       { 
284 #ifdef G4DEBUG_FIELD
285         if (dbg)
286         {
287           WarnEndPointTooFar ( endPointDist, hdid, eps, dbg ); 
288           G4cerr << "  Total steps:  bad " << fNoBadSteps
289                  << " current h= " << hdid << G4endl;
290           PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);  
291         }
292         ++no_warnings;
293 #endif
294       }
295     }
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     succeeded = false;
371 #ifdef G4DEBUG_FIELD
372     ++no_warnings;
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 G4OldMagIntDriver::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("G4OldMagIntDriver::WarnSmallStepSize()", "GeomField1001",
420               JustWarning, message);
421   ++noWarningsIssued;
422 }
423 
424 // ---------------------------------------------------------
425 
426 void
427 G4OldMagIntDriver::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("G4OldMagIntDriver::WarnTooManyStep()", "GeomField1001",
438                JustWarning, message);
439 }
440 
441 // ---------------------------------------------------------
442 
443 void
444 G4OldMagIntDriver::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 != 0) && (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("G4OldMagIntDriver::WarnEndPointTooFar()", "GeomField1001",
473                 JustWarning, message);
474   }
475 }
476 
477 // ---------------------------------------------------------
478 
479 void
480 G4OldMagIntDriver::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   const G4int max_trials=100; 
516 
517   G4ThreeVector Spin(y[9],y[10],y[11]);
518   G4double spin_mag2 = Spin.mag2();
519   G4bool hasSpin = (spin_mag2 > 0.0); 
520 
521   for (G4int iter=0; iter<max_trials; ++iter)
522   {
523     pIntStepper-> Stepper(y,dydx,h,ytemp,yerr); 
524     //            *******
525     G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep); 
526     G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos); 
527 
528     // Evaluate accuracy
529     //
530     errpos_sq =  sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
531     errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
532 
533     // Accuracy for momentum
534     G4double magvel_sq=  sqr(y[3]) + sqr(y[4]) + sqr(y[5]) ;
535     G4double sumerr_sq =  sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) ; 
536     if( magvel_sq > 0.0 )
537     { 
538        errvel_sq = sumerr_sq / magvel_sq; 
539     }
540     else
541     {
542        std::ostringstream message;
543        message << "Found case of zero momentum." << G4endl
544                << "- iteration= " << iter << "; h= " << h;
545        G4Exception("G4OldMagIntDriver::OneGoodStep()",
546                    "GeomField1001", JustWarning, message);
547        errvel_sq = sumerr_sq; 
548     }
549     errvel_sq *= inv_eps_vel_sq;
550     errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
551 
552     if( hasSpin )
553     { 
554       // Accuracy for spin
555       errspin_sq =  ( sqr(yerr[9]) + sqr(yerr[10]) + sqr(yerr[11]) )
556                     /  spin_mag2; // ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
557       errspin_sq *= inv_eps_vel_sq;
558       errmax_sq = std::max( errmax_sq, errspin_sq ); 
559     }
560 
561     if ( errmax_sq <= 1.0 )  { break; } // Step succeeded. 
562 
563     // Step failed; compute the size of retrial Step.
564     htemp = GetSafety() * h * std::pow( errmax_sq, 0.5*GetPshrnk() );
565 
566     if (htemp >= 0.1*h)  { h = htemp; }  // Truncation error too large,
567     else  { h = 0.1*h; }                 // reduce stepsize, but no more
568                                          // than a factor of 10
569     xnew = x + h;
570     if(xnew == x)
571     {
572       std::ostringstream message;
573       message << "Stepsize underflow in Stepper !" << G4endl
574               << "- Step's start x=" << x << " and end x= " << xnew 
575               << " are equal !! " << G4endl
576               << "  Due to step-size= " << h 
577               << ". Note that input step was " << htry;
578       G4Exception("G4OldMagIntDriver::OneGoodStep()",
579                   "GeomField1001", JustWarning, message);
580       break;
581     }
582   }
583 
584   // Compute size of next Step
585   if (errmax_sq > errcon*errcon)
586   { 
587     hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow());
588   }
589   else
590   {
591     hnext = max_stepping_increase*h ; // No more than a factor of 5 increase
592   }
593   x += (hdid = h);
594 
595   for(G4int k=0; k<fNoIntegrationVariables; ++k) { y[k] = ytemp[k]; }
596 
597   return;
598 }
599 
600 //----------------------------------------------------------------------
601 
602 // QuickAdvance just tries one Step - it does not ensure accuracy
603 //
604 G4bool G4OldMagIntDriver::QuickAdvance(G4FieldTrack& y_posvel,    // INOUT
605                                const G4double      dydx[],  
606                                      G4double      hstep,       // In
607                                      G4double&     dchord_step,
608                                      G4double&     dyerr_pos_sq,
609                                      G4double&     dyerr_mom_rel_sq )  
610 {
611   G4Exception("G4OldMagIntDriver::QuickAdvance()", "GeomField0001",
612               FatalException, "Not yet implemented."); 
613 
614   // Use the parameters of this method, to please compiler
615   //
616   dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0]; 
617   dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
618   return true;
619 }
620 
621 //----------------------------------------------------------------------
622 
623 G4bool G4OldMagIntDriver::QuickAdvance(G4FieldTrack& y_posvel,    // INOUT
624                                const G4double      dydx[],  
625                                      G4double      hstep,       // In
626                                      G4double&     dchord_step,
627                                      G4double&     dyerr )
628 {
629   G4double dyerr_pos_sq, dyerr_mom_rel_sq;  
630   G4double yerr_vec[G4FieldTrack::ncompSVEC],
631            yarrin[G4FieldTrack::ncompSVEC], yarrout[G4FieldTrack::ncompSVEC]; 
632   G4double s_start;
633   G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
634 
635 #ifdef  G4DEBUG_FIELD  
636   G4FieldTrack startTrack( y_posvel );  // For debugging
637 #endif
638 
639   // Move data into array
640   y_posvel.DumpToArray( yarrin );      //  yarrin  <== y_posvel 
641   s_start = y_posvel.GetCurveLength();
642 
643   // Do an Integration Step
644   pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ; 
645 
646   // Estimate curve-chord distance
647   dchord_step= pIntStepper-> DistChord();
648 
649   // Put back the values.  yarrout ==> y_posvel
650   y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );
651   y_posvel.SetCurveLength( s_start + hstep );
652 
653 #ifdef  G4DEBUG_FIELD
654   if(fVerboseLevel>2)
655   {
656     G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
657     PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep,  1); 
658   }
659 #endif
660 
661   // A single measure of the error   
662   //      TO-DO :  account for  energy,  spin, ... ? 
663   vel_mag_sq   = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
664   inv_vel_mag_sq = 1.0 / vel_mag_sq; 
665   dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
666   dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
667   dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
668 
669   // Calculate also the change in the momentum squared also ???
670   // G4double veloc_square = y_posvel.GetVelocity().mag2();
671   // ...
672 
673 #ifdef RETURN_A_NEW_STEP_LENGTH
674   // The following step cannot be done here because "eps" is not known.
675   dyerr_len = std::sqrt( dyerr_len_sq ); 
676   dyerr_len_sq /= eps ;
677 
678   // Look at the velocity deviation ?
679   //  sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
680 
681   // Set suggested new step
682   hstep = ComputeNewStepSize( dyerr_len, hstep);
683 #endif
684 
685   if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) )
686   {
687     dyerr = std::sqrt(dyerr_pos_sq);
688   }
689   else
690   {
691     // Scale it to the current step size - for now
692     dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
693   }
694 
695 #ifdef  G4DEBUG_FIELD  
696   // For debugging
697   G4cout // << "G4MagInt_Driver::"
698          << "QuickAdvance" << G4endl
699          << " Input:  hstep= " << hstep       << G4endl
700          << "         track= " << startTrack  << G4endl
701          << " Output: track= " << y_posvel    << G4endl
702          << "         d_chord = " << dchord_step
703          << " dyerr = " << dyerr << G4endl;
704 #endif
705 
706   return true;
707 }
708 
709 // --------------------------------------------------------------------------
710 
711 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT
712 G4bool  G4OldMagIntDriver::QuickAdvance(G4double  yarrin[],    // In
713                                 const G4double  dydx[],  
714                                       G4double  hstep,       // In
715                                       G4double  yarrout[],
716                                       G4double& dchord_step,
717                                       G4double& dyerr )      // In length
718 {
719   G4Exception("G4OldMagIntDriver::QuickAdvance()", "GeomField0001",
720               FatalException, "Not yet implemented.");
721   dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
722   yarrout[0]= yarrin[0];
723 }
724 #endif 
725 
726 // --------------------------------------------------------------------------
727 
728 // This method computes new step sizes - but does not limit changes to
729 // within  certain factors
730 // 
731 G4double G4OldMagIntDriver::
732 ComputeNewStepSize(G4double  errMaxNorm,    // max error  (normalised)
733                    G4double  hstepCurrent)  // current step size
734 {
735   G4double hnew;
736 
737   // Compute size of next Step for a failed step
738   if(errMaxNorm > 1.0 )
739   {
740     // Step failed; compute the size of retrial Step.
741     hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
742   }
743   else if(errMaxNorm > 0.0 )
744   {
745     // Compute size of next Step for a successful step
746     hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
747   }
748   else
749   {
750     // if error estimate is zero (possible) or negative (dubious)
751     hnew = max_stepping_increase * hstepCurrent; 
752   }
753 
754   return hnew;
755 }
756 
757 // ---------------------------------------------------------------------------
758 
759 // This method computes new step sizes limiting changes within certain factors
760 // 
761 // It shares its logic with AccurateAdvance.
762 // They are kept separate currently for optimisation.
763 //
764 G4double 
765 G4OldMagIntDriver::ComputeNewStepSize_WithinLimits( 
766                           G4double  errMaxNorm,    // max error  (normalised)
767                           G4double  hstepCurrent)  // current step size
768 {
769   G4double hnew;
770 
771   // Compute size of next Step for a failed step
772   if (errMaxNorm > 1.0 )
773   {
774     // Step failed; compute the size of retrial Step.
775     hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
776   
777     if (hnew < max_stepping_decrease*hstepCurrent)
778     {
779       hnew = max_stepping_decrease*hstepCurrent ;
780                          // reduce stepsize, but no more
781                          // than this factor (value= 1/10)
782     }
783   }
784   else
785   {
786     // Compute size of next Step for a successful step
787     if (errMaxNorm > errcon)
788      { hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()); }
789     else  // No more than a factor of 5 increase
790      { hnew = max_stepping_increase * hstepCurrent; }
791   }
792   return hnew;
793 }
794 
795 // ---------------------------------------------------------------------------
796 
797 void G4OldMagIntDriver::PrintStatus( const G4double* StartArr,  
798                                          G4double  xstart,
799                                    const G4double* CurrentArr, 
800                                          G4double  xcurrent,
801                                          G4double  requestStep, 
802                                          G4int     subStepNo )
803   // Potentially add as arguments:  
804   //                                 <dydx>           - as Initial Force
805   //                                 stepTaken(hdid)  - last step taken
806   //                                 nextStep (hnext) - proposal for size
807 {
808    G4FieldTrack  StartFT(G4ThreeVector(0,0,0),
809                  G4ThreeVector(0,0,0), 0., 0., 0., 0. );
810    G4FieldTrack  CurrentFT (StartFT);
811 
812    StartFT.LoadFromArray( StartArr, fNoIntegrationVariables); 
813    StartFT.SetCurveLength( xstart);
814    CurrentFT.LoadFromArray( CurrentArr, fNoIntegrationVariables); 
815    CurrentFT.SetCurveLength( xcurrent );
816 
817    PrintStatus(StartFT, CurrentFT, requestStep, subStepNo ); 
818 }
819 
820 // ---------------------------------------------------------------------------
821 
822 void G4OldMagIntDriver::PrintStatus(const G4FieldTrack& StartFT,
823                                   const G4FieldTrack& CurrentFT, 
824                                         G4double      requestStep, 
825                                         G4int         subStepNo)
826 {
827     G4int verboseLevel= fVerboseLevel;
828     const G4int noPrecision = 5;
829     G4long oldPrec= G4cout.precision(noPrecision);
830     // G4cout.setf(ios_base::fixed,ios_base::floatfield);
831 
832     const G4ThreeVector StartPosition=       StartFT.GetPosition();
833     const G4ThreeVector StartUnitVelocity=   StartFT.GetMomentumDir();
834     const G4ThreeVector CurrentPosition=     CurrentFT.GetPosition();
835     const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir();
836 
837     G4double  DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
838 
839     G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
840     G4double subStepSize = step_len;
841      
842     if( (subStepNo <= 1) || (verboseLevel > 3) )
843     {
844        subStepNo = - subStepNo;        // To allow printing banner
845 
846        G4cout << std::setw( 6)  << " " << std::setw( 25)
847               << " G4OldMagIntDriver: Current Position  and  Direction" << " "
848               << G4endl; 
849        G4cout << std::setw( 5) << "Step#" << " "
850               << std::setw( 7) << "s-curve" << " "
851               << std::setw( 9) << "X(mm)" << " "
852               << std::setw( 9) << "Y(mm)" << " "  
853               << std::setw( 9) << "Z(mm)" << " "
854               << std::setw( 8) << " N_x " << " "
855               << std::setw( 8) << " N_y " << " "
856               << std::setw( 8) << " N_z " << " "
857               << std::setw( 8) << " N^2-1 " << " "
858               << std::setw(10) << " N(0).N " << " "
859               << std::setw( 7) << "KinEner " << " "
860               << std::setw(12) << "Track-l" << " "   // Add the Sub-step ??
861               << std::setw(12) << "Step-len" << " " 
862               << std::setw(12) << "Step-len" << " " 
863               << std::setw( 9) << "ReqStep" << " "  
864               << G4endl;
865     }
866 
867     if( (subStepNo <= 0) )
868     {
869       PrintStat_Aux( StartFT,  requestStep, 0., 
870                        0,        0.0,         1.0);
871     }
872 
873     if( verboseLevel <= 3 )
874     {
875       G4cout.precision(noPrecision);
876       PrintStat_Aux( CurrentFT, requestStep, step_len, 
877                      subStepNo, subStepSize, DotStartCurrentVeloc );
878     }
879 
880     G4cout.precision(oldPrec);
881 }
882 
883 // ---------------------------------------------------------------------------
884 
885 void G4OldMagIntDriver::PrintStat_Aux(const G4FieldTrack& aFieldTrack,
886                                           G4double      requestStep, 
887                                           G4double      step_len,
888                                           G4int         subStepNo,
889                                           G4double      subStepSize,
890                                           G4double      dotVeloc_StartCurr)
891 {
892     const G4ThreeVector Position = aFieldTrack.GetPosition();
893     const G4ThreeVector UnitVelocity = aFieldTrack.GetMomentumDir();
894  
895     if( subStepNo >= 0)
896     {
897        G4cout << std::setw( 5) << subStepNo << " ";
898     }
899     else
900     {
901        G4cout << std::setw( 5) << "Start" << " ";
902     }
903     G4double curveLen= aFieldTrack.GetCurveLength();
904     G4cout << std::setw( 7) << curveLen;
905     G4cout << std::setw( 9) << Position.x() << " "
906            << std::setw( 9) << Position.y() << " "
907            << std::setw( 9) << Position.z() << " "
908            << std::setw( 8) << UnitVelocity.x() << " "
909            << std::setw( 8) << UnitVelocity.y() << " "
910            << std::setw( 8) << UnitVelocity.z() << " ";
911     G4long oldprec= G4cout.precision(3);
912     G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 << " ";
913     G4cout.precision(6);
914     G4cout << std::setw(10) << dotVeloc_StartCurr << " ";
915     G4cout.precision(oldprec);
916     G4cout << std::setw( 7) << aFieldTrack.GetKineticEnergy();
917     G4cout << std::setw(12) << step_len << " ";
918 
919     static G4ThreadLocal G4double oldCurveLength = 0.0;
920     static G4ThreadLocal G4double oldSubStepLength = 0.0;
921     static G4ThreadLocal G4int oldSubStepNo = -1;
922 
923     G4double subStep_len = 0.0;
924     if( curveLen > oldCurveLength )
925     {
926       subStep_len= curveLen - oldCurveLength;
927     }
928     else if (subStepNo == oldSubStepNo)
929     {
930       subStep_len= oldSubStepLength;
931     }
932     oldCurveLength= curveLen;
933     oldSubStepLength= subStep_len;
934 
935     G4cout << std::setw(12) << subStep_len << " "; 
936     G4cout << std::setw(12) << subStepSize << " "; 
937     if( requestStep != -1.0 )
938     {
939       G4cout << std::setw( 9) << requestStep << " ";
940     }
941     else
942     {
943        G4cout << std::setw( 9) << " InitialStep " << " ";
944     }
945     G4cout << G4endl;
946 }
947 
948 // ---------------------------------------------------------------------------
949 
950 void G4OldMagIntDriver::PrintStatisticsReport()
951 {
952   G4int noPrecBig = 6;
953   G4long oldPrec = G4cout.precision(noPrecBig);
954 
955   G4cout << "G4OldMagIntDriver Statistics of steps undertaken. " << G4endl;
956   G4cout << "G4OldMagIntDriver: Number of Steps: "
957          << " Total= " <<  fNoTotalSteps
958          << " Bad= "   <<  fNoBadSteps 
959          << " Small= " <<  fNoSmallSteps 
960          << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
961          << G4endl;
962  G4cout.precision(oldPrec);
963 }
964  
965 // ---------------------------------------------------------------------------
966 
967 void G4OldMagIntDriver::SetSmallestFraction(G4double newFraction)
968 {
969   if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
970   {
971     fSmallestFraction= newFraction;
972   }
973   else
974   { 
975     std::ostringstream message;
976     message << "Smallest Fraction not changed. " << G4endl
977             << "  Proposed value was " << newFraction << G4endl
978             << "  Value must be between 1.e-8 and 1.e-16";
979     G4Exception("G4OldMagIntDriver::SetSmallestFraction()",
980                 "GeomField1001", JustWarning, message);
981   }
982 }
983 
984 void G4OldMagIntDriver::
985 GetDerivatives(const G4FieldTrack& y_curr, G4double* dydx) const
986 {
987     G4double ytemp[G4FieldTrack::ncompSVEC];
988     y_curr.DumpToArray(ytemp);
989     pIntStepper->RightHandSide(ytemp, dydx);
990       // Avoid virtual call for GetStepper
991       // Was: GetStepper()->ComputeRightHandSide(ytemp, dydx);
992 }
993 
994 void G4OldMagIntDriver::GetDerivatives(const G4FieldTrack& track,
995                                      G4double dydx[],
996                                      G4double field[]) const
997 {
998     G4double ytemp[G4FieldTrack::ncompSVEC];
999     track.DumpToArray(ytemp);
1000     pIntStepper->RightHandSide(ytemp, dydx, field);
1001 }
1002 
1003 G4EquationOfMotion* G4OldMagIntDriver::GetEquationOfMotion()
1004 {
1005     return pIntStepper->GetEquationOfMotion();
1006 }
1007 
1008 void G4OldMagIntDriver::SetEquationOfMotion(G4EquationOfMotion *equation)
1009 {
1010     pIntStepper->SetEquationOfMotion(equation);
1011 }
1012 
1013 const G4MagIntegratorStepper* G4OldMagIntDriver::GetStepper() const
1014 {
1015     return pIntStepper;
1016 }
1017 
1018 G4MagIntegratorStepper* G4OldMagIntDriver::GetStepper()
1019 {
1020     return pIntStepper;
1021 }
1022 
1023 void G4OldMagIntDriver::
1024 RenewStepperAndAdjust(G4MagIntegratorStepper* pItsStepper)
1025 {  
1026     pIntStepper = pItsStepper; 
1027     ReSetParameters();
1028 }
1029 
1030 void G4OldMagIntDriver::StreamInfo( std::ostream& os ) const
1031 {
1032     os << "State of G4OldMagIntDriver: " << std::endl;
1033     os << "  Max number of Steps = " << fMaxNoSteps 
1034        << "    (base # = " << fMaxStepBase << " )" << std::endl;
1035     os << "  Safety factor       = " << safety  << std::endl;
1036     os << "  Power - shrink      = " << pshrnk << std::endl;
1037     os << "  Power - grow        = " << pgrow << std::endl;
1038     os << "  threshold (errcon)  = " << errcon << std::endl;
1039 
1040     os << "    fMinimumStep =      " << fMinimumStep  << std::endl;
1041     os << "    Smallest Fraction = " << fSmallestFraction << std::endl;
1042 
1043     os << "    No Integrat Vars  = " << fNoIntegrationVariables << std::endl;
1044     os << "    Min No Vars       = " << fMinNoVars << std::endl;
1045     os << "    Num-Vars          = " << fNoVars << std::endl;
1046     
1047     os << "    verbose level     = " << fVerboseLevel << std::endl;
1048 
1049     // auto p= const_cast<G4OldMagIntDriver*>(this);
1050     G4bool does = // p->DoesReIntegrate();
1051           const_cast<G4OldMagIntDriver*>(this)->DoesReIntegrate();    
1052     os << "    Reintegrates      = " << does  << std::endl;
1053 }
1054