61 #ifndef G4NO_FIELD_STATISTICS
72 G4int statisticsVerbose)
73 : fSmallestFraction( 1.0e-12 ),
74 fNoIntegrationVariables(numComponents),
76 fNoVars(
std::
max( fNoIntegrationVariables, fMinNoVars )),
77 fStatisticsVerboseLevel(statisticsVerbose),
78 fNoTotalSteps(0), fNoBadSteps(0), fNoSmallSteps(0),
79 fNoInitialSmallSteps(0),
80 fDyerr_max(0.0), fDyerr_mx2(0.0),
81 fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0),
82 fSumH_sm(0.0), fSumH_lg(0.0),
97 G4cout <<
"MagIntDriver version: Accur-Adv: "
98 <<
"invE_nS, QuickAdv-2sqrt with Statistics "
138 G4int nstp, i, no_warnings=0;
143 static G4int nStpPr=50;
151 G4bool succeeded =
true, lastStepSucceeded;
155 G4int noFullIntegr=0, noSmallIntegr = 0 ;
167 std::ostringstream message;
168 message <<
"Proposed step is zero; hstep = " << hstep <<
" !";
175 std::ostringstream message;
176 message <<
"Invalid run condition." <<
G4endl
177 <<
"Proposed step is negative; hstep = " << hstep <<
"." <<
G4endl
178 <<
"Requested step cannot be negative! Aborting event.";
188 x1= startCurveLength;
191 if ( (hinitial > 0.0) && (hinitial < hstep)
203 for (i=0;i<nvar;i++) { y[i] = ystart[i]; }
214 for (i=0;i<nvar;i++) { ySubStepStart[i] = y[i]; }
231 lastStepSucceeded= (hdid == h);
234 PrintStatus( ySubStepStart, xSubStepStart, y, x, h, nstp);
242 G4double dchord_step, dyerr, dyerr_len;
246 QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
265 G4cout <<
" dyerr= " << dyerr_len <<
" relative = " << dyerr_len / h
266 <<
" epsilon= " << eps <<
" hstep= " << hstep
274 "Integration Step became Zero!");
276 dyerr = dyerr_len / h;
284 lastStepSucceeded= (dyerr<=
eps);
287 if (lastStepSucceeded) { noFullIntegr++; }
288 else { noSmallIntegr++; }
293 if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
295 if( nstp==nStpPr ) {
G4cout <<
"***** Many steps ****" <<
G4endl; }
297 G4cout <<
"hdid=" << std::setw(12) << hdid <<
" "
298 <<
"hnext=" << std::setw(12) << hnext <<
" "
299 <<
"hstep=" << std::setw(12) << hstep <<
" (requested) "
301 PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
306 G4double endPointDist= (EndPos-StartPos).mag();
320 <<
" good " << noGoodSteps <<
" current h= " << hdid
322 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
343 if(std::fabs(hnext) <=
Hmin())
347 if( (x < x2 * (1-eps) ) &&
348 (std::fabs(hstep) >
Hmin()) )
353 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
380 G4cout <<
"Warning: G4MagIntegratorDriver::AccurateAdvance"
382 <<
" Integration step 'h' became "
383 << h <<
" due to roundoff. " <<
G4endl
384 <<
" Calculated as difference of x2= "<< x2 <<
" and x=" << x
385 <<
" Forcing termination of advance." <<
G4endl;
391 }
while ( ((nstp++)<=
fMaxNoSteps) && (x < x2) && (!lastStep) );
397 for (i=0;i<nvar;i++) { yEnd[i] = y[i]; }
417 if( dbg && no_warnings )
419 G4cerr <<
"G4MagIntegratorDriver exit status: no-steps " << nstp <<
G4endl;
435 const G4int maxNoWarnings = 10;
436 std::ostringstream message;
437 if( (noWarningsIssued < maxNoWarnings) ||
fVerboseLevel > 10 )
439 message <<
"The stepsize for the next iteration, " << hnext
440 <<
", is too small - in Step number " << nstp <<
"." <<
G4endl
441 <<
"The minimum for the driver is " <<
Hmin() <<
G4endl
442 <<
"Requested integr. length was " << hstep <<
" ." <<
G4endl
443 <<
"The size of this sub-step was " << h <<
" ." <<
G4endl
444 <<
"The integrations has already gone " << xDone;
448 message <<
"Too small 'next' step " << hnext
449 <<
", step-no: " << nstp <<
G4endl
450 <<
", this sub-step: " << h
451 <<
", req_tot_len: " << hstep
452 <<
", done: " << xDone <<
", min: " <<
Hmin();
454 G4Exception(
"G4MagInt_Driver::WarnSmallStepSize()",
"GeomField1001",
466 std::ostringstream message;
467 message <<
"The number of steps used in the Integration driver"
468 <<
" (Runge-Kutta) is too many." <<
G4endl
469 <<
"Integration of the interval was not completed !" <<
G4endl
470 <<
"Only a " << (xCurrent-x1start)*100/(x2end-x1start)
471 <<
" % fraction of it was done.";
472 G4Exception(
"G4MagInt_Driver::WarnTooManyStep()",
"GeomField1001",
485 G4bool isNewMax, prNewMax;
487 isNewMax = endPointDist > (1.0 + maxRelError) * h;
488 prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
489 if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
492 && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+
eps) ) ) )
495 std::ostringstream message;
496 if( (noWarnings ++ < 10) || (dbg>2) )
498 message <<
"The integration produced an end-point which " <<
G4endl
499 <<
"is further from the start-point than the curve length."
502 message <<
" Distance of endpoints = " << endPointDist
503 <<
", curve length = " << h <<
G4endl
504 <<
" Difference (curveLen-endpDist)= " << (h - endPointDist)
505 <<
", relative = " << (h-endPointDist) / h
506 <<
", epsilon = " <<
eps;
507 G4Exception(
"G4MagInt_Driver::WarnEndPointTooFar()",
"GeomField1001",
544 G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
553 const G4int max_trials=100;
557 G4bool hasSpin= (spin_mag2 > 0.0);
559 for (iter=0; iter<max_trials ;iter++)
565 G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos);
569 errpos_sq =
sqr(yerr[0]) +
sqr(yerr[1]) +
sqr(yerr[2]) ;
570 errpos_sq *= inv_eps_pos_sq;
575 if( magvel_sq > 0.0 ) {
576 errvel_sq = sumerr_sq / magvel_sq;
578 G4cerr <<
"** G4MagIntegrationDriver: found case of zero momentum."
579 <<
" iteration= " << iter <<
" h= " << h <<
G4endl;
580 errvel_sq = sumerr_sq;
582 errvel_sq *= inv_eps_vel_sq;
583 errmax_sq =
std::max( errpos_sq, errvel_sq );
588 errspin_sq = (
sqr(yerr[9]) +
sqr(yerr[10]) +
sqr(yerr[11]) )
590 errspin_sq *= inv_eps_vel_sq;
591 errmax_sq =
std::max( errmax_sq, errspin_sq );
594 if ( errmax_sq <= 1.0 ) {
break; }
599 if (htemp >= 0.1*h) { h = htemp; }
605 G4cerr <<
"G4MagIntegratorDriver::OneGoodStep:" <<
G4endl
606 <<
" Stepsize underflow in Stepper " <<
G4endl ;
607 G4cerr <<
" Step's start x=" << x <<
" and end x= " << xnew
608 <<
" are equal !! " << G4endl
609 <<
" Due to step-size= " << h
610 <<
" . Note that input step was " << htry <<
G4endl;
650 G4Exception(
"G4MagInt_Driver::QuickAdvance()",
"GeomField0001",
654 dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0];
668 G4double dyerr_pos_sq, dyerr_mom_rel_sq;
672 G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
682 pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ;
697 PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep, 1);
703 vel_mag_sq = (
sqr(yarrout[3])+
sqr(yarrout[4])+
sqr(yarrout[5]) );
704 inv_vel_mag_sq = 1.0 / vel_mag_sq;
705 dyerr_pos_sq = (
sqr(yerr_vec[0])+
sqr(yerr_vec[1])+
sqr(yerr_vec[2]));
706 dyerr_mom_sq = (
sqr(yerr_vec[3])+
sqr(yerr_vec[4])+
sqr(yerr_vec[5]));
707 dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
713 #ifdef RETURN_A_NEW_STEP_LENGTH
715 dyerr_len = std::sqrt( dyerr_len_sq );
716 dyerr_len_sq /=
eps ;
725 if( dyerr_pos_sq > ( dyerr_mom_rel_sq *
sqr(hstep) ) )
727 dyerr = std::sqrt(dyerr_pos_sq);
732 dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
740 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT
749 G4Exception(
"G4MagInt_Driver::QuickAdvance()",
"GeomField0001",
751 dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
752 yarrout[0]= yarrin[0];
769 if(errMaxNorm > 1.0 )
773 }
else if(errMaxNorm > 0.0 ) {
799 if (errMaxNorm > 1.0 )
844 PrintStatus(StartFT, CurrentFT, requestStep, subStepNo );
865 G4double DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
870 if( (subStepNo <= 1) || (verboseLevel > 3) )
872 subStepNo = - subStepNo;
874 G4cout << std::setw( 6) <<
" " << std::setw( 25)
875 <<
" G4MagInt_Driver: Current Position and Direction" <<
" "
877 G4cout << std::setw( 5) <<
"Step#" <<
" "
878 << std::setw( 7) <<
"s-curve" <<
" "
879 << std::setw( 9) <<
"X(mm)" <<
" "
880 << std::setw( 9) <<
"Y(mm)" <<
" "
881 << std::setw( 9) <<
"Z(mm)" <<
" "
882 << std::setw( 8) <<
" N_x " <<
" "
883 << std::setw( 8) <<
" N_y " <<
" "
884 << std::setw( 8) <<
" N_z " <<
" "
885 << std::setw( 8) <<
" N^2-1 " <<
" "
886 << std::setw(10) <<
" N(0).N " <<
" "
887 << std::setw( 7) <<
"KinEner " <<
" "
888 << std::setw(12) <<
"Track-l" <<
" "
889 << std::setw(12) <<
"Step-len" <<
" "
890 << std::setw(12) <<
"Step-len" <<
" "
891 << std::setw( 9) <<
"ReqStep" <<
" "
895 if( (subStepNo <= 0) )
902 if( verboseLevel <= 3 )
904 G4cout.precision(noPrecision);
906 subStepNo, subStepSize, DotStartCurrentVeloc );
922 G4cout.precision(oldPrec);
940 G4cout << std::setw( 5) << subStepNo <<
" ";
944 G4cout << std::setw( 5) <<
"Start" <<
" ";
947 G4cout << std::setw( 7) << curveLen;
948 G4cout << std::setw( 9) << Position.x() <<
" "
949 << std::setw( 9) << Position.y() <<
" "
950 << std::setw( 9) << Position.z() <<
" "
951 << std::setw( 8) << UnitVelocity.x() <<
" "
952 << std::setw( 8) << UnitVelocity.y() <<
" "
953 << std::setw( 8) << UnitVelocity.z() <<
" ";
955 G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 <<
" ";
957 G4cout << std::setw(10) << dotVeloc_StartCurr <<
" ";
958 G4cout.precision(oldprec);
960 G4cout << std::setw(12) << step_len <<
" ";
967 if( curveLen > oldCurveLength )
969 subStep_len= curveLen - oldCurveLength;
971 else if (subStepNo == oldSubStepNo)
973 subStep_len= oldSubStepLength;
975 oldCurveLength= curveLen;
976 oldSubStepLength= subStep_len;
978 G4cout << std::setw(12) << subStep_len <<
" ";
979 G4cout << std::setw(12) << subStepSize <<
" ";
980 if( requestStep != -1.0 )
982 G4cout << std::setw( 9) << requestStep <<
" ";
986 G4cout << std::setw( 9) <<
" InitialStep " <<
" ";
998 G4cout <<
"G4MagInt_Driver Statistics of steps undertaken. " <<
G4endl;
999 G4cout <<
"G4MagInt_Driver: Number of Steps: "
1012 <<
" Total h-distance: small= " <<
fSumH_sm
1017 G4int noPrecSmall=4;
1019 G4cout.precision(noPrecSmall);
1036 G4cout.precision(oldPrec);
1043 if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
1049 G4cerr <<
"Warning: SmallestFraction not changed. " <<
G4endl
1050 <<
" Proposed value was " << newFraction <<
G4endl
1051 <<
" Value must be between 1.e-8 and 1.e-16" <<
G4endl;
virtual void ComputeRightHandSide(const G4double y[], G4double dydx[])
G4double GetCurveLength() const
CLHEP::Hep3Vector G4ThreeVector
G4MagIntegratorStepper * pIntStepper
G4double GetKineticEnergy() const
void SetCurveLength(G4double nCurve_s)
void RenewStepperAndAdjust(G4MagIntegratorStepper *pItsStepper)
void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent)
void WarnEndPointTooFar(G4double endPointDist, G4double hStepSize, G4double epsilonRelative, G4int debugFlag)
G4double GetSurfaceTolerance() const
const G4ThreeVector & GetMomentumDir() const
static const G4double eps
static const double perThousand
G4MagInt_Driver(G4double hminimum, G4MagIntegratorStepper *pItsStepper, G4int numberOfComponents=6, G4int statisticsVerbosity=1)
void OneGoodStep(G4double ystart[], const G4double dydx[], G4double &x, G4double htry, G4double eps, G4double &hdid, G4double &hnext)
G4double GetSafety() const
static const G4double max_stepping_increase
G4int fStatisticsVerboseLevel
void PrintStatisticsReport()
virtual G4int IntegratorOrder() const =0
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)
G4double ComputeNewStepSize_WithinLimits(G4double errMaxNorm, G4double hstepCurrent)
G4double GetPgrow() const
void PrintStat_Aux(const G4FieldTrack &aFieldTrack, G4double requestStep, G4double actualStep, G4int subStepNo, G4double subStepSize, G4double dotVelocities)
G4double GetPshrnk() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static const double perMillion
G4bool QuickAdvance(G4FieldTrack &y_val, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr)
void DumpToArray(G4double valArr[ncompSVEC]) const
void WarnSmallStepSize(G4double hnext, G4double hstep, G4double h, G4double xDone, G4int noSteps)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void LoadFromArray(const G4double valArr[ncompSVEC], G4int noVarsIntegrated)
const G4double x[NPOINTSGL]
void SetSmallestFraction(G4double val)
G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent)
G4double fSmallestFraction
static const G4int fMaxStepBase
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
G4int fNoInitialSmallSteps
static G4GeometryTolerance * GetInstance()
G4GLOB_DLL std::ostream G4cerr
const G4int fNoIntegrationVariables
static const G4double max_stepping_decrease