#include <G4FSALIntegrationDriver.hh>

Public Member Functions
G4bool	AccurateAdvance (G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)

G4bool	QuickAdvance (G4FieldTrack &y_val, G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr)

G4bool	QuickAdvance (G4FieldTrack &y_posvel, G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr_pos_sq, G4double &dyerr_mom_rel_sq)

	G4FSALIntegrationDriver (G4double hminimum, G4VFSALIntegrationStepper *pItsStepper, G4int numberOfComponents=6, G4int statisticsVerbosity=1)

	~G4FSALIntegrationDriver ()

G4double	GetHmin () const

G4double	Hmin () const

G4double	GetSafety () const

G4double	GetPshrnk () const

G4double	GetPgrow () const

G4double	GetErrcon () const

G4int	GetNoTotalSteps () const

void	GetDerivatives (const G4FieldTrack &y_curr, G4double dydx[])

void	RenewStepperAndAdjust (G4VFSALIntegrationStepper *pItsStepper)

void	ReSetParameters (G4double new_safety=0.9)

void	SetSafety (G4double valS)

void	SetPshrnk (G4double valPs)

void	SetPgrow (G4double valPg)

void	SetErrcon (G4double valEc)

G4double	ComputeAndSetErrcon ()

const G4VFSALIntegrationStepper *	GetStepper () const

G4VFSALIntegrationStepper *	GetStepper ()

void	OneGoodStep (G4double ystart[], G4double dydx[], G4double &x, G4double htry, G4double eps, G4double &hdid, G4double &hnext)

G4double	ComputeNewStepSize (G4double errMaxNorm, G4double hstepCurrent)

G4double	ComputeNewStepSize_WithinLimits (G4double errMaxNorm, G4double hstepCurrent)

G4int	GetMaxNoSteps () const

void	SetMaxNoSteps (G4int val)

G4int	GetTotalNoStepperCalls () const

void	SetHmin (G4double newval)

void	SetVerboseLevel (G4int newLevel)

G4double	GetVerboseLevel () const

G4double	GetSmallestFraction () const

void	SetSmallestFraction (G4double val)

Protected Member Functions
void	WarnSmallStepSize (G4double hnext, G4double hstep, G4double h, G4double xDone, G4int noSteps)

void	WarnTooManyStep (G4double x1start, G4double x2end, G4double xCurrent)

void	WarnEndPointTooFar (G4double endPointDist, G4double hStepSize, G4double epsilonRelative, G4int debugFlag)

void	PrintStatus (const G4double StartArr, G4double xstart, const G4double CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)

void	PrintStatus (const G4FieldTrack &StartFT, const G4FieldTrack &CurrentFT, G4double requestStep, G4int subStepNo)

void	PrintStat_Aux (const G4FieldTrack &aFieldTrack, G4double requestStep, G4double actualStep, G4int subStepNo, G4double subStepSize, G4double dotVelocities)

void	PrintStatisticsReport ()

Detailed Description

Definition at line 48 of file G4FSALIntegrationDriver.hh.

Constructor & Destructor Documentation

G4FSALIntegrationDriver::G4FSALIntegrationDriver	(	G4double	hminimum,
		G4VFSALIntegrationStepper *	pItsStepper,
		G4int	numberOfComponents = `6`,
		G4int	statisticsVerbosity = `1`
	)

Definition at line 67 of file G4FSALIntegrationDriver.cc.

 : fSmallestFraction( 1.0e-12 ),
 fNoIntegrationVariables(numComponents),
 fMinNoVars(12),
 fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )),
 fStatisticsVerboseLevel(statisticsVerbose),
 fNoTotalSteps(0),  fNoBadSteps(0), fNoSmallSteps(0),
 fNoInitialSmallSteps(0),
 fDyerr_max(0.0), fDyerr_mx2(0.0),
 fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0),
 fSumH_sm(0.0), fSumH_lg(0.0),
 fVerboseLevel(0),
 TotalNoStepperCalls(0)
 {
     // In order to accomodate "Laboratory Time", which is [7], fMinNoVars=8
     // is required. For proper time of flight and spin,  fMinNoVars must be 12
     
     RenewStepperAndAdjust( pStepper );
     fMinimumStep= hminimum;
     fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
 #ifdef G4DEBUG_FIELD
     fVerboseLevel=2;
 #endif
     
     if( (fVerboseLevel > 0) || (fStatisticsVerboseLevel > 1) )
     {
         G4cout << "MagIntDriver version: Accur-Adv: "
         << "invE_nS, QuickAdv-2sqrt with Statistics "
 #ifdef G4FLD_STATS
         << " enabled "
 #else
         << " disabled "
 #endif
         << G4endl;
     }
 }

Here is the call graph for this function:

G4FSALIntegrationDriver::~G4FSALIntegrationDriver ( )

