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

#include <G4MultiLevelLocator.hh>

+ Inheritance diagram for G4MultiLevelLocator:

Public Member Functions

 G4MultiLevelLocator (G4Navigator *theNavigator)
 
 ~G4MultiLevelLocator () override
 
G4bool EstimateIntersectionPoint (const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin) override
 
void ReportStatistics ()
 
void SetMaxSteps (unsigned int valMax)
 
void SetWarnSteps (unsigned int valWarn)
 
- Public Member Functions inherited from G4VIntersectionLocator
 G4VIntersectionLocator (G4Navigator *theNavigator)
 
virtual ~G4VIntersectionLocator ()
 
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 ()
 
G4NavigatorGetNavigatorFor ()
 
G4ChordFinderGetChordFinderFor ()
 
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
 
G4NavigatorfiNavigator
 
G4ChordFinderfiChordFinder = nullptr
 
G4double fiEpsilonStep = -1.0
 
G4double fiDeltaIntersection = -1.0
 
G4NavigatorfHelpingNavigator
 
G4TouchableHistoryfpTouchable = 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.

40 : G4VIntersectionLocator(theNavigator)
41{
42 // In case of too slow progress in finding Intersection Point
43 // intermediates Points on the Track must be stored.
44 // Initialise the array of Pointers [max_depth+1] to do this
45
46 G4ThreeVector zeroV(0.0,0.0,0.0);
47 for (auto & idepth : ptrInterMedFT)
48 {
49 idepth = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.);
50 }
51
52 if (fCheckMode)
53 {
54 // Trial values Loose Medium Tight
55 // To happen: Infrequent Occasional Often
56 SetMaxSteps(150); // 300 85 25
57 SetWarnSteps(80); // 250 60 15
58 }
59}
void SetMaxSteps(unsigned int valMax)
void SetWarnSteps(unsigned int valWarn)
G4VIntersectionLocator(G4Navigator *theNavigator)

◆ ~G4MultiLevelLocator()

G4MultiLevelLocator::~G4MultiLevelLocator ( )
override

Definition at line 61 of file G4MultiLevelLocator.cc.

62{
63 for (auto & idepth : ptrInterMedFT)
64 {
65 delete idepth;
66 }
67#ifdef G4DEBUG_FIELD
69#endif
70}

Member Function Documentation

◆ EstimateIntersectionPoint()

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

Implements G4VIntersectionLocator.

Definition at line 115 of file G4MultiLevelLocator.cc.

