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