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