Definition at line 111 of file G4FSALIntegrationDriver.cc.

 {
     if( fStatisticsVerboseLevel > 1 )
     {
         PrintStatisticsReport();
     }
 }

Here is the call graph for this function:

Member Function Documentation

G4bool G4FSALIntegrationDriver::AccurateAdvance	(	G4FieldTrack &	y_current,
		G4double	hstep,
		G4double	eps,
		G4double	hinitial = `0.0`
	)

Definition at line 126 of file G4FSALIntegrationDriver.cc.

 {
     // Runge-Kutta driver with adaptive stepsize control. Integrate starting
     // values at y_current over hstep x2 with accuracy eps.
     // On output ystart is replaced by values at the end of the integration
     // interval. RightHandSide is the right-hand side of ODE system.
     // The source is similar to odeint routine from NRC p.721-722 .
     
     G4int nstp, i, no_warnings=0;
     G4double x, hnext, hdid, h;
     
 #ifdef G4DEBUG_FIELD
     static G4int dbg=10;
     static G4int nStpPr=50;   // For debug printing of long integrations
     G4double ySubStepStart[G4FieldTrack::ncompSVEC];
     G4FieldTrack  yFldTrkStart(y_current);
 #endif
     
     G4double y[G4FieldTrack::ncompSVEC], dydx[G4FieldTrack::ncompSVEC];
     G4double ystart[G4FieldTrack::ncompSVEC], yEnd[G4FieldTrack::ncompSVEC];
     G4double  x1, x2;
     G4bool succeeded = true, lastStepSucceeded;
     
     G4double startCurveLength;
     
     G4int  noFullIntegr=0, noSmallIntegr = 0 ;
     static G4ThreadLocal G4int  noGoodSteps =0 ;  // Bad = chord > curve-len
     const  G4int  nvar= fNoVars;
     
     G4FieldTrack yStartFT(y_current);
     
     //  Ensure that hstep > 0
     //
     if( hstep <= 0.0 )
     {
         if(hstep==0.0)
         {
             std::ostringstream message;
             message << "Proposed step is zero; hstep = " << hstep << " !";
             G4Exception("G4FSALIntegrationDriver::AccurateAdvance()",
                         "GeomField1001", JustWarning, message);
             return succeeded;
         }
         else
         {
             std::ostringstream message;
             message << "Invalid run condition." << G4endl
             << "Proposed step is negative; hstep = " << hstep << "." << G4endl
             << "Requested step cannot be negative! Aborting event.";
             G4Exception("G4FSALIntegrationDriver::AccurateAdvance()",
                         "GeomField0003", EventMustBeAborted, message);
             return false;
         }
     }
     
     y_current.DumpToArray( ystart );
     
     startCurveLength= y_current.GetCurveLength();
     x1= startCurveLength;
     x2= x1 + hstep;
     
     if ( (hinitial > 0.0) && (hinitial < hstep)
         && (hinitial > perMillion * hstep) )
     {
         h = hinitial;
     }
     else  //  Initial Step size "h" defaults to the full interval
     {
         h = hstep;
     }
     
     x = x1;
     
     for (i=0;i<nvar;i++)  { y[i] = ystart[i]; }
     
     G4bool lastStep= false;
     nstp=1;
     
     pIntStepper->ComputeRightHandSide( y, dydx );
 
     
     do
     {
         G4ThreeVector StartPos( y[0], y[1], y[2] );
         
 #ifdef G4DEBUG_FIELD
         G4double xSubStepStart= x;
         for (i=0;i<nvar;i++)  { ySubStepStart[i] = y[i]; }
         yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
         yFldTrkStart.SetCurveLength(x);
 #endif
         
         
         
         
         fNoTotalSteps++;
         
         // Perform the Integration
         //
         if( h > fMinimumStep )
         {
             OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ;
             //--------------------------------------
             lastStepSucceeded= (hdid == h);
 #ifdef G4DEBUG_FIELD
             if (dbg>2) {
                 PrintStatus( ySubStepStart, xSubStepStart, y, x, h,  nstp); // Only
             }
 #endif
         }
         else
         {
             G4FieldTrack yFldTrk( G4ThreeVector(0,0,0),
                                  G4ThreeVector(0,0,0), 0., 0., 0., 0. );
             G4double dchord_step, dyerr, dyerr_len;   // What to do with these ?
             yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
             yFldTrk.SetCurveLength( x );
             
             QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
             //-----------------------------------------------------
             
             yFldTrk.DumpToArray(y);
             
 #ifdef G4FLD_STATS
             fNoSmallSteps++;
             if ( dyerr_len > fDyerr_max)  { fDyerr_max= dyerr_len; }
             fDyerrPos_smTot += dyerr_len;
             fSumH_sm += h;  // Length total for 'small' steps
             if (nstp<=1)  { fNoInitialSmallSteps++; }
 #endif
 #ifdef G4DEBUG_FIELD
             if (dbg>1)
             {
                 if(fNoSmallSteps<2) { PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
                 G4cout << "Another sub-min step, no " << fNoSmallSteps
                 << " of " << fNoTotalSteps << " this time " << nstp << G4endl;
                 PrintStatus( ySubStepStart, x1, y, x, h,  nstp);   // Only this
                 G4cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h
                 << " epsilon= " << eps << " hstep= " << hstep
                 << " h= " << h << " hmin= " << fMinimumStep << G4endl;
             }
 #endif
             if( h == 0.0 )
             {
                 G4Exception("G4FSALIntegrationDriver::AccurateAdvance()",
                             "GeomField0003", FatalException,
                             "Integration Step became Zero!");
             }
             dyerr = dyerr_len / h;
             hdid= h;
             x += hdid;
             
             // Compute suggested new step
             hnext= ComputeNewStepSize( dyerr/eps, h);
             
             // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
             lastStepSucceeded= (dyerr<= eps);
         }
         
         if (lastStepSucceeded)  { noFullIntegr++; }
         else                    { noSmallIntegr++; }
         
         G4ThreeVector EndPos( y[0], y[1], y[2] );
         
 #ifdef  G4DEBUG_FIELD
         if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
         {
             if( nstp==nStpPr )  { G4cout << "***** Many steps ****" << G4endl; }
             G4cout << "MagIntDrv: " ;
             G4cout << "hdid="  << std::setw(12) << hdid  << " "
             << "hnext=" << std::setw(12) << hnext << " "
             << "hstep=" << std::setw(12) << hstep << " (requested) "
             << G4endl;
             PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
         }
 #endif
         
         // Check the endpoint
         G4double endPointDist= (EndPos-StartPos).mag();
         if ( endPointDist >= hdid*(1.+perMillion) )
         {
             fNoBadSteps++;
             
             // Issue a warning only for gross differences -
             // we understand how small difference occur.
             if ( endPointDist >= hdid*(1.+perThousand) )
             {
 #ifdef G4DEBUG_FIELD
                 if (dbg)
                 {
                     WarnEndPointTooFar ( endPointDist, hdid, eps, dbg );
                     G4cerr << "  Total steps:  bad " << fNoBadSteps
                     << " good " << noGoodSteps << " current h= " << hdid
                     << G4endl;
                     PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
                 }
 #endif
                 no_warnings++;
             }
         }
         else
         {
             noGoodSteps ++;
         }
         // #endif
         
         //  Avoid numerous small last steps
         if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
         {
             // No more integration -- the next step will not happen
             lastStep = true;
         }
         else
         {
             // Check the proposed next stepsize
             if(std::fabs(hnext) <= Hmin())
             {
 #ifdef  G4DEBUG_FIELD
                 // If simply a very small interval is being integrated, do not warn
                 if( (x < x2 * (1-eps) ) &&        // The last step can be small: OK
                    (std::fabs(hstep) > Hmin()) ) // and if we are asked, it's OK
                 {
                     if(dbg>0)
                     {
                         WarnSmallStepSize( hnext, hstep, h, x-x1, nstp );
                         PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
                     }
                     no_warnings++;
                 }
 #endif
                 // Make sure that the next step is at least Hmin.
                 h = Hmin();
             }
             else
             {
                 h = hnext;
             }
             
             //  Ensure that the next step does not overshoot
             if ( x+h > x2 )
             {                // When stepsize overshoots, decrease it!
                 h = x2 - x ;   // Must cope with difficult rounding-error
             }                // issues if hstep << x2
             
             if ( h == 0.0 )
             {
                 // Cannot progress - accept this as last step - by default
                 lastStep = true;
 #ifdef G4DEBUG_FIELD
                 if (dbg>2)
                 {
                     int prec= G4cout.precision(12);
                     G4cout << "Warning: FSALMagIntegratorDriver::AccurateAdvance"
                     << G4endl
                     << "  Integration step 'h' became "
                     << h << " due to roundoff. " << G4endl
                     << " Calculated as difference of x2= "<< x2 << " and x=" << x
                     << "  Forcing termination of advance." << G4endl;
                     G4cout.precision(prec);
                 }
 #endif
             }
         }
     } while ( ((nstp++)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
     // Loop checking, 07.10.2016, J. Apostolakis
 
     // Have we reached the end ?
     // --> a better test might be x-x2 > an_epsilon
     
     succeeded=  (x>=x2);  // If it was a "forced" last step
     
     for (i=0;i<nvar;i++)  { yEnd[i] = y[i]; }
     
     // Put back the values.
     y_current.LoadFromArray( yEnd, fNoIntegrationVariables );
     y_current.SetCurveLength( x );
     
     if(nstp > fMaxNoSteps)
     {
         no_warnings++;
         succeeded = false;
 #ifdef G4DEBUG_FIELD
         if (dbg)
         {
             WarnTooManyStep( x1, x2, x );  //  Issue WARNING
             PrintStatus( yEnd, x1, y, x, hstep, -nstp);
         }
 #endif
     }
     
 #ifdef G4DEBUG_FIELD
     if( dbg && no_warnings )
     {
         G4cerr << "FSALMagIntegratorDriver exit status: no-steps " << nstp <<G4endl;
         PrintStatus( yEnd, x1, y, x, hstep, nstp);
     }
 #endif
     
     return succeeded;
 }  // end of AccurateAdvance ...........................

Here is the call graph for this function:

G4double G4FSALIntegrationDriver::ComputeAndSetErrcon ( )

inline

G4double G4FSALIntegrationDriver::ComputeNewStepSize	(	G4double	errMaxNorm,
		G4double	hstepCurrent
	)

Definition at line 770 of file G4FSALIntegrationDriver.cc.

 {
     G4double hnew;
     
     // Compute size of next Step for a failed step
     if(errMaxNorm > 1.0 )
     {
         // Step failed; compute the size of retrial Step.
         hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
     } else if(errMaxNorm > 0.0 ) {
         // Compute size of next Step for a successful step
         hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
     } else {
         // if error estimate is zero (possible) or negative (dubious)
         hnew = max_stepping_increase * hstepCurrent;
     }
     
     return hnew;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4FSALIntegrationDriver::ComputeNewStepSize_WithinLimits	(	G4double	errMaxNorm,
		G4double	hstepCurrent
	)

Definition at line 800 of file G4FSALIntegrationDriver.cc.

 {
     G4double hnew;
     
     // Compute size of next Step for a failed step
     if (errMaxNorm > 1.0 )
     {
         // Step failed; compute the size of retrial Step.
         hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
         
         if (hnew < max_stepping_decrease*hstepCurrent)
         {
             hnew = max_stepping_decrease*hstepCurrent ;
             // reduce stepsize, but no more
             // than this factor (value= 1/10)
         }
     }
     else
     {
         // Compute size of next Step for a successful step
         if (errMaxNorm > errcon)
         { hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()); }
         else  // No more than a factor of 5 increase
         { hnew = max_stepping_increase * hstepCurrent; }
     }
     return hnew;
 }

Here is the call graph for this function:

void G4FSALIntegrationDriver::GetDerivatives	(	const G4FieldTrack &	y_curr,
		G4double	dydx[]
	)

inline

G4double G4FSALIntegrationDriver::GetErrcon ( ) const

inline

G4double G4FSALIntegrationDriver::GetHmin ( ) const

inline

G4int G4FSALIntegrationDriver::GetMaxNoSteps ( ) const

inline

G4int G4FSALIntegrationDriver::GetNoTotalSteps ( ) const

inline

G4double G4FSALIntegrationDriver::GetPgrow ( ) const

inline

Here is the caller graph for this function:

G4double G4FSALIntegrationDriver::GetPshrnk ( ) const

inline

Here is the caller graph for this function:

G4double G4FSALIntegrationDriver::GetSafety ( ) const

inline

Here is the caller graph for this function:

G4double G4FSALIntegrationDriver::GetSmallestFraction ( ) const

inline

const G4VFSALIntegrationStepper* G4FSALIntegrationDriver::GetStepper ( ) const

inline

G4VFSALIntegrationStepper* G4FSALIntegrationDriver::GetStepper ( )

inline

G4int G4FSALIntegrationDriver::GetTotalNoStepperCalls ( ) const

inline

G4double G4FSALIntegrationDriver::GetVerboseLevel ( ) const

inline

G4double G4FSALIntegrationDriver::Hmin ( ) const

inline

Here is the caller graph for this function:

void G4FSALIntegrationDriver::OneGoodStep	(	G4double	ystart[],
		G4double	dydx[],
		G4double &	x,
		G4double	htry,
		G4double	eps,
		G4double &	hdid,
		G4double &	hnext
	)

Definition at line 518 of file G4FSALIntegrationDriver.cc.

 {
     G4double errmax_sq;
     G4double h, htemp, xnew ;
     
     G4double yerr[G4FieldTrack::ncompSVEC], ytemp[G4FieldTrack::ncompSVEC];
     
     h = htry ; // Set stepsize to the initial trial value
     
     G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
     
     G4double errpos_sq=0.0;    // square of displacement error
     G4double errvel_sq=0.0;    // square of momentum vector difference
     G4double errspin_sq=0.0;   // square of spin vector difference
     
     G4int iter;
     
     static G4ThreadLocal G4int tot_no_trials=0;
     //-----------------------------------------------------------
     //Made a change :
     //max_trials = 100
     const G4int max_trials=100000;
     //[hackabot]
     //------------------------------------------------------------
     G4ThreeVector Spin(y[9],y[10],y[11]);
     G4double   spin_mag2 =Spin.mag2() ;
     G4bool     hasSpin= (spin_mag2 > 0.0);
     
     for (iter=0; iter<max_trials ;iter++)
     {
         tot_no_trials++;
         pIntStepper-> Stepper(y,dydx,h,ytemp,yerr, dydx);
         //            *******
         TotalNoStepperCalls++;
         G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep);
         G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos);
         
         // Evaluate accuracy
         //
         errpos_sq =  sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
         errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
         
         // Accuracy for momentum
         G4double magvel_sq=  sqr(y[3]) + sqr(y[4]) + sqr(y[5]) ;
         G4double sumerr_sq =  sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) ;
         if( magvel_sq > 0.0 ) {
             errvel_sq = sumerr_sq / magvel_sq;
         }else{
             G4cerr << "** G4MagIntegrationDriver: found case of zero momentum."
             << " iteration=  " << iter << " h= " << h << G4endl;
             errvel_sq = sumerr_sq;
         }
         errvel_sq *= inv_eps_vel_sq;
         errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
         
         if( hasSpin )
         {
             // Accuracy for spin
             errspin_sq =  ( sqr(yerr[9]) + sqr(yerr[10]) + sqr(yerr[11]) )
             /  spin_mag2; // ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
             errspin_sq *= inv_eps_vel_sq;
             errmax_sq = std::max( errmax_sq, errspin_sq );
         }
         
         if ( errmax_sq <= 1.0 )  { break; } // Step succeeded.
         
         // Step failed; compute the size of retrial Step.
         htemp = GetSafety()*h* std::pow( errmax_sq, 0.5*GetPshrnk() );
         
         if (htemp >= 0.1*h)  { h = htemp; }  // Truncation error too large,
         else  { h = 0.1*h; }                 // reduce stepsize, but no more
         // than a factor of 10
         xnew = x + h;
         if(xnew == x)
         {
             G4cerr << "FSALMagIntegratorDriver::OneGoodStep:" << G4endl
             << "  Stepsize underflow in Stepper " << G4endl ;
             G4cerr << "  Step's start x=" << x << " and end x= " << xnew
             << " are equal !! " << G4endl
             <<"  Due to step-size= " << h
             << " . Note that input step was " << htry << G4endl;
             break;
         }
     }
     
 #ifdef G4FLD_STATS
     // Sum of squares of position error // and momentum dir (underestimated)
     fSumH_lg += h;
     fDyerrPos_lgTot += errpos_sq;
     fDyerrVel_lgTot += errvel_sq * h * h;
 #endif
     
     // Compute size of next Step
     if (errmax_sq > errcon*errcon)
     {
         hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow());
     }
     else
     {
         hnext = max_stepping_increase*h ; // No more than a factor of 5 increase
     }
     x += (hdid = h);
     
     for(G4int k=0;k<fNoIntegrationVariables;k++) { y[k] = ytemp[k]; }
     
     return;
 }   // end of  OneGoodStep .............................

