G4MultiLevelLocator Class Reference

#include <G4MultiLevelLocator.hh>

Inheritance diagram for G4MultiLevelLocator:

Public Member Functions
	G4MultiLevelLocator (G4Navigator *theNavigator)
	~G4MultiLevelLocator ()
G4bool	EstimateIntersectionPoint (const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin)
void	ReportStatistics ()
void	SetMaxSteps (unsigned int valMax)
void	SetWarnSteps (unsigned int valWarn)
Public Member Functions inherited from G4VIntersectionLocator
	G4VIntersectionLocator (G4Navigator *theNavigator)
virtual	~G4VIntersectionLocator ()
virtual G4bool	EstimateIntersectionPoint (const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin)=0
void	printStatus (const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int stepNum)
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	AddAdjustementOfFoundIntersection (G4bool UseCorrection)
G4bool	GetAdjustementOfFoundIntersection ()
void	AdjustIntersections (G4bool UseCorrection)
G4bool	AreIntersectionsAdjusted ()
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

Definition at line 47 of file G4MultiLevelLocator.hh.

Constructor & Destructor Documentation

G4MultiLevelLocator()

G4MultiLevelLocator::G4MultiLevelLocator

(

G4Navigator *

theNavigator

)

Definition at line 39 of file G4MultiLevelLocator.cc.

  : G4VIntersectionLocator(theNavigator)
{
  // In case of too slow progress in finding Intersection Point
  // intermediates Points on the Track must be stored.
  // Initialise the array of Pointers [max_depth+1] to do this  
  
  G4ThreeVector zeroV(0.0,0.0,0.0);
  for ( auto idepth=0; idepth<max_depth+1; ++idepth )
  {
    ptrInterMedFT[ idepth ] = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.);
  }
 
  if (fCheckMode)
  {
    //  Trial values       Loose         Medium      Tight
    //  To happen:         Infrequent    Occasional  Often
    SetMaxSteps(150);   //  300             85          25
    SetWarnSteps(80);   //  250             60          15
  }  
}      

~G4MultiLevelLocator()

G4MultiLevelLocator::~G4MultiLevelLocator

(

)

Definition at line 61 of file G4MultiLevelLocator.cc.

{
  for ( auto idepth=0; idepth<max_depth+1; ++idepth )
  {
    delete ptrInterMedFT[idepth];
  }
#ifdef G4DEBUG_FIELD
  ReportStatistics(); 
#endif  
}

Member Function Documentation

EstimateIntersectionPoint()

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

virtual

Implements G4VIntersectionLocator.

Definition at line 115 of file G4MultiLevelLocator.cc.

