G4SimpleLocator Class Reference

G4SimpleLocator implements the calculation of the intersection point with a boundary when G4PropagationInField is used. It is based on a linear method for finding the intersection point; the difference compared to G4MultiLevelLocator is that no 'depth' algorithm is used in case of slow progress for finding the intersection point. More...

#include <G4SimpleLocator.hh>

Inheritance diagram for G4SimpleLocator:

Public Member Functions
	G4SimpleLocator (G4Navigator *aNavigator)
	~G4SimpleLocator () override
G4bool	EstimateIntersectionPoint (const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin) override
Public Member Functions inherited from G4VIntersectionLocator
	G4VIntersectionLocator (G4Navigator *theNavigator)
virtual	~G4VIntersectionLocator ()
G4bool	IntersectChord (const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=nullptr)
void	SetEpsilonStepFor (G4double EpsilonStep)
void	SetDeltaIntersectionFor (G4double deltaIntersection)
void	SetNavigatorFor (G4Navigator *fNavigator)
void	SetChordFinderFor (G4ChordFinder *fCFinder)
void	SetVerboseFor (G4int fVerbose)
G4int	GetVerboseFor ()
G4double	GetDeltaIntersectionFor ()
G4double	GetEpsilonStepFor ()
G4Navigator *	GetNavigatorFor ()
G4ChordFinder *	GetChordFinderFor ()
void	SetSafetyParametersFor (G4bool UseSafety)
void	AdjustIntersections (G4bool UseCorrection)
G4bool	AreIntersectionsAdjusted ()
void	AddAdjustementOfFoundIntersection (G4bool UseCorrection)
G4bool	GetAdjustementOfFoundIntersection ()
void	printStatus (const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int stepNum)
void	SetCheckMode (G4bool value)
G4bool	GetCheckMode ()