Here is the call graph for this function:

Here is the caller graph for this function:

void G4FSALIntegrationDriver::PrintStat_Aux	(	const G4FieldTrack &	aFieldTrack,
		G4double	requestStep,
		G4double	actualStep,
		G4int	subStepNo,
		G4double	subStepSize,
		G4double	dotVelocities
	)

protected

Definition at line 935 of file G4FSALIntegrationDriver.cc.

 {
     const G4ThreeVector Position=      aFieldTrack.GetPosition();
     const G4ThreeVector UnitVelocity=  aFieldTrack.GetMomentumDir();
     
     if( subStepNo >= 0)
     {
         G4cout << std::setw( 5) << subStepNo << " ";
     }
     else
     {
         G4cout << std::setw( 5) << "Start" << " ";
     }
     G4double curveLen= aFieldTrack.GetCurveLength();
     G4cout << std::setw( 7) << curveLen;
     G4cout << std::setw( 9) << Position.x() << " "
     << std::setw( 9) << Position.y() << " "
     << std::setw( 9) << Position.z() << " "
     << std::setw( 8) << UnitVelocity.x() << " "
     << std::setw( 8) << UnitVelocity.y() << " "
     << std::setw( 8) << UnitVelocity.z() << " ";
     G4int oldprec= G4cout.precision(3);
     G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 << " ";
     G4cout.precision(6);
     G4cout << std::setw(10) << dotVeloc_StartCurr << " ";
     G4cout.precision(oldprec);
     G4cout << std::setw( 7) << aFieldTrack.GetKineticEnergy();
     G4cout << std::setw(12) << step_len << " ";
     
     static G4ThreadLocal G4double oldCurveLength= 0.0;
     static G4ThreadLocal G4double oldSubStepLength= 0.0;
     static G4ThreadLocal G4int oldSubStepNo= -1;
     
     G4double subStep_len=0.0;
     if( curveLen > oldCurveLength )
     {
         subStep_len= curveLen - oldCurveLength;
     }
     else if (subStepNo == oldSubStepNo)
     {
         subStep_len= oldSubStepLength;
     }
     oldCurveLength= curveLen;
     oldSubStepLength= subStep_len;
     
     G4cout << std::setw(12) << subStep_len << " ";
     G4cout << std::setw(12) << subStepSize << " ";
     if( requestStep != -1.0 )
     {
         G4cout << std::setw( 9) << requestStep << " ";
     }
     else
     {
         G4cout << std::setw( 9) << " InitialStep " << " ";
     }
     G4cout << G4endl;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void G4FSALIntegrationDriver::PrintStatisticsReport ( )

protected

Definition at line 1001 of file G4FSALIntegrationDriver.cc.

 {
     G4int noPrecBig= 6;
     G4int oldPrec= G4cout.precision(noPrecBig);
     
     G4cout << "G4FSALIntegrationDriver Statistics of steps undertaken. " << G4endl;
     G4cout << "G4FSALIntegrationDriver: Number of Steps: "
     << " Total= " <<  fNoTotalSteps
     << " Bad= "   <<  fNoBadSteps 
     << " Small= " <<  fNoSmallSteps 
     << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
     << G4endl;
     
 #ifdef G4FLD_STATS
     G4cout << "MID dyerr: " 
     << " maximum= " << fDyerr_max 
     << " Sum small= " << fDyerrPos_smTot 
     << " std::sqrt(Sum large^2): pos= " << std::sqrt(fDyerrPos_lgTot)
     << " vel= " << std::sqrt( fDyerrVel_lgTot )
     << " Total h-distance: small= " << fSumH_sm 
     << " large= " << fSumH_lg
     << G4endl;
     
 #if 0
     G4int noPrecSmall=4; 
     // Single line precis of statistics ... optional
     G4cout.precision(noPrecSmall);
     G4cout << "MIDnums: " << fMinimumStep
     << "   " << fNoTotalSteps 
     << "  "  <<  fNoSmallSteps
     << "  "  << fNoSmallSteps-fNoInitialSmallSteps
     << "  "  << fNoBadSteps         
     << "   " << fDyerr_max
     << "   " << fDyerr_mx2 
     << "   " << fDyerrPos_smTot 
     << "   " << fSumH_sm
     << "   " << fDyerrPos_lgTot
     << "   " << fDyerrVel_lgTot
     << "   " << fSumH_lg
     << G4endl;
 #endif 
 #endif 
     
     G4cout.precision(oldPrec);
 }

Here is the caller graph for this function:

void G4FSALIntegrationDriver::PrintStatus	(	const G4double *	StartArr,
		G4double	xstart,
		const G4double *	CurrentArr,
		G4double	xcurrent,
		G4double	requestStep,
		G4int	subStepNo
	)

protected

Definition at line 832 of file G4FSALIntegrationDriver.cc.

 {
     G4FieldTrack  StartFT(G4ThreeVector(0,0,0),
                           G4ThreeVector(0,0,0), 0., 0., 0., 0. );
     G4FieldTrack  CurrentFT (StartFT);
     
     StartFT.LoadFromArray( StartArr, fNoIntegrationVariables);
     StartFT.SetCurveLength( xstart);
     CurrentFT.LoadFromArray( CurrentArr, fNoIntegrationVariables);
     CurrentFT.SetCurveLength( xcurrent );
     
     PrintStatus(StartFT, CurrentFT, requestStep, subStepNo );
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void G4FSALIntegrationDriver::PrintStatus	(	const G4FieldTrack &	StartFT,
		const G4FieldTrack &	CurrentFT,
		G4double	requestStep,
		G4int	subStepNo
	)

protected

Definition at line 857 of file G4FSALIntegrationDriver.cc.

 {
     G4int verboseLevel= fVerboseLevel;
     static G4ThreadLocal G4int noPrecision= 5;
     G4int oldPrec= G4cout.precision(noPrecision);
     // G4cout.setf(ios_base::fixed,ios_base::floatfield);
     
     const G4ThreeVector StartPosition=       StartFT.GetPosition();
     const G4ThreeVector StartUnitVelocity=   StartFT.GetMomentumDir();
     const G4ThreeVector CurrentPosition=     CurrentFT.GetPosition();
     const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir();
     
     G4double  DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
     
     G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
     G4double subStepSize = step_len;
     
     if( (subStepNo <= 1) || (verboseLevel > 3) )
     {
         subStepNo = - subStepNo;        // To allow printing banner
         
         G4cout << std::setw( 6)  << " " << std::setw( 25)
         << " G4FSALIntegrationDriver: Current Position  and  Direction" << " "
         << G4endl;
         G4cout << std::setw( 5) << "Step#" << " "
         << std::setw( 7) << "s-curve" << " "
         << std::setw( 9) << "X(mm)" << " "
         << std::setw( 9) << "Y(mm)" << " "
         << std::setw( 9) << "Z(mm)" << " "
         << std::setw( 8) << " N_x " << " "
         << std::setw( 8) << " N_y " << " "
         << std::setw( 8) << " N_z " << " "
         << std::setw( 8) << " N^2-1 " << " "
         << std::setw(10) << " N(0).N " << " "
         << std::setw( 7) << "KinEner " << " "
         << std::setw(12) << "Track-l" << " "   // Add the Sub-step ??
         << std::setw(12) << "Step-len" << " "
         << std::setw(12) << "Step-len" << " "
         << std::setw( 9) << "ReqStep" << " "
         << G4endl;
     }
     
     if( (subStepNo <= 0) )
     {
         PrintStat_Aux( StartFT,  requestStep, 0.,
                       0,        0.0,         1.0);
         //*************
     }
     
     if( verboseLevel <= 3 )
     {
         G4cout.precision(noPrecision);
         PrintStat_Aux( CurrentFT, requestStep, step_len,
                       subStepNo, subStepSize, DotStartCurrentVeloc );
         //*************
     }
     
     else // if( verboseLevel > 3 )
     {
         //  Multi-line output
         
         // G4cout << "Current  Position is " << CurrentPosition << G4endl
         //    << " and UnitVelocity is " << CurrentUnitVelocity << G4endl;
         // G4cout << "Step taken was " << step_len
         //    << " out of PhysicalStep= " <<  requestStep << G4endl;
         // G4cout << "Final safety is: " << safety << G4endl;
         // G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
         //        << G4endl << G4endl;
     }
     G4cout.precision(oldPrec);
 }

Here is the call graph for this function:

G4bool G4FSALIntegrationDriver::QuickAdvance	(	G4FieldTrack &	y_val,
		G4double	dydx[],
		G4double	hstep,
		G4double &	dchord_step,
		G4double &	dyerr
	)

Definition at line 669 of file G4FSALIntegrationDriver.cc.

 {
     G4double dyerr_pos_sq, dyerr_mom_rel_sq;
     G4double yerr_vec[G4FieldTrack::ncompSVEC],
     yarrin[G4FieldTrack::ncompSVEC], yarrout[G4FieldTrack::ncompSVEC];
     G4double s_start;
     G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
     
     static G4ThreadLocal G4int no_call=0;
     no_call ++;
     
     // Move data into array
     y_posvel.DumpToArray( yarrin );      //  yarrin  <== y_posvel
     s_start = y_posvel.GetCurveLength();
     
     // Do an Integration Step
     pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec, dydx) ;
     //            *******
     
     // Estimate curve-chord distance
     dchord_step= pIntStepper-> DistChord();
     //                         *********
     
     // Put back the values.  yarrout ==> y_posvel
     y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );
     y_posvel.SetCurveLength( s_start + hstep );
     
 #ifdef  G4DEBUG_FIELD
     if(fVerboseLevel>2)
     {
         G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
         PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep,  1);
     }
 #endif
     
     // A single measure of the error
     //      TO-DO :  account for  energy,  spin, ... ?
     vel_mag_sq   = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
     inv_vel_mag_sq = 1.0 / vel_mag_sq;
     dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
     dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
     dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
     
     // Calculate also the change in the momentum squared also ???
     // G4double veloc_square = y_posvel.GetVelocity().mag2();
     // ...
     
 #ifdef RETURN_A_NEW_STEP_LENGTH
     // The following step cannot be done here because "eps" is not known.
     dyerr_len = std::sqrt( dyerr_len_sq );
     dyerr_len_sq /= eps ;
     
     // Look at the velocity deviation ?
     //  sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
     
     // Set suggested new step
     hstep= ComputeNewStepSize( dyerr_len, hstep);
 #endif
     
     if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) )
     {
         dyerr = std::sqrt(dyerr_pos_sq);
     }
     else
     {
         // Scale it to the current step size - for now
         dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
     }
     
     return true;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

