52 G4int statisticsVerbose)
53 : fNoIntegrationVariables(numComponents),
54 fNoVars(
std::
max( fNoIntegrationVariables, fMinNoVars )),
55 fStatisticsVerboseLevel(statisticsVerbose)
70 G4cout <<
"MagIntDriver version: Accur-Adv: "
71 <<
"invE_nS, QuickAdv-2sqrt with Statistics "
107 G4int nstp, i, no_warnings = 0;
111 static G4int dbg = 1;
112 static G4int nStpPr = 50;
120 G4bool succeeded =
true, lastStepSucceeded;
124 G4int noFullIntegr = 0, noSmallIntegr = 0;
136 std::ostringstream message;
137 message <<
"Proposed step is zero; hstep = " << hstep <<
" !";
144 std::ostringstream message;
145 message <<
"Invalid run condition." <<
G4endl
146 <<
"Proposed step is negative; hstep = " << hstep <<
"." <<
G4endl
147 <<
"Requested step cannot be negative! Aborting event.";
157 x1= startCurveLength;
160 if ( (hinitial > 0.0) && (hinitial < hstep)
172 for ( i=0; i<nvar; ++i) { y[i] = ystart[i]; }
183 for (i=0; i<nvar; ++i) { ySubStepStart[i] = y[i]; }
197 lastStepSucceeded = (hdid == h);
210 G4double dchord_step, dyerr, dyerr_len;
214 QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
233 G4cout <<
" dyerr= " << dyerr_len <<
" relative = " << dyerr_len / h
234 <<
" epsilon= " <<
eps <<
" hstep= " << hstep
242 "Integration Step became Zero!");
244 dyerr = dyerr_len / h;
252 lastStepSucceeded = (dyerr<=
eps);
255 if (lastStepSucceeded) { ++noFullIntegr; }
256 else { ++noSmallIntegr; }
261 if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
263 if( nstp==nStpPr ) {
G4cout <<
"***** Many steps ****" <<
G4endl; }
265 G4cout <<
"hdid=" << std::setw(12) << hdid <<
" "
266 <<
"hnext=" << std::setw(12) << hnext <<
" "
267 <<
"hstep=" << std::setw(12) << hstep <<
" (requested) "
269 PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
274 G4double endPointDist= (EndPos-StartPos).mag();
288 <<
" good " << noGoodSteps <<
" current h= " << hdid
290 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
311 if(std::fabs(hnext) <=
Hmin())
315 if( (x < x2 * (1-
eps) ) &&
316 (std::fabs(hstep) >
Hmin()) )
321 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
348 G4cout <<
"Warning: G4MagIntegratorDriver::AccurateAdvance"
350 <<
" Integration step 'h' became "
351 << h <<
" due to roundoff. " <<
G4endl
352 <<
" Calculated as difference of x2= "<< x2 <<
" and x=" << x
353 <<
" Forcing termination of advance." <<
G4endl;
359 }
while ( ((++nstp)<=
fMaxNoSteps) && (x < x2) && (!lastStep) );
367 for (i=0; i<nvar; ++i) { yEnd[i] = y[i]; }
387 if( dbg && no_warnings )
389 G4cerr <<
"G4MagIntegratorDriver exit status: no-steps " << nstp <<
G4endl;
405 const G4int maxNoWarnings = 10;
406 std::ostringstream message;
407 if( (noWarningsIssued < maxNoWarnings) ||
fVerboseLevel > 10 )
409 message <<
"The stepsize for the next iteration, " << hnext
410 <<
", is too small - in Step number " << nstp <<
"." <<
G4endl
411 <<
"The minimum for the driver is " <<
Hmin() <<
G4endl
412 <<
"Requested integr. length was " << hstep <<
" ." <<
G4endl
413 <<
"The size of this sub-step was " << h <<
" ." <<
G4endl
414 <<
"The integrations has already gone " << xDone;
418 message <<
"Too small 'next' step " << hnext
419 <<
", step-no: " << nstp <<
G4endl
420 <<
", this sub-step: " << h
421 <<
", req_tot_len: " << hstep
422 <<
", done: " << xDone <<
", min: " <<
Hmin();
424 G4Exception(
"G4MagInt_Driver::WarnSmallStepSize()",
"GeomField1001",
436 std::ostringstream message;
437 message <<
"The number of steps used in the Integration driver"
438 <<
" (Runge-Kutta) is too many." <<
G4endl
439 <<
"Integration of the interval was not completed !" <<
G4endl
440 <<
"Only a " << (xCurrent-x1start)*100/(x2end-x1start)
441 <<
" % fraction of it was done.";
442 G4Exception(
"G4MagInt_Driver::WarnTooManyStep()",
"GeomField1001",
455 G4bool isNewMax, prNewMax;
457 isNewMax = endPointDist > (1.0 + maxRelError) * h;
458 prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
459 if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
462 && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+
eps) ) ) )
465 std::ostringstream message;
466 if( (noWarnings++ < 10) || (dbg>2) )
468 message <<
"The integration produced an end-point which " <<
G4endl
469 <<
"is further from the start-point than the curve length."
472 message <<
" Distance of endpoints = " << endPointDist
473 <<
", curve length = " << h <<
G4endl
474 <<
" Difference (curveLen-endpDist)= " << (h - endPointDist)
475 <<
", relative = " << (h-endPointDist) / h
476 <<
", epsilon = " <<
eps;
477 G4Exception(
"G4MagInt_Driver::WarnEndPointTooFar()",
"GeomField1001",
514 G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
521 const G4int max_trials=100;
525 G4bool hasSpin = (spin_mag2 > 0.0);
527 for (
G4int iter=0; iter<max_trials; ++iter)
533 G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos);
537 errpos_sq =
sqr(yerr[0]) +
sqr(yerr[1]) +
sqr(yerr[2]) ;
538 errpos_sq *= inv_eps_pos_sq;
543 if( magvel_sq > 0.0 )
545 errvel_sq = sumerr_sq / magvel_sq;
549 std::ostringstream message;
550 message <<
"Found case of zero momentum." <<
G4endl
551 <<
"- iteration= " << iter <<
"; h= " << h;
554 errvel_sq = sumerr_sq;
556 errvel_sq *= inv_eps_vel_sq;
557 errmax_sq =
std::max( errpos_sq, errvel_sq );
562 errspin_sq = (
sqr(yerr[9]) +
sqr(yerr[10]) +
sqr(yerr[11]) )
564 errspin_sq *= inv_eps_vel_sq;
565 errmax_sq =
std::max( errmax_sq, errspin_sq );
568 if ( errmax_sq <= 1.0 ) {
break; }
573 if (htemp >= 0.1*h) { h = htemp; }
579 std::ostringstream message;
580 message <<
"Stepsize underflow in Stepper !" <<
G4endl
581 <<
"- Step's start x=" << x <<
" and end x= " << xnew
582 <<
" are equal !! " <<
G4endl
583 <<
" Due to step-size= " << h
584 <<
". Note that input step was " << htry;
618 G4Exception(
"G4MagInt_Driver::QuickAdvance()",
"GeomField0001",
623 dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0];
636 G4double dyerr_pos_sq, dyerr_mom_rel_sq;
640 G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
650 pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ;
663 PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep, 1);
669 vel_mag_sq = (
sqr(yarrout[3])+
sqr(yarrout[4])+
sqr(yarrout[5]) );
670 inv_vel_mag_sq = 1.0 / vel_mag_sq;
671 dyerr_pos_sq = (
sqr(yerr_vec[0])+
sqr(yerr_vec[1])+
sqr(yerr_vec[2]));
672 dyerr_mom_sq = (
sqr(yerr_vec[3])+
sqr(yerr_vec[4])+
sqr(yerr_vec[5]));
673 dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
679#ifdef RETURN_A_NEW_STEP_LENGTH
681 dyerr_len = std::sqrt( dyerr_len_sq );
682 dyerr_len_sq /=
eps ;
691 if( dyerr_pos_sq > ( dyerr_mom_rel_sq *
sqr(hstep) ) )
693 dyerr = std::sqrt(dyerr_pos_sq);
698 dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
706#ifdef QUICK_ADV_ARRAY_IN_AND_OUT
714 G4Exception(
"G4MagInt_Driver::QuickAdvance()",
"GeomField0001",
716 dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
717 yarrout[0]= yarrin[0];
733 if(errMaxNorm > 1.0 )
738 else if(errMaxNorm > 0.0 )
778 if (errMaxNorm > 1.0 )
823 PrintStatus(StartFT, CurrentFT, requestStep, subStepNo );
843 G4double DotStartCurrentVeloc= StartUnitVelocity.
dot(CurrentUnitVelocity);
848 if( (subStepNo <= 1) || (verboseLevel > 3) )
850 subStepNo = - subStepNo;
852 G4cout << std::setw( 6) <<
" " << std::setw( 25)
853 <<
" G4MagInt_Driver: Current Position and Direction" <<
" "
855 G4cout << std::setw( 5) <<
"Step#" <<
" "
856 << std::setw( 7) <<
"s-curve" <<
" "
857 << std::setw( 9) <<
"X(mm)" <<
" "
858 << std::setw( 9) <<
"Y(mm)" <<
" "
859 << std::setw( 9) <<
"Z(mm)" <<
" "
860 << std::setw( 8) <<
" N_x " <<
" "
861 << std::setw( 8) <<
" N_y " <<
" "
862 << std::setw( 8) <<
" N_z " <<
" "
863 << std::setw( 8) <<
" N^2-1 " <<
" "
864 << std::setw(10) <<
" N(0).N " <<
" "
865 << std::setw( 7) <<
"KinEner " <<
" "
866 << std::setw(12) <<
"Track-l" <<
" "
867 << std::setw(12) <<
"Step-len" <<
" "
868 << std::setw(12) <<
"Step-len" <<
" "
869 << std::setw( 9) <<
"ReqStep" <<
" "
873 if( (subStepNo <= 0) )
879 if( verboseLevel <= 3 )
883 subStepNo, subStepSize, DotStartCurrentVeloc );
886 G4cout.precision(oldPrec);
903 G4cout << std::setw( 5) << subStepNo <<
" ";
907 G4cout << std::setw( 5) <<
"Start" <<
" ";
910 G4cout << std::setw( 7) << curveLen;
911 G4cout << std::setw( 9) << Position.x() <<
" "
912 << std::setw( 9) << Position.y() <<
" "
913 << std::setw( 9) << Position.z() <<
" "
914 << std::setw( 8) << UnitVelocity.
x() <<
" "
915 << std::setw( 8) << UnitVelocity.
y() <<
" "
916 << std::setw( 8) << UnitVelocity.
z() <<
" ";
918 G4cout << std::setw( 8) << UnitVelocity.
mag2()-1.0 <<
" ";
920 G4cout << std::setw(10) << dotVeloc_StartCurr <<
" ";
921 G4cout.precision(oldprec);
923 G4cout << std::setw(12) << step_len <<
" ";
930 if( curveLen > oldCurveLength )
932 subStep_len= curveLen - oldCurveLength;
934 else if (subStepNo == oldSubStepNo)
936 subStep_len= oldSubStepLength;
938 oldCurveLength= curveLen;
939 oldSubStepLength= subStep_len;
941 G4cout << std::setw(12) << subStep_len <<
" ";
942 G4cout << std::setw(12) << subStepSize <<
" ";
943 if( requestStep != -1.0 )
945 G4cout << std::setw( 9) << requestStep <<
" ";
949 G4cout << std::setw( 9) <<
" InitialStep " <<
" ";
961 G4cout <<
"G4MagInt_Driver Statistics of steps undertaken. " <<
G4endl;
962 G4cout <<
"G4MagInt_Driver: Number of Steps: "
968 G4cout.precision(oldPrec);
975 if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
981 std::ostringstream message;
982 message <<
"Smallest Fraction not changed. " <<
G4endl
983 <<
" Proposed value was " << newFraction <<
G4endl
984 <<
" Value must be between 1.e-8 and 1.e-16";
985 G4Exception(
"G4MagInt_Driver::SetSmallestFraction()",
1038 os <<
"State of G4MagInt_Driver: " << std::endl;
1041 os <<
" Safety factor = " <<
safety << std::endl;
1042 os <<
" Power - shrink = " <<
pshrnk << std::endl;
1043 os <<
" Power - grow = " <<
pgrow << std::endl;
1044 os <<
" threshold (errcon) = " <<
errcon << std::endl;
1050 os <<
" Min No Vars = " <<
fMinNoVars << std::endl;
1051 os <<
" Num-Vars = " <<
fNoVars << std::endl;
1059 os <<
"State of G4MagInt_Driver: " << std::endl;
1062 os <<
" Safety factor = " << magDrv.
GetSafety() << std::endl;
1063 os <<
" Power - shrink = " << magDrv.
GetPshrnk() << std::endl;
1064 os <<
" Power - grow = " << magDrv.
GetPgrow() << std::endl;
1065 os <<
" threshold (errcon) = " << magDrv.
GetErrcon() << std::endl;
1067 os <<
" fMinimumStep = " << magDrv.
GetHmin() << std::endl;
static const G4double eps
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
void PrintInfo(const G4MagInt_Driver &magDrv, std::ostream &os)
static constexpr double perMillion
static constexpr double perThousand
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cerr
G4GLOB_DLL std::ostream G4cout
double dot(const Hep3Vector &) const
static void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, unsigned int subStepNo, unsigned int noIntegrationVariables)
const G4ThreeVector & GetMomentumDir() const
G4double GetCurveLength() const
void SetCurveLength(G4double nCurve_s)
G4double GetKineticEnergy() const
G4ThreeVector GetPosition() const
void DumpToArray(G4double valArr[ncompSVEC]) const
void LoadFromArray(const G4double valArr[ncompSVEC], G4int noVarsIntegrated)
G4double GetSurfaceTolerance() const
static G4GeometryTolerance * GetInstance()
G4double GetPshrnk() const
unsigned long fNoSmallSteps
G4double ComputeNewStepSize_WithinLimits(G4double errMaxNorm, G4double hstepCurrent)
virtual void GetDerivatives(const G4FieldTrack &y_curr, G4double dydx[]) const override
unsigned long fNoTotalSteps
G4MagInt_Driver(G4double hminimum, G4MagIntegratorStepper *pItsStepper, G4int numberOfComponents=6, G4int statisticsVerbosity=0)
void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)
G4MagIntegratorStepper * pIntStepper
virtual ~G4MagInt_Driver() override
void PrintStatisticsReport()
G4int fStatisticsVerboseLevel
virtual void SetEquationOfMotion(G4EquationOfMotion *equation) override
virtual void RenewStepperAndAdjust(G4MagIntegratorStepper *pItsStepper) override
G4double GetErrcon() const
void SetSmallestFraction(G4double val)
virtual const G4MagIntegratorStepper * GetStepper() const override
unsigned long fNoInitialSmallSteps
virtual G4bool QuickAdvance(G4FieldTrack &y_val, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr) override
G4double GetSafety() const
virtual G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0) override
void WarnSmallStepSize(G4double hnext, G4double hstep, G4double h, G4double xDone, G4int noSteps)
virtual G4int GetVerboseLevel() const override
void WarnEndPointTooFar(G4double endPointDist, G4double hStepSize, G4double epsilonRelative, G4int debugFlag)
G4int GetMaxNoSteps() const
void OneGoodStep(G4double ystart[], const G4double dydx[], G4double &x, G4double htry, G4double eps, G4double &hdid, G4double &hnext)
G4double ComputeNewStepSize_WithoutReductionLimit(G4double errMaxNorm, G4double hstepCurrent)
const G4int fNoIntegrationVariables
virtual G4EquationOfMotion * GetEquationOfMotion() override
void PrintStat_Aux(const G4FieldTrack &aFieldTrack, G4double requestStep, G4double actualStep, G4int subStepNo, G4double subStepSize, G4double dotVelocities)
G4double GetSmallestFraction() const
G4double GetPgrow() const
virtual G4bool DoesReIntegrate() const override
virtual G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent) override
G4double fSmallestFraction
void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent)
unsigned long fNoBadSteps
void StreamInfo(std::ostream &os) const override
void ReSetParameters(G4double new_safety=0.9)
G4EquationOfMotion * GetEquationOfMotion()
void RightHandSide(const G4double y[], G4double dydx[]) const
void SetEquationOfMotion(G4EquationOfMotion *newEquation)
virtual G4int IntegratorOrder() const =0
static constexpr G4double max_stepping_decrease
static constexpr G4double max_stepping_increase
T max(const T t1, const T t2)
brief Return the largest of the two arguments