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