123{
124 // Find Intersection Point ( A, B, E ) of true path AB - start at E.
125 const char* MethodName= "G4MultiLevelLocator::EstimateIntersectionPoint()";
126
127 G4bool found_approximate_intersection = false;
128 G4bool there_is_no_intersection = false;
129
130 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
131 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
132 G4ThreeVector CurrentE_Point = TrialPoint;
133 G4bool validNormalAtE = false;
134 G4ThreeVector NormalAtEntry;
135
136 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
137 G4bool validIntersectP= true; // Is it current ?
138 G4double NewSafety = 0.0;
139 G4bool last_AF_intersection = false;
140
141 auto integrDriver = GetChordFinderFor()->GetIntegrationDriver();
142 G4bool driverReIntegrates = integrDriver->DoesReIntegrate();
143
144 G4bool first_section = true;
145 recalculatedEndPoint = false;
146 G4bool restoredFullEndpoint = false;
147
148 unsigned int substep_no = 0;
149
150 // Statistics for substeps
151 static G4ThreadLocal unsigned int max_no_seen= 0;
152
153#ifdef G4DEBUG_FIELD
154 unsigned int trigger_substepno_print = 0;
155 const G4double tolerance = 1.0e-8 * CLHEP::mm;
156 unsigned int biggest_depth = 0;
157 // using kInitialisingCL = G4LocatorChangeRecord::kInitialisingCL;
158#endif
159
160 // Log the location, iteration of changes in A,B
161 //----------------------------------------------
162 static G4ThreadLocal G4LocatorChangeLogger endChangeA("StartPointA"),
163 endChangeB("EndPointB"), recApproxPoint("ApproxPoint"),
164 pointH_logger("Trial points 'E': position, normal");
165 unsigned int eventCount = 0;
166
167 if (fCheckMode)
168 {
169 // Clear previous call's data
170 endChangeA.clear();
171 endChangeB.clear();
172 recApproxPoint.clear();
173 pointH_logger.clear();
174
175 // Record the initialisation
176 ++eventCount;
177 endChangeA.AddRecord( G4LocatorChangeRecord::kInitialisingCL, substep_no,
178 eventCount, CurrentA_PointVelocity );
179 endChangeB.AddRecord( G4LocatorChangeRecord::kInitialisingCL, substep_no,
180 eventCount, CurrentB_PointVelocity );
181 }
182
183 //--------------------------------------------------------------------------
184 // Algorithm for the case if progress in founding intersection is too slow.
185 // Process is defined too slow if after N=param_substeps advances on the
186 // path, it will be only 'fraction_done' of the total length.
187 // In this case the remaining length is divided in two half and
188 // the loop is restarted for each half.
189 // If progress is still too slow, the division in two halfs continue
190 // until 'max_depth'.
191 //--------------------------------------------------------------------------
192
193 const G4int param_substeps = 5; // Test value for the maximum number
194 // of substeps
195 const G4double fraction_done = 0.3;
196
197 G4bool Second_half = false; // First half or second half of divided step
198
199 // We need to know this for the 'final_section':
200 // real 'final_section' or first half 'final_section'
201 // In algorithm it is considered that the 'Second_half' is true
202 // and it becomes false only if we are in the first-half of level
203 // depthness or if we are in the first section
204
205 unsigned int depth = 0; // Depth counts subdivisions of initial step made
206 ++fNumCalls;
207
208 NormalAtEntry = GetSurfaceNormal(CurrentE_Point, validNormalAtE);
209
210 if (fCheckMode)
211 {
212 pointH_logger.AddRecord( G4LocatorChangeRecord::kInitialisingCL,
213 substep_no, eventCount,
214 G4FieldTrack( CurrentE_Point,0.,NormalAtEntry,0.,
215 0., 1.,G4ThreeVector(0.),0.,0.) );
216 #if (G4DEBUG_FIELD>1)
217 G4ThreeVector StartPosition = CurveStartPointVelocity.GetPosition();
218 if( (TrialPoint - StartPosition).mag2() < sqr(tolerance))
219 {
220 ReportImmediateHit( MethodName, StartPosition, TrialPoint,
221 tolerance, fNumCalls);
222 }
223 #endif
224 }
225
226 // Intermediates Points on the Track = Subdivided Points must be stored.
227 // Give the initial values to 'InterMedFt'
228 // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
229 //
230 *ptrInterMedFT[0] = CurveEndPointVelocity;
231 for ( auto idepth=1; idepth<max_depth+1; ++idepth )
232 {
233 *ptrInterMedFT[idepth] = CurveStartPointVelocity;
234 }
235
236 // Final_section boolean store
237 //
238 G4bool fin_section_depth[max_depth];
239 for (bool & idepth : fin_section_depth)
240 {
241 idepth = true;
242 }
243 // 'SubStartPoint' is needed to calculate the length of the divided step
244 //
245 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
246
247 do // Loop checking, 07.10.2016, J.Apostolakis
248 {
249 unsigned int substep_no_p = 0;
250 G4bool sub_final_section = false; // the same as final_section,
251 // but for 'sub_section'
252 SubStart_PointVelocity = CurrentA_PointVelocity;
253
254 do // Loop checking, 07.10.2016, J.Apostolakis
255 { // REPEAT param
256 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
257 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
258
259#ifdef G4DEBUG_FIELD
260 const G4double lenA = CurrentA_PointVelocity.GetCurveLength() ;
261 const G4double lenB = CurrentB_PointVelocity.GetCurveLength() ;
262 G4double curv_lenAB = lenB - lenA;
263 G4double distAB = (Point_B - Point_A).mag();
264 if( curv_lenAB < distAB * ( 1. - 10.*fiEpsilonStep ) )
265 {
266 G4cerr << "ERROR> (Start) Point A coincides with or has gone past (end) point B"
267 << "MLL: iters = " << substep_no << G4endl;
268 G4long op=G4cerr.precision(6);
269 G4cerr << " Difference = " << distAB - curv_lenAB
270 << " exceeds limit of relative dist (10*epsilon)= " << 10*fiEpsilonStep
271 << " i.e. limit = " << 10 * fiEpsilonStep * distAB << G4endl;
272 G4cerr.precision(9);
273 G4cerr << " Len A, B = " << lenA << " " << lenB << G4endl
274 << " Position A: " << Point_A << G4endl
275 << " Position B: " << Point_B << G4endl;
276 G4cerr.precision(op);
277 // G4LocatorChangeRecord::ReportVector(G4cerr, "endPointB", endChangeB );
278 // G4cerr<<"EndPoints A(start) and B(end): combined changes " << G4endl;
279 if (fCheckMode)
280 {
281 G4LocatorChangeLogger::ReportEndChanges(G4cerr, endChangeA, endChangeB);
282 }
283 }
284#endif
285 if( !validIntersectP )
286 {
288 errmsg << "Assertion FAILURE - invalid (stale) Interection point. Substep: "
289 << substep_no << " call: " << fNumCalls << G4endl;
290 if (fCheckMode)
291 {
292 G4LocatorChangeRecord::ReportEndChanges(errmsg, endChangeA, endChangeB );
293 }
294 G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint", "GeomNav0004",
295 JustWarning, errmsg);
296 }
297
298 // F = a point on true AB path close to point E
299 // (the closest if possible)
300 //
301 ApproxIntersecPointV = GetChordFinderFor()
302 ->ApproxCurvePointV( CurrentA_PointVelocity,
303 CurrentB_PointVelocity,
304 CurrentE_Point,
306 // The above method is the key & most intuitive part ...
307
308#ifdef G4DEBUG_FIELD
310 substep_no, eventCount, ApproxIntersecPointV ) );
311 G4double lenIntsc= ApproxIntersecPointV.GetCurveLength();
312 G4double checkVsEnd= lenB - lenIntsc;
313
314 if( lenIntsc > lenB )
315 {
316 std::ostringstream errmsg;
317 errmsg.precision(17);
318 G4double ratio = checkVsEnd / lenB;
319 G4double ratioTol = std::fabs(ratio) / tolerance;
320 errmsg << "Intermediate F point is past end B point" << G4endl
321 << " l( intersection ) = " << lenIntsc << G4endl
322 << " l( endpoint ) = " << lenB << G4endl;
323 errmsg.precision(8);
324 errmsg << " l_end - l_inters = " << checkVsEnd << G4endl
325 << " / l_end = " << ratio << G4endl
326 << " ratio / tolerance = " << ratioTol << G4endl;
327 if( ratioTol < 1.0 )
328 G4Exception(MethodName, "GeomNav0003", JustWarning, errmsg );
329 else
330 G4Exception(MethodName, "GeomNav0003", FatalException, errmsg );
331 }
332#endif
333
334 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
335
336 // First check whether EF is small - then F is a good approx. point
337 // Calculate the length and direction of the chord AF
338 //
339 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
340
341 G4ThreeVector NewMomentumDir = ApproxIntersecPointV.GetMomentumDir();
342 G4double MomDir_dot_Norm = NewMomentumDir.dot( NormalAtEntry );
343
344#ifdef G4DEBUG_FIELD
345 if( fVerboseLevel > 3 )
346 {
347 G4ThreeVector ChordAB = Point_B - Point_A;
348 G4double ABchord_length = ChordAB.mag();
349 G4double MomDir_dot_ABchord;
350 MomDir_dot_ABchord = (1.0 / ABchord_length)
351 * NewMomentumDir.dot( ChordAB );
352 G4VIntersectionLocator::ReportTrialStep( substep_no, ChordAB,
353 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
354 G4cout << " dot( MomentumDir, ABchord_unit ) = "
355 << MomDir_dot_ABchord << G4endl;
356 }
357#endif
358 G4bool adequate_angle =
359 ( MomDir_dot_Norm >= 0.0 ) // Can use ( > -epsilon) instead
360 || (! validNormalAtE ); // Invalid, cannot use this criterion
361 G4double EF_dist2 = ChordEF_Vector.mag2();
362 if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
363 || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
364 {
365 found_approximate_intersection = true;
366
367 // Create the "point" return value
368 //
369 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
370 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
371
373 {
374 // Try to Get Correction of IntersectionPoint using SurfaceNormal()
375 //
376 G4ThreeVector IP;
377 G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
378 G4bool goodCorrection = AdjustmentOfFoundIntersection(Point_A,
379 CurrentE_Point, CurrentF_Point, MomentumDir,
380 last_AF_intersection, IP, NewSafety,
381 previousSafety, previousSftOrigin );
382 if ( goodCorrection )
383 {
384 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
385 IntersectedOrRecalculatedFT.SetPosition(IP);
386 }
387 }
388 // Note: in order to return a point on the boundary,
389 // we must return E. But it is F on the curve.
390 // So we must "cheat": we are using the position at point E
391 // and the velocity at point F !!!
392 //
393 // This must limit the length we can allow for displacement!
394 }
395 else // E is NOT close enough to the curve (ie point F)
396 {
397 // Check whether any volumes are encountered by the chord AF
398 // ---------------------------------------------------------
399 // First relocate to restore any Voxel etc information
400 // in the Navigator before calling ComputeStep()
401 //
403
404 G4ThreeVector PointG; // Candidate intersection point
405 G4double stepLengthAF;
406 G4bool Intersects_FB = false;
407 G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
408 NewSafety, previousSafety,
409 previousSftOrigin,
410 stepLengthAF,
411 PointG );
412 last_AF_intersection = Intersects_AF;
413 if( Intersects_AF )
414 {
415 // G is our new Candidate for the intersection point.
416 // It replaces "E" and we will see if it's good enough.
417 CurrentB_PointVelocity = ApproxIntersecPointV; // B <- F
418 CurrentE_Point = PointG; // E <- G
419
420 validIntersectP = true; // 'E' has been updated.
421
422 G4bool validNormalLast;
423 NormalAtEntry = GetSurfaceNormal( PointG, validNormalLast );
424 validNormalAtE = validNormalLast;
425
426 // As we move point B, must take care in case the current
427 // AF has no intersection to try current FB!!
428 fin_section_depth[depth] = false;
429
430 if (fCheckMode)
431 {
432 ++eventCount;
433 endChangeB.push_back(
435 substep_no, eventCount, CurrentB_PointVelocity) );
436 }
437#ifdef G4VERBOSE
438 if( fVerboseLevel > 3 )
439 {
440 G4cout << "G4PiF::LI> Investigating intermediate point"
441 << " at s=" << ApproxIntersecPointV.GetCurveLength()
442 << " on way to full s="
443 << CurveEndPointVelocity.GetCurveLength() << G4endl;
444 }
445#endif
446 }
447 else // not Intersects_AF
448 {
449 // In this case:
450 // There is NO intersection of AF with a volume boundary.
451 // We must continue the search in the segment FB!
452 //
454
455 G4double stepLengthFB;
456 G4ThreeVector PointH;
457
458 // Check whether any volumes are encountered by the chord FB
459 // ---------------------------------------------------------
460
461 Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
462 NewSafety, previousSafety,
463 previousSftOrigin,
464 stepLengthFB,
465 PointH );
466 if( Intersects_FB )
467 {
468 // There is an intersection of FB with a volume boundary
469 // H <- First Intersection of Chord FB
470
471 // H is our new Candidate for the intersection point.
472 // It replaces "E" and we will repeat the test to see if
473 // it is a good enough approximate point for us.
474
475 // Note that F must be in volume volA (the same as A)
476 // (otherwise AF would meet a volume boundary!)
477 // A <- F
478 // E <- H
479 //
480 CurrentA_PointVelocity = ApproxIntersecPointV;
481 CurrentE_Point = PointH;
482
483 validIntersectP = true; // 'E' has been updated.
484
485 G4bool validNormalH;
486 NormalAtEntry = GetSurfaceNormal( PointH, validNormalH );
487 validNormalAtE = validNormalH;
488
489 if (fCheckMode)
490 {
491 ++eventCount;
492 endChangeA.push_back(
494 substep_no, eventCount, CurrentA_PointVelocity) );
495 G4FieldTrack intersectH_pn('0'); // Point and normal
496 // nothing else will be valid
497 intersectH_pn.SetPosition( PointH );
498 intersectH_pn.SetMomentum( NormalAtEntry );
499 pointH_logger.AddRecord(G4LocatorChangeRecord::kIntersectsFB,
500 substep_no, eventCount, intersectH_pn );
501 }
502 }
503 else // not Intersects_FB
504 {
505 validIntersectP = false; // Intersections are now stale
506 if( fin_section_depth[depth] )
507 {
508 // If B is the original endpoint, this means that whatever
509 // volume(s) intersected the original chord, none touch the
510 // smaller chords we have used.
511 // The value of 'IntersectedOrRecalculatedFT' returned is
512 // likely not valid
513
514 // Check on real final_section or SubEndSection
515 //
516 if( ((Second_half)&&(depth==0)) || (first_section) )
517 {
518 there_is_no_intersection = true; // real final_section
519 }
520 else
521 {
522 // end of subsection, not real final section
523 // exit from the and go to the depth-1 level
524 substep_no_p = param_substeps+2; // exit from the loop
525
526 // but 'Second_half' is still true because we need to find
527 // the 'CurrentE_point' for the next loop
528 Second_half = true;
529 sub_final_section = true;
530 }
531 }
532 else
533 {
534 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
535 CurrentB_PointVelocity = (depth==0) ? CurveEndPointVelocity
536 : *ptrInterMedFT[depth] ;
537 SubStart_PointVelocity = CurrentA_PointVelocity;
538 restoredFullEndpoint = true;
539
540 validIntersectP = false; // 'E' has NOT been updated.
541
542 if (fCheckMode)
543 {
544 ++eventCount;
545 endChangeA.push_back(
548 substep_no, eventCount, CurrentA_PointVelocity) );
549 endChangeB.push_back(
552 substep_no, eventCount, CurrentB_PointVelocity) );
553 }
554 }
555 } // Endif (Intersects_FB)
556 } // Endif (Intersects_AF)
557
558 G4int errorEndPt = 0; // Default: no error (if not calling CheckAnd...
559
560 G4bool recalculatedB= false;
561 if( driverReIntegrates )
562 {
563 G4FieldTrack RevisedB_FT = CurrentB_PointVelocity;
564 recalculatedB= CheckAndReEstimateEndpoint(CurrentA_PointVelocity,
565 CurrentB_PointVelocity,
566 RevisedB_FT,
567 errorEndPt );
568 if( recalculatedB )
569 {
570 CurrentB_PointVelocity = RevisedB_FT; // Use it !
571 // Do not invalidate intersection F -- it is still roughly OK.
572 //
573 // The best course would be to invalidate (reset validIntersectP)
574 // BUT if we invalidate it, we must re-estimate it somewhere! E.g.
575 // validIntersectP = false; // 'E' has NOT been updated.
576
577 if ( (fin_section_depth[depth]) // real final section
578 &&( first_section || ((Second_half)&&(depth==0)) ) )
579 {
580 recalculatedEndPoint = true;
581 IntersectedOrRecalculatedFT = RevisedB_FT;
582 // So that we can return it, if it is the endpoint!
583 }
584 // else
585 // Move forward the other points
586 // - or better flag it, so that they are re-computed when next used
587 // [ Implementation: a counter for # of recomputations
588 // => avoids extra work]
589 }
590 if (fCheckMode)
591 {
592 ++eventCount;
593 endChangeB.push_back(
595 substep_no, eventCount, RevisedB_FT ) );
596 }
597 }
598 else
599 {
600 if( CurrentB_PointVelocity.GetCurveLength()
601 < CurrentA_PointVelocity.GetCurveLength() )
602 {
603 errorEndPt = 2;
604 }
605 }
606
607 if( errorEndPt > 1 ) // errorEndPt = 1 is milder, just: len(B)=len(A)
608 {
609 std::ostringstream errmsg;
610
612 CurveStartPointVelocity, CurveEndPointVelocity,
613 NewSafety, fiEpsilonStep,
614 CurrentA_PointVelocity, CurrentB_PointVelocity,
615 SubStart_PointVelocity, CurrentE_Point,
616 ApproxIntersecPointV, substep_no, substep_no_p, depth);
617
618 if (fCheckMode)
619 {
620 G4LocatorChangeRecord::ReportEndChanges(errmsg, endChangeA, endChangeB );
621 }
622
623 errmsg << G4endl << " * Location: " << MethodName
624 << "- After EndIf(Intersects_AF)" << G4endl;
625 errmsg << " * Bool flags: Recalculated = " << recalculatedB
626 << " Intersects_AF = " << Intersects_AF
627 << " Intersects_FB = " << Intersects_FB << G4endl;
628 errmsg << " * Number of calls to MLL:EIP= " << fNumCalls << G4endl;
629 G4Exception(MethodName, "GeomNav0003", FatalException, errmsg);
630 }
631 if( restoredFullEndpoint )
632 {
633 fin_section_depth[depth] = restoredFullEndpoint;
634 restoredFullEndpoint = false;
635 }
636 } // EndIf ( E is close enough to the curve, ie point F. )
637 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
638
639#ifdef G4DEBUG_FIELD
640 if( trigger_substepno_print == 0)
641 {
642 trigger_substepno_print= fWarnSteps - 20;
643 }
644
645 if( substep_no >= trigger_substepno_print )
646 {
647 G4cout << "Difficulty in converging in " << MethodName
648 << " Substep no = " << substep_no << G4endl;
649 if( substep_no == trigger_substepno_print )
650 {
651 G4cout << " Start: ";
652 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
653 -1.0, NewSafety, 0 );
654 if( fCheckMode ) {
655 G4LocatorChangeRecord::ReportEndChanges(G4cout, endChangeA, endChangeB );
656 } else {
657 G4cout << " ** For more information enable 'check mode' in G4MultiLevelLocator "
658 << "-- (it saves and can output change records) " << G4endl;
659 }
660 }
661 G4cout << " Point A: ";
662 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
663 -1.0, NewSafety, substep_no-1 );
664 G4cout << " Point B: ";
665 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
666 -1.0, NewSafety, substep_no );
667 }
668#endif
669 ++substep_no;
670 ++substep_no_p;
671
672 } while ( ( ! found_approximate_intersection )
673 && ( ! there_is_no_intersection )
674 && validIntersectP // New condition: must refresh intersection !!
675 && ( substep_no_p <= param_substeps) ); // UNTIL found or
676 // failed param substep
677
678 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
679 {
680 G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
681 - SubStart_PointVelocity.GetCurveLength());
682 G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
683 - SubStart_PointVelocity.GetCurveLength());
684
685 G4double distAB = -1;
686 //
687 // Is progress is too slow, and is it possible to go deeper?
688 // If so, then *halve the step*
689 // ==============
690 if( (did_len < fraction_done*all_len)
691 && (depth<max_depth) && (!sub_final_section) )
692 {
693#ifdef G4DEBUG_FIELD
694 static G4ThreadLocal unsigned int numSplits=0; // For debugging only
695 biggest_depth = std::max(depth, biggest_depth);
696 ++numSplits;
697#endif
698 Second_half = false;
699 ++depth;
700 first_section = false;
701
702 G4double Sub_len = (all_len-did_len)/(2.);
703 G4FieldTrack midPoint = CurrentA_PointVelocity;
704 G4bool fullAdvance=
705 integrDriver->AccurateAdvance(midPoint, Sub_len, fiEpsilonStep);
706
707 ++fNumAdvanceTrials;
708 if( fullAdvance ) { ++fNumAdvanceFull; }
709
710 G4double lenAchieved=
711 midPoint.GetCurveLength()-CurrentA_PointVelocity.GetCurveLength();
712
713 const G4double adequateFraction = (1.0-CLHEP::perThousand);
714 G4bool goodAdvance = (lenAchieved >= adequateFraction * Sub_len);
715 if ( goodAdvance ) { ++fNumAdvanceGood; }
716
717#ifdef G4DEBUG_FIELD
718 else // !goodAdvance
719 {
720 G4cout << "MLL> AdvanceChordLimited not full at depth=" << depth
721 << " total length achieved = " << lenAchieved << " of "
722 << Sub_len << " fraction= ";
723 if (Sub_len != 0.0 ) { G4cout << lenAchieved / Sub_len; }
724 else { G4cout << "DivByZero"; }
725 G4cout << " Good-enough fraction = " << adequateFraction;
726 G4cout << " Number of call to mll = " << fNumCalls
727 << " iteration " << substep_no
728 << " inner = " << substep_no_p << G4endl;
729 G4cout << " Epsilon of integration = " << fiEpsilonStep << G4endl;
730 G4cout << " State at start is = " << CurrentA_PointVelocity
731 << G4endl
732 << " at end (midpoint)= " << midPoint << G4endl;
733 G4cout << " Particle mass = " << midPoint.GetRestMass() << G4endl;
734
735 G4EquationOfMotion *equation = integrDriver->GetEquationOfMotion();
736 ReportFieldValue( CurrentA_PointVelocity, "start", equation );
737 ReportFieldValue( midPoint, "midPoint", equation );
738 G4cout << " Original Start = "
739 << CurveStartPointVelocity << G4endl;
740 G4cout << " Original End = "
741 << CurveEndPointVelocity << G4endl;
742 G4cout << " Original TrialPoint = "
743 << TrialPoint << G4endl;
744 G4cout << " (this is STRICT mode) "
745 << " num Splits= " << numSplits;
746 G4cout << G4endl;
747 }
748#endif
749
750 *ptrInterMedFT[depth] = midPoint;
751 CurrentB_PointVelocity = midPoint;
752
753 if (fCheckMode)
754 {
755 ++eventCount;
756 endChangeB.push_back(
758 substep_no, eventCount, midPoint) );
759 }
760
761 // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
762 //
763 SubStart_PointVelocity = CurrentA_PointVelocity;
764
765 // Find new trial intersection point needed at start of the loop
766 //
767 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
768 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
769
770 G4ThreeVector PointGe;
772 G4bool Intersects_AB = IntersectChord(Point_A, Point_B,
773 NewSafety, previousSafety,
774 previousSftOrigin, distAB,
775 PointGe);
776 if( Intersects_AB )
777 {
778 last_AF_intersection = Intersects_AB;
779 CurrentE_Point = PointGe;
780 fin_section_depth[depth] = true;
781
782 validIntersectP = true; // 'E' has been updated.
783
784 G4bool validNormalAB;
785 NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
786 validNormalAtE = validNormalAB;
787 }
788 else
789 {
790 // No intersection found for first part of curve
791 // (CurrentA,InterMedPoint[depth]). Go to the second part
792 //
793 Second_half = true;
794
795 validIntersectP= false; // No new 'E' chord intersection found
796 }
797 } // if did_len
798
799 G4bool unfinished = Second_half;
800 while ( unfinished && (depth>0) ) // Loop checking, 07.10.2016, JA
801 {
802 // Second part of curve (InterMed[depth],Intermed[depth-1]))
803 // On the depth-1 level normally we are on the 'second_half'
804
805 // Find new trial intersection point needed at start of the loop
806 //
807 SubStart_PointVelocity = *ptrInterMedFT[depth];
808 CurrentA_PointVelocity = *ptrInterMedFT[depth];
809 CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
810
811 if (fCheckMode)
812 {
813 ++eventCount;
815 substep_no, eventCount, CurrentA_PointVelocity);
816 endChangeA.push_back( chngPop_a );
818 substep_no, eventCount, CurrentB_PointVelocity);
819 endChangeB.push_back( chngPop_b );
820 }
821
822 // Ensure that the new endpoints are not further apart in space
823 // than on the curve due to different errors in the integration
824 //
825 G4int errorEndPt = -1;
826 G4bool recalculatedB= false;
827 if( driverReIntegrates )
828 {
829 // Ensure that the new endpoints are not further apart in space
830 // than on the curve due to different errors in the integration
831 //
832 G4FieldTrack RevisedEndPointFT = CurrentB_PointVelocity;
833 recalculatedB =
834 CheckAndReEstimateEndpoint( CurrentA_PointVelocity,
835 CurrentB_PointVelocity,
836 RevisedEndPointFT,
837 errorEndPt );
838 if( recalculatedB )
839 {
840 CurrentB_PointVelocity = RevisedEndPointFT; // Use it !
841
842 if ( depth == 1 )
843 {
844 recalculatedEndPoint = true;
845 IntersectedOrRecalculatedFT = RevisedEndPointFT;
846 // So that we can return it, if it is the endpoint!
847 }
848 }
849 else
850 {
851 if( CurrentB_PointVelocity.GetCurveLength()
852 < CurrentA_PointVelocity.GetCurveLength() )
853 {
854 errorEndPt = 2;
855 }
856 }
857
858 if (fCheckMode)
859 {
860 ++eventCount;
861 endChangeB.push_back(
863 substep_no, eventCount, RevisedEndPointFT));
864 }
865 }
866 if( errorEndPt > 1 ) // errorEndPt = 1 is milder, just: len(B)=len(A)
867 {
868 std::ostringstream errmsg;
869 ReportReversedPoints(errmsg,
870 CurveStartPointVelocity, CurveEndPointVelocity,
871 NewSafety, fiEpsilonStep,
872 CurrentA_PointVelocity, CurrentB_PointVelocity,
873 SubStart_PointVelocity, CurrentE_Point,
874 ApproxIntersecPointV, substep_no, substep_no_p, depth);
875 errmsg << " * Location: " << MethodName << "- Second-Half" << G4endl;
876 errmsg << " * Recalculated = " << recalculatedB << G4endl; // false
877 G4Exception(MethodName, "GeomNav0003", FatalException, errmsg);
878 }
879
880 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
881 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
882 G4ThreeVector PointGi;
884 G4bool Intersects_AB = IntersectChord(Point_A, Point_B, NewSafety,
885 previousSafety,
886 previousSftOrigin, distAB,
887 PointGi);
888 if( Intersects_AB )
889 {
890 last_AF_intersection = Intersects_AB;
891 CurrentE_Point = PointGi;
892
893 validIntersectP = true; // 'E' has been updated.
894 NormalAtEntry = GetSurfaceNormal( PointGi, validNormalAtE );
895 }
896 else
897 {
898 validIntersectP = false; // No new 'E' chord intersection found
899 if( depth == 1)
900 {
901 there_is_no_intersection = true;
902 }
903 }
904 depth--;
905 fin_section_depth[depth] = true;
906 unfinished = !validIntersectP;
907 }
908#ifdef G4DEBUG_FIELD
909 if( ! ( validIntersectP || there_is_no_intersection ) )
910 {
911 // What happened ??
912 G4cout << "MLL - WARNING Potential FAILURE: Conditions not met!"
913 << G4endl
914 << " Depth = " << depth << G4endl
915 << " Num Substeps= " << substep_no << G4endl;
916 G4cout << " Found intersection= " << found_approximate_intersection
917 << G4endl;
918 G4cout << " Progress report: -- " << G4endl;
920 CurveStartPointVelocity, CurveEndPointVelocity,
921 substep_no, CurrentA_PointVelocity,
922 CurrentB_PointVelocity,
923 NewSafety, depth );
924 G4cout << G4endl;
925 }
926#endif
927 } // if(!found_aproximate_intersection)
928
929 assert( validIntersectP || there_is_no_intersection
930 || found_approximate_intersection);
931
932 } while ( ( ! found_approximate_intersection )
933 && ( ! there_is_no_intersection )
934 && ( substep_no <= fMaxSteps) ); // UNTIL found or failed
935
936 if( substep_no > max_no_seen )
937 {
938 max_no_seen = substep_no;
939#ifdef G4DEBUG_FIELD
940 if( max_no_seen > fWarnSteps )
941 {
942 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
943 }
944#endif
945 }
946
947 if( !there_is_no_intersection && !found_approximate_intersection )
948 {
949 if( substep_no >= fMaxSteps)
950 {
951 // Since we cannot go further (yet), we return as far as we have gone
952
953 recalculatedEndPoint = true;
954 IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
955 found_approximate_intersection = false;
956
957 std::ostringstream message;
958 message << G4endl;
959 message << "Convergence is requiring too many substeps: "
960 << substep_no << " ( limit = "<< fMaxSteps << ")"
961 << G4endl
962 << " Abandoning effort to intersect. " << G4endl << G4endl;
963 message << " Number of calls to MLL: " << fNumCalls;
964 message << " iteration = " << substep_no <<G4endl << G4endl;
965
966 message.precision( 12 );
967 G4double done_len = CurrentA_PointVelocity.GetCurveLength();
968 G4double full_len = CurveEndPointVelocity.GetCurveLength();
969 message << " Undertaken only length: " << done_len
970 << " out of " << full_len << " required." << G4endl
971 << " Remaining length = " << full_len - done_len;
972
973 message << " Start and end-point of requested Step:" << G4endl;
974 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
975 -1.0, NewSafety, 0, message, -1 );
976 message << " Start and end-point of current Sub-Step:" << G4endl;
977 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
978 -1.0, NewSafety, substep_no-1, message, -1 );
979 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
980 -1.0, NewSafety, substep_no, message, -1 );
981
982 G4Exception(MethodName, "GeomNav0003", JustWarning, message);
983 }
984 else if( substep_no >= fWarnSteps)
985 {
986 std::ostringstream message;
987 message << "Many substeps while trying to locate intersection."
988 << G4endl
989 << " Undertaken length: "
990 << CurrentB_PointVelocity.GetCurveLength()
991 << " - Needed: " << substep_no << " substeps." << G4endl
992 << " Warning number = " << fWarnSteps
993 << " and maximum substeps = " << fMaxSteps;
994 G4Exception(MethodName, "GeomNav1002", JustWarning, message);
995 }
996 }
997
998 return (!there_is_no_intersection) && found_approximate_intersection;
999 // Success or failure
1000}
@ JustWarning
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
const size_t max_depth
Definition G4Octree.hh:54
CLHEP::Hep3Vector G4ThreeVector
double G4double
Definition G4Types.hh:83
long G4long
Definition G4Types.hh:87
bool G4bool
Definition G4Types.hh:86
int G4int
Definition G4Types.hh:85
G4GLOB_DLL std::ostream G4cerr
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout
double mag2() const
double dot(const Hep3Vector &) const
double mag() const
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector &currentEPoint, G4double epsStep)
G4VIntegrationDriver * GetIntegrationDriver()
G4double GetCurveLength() const
G4ThreeVector GetPosition() const
G4double GetRestMass() const
static std::ostream & ReportEndChanges(std::ostream &os, const G4LocatorChangeLogger &startA, const G4LocatorChangeLogger &endB)
static std::ostream & ReportEndChanges(std::ostream &os, const std::vector< G4LocatorChangeRecord > &startA, const std::vector< G4LocatorChangeRecord > &endB)
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
virtual G4bool DoesReIntegrate() const =0
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=nullptr)
G4bool CheckAndReEstimateEndpoint(const G4FieldTrack &CurrentStartA, const G4FieldTrack &EstimatedEndB, G4FieldTrack &RevisedEndPoint, G4int &errorCode)
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)
G4double GetEpsilonStepFor()
void ReportImmediateHit(const char *MethodName, const G4ThreeVector &StartPosition, const G4ThreeVector &TrialPoint, G4double tolerance, unsigned long int numCalls)
G4bool GetAdjustementOfFoundIntersection()
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int stepNum)
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 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)
T sqr(const T &x)
Definition templates.hh:128
#define G4ThreadLocal
Definition tls.hh:77

◆ ReportStatistics()

void G4MultiLevelLocator::ReportStatistics ( )

Definition at line 1002 of file G4MultiLevelLocator.cc.

1003{
1004 G4cout << " Number of calls = " << fNumCalls << G4endl;
1005 G4cout << " Number of split level ('advances'): "
1006 << fNumAdvanceTrials << G4endl;
1007 G4cout << " Number of full advances: "
1008 << fNumAdvanceGood << G4endl;
1009 G4cout << " Number of good advances: "
1010 << fNumAdvanceFull << G4endl;
1011}

Referenced by ~G4MultiLevelLocator().

◆ SetMaxSteps()

void G4MultiLevelLocator::SetMaxSteps ( unsigned int valMax)
inline

Definition at line 71 of file G4MultiLevelLocator.hh.

71{ fMaxSteps= valMax; }

Referenced by G4MultiLevelLocator().

◆ SetWarnSteps()

void G4MultiLevelLocator::SetWarnSteps ( unsigned int valWarn)
inline

Definition at line 72 of file G4MultiLevelLocator.hh.

72{ fWarnSteps= valWarn; }

Referenced by G4MultiLevelLocator().


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