{
  // Find Intersection Point ( A, B, E )  of true path AB - start at E.
  const char* MethodName= "G4MultiLevelLocator::EstimateIntersectionPoint()";
  
  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
  G4bool        validIntersectP= true;   // Is it current ?
  G4double      NewSafety = 0.0;
  G4bool        last_AF_intersection = false;   
 
  auto    integrDriver = GetChordFinderFor()->GetIntegrationDriver();
  G4bool  driverReIntegrates = integrDriver->DoesReIntegrate();
  
  G4bool first_section = true;
  recalculatedEndPoint = false; 
  G4bool restoredFullEndpoint = false;
 
  unsigned int substep_no = 0;
 
  // Statistics for substeps
  static G4ThreadLocal unsigned int max_no_seen= 0;
 
#ifdef G4DEBUG_FIELD
  unsigned int trigger_substepno_print = 0;
  const G4double tolerance = 1.0e-8 * CLHEP::mm;
  unsigned int biggest_depth = 0;
    // using kInitialisingCL = G4LocatorChangeRecord::kInitialisingCL;
#endif
 
  // Log the location, iteration of changes in A,B 
  //----------------------------------------------
  static G4ThreadLocal G4LocatorChangeLogger endChangeA("StartPointA"),
    endChangeB("EndPointB"), recApproxPoint("ApproxPoint"),
    pointH_logger("Trial points 'E': position, normal");
  unsigned int eventCount = 0;
 
  if (fCheckMode)
  {
    // Clear previous call's data
    endChangeA.clear();
    endChangeB.clear();
    recApproxPoint.clear();
    pointH_logger.clear();
 
    // Record the initialisation
    ++eventCount;
    endChangeA.AddRecord( G4LocatorChangeRecord::kInitialisingCL, substep_no,
                          eventCount, CurrentA_PointVelocity );
    endChangeB.AddRecord( G4LocatorChangeRecord::kInitialisingCL, substep_no,
                          eventCount, CurrentB_PointVelocity );
  }
 
  //--------------------------------------------------------------------------  
  //  Algorithm for the case if progress in founding intersection is too slow.
  //  Process is defined too slow if after N=param_substeps advances on the
  //  path, it will be only 'fraction_done' of the total length.
  //  In this case the remaining length is divided in two half and 
  //  the loop is restarted for each half.  
  //  If progress is still too slow, the division in two halfs continue
  //  until 'max_depth'.
  //--------------------------------------------------------------------------
 
  const G4int param_substeps = 5;  // Test value for the maximum number
                                   // of substeps
  const G4double fraction_done = 0.3;
 
  G4bool Second_half = false;    // First half or second half of divided step
 
  // We need to know this for the 'final_section':
  // real 'final_section' or first half 'final_section'
  // In algorithm it is considered that the 'Second_half' is true
  // and it becomes false only if we are in the first-half of level
  // depthness or if we are in the first section
 
  unsigned int depth = 0; // Depth counts subdivisions of initial step made
  ++fNumCalls;
  
  NormalAtEntry = GetSurfaceNormal(CurrentE_Point, validNormalAtE);
 
  if (fCheckMode)
  {
    pointH_logger.AddRecord( G4LocatorChangeRecord::kInitialisingCL,
                             substep_no, eventCount,
                             G4FieldTrack( CurrentE_Point,0.,NormalAtEntry,0.,
                                           0., 1.,G4ThreeVector(0.),0.,0.) );
  #if (G4DEBUG_FIELD>1) 
    G4ThreeVector  StartPosition = CurveStartPointVelocity.GetPosition(); 
    if( (TrialPoint - StartPosition).mag2() < sqr(tolerance)) 
    {
       ReportImmediateHit( MethodName, StartPosition, TrialPoint,
                           tolerance, fNumCalls); 
    }
  #endif
  }
 
  // Intermediates Points on the Track = Subdivided Points must be stored.
  // Give the initial values to 'InterMedFt'
  // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
  //
  *ptrInterMedFT[0] = CurveEndPointVelocity;
  for ( auto idepth=1; idepth<max_depth+1; ++idepth )
  {
    *ptrInterMedFT[idepth] = CurveStartPointVelocity;
  }
 
  // Final_section boolean store
  //
  G4bool fin_section_depth[max_depth];
  for ( auto idepth=0; idepth<max_depth; ++idepth )
  {
    fin_section_depth[idepth] = true;
  }
  // 'SubStartPoint' is needed to calculate the length of the divided step
  //
  G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
   
  do  // Loop checking, 07.10.2016, J.Apostolakis
  {
    unsigned int substep_no_p = 0;
    G4bool sub_final_section = false; // the same as final_section,
                                      // but for 'sub_section'
    SubStart_PointVelocity = CurrentA_PointVelocity;
 
    do // Loop checking, 07.10.2016, J.Apostolakis
    { // REPEAT param
 
#ifdef G4DEBUG_FIELD       
      if( CurrentA_PointVelocity.GetCurveLength() >=
          CurrentB_PointVelocity.GetCurveLength() )
      {
        G4cerr << "ERROR> (Start) Point A coincides with or has gone past (end) point B"           
               << "MLL: iters = " << substep_no << G4endl;
     // G4LocatorChangeRecord::ReportVector(G4cerr, "endPointB", endChangeB );
     // G4cerr<<"EndPoints A(start) and B(end): combined changes " << G4endl;
        G4LocatorChangeLogger::ReportEndChanges(G4cerr, endChangeA, endChangeB);
      }
#endif    
      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      
      recApproxPoint.push_back(G4LocatorChangeRecord(G4LocatorChangeRecord::kInvalidCL,
                                                  substep_no, eventCount, ApproxIntersecPointV ) );
      G4double lenIntsc= ApproxIntersecPointV.GetCurveLength();
      G4double lenB    = CurrentB_PointVelocity.GetCurveLength();
      G4double checkVsEnd= lenB - lenIntsc;
 
      if( lenIntsc > lenB ) 
      {
         std::ostringstream errmsg;
         errmsg.precision(17);
         G4double ratio    = checkVsEnd / lenB;         
         G4double ratioTol = std::fabs(ratio) / tolerance;
         errmsg << "Intermediate F point is past end B point" << G4endl
                << "   l( intersection ) = " << lenIntsc << G4endl
                << "   l( endpoint     ) = " << lenB     << G4endl;
         errmsg.precision(8);
         errmsg << "   l_end - l_inters  = " << checkVsEnd << G4endl
                << "          / l_end      = " << ratio << G4endl
                << "   ratio  / tolerance  = " << ratioTol  << G4endl;
         if( ratioTol < 1.0 ) 
            G4Exception(MethodName, "GeomNav0003", JustWarning, errmsg );
         else
            G4Exception(MethodName, "GeomNav0003", FatalException, errmsg );
      }
#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 );
 
#ifdef G4DEBUG_FIELD
      if( fVerboseLevel > 3 )
      { 
         G4ThreeVector  ChordAB           = Point_B - Point_A;
         G4double       ABchord_length    = ChordAB.mag(); 
         G4double       MomDir_dot_ABchord;
         MomDir_dot_ABchord = (1.0 / ABchord_length)
                            * NewMomentumDir.dot( ChordAB );
         G4VIntersectionLocator::ReportTrialStep( substep_no, ChordAB,
              ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE ); 
         G4cout << " dot( MomentumDir, ABchord_unit ) = "
                << MomDir_dot_ABchord << G4endl;
      }
#endif
      G4bool adequate_angle =
             ( MomDir_dot_Norm >= 0.0 ) // Can use ( > -epsilon) instead
          || (! validNormalAtE );       // Invalid, cannot use this criterion
      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,
                                    previousSafety, previousSftOrigin );
          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 Intersects_FB = false;
        G4bool Intersects_AF = IntersectChord( Point_A,   CurrentF_Point,
                                               NewSafety, previousSafety,
                                               previousSftOrigin,
                                               stepLengthAF,
                                               PointG );
        last_AF_intersection = Intersects_AF;
        if( Intersects_AF )
        {
          // G is our new Candidate for the intersection point.
          // It replaces  "E" and we will see if it's good enough.
          CurrentB_PointVelocity = ApproxIntersecPointV;  // B  <- F
          CurrentE_Point = PointG;                        // E  <- G
 
          validIntersectP = true;    // 'E' has been updated.
 
          G4bool validNormalLast; 
          NormalAtEntry  = GetSurfaceNormal( PointG, validNormalLast ); 
          validNormalAtE = validNormalLast; 
 
          // As we move point B, must take care in case the current
          // AF has no intersection to try current FB!!
          fin_section_depth[depth] = false;
 
          if (fCheckMode)
          {
            ++eventCount;
            endChangeB.push_back(
              G4LocatorChangeRecord(G4LocatorChangeRecord::kIntersectsAF,
                             substep_no, eventCount, CurrentB_PointVelocity) );
          }
#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;
 
          // Check whether any volumes are encountered by the chord FB
          // ---------------------------------------------------------
 
          Intersects_FB = IntersectChord( CurrentF_Point, Point_B, 
                                          NewSafety, previousSafety,
                                          previousSftOrigin,
                                          stepLengthFB,
                                          PointH );
          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;
 
            validIntersectP = true;    // 'E' has been updated.
 
            G4bool validNormalH;
            NormalAtEntry  = GetSurfaceNormal( PointH, validNormalH ); 
            validNormalAtE = validNormalH;
 
            if (fCheckMode)
            {
              ++eventCount;
              endChangeA.push_back( 
                 G4LocatorChangeRecord(G4LocatorChangeRecord::kIntersectsFB,
                             substep_no, eventCount, CurrentA_PointVelocity) );
              G4FieldTrack intersectH_pn('0');  // Point and normal
                                                // nothing else will be valid
              intersectH_pn.SetPosition( PointH );
              intersectH_pn.SetMomentum( NormalAtEntry ); 
              pointH_logger.AddRecord(G4LocatorChangeRecord::kIntersectsFB,
                                    substep_no, eventCount, intersectH_pn );
            }
          }
          else  // not Intersects_FB
          {
            if( fin_section_depth[depth] )
            {
              // 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 
 
              // Check on real final_section or SubEndSection
              //
              if( ((Second_half)&&(depth==0)) || (first_section) )
              {
                there_is_no_intersection = true;   // real final_section
              }
              else
              {
                // end of subsection, not real final section 
                // exit from the and go to the depth-1 level 
                substep_no_p = param_substeps+2;  // exit from the loop
 
                // but 'Second_half' is still true because we need to find
                // the 'CurrentE_point' for the next loop
                Second_half = true; 
                sub_final_section = true;
              }
            }
            else
            {
              CurrentA_PointVelocity = CurrentB_PointVelocity;  // Got to B
              CurrentB_PointVelocity = (depth==0) ? CurveEndPointVelocity 
                                                  : *ptrInterMedFT[depth] ;
              SubStart_PointVelocity = CurrentA_PointVelocity;
              restoredFullEndpoint = true;
 
              validIntersectP = false;  // 'E' has NOT been updated.
 
              if (fCheckMode)
              {
                ++eventCount;
                endChangeA.push_back(
                  G4LocatorChangeRecord(
                    G4LocatorChangeRecord::kNoIntersectAForFB,
                             substep_no, eventCount, CurrentA_PointVelocity) ); 
                endChangeB.push_back(
                  G4LocatorChangeRecord (
                    G4LocatorChangeRecord::kNoIntersectAForFB,
                             substep_no, eventCount, CurrentB_PointVelocity) );
              }
            }
          } // Endif (Intersects_FB)
        } // Endif (Intersects_AF)
 
        G4int errorEndPt = 0; // Default: no error (if not calling CheckAnd...
 
        G4bool recalculatedB= false;
        if( driverReIntegrates )
        {
          G4FieldTrack RevisedB_FT = CurrentB_PointVelocity;
          recalculatedB= CheckAndReEstimateEndpoint(CurrentA_PointVelocity,
                                                    CurrentB_PointVelocity,
                                                    RevisedB_FT,
                                                    errorEndPt );
          if( recalculatedB )
          {
            CurrentB_PointVelocity = RevisedB_FT;  // Use it !
            // Do not invalidate intersection F -- it is still roughly OK.
            //
            // The best course would be to invalidate (reset validIntersectP)
            // BUT if we invalidate it, we must re-estimate it somewhere! E.g.
            //   validIntersectP = false;  // 'E' has NOT been updated.
 
            if ( (fin_section_depth[depth])           // real final section
               &&( first_section || ((Second_half)&&(depth==0)) ) )
            {
              recalculatedEndPoint = true;
              IntersectedOrRecalculatedFT = RevisedB_FT;
              // So that we can return it, if it is the endpoint!
            }
            // else
            //  Move forward the other points 
            //   - or better flag it, so that they are re-computed when next used
            //     [ Implementation: a counter for # of recomputations
            //       => avoids extra work]
          }
          // else
          //  Move forward the other points 
          //   - or better flag it, so that they are re-computed when next used
          //     [ Implementation: a counter for # of recomputations
          //       => avoids extra work]
          if (fCheckMode)
          {
            ++eventCount;
            endChangeB.push_back(
              G4LocatorChangeRecord( G4LocatorChangeRecord::kRecalculatedB,
                                     substep_no, eventCount, RevisedB_FT ) );
          }
        }
        else
        {
          if( CurrentB_PointVelocity.GetCurveLength() < CurrentA_PointVelocity.GetCurveLength() )
            errorEndPt = 2;
        }
        
        if( errorEndPt > 1 )  // errorEndPt = 1 is milder, just: len(B)=len(A)
        {
           std::ostringstream errmsg; 
           ReportReversedPoints(errmsg, 
                     CurveStartPointVelocity, CurveEndPointVelocity,
                     NewSafety, fiEpsilonStep, 
                     CurrentA_PointVelocity, CurrentB_PointVelocity,
                     SubStart_PointVelocity, CurrentE_Point,
                     ApproxIntersecPointV, substep_no, substep_no_p, depth);
           errmsg << G4endl << " * Location: " << MethodName
                  << "- After EndIf(Intersects_AF)" << G4endl;
           errmsg << " * Bool flags:  Recalculated = " << recalculatedB
                  << "   Intersects_AF = " << Intersects_AF
                  << "   Intersects_FB = " << Intersects_FB << G4endl;
           errmsg << " * Number of calls to MLL:EIP= " << fNumCalls << G4endl;
           G4Exception(MethodName, "GeomNav0003", FatalException, errmsg); 
        }
        if( restoredFullEndpoint )
        {
          fin_section_depth[depth] = restoredFullEndpoint;
          restoredFullEndpoint = false;
        }
      } // EndIf ( E is close enough to the curve, ie point F. )
        // tests ChordAF_Vector.mag() <= maximum_lateral_displacement 
 
