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