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