#ifdef G4DEBUG_FIELD
      if( trigger_substepno_print == 0)
      {
        trigger_substepno_print= fWarnSteps - 20;
      }
 
      if( substep_no >= trigger_substepno_print )
      {
        G4cout << "Difficulty in converging in " << MethodName
               << "    Substep no = " << substep_no << G4endl;
        if( substep_no == trigger_substepno_print )
        {
          G4cout << " Start:   ";
          printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
                       -1.0, NewSafety, 0 );
 
          G4cout << " ** Change records: " << G4endl;
          G4cout <<  "endPoints  A (start) and B (end): combined changes of AB intervals" << G4endl;          
          G4LocatorChangeRecord::ReportEndChanges(G4cout, endChangeA, endChangeB );
        }
        G4cout << " Point A: "; 
        printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
                     -1.0, NewSafety, substep_no-1 );
        G4cout << " Point B: ";
        printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
                     -1.0, NewSafety, substep_no );
      }
#endif
      ++substep_no; 
      ++substep_no_p;
 
    } while (  ( ! found_approximate_intersection )
            && ( ! there_is_no_intersection )     
            && ( substep_no_p <= param_substeps) );  // UNTIL found or
                                                     // failed param substep
 
    if( (!found_approximate_intersection) && (!there_is_no_intersection) )
    {
      G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
                       - SubStart_PointVelocity.GetCurveLength()); 
      G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
                       - SubStart_PointVelocity.GetCurveLength());
   
      G4double distAB = -1;
      //
      // Is progress is too slow, and is it possible to go deeper?
      // If so, then *halve the step*
      //              ==============
      if(   (did_len < fraction_done*all_len)
         && (depth<max_depth) && (!sub_final_section) )
      {
#ifdef G4DEBUG_FIELD         
        static G4ThreadLocal unsigned int numSplits=0;   // For debugging only 
        biggest_depth = std::max(depth, biggest_depth);
        ++numSplits;
#endif
        Second_half = false;
        ++depth;
        first_section = false;
 
        G4double Sub_len = (all_len-did_len)/(2.);
        G4FieldTrack midPoint = CurrentA_PointVelocity;
        G4bool fullAdvance=              
           integrDriver->AccurateAdvance(midPoint, Sub_len, fiEpsilonStep);
                         
        ++fNumAdvanceTrials;
        if( fullAdvance )  { ++fNumAdvanceFull; }
 
        G4double lenAchieved=
           midPoint.GetCurveLength()-CurrentA_PointVelocity.GetCurveLength();
        
        const G4double adequateFraction = (1.0-CLHEP::perThousand);
        G4bool goodAdvance = (lenAchieved >= adequateFraction * Sub_len);
        if ( goodAdvance )  { ++fNumAdvanceGood; }
 
#ifdef G4DEBUG_FIELD
        else  //  !goodAdvance
        {
           G4cout << "MLL> AdvanceChordLimited not full at depth=" << depth 
                  << "  total length achieved = " << lenAchieved << " of "
                  << Sub_len << "  fraction= ";
           if (Sub_len != 0.0 ) { G4cout << lenAchieved / Sub_len; }
           else                 { G4cout << "DivByZero"; }
           G4cout << "  Good-enough fraction = " << adequateFraction;
           G4cout << "  Number of call to mll = " << fNumCalls
                  << "   iteration " << substep_no
                  << "  inner =  " << substep_no_p << G4endl;
           G4cout << "  Epsilon of integration = " << fiEpsilonStep << G4endl;
           G4cout << "  State at start is      = " << CurrentA_PointVelocity
                  << G4endl
                  << "        at end (midpoint)= " << midPoint << G4endl;
           G4cout << "  Particle mass = " << midPoint.GetRestMass() << G4endl;
 
           G4EquationOfMotion *equation = integrDriver->GetEquationOfMotion();
           ReportFieldValue( CurrentA_PointVelocity, "start", equation );
           ReportFieldValue( midPoint, "midPoint", equation );            
           G4cout << "  Original Start = "
                  << CurveStartPointVelocity << G4endl;
           G4cout << "  Original   End = "
                  << CurveEndPointVelocity << G4endl;
           G4cout << "  Original TrialPoint = "
                  << TrialPoint << G4endl;
           G4cout << "  (this is STRICT mode) "
                  << "  num Splits= " << numSplits;           
           G4cout << G4endl;
        }
#endif
 
        *ptrInterMedFT[depth] =  midPoint;
        CurrentB_PointVelocity = midPoint; 
 
        if (fCheckMode)
        {
          ++eventCount;
          endChangeB.push_back(
            G4LocatorChangeRecord( G4LocatorChangeRecord::kInsertingMidPoint,
                                   substep_no, eventCount, midPoint) );
        }        
        
        // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
        //
        SubStart_PointVelocity = CurrentA_PointVelocity;
 
        // Find new trial intersection point needed at start of the loop
        //
        G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
        G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();   
 
        G4ThreeVector PointGe;
        GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
        G4bool Intersects_AB = IntersectChord(Point_A, Point_B,
                                              NewSafety, previousSafety,
                                              previousSftOrigin, distAB,
                                              PointGe);
        if( Intersects_AB )
        {
          last_AF_intersection = Intersects_AB;
          CurrentE_Point = PointGe;
          fin_section_depth[depth] = true;
 
          validIntersectP = true;  // 'E' has been updated.
 
          G4bool validNormalAB; 
          NormalAtEntry  = GetSurfaceNormal( PointGe, validNormalAB ); 
          validNormalAtE = validNormalAB;
        }
        else
        {
          // No intersection found for first part of curve
          // (CurrentA,InterMedPoint[depth]). Go to the second part
          //
          Second_half = true;
 
          validIntersectP=  false;  // No new 'E' chord intersection found
        }
      } // if did_len
 
      unsigned int levelPops = 0;
      
      G4bool unfinished = Second_half;
      while ( unfinished && (depth>0) )  // Loop checking, 07.10.2016, JA
      {
        // Second part of curve (InterMed[depth],Intermed[depth-1])) 
        // On the depth-1 level normally we are on the 'second_half'
 
        ++levelPops;
 
        //  Find new trial intersection point needed at start of the loop
        //
        SubStart_PointVelocity = *ptrInterMedFT[depth];
        CurrentA_PointVelocity = *ptrInterMedFT[depth];
        CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
 
        if (fCheckMode)
        {
          ++eventCount;
          G4LocatorChangeRecord chngPop_a( G4LocatorChangeRecord::kLevelPop,
                              substep_no, eventCount, CurrentA_PointVelocity);
          endChangeA.push_back( chngPop_a );
          G4LocatorChangeRecord chngPop_b( G4LocatorChangeRecord::kLevelPop,
                              substep_no, eventCount, CurrentB_PointVelocity);
          endChangeB.push_back( chngPop_b );
        }
        
        // Ensure that the new endpoints are not further apart in space
        // than on the curve due to different errors in the integration
        //
        G4int  errorEndPt = -1;
        G4bool recalculatedB= false;
        if( driverReIntegrates )
        {
          // Ensure that the new endpoints are not further apart in space
          // than on the curve due to different errors in the integration
          //
          G4FieldTrack RevisedEndPointFT = CurrentB_PointVelocity;
          recalculatedB = 
             CheckAndReEstimateEndpoint( CurrentA_PointVelocity,  
                                         CurrentB_PointVelocity,
                                         RevisedEndPointFT,
                                         errorEndPt );
          if( recalculatedB )
          {
            CurrentB_PointVelocity = RevisedEndPointFT;  // Use it !
            
            if ( depth == 1 )
            {
               recalculatedEndPoint = true;
               IntersectedOrRecalculatedFT = RevisedEndPointFT;
               // So that we can return it, if it is the endpoint!
            }
          }
          else
          {          
            if( CurrentB_PointVelocity.GetCurveLength() < CurrentA_PointVelocity.GetCurveLength() )
              errorEndPt = 2;
          }
          
          if (fCheckMode)
          {
            ++eventCount;
            endChangeB.push_back(
              G4LocatorChangeRecord(G4LocatorChangeRecord::kRecalculatedBagn,
                                    substep_no, eventCount, RevisedEndPointFT));
          }
        }
        if( errorEndPt > 1 )  // errorEndPt = 1 is milder, just: len(B)=len(A)
        {
          std::ostringstream errmsg; 
          ReportReversedPoints(errmsg, 
                    CurveStartPointVelocity, CurveEndPointVelocity,
                    NewSafety, fiEpsilonStep, 
                    CurrentA_PointVelocity, CurrentB_PointVelocity,
                    SubStart_PointVelocity, CurrentE_Point,
                    ApproxIntersecPointV, substep_no, substep_no_p, depth);
          errmsg << " * Location:  " << MethodName << "- Second-Half" << G4endl;
          errmsg << " * Recalculated = " << recalculatedB << G4endl; // false
          G4Exception(MethodName, "GeomNav0003", FatalException, errmsg);
        }
        
        G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
        G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
        G4ThreeVector PointGi;     
        GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
        G4bool Intersects_AB = IntersectChord(Point_A, Point_B, NewSafety,
                                              previousSafety,
                                              previousSftOrigin, distAB,
                                              PointGi);
        if( Intersects_AB )
        {
          last_AF_intersection = Intersects_AB;
          CurrentE_Point = PointGi;
 
          validIntersectP = true;  // 'E' has been updated.
          NormalAtEntry  = GetSurfaceNormal( PointGi, validNormalAtE ); 
        }
        else
        {
          validIntersectP = false;  // No new 'E' chord intersection found
          if( depth == 1)
          {
            there_is_no_intersection = true;
          }
        }
        depth--;
        fin_section_depth[depth] = true;
        unfinished = !validIntersectP;
      }