G4bool G4FSALIntegrationDriver::QuickAdvance	(	G4FieldTrack &	y_posvel,
		G4double	dydx[],
		G4double	hstep,
		G4double &	dchord_step,
		G4double &	dyerr_pos_sq,
		G4double &	dyerr_mom_rel_sq
	)

Definition at line 650 of file G4FSALIntegrationDriver.cc.

 {
     G4Exception("G4FSALIntegrationDriver::QuickAdvance()", "GeomField0001",
                 FatalException, "Not yet implemented.");
     
     // Use the parameters of this method, to please compiler
     dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0];
     dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
     return true;
 }

Here is the call graph for this function:

void G4FSALIntegrationDriver::RenewStepperAndAdjust ( G4VFSALIntegrationStepper * pItsStepper )

inline

Here is the caller graph for this function:

void G4FSALIntegrationDriver::ReSetParameters ( G4double new_safety = 0.9 )

inline

void G4FSALIntegrationDriver::SetErrcon ( G4double valEc )

inline

void G4FSALIntegrationDriver::SetHmin ( G4double newval )

inline

void G4FSALIntegrationDriver::SetMaxNoSteps ( G4int val )

inline

void G4FSALIntegrationDriver::SetPgrow ( G4double valPg )

inline

void G4FSALIntegrationDriver::SetPshrnk ( G4double valPs )

