86{
87
88
89 G4bool found_approximate_intersection =
false;
90 G4bool there_is_no_intersection =
false;
91
92 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
93 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
95 G4bool validNormalAtE =
false;
97
98 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity);
100 G4bool last_AF_intersection =
false;
101 G4bool final_section =
true;
102
103 recalculatedEndPoint = false;
104
105 G4bool restoredFullEndpoint =
false;
106
107 G4int substep_no = 0;
108
109
110
111 const G4int max_substeps = 100000000;
112 const G4int warn_substeps = 1000;
113
114
115
117
119
120#ifdef G4DEBUG_FIELD
122 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
123 if( (TrialPoint - StartPosition).mag() < tolerance * CLHEP::mm )
124 {
125 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
127 "Intersection point F is exactly at start point A." );
128 }
129#endif
130
131 do
132 {
135
136
137
138
141 CurrentB_PointVelocity,
142 CurrentE_Point,
144
145
146#ifdef G4DEBUG_FIELD
147 if( ApproxIntersecPointV.GetCurveLength() >
149 {
150 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
152 "Intermediate F point is past end B point!" );
153 }
154#endif
155
156 G4ThreeVector CurrentF_Point = ApproxIntersecPointV.GetPosition();
157
158
159
160
161 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
162
163 G4ThreeVector NewMomentumDir = ApproxIntersecPointV.GetMomentumDir();
164 G4double MomDir_dot_Norm = NewMomentumDir.
dot( NormalAtEntry ) ;
165
167
168#ifdef G4DEBUG_FIELD
171 NewMomentumDir, NormalAtEntry, validNormalAtE );
172#endif
173
174
175
176 G4bool adequate_angle = ( MomDir_dot_Norm >= 0.0 )
177 || (! validNormalAtE );
181 {
182 found_approximate_intersection = true;
183
184
185
186 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
187 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
188
190 {
191
192
194 G4ThreeVector MomentumDir= ApproxIntersecPointV.GetMomentumDirection();
196 CurrentE_Point, CurrentF_Point, MomentumDir,
197 last_AF_intersection, IP, NewSafety,
199
200 if ( goodCorrection )
201 {
202 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
203 IntersectedOrRecalculatedFT.SetPosition(IP);
204 }
205 }
206
207
208
209
210
211
212
213 }
214 else
215 {
216
217
218
219
220
222
225 G4bool usedNavigatorAF =
false;
227 CurrentF_Point,
228 NewSafety,
231 stepLengthAF,
232 PointG,
233 &usedNavigatorAF );
234 last_AF_intersection = Intersects_AF;
235 if( Intersects_AF )
236 {
237
238
239
240
241
242
243 CurrentB_PointVelocity = ApproxIntersecPointV;
244 CurrentE_Point = PointG;
245
246
247
248
249
250
251
254 validNormalAtE = validNormalLast;
255
256
257
258
259 final_section = false;
260
261#ifdef G4VERBOSE
263 {
264 G4cout <<
"G4PiF::LI> Investigating intermediate point"
265 << " at s=" << ApproxIntersecPointV.GetCurveLength()
266 << " on way to full s="
267 << CurveEndPointVelocity.GetCurveLength() <<
G4endl;
268 }
269#endif
270 }
271 else
272 {
273
274
275
276
278
281 G4bool usedNavigatorFB =
false;
282
283
284
285
289 stepLengthFB,
290 PointH, &usedNavigatorFB );
291 if( Intersects_FB )
292 {
293
294
295
296
297
298
299
300
301
302
303
304
305 CurrentA_PointVelocity = ApproxIntersecPointV;
306 CurrentE_Point = PointH;
307
308
309
310
311
312
313
316 validNormalAtE = validNormalLast;
317 }
318 else
319 {
320
321
322 if( final_section )
323 {
324
325
326
327
328
329
330 there_is_no_intersection = true;
331 }
332 else
333 {
334
335
336 CurrentA_PointVelocity = CurrentB_PointVelocity;
337 CurrentB_PointVelocity = CurveEndPointVelocity;
338 restoredFullEndpoint = true;
339 }
340 }
341 }
342
343
344
345
351
352
353
355 {
356
357
360 CurrentB_PointVelocity,
361 linDistSq,
362 curveDist );
364 CurrentB_PointVelocity = newEndPointFT;
365
366 if( (final_section))
367 {
368 recalculatedEndPoint = true;
369 IntersectedOrRecalculatedFT = newEndPointFT;
370
371 }
372 }
373 if( curveDist < 0.0 )
374 {
376 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
377 -1.0, NewSafety, substep_no );
378 std::ostringstream message;
379 message <<
"Error in advancing propagation." <<
G4endl
380 << " Point A (start) is " << CurrentA_PointVelocity
382 << " Point B (end) is " << CurrentB_PointVelocity
384 <<
" Curve distance is " << curveDist <<
G4endl
386 << "The final curve point is not further along"
387 <<
" than the original!" <<
G4endl;
388
389 if( recalculatedEndPoint )
390 {
391 message << "Recalculation of EndPoint was called with fEpsStep= "
393 }
394 message.precision(20);
395 message << " Point A (Curve start) is " << CurveStartPointVelocity
397 << " Point B (Curve end) is " << CurveEndPointVelocity
399 << " Point A (Current start) is " << CurrentA_PointVelocity
401 << " Point B (Current end) is " << CurrentB_PointVelocity
403 << " Point E (Trial Point) is " << CurrentE_Point
405 << " Point F (Intersection) is " << ApproxIntersecPointV
407 << " LocateIntersection parameters are : Substep no= "
408 << substep_no;
409
410 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
412 }
413
414 if ( restoredFullEndpoint )
415 {
416 final_section = restoredFullEndpoint;
417 restoredFullEndpoint = false;
418 }
419 }
420
421
422#ifdef G4DEBUG_LOCATE_INTERSECTION
423 G4int trigger_substepno_print= warn_substeps - 20;
424
425 if( substep_no >= trigger_substepno_print )
426 {
427 G4cout <<
"Difficulty in converging in "
428 << "G4SimpleLocator::EstimateIntersectionPoint():"
430 <<
" Substep no = " << substep_no <<
G4endl;
431 if( substep_no == trigger_substepno_print )
432 {
433 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
434 -1.0, NewSafety, 0);
435 }
436 G4cout <<
" State of point A: ";
437 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
438 -1.0, NewSafety, substep_no-1, 0);
439 G4cout <<
" State of point B: ";
440 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
441 -1.0, NewSafety, substep_no);
442 }
443#endif
444 ++substep_no;
445
446 } while ( ( ! found_approximate_intersection )
447 && ( ! there_is_no_intersection )
448 && ( substep_no <= max_substeps) );
449
450 if( substep_no > max_no_seen )
451 {
452 max_no_seen = substep_no;
453#ifdef G4DEBUG_LOCATE_INTERSECTION
454 if( max_no_seen > warn_substeps )
455 {
456 trigger_substepno_print = max_no_seen-20;
457 }
458#endif
459 }
460
461 if( ( substep_no >= max_substeps)
462 && !there_is_no_intersection
463 && !found_approximate_intersection )
464 {
465 G4cout <<
"ERROR - G4SimpleLocator::EstimateIntersectionPoint()" <<
G4endl
466 <<
" Start and Endpoint of Requested Step:" <<
G4endl;
467 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
468 -1.0, NewSafety, 0);
470 <<
" Start and end-point of current Sub-Step:" <<
G4endl;
471 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
472 -1.0, NewSafety, substep_no-1);
473 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
474 -1.0, NewSafety, substep_no);
475
476 std::ostringstream message;
477 message << "Convergence is requiring too many substeps: "
479 <<
" Abandoning effort to intersect." <<
G4endl
480 << " Found intersection = "
481 << found_approximate_intersection <<
G4endl
482 << " Intersection exists = "
483 << !there_is_no_intersection <<
G4endl;
484 message.precision(10);
486 G4double full_len = CurveEndPointVelocity.GetCurveLength();
487 message << " Undertaken only length: " << done_len
488 <<
" out of " << full_len <<
" required." <<
G4endl
489 << " Remaining length = " << full_len-done_len;
490
491 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
493 }
494 else if( substep_no >= warn_substeps )
495 {
496 std::ostringstream message;
497 message.precision(10);
498
499 message <<
"Many substeps while trying to locate intersection." <<
G4endl
500 << " Undertaken length: "
502 <<
" - Needed: " << substep_no <<
" substeps." <<
G4endl
503 << " Warning level = " << warn_substeps
504 << " and maximum substeps = " << max_substeps;
505 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
507 }
508 return !there_is_no_intersection;
509}
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)