Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4BrentLocator Class Reference

#include <G4BrentLocator.hh>

+ Inheritance diagram for G4BrentLocator:

Public Member Functions

 G4BrentLocator (G4Navigator *theNavigator)
 
 ~G4BrentLocator ()
 
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 ()
 
G4Navigator * GetNavigatorFor ()
 
G4ChordFinder * GetChordFinderFor ()
 
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
 
G4Navigator * fiNavigator
 
G4ChordFinder * fiChordFinder
 
G4double fiEpsilonStep
 
G4double fiDeltaIntersection
 
G4bool fiUseSafety
 
G4Navigator * fHelpingNavigator
 
G4TouchableHistory * fpTouchable
 

Detailed Description

Definition at line 50 of file G4BrentLocator.hh.

Constructor & Destructor Documentation

◆ G4BrentLocator()

G4BrentLocator::G4BrentLocator ( G4Navigator *  theNavigator)

Definition at line 38 of file G4BrentLocator.cc.

39 : G4VIntersectionLocator(theNavigator)
40{
41 // In case of too slow progress in finding Intersection Point
42 // intermediates Points on the Track must be stored.
43 // Initialise the array of Pointers [max_depth+1] to do this
44
45 G4ThreeVector zeroV(0.0,0.0,0.0);
46 for (G4int idepth=0; idepth<max_depth+1; idepth++ )
47 {
48 ptrInterMedFT[ idepth ] = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.);
49 }
50
51 // Counters for Locator
52
53 // Counter for Maximum Number Of Trial before Intersection Found
54 //
55 maxNumberOfStepsForIntersection=0;
56
57 // Counter for Number Of Calls to ReIntegrationEndPoint Method
58 //
59 maxNumberOfCallsToReIntegration=0;
60 maxNumberOfCallsToReIntegration_depth=0;
61}
G4int
int G4int
Definition: G4Types.hh:66

◆ ~G4BrentLocator()

G4BrentLocator::~G4BrentLocator ( )

Definition at line 63 of file G4BrentLocator.cc.

64{
65 for ( G4int idepth=0; idepth<max_depth+1; idepth++)
66 {
67 delete ptrInterMedFT[idepth];
68 }
69#ifdef G4DEBUG_FIELD
70 if(fVerboseLevel>0)
71 {
72 G4cout << "G4BrentLocator::Location with Max Number of Steps="
73 << maxNumberOfStepsForIntersection<<G4endl;
74 G4cout << "G4BrentLocator::ReIntegrateEndPoint was called "
75 << maxNumberOfCallsToReIntegration
76 << " times and for depth algorithm "
77 << maxNumberOfCallsToReIntegration_depth << " times." << G4endl;
78 }
79#endif
80}
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout

Member Function Documentation

◆ EstimateIntersectionPoint()

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

Implements G4VIntersectionLocator.

Definition at line 115 of file G4BrentLocator.cc.

124{
125 // Find Intersection Point ( A, B, E ) of true path AB - start at E.
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 G4double NewSafety = 0.0;
138 G4bool last_AF_intersection = false;
139
140 // G4bool final_section= true; // Shows whether current section is last
141 // (i.e. B=full end)
142 G4bool first_section = true;
143 recalculatedEndPoint = false;
144
145 G4bool restoredFullEndpoint = false;
146
147 G4int oldprc; // cout, cerr precision
148 G4int substep_no = 0;
149
150 // Limits for substep number
151 //
152 const G4int max_substeps= 10000; // Test 120 (old value 100 )
153 const G4int warn_substeps= 1000; // 100
154
155 // Statistics for substeps
156 //
157 static G4int max_no_seen= -1;
158
159 // Counter for restarting Bintermed
160 //
161 G4int restartB = 0;
162
163 //--------------------------------------------------------------------------
164 // Algorithm for the case if progress in founding intersection is too slow.
165 // Process is defined too slow if after N=param_substeps advances on the
166 // path, it will be only 'fraction_done' of the total length.
167 // In this case the remaining length is divided in two half and
168 // the loop is restarted for each half.
169 // If progress is still too slow, the division in two halfs continue
170 // until 'max_depth'.
171 //--------------------------------------------------------------------------
172
173 const G4int param_substeps=50; // Test value for the maximum number
174 // of substeps
175 const G4double fraction_done=0.3;
176
177 G4bool Second_half = false; // First half or second half of divided step
178
179 NormalAtEntry = GetSurfaceNormal(CurrentE_Point, validNormalAtE);
180
181 // We need to know this for the 'final_section':
182 // real 'final_section' or first half 'final_section'
183 // In algorithm it is considered that the 'Second_half' is true
184 // and it becomes false only if we are in the first-half of level
185 // depthness or if we are in the first section
186
187 G4int depth=0; // Depth counts how many subdivisions of initial step made
188
189#ifdef G4DEBUG_FIELD
190 static G4double tolerance= 1.0e-8;
191 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
192 if( (TrialPoint - StartPosition).mag() < tolerance * mm )
193 {
194 G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
195 "GeomNav1002", JustWarning,
196 "Intersection point F is exactly at start point A." );
197 }
198#endif
199
200 // Intermediates Points on the Track = Subdivided Points must be stored.
201 // Give the initial values to 'InterMedFt'
202 // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
203 //
204 *ptrInterMedFT[0] = CurveEndPointVelocity;
205 for (G4int idepth=1; idepth<max_depth+1; idepth++ )
206 {
207 *ptrInterMedFT[idepth]=CurveStartPointVelocity;
208 }
209
210 //Final_section boolean store
211 G4bool fin_section_depth[max_depth];
212 for (G4int idepth=0; idepth<max_depth; idepth++ )
213 {
214 fin_section_depth[idepth]=true;
215 }
216
217 // 'SubStartPoint' is needed to calculate the length of the divided step
218 //
219 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
220
221 do
222 {
223 G4int substep_no_p = 0;
224 G4bool sub_final_section = false; // the same as final_section,
225 // but for 'sub_section'
226 SubStart_PointVelocity = CurrentA_PointVelocity;
227 do // REPEAT param
228 {
229 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
230 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
231
232 // F = a point on true AB path close to point E
233 // (the closest if possible)
234 //
235 if(substep_no_p==0)
236 {
237 ApproxIntersecPointV = GetChordFinderFor()
238 ->ApproxCurvePointV( CurrentA_PointVelocity,
239 CurrentB_PointVelocity,
240 CurrentE_Point,
241 GetEpsilonStepFor());
242 // The above method is the key & most intuitive part ...
243 }
244#ifdef G4DEBUG_FIELD
245 if( ApproxIntersecPointV.GetCurveLength() >
246 CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
247 {
248 G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
249 "GeomNav0003", FatalException,
250 "Intermediate F point is past end B point" );
251 }
252#endif
253
254 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
255
256 // First check whether EF is small - then F is a good approx. point
257 // Calculate the length and direction of the chord AF
258 //
259 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
260 G4ThreeVector NewMomentumDir= ApproxIntersecPointV.GetMomentumDir();
261 G4double MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ;
262
263#ifdef G4DEBUG_FIELD
264 G4ThreeVector ChordAB = Point_B - Point_A;
265
266 G4VIntersectionLocator::ReportTrialStep( substep_no, ChordAB,
267 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
268#endif
269
270 G4bool adequate_angle;
271 adequate_angle = ( MomDir_dot_Norm >= 0.0 ) // Can use -epsilon instead.
272 || (! validNormalAtE ); // Makes criterion invalid
273 G4double EF_dist2 = ChordEF_Vector.mag2();
274 if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
275 || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
276 {
277 found_approximate_intersection = true;
278
279 // Create the "point" return value
280 //
281 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
282 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
283
284 if ( GetAdjustementOfFoundIntersection() )
285 {
286 // Try to Get Correction of IntersectionPoint using SurfaceNormal()
287 //
288 G4ThreeVector IP;
289 G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
290 G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A,
291 CurrentE_Point, CurrentF_Point, MomentumDir,
292 last_AF_intersection, IP, NewSafety,
293 fPreviousSafety, fPreviousSftOrigin );
294 if ( goodCorrection )
295 {
296 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
297 IntersectedOrRecalculatedFT.SetPosition(IP);
298 }
299 }
300
301 // Note: in order to return a point on the boundary,
302 // we must return E. But it is F on the curve.
303 // So we must "cheat": we are using the position at point E
304 // and the velocity at point F !!!
305 //
306 // This must limit the length we can allow for displacement!
307 }
308 else // E is NOT close enough to the curve (ie point F)
309 {
310 // Check whether any volumes are encountered by the chord AF
311 // ---------------------------------------------------------
312 // First relocate to restore any Voxel etc information
313 // in the Navigator before calling ComputeStep()
314 //
315 GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A );
316
317 G4ThreeVector PointG; // Candidate intersection point
318 G4double stepLengthAF;
319 G4bool usedNavigatorAF = false;
320 G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
321 NewSafety,fPreviousSafety,
322 fPreviousSftOrigin,
323 stepLengthAF,
324 PointG,
325 &usedNavigatorAF);
326 last_AF_intersection = Intersects_AF;
327 if( Intersects_AF )
328 {
329 // G is our new Candidate for the intersection point.
330 // It replaces "E" and we will repeat the test to see if
331 // it is a good enough approximate point for us.
332 // B <- F
333 // E <- G
334 //
335 G4FieldTrack EndPoint = ApproxIntersecPointV;
336 ApproxIntersecPointV = GetChordFinderFor()->ApproxCurvePointS(
337 CurrentA_PointVelocity, CurrentB_PointVelocity,
338 EndPoint,CurrentE_Point, CurrentF_Point,PointG,
339 true, GetEpsilonStepFor() );
340 CurrentB_PointVelocity = EndPoint;
341 CurrentE_Point = PointG;
342
343 // Need to recalculate the Exit Normal at the new PointG
344 // Know that a call was made to Navigator::ComputeStep in
345 // IntersectChord above.
346 //
347 G4bool validNormalLast;
348 NormalAtEntry = GetSurfaceNormal( PointG, validNormalLast );
349 validNormalAtE = validNormalLast;
350
351 // By moving point B, must take care if current
352 // AF has no intersection to try current FB!!
353 //
354 fin_section_depth[depth] = false;
355#ifdef G4VERBOSE
356 if( fVerboseLevel > 3 )
357 {
358 G4cout << "G4PiF::LI> Investigating intermediate point"
359 << " at s=" << ApproxIntersecPointV.GetCurveLength()
360 << " on way to full s="
361 << CurveEndPointVelocity.GetCurveLength() << G4endl;
362 }
363#endif
364 }
365 else // not Intersects_AF
366 {
367 // In this case:
368 // There is NO intersection of AF with a volume boundary.
369 // We must continue the search in the segment FB!
370 //
371 GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
372
373 G4double stepLengthFB;
374 G4ThreeVector PointH;
375 G4bool usedNavigatorFB = false;
376
377 // Check whether any volumes are encountered by the chord FB
378 // ---------------------------------------------------------
379
380 G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
381 NewSafety,fPreviousSafety,
382 fPreviousSftOrigin,
383 stepLengthFB,
384 PointH,
385 &usedNavigatorFB);
386 if( Intersects_FB )
387 {
388 // There is an intersection of FB with a volume boundary
389 // H <- First Intersection of Chord FB
390
391 // H is our new Candidate for the intersection point.
392 // It replaces "E" and we will repeat the test to see if
393 // it is a good enough approximate point for us.
394
395 // Note that F must be in volume volA (the same as A)
396 // (otherwise AF would meet a volume boundary!)
397 // A <- F
398 // E <- H
399 //
400 G4FieldTrack InterMed=ApproxIntersecPointV;
401 ApproxIntersecPointV = GetChordFinderFor()->ApproxCurvePointS(
402 CurrentA_PointVelocity,CurrentB_PointVelocity,
403 InterMed,CurrentE_Point,CurrentF_Point,PointH,
404 false,GetEpsilonStepFor());
405 CurrentA_PointVelocity = InterMed;
406 CurrentE_Point = PointH;
407
408 // Need to recalculate the Exit Normal at the new PointG
409 //
410 G4bool validNormalLast;
411 NormalAtEntry = GetSurfaceNormal( PointH, validNormalLast );
412 validNormalAtE= validNormalLast;
413 }
414 else // not Intersects_FB
415 {
416 // There is NO intersection of FB with a volume boundary
417
418 if( fin_section_depth[depth] )
419 {
420 // If B is the original endpoint, this means that whatever
421 // volume(s) intersected the original chord, none touch the
422 // smaller chords we have used.
423 // The value of 'IntersectedOrRecalculatedFT' returned is
424 // likely not valid
425
426 // Check on real final_section or SubEndSection
427 //
428 if( ((Second_half)&&(depth==0)) || (first_section) )
429 {
430 there_is_no_intersection = true; // real final_section
431 }
432 else
433 {
434 // end of subsection, not real final section
435 // exit from the and go to the depth-1 level
436
437 substep_no_p = param_substeps+2; // exit from the loop
438
439 // but 'Second_half' is still true because we need to find
440 // the 'CurrentE_point' for the next loop
441 //
442 Second_half = true;
443 sub_final_section = true;
444 }
445 }
446 else
447 {
448 if(depth==0)
449 {
450 // We must restore the original endpoint
451 //
452 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
453 CurrentB_PointVelocity = CurveEndPointVelocity;
454 SubStart_PointVelocity = CurrentA_PointVelocity;
455 ApproxIntersecPointV = GetChordFinderFor()
456 ->ApproxCurvePointV( CurrentA_PointVelocity,
457 CurrentB_PointVelocity,
458 CurrentE_Point,
459 GetEpsilonStepFor());
460
461 restoredFullEndpoint = true;
462 restartB++; // counter
463 }
464 else
465 {
466 // We must restore the depth endpoint
467 //
468 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
469 CurrentB_PointVelocity = *ptrInterMedFT[depth];
470 SubStart_PointVelocity = CurrentA_PointVelocity;
471 ApproxIntersecPointV = GetChordFinderFor()
472 ->ApproxCurvePointV( CurrentA_PointVelocity,
473 CurrentB_PointVelocity,
474 CurrentE_Point,
475 GetEpsilonStepFor());
476 restoredFullEndpoint = true;
477 restartB++; // counter
478 }
479 }
480 } // Endif (Intersects_FB)
481 } // Endif (Intersects_AF)
482
483 // Ensure that the new endpoints are not further apart in space
484 // than on the curve due to different errors in the integration
485 //
486 G4double linDistSq, curveDist;
487 linDistSq = ( CurrentB_PointVelocity.GetPosition()
488 - CurrentA_PointVelocity.GetPosition() ).mag2();
489 curveDist = CurrentB_PointVelocity.GetCurveLength()
490 - CurrentA_PointVelocity.GetCurveLength();
491
492 // Change this condition for very strict parameters of propagation
493 //
494 if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq )
495 {
496 // Re-integrate to obtain a new B
497 //
498 G4FieldTrack newEndPointFT=
499 ReEstimateEndpoint( CurrentA_PointVelocity,
500 CurrentB_PointVelocity,
501 linDistSq, // to avoid recalculation
502 curveDist );
503 G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
504 CurrentB_PointVelocity = newEndPointFT;
505
506 if ( (fin_section_depth[depth]) // real final section
507 &&( first_section || ((Second_half)&&(depth==0)) ) )
508 {
509 recalculatedEndPoint = true;
510 IntersectedOrRecalculatedFT = newEndPointFT;
511 // So that we can return it, if it is the endpoint!
512 }
513 }
514 if( curveDist < 0.0 )
515 {
516 fVerboseLevel = 5; // Print out a maximum of information
517 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
518 -1.0, NewSafety, substep_no );
519 std::ostringstream message;
520 message << "Error in advancing propagation." << G4endl
521 << " Error in advancing propagation." << G4endl
522 << " Point A (start) is " << CurrentA_PointVelocity
523 << G4endl
524 << " Point B (end) is " << CurrentB_PointVelocity
525 << G4endl
526 << " Curve distance is " << curveDist << G4endl
527 << G4endl
528 << "The final curve point is not further along"
529 << " than the original!" << G4endl;
530
531 if( recalculatedEndPoint )
532 {
533 message << "Recalculation of EndPoint was called with fEpsStep= "
534 << GetEpsilonStepFor() << G4endl;
535 }
536 oldprc = G4cerr.precision(20);
537 message << " Point A (Curve start) is " << CurveStartPointVelocity
538 << G4endl
539 << " Point B (Curve end) is " << CurveEndPointVelocity
540 << G4endl
541 << " Point A (Current start) is " << CurrentA_PointVelocity
542 << G4endl
543 << " Point B (Current end) is " << CurrentB_PointVelocity
544 << G4endl
545 << " Point S (Sub start) is " << SubStart_PointVelocity
546 << G4endl
547 << " Point E (Trial Point) is " << CurrentE_Point
548 << G4endl
549 << " Old Point F(Intersection) is " << CurrentF_Point
550 << G4endl
551 << " New Point F(Intersection) is " << ApproxIntersecPointV
552 << G4endl
553 << " LocateIntersection parameters are : Substep no= "
554 << substep_no << G4endl
555 << " Substep depth no= "<< substep_no_p << " Depth= "
556 << depth << G4endl
557 << " Restarted no= "<< restartB << " Epsilon= "
558 << GetEpsilonStepFor() <<" DeltaInters= "
559 << GetDeltaIntersectionFor();
560 G4cerr.precision( oldprc );
561
562 G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
563 "GeomNav0003", FatalException, message);
564 }
565
566 if(restoredFullEndpoint)
567 {
568 fin_section_depth[depth] = restoredFullEndpoint;
569 restoredFullEndpoint = false;
570 }
571 } // EndIf ( E is close enough to the curve, ie point F. )
572 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
573
574#ifdef G4DEBUG_LOCATE_INTERSECTION
575 static G4int trigger_substepno_print= warn_substeps - 20 ;
576
577 if( substep_no >= trigger_substepno_print )
578 {
579 G4cout << "Difficulty in converging in "
580 << "G4BrentLocator::EstimateIntersectionPoint()"
581 << G4endl
582 << " Substep no = " << substep_no << G4endl;
583 if( substep_no == trigger_substepno_print )
584 {
585 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
586 -1.0, NewSafety, 0);
587 }
588 G4cout << " State of point A: ";
589 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
590 -1.0, NewSafety, substep_no-1, 0);
591 G4cout << " State of point B: ";
592 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
593 -1.0, NewSafety, substep_no);
594 }
595#endif
596 substep_no++;
597 substep_no_p++;
598
599 } while ( ( ! found_approximate_intersection )
600 && ( ! there_is_no_intersection )
601 && ( substep_no_p <= param_substeps) ); // UNTIL found or
602 // failed param substep
603 first_section = false;
604
605 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
606 {
607 G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
608 - SubStart_PointVelocity.GetCurveLength());
609 G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
610 - SubStart_PointVelocity.GetCurveLength());
611
612 G4double stepLengthAB;
613 G4ThreeVector PointGe;
614
615 // Check if progress is too slow and if it possible to go deeper,
616 // then halve the step if so
617 //
618 if ( ( did_len < fraction_done*all_len )
619 && (depth<max_depth) && (!sub_final_section) )
620 {
621 Second_half=false;
622 depth++;
623
624 G4double Sub_len = (all_len-did_len)/(2.);
625 G4FieldTrack start = CurrentA_PointVelocity;
626 G4MagInt_Driver* integrDriver =
627 GetChordFinderFor()->GetIntegrationDriver();
628 integrDriver->AccurateAdvance(start, Sub_len, GetEpsilonStepFor());
629 *ptrInterMedFT[depth] = start;
630 CurrentB_PointVelocity = *ptrInterMedFT[depth];
631
632 // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
633 //
634 SubStart_PointVelocity = CurrentA_PointVelocity;
635
636 // Find new trial intersection point needed at start of the loop
637 //
638 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
639 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
640
641 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
642 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point,
643 NewSafety, fPreviousSafety,
644 fPreviousSftOrigin,stepLengthAB,
645 PointGe);
646 if( Intersects_AB )
647 {
648 last_AF_intersection = Intersects_AB;
649 CurrentE_Point = PointGe;
650 fin_section_depth[depth]=true;
651
652 // Need to recalculate the Exit Normal at the new PointG
653 //
654 G4bool validNormalAB;
655 NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
656 validNormalAtE= validNormalAB;
657 }
658 else
659 {
660 // No intersection found for first part of curve
661 // (CurrentA,InterMedPoint[depth]). Go to the second part
662 //
663 Second_half = true;
664 }
665 } // if did_len
666
667 if( (Second_half)&&(depth!=0) )
668 {
669 // Second part of curve (InterMed[depth],Intermed[depth-1]) )
670 // On the depth-1 level normally we are on the 'second_half'
671
672 Second_half = true;
673
674 // Find new trial intersection point needed at start of the loop
675 //
676 SubStart_PointVelocity = *ptrInterMedFT[depth];
677 CurrentA_PointVelocity = *ptrInterMedFT[depth];
678 CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
679 // Ensure that the new endpoints are not further apart in space
680 // than on the curve due to different errors in the integration
681 //
682 G4double linDistSq, curveDist;
683 linDistSq = ( CurrentB_PointVelocity.GetPosition()
684 - CurrentA_PointVelocity.GetPosition() ).mag2();
685 curveDist = CurrentB_PointVelocity.GetCurveLength()
686 - CurrentA_PointVelocity.GetCurveLength();
687 if( curveDist*curveDist*(1+2*GetEpsilonStepFor() ) < linDistSq )
688 {
689 // Re-integrate to obtain a new B
690 //
691 G4FieldTrack newEndPointFT=
692 ReEstimateEndpoint( CurrentA_PointVelocity,
693 CurrentB_PointVelocity,
694 linDistSq, // to avoid recalculation
695 curveDist );
696 G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
697 CurrentB_PointVelocity = newEndPointFT;
698 if (depth==1)
699 {
700 recalculatedEndPoint = true;
701 IntersectedOrRecalculatedFT = newEndPointFT;
702 // So that we can return it, if it is the endpoint!
703 }
704 }
705
706
707 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
708 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
709 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
710 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point, NewSafety,
711 fPreviousSafety,
712 fPreviousSftOrigin,stepLengthAB, PointGe);
713 if( Intersects_AB )
714 {
715 last_AF_intersection = Intersects_AB;
716 CurrentE_Point = PointGe;
717
718 G4bool validNormalAB;
719 NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
720 validNormalAtE = validNormalAB;
721 }
722
723 depth--;
724 fin_section_depth[depth]=true;
725 }
726 } // if(!found_aproximate_intersection)
727
728 } while ( ( ! found_approximate_intersection )
729 && ( ! there_is_no_intersection )
730 && ( substep_no <= max_substeps) ); // UNTIL found or failed
731
732 if( substep_no > max_no_seen )
733 {
734 max_no_seen = substep_no;
735#ifdef G4DEBUG_LOCATE_INTERSECTION
736 if( max_no_seen > warn_substeps )
737 {
738 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
739 }
740#endif
741 }
742
743 if( ( substep_no >= max_substeps)
744 && !there_is_no_intersection
745 && !found_approximate_intersection )
746 {
747 G4cout << "ERROR - G4BrentLocator::EstimateIntersectionPoint()" << G4endl
748 << " Start and end-point of requested Step:" << G4endl;
749 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
750 -1.0, NewSafety, 0);
751 G4cout << " Start and end-point of current Sub-Step:" << G4endl;
752 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
753 -1.0, NewSafety, substep_no-1);
754 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
755 -1.0, NewSafety, substep_no);
756 std::ostringstream message;
757 message << "Too many substeps!" << G4endl
758 << " Convergence is requiring too many substeps: "
759 << substep_no << G4endl
760 << " Abandoning effort to intersect. " << G4endl
761 << " Found intersection = "
762 << found_approximate_intersection << G4endl
763 << " Intersection exists = "
764 << !there_is_no_intersection << G4endl;
765 oldprc = G4cout.precision( 10 );
766 G4double done_len = CurrentA_PointVelocity.GetCurveLength();
767 G4double full_len = CurveEndPointVelocity.GetCurveLength();
768 message << " Undertaken only length: " << done_len
769 << " out of " << full_len << " required." << G4endl
770 << " Remaining length = " << full_len - done_len;
771 G4cout.precision( oldprc );
772
773 G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
774 "GeomNav0003", FatalException, message);
775 }
776 else if( substep_no >= warn_substeps )
777 {
778 oldprc= G4cout.precision( 10 );
779 std::ostringstream message;
780 message << "Many substeps while trying to locate intersection."
781 << G4endl
782 << " Undertaken length: "
783 << CurrentB_PointVelocity.GetCurveLength()
784 << " - Needed: " << substep_no << " substeps." << G4endl
785 << " Warning level = " << warn_substeps
786 << " and maximum substeps = " << max_substeps;
787 G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
788 "GeomNav1002", JustWarning, message);
789 G4cout.precision( oldprc );
790 }
791 return !there_is_no_intersection; // Success or failure
792}
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:64
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60
G4double
double G4double
Definition: G4Types.hh:64
G4bool
bool G4bool
Definition: G4Types.hh:67
G4cerr
G4DLLIMPORT std::ostream G4cerr
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
G4ChordFinder::ApproxCurvePointV
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector &currentEPoint, G4double epsStep)
Definition: G4ChordFinder.cc:525
G4ChordFinder::ApproxCurvePointS
G4FieldTrack ApproxCurvePointS(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4FieldTrack &ApproxCurveV, const G4ThreeVector &currentEPoint, const G4ThreeVector &currentFPoint, const G4ThreeVector &PointG, G4bool first, G4double epsStep)
Definition: G4ChordFinder.cc:413
G4ChordFinder::GetIntegrationDriver
G4MagInt_Driver * GetIntegrationDriver()
G4FieldTrack::GetCurveLength
G4double GetCurveLength() const
G4FieldTrack::GetPosition
G4ThreeVector GetPosition() const
G4MagInt_Driver
Definition: G4MagIntegratorDriver.hh:49
G4MagInt_Driver::AccurateAdvance
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
Definition: G4MagIntegratorDriver.cc:127
G4Navigator::LocateGlobalPointWithinVolume
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:548
G4VIntersectionLocator::GetNavigatorFor
G4Navigator * GetNavigatorFor()
G4VIntersectionLocator::GetChordFinderFor
G4ChordFinder * GetChordFinderFor()
G4VIntersectionLocator::GetSurfaceNormal
G4ThreeVector GetSurfaceNormal(const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
Definition: G4VIntersectionLocator.cc:424
G4VIntersectionLocator::ReportTrialStep
void ReportTrialStep(G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
Definition: G4VIntersectionLocator.cc:516
G4VIntersectionLocator::IntersectChord
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=0)
G4VIntersectionLocator::GetEpsilonStepFor
G4double GetEpsilonStepFor()
G4VIntersectionLocator::ReEstimateEndpoint
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
Definition: G4VIntersectionLocator.cc:171
G4VIntersectionLocator::GetAdjustementOfFoundIntersection
G4bool GetAdjustementOfFoundIntersection()
G4VIntersectionLocator::GetDeltaIntersectionFor
G4double GetDeltaIntersectionFor()
G4VIntersectionLocator::AdjustmentOfFoundIntersection
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)
Definition: G4VIntersectionLocator.cc:348
G4VIntersectionLocator::printStatus
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int step)
Definition: G4VIntersectionLocator.cc:75
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
sqr
T sqr(const T &x)
Definition: templates.hh:145
The documentation for this class was generated from the following files: