Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/navigation/src/G4SimpleLocator.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/navigation/src/G4SimpleLocator.cc (Version 11.3.0) and /geometry/navigation/src/G4SimpleLocator.cc (Version 9.6.p3)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id$
                                                   >>  27 //
 26 // Class G4SimpleLocator implementation            28 // Class G4SimpleLocator implementation
 27 //                                                 29 //
 28 // 27.10.08 - Tatiana Nikitina, extracted from     30 // 27.10.08 - Tatiana Nikitina, extracted from G4PropagatorInField class
 29 // 04.10.11 - John Apostolakis, revised conver     31 // 04.10.11 - John Apostolakis, revised convergence to use Surface Normal
 30 // -------------------------------------------     32 // ---------------------------------------------------------------------------
 31                                                    33 
 32 #include <iomanip>                                 34 #include <iomanip>
 33                                                    35 
 34 #include "G4ios.hh"                                36 #include "G4ios.hh"
 35 #include "G4SimpleLocator.hh"                      37 #include "G4SimpleLocator.hh"
 36                                                    38 
 37 G4SimpleLocator::G4SimpleLocator(G4Navigator *     39 G4SimpleLocator::G4SimpleLocator(G4Navigator *theNavigator)
 38   : G4VIntersectionLocator(theNavigator)           40   : G4VIntersectionLocator(theNavigator)
 39 {                                                  41 {
 40 }                                                  42 }      
 41                                                    43 
 42 G4SimpleLocator::~G4SimpleLocator() = default; <<  44 G4SimpleLocator::~G4SimpleLocator()
                                                   >>  45 {
                                                   >>  46 }
 43                                                    47 
 44 // -------------------------------------------     48 // --------------------------------------------------------------------------
 45 // G4bool G4PropagatorInField::LocateIntersect     49 // G4bool G4PropagatorInField::LocateIntersectionPoint( 
 46 //        const G4FieldTrack&       CurveStart     50 //        const G4FieldTrack&       CurveStartPointVelocity,   // A
 47 //        const G4FieldTrack&       CurveEndPo     51 //        const G4FieldTrack&       CurveEndPointVelocity,     // B
 48 //        const G4ThreeVector&      TrialPoint     52 //        const G4ThreeVector&      TrialPoint,                // E
 49 //              G4FieldTrack&       Intersecte     53 //              G4FieldTrack&       IntersectedOrRecalculated  // Output
 50 //              G4bool&             recalculat     54 //              G4bool&             recalculated )             // Out
 51 // -------------------------------------------     55 // --------------------------------------------------------------------------
 52 //                                                 56 //
 53 // Function that returns the intersection of t     57 // Function that returns the intersection of the true path with the surface
 54 // of the current volume (either the external      58 // of the current volume (either the external one or the inner one with one
 55 // of the daughters:                               59 // of the daughters:
 56 //                                                 60 //
 57 //     A = Initial point                           61 //     A = Initial point
 58 //     B = another point                           62 //     B = another point 
 59 //                                                 63 //
 60 // Both A and B are assumed to be on the true      64 // Both A and B are assumed to be on the true path:
 61 //                                                 65 //
 62 //     E is the first point of intersection of     66 //     E is the first point of intersection of the chord AB with 
 63 //     a volume other than A  (on the surface      67 //     a volume other than A  (on the surface of A or of a daughter)
 64 //                                                 68 //
 65 // Convention of Use :                             69 // Convention of Use :
 66 //     i) If it returns "true", then Intersect     70 //     i) If it returns "true", then IntersectionPointVelocity is set
 67 //       to the approximate intersection point     71 //       to the approximate intersection point.
 68 //    ii) If it returns "false", no intersecti     72 //    ii) If it returns "false", no intersection was found.
 69 //          The validity of IntersectedOrRecal     73 //          The validity of IntersectedOrRecalculated depends on 'recalculated'
 70 //        a) if latter is false, then Intersec     74 //        a) if latter is false, then IntersectedOrRecalculated is invalid. 
 71 //        b) if latter is true,  then Intersec     75 //        b) if latter is true,  then IntersectedOrRecalculated is
 72 //             the new endpoint, due to a re-i     76 //             the new endpoint, due to a re-integration.
 73 // -------------------------------------------     77 // --------------------------------------------------------------------------
 74 // NOTE: implementation taken from G4Propagato     78 // NOTE: implementation taken from G4PropagatorInField
 75 //                                                 79 //
 76 G4bool G4SimpleLocator::EstimateIntersectionPo     80 G4bool G4SimpleLocator::EstimateIntersectionPoint( 
 77          const  G4FieldTrack&       CurveStart     81          const  G4FieldTrack&       CurveStartPointVelocity,     // A
 78          const  G4FieldTrack&       CurveEndPo     82          const  G4FieldTrack&       CurveEndPointVelocity,       // B
 79          const  G4ThreeVector&      TrialPoint     83          const  G4ThreeVector&      TrialPoint,                  // E
 80                 G4FieldTrack&       Intersecte     84                 G4FieldTrack&       IntersectedOrRecalculatedFT, // Output
 81                 G4bool&             recalculat     85                 G4bool&             recalculatedEndPoint,        // Out
 82                 G4double            &fPrevious     86                 G4double            &fPreviousSafety,            //In/Out
 83                 G4ThreeVector       &fPrevious     87                 G4ThreeVector       &fPreviousSftOrigin)         //In/Out
 84 {                                                  88 {
 85   // Find Intersection Point ( A, B, E )  of t     89   // Find Intersection Point ( A, B, E )  of true path AB - start at E.
 86                                                    90 
 87   G4bool found_approximate_intersection = fals     91   G4bool found_approximate_intersection = false;
 88   G4bool there_is_no_intersection       = fals     92   G4bool there_is_no_intersection       = false;
 89                                                    93   
 90   G4FieldTrack  CurrentA_PointVelocity = Curve     94   G4FieldTrack  CurrentA_PointVelocity = CurveStartPointVelocity; 
 91   G4FieldTrack  CurrentB_PointVelocity = Curve     95   G4FieldTrack  CurrentB_PointVelocity = CurveEndPointVelocity;
 92   G4ThreeVector CurrentE_Point = TrialPoint;       96   G4ThreeVector CurrentE_Point = TrialPoint;
 93   G4bool        validNormalAtE = false;            97   G4bool        validNormalAtE = false;
 94   G4ThreeVector NormalAtEntry;                     98   G4ThreeVector NormalAtEntry;
 95                                                    99 
 96   G4FieldTrack  ApproxIntersecPointV(CurveEndP    100   G4FieldTrack  ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
 97   G4double      NewSafety = 0.0;                  101   G4double      NewSafety = 0.0;
 98   G4bool last_AF_intersection = false;            102   G4bool last_AF_intersection = false;
 99   G4bool final_section = true;  // Shows wheth    103   G4bool final_section = true;  // Shows whether current section is last
100                                 // (i.e. B=ful    104                                 // (i.e. B=full end)
101   recalculatedEndPoint = false;                   105   recalculatedEndPoint = false; 
102                                                   106 
103   G4bool restoredFullEndpoint = false;            107   G4bool restoredFullEndpoint = false;
104                                                   108 
105   G4int substep_no = 0;                           109   G4int substep_no = 0;
106                                                   110 
107   // Limits for substep number                    111   // Limits for substep number
108   //                                              112   //
109   const G4int max_substeps  = 100000000;  // T    113   const G4int max_substeps  = 100000000;  // Test 120  (old value 100 )
110   const G4int warn_substeps = 1000;       //      114   const G4int warn_substeps = 1000;       //      100  
111                                                   115 
112   // Statistics for substeps                      116   // Statistics for substeps
113   //                                              117   //
114   static G4ThreadLocal G4int max_no_seen= -1;  << 118   static G4int max_no_seen= -1; 
115                                                   119 
116   NormalAtEntry = GetSurfaceNormal( CurrentE_P    120   NormalAtEntry = GetSurfaceNormal( CurrentE_Point, validNormalAtE); 
117                                                   121 
118 #ifdef G4DEBUG_FIELD                              122 #ifdef G4DEBUG_FIELD
119   const G4double tolerance = 1.0e-8;           << 123   static G4double tolerance = 1.0e-8; 
120   G4ThreeVector  StartPosition= CurveStartPoin    124   G4ThreeVector  StartPosition= CurveStartPointVelocity.GetPosition(); 
121   if( (TrialPoint - StartPosition).mag() < tol << 125   if( (TrialPoint - StartPosition).mag() < tolerance * mm ) 
122   {                                               126   {
123      G4Exception("G4SimpleLocator::EstimateInt    127      G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 
124                  "GeomNav1002", JustWarning,      128                  "GeomNav1002", JustWarning,
125                  "Intersection point F is exac    129                  "Intersection point F is exactly at start point A." ); 
126   }                                               130   }
127 #endif                                            131 #endif
128                                                   132 
129   do                                              133   do 
130   {                                               134   {
131      G4ThreeVector Point_A = CurrentA_PointVel    135      G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();  
132      G4ThreeVector Point_B = CurrentB_PointVel    136      G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
133                                                   137        
134      // F = a point on true AB path close to p    138      // F = a point on true AB path close to point E 
135      // (the closest if possible)                 139      // (the closest if possible)
136      //                                           140      //
137      ApproxIntersecPointV = GetChordFinderFor(    141      ApproxIntersecPointV = GetChordFinderFor()
138                             ->ApproxCurvePoint    142                             ->ApproxCurvePointV( CurrentA_PointVelocity, 
139                                                   143                                                  CurrentB_PointVelocity, 
140                                                   144                                                  CurrentE_Point,
141                                                   145                                                  GetEpsilonStepFor());
142        //  The above method is the key & most     146        //  The above method is the key & most intuitive part ...
143                                                   147       
144 #ifdef G4DEBUG_FIELD                              148 #ifdef G4DEBUG_FIELD
145      if( ApproxIntersecPointV.GetCurveLength()    149      if( ApproxIntersecPointV.GetCurveLength() > 
146          CurrentB_PointVelocity.GetCurveLength    150          CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
147      {                                            151      {
148        G4Exception("G4SimpleLocator::EstimateI    152        G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 
149                    "GeomNav0003", FatalExcepti    153                    "GeomNav0003", FatalException,
150                    "Intermediate F point is pa    154                    "Intermediate F point is past end B point!" ); 
151      }                                            155      }
152 #endif                                            156 #endif
153                                                   157 
154      G4ThreeVector CurrentF_Point = ApproxInte << 158      G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
155                                                   159 
156      // First check whether EF is small - then    160      // First check whether EF is small - then F is a good approx. point 
157      // Calculate the length and direction of     161      // Calculate the length and direction of the chord AF
158      //                                           162      //
159      G4ThreeVector ChordEF_Vector = CurrentF_P << 163      G4ThreeVector  ChordEF_Vector = CurrentF_Point - CurrentE_Point;
160                                                   164 
161      G4ThreeVector NewMomentumDir = ApproxInte << 165      G4ThreeVector  NewMomentumDir= ApproxIntersecPointV.GetMomentumDir(); 
162      G4double MomDir_dot_Norm = NewMomentumDir << 166      G4double       MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ;
163                                                   167 
164      G4ThreeVector ChordAB = Point_B - Point_A << 168      G4ThreeVector  ChordAB           = Point_B - Point_A;
165                                                   169 
166 #ifdef G4DEBUG_FIELD                              170 #ifdef G4DEBUG_FIELD
167      G4VIntersectionLocator::                     171      G4VIntersectionLocator::
168        ReportTrialStep( substep_no, ChordAB, C    172        ReportTrialStep( substep_no, ChordAB, ChordEF_Vector, 
169                         NewMomentumDir, Normal << 173                       NewMomentumDir, NormalAtEntry, validNormalAtE ); 
170 #endif                                            174 #endif
171      // Check Sign is always exiting !! TODO      175      // Check Sign is always exiting !! TODO
172      // Could ( > -epsilon) be used instead?      176      // Could ( > -epsilon) be used instead?
173      //                                           177      //
174      G4bool adequate_angle = ( MomDir_dot_Norm    178      G4bool adequate_angle = ( MomDir_dot_Norm >= 0.0 ) 
175                           || (! validNormalAtE    179                           || (! validNormalAtE );  // Invalid
176      G4double EF_dist2= ChordEF_Vector.mag2();    180      G4double EF_dist2= ChordEF_Vector.mag2();
177      if ( ( EF_dist2  <= sqr(fiDeltaIntersecti    181      if ( ( EF_dist2  <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
178        || ( EF_dist2 <= kCarTolerance*kCarTole    182        || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
179      {                                            183      {
180        found_approximate_intersection = true;     184        found_approximate_intersection = true;
181                                                   185         
182        // Create the "point" return value         186        // Create the "point" return value
183        //                                         187        //
184        IntersectedOrRecalculatedFT = ApproxInt    188        IntersectedOrRecalculatedFT = ApproxIntersecPointV;
185        IntersectedOrRecalculatedFT.SetPosition    189        IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
186                                                   190 
187        if ( GetAdjustementOfFoundIntersection(    191        if ( GetAdjustementOfFoundIntersection() )
188        {                                          192        {
189          // Try to Get Correction of Intersect    193          // Try to Get Correction of IntersectionPoint using SurfaceNormal()
190          //                                       194          //  
191          G4ThreeVector IP;                        195          G4ThreeVector IP;
192          G4ThreeVector MomentumDir= ApproxInte    196          G4ThreeVector MomentumDir= ApproxIntersecPointV.GetMomentumDirection();
193          G4bool goodCorrection = AdjustmentOfF    197          G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A,
194                                    CurrentE_Po    198                                    CurrentE_Point, CurrentF_Point, MomentumDir,
195                                    last_AF_int    199                                    last_AF_intersection, IP, NewSafety,
196                                    fPreviousSa    200                                    fPreviousSafety, fPreviousSftOrigin );
197                                                   201 
198          if ( goodCorrection )                 << 202          if(goodCorrection)
199          {                                        203          {
200            IntersectedOrRecalculatedFT = Appro    204            IntersectedOrRecalculatedFT = ApproxIntersecPointV;
201            IntersectedOrRecalculatedFT.SetPosi    205            IntersectedOrRecalculatedFT.SetPosition(IP);
202          }                                        206          }
203        }                                          207        }
204                                                   208 
205        // Note: in order to return a point on     209        // Note: in order to return a point on the boundary, 
206        //       we must return E. But it is F     210        //       we must return E. But it is F on the curve.
207        //       So we must "cheat": we are usi    211        //       So we must "cheat": we are using the position at point E
208        //       and the velocity at point F !!    212        //       and the velocity at point F !!!
209        //                                         213        //
210        // This must limit the length we can al    214        // This must limit the length we can allow for displacement!
211      }                                            215      }
212      else  // E is NOT close enough to the cur    216      else  // E is NOT close enough to the curve (ie point F)
213      {                                            217      {
214        // Check whether any volumes are encoun    218        // Check whether any volumes are encountered by the chord AF
215        // ------------------------------------    219        // ---------------------------------------------------------
216        // First relocate to restore any Voxel     220        // First relocate to restore any Voxel etc information
217        // in the Navigator before calling Comp    221        // in the Navigator before calling ComputeStep()
218        //                                         222        //
219        GetNavigatorFor()->LocateGlobalPointWit    223        GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A );
220                                                   224 
221        G4ThreeVector PointG;   // Candidate in    225        G4ThreeVector PointG;   // Candidate intersection point
222        G4double stepLengthAF;                     226        G4double stepLengthAF; 
223        G4bool usedNavigatorAF = false;            227        G4bool usedNavigatorAF = false; 
224        G4bool Intersects_AF = IntersectChord(     228        G4bool Intersects_AF = IntersectChord( Point_A,   
225                                                   229                                               CurrentF_Point,
226                                                   230                                               NewSafety,
227                                                   231                                               fPreviousSafety,
228                                                   232                                               fPreviousSftOrigin,
229                                                   233                                               stepLengthAF,
230                                                   234                                               PointG,
231                                                   235                                               &usedNavigatorAF );
232        last_AF_intersection = Intersects_AF;      236        last_AF_intersection = Intersects_AF;
233        if( Intersects_AF )                        237        if( Intersects_AF )
234        {                                          238        {
235          // G is our new Candidate for the int    239          // G is our new Candidate for the intersection point.
236          // It replaces  "E" and we will repea    240          // It replaces  "E" and we will repeat the test to see if
237          // it is a good enough approximate po    241          // it is a good enough approximate point for us.
238          //       B    <- F                       242          //       B    <- F
239          //       E    <- G                       243          //       E    <- G
240                                                   244 
241          CurrentB_PointVelocity = ApproxInters    245          CurrentB_PointVelocity = ApproxIntersecPointV;
242          CurrentE_Point = PointG;                 246          CurrentE_Point = PointG;
243                                                   247 
244          // Need to recalculate the Exit Norma    248          // Need to recalculate the Exit Normal at the new PointG 
245          // Relies on a call to Navigator::Com    249          // Relies on a call to Navigator::ComputeStep in IntersectChord above
246          // If the safety was adequate (for th    250          // If the safety was adequate (for the step) this would NOT be called!
247          // But this must not occur, no inters    251          // But this must not occur, no intersection can be found in that case,
248          // so this branch, ie if( Intersects_    252          // so this branch, ie if( Intersects_AF ) would not be reached!
249          //                                       253          //
250          G4bool validNormalLast;                  254          G4bool validNormalLast; 
251          NormalAtEntry  = GetSurfaceNormal( Po    255          NormalAtEntry  = GetSurfaceNormal( PointG, validNormalLast ); 
252          validNormalAtE = validNormalLast;        256          validNormalAtE = validNormalLast; 
253                                                   257 
254          // By moving point B, must take care     258          // By moving point B, must take care if current
255          // AF has no intersection to try curr    259          // AF has no intersection to try current FB!!
256          //                                       260          //
257          final_section = false;                << 261          final_section= false; 
258                                                   262           
259 #ifdef G4VERBOSE                                  263 #ifdef G4VERBOSE
260          if( fVerboseLevel > 3 )                  264          if( fVerboseLevel > 3 )
261          {                                        265          {
262            G4cout << "G4PiF::LI> Investigating    266            G4cout << "G4PiF::LI> Investigating intermediate point"
263                   << " at s=" << ApproxInterse    267                   << " at s=" << ApproxIntersecPointV.GetCurveLength()
264                   << " on way to full s="         268                   << " on way to full s="
265                   << CurveEndPointVelocity.Get    269                   << CurveEndPointVelocity.GetCurveLength() << G4endl;
266          }                                        270          }
267 #endif                                            271 #endif
268        }                                          272        }
269        else  // not Intersects_AF                 273        else  // not Intersects_AF
270        {                                          274        {  
271          // In this case:                         275          // In this case:
272          // There is NO intersection of AF wit    276          // There is NO intersection of AF with a volume boundary.
273          // We must continue the search in the    277          // We must continue the search in the segment FB!
274          //                                       278          //
275          GetNavigatorFor()->LocateGlobalPointW    279          GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
276                                                   280 
277          G4double stepLengthFB;                   281          G4double stepLengthFB;
278          G4ThreeVector PointH;                    282          G4ThreeVector PointH;
279          G4bool usedNavigatorFB = false;       << 283          G4bool usedNavigatorFB=false; 
280                                                   284 
281          // Check whether any volumes are enco    285          // Check whether any volumes are encountered by the chord FB
282          // ----------------------------------    286          // ---------------------------------------------------------
283                                                   287 
284          G4bool Intersects_FB = IntersectChord    288          G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B, 
285                                                   289                                                 NewSafety,fPreviousSafety,
286                                                   290                                                 fPreviousSftOrigin,
287                                                   291                                                 stepLengthFB,
288                                                   292                                                 PointH, &usedNavigatorFB );
289          if( Intersects_FB )                      293          if( Intersects_FB )
290          {                                        294          { 
291            // There is an intersection of FB w    295            // There is an intersection of FB with a volume boundary
292            // H <- First Intersection of Chord    296            // H <- First Intersection of Chord FB 
293                                                   297 
294            // H is our new Candidate for the i    298            // H is our new Candidate for the intersection point.
295            // It replaces  "E" and we will rep    299            // It replaces  "E" and we will repeat the test to see if
296            // it is a good enough approximate     300            // it is a good enough approximate point for us.
297                                                   301 
298            // Note that F must be in volume vo    302            // Note that F must be in volume volA  (the same as A)
299            // (otherwise AF would meet a volum    303            // (otherwise AF would meet a volume boundary!)
300            //   A    <- F                         304            //   A    <- F 
301            //   E    <- H                         305            //   E    <- H
302            //                                     306            //
303            CurrentA_PointVelocity = ApproxInte    307            CurrentA_PointVelocity = ApproxIntersecPointV;
304            CurrentE_Point = PointH;               308            CurrentE_Point = PointH;
305                                                   309 
306            // Need to recalculate the Exit Nor    310            // Need to recalculate the Exit Normal at the new PointG
307            // Relies on call to Navigator::Com    311            // Relies on call to Navigator::ComputeStep in IntersectChord above
308            // If safety was adequate (for the     312            // If safety was adequate (for the step) this would NOT be called!
309            // But this must not occur, no inte    313            // But this must not occur, no intersection found in that case,
310            // so this branch, i.e. if( Interse    314            // so this branch, i.e. if( Intersects_AF ) would not be reached!
311            //                                     315            //
312            G4bool validNormalLast;                316            G4bool validNormalLast; 
313            NormalAtEntry  = GetSurfaceNormal(     317            NormalAtEntry  = GetSurfaceNormal( PointH, validNormalLast ); 
314            validNormalAtE = validNormalLast;      318            validNormalAtE = validNormalLast;
315          }                                        319          }
316          else  // not Intersects_FB               320          else  // not Intersects_FB
317          {                                        321          {
318            // There is NO intersection of FB w    322            // There is NO intersection of FB with a volume boundary
319                                                   323 
320            if( final_section  )                   324            if( final_section  )
321            {                                      325            {
322              // If B is the original endpoint,    326              // If B is the original endpoint, this means that whatever
323              // volume(s) intersected the orig    327              // volume(s) intersected the original chord, none touch the
324              // smaller chords we have used.      328              // smaller chords we have used.
325              // The value of 'IntersectedOrRec    329              // The value of 'IntersectedOrRecalculatedFT' returned is
326              // likely not valid                  330              // likely not valid 
327                                                   331               
328              there_is_no_intersection = true;     332              there_is_no_intersection = true;   // real final_section
329            }                                      333            }
330            else                                   334            else
331            {                                      335            {
332              // We must restore the original e    336              // We must restore the original endpoint
333                                                   337 
334              CurrentA_PointVelocity = CurrentB    338              CurrentA_PointVelocity = CurrentB_PointVelocity;  // Got to B
335              CurrentB_PointVelocity = CurveEnd    339              CurrentB_PointVelocity = CurveEndPointVelocity;
336              restoredFullEndpoint = true;         340              restoredFullEndpoint = true;
337            }                                      341            }
338          } // Endif (Intersects_FB)               342          } // Endif (Intersects_FB)
339        } // Endif (Intersects_AF)                 343        } // Endif (Intersects_AF)
340                                                   344 
341        // Ensure that the new endpoints are no    345        // Ensure that the new endpoints are not further apart in space
342        // than on the curve due to different e    346        // than on the curve due to different errors in the integration
343        //                                         347        //
344        G4double linDistSq, curveDist;             348        G4double linDistSq, curveDist; 
345        linDistSq = ( CurrentB_PointVelocity.Ge    349        linDistSq = ( CurrentB_PointVelocity.GetPosition() 
346                    - CurrentA_PointVelocity.Ge    350                    - CurrentA_PointVelocity.GetPosition() ).mag2(); 
347        curveDist = CurrentB_PointVelocity.GetC    351        curveDist = CurrentB_PointVelocity.GetCurveLength()
348                    - CurrentA_PointVelocity.Ge    352                    - CurrentA_PointVelocity.GetCurveLength();
349                                                   353 
350        // Change this condition for very stric    354        // Change this condition for very strict parameters of propagation 
351        //                                         355        //
352        if( curveDist*curveDist*(1+2* GetEpsilo    356        if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq )
353        {                                          357        {
354          // Re-integrate to obtain a new B        358          // Re-integrate to obtain a new B
355          //                                       359          //
356          G4FieldTrack newEndPointFT =             360          G4FieldTrack newEndPointFT =
357                  ReEstimateEndpoint( CurrentA_    361                  ReEstimateEndpoint( CurrentA_PointVelocity,
358                                      CurrentB_    362                                      CurrentB_PointVelocity,
359                                      linDistSq    363                                      linDistSq,    // to avoid recalculation
360                                      curveDist    364                                      curveDist );
361          G4FieldTrack oldPointVelB = CurrentB_    365          G4FieldTrack oldPointVelB = CurrentB_PointVelocity; 
362          CurrentB_PointVelocity = newEndPointF    366          CurrentB_PointVelocity = newEndPointFT;
363                                                   367 
364          if( (final_section)) // real final se    368          if( (final_section)) // real final section
365          {                                        369          {
366            recalculatedEndPoint = true;           370            recalculatedEndPoint = true;
367            IntersectedOrRecalculatedFT = newEn    371            IntersectedOrRecalculatedFT = newEndPointFT;
368              // So that we can return it, if i    372              // So that we can return it, if it is the endpoint!
369          }                                        373          }
370        }                                          374        }
371        if( curveDist < 0.0 )                      375        if( curveDist < 0.0 )
372        {                                          376        {
373          fVerboseLevel = 5; // Print out a max    377          fVerboseLevel = 5; // Print out a maximum of information
374          printStatus( CurrentA_PointVelocity,     378          printStatus( CurrentA_PointVelocity,  CurrentB_PointVelocity,
375                       -1.0, NewSafety,  subste    379                       -1.0, NewSafety,  substep_no );
376          std::ostringstream message;              380          std::ostringstream message;
377          message << "Error in advancing propag    381          message << "Error in advancing propagation." << G4endl
378                  << "        Point A (start) i    382                  << "        Point A (start) is " << CurrentA_PointVelocity
379                  << G4endl                        383                  << G4endl
380                  << "        Point B (end)   i    384                  << "        Point B (end)   is " << CurrentB_PointVelocity
381                  << G4endl                        385                  << G4endl
382                  << "        Curve distance is    386                  << "        Curve distance is " << curveDist << G4endl
383                  << G4endl                        387                  << G4endl
384                  << "The final curve point is     388                  << "The final curve point is not further along"
385                  << " than the original!" << G    389                  << " than the original!" << G4endl;
386                                                   390 
387          if( recalculatedEndPoint )               391          if( recalculatedEndPoint )
388          {                                        392          {
389            message << "Recalculation of EndPoi    393            message << "Recalculation of EndPoint was called with fEpsStep= "
390                    << GetEpsilonStepFor() << G    394                    << GetEpsilonStepFor() << G4endl;
391          }                                        395          }
392          message.precision(20);                   396          message.precision(20);
393          message << " Point A (Curve start)       397          message << " Point A (Curve start)   is " << CurveStartPointVelocity
394                  << G4endl                        398                  << G4endl
395                  << " Point B (Curve   end)       399                  << " Point B (Curve   end)   is " << CurveEndPointVelocity
396                  << G4endl                        400                  << G4endl
397                  << " Point A (Current start)     401                  << " Point A (Current start) is " << CurrentA_PointVelocity
398                  << G4endl                        402                  << G4endl
399                  << " Point B (Current end)       403                  << " Point B (Current end)   is " << CurrentB_PointVelocity
400                  << G4endl                        404                  << G4endl
401                  << " Point E (Trial Point)       405                  << " Point E (Trial Point)   is " << CurrentE_Point
402                  << G4endl                        406                  << G4endl
403                  << " Point F (Intersection)      407                  << " Point F (Intersection)  is " << ApproxIntersecPointV
404                  << G4endl                        408                  << G4endl
405                  << "        LocateIntersectio    409                  << "        LocateIntersection parameters are : Substep no= "
406                  << substep_no;                   410                  << substep_no;
407                                                   411 
408          G4Exception("G4SimpleLocator::Estimat    412          G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
409                      "GeomNav0003", FatalExcep    413                      "GeomNav0003", FatalException, message);
410        }                                          414        }
411                                                   415 
412        if ( restoredFullEndpoint )             << 416        if(restoredFullEndpoint)
413        {                                          417        {
414          final_section = restoredFullEndpoint;    418          final_section = restoredFullEndpoint;
415          restoredFullEndpoint = false;            419          restoredFullEndpoint = false;
416        }                                          420        }
417      } // EndIf ( E is close enough to the cur    421      } // EndIf ( E is close enough to the curve, ie point F. )
418        // tests ChordAF_Vector.mag() <= maximu    422        // tests ChordAF_Vector.mag() <= maximum_lateral_displacement 
419                                                   423 
420 #ifdef G4DEBUG_LOCATE_INTERSECTION                424 #ifdef G4DEBUG_LOCATE_INTERSECTION  
421      G4int trigger_substepno_print= warn_subst << 425      static G4int trigger_substepno_print= warn_substeps - 20;
422                                                   426 
423      if( substep_no >= trigger_substepno_print    427      if( substep_no >= trigger_substepno_print )
424      {                                            428      {
425        G4cout << "Difficulty in converging in     429        G4cout << "Difficulty in converging in "
426               << "G4SimpleLocator::EstimateInt    430               << "G4SimpleLocator::EstimateIntersectionPoint():"
427               << G4endl                           431               << G4endl
428               << "    Substep no = " << subste    432               << "    Substep no = " << substep_no << G4endl;
429        if( substep_no == trigger_substepno_pri    433        if( substep_no == trigger_substepno_print )
430        {                                          434        {
431          printStatus( CurveStartPointVelocity,    435          printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
432                       -1.0, NewSafety, 0);        436                       -1.0, NewSafety, 0);
433        }                                          437        }
434        G4cout << " State of point A: ";           438        G4cout << " State of point A: "; 
435        printStatus( CurrentA_PointVelocity, Cu    439        printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
436                     -1.0, NewSafety, substep_n    440                     -1.0, NewSafety, substep_no-1, 0);
437        G4cout << " State of point B: ";           441        G4cout << " State of point B: "; 
438        printStatus( CurrentA_PointVelocity, Cu    442        printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
439                     -1.0, NewSafety, substep_n    443                     -1.0, NewSafety, substep_no);
440      }                                            444      }
441 #endif                                            445 #endif
442      ++substep_no;                             << 446      substep_no++; 
443                                                   447 
444   } while ( ( ! found_approximate_intersection    448   } while ( ( ! found_approximate_intersection )
445            && ( ! there_is_no_intersection )      449            && ( ! there_is_no_intersection )     
446            && ( substep_no <= max_substeps) );    450            && ( substep_no <= max_substeps) ); // UNTIL found or failed
447                                                   451 
448   if( substep_no > max_no_seen )                  452   if( substep_no > max_no_seen )
449   {                                               453   {
450     max_no_seen = substep_no;                     454     max_no_seen = substep_no; 
451 #ifdef G4DEBUG_LOCATE_INTERSECTION                455 #ifdef G4DEBUG_LOCATE_INTERSECTION  
452     if( max_no_seen > warn_substeps )             456     if( max_no_seen > warn_substeps )
453     {                                             457     {
454       trigger_substepno_print = max_no_seen-20    458       trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps 
455     }                                             459     } 
456 #endif                                            460 #endif
457   }                                               461   }
458                                                   462 
459   if(  ( substep_no >= max_substeps)              463   if(  ( substep_no >= max_substeps)
460       && !there_is_no_intersection                464       && !there_is_no_intersection
461       && !found_approximate_intersection )        465       && !found_approximate_intersection )
462   {                                               466   {
463     G4cout << "ERROR - G4SimpleLocator::Estima    467     G4cout << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()" << G4endl
464            << "        Start and Endpoint of R    468            << "        Start and Endpoint of Requested Step:" << G4endl;
465     printStatus( CurveStartPointVelocity, Curv    469     printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
466                  -1.0, NewSafety, 0);             470                  -1.0, NewSafety, 0);
467     G4cout << G4endl                              471     G4cout << G4endl
468            << "        Start and end-point of     472            << "        Start and end-point of current Sub-Step:" << G4endl;
469     printStatus( CurrentA_PointVelocity, Curre    473     printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
470                  -1.0, NewSafety, substep_no-1    474                  -1.0, NewSafety, substep_no-1);
471     printStatus( CurrentA_PointVelocity, Curre    475     printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
472                  -1.0, NewSafety, substep_no);    476                  -1.0, NewSafety, substep_no);
473                                                   477 
474     std::ostringstream message;                   478     std::ostringstream message;
475     message << "Convergence is requiring too m    479     message << "Convergence is requiring too many substeps: "
476             << substep_no << G4endl               480             << substep_no << G4endl
477             << "          Abandoning effort to    481             << "          Abandoning effort to intersect." << G4endl
478             << "          Found intersection =    482             << "          Found intersection = "
479             << found_approximate_intersection     483             << found_approximate_intersection << G4endl
480             << "          Intersection exists     484             << "          Intersection exists = "
481             << !there_is_no_intersection << G4    485             << !there_is_no_intersection << G4endl;
482     message.precision(10);                        486     message.precision(10); 
483     G4double done_len = CurrentA_PointVelocity    487     G4double done_len = CurrentA_PointVelocity.GetCurveLength(); 
484     G4double full_len = CurveEndPointVelocity.    488     G4double full_len = CurveEndPointVelocity.GetCurveLength();
485     message << "          Undertaken only leng    489     message << "          Undertaken only length: " << done_len
486             << " out of " << full_len << " req    490             << " out of " << full_len << " required." << G4endl
487             << "          Remaining length = "    491             << "          Remaining length = " << full_len-done_len; 
488                                                   492 
489     G4Exception("G4SimpleLocator::EstimateInte    493     G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
490                 "GeomNav0003", FatalException,    494                 "GeomNav0003", FatalException, message);
491   }                                               495   }
492   else if( substep_no >= warn_substeps )          496   else if( substep_no >= warn_substeps )
493   {                                               497   {  
494     std::ostringstream message;                   498     std::ostringstream message;
495     message.precision(10);                        499     message.precision(10); 
496                                                   500 
497     message << "Many substeps while trying to     501     message << "Many substeps while trying to locate intersection." << G4endl
498             << "          Undertaken length: "    502             << "          Undertaken length: "  
499             << CurrentB_PointVelocity.GetCurve    503             << CurrentB_PointVelocity.GetCurveLength() 
500             << " - Needed: "  << substep_no <<    504             << " - Needed: "  << substep_no << " substeps." << G4endl
501             << "          Warning level = " <<    505             << "          Warning level = " << warn_substeps
502             << " and maximum substeps = " << m    506             << " and maximum substeps = " << max_substeps;
503     G4Exception("G4SimpleLocator::EstimateInte    507     G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
504                 "GeomNav1002", JustWarning, me    508                 "GeomNav1002", JustWarning, message);
505   }                                               509   }
506   return  !there_is_no_intersection; //  Succe    510   return  !there_is_no_intersection; //  Success or failure
507 }                                                 511 }
508                                                   512