| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Class G4VIntersectionLocator implementation
27 //
28 // 27.10.08 - John Apostolakis, Tatiana Nikiti
29 // -------------------------------------------
30
31 #include <iomanip>
32 #include <sstream>
33
34 #include "globals.hh"
35 #include "G4ios.hh"
36 #include "G4AutoDelete.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4VIntersectionLocator.hh"
39 #include "G4TouchableHandle.hh"
40 #include "G4GeometryTolerance.hh"
41
42 //////////////////////////////////////////////
43 //
44 // Constructor
45 //
46 G4VIntersectionLocator::G4VIntersectionLocator
47 : fiNavigator(theNavigator)
48 {
49 kCarTolerance = G4GeometryTolerance::GetInst
50
51 if( fiNavigator->GetExternalNavigation() ==
52 {
53 fHelpingNavigator = new G4Navigator();
54 }
55 else // Must clone the navigator, together w
56 {
57 fHelpingNavigator = fiNavigator->Clone();
58 }
59 }
60
61 //////////////////////////////////////////////
62 //
63 // Destructor.
64 //
65 G4VIntersectionLocator::~G4VIntersectionLocato
66 {
67 delete fHelpingNavigator;
68 delete fpTouchable;
69 }
70
71 //////////////////////////////////////////////
72 //
73 // Dump status of propagator to cout (old meth
74 //
75 void
76 G4VIntersectionLocator::printStatus( const G4F
77 const G4F
78 G4d
79 G4d
80 G4i
81 {
82 std::ostringstream os;
83 printStatus( StartFT,CurrentFT,requestStep,s
84 G4cout << os.str();
85 }
86
87 //////////////////////////////////////////////
88 //
89 // Dumps status of propagator.
90 //
91 void
92 G4VIntersectionLocator::printStatus( const G4F
93 const G4F
94 G4d
95 G4d
96 G4i
97 std
98 G4i
99 {
100 // const G4int verboseLevel= fVerboseLevel;
101 const G4ThreeVector StartPosition = St
102 const G4ThreeVector StartUnitVelocity = St
103 const G4ThreeVector CurrentPosition = Cu
104 const G4ThreeVector CurrentUnitVelocity = Cu
105
106 G4double step_len = CurrentFT.GetCurveLength
107 G4long oldprc; // cout/cerr precision setti
108
109 if( ((stepNo == 0) && (verboseLevel <3)) ||
110 {
111 oldprc = os.precision(4);
112 os << std::setw( 6) << " "
113 << std::setw( 25) << " Current Posi
114 << G4endl;
115 os << std::setw( 5) << "Step#"
116 << std::setw(10) << " s " << " "
117 << std::setw(10) << "X(mm)" << " "
118 << std::setw(10) << "Y(mm)" << " "
119 << std::setw(10) << "Z(mm)" << " "
120 << std::setw( 7) << " N_x " << " "
121 << std::setw( 7) << " N_y " << " "
122 << std::setw( 7) << " N_z " << " "
123 os << std::setw( 7) << " Delta|N|" << " "
124 << std::setw( 9) << "StepLen" << "
125 << std::setw(12) << "StartSafety" <
126 << std::setw( 9) << "PhsStep" << "
127 os << G4endl;
128 os.precision(oldprc);
129 }
130 if((stepNo == 0) && (verboseLevel <=3))
131 {
132 // Recurse to print the start values
133 //
134 printStatus( StartFT, StartFT, -1.0, safet
135 }
136 if( verboseLevel <= 3 )
137 {
138 if( stepNo >= 0)
139 {
140 os << std::setw( 4) << stepNo << " ";
141 }
142 else
143 {
144 os << std::setw( 5) << "Start" ;
145 }
146 oldprc = os.precision(8);
147 os << std::setw(10) << CurrentFT.GetCurveL
148 os << std::setw(10) << CurrentPosition.x()
149 << std::setw(10) << CurrentPosition
150 << std::setw(10) << CurrentPosition
151 os.precision(4);
152 os << std::setw( 7) << CurrentUnitVelocity
153 << std::setw( 7) << CurrentUnitVelo
154 << std::setw( 7) << CurrentUnitVelo
155 os.precision(3);
156 os << std::setw( 7)
157 << CurrentFT.GetMomentum().mag()- S
158 << " ";
159 os << std::setw( 9) << step_len << " ";
160 os << std::setw(12) << safety << " ";
161 if( requestStep != -1.0 )
162 {
163 os << std::setw( 9) << requestStep << "
164 }
165 else
166 {
167 os << std::setw( 9) << "Init/NotKnown" <
168 }
169 os << G4endl;
170 os.precision(oldprc);
171 }
172 else // if( verboseLevel > 3 )
173 {
174 // Multi-line output
175
176 os << "Step taken was " << step_len
177 << " out of PhysicalStep= " << req
178 os << "Final safety is: " << safety << G4e
179 os << "Chord length = " << (CurrentPositio
180 << G4endl;
181 os << G4endl;
182 }
183 }
184
185 //////////////////////////////////////////////
186 //
187 // ReEstimateEndPoint.
188 //
189 G4FieldTrack G4VIntersectionLocator::
190 ReEstimateEndpoint( const G4FieldTrack& Curren
191 const G4FieldTrack& Estima
192 G4double , // l
193 G4double
194 #ifdef G4DEBUG_FIELD
195 curveDist
196 #endif
197 )
198 {
199 G4FieldTrack newEndPoint( CurrentStateA );
200 auto integrDriver = GetChordFinderFor()->Get
201
202 G4FieldTrack retEndPoint( CurrentStateA );
203 G4bool goodAdvance;
204 G4int itrial = 0;
205 const G4int no_trials = 20;
206
207
208 G4double endCurveLen= EstimatedEndStateB.Get
209
210 do // Loop checking, 07.10.2016, JA
211 {
212 G4double currentCurveLen = newEndPoint.Get
213 G4double advanceLength = endCurveLen - cur
214 if (std::abs(advanceLength)<kCarTolerance)
215 {
216 goodAdvance=true;
217 }
218 else
219 {
220 goodAdvance = integrDriver->AccurateAdva
221
222 }
223 }
224 while( !goodAdvance && (++itrial < no_trials
225
226 if( goodAdvance )
227 {
228 retEndPoint = newEndPoint;
229 }
230 else
231 {
232 retEndPoint = EstimatedEndStateB; // Could
233 }
234
235 // All the work is done
236 // below are some diagnostics only -- befor
237 //
238 const G4String MethodName("G4VIntersectionLo
239
240 #ifdef G4VERBOSE
241 G4int latest_good_trials = 0;
242 if( itrial > 1)
243 {
244 if( fVerboseLevel > 0 )
245 {
246 G4cout << MethodName << " called - goodA
247 << " trials = " << itrial
248 << " previous good= " << latest_g
249 << G4endl;
250 }
251 latest_good_trials = 0;
252 }
253 else
254 {
255 ++latest_good_trials;
256 }
257 #endif
258
259 #ifdef G4DEBUG_FIELD
260 G4double lengthDone = newEndPoint.GetCurveLe
261 - CurrentStateA.GetCurve
262 if( !goodAdvance )
263 {
264 if( fVerboseLevel >= 3 )
265 {
266 G4cout << MethodName << "> AccurateAdvan
267 G4cout << " in " << itrial << " integrat
268 G4cout << " It went only " << lengthDone
269 << " -- a difference of " << curv
270 G4cout << " ReEstimateEndpoint> Reset en
271 << G4endl;
272 }
273 }
274 G4double linearDist = ( EstimatedEndStateB.G
275 - CurrentStateA.GetPosit
276 static G4int noInaccuracyWarnings = 0;
277 G4int maxNoWarnings = 10;
278 if ( (noInaccuracyWarnings < maxNoWarnings
279 || (fVerboseLevel > 1) )
280 {
281 G4ThreeVector move = newEndPoint.GetPosi
282 - EstimatedEndStateB.
283 std::ostringstream message;
284 message.precision(12);
285 message << " Integration inaccuracy requ
286 << " an adjustment in the step's
287 << " Two mid-points are furthe
288 << " curve length difference"
289 << " Dist = " << linearD
290 << " curve length = " << curveDi
291 message << " Correction applied is " <<
292 << " Old Estimated B position=
293 << EstimatedEndStateB.GetPositio
294 << " Recalculated Position=
295 << newEndPoint.GetPosition() <<
296 << " Change ( new - old ) =
297 G4Exception("G4VIntersectionLocator::ReE
298 "GeomNav1002", JustWarning,
299 }
300 /*
301 #else
302 // Statistics on the RMS value of the correc
303
304 static G4ThreadLocal G4int noCorrections = 0
305 ++noCorrections;
306 if( goodAdvance )
307 {
308 static G4ThreadLocal G4double sumCorrectio
309 sumCorrectionsSq += (EstimatedEndStateB.Ge
310 newEndPoint.GetPositi
311 }
312 */
313 #endif
314
315 return retEndPoint;
316 }
317
318
319 //////////////////////////////////////////////
320 //
321 // ReEstimateEndPoint.
322 //
323 // The following values are returned in curv
324 // 0: Normal - no problem
325 // 1: Unexpected co-incidence - milder m
326 // 2: Real mixup - EndB is NOT past Star
327 // ( ie. StartA's curve-lengh is bi
328
329
330 G4bool G4VIntersectionLocator::
331 CheckAndReEstimateEndpoint( const G4FieldTrack
332 const G4FieldTrack
333 G4FieldTrack
334 G4int&
335 {
336 G4double linDistSq, curveDist;
337
338 G4bool recalculated = false;
339 curveError= 0;
340
341 linDistSq = ( EstimatedEndB.GetPosition()
342 - CurrentStartA.GetPosition() ).
343 curveDist = EstimatedEndB.GetCurveLength()
344 - CurrentStartA.GetCurveLength();
345 if( (curveDist>=0.0)
346 && (curveDist*curveDist *(1.0+2.0*fiEpsil
347 {
348 G4FieldTrack newEndPointFT = EstimatedEndB
349
350 if (curveDist>0.0)
351 {
352 // Re-integrate to obtain a new B
353 RevisedEndPoint = ReEstimateEndpoint( C
354 E
355 l
356 c
357 recalculated = true;
358 }
359 else
360 {
361 // Zero length -> no advance!
362 newEndPointFT = CurrentStartA;
363 recalculated = true;
364 curveError = 1; // Unexpected co-incid
365
366 G4Exception("G4MultiLevelLocator::Estim
367 "GeomNav1002", JustWarning,
368 "A & B are at equal distance in 2nd
369 }
370 }
371
372 // Sanity check
373 //
374 if( curveDist < 0.0 )
375 {
376 curveError = 2; // Real mixup
377 }
378 return recalculated;
379 }
380
381 //////////////////////////////////////////////
382 //
383 // Method for finding SurfaceNormal of Interse
384 //
385 G4ThreeVector G4VIntersectionLocator::
386 GetLocalSurfaceNormal(const G4ThreeVector& Cur
387 {
388 G4ThreeVector Normal(G4ThreeVector(0.0,0.0,0
389 G4VPhysicalVolume* located;
390
391 validNormal = false;
392 fHelpingNavigator->SetWorldVolume(GetNavigat
393 located = fHelpingNavigator->LocateGlobalPoi
394
395 delete fpTouchable;
396 fpTouchable = fHelpingNavigator->CreateTouch
397
398 // To check if we can use GetGlobalExitNorma
399 //
400 G4ThreeVector localPosition = fpTouchable->G
401 ->GetTopTransform().TransformP
402
403 // Issue: in the case of coincident surfaces
404 // which side you are located onto (c
405 // TO-DO: use direction (of chord) to identi
406
407 if( located != nullptr)
408 {
409 G4LogicalVolume* pLogical= located->GetLog
410 G4VSolid* pSolid;
411
412 if( (pLogical != nullptr) && ( (pSolid=pLo
413 {
414 if ( ( pSolid->Inside(localPosition)==kS
415 || ( pSolid->DistanceToOut(localPositi
416 {
417 Normal = pSolid->SurfaceNormal(localPo
418 validNormal = true;
419
420 #ifdef G4DEBUG_FIELD
421 if( std::fabs(Normal.mag2() - 1.0 ) >
422 {
423 G4cerr << "PROBLEM in G4VIntersectio
424 << G4endl;
425 G4cerr << " Normal is not unit - ma
426 G4cerr << " at trial local point "
427 G4cerr << " Solid is " << *pSolid
428 }
429 #endif
430 }
431 }
432 }
433 return Normal;
434 }
435
436 //////////////////////////////////////////////
437 //
438 // Adjustment of Found Intersection
439 //
440 G4bool G4VIntersectionLocator::
441 AdjustmentOfFoundIntersection( const G4ThreeVe
442 const G4ThreeVe
443 const G4ThreeVe
444 const G4ThreeVe
445 const G4bool
446 G4ThreeVe
447 G4double&
448 G4double&
449 G4ThreeVe
450 {
451 G4double dist,lambda;
452 G4ThreeVector Normal, NewPoint, Point_G;
453 G4bool goodAdjust = false, Intersects_FP = f
454
455 // Get SurfaceNormal of Intersecting Solid
456 //
457 Normal = GetGlobalSurfaceNormal(CurrentE_Poi
458 if(!validNormal) { return false; }
459
460 // Intersection between Line and Plane
461 //
462 G4double n_d_m = Normal.dot(MomentumDir);
463 if ( std::abs(n_d_m)>kCarTolerance )
464 {
465 #ifdef G4VERBOSE
466 if ( fVerboseLevel>1 )
467 {
468 G4Exception("G4VIntersectionLocator::Adj
469 "GeomNav0003", JustWarning,
470 "No intersection. Parallels
471 }
472 #endif
473 lambda =- Normal.dot(CurrentF_Point-Curren
474
475 // New candidate for Intersection
476 //
477 NewPoint = CurrentF_Point+lambda*MomentumD
478
479 // Distance from CurrentF to Calculated In
480 //
481 dist = std::abs(lambda);
482
483 if ( dist<kCarTolerance*0.001 ) { return
484
485 // Calculation of new intersection point o
486 //
487 if ( IntersectAF ) // First part interse
488 {
489 G4double stepLengthFP;
490 G4ThreeVector Point_P = CurrentA_Point;
491 GetNavigatorFor()->LocateGlobalPointWith
492 Intersects_FP = IntersectChord( Point_P,
493 fPreviou
494 stepLeng
495
496 }
497 else // Second part intersects
498 {
499 G4double stepLengthFP;
500 GetNavigatorFor()->LocateGlobalPointWith
501 Intersects_FP = IntersectChord( CurrentF
502 fPreviou
503 stepLeng
504 }
505 if ( Intersects_FP )
506 {
507 goodAdjust = true;
508 IntersectionPoint = Point_G;
509 }
510 }
511
512 return goodAdjust;
513 }
514
515 //////////////////////////////////////////////
516 //
517 // GetSurfaceNormal.
518 //
519 G4ThreeVector G4VIntersectionLocator::
520 GetSurfaceNormal(const G4ThreeVector& CurrentI
521 G4bool& validNormal)
522 {
523 G4ThreeVector NormalAtEntry; // ( -10. , -10
524
525 G4ThreeVector NormalAtEntryLast, NormalAtEnt
526 G4bool validNormalLast;
527
528 // Relies on a call to Navigator::ComputeSte
529 // this call
530 //
531 NormalAtEntryLast = GetLastSurfaceNormal( Cu
532 // May return valid=false in cases, includ
533 // - if the candidate volume was not foun
534 // - a replica was involved -- determined
535 // (This list is not complete.)
536
537 #ifdef G4DEBUG_FIELD
538 if ( validNormalLast
539 && ( std::fabs(NormalAtEntryLast.mag2() - 1
540 {
541 std::ostringstream message;
542 message << "PROBLEM: Normal is not unit -
543 << NormalAtEntryLast.mag() << G4en
544 message << " at trial intersection point
545 message << " Obtained from Get *Last* Su
546 G4Exception("G4VIntersectionLocator::GetSu
547 "GeomNav1002", JustWarning, me
548 }
549 #endif
550
551 if( validNormalLast )
552 {
553 NormalAtEntry = NormalAtEntryLast;
554 }
555 validNormal = validNormalLast;
556
557 return NormalAtEntry;
558 }
559
560 //////////////////////////////////////////////
561 //
562 // GetGlobalSurfaceNormal.
563 //
564 G4ThreeVector G4VIntersectionLocator::
565 GetGlobalSurfaceNormal(const G4ThreeVector& Cu
566 G4bool& validNorm
567 {
568 G4ThreeVector localNormal = GetLocalSurfaceN
569 G4AffineTransform localToGlobal = /
570 fHelpingNavigator->GetLocalToGlobalTrans
571 G4ThreeVector globalNormal = localToGlobal.T
572
573 #ifdef G4DEBUG_FIELD
574 if( validNormal && ( std::fabs(globalNormal.
575 {
576 std::ostringstream message;
577 message << "******************************
578 << G4endl;
579 message << " Bad Normal in G4VIntersection
580 << G4endl;
581 message << " * Constituents: " << G4endl;
582 message << " Local Normal= " << localN
583 message << " Transform: " << G4endl
584 << " Net Translation= " << lo
585 << G4endl
586 << " Net Rotation = " << lo
587 << G4endl;
588 message << " * Result: " << G4endl;
589 message << " Global Normal= " << local
590 message << "******************************
591 G4Exception("G4VIntersectionLocator::GetGl
592 "GeomNav1002", JustWarning, me
593 }
594 #endif
595
596 return globalNormal;
597 }
598
599 //////////////////////////////////////////////
600 //
601 // GetLastSurfaceNormal.
602 //
603 G4ThreeVector G4VIntersectionLocator::
604 GetLastSurfaceNormal( const G4ThreeVector& int
605 G4bool& normalIsVa
606 {
607 G4ThreeVector normalVec;
608 G4bool validNorm;
609 normalVec = fiNavigator->GetGlobalExitNormal
610 normalIsValid = validNorm;
611
612 return normalVec;
613 }
614
615 //////////////////////////////////////////////
616 //
617 // ReportTrialStep.
618 //
619 void G4VIntersectionLocator::ReportTrialStep(
620 const
621 const
622 const
623 const
624
625 {
626 G4double ABchord_length = ChordAB_v.mag();
627 G4double MomDir_dot_Norm = NewMomentumDir.do
628 G4double MomDir_dot_ABchord;
629 MomDir_dot_ABchord = (1.0 / ABchord_length)
630
631 std::ostringstream outStream;
632 outStream << std::setw(6) << " Step# "
633 << std::setw(17) << " |ChordEF|(mag)" << "
634 << std::setw(18) << " uMomentum.Normal" <<
635 << std::setw(18) << " uMomentum.ABdir " <<
636 << std::setw(16) << " AB-dist " <<
637 << " Chord Vector (EF) "
638 << G4endl;
639 outStream.precision(7);
640 outStream << " " << std::setw(5) << step_no
641 << " " << std::setw(18) << ChordEF_v.mag()
642 << " " << std::setw(18) << MomDir_dot_Norm
643 << " " << std::setw(18) << MomDir_dot_ABch
644 << " " << std::setw(12) << ABchord_length
645 << " " << ChordEF_v
646 << G4endl;
647 outStream << " MomentumDir= " << " " << NewM
648 << " Normal at Entry E= " << NormalAtEntry
649 << " AB chord = " << ChordAB_v
650 << G4endl;
651 G4cout << outStream.str();
652
653 if( ( std::fabs(NormalAtEntry.mag2() - 1.0)
654 {
655 std::ostringstream message;
656 message << "Normal is not unit - mag= " <<
657 << " ValidNormalAtE = " <<
658 G4Exception("G4VIntersectionLocator::Repor
659 "GeomNav1002", JustWarning, me
660 }
661 return;
662 }
663
664 //////////////////////////////////////////////
665 //
666 // LocateGlobalPointWithinVolumeAndCheck.
667 //
668 // Locate point using navigator: updates state
669 // By default, it assumes that the point is in
670 // and returns true.
671 // In check mode, checks whether the point is
672 // If it is inside, it returns true
673 // If not, issues a warning and returns fals
674 //
675 G4bool G4VIntersectionLocator::
676 LocateGlobalPointWithinVolumeAndCheck( const G
677 {
678 G4bool good = true;
679 G4Navigator* nav = GetNavigatorFor();
680 const G4String
681 MethodName("G4VIntersectionLocator::LocateGl
682
683 if( fCheckMode )
684 {
685 G4bool navCheck= nav->IsCheckModeActive();
686 nav->CheckMode(true);
687
688 // Identify the current volume
689
690 G4TouchableHandle startTH= nav->CreateTouc
691 G4VPhysicalVolume* motherPhys = startTH->G
692 G4VSolid* motherSolid = startTH->
693 G4AffineTransform transform = nav->GetGlob
694 G4int motherCopyNo = motherPhys->GetCopyNo
695
696 // Let's check that the point is inside th
697 G4ThreeVector localPosition = transform.T
698 EInside inMother = motherSolid->Ins
699 if( inMother != kInside )
700 {
701 std::ostringstream message;
702 message << "Position located "
703 << ( inMother == kSurface ? " on
704 << "expected volume" << G4endl
705 << " Safety (from Outside) = "
706 << motherSolid->DistanceToIn(loc
707 G4Exception("G4VIntersectionLocator::Loc
708 "GeomNav1002", JustWarning,
709 }
710
711 // 1. Simple next step - quick relocation
712 // nav->LocateGlobalPointWithinVolume( pos
713
714 // 2. Full relocation - to cross-check ans
715 G4VPhysicalVolume* nextPhysical= nav->Loca
716 if( (nextPhysical != motherPhys)
717 || (nextPhysical->GetCopyNo() != mothe
718 )
719 {
720 G4Exception("G4VIntersectionLocator::Loc
721 "GeomNav1002", JustWarning,
722 "Position located outside ex
723 }
724 nav->CheckMode(navCheck); // Recover orig
725 }
726 else
727 {
728 nav->LocateGlobalPointWithinVolume( positi
729 }
730 return good;
731 }
732
733 //////////////////////////////////////////////
734 //
735 // LocateGlobalPointWithinVolumeCheckAndReport
736 //
737 void G4VIntersectionLocator::
738 LocateGlobalPointWithinVolumeCheckAndReport( c
739 c
740 G
741 {
742 // Save value of Check mode first
743 G4bool oldCheck = GetCheckMode();
744
745 G4bool ok = LocateGlobalPointWithinVolumeAnd
746 if( !ok )
747 {
748 std::ostringstream message;
749 message << "Failed point location." << G4e
750 << " Code Location info: " << Co
751 G4Exception("G4VIntersectionLocator::Locat
752 "GeomNav1002", JustWarning, me
753 }
754
755 SetCheckMode( oldCheck );
756 }
757
758 //////////////////////////////////////////////
759 //
760 // ReportReversedPoints.
761 //
762 void G4VIntersectionLocator::
763 ReportReversedPoints( std::ostringstream& msg,
764 const G4FieldTrack& Star
765 const G4FieldTrack& EndP
766 G4double NewSafety
767 const G4FieldTrack& A_Pt
768 const G4FieldTrack& B_Pt
769 const G4FieldTrack& SubS
770 const G4ThreeVector& E_P
771 const G4FieldTrack& Appr
772 G4int substep_no,
773 {
774 // Expect that 'msg' can hold the name of t
775
776 // FieldTrack 'points' A and B have been ta
777 // Whereas A should be before B, it is foun
778 G4int verboseLevel= 5;
779 G4double curveDist = B_PtVel.GetCurveLength
780 G4VIntersectionLocator::printStatus( A_PtVe
781 -1.0, NewSafety, s
782 msg << "Error in advancing propagation." <<
783 << " The final curve point is NOT fur
784 << " than the original!" << G4endl
785 << " Going *backwards* from len(A) =
786 << " to len(B) = " << B_PtVel.GetCurve
787 << " Curve distance is " << curveD
788 << G4endl
789 << " Point A' (start) is " << A_Pt
790 << " Point B' (end) is " << B_Pt
791 msg << " fEpsStep= " << epsStep << G4e
792
793 G4long oldprc = msg.precision(20);
794 msg << " In full precision, the position, m
795 << " ... are: " << G4endl;
796 msg << " Point A[0] (Curve start) is " <<
797 << " Point S (Sub start) is " <<
798 << " Point A' (Current start) is " <<
799 << " Point E (Trial Point) is " <<
800 << " Point F (Intersection) is " <<
801 << " Point B' (Current end) is " <<
802 << " Point B[0] (Curve end) is " <<
803 << G4endl
804 << " LocateIntersection parameters are
805 << " Substep no (total) = " << su
806 << " Substep no = " << su
807 msg.precision(oldprc);
808 }
809
810 //////////////////////////////////////////////
811 //
812 // ReportProgress.
813 //
814 void G4VIntersectionLocator::ReportProgress( s
815 const G4Fi
816 const G4Fi
817 G4in
818 const G4Fi
819 const G4Fi
820 G4do
821 G4in
822
823 {
824 oss << "ReportProgress: Current status of in
825 if( depth > 0 ) { oss << " Depth= " << depth
826 oss << " Substep no = " << substep_no << G4e
827 G4int verboseLevel = 5;
828 G4double safetyPrev = -1.0; // Add as argum
829
830 printStatus( StartPointVel, EndPointVel, -1.
831 oss, verboseLevel);
832 oss << " * Start and end-point of requested
833 oss << " ** State of point A: ";
834 printStatus( A_PtVel, A_PtVel, -1.0, safetyP
835 oss, verboseLevel);
836 oss << " ** State of point B: ";
837 printStatus( A_PtVel, B_PtVel, -1.0, safetyL
838 oss, verboseLevel);
839 }
840
841 //////////////////////////////////////////////
842 //
843 // ReportImmediateHit.
844 //
845 void
846 G4VIntersectionLocator::ReportImmediateHit( co
847 co
848 co
849
850 un
851 {
852 static G4ThreadLocal unsigned int occurredOn
853 static G4ThreadLocal G4ThreeVector* ptrLast
854 if( ptrLast == nullptr )
855 {
856 ptrLast= new G4ThreeVector( DBL_MAX, DBL_
857 G4AutoDelete::Register(ptrLast);
858 }
859 G4ThreeVector &lastStart= *ptrLast;
860
861 if( (TrialPoint - StartPosition).mag2() < to
862 {
863 static G4ThreadLocal unsigned int numUnmo
864 static G4ThreadLocal unsigned int numStil
865
866 G4cout << "Intersection F == start A in "
867 G4cout << "Start Point: " << StartPositio
868 G4cout << " Start-Trial: " << TrialPoint
869 G4cout << " Start-last: " << StartPositio
870
871 if( (StartPosition - lastStart).mag() < t
872 {
873 // We are at position of last 'Start'
874 ++numUnmoved;
875 ++numStill;
876 G4cout << " { Unmoved: " << " still#=
877 << " total # = " << numUnmoved
878 }
879 else
880 {
881 numStill = 0;
882 }
883 G4cout << " Occurred: " << ++occurredOnTo
884 G4cout << " out of total calls= " << num
885 G4cout << G4endl;
886 lastStart = StartPosition;
887 }
888 } // End of ReportImmediateHit()
889