Additional Inherited Members
Static Public Member Functions inherited from G4VIntersectionLocator
static void	printStatus (const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int stepNum, std::ostream &oss, G4int verboseLevel)
Protected Member Functions inherited from G4VIntersectionLocator
G4FieldTrack	ReEstimateEndpoint (const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
G4bool	CheckAndReEstimateEndpoint (const G4FieldTrack &CurrentStartA, const G4FieldTrack &EstimatedEndB, G4FieldTrack &RevisedEndPoint, G4int &errorCode)
G4ThreeVector	GetSurfaceNormal (const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
G4ThreeVector	GetGlobalSurfaceNormal (const G4ThreeVector &CurrentE_Point, G4bool &validNormal)
G4bool	AdjustmentOfFoundIntersection (const G4ThreeVector &A, const G4ThreeVector &CurrentE_Point, const G4ThreeVector &CurrentF_Point, const G4ThreeVector &MomentumDir, const G4bool IntersectAF, G4ThreeVector &IntersectionPoint, G4double &NewSafety, G4double &fPrevSafety, G4ThreeVector &fPrevSftOrigin)
void	ReportTrialStep (G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
G4bool	LocateGlobalPointWithinVolumeAndCheck (const G4ThreeVector &pos)
void	LocateGlobalPointWithinVolumeCheckAndReport (const G4ThreeVector &pos, const G4String &CodeLocationInfo, G4int CheckMode)
void	ReportReversedPoints (std::ostringstream &ossMsg, const G4FieldTrack &StartPointVel, const G4FieldTrack &EndPointVel, G4double NewSafety, G4double epsStep, const G4FieldTrack &CurrentA_PointVelocity, const G4FieldTrack &CurrentB_PointVelocity, const G4FieldTrack &SubStart_PointVelocity, const G4ThreeVector &CurrentE_Point, const G4FieldTrack &ApproxIntersecPointV, G4int sbstp_no, G4int sbstp_no_p, G4int depth)
void	ReportProgress (std::ostream &oss, const G4FieldTrack &StartPointVel, const G4FieldTrack &EndPointVel, G4int substep_no, const G4FieldTrack &A_PtVel, const G4FieldTrack &B_PtVel, G4double safetyLast, G4int depth=-1)
void	ReportImmediateHit (const char *MethodName, const G4ThreeVector &StartPosition, const G4ThreeVector &TrialPoint, G4double tolerance, unsigned long int numCalls)
Protected Attributes inherited from G4VIntersectionLocator
G4double	kCarTolerance
G4int	fVerboseLevel = 0
G4bool	fUseNormalCorrection = false
G4bool	fCheckMode = false
G4bool	fiUseSafety = false
G4Navigator *	fiNavigator
G4ChordFinder *	fiChordFinder = nullptr
G4double	fiEpsilonStep = -1.0
G4double	fiDeltaIntersection = -1.0
G4Navigator *	fHelpingNavigator
G4TouchableHistory *	fpTouchable = nullptr

Detailed Description

G4SimpleLocator implements the calculation of the intersection point with a boundary when G4PropagationInField is used. It is based on a linear method for finding the intersection point; the difference compared to G4MultiLevelLocator is that no 'depth' algorithm is used in case of slow progress for finding the intersection point.

Definition at line 52 of file G4SimpleLocator.hh.

Constructor & Destructor Documentation

G4SimpleLocator()

G4SimpleLocator::G4SimpleLocator

(

G4Navigator *

aNavigator

)

Constructor and default Destructor.

Definition at line 36 of file G4SimpleLocator.cc.

  : G4VIntersectionLocator(theNavigator)
{
}      

~G4SimpleLocator()

G4SimpleLocator::~G4SimpleLocator ( )

overridedefault

Member Function Documentation

EstimateIntersectionPoint()

G4bool G4SimpleLocator::EstimateIntersectionPoint ( const G4FieldTrack & curveStartPointTangent, const G4FieldTrack & curveEndPointTangent, const G4ThreeVector & trialPoint, G4FieldTrack & intersectPointTangent, G4bool & recalculatedEndPoint, G4double & fPreviousSafety, G4ThreeVector & fPreviousSftOrigin )

overridevirtual

If such an intersection exists, this function calculates the intersection point of the true path of the particle with the surface of the current volume (or of one of its daughters).

Note

Should use lateral displacement as measure of convergence.

Parameters

[in]	curveStartPointTangent	Start point tangent track.
[in]	curveEndPointTangent	End point tangent track.
[in]	trialPoint	Trial point.
[out]	intersectPointTangent	Intersection point tangent track.
[out]	recalculatedEndPoint	Flagging if end point was recomputed.
[in,out]	fPreviousSafety	Previous safety distance.
[in,out]	fPreviousSftOrigin	Previous safety point origin.

Returns

Whether intersection exists or not.

Implements G4VIntersectionLocator.

Definition at line 75 of file G4SimpleLocator.cc.

{
  // Find Intersection Point ( A, B, E )  of true path AB - start at E.
 
  G4bool found_approximate_intersection = false;
  G4bool there_is_no_intersection       = false;
  
  G4FieldTrack  CurrentA_PointVelocity = CurveStartPointVelocity; 
  G4FieldTrack  CurrentB_PointVelocity = CurveEndPointVelocity;
  G4ThreeVector CurrentE_Point = TrialPoint;
  G4bool        validNormalAtE = false;
  G4ThreeVector NormalAtEntry;
 
  G4FieldTrack  ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
  G4double      NewSafety = 0.0;
  G4bool last_AF_intersection = false;
  G4bool final_section = true;  // Shows whether current section is last
                                // (i.e. B=full end)
  recalculatedEndPoint = false; 
 
  G4bool restoredFullEndpoint = false;
 
  G4int substep_no = 0;
 
  // Limits for substep number
  //
  const G4int max_substeps  = 100000000;  // Test 120  (old value 100 )
  const G4int warn_substeps = 1000;       //      100  
 
  // Statistics for substeps
  //
  static G4ThreadLocal G4int max_no_seen= -1; 
 
  NormalAtEntry = GetSurfaceNormal( CurrentE_Point, validNormalAtE); 
 
#ifdef G4DEBUG_FIELD
  const G4double tolerance = 1.0e-8; 
  G4ThreeVector  StartPosition= CurveStartPointVelocity.GetPosition(); 
  if( (TrialPoint - StartPosition).mag() < tolerance * CLHEP::mm ) 
  {
     G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 
                 "GeomNav1002", JustWarning,
                 "Intersection point F is exactly at start point A." ); 
  }
#endif
 
  do 
  {
     G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();  
     G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
       
     // F = a point on true AB path close to point E 
     // (the closest if possible)
     //
     ApproxIntersecPointV = GetChordFinderFor()
                            ->ApproxCurvePointV( CurrentA_PointVelocity, 
                                                 CurrentB_PointVelocity, 
                                                 CurrentE_Point,
                                                 GetEpsilonStepFor());
       //  The above method is the key & most intuitive part ...
      
#ifdef G4DEBUG_FIELD
     if( ApproxIntersecPointV.GetCurveLength() > 
         CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
     {
       G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 
                   "GeomNav0003", FatalException,
                   "Intermediate F point is past end B point!" ); 
     }
#endif
 
     G4ThreeVector CurrentF_Point = ApproxIntersecPointV.GetPosition();
 
     // First check whether EF is small - then F is a good approx. point 
     // Calculate the length and direction of the chord AF
     //
     G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
 
     G4ThreeVector NewMomentumDir = ApproxIntersecPointV.GetMomentumDir(); 
     G4double MomDir_dot_Norm = NewMomentumDir.dot( NormalAtEntry ) ;
 
     G4ThreeVector ChordAB = Point_B - Point_A;
 
#ifdef G4DEBUG_FIELD
     G4VIntersectionLocator::
       ReportTrialStep( substep_no, ChordAB, ChordEF_Vector, 
                        NewMomentumDir, NormalAtEntry, validNormalAtE ); 
#endif
     // Check Sign is always exiting !! TODO
     // Could ( > -epsilon) be used instead?
     //
     G4bool adequate_angle = ( MomDir_dot_Norm >= 0.0 ) 
                          || (! validNormalAtE );  // Invalid
     G4double EF_dist2= ChordEF_Vector.mag2();
     if ( ( EF_dist2  <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
       || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
     {
       found_approximate_intersection = true;
        
       // Create the "point" return value
       //
       IntersectedOrRecalculatedFT = ApproxIntersecPointV;
       IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
 
       if ( GetAdjustementOfFoundIntersection() )
       {
         // Try to Get Correction of IntersectionPoint using SurfaceNormal()
         //  
         G4ThreeVector IP;
         G4ThreeVector MomentumDir= ApproxIntersecPointV.GetMomentumDirection();
         G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A,
                                   CurrentE_Point, CurrentF_Point, MomentumDir,
                                   last_AF_intersection, IP, NewSafety,
                                   fPreviousSafety, fPreviousSftOrigin );
 
         if ( goodCorrection )
         {
           IntersectedOrRecalculatedFT = ApproxIntersecPointV;
           IntersectedOrRecalculatedFT.SetPosition(IP);
         }
       }
 
       // Note: in order to return a point on the boundary, 
       //       we must return E. But it is F on the curve.
       //       So we must "cheat": we are using the position at point E
       //       and the velocity at point F !!!
       //
       // This must limit the length we can allow for displacement!
     }
     else  // E is NOT close enough to the curve (ie point F)
     {
       // Check whether any volumes are encountered by the chord AF
       // ---------------------------------------------------------
       // First relocate to restore any Voxel etc information
       // in the Navigator before calling ComputeStep()
       //
       GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A );
 
       G4ThreeVector PointG;   // Candidate intersection point
       G4double stepLengthAF; 
       G4bool usedNavigatorAF = false; 
       G4bool Intersects_AF = IntersectChord( Point_A,   
                                              CurrentF_Point,
                                              NewSafety,
                                              fPreviousSafety,
                                              fPreviousSftOrigin,
                                              stepLengthAF,
                                              PointG,
                                              &usedNavigatorAF );
       last_AF_intersection = Intersects_AF;
       if( Intersects_AF )
       {
         // G is our new Candidate for the intersection point.
         // It replaces  "E" and we will repeat the test to see if
         // it is a good enough approximate point for us.
         //       B    <- F
         //       E    <- G
 
         CurrentB_PointVelocity = ApproxIntersecPointV;
         CurrentE_Point = PointG;
 
         // Need to recalculate the Exit Normal at the new PointG 
         // Relies on a call to Navigator::ComputeStep in IntersectChord above
         // If the safety was adequate (for the step) this would NOT be called!
         // But this must not occur, no intersection can be found in that case,
         // so this branch, ie if( Intersects_AF ) would not be reached!
         //
         G4bool validNormalLast; 
         NormalAtEntry  = GetSurfaceNormal( PointG, validNormalLast ); 
         validNormalAtE = validNormalLast; 
 
         // By moving point B, must take care if current
         // AF has no intersection to try current FB!!
         //
         final_section = false; 
          
#ifdef G4VERBOSE
         if( fVerboseLevel > 3 )
         {
           G4cout << "G4PiF::LI> Investigating intermediate point"
                  << " at s=" << ApproxIntersecPointV.GetCurveLength()
                  << " on way to full s="
                  << CurveEndPointVelocity.GetCurveLength() << G4endl;
         }
#endif
       }
       else  // not Intersects_AF
       {  
         // In this case:
         // There is NO intersection of AF with a volume boundary.
         // We must continue the search in the segment FB!
         //
         GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
 
         G4double stepLengthFB;
         G4ThreeVector PointH;
         G4bool usedNavigatorFB = false; 
 
         // Check whether any volumes are encountered by the chord FB
         // ---------------------------------------------------------
 
         G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B, 
                                                NewSafety,fPreviousSafety,
                                                fPreviousSftOrigin,
                                                stepLengthFB,
                                                PointH, &usedNavigatorFB );
         if( Intersects_FB )
         { 
           // There is an intersection of FB with a volume boundary
           // H <- First Intersection of Chord FB 
 
           // H is our new Candidate for the intersection point.
           // It replaces  "E" and we will repeat the test to see if
           // it is a good enough approximate point for us.
 
           // Note that F must be in volume volA  (the same as A)
           // (otherwise AF would meet a volume boundary!)
           //   A    <- F 
           //   E    <- H
           //
           CurrentA_PointVelocity = ApproxIntersecPointV;
           CurrentE_Point = PointH;
 
           // Need to recalculate the Exit Normal at the new PointG
           // Relies on call to Navigator::ComputeStep in IntersectChord above
           // If safety was adequate (for the step) this would NOT be called!
           // But this must not occur, no intersection found in that case,
           // so this branch, i.e. if( Intersects_AF ) would not be reached!
           //
           G4bool validNormalLast; 
           NormalAtEntry  = GetSurfaceNormal( PointH, validNormalLast ); 
           validNormalAtE = validNormalLast;
         }
         else  // not Intersects_FB
         {
           // There is NO intersection of FB with a volume boundary
 
           if( final_section  )
           {
             // If B is the original endpoint, this means that whatever
             // volume(s) intersected the original chord, none touch the
             // smaller chords we have used.
             // The value of 'IntersectedOrRecalculatedFT' returned is
             // likely not valid 
              
             there_is_no_intersection = true;   // real final_section
           }
           else
           {
             // We must restore the original endpoint
 
             CurrentA_PointVelocity = CurrentB_PointVelocity;  // Got to B
             CurrentB_PointVelocity = CurveEndPointVelocity;
             restoredFullEndpoint = true;
           }
         } // Endif (Intersects_FB)
       } // Endif (Intersects_AF)
 
       // Ensure that the new endpoints are not further apart in space
       // than on the curve due to different errors in the integration
       //
       G4double linDistSq, curveDist; 
       linDistSq = ( CurrentB_PointVelocity.GetPosition() 
                   - CurrentA_PointVelocity.GetPosition() ).mag2(); 
       curveDist = CurrentB_PointVelocity.GetCurveLength()
                   - CurrentA_PointVelocity.GetCurveLength();
 
       // Change this condition for very strict parameters of propagation 
       //
       if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq )
       {
         // Re-integrate to obtain a new B
         //
         G4FieldTrack newEndPointFT =
                 ReEstimateEndpoint( CurrentA_PointVelocity,
                                     CurrentB_PointVelocity,
                                     linDistSq,    // to avoid recalculation
                                     curveDist );
         G4FieldTrack oldPointVelB = CurrentB_PointVelocity; 
         CurrentB_PointVelocity = newEndPointFT;
 
         if( (final_section)) // real final section
         {
           recalculatedEndPoint = true;
           IntersectedOrRecalculatedFT = newEndPointFT;
             // So that we can return it, if it is the endpoint!
         }
       }
       if( curveDist < 0.0 )
       {
         fVerboseLevel = 5; // Print out a maximum of information
         printStatus( CurrentA_PointVelocity,  CurrentB_PointVelocity,
                      -1.0, NewSafety,  substep_no );
         std::ostringstream message;
         message << "Error in advancing propagation." << G4endl
                 << "        Point A (start) is " << CurrentA_PointVelocity
                 << G4endl
                 << "        Point B (end)   is " << CurrentB_PointVelocity
                 << G4endl
                 << "        Curve distance is " << curveDist << G4endl
                 << G4endl
                 << "The final curve point is not further along"
                 << " than the original!" << G4endl;
 
         if( recalculatedEndPoint )
         {
           message << "Recalculation of EndPoint was called with fEpsStep= "
                   << GetEpsilonStepFor() << G4endl;
         }
         message.precision(20);
         message << " Point A (Curve start)   is " << CurveStartPointVelocity
                 << G4endl
                 << " Point B (Curve   end)   is " << CurveEndPointVelocity
                 << G4endl
                 << " Point A (Current start) is " << CurrentA_PointVelocity
                 << G4endl
                 << " Point B (Current end)   is " << CurrentB_PointVelocity
                 << G4endl
                 << " Point E (Trial Point)   is " << CurrentE_Point
                 << G4endl
                 << " Point F (Intersection)  is " << ApproxIntersecPointV
                 << G4endl
                 << "        LocateIntersection parameters are : Substep no= "
                 << substep_no;
 
         G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
                     "GeomNav0003", FatalException, message);
       }
 
       if ( restoredFullEndpoint )
       {
         final_section = restoredFullEndpoint;
         restoredFullEndpoint = false;
       }
     } // EndIf ( E is close enough to the curve, ie point F. )
       // tests ChordAF_Vector.mag() <= maximum_lateral_displacement 
 
#ifdef G4DEBUG_LOCATE_INTERSECTION  
     G4int trigger_substepno_print= warn_substeps - 20;
 
     if( substep_no >= trigger_substepno_print )
     {
       G4cout << "Difficulty in converging in "
              << "G4SimpleLocator::EstimateIntersectionPoint():"
              << G4endl
              << "    Substep no = " << substep_no << G4endl;
       if( substep_no == trigger_substepno_print )
       {
         printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
                      -1.0, NewSafety, 0);
       }
       G4cout << " State of point A: "; 
       printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
                    -1.0, NewSafety, substep_no-1, 0);
       G4cout << " State of point B: "; 
       printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
                    -1.0, NewSafety, substep_no);
     }