#ifdef G4DEBUG_FIELD
      if( ! ( validIntersectP || there_is_no_intersection ) )
      {
         // What happened ??
         G4cout << "MLL - WARNING Potential FAILURE: Conditions not met!"
                << G4endl
                << " Depth = " << depth << G4endl
                << " Levels popped = " << levelPops
                << " Num Substeps= " << substep_no << G4endl;
         G4cout << " Found intersection= " << found_approximate_intersection
                << G4endl;
         G4cout << " Progress report: -- " << G4endl;
         ReportProgress(G4cout, 
                        CurveStartPointVelocity, CurveEndPointVelocity,
                        substep_no, CurrentA_PointVelocity,
                        CurrentB_PointVelocity,
                        NewSafety, depth ); 
         G4cout << G4endl; 
      }
#endif      
    }  // if(!found_aproximate_intersection)
 
    assert( validIntersectP || there_is_no_intersection
                            || found_approximate_intersection);
 
  } while ( ( ! found_approximate_intersection )
            && ( ! there_is_no_intersection )     
            && ( substep_no <= fMaxSteps) ); // UNTIL found or failed
 
  if( substep_no > max_no_seen )
  {
    max_no_seen = substep_no; 
#ifdef G4DEBUG_FIELD
    if( max_no_seen > fWarnSteps )
    {
      trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps 
    } 
#endif
  }
 
  if( !there_is_no_intersection && !found_approximate_intersection )
  {
     if( substep_no >= fMaxSteps)
     {
        // Since we cannot go further (yet), we return as far as we have gone
 
        recalculatedEndPoint = true;
        IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
        found_approximate_intersection = false; 
     
        std::ostringstream message;
        message << G4endl;
        message << "Convergence is requiring too many substeps: "
                << substep_no << "  ( limit = "<<  fMaxSteps << ")"
                << G4endl 
                << "  Abandoning effort to intersect. " << G4endl << G4endl;
        message << "    Number of calls to MLL: " <<   fNumCalls;
        message << "  iteration = " << substep_no <<G4endl << G4endl;
        
        message.precision( 12 ); 
        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; 
        
        message << "     Start and end-point of requested Step:" << G4endl;
        printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
                     -1.0, NewSafety, 0,     message, -1 );
        message << "        Start and end-point of current Sub-Step:" << G4endl;
        printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
                     -1.0, NewSafety, substep_no-1, message, -1 );
        printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
                     -1.0, NewSafety, substep_no, message, -1 );
        
        G4Exception(MethodName, "GeomNav0003", JustWarning, message);
     }
     else if( substep_no >= fWarnSteps)
     {  
        std::ostringstream message;
        message << "Many substeps while trying to locate intersection."
                << G4endl
                << "          Undertaken length: "  
                << CurrentB_PointVelocity.GetCurveLength()
                << " - Needed: "  << substep_no << " substeps." << G4endl
                << "          Warning number = " << fWarnSteps
                << " and maximum substeps = " << fMaxSteps;
        G4Exception(MethodName, "GeomNav1002", JustWarning, message);
     }
  }
  
  return  (!there_is_no_intersection) && found_approximate_intersection;
    //  Success or failure
}

ReportStatistics()

void G4MultiLevelLocator::ReportStatistics

(

)

Definition at line 969 of file G4MultiLevelLocator.cc.

{
   G4cout << " Number of calls = " << fNumCalls << G4endl;
   G4cout << " Number of split level ('advances'):  "
          << fNumAdvanceTrials << G4endl;
   G4cout << " Number of full advances:             "
          << fNumAdvanceGood << G4endl;
   G4cout << " Number of good advances:             "
          << fNumAdvanceFull << G4endl;
}

Referenced by ~G4MultiLevelLocator().

SetMaxSteps()

void G4MultiLevelLocator::SetMaxSteps ( unsigned int  valMax )

inline

Definition at line 71 of file G4MultiLevelLocator.hh.

{ fMaxSteps= valMax; }

Referenced by G4MultiLevelLocator().

SetWarnSteps()

void G4MultiLevelLocator::SetWarnSteps ( unsigned int  valWarn )

inline

Definition at line 72 of file G4MultiLevelLocator.hh.

{ fWarnSteps= valWarn; }

Referenced by G4MultiLevelLocator().

---

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

• G4MultiLevelLocator.hh
• G4MultiLevelLocator.cc