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 9.4.p4)


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