101{
102
103
104 G4bool found_approximate_intersection =
false;
105 G4bool there_is_no_intersection =
false;
106
107 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
108 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
110 G4bool validNormalAtE =
false;
112
113 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity);
115 G4bool last_AF_intersection =
false;
116
117
118
119 G4bool first_section =
true;
120 recalculatedEndPoint = false;
121
122 G4bool restoredFullEndpoint =
false;
123
124 G4int substep_no = 0;
125
127
128
129
130 const G4int max_substeps= 10000;
131 const G4int warn_substeps= 1000;
132
133
134
135 static G4int max_no_seen= -1;
136
137
138
139
140
141
142
143
144
145
146
147 const G4int param_substeps=5;
148
150
151 G4bool Second_half =
false;
152
153
154
155
156
157
158
160
161#ifdef G4DEBUG_FIELD
162 static const G4double tolerance = 1.0e-8 * mm;
163 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
164 if( (TrialPoint - StartPosition).mag() < tolerance)
165 {
166 G4Exception(
"G4MultiLevelLocator::EstimateIntersectionPoint()",
168 "Intersection point F is exactly at start point A." );
169 }
170#endif
171
173
174
175
176
177
178 *ptrInterMedFT[0] = CurveEndPointVelocity;
179 for (
G4int idepth=1; idepth<max_depth+1; idepth++ )
180 {
181 *ptrInterMedFT[idepth]=CurveStartPointVelocity;
182 }
183
184
185
186 G4bool fin_section_depth[max_depth];
187 for (
G4int idepth=0; idepth<max_depth; idepth++ )
188 {
189 fin_section_depth[idepth]=true;
190 }
191
192
193 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
194
195 do
196 {
197 G4int substep_no_p = 0;
198 G4bool sub_final_section =
false;
199
200 SubStart_PointVelocity = CurrentA_PointVelocity;
201 do
202 {
205
206
207
208
211 CurrentB_PointVelocity,
212 CurrentE_Point,
214
215
216#ifdef G4DEBUG_FIELD
217 if( ApproxIntersecPointV.GetCurveLength() >
219 {
220 G4Exception(
"G4multiLevelLocator::EstimateIntersectionPoint()",
222 "Intermediate F point is past end B point" );
223 }
224#endif
225
226 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
227
228
229
230
231 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
232
233 G4ThreeVector NewMomentumDir= ApproxIntersecPointV.GetMomentumDir();
234 G4double MomDir_dot_Norm= NewMomentumDir.
dot( NormalAtEntry ) ;
235
236#ifdef G4DEBUG_FIELD
237 if( VerboseLevel > 3 )
238 {
242 MomDir_dot_ABchord = (1.0 / ABchord_length)
243 * NewMomentumDir.
dot( ChordAB );
245 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
246 G4cout <<
" dot( MomentumDir, ABchord_unit ) = "
247 << MomDir_dot_ABchord <<
G4endl;
248 }
249#endif
251 ( MomDir_dot_Norm >= 0.0 )
252 || (! validNormalAtE );
256 {
257 found_approximate_intersection = true;
258
259
260
261 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
262 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
263
265 {
266
267
269 G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
271 CurrentE_Point, CurrentF_Point, MomentumDir,
272 last_AF_intersection, IP, NewSafety,
273 previousSafety, previousSftOrigin );
274 if ( goodCorrection )
275 {
276 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
277 IntersectedOrRecalculatedFT.SetPosition(IP);
278 }
279 }
280
281
282
283
284
285
286 }
287 else
288 {
289
290
291
292
293
295
299 NewSafety, previousSafety,
300 previousSftOrigin,
301 stepLengthAF,
302 PointG );
303 last_AF_intersection = Intersects_AF;
304 if( Intersects_AF )
305 {
306
307
308
309
310
311
312 CurrentB_PointVelocity = ApproxIntersecPointV;
313 CurrentE_Point = PointG;
314
317 validNormalAtE = validNormalLast;
318
319
320
321
322 fin_section_depth[depth]=false;
323
324
325#ifdef G4VERBOSE
327 {
328 G4cout <<
"G4PiF::LI> Investigating intermediate point"
329 << " at s=" << ApproxIntersecPointV.GetCurveLength()
330 << " on way to full s="
331 << CurveEndPointVelocity.GetCurveLength() <<
G4endl;
332 }
333#endif
334 }
335 else
336 {
337
338
339
340
342
345
346
347
348
350 NewSafety, previousSafety,
351 previousSftOrigin,
352 stepLengthFB,
353 PointH );
354 if( Intersects_FB )
355 {
356
357
358
359
360
361
362
363
364
365
366
367
368 CurrentA_PointVelocity = ApproxIntersecPointV;
369 CurrentE_Point = PointH;
370
373 validNormalAtE = validNormalH;
374
375 }
376 else
377 {
378
379
380 if(fin_section_depth[depth])
381 {
382
383
384
385
386
387
388
389
390 if( ((Second_half)&&(depth==0)) || (first_section) )
391 {
392 there_is_no_intersection = true;
393 }
394 else
395 {
396
397
398
399 substep_no_p = param_substeps+2;
400
401
402
403
404 Second_half = true;
405 sub_final_section = true;
406
407 }
408 }
409 else
410 {
411 if(depth==0)
412 {
413
414
415 CurrentA_PointVelocity = CurrentB_PointVelocity;
416 CurrentB_PointVelocity = CurveEndPointVelocity;
417 SubStart_PointVelocity = CurrentA_PointVelocity;
418 restoredFullEndpoint = true;
419 }
420 else
421 {
422
423
424 CurrentA_PointVelocity = CurrentB_PointVelocity;
425 CurrentB_PointVelocity = *ptrInterMedFT[depth];
426 SubStart_PointVelocity = CurrentA_PointVelocity;
427 restoredFullEndpoint = true;
428 }
429 }
430 }
431 }
432
433
434
435
441
442
443
445 {
446
447
450 CurrentB_PointVelocity,
451 linDistSq,
452 curveDist );
454 CurrentB_PointVelocity = newEndPointFT;
455
456 if ( (fin_section_depth[depth])
457 &&( first_section || ((Second_half)&&(depth==0)) ) )
458 {
459 recalculatedEndPoint = true;
460 IntersectedOrRecalculatedFT = newEndPointFT;
461
462 }
463 }
464 if( curveDist < 0.0 )
465 {
467 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
468 -1.0, NewSafety, substep_no );
469 std::ostringstream message;
470 message <<
"Error in advancing propagation." <<
G4endl
471 << " Point A (start) is " << CurrentA_PointVelocity
473 << " Point B (end) is " << CurrentB_PointVelocity
475 <<
" Curve distance is " << curveDist <<
G4endl
477 << "The final curve point is not further along"
478 <<
" than the original!" <<
G4endl;
479
480 if( recalculatedEndPoint )
481 {
482 message << "Recalculation of EndPoint was called with fEpsStep= "
484 }
485 oldprc =
G4cerr.precision(20);
486 message << " Point A (Curve start) is " << CurveStartPointVelocity
488 << " Point B (Curve end) is " << CurveEndPointVelocity
490 << " Point A (Current start) is " << CurrentA_PointVelocity
492 << " Point B (Current end) is " << CurrentB_PointVelocity
494 << " Point S (Sub start) is " << SubStart_PointVelocity
496 << " Point E (Trial Point) is " << CurrentE_Point
498 << " Point F (Intersection) is " << ApproxIntersecPointV
500 << " LocateIntersection parameters are : Substep no= "
502 << " Substep depth no= "<< substep_no_p << " Depth= "
503 << depth;
505
506 G4Exception(
"G4MultiLevelLocator::EstimateIntersectionPoint()",
508 }
509 if(restoredFullEndpoint)
510 {
511 fin_section_depth[depth] = restoredFullEndpoint;
512 restoredFullEndpoint = false;
513 }
514 }
515
516
517#ifdef G4DEBUG_LOCATE_INTERSECTION
518 static G4int trigger_substepno_print= warn_substeps - 20 ;
519
520 if( substep_no >= trigger_substepno_print )
521 {
522 G4cout <<
"Difficulty in converging in "
523 << "G4MultiLevelLocator::EstimateIntersectionPoint():"
525 <<
" Substep no = " << substep_no <<
G4endl;
526 if( substep_no == trigger_substepno_print )
527 {
528 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
529 -1.0, NewSafety, 0);
530 }
531 G4cout <<
" State of point A: ";
532 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
533 -1.0, NewSafety, substep_no-1, 0);
534 G4cout <<
" State of point B: ";
535 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
536 -1.0, NewSafety, substep_no);
537 }
538#endif
539 substep_no++;
540 substep_no_p++;
541
542 } while ( ( ! found_approximate_intersection )
543 && ( ! there_is_no_intersection )
544 && ( substep_no_p <= param_substeps) );
545
546 first_section = false;
547
548 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
549 {
554
557
558
559
560 if( ( ( did_len )<fraction_done*all_len)
561 && (depth<max_depth) && (!sub_final_section) )
562 {
563 Second_half=false;
564 depth++;
565
566 G4double Sub_len = (all_len-did_len)/(2.);
571 *ptrInterMedFT[depth] = start;
572 CurrentB_PointVelocity = *ptrInterMedFT[depth];
573
574
575
576 SubStart_PointVelocity = CurrentA_PointVelocity;
577
578
579
582
585 NewSafety, previousSafety,
586 previousSftOrigin, stepLengthAB,
587 PointGe);
588 if( Intersects_AB )
589 {
590 last_AF_intersection = Intersects_AB;
591 CurrentE_Point = PointGe;
592 fin_section_depth[depth]=true;
593
596 validNormalAtE = validNormalAB;
597 }
598 else
599 {
600
601
602
603 Second_half = true;
604 }
605 }
606
607 if( (Second_half)&&(depth!=0) )
608 {
609
610
611
612 Second_half = true;
613
614
615 SubStart_PointVelocity = *ptrInterMedFT[depth];
616 CurrentA_PointVelocity = *ptrInterMedFT[depth];
617 CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
618
619
620
627 {
628
629
632 CurrentB_PointVelocity,
633 linDistSq,
634 curveDist );
636 CurrentB_PointVelocity = newEndPointFT;
637 if (depth==1)
638 {
639 recalculatedEndPoint = true;
640 IntersectedOrRecalculatedFT = newEndPointFT;
641
642 }
643 }
644
649 previousSafety,
650 previousSftOrigin, stepLengthAB,
651 PointGe);
652 if( Intersects_AB )
653 {
654 last_AF_intersection = Intersects_AB;
655 CurrentE_Point = PointGe;
656
659 validNormalAtE = validNormalAB;
660 }
661 depth--;
662 fin_section_depth[depth]=true;
663 }
664 }
665
666 } while ( ( ! found_approximate_intersection )
667 && ( ! there_is_no_intersection )
668 && ( substep_no <= max_substeps) );
669
670 if( substep_no > max_no_seen )
671 {
672 max_no_seen = substep_no;
673#ifdef G4DEBUG_LOCATE_INTERSECTION
674 if( max_no_seen > warn_substeps )
675 {
676 trigger_substepno_print = max_no_seen-20;
677 }
678#endif
679 }
680
681 if( ( substep_no >= max_substeps)
682 && !there_is_no_intersection
683 && !found_approximate_intersection )
684 {
685 G4cout <<
"ERROR - G4MultiLevelLocator::EstimateIntersectionPoint()"
687 <<
" Start and end-point of requested Step:" <<
G4endl;
688 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
689 -1.0, NewSafety, 0);
690 G4cout <<
" Start and end-point of current Sub-Step:" <<
G4endl;
691 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
692 -1.0, NewSafety, substep_no-1);
693 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
694 -1.0, NewSafety, substep_no);
695 std::ostringstream message;
696 message <<
"Too many substeps!" <<
G4endl
697 << " Convergence is requiring too many substeps: "
699 <<
" Abandoning effort to intersect. " <<
G4endl
700 << " Found intersection = "
701 << found_approximate_intersection <<
G4endl
702 << " Intersection exists = "
703 << !there_is_no_intersection <<
G4endl;
704
705#ifdef FUTURE_CORRECTION
706
707
708 if ( ! found_approximate_intersection )
709 {
710 recalculatedEndPoint = true;
711
712
713 IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
714
715 message << " Did not convergence after " << substep_no
717 << " The endpoint was adjused to pointA resulting"
719 << " from the last substep: " << CurrentA_PointVelocity
721 }
722#endif
723
724 oldprc =
G4cout.precision( 10 );
726 G4double full_len = CurveEndPointVelocity.GetCurveLength();
727 message << " Undertaken only length: " << done_len
728 <<
" out of " << full_len <<
" required." <<
G4endl
729 << " Remaining length = " << full_len - done_len;
730 G4cout.precision( oldprc );
731
732 G4Exception(
"G4MultiLevelLocator::EstimateIntersectionPoint()",
734 }
735 else if( substep_no >= warn_substeps )
736 {
737 oldprc =
G4cout.precision( 10 );
738 std::ostringstream message;
739 message << "Many substeps while trying to locate intersection."
741 << " Undertaken length: "
743 <<
" - Needed: " << substep_no <<
" substeps." <<
G4endl
744 << " Warning level = " << warn_substeps
745 << " and maximum substeps = " << max_substeps;
746 G4Exception(
"G4MultiLevelLocator::EstimateIntersectionPoint()",
748 G4cout.precision( oldprc );
749 }
750 return !there_is_no_intersection;
751}
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
double dot(const Hep3Vector &) const
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector ¤tEPoint, G4double epsStep)
G4MagInt_Driver * GetIntegrationDriver()
G4double GetCurveLength() const
G4ThreeVector GetPosition() const
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
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)