Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4SimpleLocator Class Reference

#include <G4SimpleLocator.hh>

+ Inheritance diagram for G4SimpleLocator:

Public Member Functions

 G4SimpleLocator (G4Navigator *theNavigator)
 
 ~G4SimpleLocator ()
 
G4bool EstimateIntersectionPoint (const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin)
 
- 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 step)
 
G4bool IntersectChord (const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=0)
 
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 ()
 
G4NavigatorGetNavigatorFor ()
 
G4ChordFinderGetChordFinderFor ()
 
void SetSafetyParametersFor (G4bool UseSafety)
 
void AddAdjustementOfFoundIntersection (G4bool UseCorrection)
 
G4bool GetAdjustementOfFoundIntersection ()
 
void AdjustIntersections (G4bool UseCorrection)
 
G4bool AreIntersectionsAdjusted ()
 

Additional Inherited Members

- Protected Member Functions inherited from G4VIntersectionLocator
G4FieldTrack ReEstimateEndpoint (const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
 
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)
 
- Protected Attributes inherited from G4VIntersectionLocator
G4double kCarTolerance
 
G4int fVerboseLevel
 
G4bool fUseNormalCorrection
 
G4NavigatorfiNavigator
 
G4ChordFinderfiChordFinder
 
G4double fiEpsilonStep
 
G4double fiDeltaIntersection
 
G4bool fiUseSafety
 
G4NavigatorfHelpingNavigator
 
G4TouchableHistoryfpTouchable
 

Detailed Description

Definition at line 52 of file G4SimpleLocator.hh.

Constructor & Destructor Documentation

◆ G4SimpleLocator()

G4SimpleLocator::G4SimpleLocator ( G4Navigator theNavigator)

Definition at line 39 of file G4SimpleLocator.cc.

◆ ~G4SimpleLocator()

G4SimpleLocator::~G4SimpleLocator ( )

Definition at line 44 of file G4SimpleLocator.cc.

45{
46}

Member Function Documentation

◆ EstimateIntersectionPoint()

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

Implements G4VIntersectionLocator.

Definition at line 80 of file G4SimpleLocator.cc.