#endif
     ++substep_no; 
 
  } while ( ( ! found_approximate_intersection )
           && ( ! there_is_no_intersection )     
           && ( substep_no <= max_substeps) ); // UNTIL found or failed
 
  if( substep_no > max_no_seen )
  {
    max_no_seen = substep_no; 
#ifdef G4DEBUG_LOCATE_INTERSECTION  
    if( max_no_seen > warn_substeps )
    {
      trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps 
    } 
#endif
  }
 
  if(  ( substep_no >= max_substeps)
      && !there_is_no_intersection
      && !found_approximate_intersection )
  {
    G4cout << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()" << G4endl
           << "        Start and Endpoint of Requested Step:" << G4endl;
    printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
                 -1.0, NewSafety, 0);
    G4cout << G4endl
           << "        Start and end-point of current Sub-Step:" << G4endl;
    printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
                 -1.0, NewSafety, substep_no-1);
    printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
                 -1.0, NewSafety, substep_no);
 
    std::ostringstream message;
    message << "Convergence is requiring too many substeps: "
            << substep_no << G4endl
            << "          Abandoning effort to intersect." << G4endl
            << "          Found intersection = "
            << found_approximate_intersection << G4endl
            << "          Intersection exists = "
            << !there_is_no_intersection << G4endl;
    message.precision(10); 
    G4double done_len = CurrentA_PointVelocity.GetCurveLength(); 
    G4double full_len = CurveEndPointVelocity.GetCurveLength();
    message << "          Undertaken only length: " << done_len
            << " out of " << full_len << " required." << G4endl
            << "          Remaining length = " << full_len-done_len; 
 
    G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
                "GeomNav0003", FatalException, message);
  }
  else if( substep_no >= warn_substeps )
  {  
    std::ostringstream message;
    message.precision(10); 
 
    message << "Many substeps while trying to locate intersection." << G4endl
            << "          Undertaken length: "  
            << CurrentB_PointVelocity.GetCurveLength() 
            << " - Needed: "  << substep_no << " substeps." << G4endl
            << "          Warning level = " << warn_substeps
            << " and maximum substeps = " << max_substeps;
    G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
                "GeomNav1002", JustWarning, message);
  }
  return  !there_is_no_intersection; //  Success or failure
}

---

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

• G4SimpleLocator.hh
• G4SimpleLocator.cc