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