inline

void G4FSALIntegrationDriver::SetSafety ( G4double valS )

inline

void G4FSALIntegrationDriver::SetSmallestFraction ( G4double val )

Definition at line 1049 of file G4FSALIntegrationDriver.cc.

 {
     if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
     {
         fSmallestFraction= newFraction;
     }
     else
     { 
         G4cerr << "Warning: SmallestFraction not changed. " << G4endl
         << "  Proposed value was " << newFraction << G4endl
         << "  Value must be between 1.e-8 and 1.e-16" << G4endl;
     }
 }

void G4FSALIntegrationDriver::SetVerboseLevel ( G4int newLevel )

inline

void G4FSALIntegrationDriver::WarnEndPointTooFar	(	G4double	endPointDist,
		G4double	hStepSize,
		G4double	epsilonRelative,
		G4int	debugFlag
	)

protected

Definition at line 482 of file G4FSALIntegrationDriver.cc.

 {
     static G4ThreadLocal G4double maxRelError=0.0;
     G4bool isNewMax, prNewMax;
     
     isNewMax = endPointDist > (1.0 + maxRelError) * h;
     prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
     if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
     
     if( dbg && (h > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance())
        && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+eps) ) ) )
     {
         static G4ThreadLocal G4int noWarnings = 0;
         std::ostringstream message;
         if( (noWarnings ++ < 10) || (dbg>2) )
         {
             message << "The integration produced an end-point which " << G4endl
             << "is further from the start-point than the curve length."
             << G4endl;
         }
         message << "  Distance of endpoints = " << endPointDist
         << ", curve length = " << h << G4endl
         << "  Difference (curveLen-endpDist)= " << (h - endPointDist)
         << ", relative = " << (h-endPointDist) / h
         << ", epsilon =  " << eps;
         G4Exception("G4FSALIntegrationDriver::WarnEndPointTooFar()", "GeomField1001",
                     JustWarning, message);
     }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void G4FSALIntegrationDriver::WarnSmallStepSize	(	G4double	hnext,
		G4double	hstep,
		G4double	h,
		G4double	xDone,
		G4int	noSteps
	)