88{
89 // Find Intersection Point ( A, B, E ) of true path AB - start at E.
90
91 G4bool found_approximate_intersection = false;
92 G4bool there_is_no_intersection = false;
93
94 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
95 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
96 G4ThreeVector CurrentE_Point = TrialPoint;
97 G4bool validNormalAtE = false;
98 G4ThreeVector NormalAtEntry;
99
100 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
101 G4double NewSafety = 0.0;
102 G4bool last_AF_intersection = false;
103 G4bool final_section = true; // Shows whether current section is last
104 // (i.e. B=full end)
105 recalculatedEndPoint = false;
106
107 G4bool restoredFullEndpoint = false;
108
109 G4int substep_no = 0;
110
111 // Limits for substep number
112 //
113 const G4int max_substeps = 100000000; // Test 120 (old value 100 )
114 const G4int warn_substeps = 1000; // 100
115
116 // Statistics for substeps
117 //
118 static G4int max_no_seen= -1;
119
120 NormalAtEntry = GetSurfaceNormal( CurrentE_Point, validNormalAtE);
121
122#ifdef G4DEBUG_FIELD
123 static G4double tolerance = 1.0e-8;
124 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
125 if( (TrialPoint - StartPosition).mag() < tolerance * mm )
126 {
127 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
128 "GeomNav1002", JustWarning,
129 "Intersection point F is exactly at start point A." );
130 }
131#endif
132
133 do
134 {
135 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
136 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
137
138 // F = a point on true AB path close to point E
139 // (the closest if possible)
140 //
141 ApproxIntersecPointV = GetChordFinderFor()
142 ->ApproxCurvePointV( CurrentA_PointVelocity,
143 CurrentB_PointVelocity,
144 CurrentE_Point,
146 // The above method is the key & most intuitive part ...
147
148#ifdef G4DEBUG_FIELD
149 if( ApproxIntersecPointV.GetCurveLength() >
150 CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
151 {
152 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
153 "GeomNav0003", FatalException,
154 "Intermediate F point is past end B point!" );
155 }
156#endif
157
158 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
159
160 // First check whether EF is small - then F is a good approx. point
161 // Calculate the length and direction of the chord AF
162 //
163 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
164
165 G4ThreeVector NewMomentumDir= ApproxIntersecPointV.GetMomentumDir();
166 G4double MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ;
167
168 G4ThreeVector ChordAB = Point_B - Point_A;
169
170#ifdef G4DEBUG_FIELD
172 ReportTrialStep( substep_no, ChordAB, ChordEF_Vector,
173 NewMomentumDir, NormalAtEntry, validNormalAtE );
174#endif
175 // Check Sign is always exiting !! TODO
176 // Could ( > -epsilon) be used instead?
177 //
178 G4bool adequate_angle = ( MomDir_dot_Norm >= 0.0 )
179 || (! validNormalAtE ); // Invalid
180 G4double EF_dist2= ChordEF_Vector.mag2();
181 if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
182 || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
183 {
184 found_approximate_intersection = true;
185
186 // Create the "point" return value
187 //
188 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
189 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
190
192 {
193 // Try to Get Correction of IntersectionPoint using SurfaceNormal()
194 //
195 G4ThreeVector IP;
196 G4ThreeVector MomentumDir= ApproxIntersecPointV.GetMomentumDirection();
197 G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A,
198 CurrentE_Point, CurrentF_Point, MomentumDir,
199 last_AF_intersection, IP, NewSafety,
200 fPreviousSafety, fPreviousSftOrigin );
201
202 if(goodCorrection)
203 {
204 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
205 IntersectedOrRecalculatedFT.SetPosition(IP);
206 }
207 }
208
209 // Note: in order to return a point on the boundary,
210 // we must return E. But it is F on the curve.
211 // So we must "cheat": we are using the position at point E
212 // and the velocity at point F !!!
213 //
214 // This must limit the length we can allow for displacement!
215 }
216 else // E is NOT close enough to the curve (ie point F)
217 {
218 // Check whether any volumes are encountered by the chord AF
219 // ---------------------------------------------------------
220 // First relocate to restore any Voxel etc information
221 // in the Navigator before calling ComputeStep()
222 //
224
225 G4ThreeVector PointG; // Candidate intersection point
226 G4double stepLengthAF;
227 G4bool usedNavigatorAF = false;
228 G4bool Intersects_AF = IntersectChord( Point_A,
229 CurrentF_Point,
230 NewSafety,
231 fPreviousSafety,
232 fPreviousSftOrigin,
233 stepLengthAF,
234 PointG,
235 &usedNavigatorAF );
236 last_AF_intersection = Intersects_AF;
237 if( Intersects_AF )
238 {
239 // G is our new Candidate for the intersection point.
240 // It replaces "E" and we will repeat the test to see if
241 // it is a good enough approximate point for us.
242 // B <- F
243 // E <- G
244
245 CurrentB_PointVelocity = ApproxIntersecPointV;
246 CurrentE_Point = PointG;
247
248 // Need to recalculate the Exit Normal at the new PointG
249 // Relies on a call to Navigator::ComputeStep in IntersectChord above
250 // If the safety was adequate (for the step) this would NOT be called!
251 // But this must not occur, no intersection can be found in that case,
252 // so this branch, ie if( Intersects_AF ) would not be reached!
253 //
254 G4bool validNormalLast;
255 NormalAtEntry = GetSurfaceNormal( PointG, validNormalLast );
256 validNormalAtE = validNormalLast;
257
258 // By moving point B, must take care if current
259 // AF has no intersection to try current FB!!
260 //
261 final_section= false;
262
263#ifdef G4VERBOSE
264 if( fVerboseLevel > 3 )
265 {
266 G4cout << "G4PiF::LI> Investigating intermediate point"
267 << " at s=" << ApproxIntersecPointV.GetCurveLength()
268 << " on way to full s="
269 << CurveEndPointVelocity.GetCurveLength() << G4endl;
270 }
271#endif
272 }
273 else // not Intersects_AF
274 {
275 // In this case:
276 // There is NO intersection of AF with a volume boundary.
277 // We must continue the search in the segment FB!
278 //
280
281 G4double stepLengthFB;
282 G4ThreeVector PointH;
283 G4bool usedNavigatorFB=false;
284
285 // Check whether any volumes are encountered by the chord FB
286 // ---------------------------------------------------------
287
288 G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
289 NewSafety,fPreviousSafety,
290 fPreviousSftOrigin,
291 stepLengthFB,
292 PointH, &usedNavigatorFB );
293 if( Intersects_FB )
294 {
295 // There is an intersection of FB with a volume boundary
296 // H <- First Intersection of Chord FB
297
298 // H is our new Candidate for the intersection point.
299 // It replaces "E" and we will repeat the test to see if
300 // it is a good enough approximate point for us.
301
302 // Note that F must be in volume volA (the same as A)
303 // (otherwise AF would meet a volume boundary!)
304 // A <- F
305 // E <- H
306 //
307 CurrentA_PointVelocity = ApproxIntersecPointV;
308 CurrentE_Point = PointH;
309
310 // Need to recalculate the Exit Normal at the new PointG
311 // Relies on call to Navigator::ComputeStep in IntersectChord above
312 // If safety was adequate (for the step) this would NOT be called!
313 // But this must not occur, no intersection found in that case,
314 // so this branch, i.e. if( Intersects_AF ) would not be reached!
315 //
316 G4bool validNormalLast;
317 NormalAtEntry = GetSurfaceNormal( PointH, validNormalLast );
318 validNormalAtE = validNormalLast;
319 }
320 else // not Intersects_FB
321 {
322 // There is NO intersection of FB with a volume boundary
323
324 if( final_section )
325 {
326 // If B is the original endpoint, this means that whatever
327 // volume(s) intersected the original chord, none touch the
328 // smaller chords we have used.
329 // The value of 'IntersectedOrRecalculatedFT' returned is
330 // likely not valid
331
332 there_is_no_intersection = true; // real final_section
333 }
334 else
335 {
336 // We must restore the original endpoint
337
338 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
339 CurrentB_PointVelocity = CurveEndPointVelocity;
340 restoredFullEndpoint = true;
341 }
342 } // Endif (Intersects_FB)
343 } // Endif (Intersects_AF)
344
345 // Ensure that the new endpoints are not further apart in space
346 // than on the curve due to different errors in the integration
347 //
348 G4double linDistSq, curveDist;
349 linDistSq = ( CurrentB_PointVelocity.GetPosition()
350 - CurrentA_PointVelocity.GetPosition() ).mag2();
351 curveDist = CurrentB_PointVelocity.GetCurveLength()
352 - CurrentA_PointVelocity.GetCurveLength();
353
354 // Change this condition for very strict parameters of propagation
355 //
356 if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq )
357 {
358 // Re-integrate to obtain a new B
359 //
360 G4FieldTrack newEndPointFT =
361 ReEstimateEndpoint( CurrentA_PointVelocity,
362 CurrentB_PointVelocity,
363 linDistSq, // to avoid recalculation
364 curveDist );
365 G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
366 CurrentB_PointVelocity = newEndPointFT;
367
368 if( (final_section)) // real final section
369 {
370 recalculatedEndPoint = true;
371 IntersectedOrRecalculatedFT = newEndPointFT;
372 // So that we can return it, if it is the endpoint!
373 }
374 }
375 if( curveDist < 0.0 )
376 {
377 fVerboseLevel = 5; // Print out a maximum of information
378 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
379 -1.0, NewSafety, substep_no );
380 std::ostringstream message;
381 message << "Error in advancing propagation." << G4endl
382 << " Point A (start) is " << CurrentA_PointVelocity
383 << G4endl
384 << " Point B (end) is " << CurrentB_PointVelocity
385 << G4endl
386 << " Curve distance is " << curveDist << G4endl
387 << G4endl
388 << "The final curve point is not further along"
389 << " than the original!" << G4endl;
390
391 if( recalculatedEndPoint )
392 {
393 message << "Recalculation of EndPoint was called with fEpsStep= "
395 }
396 message.precision(20);
397 message << " Point A (Curve start) is " << CurveStartPointVelocity
398 << G4endl
399 << " Point B (Curve end) is " << CurveEndPointVelocity
400 << G4endl
401 << " Point A (Current start) is " << CurrentA_PointVelocity
402 << G4endl
403 << " Point B (Current end) is " << CurrentB_PointVelocity
404 << G4endl
405 << " Point E (Trial Point) is " << CurrentE_Point
406 << G4endl
407 << " Point F (Intersection) is " << ApproxIntersecPointV
408 << G4endl
409 << " LocateIntersection parameters are : Substep no= "
410 << substep_no;
411
412 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
413 "GeomNav0003", FatalException, message);
414 }
415
416 if(restoredFullEndpoint)
417 {
418 final_section = restoredFullEndpoint;
419 restoredFullEndpoint = false;
420 }
421 } // EndIf ( E is close enough to the curve, ie point F. )
422 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
423
424#ifdef G4DEBUG_LOCATE_INTERSECTION
425 static G4int trigger_substepno_print= warn_substeps - 20;
426
427 if( substep_no >= trigger_substepno_print )
428 {
429 G4cout << "Difficulty in converging in "
430 << "G4SimpleLocator::EstimateIntersectionPoint():"
431 << G4endl
432 << " Substep no = " << substep_no << G4endl;
433 if( substep_no == trigger_substepno_print )
434 {
435 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
436 -1.0, NewSafety, 0);
437 }
438 G4cout << " State of point A: ";
439 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
440 -1.0, NewSafety, substep_no-1, 0);
441 G4cout << " State of point B: ";
442 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
443 -1.0, NewSafety, substep_no);
444 }
445#endif
446 substep_no++;
447
448 } while ( ( ! found_approximate_intersection )
449 && ( ! there_is_no_intersection )
450 && ( substep_no <= max_substeps) ); // UNTIL found or failed
451
452 if( substep_no > max_no_seen )
453 {
454 max_no_seen = substep_no;
455#ifdef G4DEBUG_LOCATE_INTERSECTION
456 if( max_no_seen > warn_substeps )
457 {
458 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
459 }
460#endif
461 }
462
463 if( ( substep_no >= max_substeps)
464 && !there_is_no_intersection
465 && !found_approximate_intersection )
466 {
467 G4cout << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()" << G4endl
468 << " Start and Endpoint of Requested Step:" << G4endl;
469 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
470 -1.0, NewSafety, 0);
471 G4cout << G4endl
472 << " Start and end-point of current Sub-Step:" << G4endl;
473 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
474 -1.0, NewSafety, substep_no-1);
475 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
476 -1.0, NewSafety, substep_no);
477
478 std::ostringstream message;
479 message << "Convergence is requiring too many substeps: "
480 << substep_no << G4endl
481 << " Abandoning effort to intersect." << G4endl
482 << " Found intersection = "
483 << found_approximate_intersection << G4endl
484 << " Intersection exists = "
485 << !there_is_no_intersection << G4endl;
486 message.precision(10);
487 G4double done_len = CurrentA_PointVelocity.GetCurveLength();
488 G4double full_len = CurveEndPointVelocity.GetCurveLength();
489 message << " Undertaken only length: " << done_len
490 << " out of " << full_len << " required." << G4endl
491 << " Remaining length = " << full_len-done_len;
492
493 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
494 "GeomNav0003", FatalException, message);
495 }
496 else if( substep_no >= warn_substeps )
497 {
498 std::ostringstream message;
499 message.precision(10);
500
501 message << "Many substeps while trying to locate intersection." << G4endl
502 << " Undertaken length: "
503 << CurrentB_PointVelocity.GetCurveLength()
504 << " - Needed: " << substep_no << " substeps." << G4endl
505 << " Warning level = " << warn_substeps
506 << " and maximum substeps = " << max_substeps;
507 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
508 "GeomNav1002", JustWarning, message);
509 }
510 return !there_is_no_intersection; // Success or failure
511}
@ JustWarning
@ FatalException
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
bool G4bool
Definition: G4Types.hh:67
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
double mag2() const
double dot(const Hep3Vector &) const
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector &currentEPoint, G4double epsStep)
G4double GetCurveLength() const
G4ThreeVector GetPosition() const
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:548
G4Navigator * GetNavigatorFor()
G4ChordFinder * GetChordFinderFor()
G4ThreeVector GetSurfaceNormal(const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
void ReportTrialStep(G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=0)
G4double GetEpsilonStepFor()
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
G4bool GetAdjustementOfFoundIntersection()
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 printStatus(const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int step)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
T sqr(const T &x)
Definition: templates.hh:145

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