88{
89
90
91 G4bool found_approximate_intersection =
false;
92 G4bool there_is_no_intersection =
false;
93
94 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
95 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
97 G4bool validNormalAtE =
false;
99
100 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity);
102 G4bool last_AF_intersection =
false;
103 G4bool final_section =
true;
104
105 recalculatedEndPoint = false;
106
107 G4bool restoredFullEndpoint =
false;
108
109 G4int substep_no = 0;
110
111
112
113 const G4int max_substeps = 100000000;
114 const G4int warn_substeps = 1000;
115
116
117
118 static G4int max_no_seen= -1;
119
121
122#ifdef G4DEBUG_FIELD
124 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
125 if( (TrialPoint - StartPosition).mag() < tolerance * mm )
126 {
127 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
129 "Intersection point F is exactly at start point A." );
130 }
131#endif
132
133 do
134 {
137
138
139
140
143 CurrentB_PointVelocity,
144 CurrentE_Point,
146
147
148#ifdef G4DEBUG_FIELD
149 if( ApproxIntersecPointV.GetCurveLength() >
151 {
152 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
154 "Intermediate F point is past end B point!" );
155 }
156#endif
157
158 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
159
160
161
162
163 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
164
165 G4ThreeVector NewMomentumDir= ApproxIntersecPointV.GetMomentumDir();
166 G4double MomDir_dot_Norm= NewMomentumDir.
dot( NormalAtEntry ) ;
167
169
170#ifdef G4DEBUG_FIELD
173 NewMomentumDir, NormalAtEntry, validNormalAtE );
174#endif
175
176
177
178 G4bool adequate_angle = ( MomDir_dot_Norm >= 0.0 )
179 || (! validNormalAtE );
183 {
184 found_approximate_intersection = true;
185
186
187
188 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
189 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
190
192 {
193
194
196 G4ThreeVector MomentumDir= ApproxIntersecPointV.GetMomentumDirection();
198 CurrentE_Point, CurrentF_Point, MomentumDir,
199 last_AF_intersection, IP, NewSafety,
200 fPreviousSafety, fPreviousSftOrigin );
201
202 if(goodCorrection)
203 {
204 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
205 IntersectedOrRecalculatedFT.SetPosition(IP);
206 }
207 }
208
209
210
211
212
213
214
215 }
216 else
217 {
218
219
220
221
222
224
227 G4bool usedNavigatorAF =
false;
229 CurrentF_Point,
230 NewSafety,
231 fPreviousSafety,
232 fPreviousSftOrigin,
233 stepLengthAF,
234 PointG,
235 &usedNavigatorAF );
236 last_AF_intersection = Intersects_AF;
237 if( Intersects_AF )
238 {
239
240
241
242
243
244
245 CurrentB_PointVelocity = ApproxIntersecPointV;
246 CurrentE_Point = PointG;
247
248
249
250
251
252
253
256 validNormalAtE = validNormalLast;
257
258
259
260
261 final_section= false;
262
263#ifdef G4VERBOSE
265 {
266 G4cout <<
"G4PiF::LI> Investigating intermediate point"
267 << " at s=" << ApproxIntersecPointV.GetCurveLength()
268 << " on way to full s="
269 << CurveEndPointVelocity.GetCurveLength() <<
G4endl;
270 }
271#endif
272 }
273 else
274 {
275
276
277
278
280
283 G4bool usedNavigatorFB=
false;
284
285
286
287
289 NewSafety,fPreviousSafety,
290 fPreviousSftOrigin,
291 stepLengthFB,
292 PointH, &usedNavigatorFB );
293 if( Intersects_FB )
294 {
295
296
297
298
299
300
301
302
303
304
305
306
307 CurrentA_PointVelocity = ApproxIntersecPointV;
308 CurrentE_Point = PointH;
309
310
311
312
313
314
315
318 validNormalAtE = validNormalLast;
319 }
320 else
321 {
322
323
324 if( final_section )
325 {
326
327
328
329
330
331
332 there_is_no_intersection = true;
333 }
334 else
335 {
336
337
338 CurrentA_PointVelocity = CurrentB_PointVelocity;
339 CurrentB_PointVelocity = CurveEndPointVelocity;
340 restoredFullEndpoint = true;
341 }
342 }
343 }
344
345
346
347
353
354
355
357 {
358
359
362 CurrentB_PointVelocity,
363 linDistSq,
364 curveDist );
366 CurrentB_PointVelocity = newEndPointFT;
367
368 if( (final_section))
369 {
370 recalculatedEndPoint = true;
371 IntersectedOrRecalculatedFT = newEndPointFT;
372
373 }
374 }
375 if( curveDist < 0.0 )
376 {
378 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
379 -1.0, NewSafety, substep_no );
380 std::ostringstream message;
381 message <<
"Error in advancing propagation." <<
G4endl
382 << " Point A (start) is " << CurrentA_PointVelocity
384 << " Point B (end) is " << CurrentB_PointVelocity
386 <<
" Curve distance is " << curveDist <<
G4endl
388 << "The final curve point is not further along"
389 <<
" than the original!" <<
G4endl;
390
391 if( recalculatedEndPoint )
392 {
393 message << "Recalculation of EndPoint was called with fEpsStep= "
395 }
396 message.precision(20);
397 message << " Point A (Curve start) is " << CurveStartPointVelocity
399 << " Point B (Curve end) is " << CurveEndPointVelocity
401 << " Point A (Current start) is " << CurrentA_PointVelocity
403 << " Point B (Current end) is " << CurrentB_PointVelocity
405 << " Point E (Trial Point) is " << CurrentE_Point
407 << " Point F (Intersection) is " << ApproxIntersecPointV
409 << " LocateIntersection parameters are : Substep no= "
410 << substep_no;
411
412 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
414 }
415
416 if(restoredFullEndpoint)
417 {
418 final_section = restoredFullEndpoint;
419 restoredFullEndpoint = false;
420 }
421 }
422
423
424#ifdef G4DEBUG_LOCATE_INTERSECTION
425 static G4int trigger_substepno_print= warn_substeps - 20;
426
427 if( substep_no >= trigger_substepno_print )
428 {
429 G4cout <<
"Difficulty in converging in "
430 << "G4SimpleLocator::EstimateIntersectionPoint():"
432 <<
" Substep no = " << substep_no <<
G4endl;
433 if( substep_no == trigger_substepno_print )
434 {
435 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
436 -1.0, NewSafety, 0);
437 }
438 G4cout <<
" State of point A: ";
439 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
440 -1.0, NewSafety, substep_no-1, 0);
441 G4cout <<
" State of point B: ";
442 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
443 -1.0, NewSafety, substep_no);
444 }
445#endif
446 substep_no++;
447
448 } while ( ( ! found_approximate_intersection )
449 && ( ! there_is_no_intersection )
450 && ( substep_no <= max_substeps) );
451
452 if( substep_no > max_no_seen )
453 {
454 max_no_seen = substep_no;
455#ifdef G4DEBUG_LOCATE_INTERSECTION
456 if( max_no_seen > warn_substeps )
457 {
458 trigger_substepno_print = max_no_seen-20;
459 }
460#endif
461 }
462
463 if( ( substep_no >= max_substeps)
464 && !there_is_no_intersection
465 && !found_approximate_intersection )
466 {
467 G4cout <<
"ERROR - G4SimpleLocator::EstimateIntersectionPoint()" <<
G4endl
468 <<
" Start and Endpoint of Requested Step:" <<
G4endl;
469 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
470 -1.0, NewSafety, 0);
472 <<
" Start and end-point of current Sub-Step:" <<
G4endl;
473 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
474 -1.0, NewSafety, substep_no-1);
475 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
476 -1.0, NewSafety, substep_no);
477
478 std::ostringstream message;
479 message << "Convergence is requiring too many substeps: "
481 <<
" Abandoning effort to intersect." <<
G4endl
482 << " Found intersection = "
483 << found_approximate_intersection <<
G4endl
484 << " Intersection exists = "
485 << !there_is_no_intersection <<
G4endl;
486 message.precision(10);
488 G4double full_len = CurveEndPointVelocity.GetCurveLength();
489 message << " Undertaken only length: " << done_len
490 <<
" out of " << full_len <<
" required." <<
G4endl
491 << " Remaining length = " << full_len-done_len;
492
493 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
495 }
496 else if( substep_no >= warn_substeps )
497 {
498 std::ostringstream message;
499 message.precision(10);
500
501 message <<
"Many substeps while trying to locate intersection." <<
G4endl
502 << " Undertaken length: "
504 <<
" - Needed: " << substep_no <<
" substeps." <<
G4endl
505 << " Warning level = " << warn_substeps
506 << " and maximum substeps = " << max_substeps;
507 G4Exception(
"G4SimpleLocator::EstimateIntersectionPoint()",
509 }
510 return !there_is_no_intersection;
511}
G4DLLIMPORT 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=0)
G4double GetEpsilonStepFor()
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
G4bool GetAdjustementOfFoundIntersection()
G4bool AdjustmentOfFoundIntersection(const G4ThreeVector &A, const G4ThreeVector &CurrentE_Point, const G4ThreeVector &CurrentF_Point, const G4ThreeVector &MomentumDir, const G4bool IntersectAF, G4ThreeVector &IntersectionPoint, G4double &NewSafety, G4double &fPrevSafety, G4ThreeVector &fPrevSftOrigin)
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, G4double requestStep, G4double safety, G4int step)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)