45 : fDefaultDeltaChord( 0.25 *
mm ),
46 fDeltaChord( fDefaultDeltaChord ),
47 fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
48 fMultipleRadius(15.0),
51 fAllocatedStepper(false),
53 fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0)
56 fIntgrDriver= pIntegrationDriver;
57 fAllocatedStepper=
false;
59 fLastStepEstimate_Unconstrained =
DBL_MAX;
71 : fDefaultDeltaChord( 0.25 *
mm ),
72 fDeltaChord( fDefaultDeltaChord ),
73 fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
74 fMultipleRadius(15.0),
77 fAllocatedStepper(false),
79 fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0)
85 fEquation = pEquation;
86 fLastStepEstimate_Unconstrained =
DBL_MAX;
95 if( pItsStepper == 0 )
98 fAllocatedStepper=
true;
102 fAllocatedStepper=
false;
114 if( fAllocatedStepper)
116 delete fDriversStepper;
130 if( fractLast == -1.0 ) fractLast = 1.0;
131 if( fractNext == -1.0 ) fractNext = 0.98;
138 G4cout <<
" ChordFnd> Trying to set fractions: "
139 <<
" first " << fFirstFraction
140 <<
" last " << fractLast
141 <<
" next " << fractNext
142 <<
" and multiple " << fMultipleRadius
146 if( (fractLast > 0.0) && (fractLast <=1.0) )
148 fFractionLast= fractLast;
152 G4cerr <<
"G4ChordFinder::SetFractions_Last_Next: Invalid "
153 <<
" fraction Last = " << fractLast
154 <<
" must be 0 < fractionLast <= 1 " <<
G4endl;
156 if( (fractNext > 0.0) && (fractNext <1.0) )
158 fFractionNextEstimate = fractNext;
162 G4cerr <<
"G4ChordFinder:: SetFractions_Last_Next: Invalid "
163 <<
" fraction Next = " << fractNext
164 <<
" must be 0 < fractionNext < 1 " <<
G4endl;
184 stepPossible=
FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep,
185 &nextStep, latestSafetyOrigin, latestSafetyRadius
191 if ( dyErr < epsStep * stepPossible )
204 if ( ! good_advance )
230 G4double stepTrial, stepForAccuracy;
237 G4bool validEndPoint=
false;
238 G4double dChordStep, lastStepLength;
240 fIntgrDriver-> GetDerivatives( yCurrent, dydx );
243 const G4double safetyFactor= fFirstFraction;
245 stepTrial =
std::min( stepMax, safetyFactor*fLastStepEstimate_Unconstrained );
255 dChordStep, dyErrPos);
261 lastStepLength = stepTrial;
264 stepForChord =
NewStep(stepTrial, dChordStep, newStepEst_Uncons );
266 if( ! validEndPoint )
270 stepTrial = stepForChord;
272 else if (stepForChord <= stepTrial)
275 stepTrial =
std::min( stepForChord, fFractionLast * stepTrial);
284 while( ! validEndPoint );
286 if( newStepEst_Uncons > 0.0 )
288 fLastStepEstimate_Unconstrained= newStepEst_Uncons;
293 if( pStepForAccuracy )
297 G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
298 if( dyErr_relative > 1.0 )
305 stepForAccuracy = 0.0;
307 *pStepForAccuracy = stepForAccuracy;
310 #ifdef TEST_CHORD_PRINT
314 G4cout <<
"ChordF/FindNextChord: NoTrials= " << noTrials
315 <<
" StepForGoodChord=" << std::setw(10) << stepTrial <<
G4endl;
327 G4double& stepEstimate_Unconstrained )
337 if (dChordStep > 0.0)
339 stepEstimate_Unconstrained =
340 stepTrialOld*std::sqrt( fDeltaChord / dChordStep );
341 stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
346 stepTrial = stepTrialOld * 2.;
349 if( stepTrial <= 0.001 * stepTrialOld)
351 if ( dChordStep > 1000.0 * fDeltaChord )
353 stepTrial= stepTrialOld * 0.03;
357 if ( dChordStep > 100. * fDeltaChord )
359 stepTrial= stepTrialOld * 0.1;
363 stepTrial= stepTrialOld * 0.5;
367 else if (stepTrial > 1000.0 * stepTrialOld)
369 stepTrial= 1000.0 * stepTrialOld;
372 if( stepTrial == 0.0 )
379 if ( dChordStep > 1000. * fDeltaChord )
381 stepTrial= stepTrialOld * 0.03;
385 if ( dChordStep > 100. * fDeltaChord )
387 stepTrial= stepTrialOld * 0.1;
391 stepTrial= stepTrialOld * 0.5;
431 if(!first){EndPoint= ApproxCurveV;}
444 ya=(PointG-Point_A).mag();
445 xb=(Point_A-CurrentF_Point).mag();
446 yb=-(PointG-CurrentF_Point).mag();
447 xc=(Point_A-Point_B).mag();
448 yc=-(CurrentE_Point-Point_B).mag();
453 ya=(Point_A-CurrentE_Point).mag();
454 xb=(Point_A-CurrentF_Point).mag();
455 yb=(PointG-CurrentF_Point).mag();
456 xc=(Point_A-Point_B).mag();
457 yc=-(Point_B-PointG).mag();
462 CurrentE_Point, eps_step);
468 if(std::abs(ya)<=tolerance||std::abs(yc)<=tolerance)
484 test_step=(test_step-xb);
487 xb=(CurrentF_Point-Point_B).mag();
490 if(test_step<=0) { test_step=0.1*xb; }
491 if(test_step>=xb) { test_step=0.5*xb; }
492 if(test_step>=curve){ test_step=0.5*curve; }
494 if(curve*(1.+eps_step)<xb)
504 G4cout <<
"G4ChordFinder::ApproxCurvePointS() - test-step ShF = "
505 << test_step <<
" EndPoint = " << EndPoint <<
G4endl;
511 CurveB_PointVelocity,
512 CurrentE_Point, eps_step );
514 G4cout <<
"G4ChordFinder::BrentApprox = " << EndPoint <<
G4endl;
515 G4cout <<
"G4ChordFinder::LinearApprox= " << TestTrack <<
G4endl;
533 G4FieldTrack Current_PointVelocity = CurveA_PointVelocity;
549 if( curve_length < ABdist * (1. - integrationInaccuracyLimit) )
552 G4cerr <<
" Warning in G4ChordFinder::ApproxCurvePoint: "
554 <<
" The two points are further apart than the curve length "
556 <<
" Dist = " << ABdist
557 <<
" curve length = " << curve_length
558 <<
" relativeDiff = " << (curve_length-ABdist)/ABdist
560 if( curve_length < ABdist * (1. - 10*eps_step) )
562 std::ostringstream message;
563 message <<
"Unphysical curve length." <<
G4endl
564 <<
"The size of the above difference exceeds allowed limits."
567 G4Exception(
"G4ChordFinder::ApproxCurvePointV()",
"GeomField0003",
580 AE_fraction = ChordAE_Vector.
mag() / ABdist;
586 G4cout <<
"Warning in G4ChordFinder::ApproxCurvePointV():"
587 <<
" A and B are the same point!" <<
G4endl
588 <<
" Chord AB length = " << ChordAE_Vector.
mag() <<
G4endl
593 if( (AE_fraction> 1.0 +
perMillion) || (AE_fraction< 0.) )
596 G4cerr <<
" G4ChordFinder::ApproxCurvePointV() - Warning:"
597 <<
" Anomalous condition:AE > AB or AE/AB <= 0 " <<
G4endl
598 <<
" AE_fraction = " << AE_fraction <<
G4endl
599 <<
" Chord AE length = " << ChordAE_Vector.
mag() <<
G4endl
601 G4cerr <<
" OK if this condition occurs after a recalculation of 'B'"
602 << G4endl <<
" Otherwise it is an error. " <<
G4endl ;
612 new_st_length= AE_fraction * curve_length;
614 if ( AE_fraction > 0.0 )
617 new_st_length, eps_step );
627 return Current_PointVelocity;
638 G4cout <<
"G4ChordFinder statistics report: " <<
G4endl;
640 <<
" No trials: " << fTotalNoTrials_FNC
641 <<
" No Calls: " << fNoCalls_FNC
642 <<
" Max-trial: " << fmaxTrials_FNC
646 <<
" fFirstFraction " << fFirstFraction
647 <<
" fFractionLast " << fFractionLast
648 <<
" fFractionNextEstimate " << fFractionNextEstimate
661 G4cout <<
" ChF/fnc: notrial " << std::setw( 3) << noTrials
662 <<
" this_step= " << std::setw(10) << lastStepTrial;
663 if( std::fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 )
671 G4cout <<
" dChordStep= " << std::setw(12) << dChordStep;
672 if( dChordStep > fDeltaChord ) {
G4cout <<
" d+"; }
675 G4cout <<
" new_step= " << std::setw(10)
676 << fLastStepEstimate_Unconstrained
677 <<
" new_step_constr= " << std::setw(10)
678 << lastStepTrial <<
G4endl;
679 G4cout <<
" nextStepTrial = " << std::setw(10) << nextStepTrial <<
G4endl;
680 G4cout.precision(oldprec);
G4double InvParabolic(const G4double xa, const G4double ya, const G4double xb, const G4double yb, const G4double xc, const G4double yc)
G4double GetCurveLength() const
G4ChordFinder(G4MagInt_Driver *pIntegrationDriver)
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector ¤tEPoint, G4double epsStep)
virtual void PrintStatistics()
G4int GetNumberOfVariables() const
virtual G4double FindNextChord(const G4FieldTrack &yStart, G4double stepMax, G4FieldTrack &yEnd, G4double &dyErr, G4double epsStep, G4double *pNextStepForAccuracy, const G4ThreeVector latestSafetyOrigin, G4double latestSafetyRadius)
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
void TestChordPrint(G4int noTrials, G4int lastStepTrial, G4double dChordStep, G4double nextStepTrial)
G4bool AcceptableMissDist(G4double dChordStep) const
void AccumulateStatistics(G4int noTrials)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double AdvanceChordLimited(G4FieldTrack &yCurrent, G4double stepInitial, G4double epsStep_Relative, const G4ThreeVector latestSafetyOrigin, G4double lasestSafetyRadius)
G4bool QuickAdvance(G4FieldTrack &y_val, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void SetFractions_Last_Next(G4double fractLast=0.90, G4double fractNext=0.95)
G4double NewStep(G4double stepTrialOld, G4double dChordStep, G4double &stepEstimate_Unconstrained)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent)
G4FieldTrack ApproxCurvePointS(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4FieldTrack &ApproxCurveV, const G4ThreeVector ¤tEPoint, const G4ThreeVector ¤tFPoint, const G4ThreeVector &PointG, G4bool first, G4double epsStep)
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
G4GLOB_DLL std::ostream G4cerr