Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4MagIntegratorDriver.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 ]

Diff markup

Differences between /geometry/magneticfield/src/G4MagIntegratorDriver.cc (Version 11.3.0) and /geometry/magneticfield/src/G4MagIntegratorDriver.cc (Version 7.0.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
 26 // G4MagInt_Driver implementation              << 
 27 //                                                 23 //
 28 // V.Grichine, 07.10.1996 - Created            <<  24 // $Id: G4MagIntegratorDriver.cc,v 1.42 2004/12/02 09:55:20 gcosmo Exp $
 29 // W.Wander, 28.01.1998 - Added ability for lo <<  25 // GEANT4 tag $Name: geant4-07-00-patch-01 $
 30 // J.Apostolakis, 08.11.2001 - Respect minimum <<  26 //
                                                   >>  27 // 
                                                   >>  28 //
                                                   >>  29 // Implementation for class G4MagInt_Driver
                                                   >>  30 // Tracking in space dependent magnetic field
                                                   >>  31 //
                                                   >>  32 // History of major changes:
                                                   >>  33 //  8 Nov 01  J. Apostolakis:   Respect minimum step in AccurateAdvance
                                                   >>  34 // 27 Jul 99  J. Apostolakis:   Ensured that AccurateAdvance does not loop 
                                                   >>  35 //                              due to very small eps & step size (precision)
                                                   >>  36 // 28 Jan 98  W. Wander:        Added ability for low order integrators
                                                   >>  37 //  7 Oct 96  V. Grichine       First version
 31 // -------------------------------------------     38 // --------------------------------------------------------------------
 32                                                    39 
 33 #include <iomanip>                             << 
 34                                                << 
 35 #include "globals.hh"                              40 #include "globals.hh"
 36 #include "G4SystemOfUnits.hh"                  <<  41 #include "geomdefs.hh"         //  for kCarTolerance
 37 #include "G4GeometryTolerance.hh"              <<  42 #include <iomanip>
 38 #include "G4MagIntegratorDriver.hh"                43 #include "G4MagIntegratorDriver.hh"
 39 #include "G4FieldTrack.hh"                         44 #include "G4FieldTrack.hh"
 40                                                    45 
 41 #ifdef   G4DEBUG_FIELD                         <<  46 //  Stepsize can increase by no more than 5.0
 42 #include "G4DriverReporter.hh"                 <<  47 //           and decrease by no more than 1/10. = 0.1
 43 #endif                                         <<  48 //
                                                   >>  49 const G4double G4MagInt_Driver::max_stepping_increase = 5.0;
                                                   >>  50 const G4double G4MagInt_Driver::max_stepping_decrease = 0.1;
 44                                                    51 
 45 // ------------------------------------------- <<  52 //  The (default) maximum number of steps is Base divided by the order of Stepper
                                                   >>  53 //
                                                   >>  54 const G4int  G4MagInt_Driver::fMaxStepBase = 250;  // Was 5000
                                                   >>  55 
                                                   >>  56 #ifndef G4NO_FIELD_STATISTICS
                                                   >>  57 #define G4FLD_STATS  1
                                                   >>  58 #endif
 46                                                    59 
 47 //  Constructor                                    60 //  Constructor
 48 //                                                 61 //
 49 G4MagInt_Driver::G4MagInt_Driver( G4double         62 G4MagInt_Driver::G4MagInt_Driver( G4double                hminimum, 
 50                                   G4MagIntegra <<  63           G4MagIntegratorStepper *pItsStepper,
 51                                   G4int        <<  64           G4int                   numComponents,
 52                                   G4int        <<  65           G4int                   statisticsVerbose)
 53   : fNoIntegrationVariables(numComponents),        66   : fNoIntegrationVariables(numComponents), 
                                                   >>  67     fMinNoVars(12), 
 54     fNoVars( std::max( fNoIntegrationVariables     68     fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )),
                                                   >>  69     fVerboseLevel(0),
                                                   >>  70     fNoTotalSteps(0),  fNoBadSteps(0), fNoSmallSteps(0), fNoInitialSmallSteps(0),
                                                   >>  71     fDyerr_max(0.0), fDyerr_mx2(0.0), 
                                                   >>  72     fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0), 
                                                   >>  73     fSumH_sm(0.0),   fSumH_lg(0.0),
 55     fStatisticsVerboseLevel(statisticsVerbose)     74     fStatisticsVerboseLevel(statisticsVerbose)
 56 {                                                  75 {  
 57   // In order to accomodate "Laboratory Time", <<  76 // In order to accomodate "Laboratory Time", which is [7], fMinNoVars=8 is required.
 58   // is required. For proper time of flight an <<  77 // For proper time of flight and spin,  fMinNoVars must be 12
 59                                                << 
 60   RenewStepperAndAdjust( pStepper );           << 
 61   fMinimumStep = hminimum;                     << 
 62                                                    78 
 63   fMaxNoSteps = fMaxStepBase / pIntStepper->In <<  79       // fNoVars= std::max( fNoVars, fMinNoVars );  
                                                   >>  80       RenewStepperAndAdjust( pItsStepper );
                                                   >>  81       fMinimumStep= hminimum;
                                                   >>  82       fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
 64 #ifdef G4DEBUG_FIELD                               83 #ifdef G4DEBUG_FIELD
 65   fVerboseLevel=2;                             <<  84       fVerboseLevel=2;
 66 #endif                                             85 #endif
 67                                                    86 
 68   if( (fVerboseLevel > 0) || (fStatisticsVerbo <<  87       if( (fVerboseLevel > 0) || (fStatisticsVerboseLevel > 1) ){
 69   {                                            <<  88    G4cout << "MagIntDriver version: Accur-Adv: invE_nS, QuickAdv-2sqrt with Statistics "
 70     G4cout << "MagIntDriver version: Accur-Adv << 
 71            << "invE_nS, QuickAdv-2sqrt with St << 
 72 #ifdef G4FLD_STATS                                 89 #ifdef G4FLD_STATS
 73            << " enabled "                      <<  90     << " enabled "
 74 #else                                              91 #else
 75            << " disabled "                     <<  92     << " disabled "
 76 #endif                                             93 #endif
 77            << G4endl;                          <<  94     << G4endl;
 78   }                                            <<  95       }
 79 }                                                  96 }
 80                                                    97 
 81 // ------------------------------------------- << 
 82                                                << 
 83 //  Destructor                                     98 //  Destructor
 84 //                                                 99 //
 85 G4MagInt_Driver::~G4MagInt_Driver()               100 G4MagInt_Driver::~G4MagInt_Driver()
 86 {                                                 101 { 
 87   if( fStatisticsVerboseLevel > 1 )            << 102      if( fStatisticsVerboseLevel > 1 ){
 88   {                                            << 103         PrintStatisticsReport() ;
 89     PrintStatisticsReport();                   << 104      }
 90   }                                            << 105      // Future: for default verbose level, print an understandable summary
 91 }                                                 106 }
 92                                                   107 
 93 // ------------------------------------------- << 108 
                                                   >> 109 // To add much printing for debugging purposes, uncomment this:
                                                   >> 110 // #define  G4DEBUG_FIELD 1    
 94                                                   111 
 95 G4bool                                            112 G4bool
 96 G4MagInt_Driver::AccurateAdvance(G4FieldTrack&    113 G4MagInt_Driver::AccurateAdvance(G4FieldTrack& y_current,
 97                                  G4double      << 114                      G4double     hstep,
 98                                  G4double      << 115                      G4double     eps,
 99                                  G4double      << 116          G4double hinitial )
100 {                                              << 117 //         const G4double dydx[6],    // We could may add this ??
101   // Runge-Kutta driver with adaptive stepsize << 118 
102   // values at y_current over hstep x2 with ac << 119 // Runge-Kutta driver with adaptive stepsize control. Integrate starting
103   // On output ystart is replaced by values at << 120 // values at y_current over hstep x2 with accuracy eps. 
104   // interval. RightHandSide is the right-hand << 121 // On output ystart is replaced by values at the end of the integration 
105   // The source is similar to odeint routine f << 122 // interval.  
                                                   >> 123 // RightHandSide is the right-hand side of ODE system. 
                                                   >> 124 // The source is similar to odeint routine from NRC p.721-722 .
106                                                   125 
107   G4int nstp, i;                               << 126 {
108   G4double x, hnext, hdid, h;                  << 127   G4int nstp, i, no_warnings=0;; 
                                                   >> 128   G4double x, hnext, hdid, h ;
109                                                   129 
110 #ifdef G4DEBUG_FIELD                              130 #ifdef G4DEBUG_FIELD
111   G4int no_warnings = 0;                       << 131   static G4int dbg=1;
112   static G4int dbg = 1;                        << 132   static G4int nStpPr=50;   // For debug printing of long integrations
113   static G4int nStpPr = 50;   // For debug pri << 
114   G4double ySubStepStart[G4FieldTrack::ncompSV    133   G4double ySubStepStart[G4FieldTrack::ncompSVEC];
115   G4FieldTrack  yFldTrkStart(y_current);          134   G4FieldTrack  yFldTrkStart(y_current);
116 #endif                                            135 #endif
117                                                   136 
118   G4double y[G4FieldTrack::ncompSVEC] = {0., 0 << 137   G4double y[G4FieldTrack::ncompSVEC], dydx[G4FieldTrack::ncompSVEC];
119   G4double dydx[G4FieldTrack::ncompSVEC] = {0. << 138   G4double ystart[G4FieldTrack::ncompSVEC], yEnd[G4FieldTrack::ncompSVEC]; 
120   G4double ystart[G4FieldTrack::ncompSVEC] = { << 
121   G4double yEnd[G4FieldTrack::ncompSVEC] = {0. << 
122   G4double  x1, x2;                               139   G4double  x1, x2;
123   G4bool succeeded = true;                     << 140   G4bool succeeded = true, lastStepSucceeded;
124                                                << 
125   G4double startCurveLength;                   << 
126                                                   141 
127   const G4int nvar = fNoVars;                  << 142   G4int  noFullIntegr=0, noSmallIntegr = 0 ;
                                                   >> 143   static G4int  noGoodSteps =0 ;  // Bad = chord > curve-len 
                                                   >> 144   const  int    nvar= fNoVars;
128                                                   145 
129   G4FieldTrack yStartFT(y_current);               146   G4FieldTrack yStartFT(y_current);
130                                                   147 
131   //  Ensure that hstep > 0                    << 148   //  Assume that hstep > 0 
132   //                                           << 
133   if( hstep <= 0.0 )                           << 
134   {                                            << 
135     if( hstep == 0.0 )                         << 
136     {                                          << 
137       std::ostringstream message;              << 
138       message << "Proposed step is zero; hstep << 
139       G4Exception("G4MagInt_Driver::AccurateAd << 
140                   "GeomField1001", JustWarning << 
141       return succeeded;                        << 
142     }                                          << 
143                                                << 
144     std::ostringstream message;                << 
145     message << "Invalid run condition." << G4e << 
146             << "Proposed step is negative; hst << 
147             << "Requested step cannot be negat << 
148     G4Exception("G4MagInt_Driver::AccurateAdva << 
149                 "GeomField0003", EventMustBeAb << 
150     return false;                              << 
151   }                                            << 
152                                                   149 
153   y_current.DumpToArray( ystart );                150   y_current.DumpToArray( ystart );
154                                                << 151   x1= y_current.GetCurveLength();
155   startCurveLength= y_current.GetCurveLength() << 
156   x1= startCurveLength;                        << 
157   x2= x1 + hstep;                                 152   x2= x1 + hstep;
158                                                   153 
159   if ( (hinitial > 0.0) && (hinitial < hstep)  << 154   if( (hinitial > 0.0) 
160     && (hinitial > perMillion * hstep) )       << 155       && (hinitial < hstep)
161   {                                            << 156       && (hinitial > perMillion * hstep) ){
162      h = hinitial;                                157      h = hinitial;
163   }                                            << 158   }else{
164   else  //  Initial Step size "h" defaults to  << 159      //  Initial Step size "h" defaults to the full interval
165   {                                            << 
166      h = hstep;                                   160      h = hstep;
167   }                                               161   }
168                                                   162 
169   x = x1;                                         163   x = x1;
170                                                   164 
171   for ( i=0; i<nvar; ++i)  { y[i] = ystart[i]; << 165   for(i=0;i<nvar;i++) y[i] = ystart[i] ;
172                                                << 
173   G4bool lastStep= false;                      << 
174   nstp = 1;                                    << 
175                                                   166 
176   do                                           << 167   G4bool  lastStep= false;
177   {                                            << 168   nstp=1;
178     G4ThreeVector StartPos( y[0], y[1], y[2] ) << 169   // G4double  lastStepThreshold = std::min( eps * hstep, Hmin() ); 
179                                                << 170 
180 #ifdef G4DEBUG_FIELD                           << 171   do{
181     G4double xSubStepStart= x;                 << 172      G4ThreeVector StartPos( y[0], y[1], y[2] );   
182     for (i=0; i<nvar; ++i)  { ySubStepStart[i] << 173 #    ifdef G4DEBUG_FIELD
183     yFldTrkStart.LoadFromArray(y, fNoIntegrati << 174        for(i=0;i<nvar;i++) ySubStepStart[i] = y[i] ;
184     yFldTrkStart.SetCurveLength(x);            << 175        yFldTrkStart.LoadFromArray(y);
185 #endif                                         << 176        yFldTrkStart.SetCurveLength(x);
186                                                << 177 #    endif
187     pIntStepper->RightHandSide( y, dydx );     << 178 
188     ++fNoTotalSteps;                           << 179      pIntStepper->RightHandSide( y, dydx );
189                                                << 180 
190     // Perform the Integration                 << 181      if( x+h > x2 ) {
191     //                                         << 182         h = x2 - x ;     // When stepsize overshoots, decrease it!
192     if( h > fMinimumStep )                     << 183         if( h < eps * hstep) {
193     {                                          << 184       lastStep = true;   //  Avoid numerous small last steps
194       OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ; << 185   }
195       //-------------------------------------- << 186      }
196 #ifdef G4DEBUG_FIELD                           << 187 
197       if (dbg>2)                               << 188      fNoTotalSteps++;
198       {                                        << 189      // Perform the Integration
199          // PrintStatus( ySubStepStart, xSubSt << 190      //      
200         G4DriverReporter::PrintStatus( ySubSte << 191      if( h > fMinimumStep ){ 
201       }                                        << 192         OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ;
                                                   >> 193         //--------------------------------------
                                                   >> 194         lastStepSucceeded= (hdid == h);   
                                                   >> 195 #ifdef  G4DEBUG_FIELD
                                                   >> 196       if(dbg>2) PrintStatus( ySubStepStart, x1, y, x, h,  nstp); // Only
202 #endif                                            197 #endif
203     }                                          << 198      }else{
204     else                                       << 199         G4FieldTrack yFldTrk( G4ThreeVector(0,0,0), 
205     {                                          << 200             G4ThreeVector(0,0,0), 0., 0., 0., 0. );
206       G4FieldTrack yFldTrk( G4ThreeVector(0,0, << 201         G4double dchord_step, dyerr, dyerr_len;  //  Must figure what to do with these
207                             G4ThreeVector(0,0, << 202   yFldTrk.LoadFromArray(y); 
208       G4double dchord_step, dyerr, dyerr_len;  << 203         yFldTrk.SetCurveLength( x );
209       yFldTrk.LoadFromArray(y, fNoIntegrationV << 204 
210       yFldTrk.SetCurveLength( x );             << 205         QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len ); 
211                                                << 206         //-----------------------------------------------------
212       QuickAdvance( yFldTrk, dydx, h, dchord_s << 207         // #ifdef G4DEBUG_FIELD
213       //-------------------------------------- << 208   // if(dbg>1) OneGoodStep(y,dydx,x,h,2*eps,hdid,hnext) ;
                                                   >> 209   // if(dbg>1) PrintStatus( ystart, x1, y, x, h, -nstp);  
214                                                   210 
215       yFldTrk.DumpToArray(y);                  << 211         yFldTrk.DumpToArray(y);    
216                                                   212 
217 #ifdef G4FLD_STATS                                213 #ifdef G4FLD_STATS
218       ++fNoSmallSteps;                         << 214   fNoSmallSteps++; 
219       if ( dyerr_len > fDyerr_max )  { fDyerr_ << 215         if( dyerr_len > fDyerr_max) fDyerr_max= dyerr_len;
220       fDyerrPos_smTot += dyerr_len;            << 216         fDyerrPos_smTot += dyerr_len;
221       fSumH_sm += h;  // Length total for 'sma << 217   fSumH_sm += h;  // Length total for 'small' steps
222       if (nstp==1)  { ++fNoInitialSmallSteps;  << 218         if(nstp<=1) fNoInitialSmallSteps++;
223 #endif                                            219 #endif
224 #ifdef G4DEBUG_FIELD                              220 #ifdef G4DEBUG_FIELD
225       if (dbg>1)                               << 221   if(dbg>1) {
226       {                                        << 222      if(fNoSmallSteps<2) PrintStatus( ySubStepStart, x1, y, x, h, -nstp);
227         if(fNoSmallSteps<2) { PrintStatus(ySub << 223      G4cout << "Another sub-min step, no " << fNoSmallSteps 
228         G4cout << "Another sub-min step, no "  << 224       << " of " << fNoTotalSteps << " this time " << nstp << G4endl; 
229                << " of " << fNoTotalSteps << " << 225      PrintStatus( ySubStepStart, x1, y, x, h,  nstp);   // Only this
230         PrintStatus( ySubStepStart, x1, y, x,  << 226      G4cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h 
231         G4cout << " dyerr= " << dyerr_len << " << 227       << " epsilon= " << eps << " hstep= " << hstep 
232                << " epsilon= " << eps << " hst << 228       << " h= " << h << " hmin= " << fMinimumStep
233                << " h= " << h << " hmin= " <<  << 229       << G4endl;
234       }                                        << 230         }
235 #endif                                         << 231 #endif  
236       if( h == 0.0 )                           << 232   dyerr = dyerr_len / h;
237       {                                        << 233   hdid= h;
238         G4Exception("G4MagInt_Driver::Accurate << 234         x += hdid;
239                     "GeomField0003", FatalExce << 235         // Compute suggested new step
240                     "Integration Step became Z << 236   hnext= ComputeNewStepSize( dyerr/eps, h);
241       }                                        << 237   // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
242       dyerr = dyerr_len / h;                   << 238   lastStepSucceeded= (dyerr<= eps);
243       hdid = h;                                << 239      }
244       x += hdid;                               << 
245                                                << 
246       // Compute suggested new step            << 
247       hnext = ComputeNewStepSize( dyerr/eps, h << 
248     }                                          << 
249                                                   240 
250     G4ThreeVector EndPos( y[0], y[1], y[2] );  << 241      if(lastStepSucceeded) noFullIntegr++ ; else noSmallIntegr++ ;
                                                   >> 242      G4ThreeVector EndPos( y[0], y[1], y[2] );
251                                                   243 
252 #ifdef  G4DEBUG_FIELD                             244 #ifdef  G4DEBUG_FIELD
253     if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))  << 245      if(dbg && (nstp>nStpPr)) {
254     {                                          << 246        G4cout << "hdid="  << std::setw(12) << hdid  << " "
255       if( nstp==nStpPr )  { G4cout << "***** M << 247         << "hnext=" << std::setw(12) << hnext << " " << G4endl;
256       G4cout << "MagIntDrv: " ;                << 248        PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp); 
257       G4cout << "hdid="  << std::setw(12) << h << 249      }
258              << "hnext=" << std::setw(12) << h << 
259              << "hstep=" << std::setw(12) << h << 
260              << G4endl;                        << 
261       PrintStatus( ystart, x1, y, x, h, (nstp= << 
262     }                                          << 
263 #endif                                            250 #endif
264                                                   251 
265     // Check the endpoint                      << 252      // Check the endpoint
266     G4double endPointDist= (EndPos-StartPos).m << 253      G4double endPointDist= (EndPos-StartPos).mag(); 
267     if ( endPointDist >= hdid*(1.+perMillion)  << 254      if( endPointDist >= hdid*(1.+perMillion) ){
268     {                                          << 255         fNoBadSteps  ++;
269       ++fNoBadSteps;                           << 256         // Issue a warning only for gross differences -
270                                                << 257         //   we understand how small difference occur.
271       // Issue a warning only for gross differ << 258         if( endPointDist >= hdid*(1.+perThousand) ){ 
272       // we understand how small difference oc << 259 #ifdef  G4DEBUG_FIELD
273       if ( endPointDist >= hdid*(1.+perThousan << 260            if(dbg){
274       {                                        << 261         WarnEndPointTooFar ( endPointDist, hdid, eps, dbg ); 
275 #ifdef G4DEBUG_FIELD                           << 262         G4cerr << "  Total steps:  bad " << fNoBadSteps << " good " << noGoodSteps << " current h= " << hdid << G4endl;
276         if (dbg)                               << 263         // G4cerr << "Mid:EndPtFar> ";
277         {                                      << 264         PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);  
278           WarnEndPointTooFar ( endPointDist, h << 265                 // Potentially add as arguments:  <dydx> - as Initial Force
279           G4cerr << "  Total steps:  bad " <<  << 266            }
280                  << " current h= " << hdid <<  << 267 #endif
281           PrintStatus( ystart, x1, y, x, hstep << 268      no_warnings++;
282         }                                      << 269   }
283         ++no_warnings;                         << 270      } else { // ie (!dbg)
284 #endif                                         << 271         noGoodSteps ++;
285       }                                        << 272      } 
286     }                                          << 273 // #endif
287                                                << 274 
288     //  Avoid numerous small last steps        << 275      // Check the proposed next stepsize
289     if( (h < eps * hstep) || (h < fSmallestFra << 276      if(std::fabs(hnext) <= Hmin())
290     {                                          << 277      {
291       // No more integration -- the next step  << 
292       lastStep = true;                         << 
293     }                                          << 
294     else                                       << 
295     {                                          << 
296       // Check the proposed next stepsize      << 
297       if(std::fabs(hnext) <= Hmin())           << 
298       {                                        << 
299 #ifdef  G4DEBUG_FIELD                             278 #ifdef  G4DEBUG_FIELD
300         // If simply a very small interval is     279         // If simply a very small interval is being integrated, do not warn
301         if( (x < x2 * (1-eps) ) &&        // T << 280         if( (x < x2 * (1-eps) ) &&     //  The last step can be small: it's OK
302             (std::fabs(hstep) > Hmin()) ) // a << 281             (std::fabs(hstep) > Hmin())     //   and if we are asked, it's OK
303         {                                      << 282             //   && (hnext < hstep * PerThousand ) 
304           if(dbg>0)                            << 283     ){
305           {                                    << 284              if(dbg>0){ 
306             WarnSmallStepSize( hnext, hstep, h << 285     //  Issue WARNING
307             PrintStatus( ystart, x1, y, x, hst << 286     WarnSmallStepSize( hnext, hstep, h, x-x1, nstp ); 
308           }                                    << 287     // G4cerr << "Mid:SmallStep> ";
309           ++no_warnings;                       << 288     PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
310         }                                      << 289        }
311 #endif                                         << 290        no_warnings++;
                                                   >> 291      }
                                                   >> 292 #endif
                                                   >> 293         // else 
                                                   >> 294   //   succeeded = false;  // Meaningful only if we break out of the loop.
                                                   >> 295   // 
                                                   >> 296         // lastStep = true;   //  Make this the last step ... Dubious now
                                                   >> 297 
312         // Make sure that the next step is at     298         // Make sure that the next step is at least Hmin.
313         h = Hmin();                               299         h = Hmin();
314       }                                        << 300      }else{
315       else                                     << 301         h = hnext ;
316       {                                        << 302      }
317         h = hnext;                             << 303 
318       }                                        << 304   }while ( ((nstp++)<=fMaxNoSteps)
319                                                << 305           && (x < x2)           //  Have we reached the end ?
320       //  Ensure that the next step does not o << 306                                 //   --> a better test might be x-x2 > an_epsilon
321       if ( x+h > x2 )                          << 307           && (!lastStep)
322       {                // When stepsize oversh << 308          );
323         h = x2 - x ;   // Must cope with diffi << 
324       }                // issues if hstep << x << 
325                                                << 
326       if ( h == 0.0 )                          << 
327       {                                        << 
328         // Cannot progress - accept this as la << 
329         lastStep = true;                       << 
330 #ifdef G4DEBUG_FIELD                           << 
331         if (dbg>2)                             << 
332         {                                      << 
333           int prec= G4cout.precision(12);      << 
334           G4cout << "Warning: G4MagIntegratorD << 
335                  << G4endl                     << 
336                  << "  Integration step 'h' be << 
337                  << h << " due to roundoff. "  << 
338                  << " Calculated as difference << 
339                  << "  Forcing termination of  << 
340           G4cout.precision(prec);              << 
341         }                                      << 
342 #endif                                         << 
343       }                                        << 
344     }                                          << 
345   } while ( ((++nstp)<=fMaxNoSteps) && (x < x2 << 
346   // Loop checking, 07.10.2016, J. Apostolakis << 
347                                                   309 
348      // Have we reached the end ?              << 310   succeeded=  (x>=x2);  // If it was a "forced" last step
349      // --> a better test might be x-x2 > an_e << 
350                                                   311 
351   succeeded = (x>=x2);  // If it was a "forced << 312   for(i=0;i<nvar;i++)  yEnd[i] = y[i] ;
352                                                << 
353   for (i=0; i<nvar; ++i)  { yEnd[i] = y[i]; }  << 
354                                                   313 
355   // Put back the values.                         314   // Put back the values.
356   y_current.LoadFromArray( yEnd, fNoIntegratio << 315   y_current.LoadFromArray( yEnd );
357   y_current.SetCurveLength( x );                  316   y_current.SetCurveLength( x );
358                                                   317 
359   if(nstp > fMaxNoSteps)                       << 318   if(nstp > fMaxNoSteps){
360   {                                            << 319      no_warnings++;
361     succeeded = false;                         << 320      succeeded = false;
362 #ifdef G4DEBUG_FIELD                           << 321 #ifdef  G4DEBUG_FIELD
363     ++no_warnings;                             << 322         if(dbg){
364     if (dbg)                                   << 323      WarnTooManyStep( x1, x2, x );  //  Issue WARNING
365     {                                          << 324      PrintStatus( yEnd, x1, y, x, hstep, -nstp);
366       WarnTooManyStep( x1, x2, x );  //  Issue << 325         }
367       PrintStatus( yEnd, x1, y, x, hstep, -nst << 
368     }                                          << 
369 #endif                                            326 #endif
370   }                                               327   }
371                                                   328 
372 #ifdef G4DEBUG_FIELD                              329 #ifdef G4DEBUG_FIELD
373   if( dbg && no_warnings )                     << 330   if( dbg && no_warnings ){
374   {                                            << 331      G4cerr << "G4MagIntegratorDriver exit status: no-steps " << nstp <<G4endl;
375     G4cerr << "G4MagIntegratorDriver exit stat << 332      PrintStatus( yEnd, x1, y, x, hstep, nstp);
376     PrintStatus( yEnd, x1, y, x, hstep, nstp); << 
377   }                                               333   }
378 #endif                                            334 #endif
379                                                   335 
380   return succeeded;                               336   return succeeded;
381 }  // end of AccurateAdvance ................. << 
382                                                   337 
383 // ------------------------------------------- << 338 }  // end of AccurateAdvance ...........................
384                                                   339 
385 void                                              340 void
386 G4MagInt_Driver::WarnSmallStepSize( G4double h    341 G4MagInt_Driver::WarnSmallStepSize( G4double hnext, G4double hstep, 
387                                     G4double h << 342             G4double h, G4double xDone,
388                                     G4int nstp << 343             G4int nstp)
389 {                                                 344 {
390   static G4ThreadLocal G4int noWarningsIssued  << 345   static G4int noWarningsIssued =0;
391   const  G4int maxNoWarnings = 10;   // Number << 346   const  G4int maxNoWarnings =  10;   // Number of verbose warnings
392   std::ostringstream message;                  << 347   if( (noWarningsIssued < maxNoWarnings) || fVerboseLevel > 10 ){
393   if( (noWarningsIssued < maxNoWarnings) || fV << 348     G4cerr<< " Warning (G4MagIntegratorDriver::AccurateAdvance): The stepsize for the " 
394   {                                            << 349     << " next iteration=" << hnext << " is too small " 
395     message << "The stepsize for the next iter << 350     << "- in Step number " << nstp << "." << G4endl;
396             << ", is too small - in Step numbe << 351     G4cerr << "     The minimum for the driver is " << Hmin()  << G4endl ;
397             << "The minimum for the driver is  << 352     G4cerr << "     Requested integr. length was " << hstep << " ." << G4endl ;
398             << "Requested integr. length was " << 353     G4cerr << "     The size of this sub-step was " << h     << " ." << G4endl ;
399             << "The size of this sub-step was  << 354     G4cerr << "     The integrations has already gone " << xDone << G4endl ;
400             << "The integrations has already g << 355   }else{
401   }                                            << 356     G4cerr<< " G4MagInt_Driver: Too small 'next' step " << hnext 
402   else                                         << 357     << " step-no "  << nstp ;                     // << G4setw(4)
403   {                                            << 358     G4cerr << "  this sub-step " << h     
404     message << "Too small 'next' step " << hne << 359      << "  req_tot_len " << hstep 
405             << ", step-no: " << nstp << G4endl << 360      << "  done " << xDone 
406             << ", this sub-step: " << h        << 361      << "  min " << Hmin() 
407             << ",  req_tot_len: " << hstep     << 362      << G4endl ;
408             << ", done: " << xDone << ", min:  << 
409   }                                               363   }
410   G4Exception("G4MagInt_Driver::WarnSmallStepS << 364   noWarningsIssued++;
411               JustWarning, message);           << 
412   ++noWarningsIssued;                          << 
413 }                                                 365 }
414                                                   366 
415 // ------------------------------------------- << 
416                                                << 
417 void                                              367 void
418 G4MagInt_Driver::WarnTooManyStep( G4double x1s    368 G4MagInt_Driver::WarnTooManyStep( G4double x1start, 
419                                   G4double x2e << 369           G4double x2end, 
420                                   G4double xCu << 370           G4double xCurrent)
421 {                                                 371 {
422    std::ostringstream message;                 << 372     G4cerr << " Warning (G4MagIntegratorDriver): The number of steps " 
423    message << "The number of steps used in the << 373          << "used in the Integration driver (Runge-Kutta) is too many.  "
424            << " (Runge-Kutta) is too many." << << 374    << G4endl ;
425            << "Integration of the interval was << 375     G4cerr << "Integration of the interval was not completed - only a " 
426            << "Only a " << (xCurrent-x1start)* << 376          << (xCurrent-x1start)*100/(x2end-x1start)
427            << " % fraction of it was done.";   << 377      <<" % fraction of it was Done." << G4endl;
428    G4Exception("G4MagInt_Driver::WarnTooManySt << 
429                JustWarning, message);          << 
430 }                                                 378 }
431                                                   379 
432 // ------------------------------------------- << 
433                                                << 
434 void                                              380 void
435 G4MagInt_Driver::WarnEndPointTooFar (G4double     381 G4MagInt_Driver::WarnEndPointTooFar (G4double endPointDist, 
436                                      G4double  << 382              G4double   h , 
437                                      G4double  << 383              G4double  eps,
438                                      G4int     << 384              G4int     dbg)
439 {                                              << 385 {
440   static G4ThreadLocal G4double maxRelError =  << 386   static G4double maxRelError= 0.0, maxRelError_last_printed=0.0;
441   G4bool isNewMax, prNewMax;                   << 387   G4bool isNewMax, prNewMax;
442                                                << 388 
443   isNewMax = endPointDist > (1.0 + maxRelError << 389         isNewMax = endPointDist > (1.0 + maxRelError) * h;
444   prNewMax = endPointDist > (1.0 + 1.05 * maxR << 390         prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
445   if( isNewMax ) { maxRelError= endPointDist / << 391   if( isNewMax )
446                                                << 392      maxRelError= endPointDist / h - 1.0; 
447   if( (dbg != 0) && (h > G4GeometryTolerance:: << 393         if( prNewMax )
448           && ( (dbg>1) || prNewMax || (endPoin << 394      maxRelError_last_printed = maxRelError;
449   {                                            << 395 
450     static G4ThreadLocal G4int noWarnings = 0; << 396         if(    dbg 
451     std::ostringstream message;                << 397       && (h > kCarTolerance) 
452     if( (noWarnings++ < 10) || (dbg>2) )       << 398       && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+eps) ) ) 
453     {                                          << 399           ){ 
454       message << "The integration produced an  << 400            static G4int noWarnings = 0;
455               << "is further from the start-po << 401            if( (noWarnings ++ < 10) || (dbg>2) ){
456               << G4endl;                       << 402         G4cerr << " Warning (G4MagIntegratorDriver): "
457     }                                          << 403          << " The integration produced an endpoint which " << G4endl
458     message << "  Distance of endpoints = " << << 404          << "   is further from the startpoint than the curve length." << G4endl; 
459             << ", curve length = " << h << G4e << 405           
460             << "  Difference (curveLen-endpDis << 406         G4cerr << "   Distance of endpoints = " << endPointDist
461             << ", relative = " << (h-endPointD << 407          << "  curve length = " <<  h
462             << ", epsilon =  " << eps;         << 408          << "  Difference (curveLen-endpDist)= " << (h - endPointDist)
463     G4Exception("G4MagInt_Driver::WarnEndPoint << 409          << "  relative = " << (h-endPointDist) / h 
464                 JustWarning, message);         << 410          << "  epsilon =  " << eps 
465   }                                            << 411          << G4endl;
                                                   >> 412      }else{
                                                   >> 413         G4cerr << " Warning:" 
                                                   >> 414          << "  dist_e= " << endPointDist
                                                   >> 415          << "  h_step = " <<  h
                                                   >> 416          << "  Diff (hs-ed)= " << (h - endPointDist)
                                                   >> 417          << "  rel = " << (h-endPointDist) / h 
                                                   >> 418          << "  eps = " << eps
                                                   >> 419          << " (from G4MagInt_Driver)" << G4endl;
                                                   >> 420      }
                                                   >> 421   }
466 }                                                 422 }
467                                                << 
468 // -------------------------------------------    423 // ---------------------------------------------------------
469                                                   424 
470 void                                              425 void
471 G4MagInt_Driver::OneGoodStep(      G4double y[    426 G4MagInt_Driver::OneGoodStep(      G4double y[],        // InOut
472                              const G4double dy << 427            const G4double dydx[],
473                                    G4double& x << 428            G4double& x,         // InOut
474                                    G4double ht << 429                  G4double htry,
475                                    G4double ep << 430                  G4double eps_rel_max,
476                                    G4double& h << 431            G4double& hdid,      // Out
477                                    G4double& h << 432            G4double& hnext )    // Out
478                                                   433 
479 // Driver for one Runge-Kutta Step with monito    434 // Driver for one Runge-Kutta Step with monitoring of local truncation error
480 // to ensure accuracy and adjust stepsize. Inp    435 // to ensure accuracy and adjust stepsize. Input are dependent variable
481 // array y[0,...,5] and its derivative dydx[0,    436 // array y[0,...,5] and its derivative dydx[0,...,5] at the
482 // starting value of the independent variable     437 // starting value of the independent variable x . Also input are stepsize
483 // to be attempted htry, and the required accu    438 // to be attempted htry, and the required accuracy eps. On output y and x
484 // are replaced by their new values, hdid is t    439 // are replaced by their new values, hdid is the stepsize that was actually
485 // accomplished, and hnext is the estimated ne    440 // accomplished, and hnext is the estimated next stepsize. 
486 // This is similar to the function rkqs from t    441 // This is similar to the function rkqs from the book:
487 // Numerical Recipes in C: The Art of Scientif    442 // Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
488 // Edition, by William H. Press, Saul A. Teuko    443 // Edition, by William H. Press, Saul A. Teukolsky, William T.
489 // Vetterling, and Brian P. Flannery (Cambridg    444 // Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
490 // 16.2 Adaptive StepSize Control for Runge-Ku    445 // 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
491                                                   446 
492 {                                                 447 {
493   G4double errmax_sq;                          << 448       // G4double errpos_rel_sq, errvel_rel_sq
494   G4double h, htemp, xnew ;                    << 449       G4double errmax_sq;
                                                   >> 450       // G4double errmax; 
                                                   >> 451       G4double h, htemp, xnew ;
495                                                   452 
496   G4double yerr[G4FieldTrack::ncompSVEC], ytem << 453       G4double yerr[G4FieldTrack::ncompSVEC], ytemp[G4FieldTrack::ncompSVEC];
497                                                   454 
498   h = htry ; // Set stepsize to the initial tr << 455       h = htry ; // Set stepsize to the initial trial value
499                                                   456 
500   G4double inv_eps_vel_sq = 1.0 / (eps_rel_max << 457       // G4double inv_epspos_sq = 1.0 / eps * eps; 
                                                   >> 458       G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max); 
501                                                   459 
502   G4double errpos_sq = 0.0;    // square of di << 460       G4double errpos_sq=0.0;    // square of displacement error
503   G4double errvel_sq = 0.0;    // square of mo << 461       G4double errvel_sq=0.0;    // square of momentum vector difference
504   G4double errspin_sq = 0.0;   // square of sp << 
505                                                   462 
506   const G4int max_trials=100;                  << 463       G4int iter;
507                                                   464 
508   G4ThreeVector Spin(y[9],y[10],y[11]);        << 465       static G4int tot_no_trials=0; 
509   G4double spin_mag2 = Spin.mag2();            << 466       const G4int max_trials=100; 
510   G4bool hasSpin = (spin_mag2 > 0.0);          << 
511                                                << 
512   for (G4int iter=0; iter<max_trials; ++iter)  << 
513   {                                            << 
514     pIntStepper-> Stepper(y,dydx,h,ytemp,yerr) << 
515     //            *******                      << 
516     G4double eps_pos = eps_rel_max * std::max( << 
517     G4double inv_eps_pos_sq = 1.0 / (eps_pos*e << 
518                                                << 
519     // Evaluate accuracy                       << 
520     //                                         << 
521     errpos_sq =  sqr(yerr[0]) + sqr(yerr[1]) + << 
522     errpos_sq *= inv_eps_pos_sq; // Scale rela << 
523                                                << 
524     // Accuracy for momentum                   << 
525     G4double magvel_sq=  sqr(y[3]) + sqr(y[4]) << 
526     G4double sumerr_sq =  sqr(yerr[3]) + sqr(y << 
527     if( magvel_sq > 0.0 )                      << 
528     {                                          << 
529        errvel_sq = sumerr_sq / magvel_sq;      << 
530     }                                          << 
531     else                                       << 
532     {                                          << 
533        std::ostringstream message;             << 
534        message << "Found case of zero momentum << 
535                << "- iteration= " << iter << " << 
536        G4Exception("G4MagInt_Driver::OneGoodSt << 
537                    "GeomField1001", JustWarnin << 
538        errvel_sq = sumerr_sq;                  << 
539     }                                          << 
540     errvel_sq *= inv_eps_vel_sq;               << 
541     errmax_sq = std::max( errpos_sq, errvel_sq << 
542                                                   467 
543     if( hasSpin )                              << 468       for (iter=0; iter<max_trials ;iter++)
544     {                                          << 469       {
545       // Accuracy for spin                     << 470     tot_no_trials++;
546       errspin_sq =  ( sqr(yerr[9]) + sqr(yerr[ << 471     pIntStepper-> Stepper(y,dydx,h,ytemp,yerr); 
547                     /  spin_mag2; // ( sqr(y[9 << 472           //            *******
548       errspin_sq *= inv_eps_vel_sq;            << 473           G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep); 
549       errmax_sq = std::max( errmax_sq, errspin << 474     G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos); 
550     }                                          << 475 
                                                   >> 476     // Evaluate accuracy
                                                   >> 477     //
                                                   >> 478     errpos_sq =  sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
                                                   >> 479     // errpos_sq /= eps_pos*eps_pos; // Scale to tolerance
                                                   >> 480     errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
                                                   >> 481 
                                                   >> 482           // Accuracy for momentum
                                                   >> 483           errvel_sq =  (sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) )
                                                   >> 484                      / (sqr(y[3]) + sqr(y[4]) + sqr(y[5]) );
                                                   >> 485           // errvel_sq /= eps_rel_max*eps_rel_max; 
                                                   >> 486           errvel_sq *= inv_eps_vel_sq;
                                                   >> 487 
                                                   >> 488           errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
                                                   >> 489           // errmax = std::sqrt( errmax_sq );
                                                   >> 490     if(errmax_sq <= 1.0 ) break ; // Step succeeded. 
                                                   >> 491 
                                                   >> 492     // Step failed; compute the size of retrial Step.
                                                   >> 493     htemp = GetSafety()*h* std::pow( errmax_sq, 0.5*GetPshrnk() );
                                                   >> 494 
                                                   >> 495     if(htemp >= 0.1*h) h = htemp ;  // Truncation error too large,
                                                   >> 496     else h = 0.1*h ;                // reduce stepsize, but no more
                                                   >> 497             // than a factor of 10
                                                   >> 498     xnew = x + h ;
                                                   >> 499     if(xnew == x) {
                                                   >> 500        G4cerr<<"G4MagIntegratorDriver::OneGoodStep: Stepsize underflow in Stepper "<<G4endl ;
                                                   >> 501        G4cerr<<"  Step's start x=" << x << " and end x= " << xnew 
                                                   >> 502        << " are equal !! " << G4endl
                                                   >> 503        <<"  Due to step-size= " << h 
                                                   >> 504                    << " . Note that input step was " << htry << G4endl;
                                                   >> 505        break;
                                                   >> 506     }
                                                   >> 507       }
                                                   >> 508       // tot_no_trials+= (iter+1); 
551                                                   509 
552     if ( errmax_sq <= 1.0 )  { break; } // Ste << 510 #ifdef G4FLD_STATS
                                                   >> 511       // Sum of squares of position error // and momentum dir (underestimated)
                                                   >> 512       fSumH_lg += h; 
                                                   >> 513       fDyerrPos_lgTot += errpos_sq;  //  + errvel_last_sq * h * h ; 
                                                   >> 514       fDyerrVel_lgTot += errvel_sq * h * h; 
                                                   >> 515 #endif
553                                                   516 
554     // Step failed; compute the size of retria << 517       // Compute size of next Step
555     htemp = GetSafety() * h * std::pow( errmax << 518       if(errmax_sq > errcon*errcon) 
                                                   >> 519                hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow()) ;
                                                   >> 520       else hnext = max_stepping_increase*h ;
                                                   >> 521                      // No more than a factor of 5 increase
556                                                   522 
557     if (htemp >= 0.1*h)  { h = htemp; }  // Tr << 523       x += (hdid = h) ;
558     else  { h = 0.1*h; }                 // re << 
559                                          // th << 
560     xnew = x + h;                              << 
561     if(xnew == x)                              << 
562     {                                          << 
563       std::ostringstream message;              << 
564       message << "Stepsize underflow in Steppe << 
565               << "- Step's start x=" << x << " << 
566               << " are equal !! " << G4endl    << 
567               << "  Due to step-size= " << h   << 
568               << ". Note that input step was " << 
569       G4Exception("G4MagInt_Driver::OneGoodSte << 
570                   "GeomField1001", JustWarning << 
571       break;                                   << 
572     }                                          << 
573   }                                            << 
574                                                   524 
575   // Compute size of next Step                 << 525       int i;
576   if (errmax_sq > errcon*errcon)               << 526       const int nvar= fNoIntegrationVariables;   
577   {                                            << 527       for(i=0;i<nvar;i++) y[i] = ytemp[i] ;
578     hnext = GetSafety()*h*std::pow(errmax_sq,  << 
579   }                                            << 
580   else                                         << 
581   {                                            << 
582     hnext = max_stepping_increase*h ; // No mo << 
583   }                                            << 
584   x += (hdid = h);                             << 
585                                                   528 
586   for(G4int k=0; k<fNoIntegrationVariables; ++ << 529       return ;
587                                                   530 
588   return;                                      << 531 }   // end of  OneGoodStep .............................
589 }                                              << 
590                                                   532 
591 //--------------------------------------------    533 //----------------------------------------------------------------------
592                                                   534 
                                                   >> 535 #ifdef QUIK_ADVANCE_NEW2
593 // QuickAdvance just tries one Step - it does     536 // QuickAdvance just tries one Step - it does not ensure accuracy
594 //                                                537 //
595 G4bool G4MagInt_Driver::QuickAdvance(G4FieldTr << 538 //   This original interface does not return individual element errors
596                                const G4double  << 539 //    It is kept only for compatibility, and will be obsolete as of G4 6.0
597                                      G4double  << 540 //
598                                      G4double& << 541 G4bool  G4MagInt_Driver::QuickAdvance(       
599                                      G4double& << 542           G4FieldTrack& y_posvel,         // INOUT
600                                      G4double& << 543                 const G4double     dydx[],  
601 {                                              << 544                       G4double     hstep,       // In
602   G4Exception("G4MagInt_Driver::QuickAdvance() << 545                 G4double&    dchord_step,
603               FatalException, "Not yet impleme << 546                 G4double&    dyerr )
604                                                << 547 {
605   // Use the parameters of this method, to ple << 548     G4double dyerr_pos_sq=0.0, dyerr_mom_rel_sq=0.0; // dyerr_ener_sq=0.0;
606   //                                           << 549     G4double dyerr_pos; 
607   dchord_step = dyerr_pos_sq = hstep * hstep * << 550 
608   dyerr_mom_rel_sq = y_posvel.GetPosition().ma << 551     QuickAdvance( y_posvel,  dydx,  hstep,      
609   return true;                                 << 552       dchord_step, dyerr_pos_sq, dyerr_mom_rel_sq);
                                                   >> 553                           //  , dyerr_ener_sq ); 
                                                   >> 554     ... 
                                                   >> 555     // Calculate dyerr from the above -- as at the end of the func below
610 }                                                 556 }
                                                   >> 557 #endif
611                                                   558 
612 //--------------------------------------------    559 //----------------------------------------------------------------------
                                                   >> 560 // QuickAdvance just tries one Step - it does not ensure accuracy
                                                   >> 561 //
                                                   >> 562 G4bool  G4MagInt_Driver::QuickAdvance(       
                                                   >> 563           G4FieldTrack& y_posvel,         // INOUT
                                                   >> 564                 const G4double     dydx[],  
                                                   >> 565                       G4double     hstep,       // In
                                                   >> 566                 G4double&    dchord_step,
                                                   >> 567                   G4double&    dyerr_pos_sq,
                                                   >> 568                 G4double&    dyerr_mom_rel_sq
                                                   >> 569           //    G4double&    dyerr_ener_sq // Future
                                                   >> 570  )  
                                                   >> 571 {
                                                   >> 572       G4Exception("G4MagInt_Driver::QuickAdvance()", "NotImplemented",
                                                   >> 573                   FatalException, "Not yet implemented."); 
                                                   >> 574 
                                                   >> 575       // Use the parameters of this method, to please compiler
                                                   >> 576       dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0]; 
                                                   >> 577       dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
                                                   >> 578       return true;
                                                   >> 579 }
                                                   >> 580 
                                                   >> 581 G4bool  G4MagInt_Driver::QuickAdvance(       
                                                   >> 582           G4FieldTrack& y_posvel,         // INOUT
                                                   >> 583                 const G4double     dydx[],  
                                                   >> 584                       G4double     hstep,       // In
                                                   >> 585                 G4double&    dchord_step,
                                                   >> 586                 G4double&    dyerr )
                                                   >> 587 {
                                                   >> 588     G4double   dyerr_pos_sq,  dyerr_mom_rel_sq;  
                                                   >> 589     G4double yerr_vec[G4FieldTrack::ncompSVEC], yarrin[G4FieldTrack::ncompSVEC], yarrout[G4FieldTrack::ncompSVEC]; 
                                                   >> 590     G4double s_start;
                                                   >> 591     // G4double dyerr_len=0.0;  // , dyerr_vel, vel_mag;
                                                   >> 592     G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
                                                   >> 593 
                                                   >> 594     static G4int no_call=0; 
                                                   >> 595     no_call ++; 
                                                   >> 596 
                                                   >> 597     // Move data into array
                                                   >> 598     y_posvel.DumpToArray( yarrin );      //  yarrin  <== y_posvel 
                                                   >> 599     s_start = y_posvel.GetCurveLength();
613                                                   600 
614 G4bool G4MagInt_Driver::QuickAdvance(G4FieldTr << 601     // Do an Integration Step
615                                const G4double  << 602     pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ; 
616                                      G4double  << 603     //            *******
617                                      G4double& << 
618                                      G4double& << 
619 {                                              << 
620   G4double dyerr_pos_sq, dyerr_mom_rel_sq;     << 
621   G4double yerr_vec[G4FieldTrack::ncompSVEC],  << 
622            yarrin[G4FieldTrack::ncompSVEC], ya << 
623   G4double s_start;                            << 
624   G4double dyerr_mom_sq, vel_mag_sq, inv_vel_m << 
625                                                << 
626   // Move data into array                      << 
627   y_posvel.DumpToArray( yarrin );      //  yar << 
628   s_start = y_posvel.GetCurveLength();         << 
629                                                << 
630   // Do an Integration Step                    << 
631   pIntStepper-> Stepper(yarrin, dydx, hstep, y << 
632                                                << 
633   // Estimate curve-chord distance             << 
634   dchord_step= pIntStepper-> DistChord();      << 
635                                                << 
636   // Put back the values.  yarrout ==> y_posve << 
637   y_posvel.LoadFromArray( yarrout, fNoIntegrat << 
638   y_posvel.SetCurveLength( s_start + hstep );  << 
639                                                << 
640 #ifdef  G4DEBUG_FIELD                          << 
641   if(fVerboseLevel>2)                          << 
642   {                                            << 
643     G4cout << "G4MagIntDrv: Quick Advance" <<  << 
644     PrintStatus( yarrin, s_start, yarrout, s_s << 
645   }                                            << 
646 #endif                                         << 
647                                                   604 
648   // A single measure of the error             << 605     // Estimate curve-chord distance
649   //      TO-DO :  account for  energy,  spin, << 606     dchord_step= pIntStepper-> DistChord();
650   vel_mag_sq   = ( sqr(yarrout[3])+sqr(yarrout << 607     //                         *********
651   inv_vel_mag_sq = 1.0 / vel_mag_sq;           << 608 
652   dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_v << 609     // Put back the values.
653   dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_v << 610     y_posvel.LoadFromArray( yarrout );   //  yarrout ==> y_posvel
654   dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_ma << 611     y_posvel.SetCurveLength( s_start + hstep );
655                                                << 612 
656   // Calculate also the change in the momentum << 613     // A single measure of the error   
657   // G4double veloc_square = y_posvel.GetVeloc << 614     //      TO-DO :  account for  energy,  spin, ... ? 
658   // ...                                       << 615     vel_mag_sq   = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
                                                   >> 616     inv_vel_mag_sq = 1.0 / vel_mag_sq; 
                                                   >> 617     dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
                                                   >> 618     dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
                                                   >> 619 
                                                   >> 620     dyerr_mom_rel_sq =  dyerr_mom_sq * inv_vel_mag_sq;
                                                   >> 621 
                                                   >> 622     //// Calculate also the change in the momentum squared also ???
                                                   >> 623     // G4double veloc_square = y_posvel.GetVelocity().mag2();
                                                   >> 624     // ...
659                                                   625 
660 #ifdef RETURN_A_NEW_STEP_LENGTH                   626 #ifdef RETURN_A_NEW_STEP_LENGTH
661   // The following step cannot be done here be << 627     // The following step cannot be done here because "eps" is not known.
662   dyerr_len = std::sqrt( dyerr_len_sq );       << 628     dyerr_len = std::sqrt( dyerr_len_sq ); 
663   dyerr_len_sq /= eps ;                        << 629     dyerr_len_sq /= eps ;
664                                                   630 
665   // Look at the velocity deviation ?          << 631     // Look at the velocity deviation ?
666   //  sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(ye << 632     //  sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
667                                                   633 
668   // Set suggested new step                    << 634    // Set suggested new step
669   hstep = ComputeNewStepSize( dyerr_len, hstep << 635     hstep= ComputeNewStepSize( dyerr_len, hstep);
670 #endif                                            636 #endif
671                                                   637 
672   if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr( << 638     if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) ) {
673   {                                            << 639        dyerr = std::sqrt(dyerr_pos_sq);
674     dyerr = std::sqrt(dyerr_pos_sq);           << 640     }else{
675   }                                            << 641        // Scale it to the current step size - for now
676   else                                         << 642        dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
677   {                                            << 643     }
678     // Scale it to the current step size - for << 
679     dyerr = std::sqrt(dyerr_mom_rel_sq) * hste << 
680   }                                            << 
681                                                   644 
682   return true;                                 << 645     return true;
683 }                                                 646 }
684                                                   647 
685 // ------------------------------------------- << 
686                                                << 
687 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT                 648 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT
688 G4bool  G4MagInt_Driver::QuickAdvance(G4double << 649 G4bool  G4MagInt_Driver::QuickAdvance(       
689                                 const G4double << 650           G4double     yarrin[],        // IN
690                                       G4double << 651                 const G4double     dydx[],  
691                                       G4double << 652                       G4double     hstep,       // In
692                                       G4double << 653           G4double     yarrout[],
693                                       G4double << 654                 G4double&    dchord_step,
                                                   >> 655                 G4double&    dyerr )      // in length
694 {                                                 656 {
695   G4Exception("G4MagInt_Driver::QuickAdvance() << 657    G4Exception("G4MagInt_Driver::QuickAdvance()", "NotImplemented",
696               FatalException, "Not yet impleme << 658                 FatalException, "Not yet implemented.");
697   dyerr = dchord_step = hstep * yarrin[0] * dy << 659 
698   yarrout[0]= yarrin[0];                       << 660    dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
                                                   >> 661    yarrout[0]= yarrin[0];
699 }                                                 662 }
700 #endif                                            663 #endif 
701                                                   664 
702 // -------------------------------------------    665 // --------------------------------------------------------------------------
703                                                << 666 //  This method computes new step sizes - but does not limit changes to
704 // This method computes new step sizes - but d << 667 //   within  certain factors
705 // within  certain factors                     << 
706 //                                                668 // 
707 G4double G4MagInt_Driver::                     << 669 
708 ComputeNewStepSize_WithoutReductionLimit(G4dou << 670 G4double 
709                    G4double  hstepCurrent)  // << 671 G4MagInt_Driver::ComputeNewStepSize( 
                                                   >> 672                           G4double  errMaxNorm,    // max error  (normalised)
                                                   >> 673         G4double  hstepCurrent)  // current step size
710 {                                                 674 {
711   G4double hnew;                                  675   G4double hnew;
712                                                   676 
713   // Compute size of next Step for a failed st    677   // Compute size of next Step for a failed step
714   if(errMaxNorm > 1.0 )                        << 678   if(errMaxNorm > 1.0 ) {
715   {                                            << 679 
716     // Step failed; compute the size of retria    680     // Step failed; compute the size of retrial Step.
717     hnew = GetSafety()*hstepCurrent*std::pow(e    681     hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
718   }                                            << 682   }else if(errMaxNorm > 0.0 ){
719   else if(errMaxNorm > 0.0 )                   << 
720   {                                            << 
721     // Compute size of next Step for a success    683     // Compute size of next Step for a successful step
722     hnew = GetSafety()*hstepCurrent*std::pow(e    684     hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
723   }                                            << 685   }else {
724   else                                         << 
725   {                                            << 
726     // if error estimate is zero (possible) or    686     // if error estimate is zero (possible) or negative (dubious)
727     hnew = max_stepping_increase * hstepCurren    687     hnew = max_stepping_increase * hstepCurrent; 
728   }                                               688   }
729                                                   689 
730   return hnew;                                    690   return hnew;
731 }                                                 691 }
732                                                   692 
733 // ------------------------------------------- << 693 // -----------------------------------------------------------------------------
734                                                << 694 //  This method computes new step sizes limiting changes within certain factors
735 G4double                                       << 
736 G4MagInt_Driver::ComputeNewStepSize(           << 
737                           G4double  errMaxNorm << 
738                           G4double  hstepCurre << 
739 {                                              << 
740    // Legacy behaviour:                        << 
741    return ComputeNewStepSize_WithoutReductionL << 
742    // 'Improved' behaviour - at least more con << 
743    // return ComputeNewStepSize_WithinLimits(  << 
744 }                                              << 
745                                                << 
746 // This method computes new step sizes limitin << 
747 //                                                695 // 
748 // It shares its logic with AccurateAdvance.   << 696 //   It shares its logic with AccurateAdvance.
749 // They are kept separate currently for optimi << 697 //    They are kept separate currently for optimisation.
750 //                                             << 698 
751 G4double                                          699 G4double 
752 G4MagInt_Driver::ComputeNewStepSize_WithinLimi    700 G4MagInt_Driver::ComputeNewStepSize_WithinLimits( 
753                           G4double  errMaxNorm    701                           G4double  errMaxNorm,    // max error  (normalised)
754                           G4double  hstepCurre << 702         G4double  hstepCurrent)  // current step size
755 {                                                 703 {
756   G4double hnew;                                  704   G4double hnew;
757                                                   705 
758   // Compute size of next Step for a failed st    706   // Compute size of next Step for a failed step
759   if (errMaxNorm > 1.0 )                       << 707   if(errMaxNorm > 1.0 ) {
760   {                                            << 708 
761     // Step failed; compute the size of retria    709     // Step failed; compute the size of retrial Step.
762     hnew = GetSafety()*hstepCurrent*std::pow(e    710     hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
763                                                   711   
764     if (hnew < max_stepping_decrease*hstepCurr << 712     if(hnew < max_stepping_decrease*hstepCurrent) 
765     {                                          << 713          hnew = max_stepping_decrease*hstepCurrent ;
766       hnew = max_stepping_decrease*hstepCurren << 
767                          // reduce stepsize, b    714                          // reduce stepsize, but no more
768                          // than this factor (    715                          // than this factor (value= 1/10)
769     }                                          << 716   }else{
770   }                                            << 
771   else                                         << 
772   {                                            << 
773     // Compute size of next Step for a success    717     // Compute size of next Step for a successful step
774     if (errMaxNorm > errcon)                   << 718     if(errMaxNorm > errcon) hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
775      { hnew = GetSafety()*hstepCurrent*std::po << 719     else                    hnew = max_stepping_increase * hstepCurrent ;
776     else  // No more than a factor of 5 increa << 720       // No more than a factor of 5 increase
777      { hnew = max_stepping_increase * hstepCur << 
778   }                                               721   }
                                                   >> 722 
779   return hnew;                                    723   return hnew;
780 }                                                 724 }
781                                                   725 
782 // ------------------------------------------- << 
783                                                   726 
784 void G4MagInt_Driver::PrintStatus( const G4dou << 727 
785                                          G4dou << 728 void G4MagInt_Driver::PrintStatus( const G4double*   StartArr,  
786                                    const G4dou << 729            G4double          xstart,
787                                          G4dou << 730            const G4double*   CurrentArr, 
788                                          G4dou << 731            G4double          xcurrent,
789                                          G4int << 732            G4double          requestStep, 
                                                   >> 733            G4int             subStepNo)
790   // Potentially add as arguments:                734   // Potentially add as arguments:  
791   //                                 <dydx>       735   //                                 <dydx>           - as Initial Force
792   //                                 stepTaken    736   //                                 stepTaken(hdid)  - last step taken
793   //                                 nextStep     737   //                                 nextStep (hnext) - proposal for size
794 {                                                 738 {
795    G4FieldTrack  StartFT(G4ThreeVector(0,0,0), << 739    G4FieldTrack  StartFT(G4ThreeVector(0,0,0), G4ThreeVector(0,0,0), 0., 0., 0., 0. );
796                  G4ThreeVector(0,0,0), 0., 0., << 
797    G4FieldTrack  CurrentFT (StartFT);             740    G4FieldTrack  CurrentFT (StartFT);
798                                                   741 
799    StartFT.LoadFromArray( StartArr, fNoIntegra << 742    StartFT.LoadFromArray( StartArr); 
800    StartFT.SetCurveLength( xstart);               743    StartFT.SetCurveLength( xstart);
801    CurrentFT.LoadFromArray( CurrentArr, fNoInt << 744    CurrentFT.LoadFromArray( CurrentArr); 
802    CurrentFT.SetCurveLength( xcurrent );          745    CurrentFT.SetCurveLength( xcurrent );
803                                                   746 
804    PrintStatus(StartFT, CurrentFT, requestStep    747    PrintStatus(StartFT, CurrentFT, requestStep, subStepNo ); 
805 }                                                 748 }
806                                                   749 
807 // ------------------------------------------- << 
808                                                   750 
809 void G4MagInt_Driver::PrintStatus(const G4Fiel << 751 
810                                   const G4Fiel << 752 void G4MagInt_Driver::PrintStatus(
811                                         G4doub << 753                   const G4FieldTrack&  StartFT,
812                                         G4int  << 754       const G4FieldTrack&  CurrentFT, 
                                                   >> 755                   G4double             requestStep, 
                                                   >> 756                   // G4double             safety,
                                                   >> 757                   G4int                subStepNo)
813 {                                                 758 {
814     G4int verboseLevel= fVerboseLevel;            759     G4int verboseLevel= fVerboseLevel;
815     const G4int noPrecision = 5;               << 760     static G4int noPrecision= 5;
816     G4long oldPrec= G4cout.precision(noPrecisi << 761     G4int oldPrec= G4cout.precision(noPrecision);
817     // G4cout.setf(ios_base::fixed,ios_base::f    762     // G4cout.setf(ios_base::fixed,ios_base::floatfield);
818                                                   763 
819     const G4ThreeVector StartPosition=       S << 764     const G4ThreeVector StartPosition=      StartFT.GetPosition();
820     const G4ThreeVector StartUnitVelocity=   S << 765     const G4ThreeVector StartUnitVelocity=  StartFT.GetMomentumDir();
821     const G4ThreeVector CurrentPosition=     C << 766     const G4ThreeVector CurrentPosition=    CurrentFT.GetPosition();
822     const G4ThreeVector CurrentUnitVelocity= C << 767     const G4ThreeVector CurrentUnitVelocity=    CurrentFT.GetMomentumDir();
823                                                   768 
824     G4double  DotStartCurrentVeloc= StartUnitV    769     G4double  DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
825                                                   770 
826     G4double step_len= CurrentFT.GetCurveLengt << 771     G4double step_len= CurrentFT.GetCurveLength() 
                                                   >> 772                          - StartFT.GetCurveLength();
827     G4double subStepSize = step_len;              773     G4double subStepSize = step_len;
828                                                   774      
829     if( (subStepNo <= 1) || (verboseLevel > 3) << 775     if( (subStepNo <= 0) && (verboseLevel <= 3) )
830     {                                             776     {
831        subStepNo = - subStepNo;        // To a    777        subStepNo = - subStepNo;        // To allow printing banner
832                                                   778 
833        G4cout << std::setw( 6)  << " " << std: << 779        G4cout << std::setw( 6)  << " " 
834               << " G4MagInt_Driver: Current Po << 780         << std::setw( 25) << " G4MagInt_Driver: Current Position  and  Direction" << " "
835               << G4endl;                       << 781         << G4endl; 
836        G4cout << std::setw( 5) << "Step#" << "    782        G4cout << std::setw( 5) << "Step#" << " "
837               << std::setw( 7) << "s-curve" << << 783         << std::setw( 7) << "s-curve" << " "
838               << std::setw( 9) << "X(mm)" << " << 784         << std::setw( 9) << "X(mm)" << " "
839               << std::setw( 9) << "Y(mm)" << " << 785         << std::setw( 9) << "Y(mm)" << " "  
840               << std::setw( 9) << "Z(mm)" << " << 786         << std::setw( 9) << "Z(mm)" << " "
841               << std::setw( 8) << " N_x " << " << 787         << std::setw( 8) << " N_x " << " "
842               << std::setw( 8) << " N_y " << " << 788         << std::setw( 8) << " N_y " << " "
843               << std::setw( 8) << " N_z " << " << 789         << std::setw( 8) << " N_z " << " "
844               << std::setw( 8) << " N^2-1 " << << 790         << std::setw( 7) << " N^2-1 " << " "
845               << std::setw(10) << " N(0).N " < << 791         << std::setw(10) << " N(0).N " << " "
846               << std::setw( 7) << "KinEner " < << 792         << std::setw( 7) << "KinEner " << " "
847               << std::setw(12) << "Track-l" << << 793         << std::setw(12) << "Track-l" << " "   // Add the Sub-step ??
848               << std::setw(12) << "Step-len" < << 794         << std::setw(12) << "Step-len" << " " 
849               << std::setw(12) << "Step-len" < << 795         << std::setw(12) << "Step-len" << " " 
850               << std::setw( 9) << "ReqStep" << << 796         << std::setw( 9) << "ReqStep" << " "  
851               << G4endl;                       << 797         << G4endl;
852     }                                          << 798 
853                                                << 799         PrintStat_Aux( StartFT,  requestStep, 0., 
854     if( (subStepNo <= 0) )                     << 800            0,        0.0,         1.0);
855     {                                          << 801         //*************
856       PrintStat_Aux( StartFT,  requestStep, 0. << 
857                        0,        0.0,          << 
858     }                                             802     }
859                                                   803 
860     if( verboseLevel <= 3 )                       804     if( verboseLevel <= 3 )
861     {                                             805     {
862       G4cout.precision(noPrecision);           << 806        G4cout.precision(noPrecision);
863       PrintStat_Aux( CurrentFT, requestStep, s << 807        PrintStat_Aux( CurrentFT, requestStep, step_len, 
864                      subStepNo, subStepSize, D << 808           subStepNo, subStepSize, DotStartCurrentVeloc );
865     }                                          << 809        //*************
                                                   >> 810     }
                                                   >> 811 
                                                   >> 812     else // if( verboseLevel > 3 )
                                                   >> 813     {
                                                   >> 814        //  Multi-line output
                                                   >> 815        
                                                   >> 816        // G4cout << "Current  Position is " << CurrentPosition << G4endl 
                                                   >> 817        //    << " and UnitVelocity is " << CurrentUnitVelocity << G4endl;
                                                   >> 818        // G4cout << "Step taken was " << step_len  
                                                   >> 819        //    << " out of PhysicalStep= " <<  requestStep << G4endl;
                                                   >> 820        // G4cout << "Final safety is: " << safety << G4endl;
866                                                   821 
                                                   >> 822        // G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag() << G4endl;
                                                   >> 823        // G4cout << G4endl; 
                                                   >> 824     }
867     G4cout.precision(oldPrec);                    825     G4cout.precision(oldPrec);
868 }                                                 826 }
869                                                   827 
870 // ------------------------------------------- << 828 void G4MagInt_Driver::PrintStat_Aux(
871                                                << 829                   const G4FieldTrack&  aFieldTrack,
872 void G4MagInt_Driver::PrintStat_Aux(const G4Fi << 830                   G4double             requestStep, 
873                                           G4do << 831       G4double             step_len,
874                                           G4do << 832                   G4int                subStepNo,
875                                           G4in << 833       G4double             subStepSize,
876                                           G4do << 834       G4double             dotVeloc_StartCurr)
877                                           G4do << 
878 {                                                 835 {
879     const G4ThreeVector Position = aFieldTrack << 836     const G4ThreeVector Position=      aFieldTrack.GetPosition();
880     const G4ThreeVector UnitVelocity = aFieldT << 837     const G4ThreeVector UnitVelocity=  aFieldTrack.GetMomentumDir();
881                                                   838  
882     if( subStepNo >= 0)                           839     if( subStepNo >= 0)
883     {                                          << 
884        G4cout << std::setw( 5) << subStepNo <<    840        G4cout << std::setw( 5) << subStepNo << " ";
885     }                                          << 
886     else                                          841     else
887     {                                          << 
888        G4cout << std::setw( 5) << "Start" << "    842        G4cout << std::setw( 5) << "Start" << " ";
889     }                                          << 
890     G4double curveLen= aFieldTrack.GetCurveLen    843     G4double curveLen= aFieldTrack.GetCurveLength();
891     G4cout << std::setw( 7) << curveLen;          844     G4cout << std::setw( 7) << curveLen;
892     G4cout << std::setw( 9) << Position.x() <<    845     G4cout << std::setw( 9) << Position.x() << " "
893            << std::setw( 9) << Position.y() << << 846      << std::setw( 9) << Position.y() << " "
894            << std::setw( 9) << Position.z() << << 847      << std::setw( 9) << Position.z() << " "
895            << std::setw( 8) << UnitVelocity.x( << 848      << std::setw( 8) << UnitVelocity.x() << " "
896            << std::setw( 8) << UnitVelocity.y( << 849      << std::setw( 8) << UnitVelocity.y() << " "
897            << std::setw( 8) << UnitVelocity.z( << 850      << std::setw( 8) << UnitVelocity.z() << " ";
898     G4long oldprec= G4cout.precision(3);       << 851     G4int oldprec= G4cout.precision(3);
899     G4cout << std::setw( 8) << UnitVelocity.ma << 852     G4cout << std::setw( 7) << UnitVelocity.mag2()-1.0 << " ";
900     G4cout.precision(6);                          853     G4cout.precision(6);
901     G4cout << std::setw(10) << dotVeloc_StartC    854     G4cout << std::setw(10) << dotVeloc_StartCurr << " ";
902     G4cout.precision(oldprec);                    855     G4cout.precision(oldprec);
903     G4cout << std::setw( 7) << aFieldTrack.Get    856     G4cout << std::setw( 7) << aFieldTrack.GetKineticEnergy();
904     G4cout << std::setw(12) << step_len << " "    857     G4cout << std::setw(12) << step_len << " ";
905                                                   858 
906     static G4ThreadLocal G4double oldCurveLeng << 859     static G4double oldCurveLength= 0.0;
907     static G4ThreadLocal G4double oldSubStepLe << 860     static G4double oldSubStepLength= 0.0;
908     static G4ThreadLocal G4int oldSubStepNo =  << 861     static int oldSubStepNo= -1;
909                                                   862 
910     G4double subStep_len = 0.0;                << 863     G4double subStep_len=0.0;
911     if( curveLen > oldCurveLength )               864     if( curveLen > oldCurveLength )
912     {                                          << 865        subStep_len= curveLen - oldCurveLength;
913       subStep_len= curveLen - oldCurveLength;  << 
914     }                                          << 
915     else if (subStepNo == oldSubStepNo)           866     else if (subStepNo == oldSubStepNo)
916     {                                          << 867        subStep_len= oldSubStepLength;
917       subStep_len= oldSubStepLength;           << 868     //     else  subStepLen_NotAvail;
918     }                                          << 
919     oldCurveLength= curveLen;                     869     oldCurveLength= curveLen;
920     oldSubStepLength= subStep_len;                870     oldSubStepLength= subStep_len;
921                                                   871 
922     G4cout << std::setw(12) << subStep_len <<     872     G4cout << std::setw(12) << subStep_len << " "; 
923     G4cout << std::setw(12) << subStepSize <<     873     G4cout << std::setw(12) << subStepSize << " "; 
924     if( requestStep != -1.0 )                  << 874     if( requestStep != -1.0 ) 
925     {                                          << 875        G4cout << std::setw( 9) << requestStep << " ";
926       G4cout << std::setw( 9) << requestStep < << 
927     }                                          << 
928     else                                          876     else
929     {                                          << 877        G4cout << std::setw( 9) << " InitialStep " << " "; 
930        G4cout << std::setw( 9) << " InitialSte << 878     // G4cout << std::setw(12) << safety << " ";
931     }                                          << 
932     G4cout << G4endl;                             879     G4cout << G4endl;
933 }                                                 880 }
934                                                   881 
935 // ------------------------------------------- << 
936                                                << 
937 void G4MagInt_Driver::PrintStatisticsReport()     882 void G4MagInt_Driver::PrintStatisticsReport()
938 {                                                 883 {
939   G4int noPrecBig = 6;                         << 884   G4int noPrecBig= 6;
940   G4long oldPrec = G4cout.precision(noPrecBig) << 885   G4int oldPrec= G4cout.precision(noPrecBig);
941                                                   886 
942   G4cout << "G4MagInt_Driver Statistics of ste    887   G4cout << "G4MagInt_Driver Statistics of steps undertaken. " << G4endl;
943   G4cout << "G4MagInt_Driver: Number of Steps:    888   G4cout << "G4MagInt_Driver: Number of Steps: "
944          << " Total= " <<  fNoTotalSteps       << 889    << " Total= " <<  fNoTotalSteps
945          << " Bad= "   <<  fNoBadSteps            890          << " Bad= "   <<  fNoBadSteps 
946          << " Small= " <<  fNoSmallSteps          891          << " Small= " <<  fNoSmallSteps 
947          << " Non-initial small= " << (fNoSmal << 892    << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
948          << G4endl;                            << 893    << G4endl;
949  G4cout.precision(oldPrec);                    << 
950 }                                              << 
951                                                << 
952 // ------------------------------------------- << 
953                                                << 
954 void G4MagInt_Driver::SetSmallestFraction(G4do << 
955 {                                              << 
956   if( (newFraction > 1.e-16) && (newFraction < << 
957   {                                            << 
958     fSmallestFraction= newFraction;            << 
959   }                                            << 
960   else                                         << 
961   {                                            << 
962     std::ostringstream message;                << 
963     message << "Smallest Fraction not changed. << 
964             << "  Proposed value was " << newF << 
965             << "  Value must be between 1.e-8  << 
966     G4Exception("G4MagInt_Driver::SetSmallestF << 
967                 "GeomField1001", JustWarning,  << 
968   }                                            << 
969 }                                              << 
970                                                << 
971 void G4MagInt_Driver::                         << 
972 GetDerivatives(const G4FieldTrack& y_curr, G4d << 
973 {                                              << 
974     G4double ytemp[G4FieldTrack::ncompSVEC];   << 
975     y_curr.DumpToArray(ytemp);                 << 
976     pIntStepper->RightHandSide(ytemp, dydx);   << 
977       // Avoid virtual call for GetStepper     << 
978       // Was: GetStepper()->ComputeRightHandSi << 
979 }                                              << 
980                                                << 
981 void G4MagInt_Driver::GetDerivatives(const G4F << 
982                                      G4double  << 
983                                      G4double  << 
984 {                                              << 
985     G4double ytemp[G4FieldTrack::ncompSVEC];   << 
986     track.DumpToArray(ytemp);                  << 
987     pIntStepper->RightHandSide(ytemp, dydx, fi << 
988 }                                              << 
989                                                << 
990 G4EquationOfMotion* G4MagInt_Driver::GetEquati << 
991 {                                              << 
992     return pIntStepper->GetEquationOfMotion(); << 
993 }                                              << 
994                                                << 
995 void G4MagInt_Driver::SetEquationOfMotion(G4Eq << 
996 {                                              << 
997     pIntStepper->SetEquationOfMotion(equation) << 
998 }                                              << 
999                                                   894 
1000 const G4MagIntegratorStepper* G4MagInt_Driver << 895  #ifdef G4FLD_STATS
1001 {                                             << 896  G4cout << "MID dyerr: " 
1002     return pIntStepper;                       << 897    << " maximum= " << fDyerr_max 
1003 }                                             << 898    // << " 2nd max= " << fDyerr_mx2
                                                   >> 899    << " Sum small= " << fDyerrPos_smTot 
                                                   >> 900    << " std::sqrt(Sum large^2): pos= " << std::sqrt(fDyerrPos_lgTot)
                                                   >> 901    << " vel= " << std::sqrt( fDyerrVel_lgTot )
                                                   >> 902    << " Total h-distance: small= " << fSumH_sm 
                                                   >> 903    << " large= " << fSumH_lg
                                                   >> 904    << G4endl;
                                                   >> 905 
                                                   >> 906 #if 0
                                                   >> 907   G4int noPrecSmall=4; 
                                                   >> 908   // Single line precis of statistics ... optional
                                                   >> 909   G4cout.precision(noPrecSmall);
                                                   >> 910   G4cout << "MIDnums: " << fMinimumStep
                                                   >> 911    << "   " << fNoTotalSteps 
                                                   >> 912    << "  "  <<  fNoSmallSteps
                                                   >> 913          << "  "  << fNoSmallSteps-fNoInitialSmallSteps
                                                   >> 914    << "  "  << fNoBadSteps         
                                                   >> 915    << "   " << fDyerr_max
                                                   >> 916    << "   " << fDyerr_mx2 
                                                   >> 917    << "   " << fDyerrPos_smTot 
                                                   >> 918    << "   " << fSumH_sm
                                                   >> 919    << "   " << fDyerrPos_lgTot
                                                   >> 920    << "   " << fDyerrVel_lgTot
                                                   >> 921    << "   " << fSumH_lg
                                                   >> 922    << G4endl;
                                                   >> 923  #endif 
                                                   >> 924  #endif 
1004                                                  925 
1005 G4MagIntegratorStepper* G4MagInt_Driver::GetS << 926  G4cout.precision(oldPrec);
1006 {                                             << 
1007     return pIntStepper;                       << 
1008 }                                             << 
1009                                               << 
1010 void G4MagInt_Driver::                        << 
1011 RenewStepperAndAdjust(G4MagIntegratorStepper* << 
1012 {                                             << 
1013     pIntStepper = pItsStepper;                << 
1014     ReSetParameters();                        << 
1015 }                                             << 
1016                                               << 
1017 void G4MagInt_Driver::StreamInfo( std::ostrea << 
1018 {                                             << 
1019     os << "State of G4MagInt_Driver: " << std << 
1020     os << "  Max number of Steps = " << fMaxN << 
1021        << "    (base # = " << fMaxStepBase << << 
1022     os << "  Safety factor       = " << safet << 
1023     os << "  Power - shrink      = " << pshrn << 
1024     os << "  Power - grow        = " << pgrow << 
1025     os << "  threshold (errcon)  = " << errco << 
1026                                               << 
1027     os << "    fMinimumStep =      " << fMini << 
1028     os << "    Smallest Fraction = " << fSmal << 
1029                                               << 
1030     os << "    No Integrat Vars  = " << fNoIn << 
1031     os << "    Min No Vars       = " << fMinN << 
1032     os << "    Num-Vars          = " << fNoVa << 
1033                                               << 
1034     os << "    verbose level     = " << fVerb << 
1035     os << "    Reintegrates      = " << DoesR << 
1036 }                                             << 
1037                                               << 
1038 void PrintInfo( const G4MagInt_Driver & magDr << 
1039 {                                             << 
1040     os << "State of G4MagInt_Driver: " << std << 
1041     os << "  Max number of Steps = " << magDr << 
1042     //   << "    (base # = " << magDrv.fMaxSt << 
1043     os << "  Safety factor       = " << magDr << 
1044     os << "  Power - shrink      = " << magDr << 
1045     os << "  Power - grow        = " << magDr << 
1046     os << "  threshold (errcon)  = " << magDr << 
1047                                               << 
1048     os << "    fMinimumStep =      " << magDr << 
1049     os << "    Smallest Fraction = " << magDr << 
1050                                               << 
1051     /*****                                    << 
1052     os << "    No Integrat Vars  = " << magDr << 
1053     os << "    Min No Vars       = " << magDr << 
1054     os << "    Num-Vars          = " << magDr << 
1055     *****/                                    << 
1056     os << "    verbose level     = " << magDr << 
1057     os << "    Reintegrates      = " << magDr << 
1058 }                                                927 }
1059                                                  928