protected

Definition at line 433 of file G4FSALIntegrationDriver.cc.

 {
     static G4ThreadLocal G4int noWarningsIssued =0;
     const  G4int maxNoWarnings =  10;   // Number of verbose warnings
     std::ostringstream message;
     if( (noWarningsIssued < maxNoWarnings) || fVerboseLevel > 10 )
     {
         message << "The stepsize for the next iteration, " << hnext
         << ", is too small - in Step number " << nstp << "." << G4endl
         << "The minimum for the driver is " << Hmin()  << G4endl
         << "Requested integr. length was " << hstep << " ." << G4endl
         << "The size of this sub-step was " << h     << " ." << G4endl
         << "The integrations has already gone " << xDone;
     }
     else
     {
         message << "Too small 'next' step " << hnext
         << ", step-no: " << nstp << G4endl
         << ", this sub-step: " << h
         << ",  req_tot_len: " << hstep
         << ", done: " << xDone << ", min: " << Hmin();
     }
     G4Exception("G4FSALIntegrationDriver::WarnSmallStepSize()", "GeomField1001",
                 JustWarning, message);
     noWarningsIssued++;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void G4FSALIntegrationDriver::WarnTooManyStep	(	G4double	x1start,
		G4double	x2end,
		G4double	xCurrent
	)

protected

Definition at line 465 of file G4FSALIntegrationDriver.cc.

 {
     std::ostringstream message;
     message << "The number of steps used in the Integration driver"
     << " (Runge-Kutta) is too many." << G4endl
     << "Integration of the interval was not completed !" << G4endl
     << "Only a " << (xCurrent-x1start)*100/(x2end-x1start)
     << " % fraction of it was done.";
     G4Exception("G4FSALIntegrationDriver::WarnTooManyStep()", "GeomField1001",
                 JustWarning, message);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

The documentation for this class was generated from the following files:

source/geant4.10.03.p03/source/geometry/magneticfield/include/G4FSALIntegrationDriver.hh
source/geant4.10.03.p03/source/geometry/magneticfield/src/G4FSALIntegrationDriver.cc

Public Member Functions

Protected Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation