84{
85
86
87 G4bool found_approximate_intersection =
false;
88 G4bool there_is_no_intersection =
false;
89
90 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
91 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
93 G4bool validNormalAtE =
false;
95
96 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity);
98 G4bool last_AF_intersection =
false;
99 G4bool final_section =
true;
100
101 recalculatedEndPoint = false;
102
103 G4bool restoredFullEndpoint =
false;
104
105 G4int substep_no = 0;
106
107
108
109 const G4int max_substeps = 100000000;
110 const G4int warn_substeps = 1000;
111
112
113
115
117
118#ifdef G4DEBUG_FIELD
120 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
121 if( (TrialPoint - StartPosition).mag() < tolerance * CLHEP::mm )
122 {
123 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
125 "Intersection point F is exactly at start point A." );
126 }
127#endif
128
129 do
130 {
133
134
135
136
139 CurrentB_PointVelocity,
140 CurrentE_Point,
142
143
144#ifdef G4DEBUG_FIELD
145 if( ApproxIntersecPointV.GetCurveLength() >
147 {
148 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
150 "Intermediate F point is past end B point!" );
151 }
152#endif
153
154 G4ThreeVector CurrentF_Point = ApproxIntersecPointV.GetPosition();
155
156
157
158
159 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
160
161 G4ThreeVector NewMomentumDir = ApproxIntersecPointV.GetMomentumDir();
162 G4double MomDir_dot_Norm = NewMomentumDir.
dot( NormalAtEntry ) ;
163
165
166#ifdef G4DEBUG_FIELD
169 NewMomentumDir, NormalAtEntry, validNormalAtE );
170#endif
171
172
173
174 G4bool adequate_angle = ( MomDir_dot_Norm >= 0.0 )
175 || (! validNormalAtE );
179 {
180 found_approximate_intersection = true;
181
182
183
184 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
185 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
186
188 {
189
190
192 G4ThreeVector MomentumDir= ApproxIntersecPointV.GetMomentumDirection();
194 CurrentE_Point, CurrentF_Point, MomentumDir,
195 last_AF_intersection, IP, NewSafety,
197
198 if ( goodCorrection )
199 {
200 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
201 IntersectedOrRecalculatedFT.SetPosition(IP);
202 }
203 }
204
205
206
207
208
209
210
211 }
212 else
213 {
214
215
216
217
218
220
223 G4bool usedNavigatorAF =
false;
225 CurrentF_Point,
226 NewSafety,
229 stepLengthAF,
230 PointG,
231 &usedNavigatorAF );
232 last_AF_intersection = Intersects_AF;
233 if( Intersects_AF )
234 {
235
236
237
238
239
240
241 CurrentB_PointVelocity = ApproxIntersecPointV;
242 CurrentE_Point = PointG;
243
244
245
246
247
248
249
252 validNormalAtE = validNormalLast;
253
254
255
256
257 final_section = false;
258
259#ifdef G4VERBOSE
261 {
262 G4cout <<
"G4PiF::LI> Investigating intermediate point"
263 << " at s=" << ApproxIntersecPointV.GetCurveLength()
264 << " on way to full s="
265 << CurveEndPointVelocity.GetCurveLength() <<
G4endl;
266 }
267#endif
268 }
269 else
270 {
271
272
273
274
276
279 G4bool usedNavigatorFB =
false;
280
281
282
283
287 stepLengthFB,
288 PointH, &usedNavigatorFB );
289 if( Intersects_FB )
290 {
291
292
293
294
295
296
297
298
299
300
301
302
303 CurrentA_PointVelocity = ApproxIntersecPointV;
304 CurrentE_Point = PointH;
305
306
307
308
309
310
311
314 validNormalAtE = validNormalLast;
315 }
316 else
317 {
318
319
320 if( final_section )
321 {
322
323
324
325
326
327
328 there_is_no_intersection = true;
329 }
330 else
331 {
332
333
334 CurrentA_PointVelocity = CurrentB_PointVelocity;
335 CurrentB_PointVelocity = CurveEndPointVelocity;
336 restoredFullEndpoint = true;
337 }
338 }
339 }
340
341
342
343
349
350
351
353 {
354
355
358 CurrentB_PointVelocity,
359 linDistSq,
360 curveDist );
362 CurrentB_PointVelocity = newEndPointFT;
363
364 if( (final_section))
365 {
366 recalculatedEndPoint = true;
367 IntersectedOrRecalculatedFT = newEndPointFT;
368
369 }
370 }
371 if( curveDist < 0.0 )
372 {
374 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
375 -1.0, NewSafety, substep_no );
376 std::ostringstream message;
377 message <<
"Error in advancing propagation." <<
G4endl
378 << " Point A (start) is " << CurrentA_PointVelocity
380 << " Point B (end) is " << CurrentB_PointVelocity
382 <<
" Curve distance is " << curveDist <<
G4endl
384 << "The final curve point is not further along"
385 <<
" than the original!" <<
G4endl;
386
387 if( recalculatedEndPoint )
388 {
389 message << "Recalculation of EndPoint was called with fEpsStep= "
391 }
392 message.precision(20);
393 message << " Point A (Curve start) is " << CurveStartPointVelocity
395 << " Point B (Curve end) is " << CurveEndPointVelocity
397 << " Point A (Current start) is " << CurrentA_PointVelocity
399 << " Point B (Current end) is " << CurrentB_PointVelocity
401 << " Point E (Trial Point) is " << CurrentE_Point
403 << " Point F (Intersection) is " << ApproxIntersecPointV
405 << " LocateIntersection parameters are : Substep no= "
406 << substep_no;
407
408 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
410 }
411
412 if ( restoredFullEndpoint )
413 {
414 final_section = restoredFullEndpoint;
415 restoredFullEndpoint = false;
416 }
417 }
418
419
420#ifdef G4DEBUG_LOCATE_INTERSECTION
421 G4int trigger_substepno_print= warn_substeps - 20;
422
423 if( substep_no >= trigger_substepno_print )
424 {
425 G4cout <<
"Difficulty in converging in "
426 << "G4SimpleLocator::EstimateIntersectionPoint():"
428 <<
" Substep no = " << substep_no <<
G4endl;
429 if( substep_no == trigger_substepno_print )
430 {
431 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
432 -1.0, NewSafety, 0);
433 }
434 G4cout <<
" State of point A: ";
435 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
436 -1.0, NewSafety, substep_no-1, 0);
437 G4cout <<
" State of point B: ";
438 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
439 -1.0, NewSafety, substep_no);
440 }
441#endif
442 ++substep_no;
443
444 } while ( ( ! found_approximate_intersection )
445 && ( ! there_is_no_intersection )
446 && ( substep_no <= max_substeps) );
447
448 if( substep_no > max_no_seen )
449 {
450 max_no_seen = substep_no;
451#ifdef G4DEBUG_LOCATE_INTERSECTION
452 if( max_no_seen > warn_substeps )
453 {
454 trigger_substepno_print = max_no_seen-20;
455 }
456#endif
457 }
458
459 if( ( substep_no >= max_substeps)
460 && !there_is_no_intersection
461 && !found_approximate_intersection )
462 {
463 G4cout <<
"ERROR - G4SimpleLocator::EstimateIntersectionPoint()" <<
G4endl
464 <<
" Start and Endpoint of Requested Step:" <<
G4endl;
465 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
466 -1.0, NewSafety, 0);
468 <<
" Start and end-point of current Sub-Step:" <<
G4endl;
469 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
470 -1.0, NewSafety, substep_no-1);
471 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
472 -1.0, NewSafety, substep_no);
473
474 std::ostringstream message;
475 message << "Convergence is requiring too many substeps: "
477 <<
" Abandoning effort to intersect." <<
G4endl
478 << " Found intersection = "
479 << found_approximate_intersection <<
G4endl
480 << " Intersection exists = "
481 << !there_is_no_intersection <<
G4endl;
482 message.precision(10);
484 G4double full_len = CurveEndPointVelocity.GetCurveLength();
485 message << " Undertaken only length: " << done_len
486 <<
" out of " << full_len <<
" required." <<
G4endl
487 << " Remaining length = " << full_len-done_len;
488
489 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
491 }
492 else if( substep_no >= warn_substeps )
493 {
494 std::ostringstream message;
495 message.precision(10);
496
497 message <<
"Many substeps while trying to locate intersection." <<
G4endl
498 << " Undertaken length: "
500 <<
" - Needed: " << substep_no <<
" substeps." <<
G4endl
501 << " Warning level = " << warn_substeps
502 << " and maximum substeps = " << max_substeps;
503 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
505 }
506 return !there_is_no_intersection;
507}
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
#define fPreviousSftOrigin
G4GLOB_DLL std::ostream G4cout
double dot(const Hep3Vector &) const
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector ¤tEPoint, G4double epsStep)
G4double GetCurveLength() const
G4ThreeVector GetPosition() const
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
G4double fiDeltaIntersection
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)
G4double GetEpsilonStepFor()
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
G4bool GetAdjustementOfFoundIntersection()
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, 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)