57 G4int statisticsVerbose)
58 : fSmallestFraction( 1.0
e-12 ),
59 fNoIntegrationVariables(numComponents),
61 fNoVars( std::
max( fNoIntegrationVariables, fMinNoVars )),
62 fStatisticsVerboseLevel(statisticsVerbose),
63 fNoTotalSteps(0), fNoBadSteps(0), fNoSmallSteps(0),
64 fNoInitialSmallSteps(0), fNoCalls(0),
65 fDyerr_max(0.0), fDyerr_mx2(0.0),
66 fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0),
67 fSumH_sm(0.0), fSumH_lg(0.0),
88 G4cout <<
"MagIntDriver version: Accur-Adv: "
89 <<
"invE_nS, QuickAdv-2sqrt with Statistics "
125 G4int nstp, i, no_warnings=0;
130 static G4int nStpPr=50;
138 G4bool succeeded =
true, lastStepSucceeded;
142 G4int noFullIntegr=0, noSmallIntegr = 0 ;
155 message <<
"Proposed step is zero; hstep = " << hstep <<
" !";
163 message <<
"Invalid run condition." <<
G4endl
164 <<
"Proposed step is negative; hstep = " << hstep <<
"." <<
G4endl
165 <<
"Requested step cannot be negative! Aborting event.";
175 x1= startCurveLength;
178 if ( (hinitial > 0.0) && (hinitial < hstep)
190 for (i=0;i<nvar;i++) { y[i] = ystart[i]; }
201 for (i=0;i<nvar;i++) { ySubStepStart[i] = y[i]; }
218 lastStepSucceeded= (hdid == h);
222 PrintStatus( ySubStepStart, xSubStepStart, y, x, h, nstp);
230 G4double dchord_step, dyerr, dyerr_len;
234 QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
253 G4cout <<
" dyerr= " << dyerr_len <<
" relative = " << dyerr_len / h
254 <<
" epsilon= " << eps <<
" hstep= " << hstep
262 "Integration Step became Zero!");
264 dyerr = dyerr_len / h;
272 lastStepSucceeded= (dyerr<=
eps);
275 if (lastStepSucceeded) { noFullIntegr++; }
276 else { noSmallIntegr++; }
281 if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
283 if( nstp==nStpPr ) {
G4cout <<
"***** Many steps ****" <<
G4endl; }
285 G4cout <<
"hdid=" << std::setw(12) << hdid <<
" "
286 <<
"hnext=" << std::setw(12) << hnext <<
" "
287 <<
"hstep=" << std::setw(12) << hstep <<
" (requested) "
289 PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
294 G4double endPointDist= (EndPos-StartPos).mag();
308 <<
" good " << noGoodSteps <<
" current h= " << hdid
310 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
331 if(std::fabs(hnext) <=
Hmin())
335 if( (x < x2 * (1-eps) ) &&
336 (std::fabs(hstep) >
Hmin()) )
341 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
368 G4cout <<
"Warning: G4MagIntegratorDriver::AccurateAdvance"
370 <<
" Integration step 'h' became "
371 << h <<
" due to roundoff. " <<
G4endl
372 <<
" Calculated as difference of x2= "<< x2 <<
" and x=" << x
373 <<
" Forcing termination of advance." <<
G4endl;
379 }
while ( ((nstp++)<=
fMaxNoSteps) && (x < x2) && (!lastStep) );
387 for (i=0;i<nvar;i++) { yEnd[i] = y[i]; }
407 if( dbg && no_warnings )
409 G4cerr <<
"G4MagIntegratorDriver exit status: no-steps " << nstp <<
G4endl;
425 const G4int maxNoWarnings = 10;
427 if( (noWarningsIssued < maxNoWarnings) ||
fVerboseLevel > 10 )
429 message <<
"The stepsize for the next iteration, " << hnext
430 <<
", is too small - in Step number " << nstp <<
"." <<
G4endl
431 <<
"The minimum for the driver is " <<
Hmin() <<
G4endl
432 <<
"Requested integr. length was " << hstep <<
" ." <<
G4endl
433 <<
"The size of this sub-step was " << h <<
" ." <<
G4endl
434 <<
"The integrations has already gone " << xDone;
438 message <<
"Too small 'next' step " << hnext
439 <<
", step-no: " << nstp <<
G4endl
440 <<
", this sub-step: " << h
441 <<
", req_tot_len: " << hstep
442 <<
", done: " << xDone <<
", min: " <<
Hmin();
444 G4Exception(
"G4MagInt_Driver::WarnSmallStepSize()",
"GeomField1001",
457 message <<
"The number of steps used in the Integration driver"
458 <<
" (Runge-Kutta) is too many." <<
G4endl
459 <<
"Integration of the interval was not completed !" <<
G4endl
460 <<
"Only a " << (xCurrent-x1start)*100/(x2end-x1start)
461 <<
" % fraction of it was done.";
462 G4Exception(
"G4MagInt_Driver::WarnTooManyStep()",
"GeomField1001",
475 G4bool isNewMax, prNewMax;
477 isNewMax = endPointDist > (1.0 + maxRelError) * h;
478 prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
479 if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
482 && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+
eps) ) ) )
486 if( (noWarnings ++ < 10) || (dbg>2) )
488 message <<
"The integration produced an end-point which " <<
G4endl
489 <<
"is further from the start-point than the curve length."
492 message <<
" Distance of endpoints = " << endPointDist
493 <<
", curve length = " << h <<
G4endl
494 <<
" Difference (curveLen-endpDist)= " << (h - endPointDist)
495 <<
", relative = " << (h-endPointDist) / h
496 <<
", epsilon = " <<
eps;
497 G4Exception(
"G4MagInt_Driver::WarnEndPointTooFar()",
"GeomField1001",
534 G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
543 const G4int max_trials=100;
547 G4bool hasSpin= (spin_mag2 > 0.0);
549 for (iter=0; iter<max_trials ;iter++)
555 G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos);
559 errpos_sq =
sqr(yerr[0]) +
sqr(yerr[1]) +
sqr(yerr[2]) ;
560 errpos_sq *= inv_eps_pos_sq;
565 if( magvel_sq > 0.0 )
567 errvel_sq = sumerr_sq / magvel_sq;
572 message <<
"Found case of zero momentum." <<
G4endl
573 <<
"- iteration= " << iter <<
"; h= " << h;
576 errvel_sq = sumerr_sq;
578 errvel_sq *= inv_eps_vel_sq;
579 errmax_sq =
std::max( errpos_sq, errvel_sq );
584 errspin_sq = (
sqr(yerr[9]) +
sqr(yerr[10]) +
sqr(yerr[11]) )
586 errspin_sq *= inv_eps_vel_sq;
587 errmax_sq =
std::max( errmax_sq, errspin_sq );
590 if ( errmax_sq <= 1.0 ) {
break; }
595 if (htemp >= 0.1*h) { h = htemp; }
602 message <<
"Stepsize underflow in Stepper !" <<
G4endl
603 <<
"- Step's start x=" << x <<
" and end x= " << xnew
604 <<
" are equal !! " <<
G4endl
605 <<
" Due to step-size= " << h
606 <<
". Note that input step was " << htry;
641 G4Exception(
"G4MagInt_Driver::QuickAdvance()",
"GeomField0001",
645 dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0];
659 G4double dyerr_pos_sq, dyerr_mom_rel_sq;
663 G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
673 pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ;
686 PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep, 1);
692 vel_mag_sq = (
sqr(yarrout[3])+
sqr(yarrout[4])+
sqr(yarrout[5]) );
693 inv_vel_mag_sq = 1.0 / vel_mag_sq;
694 dyerr_pos_sq = (
sqr(yerr_vec[0])+
sqr(yerr_vec[1])+
sqr(yerr_vec[2]));
695 dyerr_mom_sq = (
sqr(yerr_vec[3])+
sqr(yerr_vec[4])+
sqr(yerr_vec[5]));
696 dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
702 #ifdef RETURN_A_NEW_STEP_LENGTH
704 dyerr_len = std::sqrt( dyerr_len_sq );
705 dyerr_len_sq /=
eps ;
714 if( dyerr_pos_sq > ( dyerr_mom_rel_sq *
sqr(hstep) ) )
716 dyerr = std::sqrt(dyerr_pos_sq);
721 dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
729 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT
738 G4Exception(
"G4MagInt_Driver::QuickAdvance()",
"GeomField0001",
740 dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
741 yarrout[0]= yarrin[0];
757 if(errMaxNorm > 1.0 )
762 else if(errMaxNorm > 0.0 )
791 if (errMaxNorm > 1.0 )
836 PrintStatus(StartFT, CurrentFT, requestStep, subStepNo );
848 const G4int noPrecision = 5;
857 G4double DotStartCurrentVeloc= StartUnitVelocity.
dot(CurrentUnitVelocity);
862 if( (subStepNo <= 1) || (verboseLevel > 3) )
864 subStepNo = - subStepNo;
866 G4cout << std::setw( 6) <<
" " << std::setw( 25)
867 <<
" G4MagInt_Driver: Current Position and Direction" <<
" "
869 G4cout << std::setw( 5) <<
"Step#" <<
" "
870 << std::setw( 7) <<
"s-curve" <<
" "
871 << std::setw( 9) <<
"X(mm)" <<
" "
872 << std::setw( 9) <<
"Y(mm)" <<
" "
873 << std::setw( 9) <<
"Z(mm)" <<
" "
874 << std::setw( 8) <<
" N_x " <<
" "
875 << std::setw( 8) <<
" N_y " <<
" "
876 << std::setw( 8) <<
" N_z " <<
" "
877 << std::setw( 8) <<
" N^2-1 " <<
" "
878 << std::setw(10) <<
" N(0).N " <<
" "
879 << std::setw( 7) <<
"KinEner " <<
" "
880 << std::setw(12) <<
"Track-l" <<
" "
881 << std::setw(12) <<
"Step-len" <<
" "
882 << std::setw(12) <<
"Step-len" <<
" "
883 << std::setw( 9) <<
"ReqStep" <<
" "
887 if( (subStepNo <= 0) )
893 if( verboseLevel <= 3 )
895 G4cout.precision(noPrecision);
897 subStepNo, subStepSize, DotStartCurrentVeloc );
900 G4cout.precision(oldPrec);
918 G4cout << std::setw( 5) << subStepNo <<
" ";
922 G4cout << std::setw( 5) <<
"Start" <<
" ";
925 G4cout << std::setw( 7) << curveLen;
926 G4cout << std::setw( 9) << Position.
x() <<
" "
927 << std::setw( 9) << Position.
y() <<
" "
928 << std::setw( 9) << Position.
z() <<
" "
929 << std::setw( 8) << UnitVelocity.
x() <<
" "
930 << std::setw( 8) << UnitVelocity.
y() <<
" "
931 << std::setw( 8) << UnitVelocity.
z() <<
" ";
933 G4cout << std::setw( 8) << UnitVelocity.
mag2()-1.0 <<
" ";
935 G4cout << std::setw(10) << dotVeloc_StartCurr <<
" ";
936 G4cout.precision(oldprec);
938 G4cout << std::setw(12) << step_len <<
" ";
945 if( curveLen > oldCurveLength )
947 subStep_len= curveLen - oldCurveLength;
949 else if (subStepNo == oldSubStepNo)
951 subStep_len= oldSubStepLength;
953 oldCurveLength= curveLen;
954 oldSubStepLength= subStep_len;
956 G4cout << std::setw(12) << subStep_len <<
" ";
957 G4cout << std::setw(12) << subStepSize <<
" ";
958 if( requestStep != -1.0 )
960 G4cout << std::setw( 9) << requestStep <<
" ";
964 G4cout << std::setw( 9) <<
" InitialStep " <<
" ";
976 G4cout <<
"G4MagInt_Driver Statistics of steps undertaken. " <<
G4endl;
977 G4cout <<
"G4MagInt_Driver: Number of Steps: "
983 G4cout.precision(oldPrec);
990 if( (newFraction > 1.
e-16) && (newFraction < 1
e-8) )
997 message <<
"Smallest Fraction not changed. " <<
G4endl
998 <<
" Proposed value was " << newFraction <<
G4endl
999 <<
" Value must be between 1.e-8 and 1.e-16";
1000 G4Exception(
"G4MagInt_Driver::SetSmallestFraction()",
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0) override
void LoadFromArray(const G4double valArr[ncompSVEC], G4int noVarsIntegrated)
CLHEP::Hep3Vector G4ThreeVector
virtual G4bool QuickAdvance(G4FieldTrack &y_val, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr) override
G4EquationOfMotion * GetEquationOfMotion()
void SetEquationOfMotion(G4EquationOfMotion *newEquation)
Float_t x1[n_points_granero]
void message(RunManager *runmanager)
void PrintStatisticsReport()
void DumpToArray(G4double valArr[ncompSVEC]) const
G4double GetPshrnk() const
static constexpr double perMillion
double dot(const Hep3Vector &) const
G4double G4MagInt_Driver::G4double hstepCurrent G4double hnew
void SetSmallestFraction(G4double val)
virtual G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent) override
void WarnSmallStepSize(G4double hnext, G4double hstep, G4double h, G4double xDone, G4int noSteps)
const G4int fNoIntegrationVariables
G4MagInt_Driver(G4double hminimum, G4MagIntegratorStepper *pItsStepper, G4int numberOfComponents=6, G4int statisticsVerbosity=1)
unsigned long fNoSmallSteps
void WarnEndPointTooFar(G4double endPointDist, G4double hStepSize, G4double epsilonRelative, G4int debugFlag)
G4ThreeVector GetPosition() const
G4double GetSafety() const
virtual void GetDerivatives(const G4FieldTrack &y_curr, G4double dydx[]) const override
G4double GetKineticEnergy() const
G4MagIntegratorStepper * pIntStepper
void SetCurveLength(G4double nCurve_s)
G4double max_stepping_decrease
void RightHandSide(const double y[], double dydx[]) const
void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)
G4double max_stepping_increase
virtual void SetEquationOfMotion(G4EquationOfMotion *equation) override
G4double fSmallestFraction
virtual void ComputeRightHandSide(const G4double y[], G4double dydx[])
G4GLOB_DLL std::ostream G4cerr
virtual G4EquationOfMotion * GetEquationOfMotion() override
G4double ComputeNewStepSize_WithinLimits(G4double errMaxNorm, G4double hstepCurrent)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4int fStatisticsVerboseLevel
G4double GetPgrow() const
virtual ~G4MagInt_Driver() override
const G4ThreeVector & GetMomentumDir() const
void ReSetParameters(G4double new_safety=0.9)
G4GLOB_DLL std::ostream G4cout
static G4GeometryTolerance * GetInstance()
virtual G4int IntegratorOrder() const =0
static const G4double eps
unsigned long fNoInitialSmallSteps
unsigned long fNoTotalSteps
unsigned long fNoBadSteps
Float_t x2[n_points_geant4]
void OneGoodStep(G4double ystart[], const G4double dydx[], G4double &x, G4double htry, G4double eps, G4double &hdid, G4double &hnext)
virtual void RenewStepperAndAdjust(G4MagIntegratorStepper *pItsStepper) override
void PrintStat_Aux(const G4FieldTrack &aFieldTrack, G4double requestStep, G4double actualStep, G4int subStepNo, G4double subStepSize, G4double dotVelocities)
void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent)
virtual const G4MagIntegratorStepper * GetStepper() const override
static constexpr double perThousand
G4double GetSurfaceTolerance() const
G4double GetCurveLength() const