| 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 // G4PathFinder Implementation
27 //
28 // Original author: John Apostolakis, April 20
29 // -------------------------------------------
30
31 #include <iomanip>
32
33 #include "G4PathFinder.hh"
34
35 #include "G4SystemOfUnits.hh"
36 #include "G4GeometryTolerance.hh"
37 #include "G4Navigator.hh"
38 #include "G4PropagatorInField.hh"
39 #include "G4TransportationManager.hh"
40 #include "G4MultiNavigator.hh"
41 #include "G4SafetyHelper.hh"
42
43 // Initialise the static instance of the singl
44 //
45 G4ThreadLocal G4PathFinder* G4PathFinder::fpPa
46
47 // -------------------------------------------
48 // GetInstance()
49 //
50 // Retrieve the static instance of the singlet
51 //
52 G4PathFinder* G4PathFinder::GetInstance()
53 {
54 if( fpPathFinder == nullptr )
55 {
56 fpPathFinder = new G4PathFinder;
57 }
58 return fpPathFinder;
59 }
60
61 // -------------------------------------------
62 // GetInstanceIfExist()
63 //
64 // Retrieve the static instance pointer of the
65 //
66 G4PathFinder* G4PathFinder::GetInstanceIfExist
67 {
68 return fpPathFinder;
69 }
70
71 // -------------------------------------------
72 // Constructor
73 //
74 G4PathFinder::G4PathFinder()
75 : fEndState( G4ThreeVector(), G4ThreeVector(
76 {
77 fpMultiNavigator = new G4MultiNavigator();
78
79 fpTransportManager= G4TransportationManager
80 fpFieldPropagator = fpTransportManager->Get
81
82 kCarTolerance = G4GeometryTolerance::GetIns
83
84 G4ThreeVector Big3Vector( kInfinity, kInfi
85 fLastLocatedPosition= Big3Vector;
86 fSafetyLocation= Big3Vector;
87 fPreSafetyLocation= Big3Vector;
88 fPreStepLocation= Big3Vector;
89
90 for( auto num=0; num<fMaxNav; ++num )
91 {
92 fpNavigator[num] = nullptr;
93 fLimitTruth[num] = false;
94 fLimitedStep[num] = kUndefLimited;
95 fCurrentStepSize[num] = -1.0;
96 fLocatedVolume[num] = nullptr;
97 fPreSafetyValues[num]= -1.0;
98 fCurrentPreStepSafety[num] = -1.0;
99 fNewSafetyComputed[num]= -1.0;
100 }
101 }
102
103 // -------------------------------------------
104 // Destructor
105 //
106 G4PathFinder::~G4PathFinder()
107 {
108 delete fpMultiNavigator;
109 fpPathFinder = nullptr;
110 }
111
112 // -------------------------------------------
113 //
114 void
115 G4PathFinder::EnableParallelNavigation(G4bool
116 {
117 G4Navigator *navigatorForPropagation = null
118
119 massNavigator = fpTransportManager->GetNavi
120 if( enableChoice )
121 {
122 navigatorForPropagation = fpMultiNavigat
123
124 // Enable SafetyHelper to use PF
125 //
126 fpTransportManager->GetSafetyHelper()->E
127 }
128 else
129 {
130 navigatorForPropagation = massNavigator;
131
132 // Disable SafetyHelper to use PF
133 //
134 fpTransportManager->GetSafetyHelper()->E
135 }
136 fpFieldPropagator->SetNavigatorForPropagati
137 }
138
139 // -------------------------------------------
140 //
141 G4double
142 G4PathFinder::ComputeStep( const G4FieldTrack&
143 G4double
144 G4int
145 G4int
146 G4double&
147 ELimited&
148 G4FieldTrack&
149 G4VPhysicalVo
150 {
151 G4double possibleStep = -1.0;
152
153 #ifdef G4DEBUG_PATHFINDER
154 if( fVerboseLevel > 2 )
155 {
156 G4cout << " -------------------------" <<
157 G4cout << " G4PathFinder::ComputeStep - en
158 G4cout << " - stepNo = " << std::setw(4
159 << " navigatorId = " << std::setw(2
160 << " proposed step len = " << propo
161 G4cout << " PF::ComputeStep requested step
162 << " from " << InitialFieldTrack.Ge
163 << " dir " << InitialFieldTrack.Ge
164 }
165 #endif
166 #ifdef G4VERBOSE
167 if( navigatorNo >= fNoActiveNavigators )
168 {
169 std::ostringstream message;
170 message << "Bad Navigator ID !" << G4endl
171 << " Requested Navigator ID
172 << " Number of active navig
173 G4Exception("G4PathFinder::ComputeStep()",
174 FatalException, message);
175 }
176 #endif
177
178 if( fNewTrack || (stepNo != fLastStepNo) )
179 {
180 // This is a new track or a new step, so w
181 // ( else we can simply retrieve its resul
182
183 G4FieldTrack currentState = InitialFieldTr
184
185 fCurrentStepNo = stepNo;
186
187 // Check whether a process shifted the pos
188 // since the last step -- by physics proce
189 //
190 G4ThreeVector newPosition = InitialFieldTr
191 G4ThreeVector moveVector = newPosition - f
192 G4double moveLenSq = moveVector.mag2();
193 if( moveLenSq > sqr(kCarTolerance) )
194 {
195 G4ThreeVector newDirection = InitialFie
196 #ifdef G4DEBUG_PATHFINDER
197 if( fVerboseLevel > 2 )
198 {
199 G4double moveLen= std::sqrt( moveLen
200 G4cout << " G4PathFinder::ComputeSte
201 << " -- at step # = " << step
202 << " by " << moveLen << " to
203 }
204 #endif
205 MovePoint(); // Unintentional changed
206
207 // Relocate to cope with this move -- e
208 //
209 Locate( newPosition, newDirection );
210 }
211
212 // Check whether the particle have an (EM)
213 //
214 G4double particleCharge = currentState.Get
215
216 G4FieldManager* fieldMgr = nullptr;
217 G4bool fieldExertsForce = false ;
218 if( particleCharge != 0.0 )
219 {
220 fieldMgr = fpFieldPropagator->FindAndS
221
222 // Protect for case where field manage
223 //
224 fieldExertsForce = (fieldMgr != nullpt
225 && (fieldMgr->GetDetec
226 }
227 fFieldExertedForce = fieldExertsForce; //
228 //
229
230 fNoGeometriesLimiting = -1; // At start o
231 if( fieldExertsForce )
232 {
233 DoNextCurvedStep( currentState, propose
234 //--------------
235 }else{
236 DoNextLinearStep( currentState, propose
237 //--------------
238 }
239 fLastStepNo = stepNo;
240 fRelocatedPoint = false;
241
242 #ifdef G4DEBUG_PATHFINDER
243 if ( (fNoGeometriesLimiting < 0)
244 || (fNoGeometriesLimiting > fNoActiveNav
245 {
246 std::ostringstream message;
247 message << "Number of geometries limitin
248 << " Number of geometries
249 << fNoGeometriesLimiting;
250 G4Exception("G4PathFinder::ComputeStep()
251 "GeomNav0002", FatalExceptio
252 }
253 #endif
254 }
255 #ifdef G4DEBUG_PATHFINDER
256 else
257 {
258 const G4double checkTolerance = 1.0e-9;
259 if( proposedStepLength < fTrueMinStep * (
260 {
261 std::ostringstream message;
262 message.precision( 12 );
263 message << "Problem in step size reques
264 << " Being requested to
265 << " than the minimum Step " <<
266 << " already computed fo
267 << " this tracking-step: " << G
268 << " This could happen d
269 << G4endl
270 << " Check that all phys
271 << " before all processe
272 << G4endl
273 << " If using pre-packag
274 << " functionality, plea
275 << G4endl << G4endl
276 << " Additional informat
277 << " Steps request/propo
278 << G4endl
279 << " MinimumStep (true)
280 << G4endl
281 << " MinimumStep (navraw
282 << G4endl
283 << " Navigator raw retur
284 << " Requested step now
285 << G4endl
286 << " Difference min-req
287 << fTrueMinStep-proposedStepLen
288 << " Relative (to max of
289 << (fTrueMinStep-proposedStepLe
290 / std::max(proposedStepLengt
291 << " -- Step info> stepNo=
292 << " last= " << fLastStepNo
293 << " newTr= " << fNewTrack << G
294 G4Exception("G4PathFinder::ComputeStep
295 "GeomNav0003", FatalExcept
296 }
297 else
298 {
299 // This is neither a new track nor a n
300 // client accessing information for th
301 // We will simply retrieve the results
302 // stepping for this Navigator Id belo
303 //
304 if( fVerboseLevel > 1 )
305 {
306 G4cout << " G4P::CS -> Not calling
307 << " stepNo= " << stepNo <<
308 << " new= " << fNewTrack <<
309 }
310 }
311 }
312 #endif
313
314 fNewTrack = false;
315
316 // Prepare the information to return
317
318 pNewSafety = fCurrentPreStepSafety[ navigat
319 limitedStep = fLimitedStep[ navigatorNo ];
320
321 possibleStep = std::min(proposedStepLength,
322 EndState = fEndState; // now corrected for
323
324 #ifdef G4DEBUG_PATHFINDER
325 if( fVerboseLevel > 0 )
326 {
327 G4cout << " G4PathFinder::ComputeStep retu
328 << fCurrentStepSize[ navigatorNo ]
329 << " for Navigator " << navigatorNo
330 << " Limited step = " << limitedSte
331 << " Safety(mm) = " << pNewSafety /
332 << G4endl;
333 }
334 #endif
335
336 return possibleStep;
337 }
338
339 // -------------------------------------------
340
341 void
342 G4PathFinder::PrepareNewTrack( const G4ThreeVe
343 const G4ThreeVe
344 G4VPhysicalVolu
345 {
346 // Key purposes:
347 // - Check and cache set of active navigat
348 // - Reset state for new track
349
350 G4int num=0;
351
352 EnableParallelNavigation(true);
353 // Switch PropagatorInField to use MultiNa
354
355 fpTransportManager->GetSafetyHelper()->Initi
356 // Reinitialise state of safety helper --
357
358 fNewTrack = true;
359 this->MovePoint(); // Signal further that
360
361 fpFieldPropagator->PrepareNewTrack(); // Inf
362
363 // Message the G4NavigatorPanel / Dispatcher
364 //
365 std::vector<G4Navigator*>::iterator pNavigat
366
367 fNoActiveNavigators = (G4int)fpTransportMana
368 if( fNoActiveNavigators > fMaxNav )
369 {
370 std::ostringstream message;
371 message << "Too many active Navigators / w
372 << " Transportation Manager
373 << fNoActiveNavigators << " active
374 << " This is more than the
375 << fMaxNav << " !";
376 G4Exception("G4PathFinder::PrepareNewTrack
377 FatalException, message);
378 }
379
380 fpMultiNavigator->PrepareNavigators();
381 //------------------------------------
382
383 pNavigatorIter = fpTransportManager->GetActi
384 for( num=0; num< fNoActiveNavigators; ++pNav
385 {
386 // Keep information in C-array ... for cr
387 //
388 fpNavigator[num] = *pNavigatorIter;
389 fLimitTruth[num] = false;
390 fLimitedStep[num] = kDoNot;
391 fCurrentStepSize[num] = 0.0;
392 fLocatedVolume[num] = nullptr;
393 }
394 fNoGeometriesLimiting = 0; // At start of t
395
396 // In case of one geometry, the tracking wil
397
398 if( fNoActiveNavigators > 1 )
399 {
400 Locate( position, direction, false );
401 }
402 else
403 {
404 // Update state -- depending on the track
405
406 fLastLocatedPosition = position;
407 fLocatedVolume[0] = massStartVol; // This
408 // by t
409 fLimitedStep[0] = kDoNot;
410 fCurrentStepSize[0] = 0.0;
411 }
412
413 // Reset Safety Information -- as in case of
414 // inconsistencies ...
415 //
416 fMinSafety_PreStepPt = fPreSafetyMinValue =
417
418 for( num=0; num<fNoActiveNavigators; ++num )
419 {
420 fPreSafetyValues[num] = 0.0;
421 fNewSafetyComputed[num] = 0.0;
422 fCurrentPreStepSafety[num] = 0.0;
423 }
424
425 // The first location for each Navigator mus
426 // or else call ResetStackAndState() for eac
427
428 fRelocatedPoint = false;
429 }
430
431
432 void G4PathFinder::EndTrack()
433 // Signal end of tracking of current track.
434 // Reset TransportationManager to use 'ordin
435 // Reset internal state, if needed
436 {
437 EnableParallelNavigation(false); // Else it
438 }
439
440 void G4PathFinder::ReportMove( const G4ThreeVe
441 const G4ThreeVe
442 const G4String&
443 {
444 G4ThreeVector moveVec = ( NewVector - OldV
445
446 std::ostringstream message;
447 message.precision(16);
448 message << "Endpoint moved between value r
449 << " and call to Locate(). " << G4
450 << " Change of " << Quant
451 << moveVec.mag() / mm << " mm long
452 << " and its vector is "
453 << (1.0/mm) * moveVec << " mm " <<
454 << " Endpoint of ComputeS
455 << G4endl
456 << " and current position
457 G4Exception("G4PathFinder::ReportMove()",
458 JustWarning, message);
459 }
460
461 void G4PathFinder::Locate( const G4ThreeVector
462 const G4ThreeVector
463 G4bool relati
464 {
465 // Locate the point in each geometry
466
467 auto pNavIter = fpTransportManager->GetActiv
468
469 G4ThreeVector lastEndPosition = fRelocatedPo
470 ? fLastLocated
471 : fEndState.Ge
472 fLastLocatedPosition = position;
473
474 #ifdef G4DEBUG_PATHFINDER
475 static const G4double movLenTol = 10*sqr(kCa
476
477 G4ThreeVector moveVec = ( position - lastEnd
478 G4double moveLenSq = moveVec.mag2();
479 if( (!fNewTrack) && ( moveLenSq > movLenTol
480 {
481 ReportMove( lastEndPosition, position,
482 " (End) Position / G4PathFind
483 }
484
485 if( fVerboseLevel > 2 )
486 {
487 G4cout << G4endl;
488 G4cout << " G4PathFinder::Locate : entered
489 G4cout << " -------------------- -------
490 G4cout << " Locating at position " << po
491 << " with direction " << direction
492 << " relative= " << relative << G4
493 if ( (fVerboseLevel > 1) || ( moveLenSq >
494 {
495 G4cout << " lastEndPosition = " << las
496 << " moveVec = " << moveVec
497 << " newTr = " << fNewTrack
498 << " relocated = " << fRelocate
499 }
500
501 G4cout << " Located at " << position ;
502 if( fNoActiveNavigators > 1 ) { G4cout <<
503 }
504 #endif
505
506 for ( auto num=0; num<fNoActiveNavigators ;
507 {
508 // ... who limited the step ....
509
510 if( fLimitTruth[num] ) { (*pNavIter)->Set
511
512 G4VPhysicalVolume *pLocated=
513 (*pNavIter)->LocateGlobalPointAndSetup( p
514 r
515 f
516 // Set the state related to the location
517 //
518 fLocatedVolume[num] = pLocated;
519
520 // Clear state related to the step
521 //
522 fLimitedStep[num] = kDoNot;
523 fCurrentStepSize[num] = 0.0;
524
525 #ifdef G4DEBUG_PATHFINDER
526 if( fVerboseLevel > 2 )
527 {
528 G4cout << " - In world " << num << " ge
529 << " gives volume= " << pLocate
530 if( pLocated )
531 {
532 G4cout << " name = '" << pLocated->G
533 G4cout << " - CopyNo= " << pLocated->
534 }
535 G4cout << G4endl;
536 }
537 #endif
538 }
539
540 fRelocatedPoint = false;
541 }
542
543 void G4PathFinder::ReLocate( const G4ThreeVect
544 {
545 // Locate the point in each geometry
546
547 auto pNavIter = fpTransportManager->GetActiv
548
549 #ifdef G4DEBUG_PATHFINDER
550
551 // Check that this relocation does not viola
552 // - at endpoint (computed from start poin
553 // - at last safety location (likely just
554
555 G4ThreeVector lastEndPosition = fEndState.Ge
556
557 // Calculate end-point safety ...
558 //
559 G4double DistanceStartEnd = (lastEndPosition
560 G4double endPointSafety_raw = fMinSafety_Pre
561 G4double endPointSafety_Est1 = std::max( 0.0
562
563 // ... and check move from endpoint against
564 //
565 G4ThreeVector moveVecEndPos = position - la
566 G4double moveLenEndPosSq = moveVecEndPo
567
568 // Check that move from endpoint of last ste
569 // -- or check against last location or relo
570 //
571 G4ThreeVector moveVecSafety = position - fSa
572 G4double moveLenSafSq = moveVecSafety.m
573
574 G4double distCheckEnd_sq = ( moveLenEndPosSq
575
576 G4double distCheckSaf_sq = ( moveLenSafSq -
577 *
578
579 G4bool longMoveEnd = distCheckEnd_sq > 0.0;
580 G4bool longMoveSaf = distCheckSaf_sq > 0.0;
581
582 G4double revisedSafety = 0.0;
583
584 if( (!fNewTrack) && ( longMoveEnd && longMov
585 {
586 // Recompute ComputeSafety for end positi
587 //
588 revisedSafety = ComputeSafety(lastEndPosi
589
590 const G4double kRadTolerance =
591 G4GeometryTolerance::GetInstance()-
592 const G4double cErrorTolerance = 1e-12;
593 // Maximum relative error from roundoff
594
595 G4double distCheckRevisedEnd = moveLenEnd
596
597 G4bool longMoveRevisedEnd = ( distCheckR
598
599 G4double moveMinusSafety = 0.0;
600 G4double moveLenEndPosition = std::sqrt(
601 moveMinusSafety = moveLenEndPosition - re
602
603 if ( longMoveRevisedEnd && ( moveMinusSaf
604 && ( revisedSafet
605 {
606 // Take into account possibility of ro
607 // this apparent move further than saf
608
609 if( fVerboseLevel > 0 )
610 {
611 G4cout << " G4PF:Relocate> Ratio to
612 << std::fabs(moveMinusSafety
613 }
614
615 G4double absMoveMinusSafety = std::fab
616 G4bool smallRatio = absMoveMinusSafety
617 G4double maxCoordPos = std::max(
618 std::max
619
620 std::fab
621 G4bool smallValue= absMoveMinusSafety
622 if( !(smallRatio || smallValue) )
623 {
624 G4cout << " G4PF:Relocate> Ratio to
625 << std::fabs(moveMinusSafety
626 G4cout << " Difference of move and
627 << G4endl;
628 }
629 else
630 {
631 moveMinusSafety = 0.0;
632 longMoveRevisedEnd = false; // Num
633
634 G4cout << " Difference of move & saf
635 << absMoveMinusSafety << G4en
636 if( smallRatio )
637 {
638 G4cout << " ratio to safety " << r
639 << " is " << absMoveMinusS
640 << "smaller than " << kRadT
641 }
642 else
643 {
644 G4cout << " as fraction " << absMo
645 << " of position vector max
646 << " smaller than " << cErr
647 }
648 G4cout << " -- reset moveMinusSafety
649 << moveMinusSafety << G4endl;
650 }
651 }
652
653 if ( longMoveEnd && longMoveSaf
654 && longMoveRevisedEnd && (moveMinusSafe
655 {
656 std::ostringstream message;
657 message.precision(9);
658 message << "ReLocation is further than
659 << " Moved from last endpoint
660 << " compared to end safety (f
661 << endPointSafety_Est1 << G4en
662 << " --> last PreSafety Locat
663 << G4endl
664 << " safety value = " <
665 << " --> last PreStep Locatio
666 << G4endl
667 << " safety value = " <
668 << " --> last EndStep Locatio
669 << G4endl
670 << " safety value = " <
671 << " raw-value = " << endPoint
672 << " --> Calling again at thi
673 << revisedSafety << " as safe
674 << " --> last position for sa
675 << G4endl
676 << " its safety value =
677 << G4endl
678 << " move from safety lo
679 << std::sqrt(moveLenSafSq) <<
680 << " again= " << moveV
681 << " safety - Move-from-
682 << revisedSafety - moveLenEndP
683 << " (negative is Bad.)" << G4
684 << " Debug: distCheckRevisedE
685 << distCheckRevisedEnd;
686 ReportMove( lastEndPosition, position,
687 G4Exception("G4PathFinder::ReLocate",
688 FatalException, message);
689 }
690 }
691
692 if( fVerboseLevel > 2 )
693 {
694 G4cout << G4endl;
695 G4cout << " G4PathFinder::ReLocate : enter
696 G4cout << " ---------------------- -----
697 G4cout << " *Re*Locating at position " <<
698 if ( (fVerboseLevel > -1) || ( moveLenEndP
699 {
700 G4cout << " lastEndPosition = " << las
701 << " moveVec from step-end = "
702 << " is new Track = " << fNewTr
703 << " relocated = " << fRelocate
704 }
705 }
706 #endif // G4DEBUG_PATHFINDER
707
708 for ( auto num=0; num< fNoActiveNavigators ;
709 {
710 // ... none limited the step
711
712 (*pNavIter)->LocateGlobalPointWithinVolum
713
714 // Clear state related to the step
715 //
716 fLimitedStep[num] = kDoNot;
717 fCurrentStepSize[num] = 0.0;
718 fLimitTruth[num] = false;
719 }
720
721 fLastLocatedPosition = position;
722 fRelocatedPoint = true;
723
724 #ifdef G4DEBUG_PATHFINDER
725 if( fVerboseLevel > 2 )
726 {
727 G4cout << " G4PathFinder::ReLocate : exiti
728 << " at position " << fLastLocated
729 }
730 #endif
731 }
732
733 // -------------------------------------------
734
735 G4double G4PathFinder::ComputeSafety( const G
736 {
737 // Recompute safety for the relevant point
738
739 G4double minSafety = kInfinity;
740
741 std::vector<G4Navigator*>::iterator pNaviga
742 pNavigatorIter = fpTransportManager->GetAct
743
744 for( auto num=0; num<fNoActiveNavigators; +
745 {
746 G4double safety = (*pNavigatorIter)->Com
747 if( safety < minSafety ) { minSafety = s
748 fNewSafetyComputed[num] = safety;
749 }
750
751 fSafetyLocation = position;
752 fMinSafety_atSafLocation = minSafety;
753
754 #ifdef G4DEBUG_PATHFINDER
755 if( fVerboseLevel > 1 )
756 {
757 G4cout << " G4PathFinder::ComputeSafety -
758 << minSafety << " at location " <<
759 }
760 #endif
761 return minSafety;
762 }
763
764
765 // -------------------------------------------
766
767 G4TouchableHandle
768 G4PathFinder::CreateTouchableHandle( G4int nav
769 {
770 #ifdef G4DEBUG_PATHFINDER
771 if( fVerboseLevel > 2 )
772 {
773 G4cout << "G4PathFinder::CreateTouchableHa
774 << navId << " -- " << GetNavigator(
775 }
776 #endif
777
778 G4TouchableHistory* touchHist;
779 touchHist = GetNavigator(navId)->CreateTouch
780
781 G4VPhysicalVolume* locatedVolume = fLocatedV
782 if( locatedVolume == nullptr )
783 {
784 // Workaround to ensure that the touchabl
785
786 touchHist->UpdateYourself( locatedVolume,
787 }
788
789 #ifdef G4DEBUG_PATHFINDER
790 if( fVerboseLevel > 2 )
791 {
792 G4String VolumeName("None");
793 if( locatedVolume ) { VolumeName = located
794 G4cout << " Touchable History created at a
795 << "; volume = " << locatedVolume <
796 << G4endl;
797 }
798 #endif
799
800 return {touchHist};
801 }
802
803 G4double
804 G4PathFinder::DoNextLinearStep( const G4FieldT
805 G4double
806 {
807 std::vector<G4Navigator*>::iterator pNavigat
808 G4double safety = 0.0, step =0.0;
809 G4double minSafety = kInfinity, minStep = kI
810
811 const G4int IdTransport = 0; // Id of Mass
812 G4int num = 0;
813
814 #ifdef G4DEBUG_PATHFINDER
815 if( fVerboseLevel > 2 )
816 {
817 G4cout << " G4PathFinder::DoNextLinearStep
818 G4cout << " Input field track= " << init
819 G4cout << " Requested step= " << propose
820 }
821 #endif
822
823 G4ThreeVector initialPosition = initialState
824 G4ThreeVector initialDirection = initialStat
825
826 G4ThreeVector OriginShift = initialPosition
827 G4double MagSqShift = OriginShift.mag2
828 G4double MagShift; // Only given value
829
830 // Potential optimisation using Maximum Valu
831 // if( MagSqShift >= sqr(fPreSafetyMaxValue
832 // MagShift= kInfinity; // Not a useful
833 // else
834 // MagShift= std::sqrt(MagSqShift) ;
835
836 MagShift= std::sqrt(MagSqShift) ;
837
838 #ifdef G4PATHFINDER_OPTIMISATION
839
840 G4double fullSafety; // For all geometries,
841
842 if( MagSqShift >= sqr(fPreSafetyMinValue) )
843 {
844 fullSafety = 0.0 ;
845 }
846 else
847 {
848 fullSafety = fPreSafetyMinValue - MagShif
849 }
850 if( proposedStepLength < fullSafety )
851 {
852 // Move is smaller than all safeties
853 // -> so we do not have to move the safe
854
855 fPreStepCenterRenewed = false;
856
857 for( num=0; num< fNoActiveNavigators; ++n
858 {
859 fCurrentStepSize[num] = kInfinity;
860 safety = std::max( 0.0, fPreSafetyVal
861 minSafety= std::min( safety, minSafety
862 fCurrentPreStepSafety[num] = safety;
863 }
864 minStep = kInfinity;
865
866 #ifdef G4DEBUG_PATHFINDER
867 if( fVerboseLevel > 2 )
868 {
869 G4cout << "G4PathFinder::DoNextLinearSt
870 << " proposedStepLength " << p
871 << " < (full) safety = " << ful
872 << " at " << initialPosition
873 << G4endl;
874 }
875 #endif
876 }
877 else
878 #endif // End of G4PATHFINDER_OPTIMISATION 1
879 {
880 // Move is larger than at least one of th
881 // -> so we must move the safety center!
882
883 fPreStepCenterRenewed = true;
884 pNavigatorIter = fpTransportManager-> Get
885
886 minStep = kInfinity; // Not proposedStep
887
888 for( num=0; num< fNoActiveNavigators; ++p
889 {
890 safety = std::max( 0.0, fPreSafetyVal
891
892 #ifdef G4PATHFINDER_OPTIMISATION
893 if( proposedStepLength <= safety ) //
894 {
895 // The Step is guaranteed to be tak
896
897 step = kInfinity; // ComputeSte
898
899 #ifdef G4DEBUG_PATHFINDER
900 G4cout.precision(8);
901 G4cout << "PathFinder::ComputeStep>
902 << proposedStepLength
903 << " <= safety = " << safet
904 << " Step fully taken. " <<
905 #endif
906 }
907 else
908 #endif // End of G4PATHFINDER_OPTIMISATION 2
909 {
910 #ifdef G4DEBUG_PATHFINDER
911 G4double previousSafety = safety;
912 #endif
913 step = (*pNavigatorIter)->ComputeSt
914
915
916
917 minStep = std::min(step, minStep);
918
919
920 #ifdef G4DEBUG_PATHFINDER
921 if( fVerboseLevel > 0)
922 {
923 G4cout.precision(8);
924 G4cout << "PathFinder::ComputeSte
925 << proposedStepLength
926 << " > safety = " << pre
927 << " for nav " << num
928 << " . New safety = " <<
929 << G4endl;
930 }
931 #endif
932 }
933 fCurrentStepSize[num] = step; // Raw
934
935 // TODO: consider whether/how to r
936 // to the latest minStep val
937 // ( If so, much change 1st
938
939 // Save safety value, must be done for
940 // (even if not recomputed using call
941 // since they share the fPreSafetyLoca
942
943 fPreSafetyValues[num] = safety;
944 fCurrentPreStepSafety[num] = safety;
945
946 minSafety = std::min( safety, minSafet
947
948 #ifdef G4DEBUG_PATHFINDER
949 if( fVerboseLevel > 2 )
950 {
951 G4cout << "G4PathFinder::DoNextLinea
952 << num << "] -- step size " <
953 }
954 #endif
955 }
956
957 // Only change these when safety is recal
958 // it is good/relevant only for safety ca
959
960 fPreSafetyLocation = initialPosition;
961 fPreSafetyMinValue = minSafety;
962 } // end of else for if( proposedStepLength
963
964 // For use in Relocation, need PreStep point
965 //
966 fPreStepLocation = initialPosition;
967 fMinSafety_PreStepPt = minSafety;
968
969 fMinStep = minStep;
970
971 if( fMinStep == kInfinity )
972 {
973 minStep = proposedStepLength; // Use t
974 }
975 fTrueMinStep = minStep;
976
977 // Set the EndState
978
979 G4ThreeVector endPosition;
980
981 fEndState = initialState;
982 endPosition = initialPosition + minStep * in
983
984 #ifdef G4DEBUG_PATHFINDER
985 if( fVerboseLevel > 1 )
986 {
987 G4int oldPrec= G4cout.precision(14);
988 G4cout << "G4PathFinder::DoNextLinearStep
989 << " initialPosition = " << initial
990 << " and endPosition = " << endPosi
991 G4cout.precision(oldPrec);
992 }
993 #endif
994
995 fEndState.SetPosition( endPosition );
996 fEndState.SetProperTimeOfFlight( -1.000 );
997
998 if( fNoActiveNavigators == 1 )
999 {
1000 G4bool transportLimited = (fMinStep!= kI
1001 fLimitTruth[IdTransport] = transportLimi
1002 fLimitedStep[IdTransport] = transportLim
1003
1004 // Set fNoGeometriesLimiting - as WhichL
1005 fNoGeometriesLimiting = transportLimited
1006 }
1007 else
1008 {
1009 WhichLimited();
1010 }
1011
1012 #ifdef G4DEBUG_PATHFINDER
1013 if( fVerboseLevel > 2 )
1014 {
1015 G4cout << " G4PathFinder::DoNextLinearSte
1016 << minStep << G4endl;
1017 G4cout << " - Endpoint values = " << fEnd
1018 G4cout << G4endl;
1019 }
1020 #endif
1021
1022 return minStep;
1023 }
1024
1025 void G4PathFinder::WhichLimited()
1026 {
1027 // Flag which processes limited the step
1028
1029 G4int num = -1, last = -1;
1030 G4int noLimited = 0;
1031 ELimited shared = kSharedOther;
1032
1033 const G4int IdTransport = 0; // Id of Mass
1034
1035 // Assume that [IdTransport] is Mass / Tran
1036 //
1037 G4bool transportLimited = (fCurrentStepSize
1038 && (fMinStep != kInf
1039
1040 if( transportLimited )
1041 {
1042 shared= kSharedTransport;
1043 }
1044
1045 for ( num = 0; num < fNoActiveNavigators; +
1046 {
1047 G4bool limitedStep;
1048
1049 G4double step = fCurrentStepSize[num];
1050
1051 limitedStep = ( std::fabs(step - fMinStep
1052 && ( step != kInfinity);
1053
1054 fLimitTruth[ num ] = limitedStep;
1055 if( limitedStep )
1056 {
1057 ++noLimited;
1058 fLimitedStep[num] = shared;
1059 last= num;
1060 }
1061 else
1062 {
1063 fLimitedStep[num] = kDoNot;
1064 }
1065 }
1066 fNoGeometriesLimiting= noLimited; // Save
1067
1068 if( (last > -1) && (noLimited == 1 ) )
1069 {
1070 fLimitedStep[ last ] = kUnique;
1071 }
1072
1073 #ifdef G4DEBUG_PATHFINDER
1074 if( fVerboseLevel > 1 )
1075 {
1076 PrintLimited(); // --> for tracing
1077 if( fVerboseLevel > 4 )
1078 {
1079 G4cout << " G4PathFinder::WhichLimited
1080 }
1081 }
1082 #endif
1083 }
1084
1085 void G4PathFinder::PrintLimited()
1086 {
1087 // Report results -- for checking
1088
1089 G4cout << "G4PathFinder::PrintLimited repor
1090 G4cout << " Minimum step (true)= " << fTru
1091 << " reported min = " << fMinStep
1092 << G4endl;
1093 if( (fCurrentStepNo <= 2) || (fVerboseLeve
1094 {
1095 G4cout << std::setw(5) << " Step#" << "
1096 << std::setw(5) << " NavId" << "
1097 << std::setw(12) << " step-size "
1098 << std::setw(12) << " raw-size "
1099 << std::setw(12) << " pre-safety "
1100 << std::setw(15) << " Limited / fl
1101 << std::setw(15) << " World " <<
1102 << G4endl;
1103 }
1104 for ( auto num = 0; num < fNoActiveNavigato
1105 {
1106 G4double rawStep = fCurrentStepSize[num];
1107 G4double stepLen = fCurrentStepSize[num];
1108 if( stepLen > fTrueMinStep )
1109 {
1110 stepLen = fTrueMinStep; // did not
1111 }
1112 G4long oldPrec = G4cout.precision(9);
1113
1114 G4cout << std::setw(5) << fCurrentStepNo
1115 << std::setw(5) << num << " "
1116 << std::setw(12) << stepLen << " "
1117 << std::setw(12) << rawStep << " "
1118 << std::setw(12) << fCurrentPreSte
1119 << std::setw(5) << (fLimitTruth[nu
1120 G4String limitedStr= LimitedString(fLimit
1121 G4cout << " " << std::setw(15) << limited
1122 G4cout.precision(oldPrec);
1123
1124 G4Navigator* pNav = GetNavigator( num );
1125 G4String WorldName( "Not-Set" );
1126 if (pNav != nullptr)
1127 {
1128 G4VPhysicalVolume *pWorld = pNav->GetW
1129 if( pWorld != nullptr )
1130 {
1131 WorldName = pWorld->GetName();
1132 }
1133 }
1134 G4cout << " " << WorldName ;
1135 G4cout << G4endl;
1136 }
1137
1138 if( fVerboseLevel > 4 )
1139 {
1140 G4cout << " G4PathFinder::PrintLimited -
1141 }
1142 }
1143
1144 G4double
1145 G4PathFinder::DoNextCurvedStep( const G4Field
1146 G4doubl
1147 G4VPhysicalVo
1148 {
1149 const G4double toleratedRelativeError = 1.0
1150 G4double minStep= kInfinity, newSafety = 0.
1151 G4int numNav;
1152 G4FieldTrack fieldTrack = initialState;
1153 G4ThreeVector startPoint = initialState.Get
1154
1155
1156 G4EquationOfMotion* equationOfMotion =
1157 fpFieldPropagator->GetCurrentEquationOfM
1158
1159 equationOfMotion->SetChargeMomentumMass( *(
1160 in
1161 in
1162
1163 #ifdef G4DEBUG_PATHFINDER
1164 G4int prc = G4cout.precision(9);
1165 if( fVerboseLevel > 2 )
1166 {
1167 G4cout << " G4PathFinder::DoNextCurvedSte
1168 G4cout << " Initial value of field track
1169 << " and proposed step= " << propo
1170 }
1171 #endif
1172
1173 fPreStepCenterRenewed = true; // Always upd
1174
1175 if( fNoActiveNavigators > 1 )
1176 {
1177 // Calculate the safety values before ma
1178
1179 G4double minSafety= kInfinity, safety;
1180 for( numNav=0; numNav < fNoActiveNavigat
1181 {
1182 safety= fpNavigator[numNav]->ComputeS
1183 fPreSafetyValues[numNav] = safety;
1184 fCurrentPreStepSafety[numNav] = safet
1185 minSafety = std::min( safety, minSafe
1186 }
1187
1188 // Save safety value, related position
1189
1190 fPreSafetyLocation = startPoint;
1191 fPreSafetyMinValue = minSafety;
1192 fPreStepLocation = startPoint;
1193 fMinSafety_PreStepPt = minSafety;
1194 }
1195
1196 // Allow Propagator In Field to do the hard
1197 //
1198 minStep = fpFieldPropagator->ComputeStep( f
1199 p
1200 n
1201 p
1202 f
1203
1204 // fieldTrack now contains the endpoint inf
1205 //
1206 fEndState = fieldTrack;
1207 fMinStep = minStep;
1208 fTrueMinStep = std::min( minStep, proposedS
1209
1210 if( fNoActiveNavigators == 1 )
1211 {
1212 // Update the 'PreSafety' sphere - as an
1213 // (must be done anyway in field)
1214
1215 fPreSafetyValues[0] = newSafety;
1216 fPreSafetyLocation = startPoint;
1217 fPreSafetyMinValue = newSafety;
1218
1219 // Update the current 'PreStep' point's
1220 //
1221 fCurrentPreStepSafety[0] = newSafety;
1222 fPreStepLocation = startPoint;
1223 fMinSafety_PreStepPt= newSafety;
1224 }
1225
1226 #ifdef G4DEBUG_PATHFINDER
1227 if( fVerboseLevel > 2 )
1228 {
1229 G4cout << "G4PathFinder::DoNextCurvedStep
1230 << " initialState = " << initialSt
1231 << " and endState = " << fEndState
1232 G4cout << "G4PathFinder::DoNextCurvedStep
1233 << " minStep = " << minStep
1234 << " proposedStepLength " << propo
1235 << " safety = " << newSafety << G4
1236 }
1237 #endif
1238 G4double currentStepSize; // = 0.0;
1239 if( minStep < proposedStepLength ) // if ==
1240 {
1241 // Recover the remaining information from
1242 // especially regarding which Navigator l
1243
1244 G4int noLimited = 0; // No geometries
1245 for( numNav=0; numNav < fNoActiveNavigato
1246 {
1247 G4double finalStep, lastPreSafety = 0.0
1248 ELimited didLimit;
1249 G4bool limited;
1250
1251 finalStep= fpMultiNavigator->ObtainFin
1252
1253
1254 // Calculate the step for this geometry
1255 // final step (the only one which can d
1256
1257 currentStepSize = fTrueMinStep;
1258 G4double diffStep = 0.0;
1259 if( (minStepLast != kInfinity) )
1260 {
1261 diffStep = (finalStep-minStepLast);
1262 if ( std::abs(diffStep) <= toleratedR
1263 {
1264 diffStep = 0.0;
1265 }
1266 currentStepSize += diffStep;
1267 }
1268 fCurrentStepSize[numNav] = currentStepS
1269
1270 // TODO: could refine the way to obtain
1271 // - for pre step safety
1272 // notify MultiNavigator about n
1273 // allow it to return this value
1274 // instead of lastPreSafety (or
1275 // - for final step start (availabl
1276 // get final Step start from Mul
1277 // and corresponding safety valu
1278 // and/or ALSO calculate ComputeSafety
1279 // endSafety= fpNavigator[numNav]->
1280
1281 fLimitedStep[numNav] = didLimit;
1282 fLimitTruth[numNav] = limited = (didLim
1283 if( limited ) { ++noLimited; }
1284
1285 #ifdef G4DEBUG_PATHFINDER
1286 G4bool StepError = (currentStepSize < 0
1287 || ( (minStepLast != kI
1288 if( StepError || (fVerboseLevel > 2) )
1289 {
1290 G4String limitedString = LimitedStr
1291
1292 G4cout << " G4PathFinder::ComputeStep
1293 << " step= " << fCurrentStepS
1294 << " from final-step= " << fin
1295 << " fTrueMinStep= " << fTrueM
1296 << " minStepLast= " << minSte
1297 << " limited = " << (fLimitTr
1298 << " ";
1299 G4cout << " status = " << limitedStr
1300 << G4endl;
1301
1302 if( StepError )
1303 {
1304 std::ostringstream message;
1305 message << "Incorrect calculation o
1306 << G4endl
1307 << " currentStepSize
1308 << ", diffStep= " << diffSt
1309 << "ERROR in computing step
1310 G4Exception("G4PathFinder::DoNextCu
1311 "GeomNav0003", FatalExc
1312 }
1313 }
1314 #endif
1315 } // for num Navigators
1316
1317 fNoGeometriesLimiting = noLimited; // Sa
1318 }
1319 else if ( (minStep == proposedStepLength)
1320 || (minStep == kInfinity)
1321 || ( std::abs(minStep-proposedSte
1322 < toleratedRelativeError * pro
1323 {
1324 // In case the step was not limited, use
1325 // --> all Navigators
1326 // Also avoid problems in case of PathFin
1327 // - it is possible that the Navigators
1328 // if the safety was already satisfact
1329 // (In that case calling ObtainFinalSt
1330
1331 currentStepSize = minStep;
1332 for( numNav=0; numNav < fNoActiveNavigato
1333 {
1334 fCurrentStepSize[numNav] = minStep;
1335 // Safety for endpoint ?? // Can event
1336 fLimitedStep[numNav] = kDoNot;
1337 fLimitTruth[numNav] = false;
1338 }
1339 fNoGeometriesLimiting = 0; // Save # pro
1340 }
1341 else // (minStep > proposedStepLength)
1342 {
1343 std::ostringstream message;
1344 message << "Incorrect calculation of step
1345 << " currentStepSize = " <
1346 << " proposed StepSize = " << pro
1347 G4Exception("G4PathFinder::DoNextCurvedSt
1348 "GeomNav0003", FatalException
1349 }
1350
1351 #ifdef G4DEBUG_PATHFINDER
1352 if( fVerboseLevel > 2 )
1353 {
1354 G4cout << " Exiting G4PathFinder::DoNextC
1355 PrintLimited();
1356 }
1357 G4cout.precision(prc);
1358 #endif
1359
1360 return minStep;
1361 }
1362
1363 G4String& G4PathFinder::LimitedString( ELimit
1364 {
1365 static G4String StrDoNot("DoNot"),
1366 StrUnique("Unique"),
1367 StrUndefined("Undefined"),
1368 StrSharedTransport("SharedT
1369 StrSharedOther("SharedOther
1370
1371 G4String* limitedStr;
1372 switch ( lim )
1373 {
1374 case kDoNot: limitedStr = &StrDoNot; br
1375 case kUnique: limitedStr = &StrUnique; b
1376 case kSharedTransport: limitedStr = &St
1377 case kSharedOther: limitedStr = &StrShar
1378 default: limitedStr = &StrUndefined; bre
1379 }
1380 return *limitedStr;
1381 }
1382
1383 void G4PathFinder::PushPostSafetyToPreSafety(
1384 {
1385 fPreSafetyLocation = fSafetyLocation;
1386 fPreSafetyMinValue = fMinSafety_atSafLocati
1387 for( auto nav=0; nav < fNoActiveNavigators;
1388 {
1389 fPreSafetyValues[nav] = fNewSafetyComput
1390 }
1391 }
1392