| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // class G4ITNavigator2 Implementation
28 //
29 // Original author: Paul Kent, July 95/96
30 //
31 // G4ITNavigator2 is a duplicate version of G4Navigator starting from Geant4.9.5
32 // initially written by Paul Kent and colleagues.
33 // The only difference resides in the way the information is saved and managed
34 //
35 // History:
36 // - Created. Paul Kent, Jul 95/96
37 // - Zero step protections J.A. / G.C., Nov 2004
38 // - Added check mode G. Cosmo, Mar 2004
39 // - Made Navigator Abstract G. Cosmo, Nov 2003
40 // - G4ITNavigator2 created M.K., Nov 2012
41 // --------------------------------------------------------------------
42
43 #include <iomanip>
44
45 //#include "G4ITNavigator2.hh"
46 #include "G4ITNavigator.hh"
47 #include "G4ios.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4GeometryTolerance.hh"
50 #include "G4VPhysicalVolume.hh"
51
52 //#define G4DEBUG_NAVIGATION 1
53 #include "G4VoxelSafety.hh"
54 #include "Randomize.hh"
55 #include "G4VoxelLimits.hh"
56
57 #define fHistory fpNavigatorState->fHistory
58 #define fLastTriedStepComputation fpNavigatorState->fLastTriedStepComputation
59 #define fEnteredDaughter fpNavigatorState->fEnteredDaughter
60 #define fExitedMother fpNavigatorState->fExitedMother
61 #define fWasLimitedByGeometry fpNavigatorState->fWasLimitedByGeometry
62 #define fStepEndPoint fpNavigatorState->fStepEndPoint
63 #define fLastStepEndPointLocal fpNavigatorState->fLastStepEndPointLocal
64 #define fPushed fpNavigatorState->fPushed
65 #define fLastTriedStepComputation fpNavigatorState->fLastTriedStepComputation
66 #define fEntering fpNavigatorState->fEntering
67 #define fExiting fpNavigatorState->fExiting
68 #define fBlockedPhysicalVolume fpNavigatorState->fBlockedPhysicalVolume
69 #define fBlockedReplicaNo fpNavigatorState->fBlockedReplicaNo
70 #define fLastLocatedPointLocal fpNavigatorState->fLastLocatedPointLocal
71 #define fLocatedOutsideWorld fpNavigatorState->fLocatedOutsideWorld
72 #define fValidExitNormal fpNavigatorState->fValidExitNormal
73 #define fExitNormal fpNavigatorState->fExitNormal
74 #define fGrandMotherExitNormal fpNavigatorState->fGrandMotherExitNormal
75 #define fChangedGrandMotherRefFrame fpNavigatorState->fChangedGrandMotherRefFrame
76 #define fExitNormalGlobalFrame fpNavigatorState->fExitNormalGlobalFrame
77 #define fCalculatedExitNormal fpNavigatorState->fCalculatedExitNormal
78 #define fLastStepWasZero fpNavigatorState->fLastStepWasZero
79 #define fLocatedOnEdge fpNavigatorState->fLocatedOnEdge
80 #define fNumberZeroSteps fpNavigatorState->fNumberZeroSteps
81 #define fPreviousSftOrigin fpNavigatorState->fPreviousSftOrigin
82 #define fPreviousSafety fpNavigatorState->fPreviousSafety
83
84 //------------------------------------------------------------------------------
85
86 namespace BoundingBox
87 {
88 enum Boundary
89 {
90 kMin,
91 kMax
92 };
93 }
94
95 // ********************************************************************
96 // Constructor
97 // ********************************************************************
98 //
99 G4ITNavigator2::G4ITNavigator2()
100 {
101 fActive= false;
102 fActionThreshold_NoZeroSteps = 1000;
103 fAbandonThreshold_NoZeroSteps = 2500;
104 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
105 fregularNav.SetNormalNavigation( &fnormalNav );
106
107 fpNavigatorState = nullptr;
108 // fSaveState = 0;
109 fpVoxelSafety= new G4VoxelSafety();
110
111 // this->SetVerboseLevel(3);
112 // this->CheckMode(true);
113 }
114
115 // ********************************************************************
116 // Destructor
117 // ********************************************************************
118 //
119 G4ITNavigator2::~G4ITNavigator2()
120 { delete fpVoxelSafety; }
121
122 // ********************************************************************
123 // ResetHierarchyAndLocate
124 // ********************************************************************
125 //
126 G4VPhysicalVolume*
127 G4ITNavigator2::ResetHierarchyAndLocate(const G4ThreeVector &p,
128 const G4ThreeVector &direction,
129 const G4TouchableHistory &h)
130 {
131 // ResetState();
132 fHistory = *h.GetHistory();
133 SetupHierarchy();
134 fLastTriedStepComputation= false; // Redundant, but best
135 return LocateGlobalPointAndSetup(p, &direction, true, false);
136 }
137
138 // ********************************************************************
139 // LocateGlobalPointAndSetup
140 //
141 // Locate the point in the hierarchy return 0 if outside
142 // The direction is required
143 // - if on an edge shared by more than two surfaces
144 // (to resolve likely looping in tracking)
145 // - at initial location of a particle
146 // (to resolve potential ambiguity at boundary)
147 //
148 // Flags on exit: (comments to be completed)
149 // fEntering - True if entering `daughter' volume (or replica)
150 // whether daughter of last mother directly
151 // or daughter of that volume's ancestor.
152 // ********************************************************************
153 //
154 G4VPhysicalVolume*
155 G4ITNavigator2::LocateGlobalPointAndSetup( const G4ThreeVector& globalPoint,
156 const G4ThreeVector* pGlobalDirection,
157 const G4bool relativeSearch,
158 const G4bool ignoreDirection )
159 {
160 CheckNavigatorStateIsValid();
161 G4bool notKnownContained=true, noResult;
162 G4VPhysicalVolume *targetPhysical;
163 G4LogicalVolume *targetLogical;
164 G4VSolid *targetSolid=nullptr;
165 G4ThreeVector localPoint, globalDirection;
166 EInside insideCode;
167
168 G4bool considerDirection = (!ignoreDirection) || fLocatedOnEdge;
169
170 fLastTriedStepComputation= false;
171 fChangedGrandMotherRefFrame= false; // For local exit normal
172
173 if( considerDirection && pGlobalDirection != nullptr )
174 {
175 globalDirection=*pGlobalDirection;
176 }
177
178
179 #ifdef G4VERBOSE
180 if( fVerbose > 2 )
181 {
182 G4long oldcoutPrec = G4cout.precision(8);
183 G4cout << "*** G4ITNavigator2::LocateGlobalPointAndSetup: ***" << G4endl;
184 G4cout << " Called with arguments: " << G4endl
185 << " Globalpoint = " << globalPoint << G4endl
186 << " RelativeSearch = " << relativeSearch << G4endl;
187 if( fVerbose == 4 )
188 {
189 G4cout << " ----- Upon entering:" << G4endl;
190 PrintState();
191 }
192 G4cout.precision(oldcoutPrec);
193 }
194 #endif
195
196 G4int noLevelsExited=0 ;
197
198 if ( !relativeSearch )
199 {
200 ResetStackAndState();
201 }
202 else
203 {
204 if ( fWasLimitedByGeometry )
205 {
206 fWasLimitedByGeometry = false;
207 fEnteredDaughter = fEntering; // Remember
208 fExitedMother = fExiting; // Remember
209 if ( fExiting )
210 {
211 noLevelsExited++; // count this first level entered too
212
213 if ( fHistory.GetDepth() )
214 {
215 fBlockedPhysicalVolume = fHistory.GetTopVolume();
216 fBlockedReplicaNo = fHistory.GetTopReplicaNo();
217 fHistory.BackLevel();
218 }
219 else
220 {
221 fLastLocatedPointLocal = localPoint;
222 fLocatedOutsideWorld = true;
223 fBlockedPhysicalVolume = nullptr; // to be sure
224 fBlockedReplicaNo = -1;
225 fEntering = false; // No longer
226 fEnteredDaughter = false;
227 fExitedMother = true; // ??
228
229 return nullptr; // Have exited world volume
230 }
231 // A fix for the case where a volume is "entered" at an edge
232 // and a coincident surface exists outside it.
233 // - This stops it from exiting further volumes and cycling
234 // - However ReplicaNavigator treats this case itself
235 //
236 // assert( fBlockedPhysicalVolume!=0 );
237
238 // Expect to be on edge => on surface
239 //
240 if ( fLocatedOnEdge && (VolumeType(fBlockedPhysicalVolume)!=kReplica ))
241 {
242 fExiting= false;
243 // Consider effect on Exit Normal !?
244 }
245 }
246 else
247 if ( fEntering )
248 {
249 // assert( fBlockedPhysicalVolume!=0 );
250
251 switch (VolumeType(fBlockedPhysicalVolume))
252 {
253 case kNormal:
254 fHistory.NewLevel(fBlockedPhysicalVolume, kNormal,
255 fBlockedPhysicalVolume->GetCopyNo());
256 break;
257 case kReplica:
258 freplicaNav.ComputeTransformation(fBlockedReplicaNo,
259 fBlockedPhysicalVolume);
260 fHistory.NewLevel(fBlockedPhysicalVolume, kReplica,
261 fBlockedReplicaNo);
262 fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
263 break;
264 case kParameterised:
265 if( fBlockedPhysicalVolume->GetRegularStructureId() == 0 )
266 {
267 G4VSolid *pSolid;
268 G4VPVParameterisation *pParam;
269 G4TouchableHistory parentTouchable( fHistory );
270 pParam = fBlockedPhysicalVolume->GetParameterisation();
271 pSolid = pParam->ComputeSolid(fBlockedReplicaNo,
272 fBlockedPhysicalVolume);
273 pSolid->ComputeDimensions(pParam, fBlockedReplicaNo,
274 fBlockedPhysicalVolume);
275 pParam->ComputeTransformation(fBlockedReplicaNo,
276 fBlockedPhysicalVolume);
277 fHistory.NewLevel(fBlockedPhysicalVolume, kParameterised,
278 fBlockedReplicaNo);
279 fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
280 //
281 // Set the correct solid and material in Logical Volume
282 //
283 G4LogicalVolume *pLogical;
284 pLogical = fBlockedPhysicalVolume->GetLogicalVolume();
285 pLogical->SetSolid( pSolid );
286 pLogical->UpdateMaterial(pParam ->
287 ComputeMaterial(fBlockedReplicaNo,
288 fBlockedPhysicalVolume,
289 &parentTouchable));
290 }
291 break;
292 case kExternal:
293 G4Exception("G4ITNavigator2::LocateGlobalPointAndSetup()",
294 "GeomNav0001", FatalException,
295 "Not applicable for external volumes.");
296 break;
297 }
298 fEntering = false;
299 fBlockedPhysicalVolume = nullptr;
300 localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
301 notKnownContained = false;
302 }
303 }
304 else
305 {
306 fBlockedPhysicalVolume = nullptr;
307 fEntering = false;
308 fEnteredDaughter = false; // Full Step was not taken, did not enter
309 fExiting = false;
310 fExitedMother = false; // Full Step was not taken, did not exit
311 }
312 }
313 //
314 // Search from top of history up through geometry until
315 // containing volume found:
316 // If on
317 // o OUTSIDE - Back up level, not/no longer exiting volumes
318 // o SURFACE and EXITING - Back up level, setting new blocking no.s
319 // else
320 // o containing volume found
321 //
322
323 while (notKnownContained)
324 {
325 if ( fHistory.GetTopVolumeType()!=kReplica )
326 {
327 targetSolid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid();
328 localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
329 insideCode = targetSolid->Inside(localPoint);
330 #ifdef G4VERBOSE
331 if(( fVerbose == 1 ) && ( fCheck ))
332 {
333 G4String solidResponse = "-kInside-";
334 if (insideCode == kOutside)
335 solidResponse = "-kOutside-";
336 else if (insideCode == kSurface)
337 solidResponse = "-kSurface-";
338 G4cout << "*** G4ITNavigator2::LocateGlobalPointAndSetup(): ***" << G4endl
339 << " Invoked Inside() for solid: " << targetSolid->GetName()
340 << ". Solid replied: " << solidResponse << G4endl
341 << " For local point p: " << localPoint << G4endl;
342 }
343 #endif
344 }
345 else
346 {
347 insideCode = freplicaNav.BackLocate(fHistory, globalPoint, localPoint,
348 fExiting, notKnownContained);
349 // !CARE! if notKnownContained returns false then the point is within
350 // the containing placement volume of the replica(s). If insidecode
351 // will result in the history being backed up one level, then the
352 // local point returned is the point in the system of this new level
353 }
354
355
356 if ( insideCode==kOutside )
357 {
358 noLevelsExited++;
359 if ( fHistory.GetDepth() )
360 {
361 fBlockedPhysicalVolume = fHistory.GetTopVolume();
362 fBlockedReplicaNo = fHistory.GetTopReplicaNo();
363 fHistory.BackLevel();
364 fExiting = false;
365
366 if( noLevelsExited > 1 )
367 {
368 // The first transformation was done by the sub-navigator
369 //
370 const G4RotationMatrix* mRot = fBlockedPhysicalVolume->GetRotation();
371 if( mRot != nullptr )
372 {
373 fGrandMotherExitNormal *= (*mRot).inverse();
374 fChangedGrandMotherRefFrame= true;
375 }
376 }
377 }
378 else
379 {
380 fLastLocatedPointLocal = localPoint;
381 fLocatedOutsideWorld = true;
382 // No extra transformation for ExitNormal - is in frame of Top Volume
383 return nullptr; // Have exited world volume
384 }
385 }
386 else
387 if ( insideCode==kSurface )
388 {
389 G4bool isExiting = fExiting;
390 if( (!fExiting)&&considerDirection )
391 {
392 // Figure out whether we are exiting this level's volume
393 // by using the direction
394 //
395 G4bool directionExiting = false;
396 G4ThreeVector localDirection =
397 fHistory.GetTopTransform().TransformAxis(globalDirection);
398
399 // Make sure localPoint in correct reference frame
400 // ( Was it already correct ? How ? )
401 //
402 localPoint= fHistory.GetTopTransform().TransformPoint(globalPoint);
403 if ( fHistory.GetTopVolumeType()!=kReplica )
404 {
405 G4ThreeVector normal = targetSolid->SurfaceNormal(localPoint);
406 directionExiting = normal.dot(localDirection) > 0.0;
407 isExiting = isExiting || directionExiting;
408 }
409 }
410 if( isExiting )
411 {
412 noLevelsExited++;
413 if ( fHistory.GetDepth() )
414 {
415 fBlockedPhysicalVolume = fHistory.GetTopVolume();
416 fBlockedReplicaNo = fHistory.GetTopReplicaNo();
417 fHistory.BackLevel();
418 //
419 // Still on surface but exited volume not necessarily convex
420 //
421 fValidExitNormal = false;
422
423 if( noLevelsExited > 1 )
424 {
425 // The first transformation was done by the sub-navigator
426 //
427 const G4RotationMatrix* mRot =
428 fBlockedPhysicalVolume->GetRotation();
429 if( mRot != nullptr )
430 {
431 fGrandMotherExitNormal *= (*mRot).inverse();
432 fChangedGrandMotherRefFrame= true;
433 }
434 }
435 }
436 else
437 {
438 fLastLocatedPointLocal = localPoint;
439 fLocatedOutsideWorld = true;
440 // No extra transformation for ExitNormal, is in frame of Top Vol
441 return nullptr; // Have exited world volume
442 }
443 }
444 else
445 {
446 notKnownContained=false;
447 }
448 }
449 else
450 {
451 notKnownContained=false;
452 }
453 } // END while (notKnownContained)
454 //
455 // Search downwards until deepest containing volume found,
456 // blocking fBlockedPhysicalVolume/BlockedReplicaNum
457 //
458 // 3 Cases:
459 //
460 // o Parameterised daughters
461 // =>Must be one G4PVParameterised daughter & voxels
462 // o Positioned daughters & voxels
463 // o Positioned daughters & no voxels
464
465 noResult = true; // noResult should be renamed to
466 // something like enteredLevel, as that is its meaning.
467 do
468 {
469 // Determine `type' of current mother volume
470 //
471 targetPhysical = fHistory.GetTopVolume();
472 if (targetPhysical == nullptr) { break; }
473 targetLogical = targetPhysical->GetLogicalVolume();
474 switch( CharacteriseDaughters(targetLogical) )
475 {
476 case kNormal:
477 if ( targetLogical->GetVoxelHeader() != nullptr ) // use optimised navigation
478 {
479 noResult = fvoxelNav.LevelLocate(fHistory,
480 fBlockedPhysicalVolume,
481 fBlockedReplicaNo,
482 globalPoint,
483 pGlobalDirection,
484 considerDirection,
485 localPoint);
486 }
487 else // do not use optimised navigation
488 {
489 noResult = fnormalNav.LevelLocate(fHistory,
490 fBlockedPhysicalVolume,
491 fBlockedReplicaNo,
492 globalPoint,
493 pGlobalDirection,
494 considerDirection,
495 localPoint);
496 }
497 break;
498 case kReplica:
499 noResult = freplicaNav.LevelLocate(fHistory,
500 fBlockedPhysicalVolume,
501 fBlockedReplicaNo,
502 globalPoint,
503 pGlobalDirection,
504 considerDirection,
505 localPoint);
506 break;
507 case kParameterised:
508 if( GetDaughtersRegularStructureId(targetLogical) != 1 )
509 {
510 noResult = fparamNav.LevelLocate(fHistory,
511 fBlockedPhysicalVolume,
512 fBlockedReplicaNo,
513 globalPoint,
514 pGlobalDirection,
515 considerDirection,
516 localPoint);
517 }
518 else // Regular structure
519 {
520 noResult = fregularNav.LevelLocate(fHistory,
521 fBlockedPhysicalVolume,
522 fBlockedReplicaNo,
523 globalPoint,
524 pGlobalDirection,
525 considerDirection,
526 localPoint);
527 }
528 break;
529 case kExternal:
530 G4Exception("G4ITNavigator2::LocateGlobalPointAndSetup()",
531 "GeomNav0001", FatalException,
532 "Not applicable for external volumes.");
533 break;
534 }
535
536 // LevelLocate returns true if it finds a daughter volume
537 // in which globalPoint is inside (or on the surface).
538
539 if ( noResult )
540 {
541 // Entering a daughter after ascending
542 //
543 // The blocked volume is no longer valid - it was for another level
544 //
545 fBlockedPhysicalVolume = nullptr;
546 fBlockedReplicaNo = -1;
547
548 // fEntering should be false -- else blockedVolume is assumed good.
549 // fEnteredDaughter is used for ExitNormal
550 //
551 fEntering = false;
552 fEnteredDaughter = true;
553
554 if( fExitedMother )
555 {
556 G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume();
557 const G4RotationMatrix* mRot = enteredPhysical->GetRotation();
558 if( mRot != nullptr )
559 {
560 // Go deeper, i.e. move 'down' in the hierarchy
561 // Apply direct rotation, not inverse
562 //
563 fGrandMotherExitNormal *= (*mRot);
564 fChangedGrandMotherRefFrame= true;
565 }
566 }
567
568 #ifdef G4DEBUG_NAVIGATION
569 if( fVerbose > 2 )
570 {
571 G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume();
572 G4cout << "*** G4ITNavigator2::LocateGlobalPointAndSetup() ***" << G4endl;
573 G4cout << " Entering volume: " << enteredPhysical->GetName()
574 << G4endl;
575 }
576 #endif
577 }
578 } while (noResult);
579
580 fLastLocatedPointLocal = localPoint;
581
582 #ifdef G4VERBOSE
583 if( fVerbose >= 4 )
584 {
585 G4long oldcoutPrec = G4cout.precision(8);
586 G4String curPhysVol_Name("None");
587 if (targetPhysical != nullptr) { curPhysVol_Name = targetPhysical->GetName(); }
588 G4cout << " Return value = new volume = " << curPhysVol_Name << G4endl;
589 G4cout << " ----- Upon exiting:" << G4endl;
590 PrintState();
591 if( fVerbose >= 5 )
592 {
593 G4cout << "Upon exiting LocateGlobalPointAndSetup():" << G4endl;
594 G4cout << " History = " << G4endl << fHistory << G4endl << G4endl;
595 }
596 G4cout.precision(oldcoutPrec);
597 }
598 #endif
599
600 fLocatedOutsideWorld= false;
601
602 return targetPhysical;
603 }
604
605 // ********************************************************************
606 // LocateGlobalPointWithinVolume
607 //
608 // -> the state information of this Navigator and its subNavigators
609 // is updated in order to start the next step at pGlobalpoint
610 // -> no check is performed whether pGlobalpoint is inside the
611 // original volume (this must be the case).
612 //
613 // Note: a direction could be added to the arguments, to aid in future
614 // optional checking (via the old code below, flagged by OLD_LOCATE).
615 // [ This would be done only in verbose mode ]
616 // ********************************************************************
617 //
618 void
619 G4ITNavigator2::LocateGlobalPointWithinVolume(const G4ThreeVector& pGlobalpoint)
620 {
621 CheckNavigatorStateIsValid();
622
623 #ifdef G4DEBUG_NAVIGATION
624 // Check: Either step was not limited by a boundary
625 // or else the full step is no longer being taken
626 assert( !fWasLimitedByGeometry );
627 #endif
628
629 fLastLocatedPointLocal = ComputeLocalPoint(pGlobalpoint);
630 fLastTriedStepComputation= false;
631 fChangedGrandMotherRefFrame= false; // Frame for Exit Normal
632
633 #ifdef G4DEBUG_NAVIGATION
634 if( fVerbose > 2 )
635 {
636 G4cout << "Entering LocateGlobalWithinVolume(): History = " << G4endl;
637 G4cout << fHistory << G4endl;
638 }
639 #endif
640
641 // For the case of Voxel (or Parameterised) volume the respective
642 // Navigator must be messaged to update its voxel information etc
643
644 // Update the state of the Sub Navigators
645 // - in particular any voxel information they store/cache
646 //
647 G4VPhysicalVolume* motherPhysical = fHistory.GetTopVolume();
648 G4LogicalVolume* motherLogical = motherPhysical->GetLogicalVolume();
649 G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
650
651 if ( fHistory.GetTopVolumeType()!=kReplica )
652 {
653 switch( CharacteriseDaughters(motherLogical) )
654 {
655 case kNormal:
656 if ( pVoxelHeader != nullptr )
657 {
658 fvoxelNav.VoxelLocate( pVoxelHeader, fLastLocatedPointLocal );
659 }
660 break;
661 case kParameterised:
662 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
663 {
664 // Resets state & returns voxel node
665 //
666 fparamNav.ParamVoxelLocate( pVoxelHeader, fLastLocatedPointLocal );
667 }
668 break;
669 case kReplica:
670 G4Exception("G4ITNavigator2::LocateGlobalPointWithinVolume()",
671 "GeomNav0001", FatalException,
672 "Not applicable for replicated volumes.");
673 break;
674 case kExternal:
675 G4Exception("G4ITNavigator2::LocateGlobalPointWithinVolume()",
676 "GeomNav0001", FatalException,
677 "Not applicable for external volumes.");
678 break;
679 }
680 }
681
682 // Reset the state variables
683 // - which would have been affected
684 // by the 'equivalent' call to LocateGlobalPointAndSetup
685 // - who's values have been invalidated by the 'move'.
686 //
687 fBlockedPhysicalVolume = nullptr;
688 fBlockedReplicaNo = -1;
689 fEntering = false;
690 fEnteredDaughter = false; // Boundary not encountered, did not enter
691 fExiting = false;
692 fExitedMother = false; // Boundary not encountered, did not exit
693 }
694
695 // !>
696
697 void G4ITNavigator2::CheckNavigatorState() const
698 {
699 if(fpNavigatorState == nullptr)
700 {
701 G4ExceptionDescription exceptionDescription;
702 exceptionDescription << "The navigator state is NULL. ";
703 exceptionDescription << "Either NewNavigatorStateAndLocate was not called ";
704 exceptionDescription << "or the provided navigator state was already NULL.";
705
706 G4Exception("G4ITNavigator::CheckNavigatorStateIsValid",
707 "NavigatorStateNotValid",FatalException,exceptionDescription);
708 return;
709 }
710 }
711
712 G4ITNavigatorState_Lock2* G4ITNavigator2::GetNavigatorState()
713 {
714 return fpNavigatorState;
715 }
716
717 void G4ITNavigator2::SetNavigatorState(G4ITNavigatorState_Lock2* navState)
718 {
719 fpNavigatorState = (G4NavigatorState*) navState;
720 if(fpNavigatorState != nullptr) SetupHierarchy(); // ?
721 // fpSaveState = (G4SaveNavigatorState*) navState;
722 // if(navState) RestoreSavedState();
723 }
724
725 void G4ITNavigator2::ResetNavigatorState()
726 {
727 fpNavigatorState = nullptr;
728 }
729
730 void G4ITNavigator2::NewNavigatorState()
731 {
732 fpNavigatorState = new G4NavigatorState();
733 if(fTopPhysical == nullptr)
734 {
735 G4ExceptionDescription exceptionDescription;
736 exceptionDescription << "No World Volume";
737
738 G4Exception("G4ITNavigator::NewNavigatorState",
739 "NoWorldVolume",FatalException,exceptionDescription);
740 return;
741 }
742
743 fHistory.SetFirstEntry(fTopPhysical );
744 SetupHierarchy();
745 }
746
747 void G4ITNavigator2::NewNavigatorState(const G4TouchableHistory &h)
748 {
749 fpNavigatorState = new G4NavigatorState();
750 if(fTopPhysical == nullptr)
751 {
752 G4ExceptionDescription exceptionDescription;
753 exceptionDescription << "No World Volume";
754
755 G4Exception("G4ITNavigator::NewNavigatorState",
756 "NoWorldVolume",FatalException,exceptionDescription);
757 return;
758 }
759
760 fHistory = *h.GetHistory();
761 fLastTriedStepComputation= false; // Redundant, but best
762 SetupHierarchy();
763 }
764
765 G4VPhysicalVolume* G4ITNavigator2::NewNavigatorStateAndLocate(const G4ThreeVector &p,
766 const G4ThreeVector &direction)
767 {
768 fpNavigatorState = new G4NavigatorState();
769
770 if(fTopPhysical == nullptr)
771 {
772 G4ExceptionDescription exceptionDescription;
773 exceptionDescription << "No World Volume";
774
775 G4Exception("G4ITNavigator::NewNavigatorStateAndLocate",
776 "NoWorldVolume",FatalException,exceptionDescription);
777 return nullptr;
778 }
779
780 fHistory.SetFirstEntry(fTopPhysical );
781 SetupHierarchy();
782 return LocateGlobalPointAndSetup(p, &direction, false, false);
783 }
784
785 // ********************************************************************
786 // SetSavedState
787 //
788 // Save the state, in case this is a parasitic call
789 // Save fValidExitNormal, fExitNormal, fExiting, fEntering,
790 // fBlockedPhysicalVolume, fBlockedReplicaNo, fLastStepWasZero;
791 // ********************************************************************
792 //
793 /*
794 void G4ITNavigator2::SetSavedState()
795 {
796 fSaveState = *fpNavigatorState;
797 }
798 */
799
800 /*
801 void G4ITNavigator2::SetSavedState()
802 {
803 // !>
804 // This check can be avoid if instead, at every first step of a track,
805 // the IT tracking uses NewNavigatorSate
806 // The normal tracking would just call once NewNavigatorState() before tracking
807
808 // if(fpSaveState == 0)
809 // fpSaveState = new G4SaveNavigatorState;
810 // <!
811
812 // fSaveExitNormal = fExitNormal;
813 fpSaveState->sExitNormal = fExitNormal;
814 fpSaveState->sValidExitNormal = fValidExitNormal;
815 fpSaveState->sExiting = fExiting;
816 fpSaveState->sEntering = fEntering;
817
818 fpSaveState->spBlockedPhysicalVolume = fBlockedPhysicalVolume;
819 fpSaveState->sBlockedReplicaNo = fBlockedReplicaNo,
820
821 fpSaveState->sLastStepWasZero = fLastStepWasZero;
822
823 // !>
824 fpSaveState->sPreviousSftOrigin = fPreviousSftOrigin;
825 fpSaveState->sPreviousSafety = fPreviousSafety;
826 fpSaveState->sNumberZeroSteps = fNumberZeroSteps;
827 fpSaveState->sLocatedOnEdge = fLocatedOnEdge;
828 fpSaveState->sWasLimitedByGeometry= fWasLimitedByGeometry;
829 fpSaveState->sPushed=fPushed;
830 fpSaveState->sNumberZeroSteps=fNumberZeroSteps;
831 fpSaveState->sEnteredDaughter = fEnteredDaughter;
832 fpSaveState->sExitedMother = fExitedMother;
833
834 fpSaveState->sLastLocatedPointLocal = fLastLocatedPointLocal;
835 fpSaveState->sLocatedOutsideWorld = fLocatedOutsideWorld;
836 // <!
837 }
838 */
839 // ********************************************************************
840 // RestoreSavedState
841 //
842 // Restore the state (in Compute Step), in case this is a parasitic call
843 // ********************************************************************
844 //
845 /*
846 void G4ITNavigator2::RestoreSavedState()
847 {
848 *fpNavigatorState = fSaveState;
849 }
850 */
851 /*
852 void G4ITNavigator2::RestoreSavedState()
853 {
854 fExitNormal = fpSaveState->sExitNormal;
855 fValidExitNormal = fpSaveState->sValidExitNormal;
856 fExiting = fpSaveState->sExiting;
857 fEntering = fpSaveState->sEntering;
858
859 fBlockedPhysicalVolume = fpSaveState->spBlockedPhysicalVolume;
860 fBlockedReplicaNo = fpSaveState->sBlockedReplicaNo,
861
862 fLastStepWasZero = fpSaveState->sLastStepWasZero;
863
864 // !>
865 fPreviousSftOrigin = fpSaveState->sPreviousSftOrigin ;
866 fPreviousSafety = fpSaveState->sPreviousSafety ;
867 fNumberZeroSteps = fpSaveState->sNumberZeroSteps ;
868 fLocatedOnEdge = fpSaveState->sLocatedOnEdge ;
869 fWasLimitedByGeometry = fpSaveState->sWasLimitedByGeometry;
870 fPushed = fpSaveState->sPushed;
871 fNumberZeroSteps = fpSaveState->sNumberZeroSteps;
872 fEnteredDaughter= fpSaveState->sEnteredDaughter ;
873 fExitedMother = fpSaveState->sExitedMother ;
874
875 fLastLocatedPointLocal = fpSaveState->sLastLocatedPointLocal ;
876 fLocatedOutsideWorld = fpSaveState->sLocatedOutsideWorld;
877 // <!
878 }
879 */
880 // <!
881
882 // ********************************************************************
883 // ComputeStep
884 //
885 // Computes the next geometric Step: intersections with current
886 // mother and `daughter' volumes.
887 //
888 // NOTE:
889 //
890 // Flags on entry:
891 // --------------
892 // fValidExitNormal - Normal of exited volume is valid (convex, not a
893 // coincident boundary)
894 // fExitNormal - Surface normal of exited volume
895 // fExiting - True if have exited solid
896 //
897 // fBlockedPhysicalVolume - Ptr to exited volume (or 0)
898 // fBlockedReplicaNo - Replication no of exited volume
899 // fLastStepWasZero - True if last Step size was zero.
900 //
901 // Flags on exit:
902 // -------------
903 // fValidExitNormal - True if surface normal of exited volume is valid
904 // fExitNormal - Surface normal of exited volume rotated to mothers
905 // reference system
906 // fExiting - True if exiting mother
907 // fEntering - True if entering `daughter' volume (or replica)
908 // fBlockedPhysicalVolume - Ptr to candidate (entered) volume
909 // fBlockedReplicaNo - Replication no of candidate (entered) volume
910 // fLastStepWasZero - True if this Step size was zero.
911 // ********************************************************************
912 //
913 G4double G4ITNavigator2::ComputeStep( const G4ThreeVector &pGlobalpoint,
914 const G4ThreeVector &pDirection,
915 const G4double pCurrentProposedStepLength,
916 G4double &pNewSafety)
917 {
918 CheckNavigatorStateIsValid();
919 G4ThreeVector localDirection = ComputeLocalAxis(pDirection);
920 G4double Step = kInfinity;
921 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
922 G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
923
924 // All state relating to exiting normals must be reset
925 //
926 fExitNormalGlobalFrame= G4ThreeVector( 0., 0., 0.);
927 // Reset value - to erase its memory
928 fChangedGrandMotherRefFrame= false;
929 // Reset - used for local exit normal
930 fGrandMotherExitNormal= G4ThreeVector( 0., 0., 0.);
931 fCalculatedExitNormal = false;
932 // Reset for new step
933
934 fLastTriedStepComputation= true;
935
936 #ifdef G4VERBOSE
937 if( fVerbose > 0 )
938 {
939 G4cout << "*** G4ITNavigator2::ComputeStep: ***" << G4endl;
940 G4cout << " Volume = " << motherPhysical->GetName()
941 << " - Proposed step length = " << pCurrentProposedStepLength
942 << G4endl;
943 #ifdef G4DEBUG_NAVIGATION
944 if( fVerbose >= 2 )
945 {
946 G4cout << " Called with the arguments: " << G4endl
947 << " Globalpoint = " << std::setw(25) << pGlobalpoint << G4endl
948 << " Direction = " << std::setw(25) << pDirection << G4endl;
949 if( fVerbose >= 4 )
950 {
951 G4cout << " ---- Upon entering : State" << G4endl;
952 PrintState();
953 }
954 }
955 #endif
956 }
957 #endif
958
959 G4ThreeVector newLocalPoint = ComputeLocalPoint(pGlobalpoint);
960 if( newLocalPoint != fLastLocatedPointLocal )
961 {
962 // Check whether the relocation is within safety
963 //
964 G4ThreeVector oldLocalPoint = fLastLocatedPointLocal;
965 G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2();
966
967 if ( moveLenSq >= kCarTolerance*kCarTolerance )
968 {
969 #ifdef G4VERBOSE
970 ComputeStepLog(pGlobalpoint, moveLenSq);
971 #endif
972 // Relocate the point within the same volume
973 //
974 LocateGlobalPointWithinVolume( pGlobalpoint );
975 fLastTriedStepComputation= true; // Ensure that this is set again !!
976 }
977 }
978 if ( fHistory.GetTopVolumeType()!=kReplica )
979 {
980 switch( CharacteriseDaughters(motherLogical) )
981 {
982 case kNormal:
983 if ( motherLogical->GetVoxelHeader() != nullptr )
984 {
985 LocateGlobalPointWithinVolume(pGlobalpoint);
986 Step = fvoxelNav.ComputeStep(fLastLocatedPointLocal,
987 localDirection,
988 pCurrentProposedStepLength,
989 pNewSafety,
990 fHistory,
991 fValidExitNormal,
992 fExitNormal,
993 fExiting,
994 fEntering,
995 &fBlockedPhysicalVolume,
996 fBlockedReplicaNo);
997
998 }
999 else
1000 {
1001 if( motherPhysical->GetRegularStructureId() == 0 )
1002 {
1003 Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
1004 localDirection,
1005 pCurrentProposedStepLength,
1006 pNewSafety,
1007 fHistory,
1008 fValidExitNormal,
1009 fExitNormal,
1010 fExiting,
1011 fEntering,
1012 &fBlockedPhysicalVolume,
1013 fBlockedReplicaNo);
1014 }
1015 else // Regular (non-voxelised) structure
1016 {
1017 LocateGlobalPointAndSetup( pGlobalpoint, &pDirection, true, true );
1018 fLastTriedStepComputation= true; // Ensure that this is set again !!
1019 //
1020 // if physical process limits the step, the voxel will not be the
1021 // one given by ComputeStepSkippingEqualMaterials() and the local
1022 // point will be wrongly calculated.
1023
1024 // There is a problem: when msc limits the step and the point is
1025 // assigned wrongly to phantom in previous step (while it is out
1026 // of the container volume). Then LocateGlobalPointAndSetup() has
1027 // reset the history topvolume to world.
1028 //
1029 if(fHistory.GetTopVolume()->GetRegularStructureId() == 0 )
1030 {
1031 G4Exception("G4ITNavigator2::ComputeStep()",
1032 "GeomNav1001", JustWarning,
1033 "Point is relocated in voxels, while it should be outside!");
1034 Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
1035 localDirection,
1036 pCurrentProposedStepLength,
1037 pNewSafety,
1038 fHistory,
1039 fValidExitNormal,
1040 fExitNormal,
1041 fExiting,
1042 fEntering,
1043 &fBlockedPhysicalVolume,
1044 fBlockedReplicaNo);
1045 }
1046 else
1047 {
1048 Step = fregularNav.
1049 ComputeStepSkippingEqualMaterials(fLastLocatedPointLocal,
1050 localDirection,
1051 pCurrentProposedStepLength,
1052 pNewSafety,
1053 fHistory,
1054 fValidExitNormal,
1055 fExitNormal,
1056 fExiting,
1057 fEntering,
1058 &fBlockedPhysicalVolume,
1059 fBlockedReplicaNo,
1060 motherPhysical);
1061 }
1062 }
1063 }
1064 break;
1065 case kParameterised:
1066 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
1067 {
1068 Step = fparamNav.ComputeStep(fLastLocatedPointLocal,
1069 localDirection,
1070 pCurrentProposedStepLength,
1071 pNewSafety,
1072 fHistory,
1073 fValidExitNormal,
1074 fExitNormal,
1075 fExiting,
1076 fEntering,
1077 &fBlockedPhysicalVolume,
1078 fBlockedReplicaNo);
1079 }
1080 else // Regular structure
1081 {
1082 Step = fregularNav.ComputeStep(fLastLocatedPointLocal,
1083 localDirection,
1084 pCurrentProposedStepLength,
1085 pNewSafety,
1086 fHistory,
1087 fValidExitNormal,
1088 fExitNormal,
1089 fExiting,
1090 fEntering,
1091 &fBlockedPhysicalVolume,
1092 fBlockedReplicaNo);
1093 }
1094 break;
1095 case kReplica:
1096 G4Exception("G4ITNavigator2::ComputeStep()", "GeomNav0001",
1097 FatalException, "Not applicable for replicated volumes.");
1098 break;
1099 case kExternal:
1100 G4Exception("G4ITNavigator2::ComputeStep()", "GeomNav0001",
1101 FatalException, "Not applicable for external volumes.");
1102 break;
1103 }
1104 }
1105 else
1106 {
1107 // In the case of a replica, it must handle the exiting
1108 // edge/corner problem by itself
1109 //
1110 G4bool exitingReplica = fExitedMother;
1111 G4bool calculatedExitNormal;
1112 Step = freplicaNav.ComputeStep(pGlobalpoint,
1113 pDirection,
1114 fLastLocatedPointLocal,
1115 localDirection,
1116 pCurrentProposedStepLength,
1117 pNewSafety,
1118 fHistory,
1119 fValidExitNormal,
1120 calculatedExitNormal,
1121 fExitNormal,
1122 exitingReplica,
1123 fEntering,
1124 &fBlockedPhysicalVolume,
1125 fBlockedReplicaNo);
1126 fExiting= exitingReplica;
1127 fCalculatedExitNormal= calculatedExitNormal;
1128 }
1129
1130
1131 // G4cout << " !!!! Step = " << Step << G4endl;
1132
1133 // Remember last safety origin & value.
1134 //
1135 fPreviousSftOrigin = pGlobalpoint;
1136 fPreviousSafety = pNewSafety;
1137
1138 // Count zero steps - one can occur due to changing momentum at a boundary
1139 // - one, two (or a few) can occur at common edges between
1140 // volumes
1141 // - more than two is likely a problem in the geometry
1142 // description or the Navigation
1143
1144 // Rule of thumb: likely at an Edge if two consecutive steps are zero,
1145 // because at least two candidate volumes must have been
1146 // checked
1147 //
1148 fLocatedOnEdge = fLastStepWasZero && (Step==0.0);
1149 fLastStepWasZero = (Step==0.0);
1150 if (fPushed) { fPushed = fLastStepWasZero; }
1151
1152 // Handle large number of consecutive zero steps
1153 //
1154 if ( fLastStepWasZero )
1155 {
1156 fNumberZeroSteps++;
1157 #ifdef G4DEBUG_NAVIGATION
1158 if( fNumberZeroSteps > 1 )
1159 {
1160 G4cout << "G4ITNavigator2::ComputeStep(): another zero step, # "
1161 << fNumberZeroSteps
1162 << " at " << pGlobalpoint
1163 << " in volume " << motherPhysical->GetName()
1164 << G4endl;
1165 }
1166 #endif
1167 if( fNumberZeroSteps > fActionThreshold_NoZeroSteps-1 )
1168 {
1169 // Act to recover this stuck track. Pushing it along direction
1170 //
1171 Step += 100*kCarTolerance;
1172 #ifdef G4VERBOSE
1173 if ((!fPushed) && (fWarnPush))
1174 {
1175 std::ostringstream message;
1176 message << "Track stuck or not moving." << G4endl
1177 << " Track stuck, not moving for "
1178 << fNumberZeroSteps << " steps" << G4endl
1179 << " in volume -" << motherPhysical->GetName()
1180 << "- at point " << pGlobalpoint << G4endl
1181 << " direction: " << pDirection << "." << G4endl
1182 << " Potential geometry or navigation problem !"
1183 << G4endl
1184 << " Trying pushing it of " << Step << " mm ...";
1185 G4Exception("G4ITNavigator2::ComputeStep()", "GeomNav1002",
1186 JustWarning, message, "Potential overlap in geometry!");
1187 }
1188 #endif
1189 fPushed = true;
1190 }
1191 if( fNumberZeroSteps > fAbandonThreshold_NoZeroSteps-1 )
1192 {
1193 // Must kill this stuck track
1194 //
1195 std::ostringstream message;
1196 message << "Stuck Track: potential geometry or navigation problem."
1197 << G4endl
1198 << " Track stuck, not moving for "
1199 << fNumberZeroSteps << " steps" << G4endl
1200 << " in volume -" << motherPhysical->GetName()
1201 << "- at point " << pGlobalpoint << G4endl
1202 << " direction: " << pDirection << ".";
1203 motherPhysical->CheckOverlaps(5000, 0.0);
1204 G4Exception("G4ITNavigator2::ComputeStep()", "GeomNav0003",
1205 EventMustBeAborted, message);
1206 }
1207 }
1208 else
1209 {
1210 if (!fPushed) fNumberZeroSteps = 0;
1211 }
1212
1213 fEnteredDaughter = fEntering; // I expect to enter a volume in this Step
1214 fExitedMother = fExiting;
1215
1216 fStepEndPoint = pGlobalpoint
1217 + std::min(Step,pCurrentProposedStepLength) * pDirection;
1218 fLastStepEndPointLocal = fLastLocatedPointLocal + Step * localDirection;
1219
1220 if( fExiting )
1221 {
1222 #ifdef G4DEBUG_NAVIGATION
1223 if( fVerbose > 2 )
1224 {
1225 G4cout << " At G4Nav CompStep End - if(exiting) - fExiting= " << fExiting
1226 << " fValidExitNormal = " << fValidExitNormal << G4endl;
1227 G4cout << " fExitNormal= " << fExitNormal << G4endl;
1228 }
1229 #endif
1230
1231 if(fValidExitNormal || fCalculatedExitNormal)
1232 {
1233 if ( fHistory.GetTopVolumeType()!=kReplica )
1234 {
1235 // Convention: fExitNormal is in the 'grand-mother' coordinate system
1236 //
1237 fGrandMotherExitNormal= fExitNormal;
1238 fCalculatedExitNormal= true;
1239 }
1240 else
1241 {
1242 fGrandMotherExitNormal = fExitNormal;
1243 }
1244 }
1245 else
1246 {
1247 // We must calculate the normal anyway (in order to have it if requested)
1248 //
1249 G4ThreeVector finalLocalPoint =
1250 fLastLocatedPointLocal + localDirection*Step;
1251
1252 if ( fHistory.GetTopVolumeType()!=kReplica )
1253 {
1254 // Find normal in the 'mother' coordinate system
1255 //
1256 G4ThreeVector exitNormalMotherFrame=
1257 motherLogical->GetSolid()->SurfaceNormal(finalLocalPoint);
1258
1259 // Transform it to the 'grand-mother' coordinate system
1260 //
1261 const G4RotationMatrix* mRot = motherPhysical->GetRotation();
1262 if( mRot != nullptr )
1263 {
1264 fChangedGrandMotherRefFrame= true;
1265 fGrandMotherExitNormal = (*mRot).inverse() * exitNormalMotherFrame;
1266 }
1267 else
1268 {
1269 fGrandMotherExitNormal = exitNormalMotherFrame;
1270 }
1271
1272 // Do not set fValidExitNormal -- this signifies
1273 // that the solid is convex!
1274 //
1275 fCalculatedExitNormal= true;
1276 }
1277 else
1278 {
1279 fCalculatedExitNormal = false;
1280 //
1281 // Nothing can be done at this stage currently - to solve this
1282 // Replica Navigation must have calculated the normal for this case
1283 // already.
1284 // Cases: mother is not convex, and exit is at previous replica level
1285
1286 #ifdef G4DEBUG_NAVIGATION
1287 G4ExceptionDescription desc;
1288
1289 desc << "Problem in ComputeStep: Replica Navigation did not provide"
1290 << " valid exit Normal. " << G4endl;
1291 desc << " Do not know how calculate it in this case." << G4endl;
1292 desc << " Location = " << finalLocalPoint << G4endl;
1293 desc << " Volume name = " << motherPhysical->GetName()
1294 << " copy/replica No = " << motherPhysical->GetCopyNo() << G4endl;
1295 G4Exception("G4ITNavigator2::ComputeStep()", "GeomNav0003",
1296 JustWarning, desc, "Normal not available for exiting.");
1297 #endif
1298 }
1299 }
1300
1301 // Now transform it to the global reference frame !!
1302 //
1303 if( fValidExitNormal || fCalculatedExitNormal )
1304 {
1305 auto depth= (G4int)fHistory.GetDepth();
1306 if( depth > 0 )
1307 {
1308 G4AffineTransform GrandMotherToGlobalTransf =
1309 fHistory.GetTransform(depth-1).Inverse();
1310 fExitNormalGlobalFrame =
1311 GrandMotherToGlobalTransf.TransformAxis( fGrandMotherExitNormal );
1312 }
1313 else
1314 {
1315 fExitNormalGlobalFrame= fGrandMotherExitNormal;
1316 }
1317 }
1318 else
1319 {
1320 fExitNormalGlobalFrame= G4ThreeVector( 0., 0., 0.);
1321 }
1322 }
1323
1324 if( (Step == pCurrentProposedStepLength) && (!fExiting) && (!fEntering) )
1325 {
1326 // This if Step is not really limited by the geometry.
1327 // The Navigator is obliged to return "infinity"
1328 //
1329 Step = kInfinity;
1330 }
1331
1332 #ifdef G4VERBOSE
1333 if( fVerbose > 1 )
1334 {
1335 if( fVerbose >= 4 )
1336 {
1337 G4cout << " ----- Upon exiting :" << G4endl;
1338 PrintState();
1339 }
1340 G4cout << " Returned step= " << Step;
1341 if( fVerbose > 5 ) G4cout << G4endl;
1342 if( Step == kInfinity )
1343 {
1344 G4cout << " Requested step= " << pCurrentProposedStepLength ;
1345 if( fVerbose > 5) G4cout << G4endl;
1346 }
1347 G4cout << " Safety = " << pNewSafety << G4endl;
1348 }
1349 #endif
1350
1351 return Step;
1352 }
1353
1354 // ********************************************************************
1355 // CheckNextStep
1356 //
1357 // Compute the step without altering the navigator state
1358 // ********************************************************************
1359 //
1360 G4double G4ITNavigator2::CheckNextStep( const G4ThreeVector& pGlobalpoint,
1361 const G4ThreeVector& pDirection,
1362 const G4double pCurrentProposedStepLength,
1363 G4double& pNewSafety)
1364 {
1365 CheckNavigatorStateIsValid();
1366 G4double step;
1367
1368 // Save the state, for this parasitic call
1369 //
1370 //SetSavedState();
1371 // G4SaveNavigatorState savedState(fpNavigatorState);
1372 G4NavigatorState savedState(*fpNavigatorState);
1373
1374 step = ComputeStep ( pGlobalpoint,
1375 pDirection,
1376 pCurrentProposedStepLength,
1377 pNewSafety );
1378
1379 // If a parasitic call, then attempt to restore the key parts of the state
1380 //
1381 *fpNavigatorState = savedState;
1382 //RestoreSavedState();
1383 // NOTE: the state of the current subnavigator is NOT restored.
1384 // ***> TODO: restore subnavigator state
1385 // if( last_located) Need Position of last location
1386 // if( last_computed step) Need Endposition of last step
1387
1388
1389 return step;
1390 }
1391
1392 // ********************************************************************
1393 // ResetState
1394 //
1395 // Resets stack and minimum of navigator state `machine'
1396 // ********************************************************************
1397 //
1398 void G4ITNavigator2::ResetState()
1399 {
1400 fWasLimitedByGeometry = false;
1401 fEntering = false;
1402 fExiting = false;
1403 fLocatedOnEdge = false;
1404 fLastStepWasZero = false;
1405 fEnteredDaughter = false;
1406 fExitedMother = false;
1407 fPushed = false;
1408
1409 fValidExitNormal = false;
1410 fChangedGrandMotherRefFrame= false;
1411 fCalculatedExitNormal = false;
1412
1413 fExitNormal = G4ThreeVector(0,0,0);
1414 fGrandMotherExitNormal = G4ThreeVector(0,0,0);
1415 fExitNormalGlobalFrame = G4ThreeVector(0,0,0);
1416
1417 fPreviousSftOrigin = G4ThreeVector(0,0,0);
1418 fPreviousSafety = 0.0;
1419
1420 fNumberZeroSteps = 0;
1421
1422 fBlockedPhysicalVolume = nullptr;
1423 fBlockedReplicaNo = -1;
1424
1425 fLastLocatedPointLocal = G4ThreeVector( kInfinity, -kInfinity, 0.0 );
1426 fLocatedOutsideWorld = false;
1427 }
1428
1429 // ********************************************************************
1430 // SetupHierarchy
1431 //
1432 // Renavigates & resets hierarchy described by current history
1433 // o Reset volumes
1434 // o Recompute transforms and/or solids of replicated/parameterised volumes
1435 // ********************************************************************
1436 //
1437 void G4ITNavigator2::SetupHierarchy()
1438 {
1439 G4int i;
1440 const auto cdepth = (G4int)fHistory.GetDepth();
1441 G4VPhysicalVolume *current;
1442 G4VSolid *pSolid;
1443 G4VPVParameterisation *pParam;
1444
1445 for ( i=1; i<=cdepth; i++ )
1446 {
1447 current = fHistory.GetVolume(i);
1448 switch ( fHistory.GetVolumeType(i) )
1449 {
1450 case kNormal:
1451 break;
1452 case kReplica:
1453 freplicaNav.ComputeTransformation(fHistory.GetReplicaNo(i), current);
1454 break;
1455 case kParameterised:
1456 {
1457 G4int replicaNo;
1458 pParam = current->GetParameterisation();
1459 replicaNo = fHistory.GetReplicaNo(i);
1460 pSolid = pParam->ComputeSolid(replicaNo, current);
1461
1462 // Set up dimensions & transform in solid/physical volume
1463 //
1464 pSolid->ComputeDimensions(pParam, replicaNo, current);
1465 pParam->ComputeTransformation(replicaNo, current);
1466
1467 G4TouchableHistory *pTouchable= nullptr;
1468 if( pParam->IsNested() )
1469 {
1470 pTouchable= new G4TouchableHistory( fHistory );
1471 pTouchable->MoveUpHistory(); // Move up to the parent level
1472 // Adequate only if Nested at the Branch level (last)
1473 // To extend to other cases:
1474 // pTouchable->MoveUpHistory(cdepth-i-1);
1475 // Move to the parent level of *Current* level
1476 // Could replace this line and constructor with a revised
1477 // c-tor for History(levels to drop)
1478 }
1479 // Set up the correct solid and material in Logical Volume
1480 //
1481 G4LogicalVolume *pLogical = current->GetLogicalVolume();
1482 pLogical->SetSolid( pSolid );
1483 pLogical->UpdateMaterial( pParam ->
1484 ComputeMaterial(replicaNo, current, pTouchable) );
1485 delete pTouchable;
1486 }
1487 break;
1488
1489 case kExternal:
1490 G4Exception("G4ITNavigator2::SetupHierarchy()",
1491 "GeomNav0001", FatalException,
1492 "Not applicable for external volumes.");
1493 break;
1494
1495 }
1496 }
1497 }
1498
1499 // ********************************************************************
1500 // GetLocalExitNormal
1501 //
1502 // Obtains the Normal vector to a surface (in local coordinates)
1503 // pointing out of previous volume and into current volume
1504 // ********************************************************************
1505 //
1506 G4ThreeVector G4ITNavigator2::GetLocalExitNormal( G4bool* valid )
1507 {
1508 CheckNavigatorStateIsValid();
1509 G4ThreeVector ExitNormal(0.,0.,0.);
1510 G4VSolid *currentSolid=nullptr;
1511 G4LogicalVolume *candidateLogical;
1512 if ( fLastTriedStepComputation )
1513 {
1514 // use fLastLocatedPointLocal and next candidate volume
1515 //
1516 G4ThreeVector nextSolidExitNormal(0.,0.,0.);
1517
1518 if( fEntering && (fBlockedPhysicalVolume!=nullptr) )
1519 {
1520 candidateLogical= fBlockedPhysicalVolume->GetLogicalVolume();
1521 if( candidateLogical != nullptr )
1522 {
1523 // fLastStepEndPointLocal is in the coordinates of the mother
1524 // we need it in the daughter's coordinate system.
1525
1526 // The following code should also work in case of Replica
1527 {
1528 // First transform fLastLocatedPointLocal to the new daughter
1529 // coordinates
1530 //
1531 G4AffineTransform MotherToDaughterTransform=
1532 GetMotherToDaughterTransform( fBlockedPhysicalVolume,
1533 fBlockedReplicaNo,
1534 VolumeType(fBlockedPhysicalVolume) );
1535 G4ThreeVector daughterPointOwnLocal=
1536 MotherToDaughterTransform.TransformPoint( fLastStepEndPointLocal );
1537
1538 // OK if it is a parameterised volume
1539 //
1540 EInside inSideIt;
1541 G4bool onSurface;
1542 G4double safety= -1.0;
1543 currentSolid= candidateLogical->GetSolid();
1544 inSideIt = currentSolid->Inside(daughterPointOwnLocal);
1545 onSurface = (inSideIt == kSurface);
1546 if( ! onSurface )
1547 {
1548 if( inSideIt == kOutside )
1549 {
1550 safety = (currentSolid->DistanceToIn(daughterPointOwnLocal));
1551 onSurface = safety < 100.0 * kCarTolerance;
1552 }
1553 else if (inSideIt == kInside )
1554 {
1555 safety = (currentSolid->DistanceToOut(daughterPointOwnLocal));
1556 onSurface = safety < 100.0 * kCarTolerance;
1557 }
1558 }
1559
1560 if( onSurface )
1561 {
1562 nextSolidExitNormal =
1563 currentSolid->SurfaceNormal(daughterPointOwnLocal);
1564
1565 // Entering the solid ==> opposite
1566 //
1567 ExitNormal = -nextSolidExitNormal;
1568 fCalculatedExitNormal= true;
1569 }
1570 else
1571 {
1572 #ifdef G4VERBOSE
1573 if(( fVerbose == 1 ) && ( fCheck ))
1574 {
1575 std::ostringstream message;
1576 message << "Point not on surface ! " << G4endl
1577 << " Point = "
1578 << daughterPointOwnLocal << G4endl
1579 << " Physical volume = "
1580 << fBlockedPhysicalVolume->GetName() << G4endl
1581 << " Logical volume = "
1582 << candidateLogical->GetName() << G4endl
1583 << " Solid = " << currentSolid->GetName()
1584 << " Type = "
1585 << currentSolid->GetEntityType() << G4endl
1586 << *currentSolid << G4endl;
1587 if( inSideIt == kOutside )
1588 {
1589 message << "Point is Outside. " << G4endl
1590 << " Safety (from outside) = " << safety << G4endl;
1591 }
1592 else // if( inSideIt == kInside )
1593 {
1594 message << "Point is Inside. " << G4endl
1595 << " Safety (from inside) = " << safety << G4endl;
1596 }
1597 G4Exception("G4ITNavigator2::GetLocalExitNormal()", "GeomNav1001",
1598 JustWarning, message);
1599 }
1600 #endif
1601 }
1602 *valid = onSurface; // was =true;
1603 }
1604 }
1605 }
1606 else if ( fExiting )
1607 {
1608 ExitNormal = fGrandMotherExitNormal;
1609 *valid = true;
1610 fCalculatedExitNormal= true; // Should be true already
1611 }
1612 else // i.e. ( fBlockedPhysicalVolume == 0 )
1613 {
1614 *valid = false;
1615 G4Exception("G4ITNavigator2::GetLocalExitNormal()",
1616 "GeomNav0003", JustWarning,
1617 "Incorrect call to GetLocalSurfaceNormal." );
1618 }
1619 }
1620 else // ( ! fLastTriedStepComputation ) ie. last call was to Locate
1621 {
1622 if ( EnteredDaughterVolume() )
1623 {
1624 G4VSolid* daughterSolid =fHistory.GetTopVolume()->GetLogicalVolume()
1625 ->GetSolid();
1626 ExitNormal= -(daughterSolid->SurfaceNormal(fLastLocatedPointLocal));
1627 if( std::fabs(ExitNormal.mag2()-1.0 ) > CLHEP::perMillion )
1628 {
1629 G4ExceptionDescription desc;
1630 desc << " Parameters of solid: " << *daughterSolid
1631 << " Point for surface = " << fLastLocatedPointLocal << std::endl;
1632 G4Exception("G4ITNavigator2::GetLocalExitNormal()",
1633 "GeomNav0003", FatalException, desc,
1634 "Surface Normal returned by Solid is not a Unit Vector." );
1635 }
1636 fCalculatedExitNormal= true;
1637 *valid = true;
1638 }
1639 else
1640 {
1641 if( fExitedMother )
1642 {
1643 ExitNormal = fGrandMotherExitNormal;
1644 *valid = true;
1645 fCalculatedExitNormal= true;
1646 }
1647 else // We are not at a boundary. ExitNormal remains (0,0,0)
1648 {
1649 *valid = false;
1650 fCalculatedExitNormal= false;
1651 G4ExceptionDescription message;
1652 message << "Function called when *NOT* at a Boundary." << G4endl;
1653 G4Exception("G4ITNavigator2::GetLocalExitNormal()",
1654 "GeomNav0003", JustWarning, message);
1655 }
1656 }
1657 }
1658 return ExitNormal;
1659 }
1660
1661 // ********************************************************************
1662 // GetMotherToDaughterTransform
1663 //
1664 // Obtains the mother to daughter affine transformation
1665 // ********************************************************************
1666 //
1667 G4AffineTransform
1668 G4ITNavigator2::GetMotherToDaughterTransform( G4VPhysicalVolume *pEnteringPhysVol, // not Const
1669 G4int enteringReplicaNo,
1670 EVolume enteringVolumeType )
1671 {
1672 CheckNavigatorStateIsValid();
1673 switch (enteringVolumeType)
1674 {
1675 case kNormal: // Nothing is needed to prepare the transformation
1676 break; // It is stored already in the physical volume (placement)
1677 case kReplica: // Sets the transform in the Replica - tbc
1678 G4Exception("G4ITNavigator2::GetMotherToDaughterTransform()",
1679 "GeomNav0001", FatalException,
1680 "Method NOT Implemented yet for replica volumes.");
1681 break;
1682 case kParameterised:
1683 if( pEnteringPhysVol->GetRegularStructureId() == 0 )
1684 {
1685 G4VPVParameterisation *pParam =
1686 pEnteringPhysVol->GetParameterisation();
1687 G4VSolid* pSolid =
1688 pParam->ComputeSolid(enteringReplicaNo, pEnteringPhysVol);
1689 pSolid->ComputeDimensions(pParam, enteringReplicaNo, pEnteringPhysVol);
1690
1691 // Sets the transform in the Parameterisation
1692 //
1693 pParam->ComputeTransformation(enteringReplicaNo, pEnteringPhysVol);
1694
1695 // Set the correct solid and material in Logical Volume
1696 //
1697 G4LogicalVolume* pLogical = pEnteringPhysVol->GetLogicalVolume();
1698 pLogical->SetSolid( pSolid );
1699 }
1700 break;
1701 case kExternal:
1702 G4Exception("G4ITNavigator2::GetMotherToDaughterTransform()",
1703 "GeomNav0001", FatalException,
1704 "Not applicable for external volumes.");
1705 break;
1706 }
1707 return G4AffineTransform(pEnteringPhysVol->GetRotation(),
1708 pEnteringPhysVol->GetTranslation()).Invert();
1709 }
1710
1711 // ********************************************************************
1712 // GetLocalExitNormalAndCheck
1713 //
1714 // Obtains the Normal vector to a surface (in local coordinates)
1715 // pointing out of previous volume and into current volume, and
1716 // checks the current point against expected 'local' value.
1717 // ********************************************************************
1718 //
1719 G4ThreeVector G4ITNavigator2::
1720 GetLocalExitNormalAndCheck(
1721 #ifdef G4DEBUG_NAVIGATION
1722 const G4ThreeVector& ExpectedBoundaryPointGlobal,
1723 #else
1724 const G4ThreeVector&,
1725 #endif
1726 G4bool* pValid)
1727 {
1728 CheckNavigatorStateIsValid();
1729 #ifdef G4DEBUG_NAVIGATION
1730 // Check Current point against expected 'local' value
1731 //
1732 if ( fLastTriedStepComputation )
1733 {
1734 G4ThreeVector ExpectedBoundaryPointLocal;
1735
1736 const G4AffineTransform& GlobalToLocal= GetGlobalToLocalTransform();
1737 ExpectedBoundaryPointLocal =
1738 GlobalToLocal.TransformPoint( ExpectedBoundaryPointGlobal );
1739
1740 // Add here: Comparison against expected position,
1741 // i.e. the endpoint of ComputeStep
1742 }
1743 #endif
1744
1745 return GetLocalExitNormal( pValid);
1746 }
1747
1748 // ********************************************************************
1749 // GetGlobalExitNormal
1750 //
1751 // Obtains the Normal vector to a surface (in global coordinates)
1752 // pointing out of previous volume and into current volume
1753 // ********************************************************************
1754 //
1755 G4ThreeVector
1756 G4ITNavigator2::GetGlobalExitNormal(const G4ThreeVector& IntersectPointGlobal,
1757 G4bool* pNormalCalculated)
1758 {
1759 CheckNavigatorStateIsValid();
1760 G4bool validNormal;
1761 G4ThreeVector localNormal, globalNormal;
1762 G4bool usingStored;
1763
1764 usingStored=
1765 fCalculatedExitNormal &&
1766 ( ( fLastTriedStepComputation && fExiting) // Just calculated it
1767 || // No locate in between
1768 (!fLastTriedStepComputation
1769 && (IntersectPointGlobal-fStepEndPoint).mag2()
1770 < (10.0*kCarTolerance*kCarTolerance)
1771 ) // Calculated it 'just' before & then called locate
1772 // but it did not move position
1773 );
1774
1775 if( usingStored )
1776 {
1777 // This was computed in ComputeStep -- and only on arrival at boundary
1778 //
1779 globalNormal = fExitNormalGlobalFrame;
1780 G4double normMag2 = globalNormal.mag2();
1781 if( std::fabs ( normMag2 - 1.0 ) < perMillion ) // Value is good
1782 {
1783 *pNormalCalculated = true; // ComputeStep always computes it if Exiting
1784 // (fExiting==true)
1785 }
1786 else
1787 {
1788 G4ExceptionDescription message;
1789
1790 message << " ERROR> Expected normal-global-frame to valid (unit vector) "
1791 << " - but |normal| = " << std::sqrt(normMag2)
1792 << " - and |normal|^ = " << normMag2
1793 << " which differs from 1.0 by " << normMag2 - 1.0 << G4endl
1794 << " n = " << fExitNormalGlobalFrame << G4endl;
1795 message << "============================================================"
1796 << G4endl;
1797 G4int oldVerbose = fVerbose;
1798 fVerbose=4;
1799 message << " State of Navigator: " << G4endl;
1800 message << *this << G4endl;
1801 fVerbose = oldVerbose;
1802 message << "============================================================"
1803 << G4endl;
1804
1805 G4Exception("G4ITNavigator2::GetGlobalExitNormal()",
1806 "GeomNav0003",JustWarning, message,
1807 "Value obtained from stored global-normal is not a unit vector.");
1808
1809 // (Re)Compute it now -- as either it was not computed, or it is wrong.
1810 //
1811 localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,
1812 &validNormal);
1813 *pNormalCalculated = fCalculatedExitNormal;
1814
1815 G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
1816 globalNormal = localToGlobal.TransformAxis( localNormal );
1817 }
1818 }
1819 else
1820 {
1821 localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,&validNormal);
1822 *pNormalCalculated = fCalculatedExitNormal;
1823
1824 #ifdef G4DEBUG_NAVIGATION
1825 usingStored= false;
1826
1827 if( (!validNormal) && !fCalculatedExitNormal)
1828 {
1829 G4ExceptionDescription edN;
1830 edN << " Calculated = " << fCalculatedExitNormal << G4endl;
1831 edN << " Entering= " << fEntering << G4endl;
1832 G4int oldVerbose= this->GetVerboseLevel();
1833 this->SetVerboseLevel(4);
1834 edN << " State of Navigator: " << G4endl;
1835 edN << *this << G4endl;
1836 this->SetVerboseLevel( oldVerbose );
1837
1838 G4Exception("G4ITNavigator2::GetGlobalExitNormal()",
1839 "GeomNav0003", JustWarning, edN,
1840 "LocalExitNormalAndCheck() did not calculate Normal.");
1841 }
1842 #endif
1843
1844 G4double localMag2= localNormal.mag2();
1845 if( validNormal && (std::fabs(localMag2-1.0)) > CLHEP::perMillion )
1846 {
1847 G4ExceptionDescription edN;
1848
1849 edN << "G4ITNavigator2::GetGlobalExitNormal: "
1850 << " Using Local Normal - from call to GetLocalExitNormalAndCheck. "
1851 << G4endl
1852 << " Local Exit Normal : " << " || = " << std::sqrt(localMag2)
1853 << " vec = " << localNormal << G4endl
1854 << " Global Exit Normal : " << " || = " << globalNormal.mag()
1855 << " vec = " << globalNormal << G4endl;
1856 edN << " Calculated It = " << fCalculatedExitNormal << G4endl;
1857
1858 G4Exception("G4ITNavigator2::GetGlobalExitNormal()",
1859 "GeomNav0003",JustWarning, edN,
1860 "Value obtained from new local *solid* is incorrect.");
1861 localNormal = localNormal.unit(); // Should we correct it ??
1862 }
1863 G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
1864 globalNormal = localToGlobal.TransformAxis( localNormal );
1865 }
1866
1867 #ifdef G4DEBUG_NAVIGATION
1868 if( usingStored )
1869 {
1870 G4ThreeVector globalNormAgn;
1871
1872 localNormal= GetLocalExitNormalAndCheck(IntersectPointGlobal, &validNormal);
1873
1874 G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
1875 globalNormAgn = localToGlobal.TransformAxis( localNormal );
1876
1877 // Check the value computed against fExitNormalGlobalFrame
1878 G4ThreeVector diffNorm = globalNormAgn - fExitNormalGlobalFrame;
1879 if( diffNorm.mag2() > perMillion*CLHEP::perMillion)
1880 {
1881 G4ExceptionDescription edDfn;
1882 edDfn << "Found difference in normals in case of exiting mother "
1883 << "- when Get is called after ComputingStep " << G4endl;
1884 edDfn << " Magnitude of diff = " << diffNorm.mag() << G4endl;
1885 edDfn << " Normal stored (Global) = " << fExitNormalGlobalFrame
1886 << G4endl;
1887 edDfn << " Global Computed from Local = " << globalNormAgn << G4endl;
1888 G4Exception("G4ITNavigator::GetGlobalExitNormal()", "GeomNav0003",
1889 JustWarning, edDfn);
1890 }
1891 }
1892 #endif
1893
1894 return globalNormal;
1895 }
1896
1897 // To make the new Voxel Safety the default, uncomment the next line
1898 #define G4NEW_SAFETY 1
1899
1900 // ********************************************************************
1901 // ComputeSafety
1902 //
1903 // It assumes that it will be
1904 // i) called at the Point in the same volume as the EndPoint of the
1905 // ComputeStep.
1906 // ii) after (or at the end of) ComputeStep OR after the relocation.
1907 // ********************************************************************
1908 //
1909 G4double G4ITNavigator2::ComputeSafety( const G4ThreeVector &pGlobalpoint,
1910 const G4double pMaxLength,
1911 const G4bool keepState)
1912 {
1913 CheckNavigatorStateIsValid();
1914 G4double newSafety = 0.0;
1915
1916 #ifdef G4DEBUG_NAVIGATION
1917 G4int oldcoutPrec = G4cout.precision(8);
1918 if( fVerbose > 0 )
1919 {
1920 G4cout << "*** G4ITNavigator2::ComputeSafety: ***" << G4endl
1921 << " Called at point: " << pGlobalpoint << G4endl;
1922
1923 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1924 G4cout << " Volume = " << motherPhysical->GetName()
1925 << " - Maximum length = " << pMaxLength << G4endl;
1926 if( fVerbose >= 4 )
1927 {
1928 G4cout << " ----- Upon entering Compute Safety:" << G4endl;
1929 PrintState();
1930 }
1931 }
1932 #endif
1933
1934 G4SaveNavigatorState* savedState(nullptr);
1935
1936 G4double distEndpointSq = (pGlobalpoint-fStepEndPoint).mag2();
1937 G4bool stayedOnEndpoint = distEndpointSq < kCarTolerance*kCarTolerance;
1938 G4bool endpointOnSurface = fEnteredDaughter || fExitedMother;
1939
1940 if( endpointOnSurface && stayedOnEndpoint )
1941 {
1942 #ifdef G4DEBUG_NAVIGATION
1943 if( fVerbose >= 2 )
1944 {
1945 G4cout << " G4Navigator::ComputeSafety() finds that point - "
1946 << pGlobalpoint << " - is on surface " << G4endl;
1947 if( fEnteredDaughter ) { G4cout << " entered new daughter volume"; }
1948 if( fExitedMother ) { G4cout << " and exited previous volume."; }
1949 G4cout << G4endl;
1950 G4cout << " EndPoint was = " << fStepEndPoint << G4endl;
1951 }
1952 #endif
1953 newSafety = 0.0;
1954 // return newSafety;
1955 }
1956 else // if( !(endpointOnSurface && stayedOnEndpoint) )
1957 {
1958
1959 if (keepState)
1960 {
1961 // SetSavedState();
1962 savedState = new G4SaveNavigatorState(fpNavigatorState);
1963 }
1964
1965 // Pseudo-relocate to this point (updates voxel information only)
1966 //
1967 LocateGlobalPointWithinVolume( pGlobalpoint );
1968 // --->> DANGER: Side effects on sub-navigator voxel information <<---
1969 // Could be replaced again by 'granular' calls to sub-navigator
1970 // locates (similar side-effects, but faster.
1971 // Solutions:
1972 // 1) Re-locate (to where?)
1973 // 2) Insure that the methods using (G4ComputeStep?)
1974 // does a relocation (if information is disturbed only ?)
1975
1976 #ifdef G4DEBUG_NAVIGATION
1977 if( fVerbose >= 2 )
1978 {
1979 G4cout << " G4ITNavigator2::ComputeSafety() relocates-in-volume to point: "
1980 << pGlobalpoint << G4endl;
1981 }
1982 #endif
1983 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1984 G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
1985 G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
1986 G4ThreeVector localPoint = ComputeLocalPoint(pGlobalpoint);
1987
1988 if ( fHistory.GetTopVolumeType()!=kReplica )
1989 {
1990 switch(CharacteriseDaughters(motherLogical))
1991 {
1992 case kNormal:
1993 if ( pVoxelHeader != nullptr )
1994 {
1995 #ifdef G4NEW_SAFETY
1996 G4double safetyTwo = fpVoxelSafety->ComputeSafety(localPoint,
1997 *motherPhysical, pMaxLength);
1998 newSafety= safetyTwo; // Faster and best available
1999 #else
2000 G4double safetyOldVoxel;
2001 LocateGlobalPointWithinVolume(pGlobalpoint);
2002 safetyOldVoxel =
2003 fvoxelNav.ComputeSafety(localPoint,fHistory,pMaxLength);
2004 newSafety= safetyOldVoxel;
2005 #endif
2006 }
2007 else
2008 {
2009 newSafety=fnormalNav.ComputeSafety(localPoint,fHistory,pMaxLength);
2010 }
2011 break;
2012 case kParameterised:
2013 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
2014 {
2015 newSafety=fparamNav.ComputeSafety(localPoint,fHistory,pMaxLength);
2016 }
2017 else // Regular structure
2018 {
2019 newSafety=fregularNav.ComputeSafety(localPoint,fHistory,pMaxLength);
2020 }
2021 break;
2022 case kReplica:
2023 G4Exception("G4ITNavigator2::ComputeSafety()", "GeomNav0001",
2024 FatalException, "Not applicable for replicated volumes.");
2025 break;
2026 case kExternal:
2027 G4Exception("G4ITNavigator2::ComputeSafety()",
2028 "GeomNav0001", FatalException,
2029 "Not applicable for external volumes.");
2030 break;
2031 }
2032 }
2033 else
2034 {
2035 newSafety = freplicaNav.ComputeSafety(pGlobalpoint, localPoint,
2036 fHistory, pMaxLength);
2037 }
2038
2039 if (keepState)
2040 {
2041 *fpNavigatorState = *savedState;
2042 delete savedState;
2043 // RestoreSavedState();
2044 // This now overwrites the values of the Safety 'sphere' (correction)
2045 }
2046
2047 // Remember last safety origin & value
2048 //
2049 // We overwrite the Safety 'sphere' - keeping old behaviour
2050 fPreviousSftOrigin = pGlobalpoint;
2051 fPreviousSafety = newSafety;
2052 }
2053
2054 #ifdef G4DEBUG_NAVIGATION
2055 if( fVerbose > 1 )
2056 {
2057 G4cout << " ---- Exiting ComputeSafety " << G4endl;
2058 if( fVerbose > 2 ) { PrintState(); }
2059 G4cout << " Returned value of Safety = " << newSafety << G4endl;
2060 }
2061 G4cout.precision(oldcoutPrec);
2062 #endif
2063
2064 return newSafety;
2065 }
2066
2067
2068 // ********************************************************************
2069 // RecheckDistanceToCurrentBoundary
2070 //
2071 // Trial method for checking potential displacement for MS
2072 // Check position aDisplacedGlobalpoint, to see whether it is in the
2073 // current volume (mother).
2074 // If in mother, check distance to boundary along aNewDirection.
2075 // If in entering daughter, check distance back to boundary.
2076 // NOTE:
2077 // Can be called only after ComputeStep() is called - before ReLocation
2078 // Deals only with current volume (and potentially entered)
2079 // Caveats
2080 // First VERSION: Does not consider daughter volumes if it remained in mother
2081 // neither for checking whether it has exited current (mother) volume,
2082 // nor for determining distance to exit this (mother) volume.
2083 // ********************************************************************
2084 //
2085 G4bool G4ITNavigator2::RecheckDistanceToCurrentBoundary(
2086 const G4ThreeVector &aDisplacedGlobalPoint,
2087 const G4ThreeVector &aNewDirection,
2088 const G4double ProposedMove,
2089 G4double *prDistance,
2090 G4double *prNewSafety) const
2091 {
2092 G4ThreeVector localPosition = ComputeLocalPoint(aDisplacedGlobalPoint);
2093 G4ThreeVector localDirection = ComputeLocalAxis(aNewDirection);
2094 // G4double Step = kInfinity;
2095
2096 G4bool validExitNormal;
2097 G4ThreeVector exitNormal;
2098 // Check against mother solid
2099 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
2100 G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
2101
2102 #ifdef CHECK_ORDER_OF_METHODS
2103 if( ! fLastTriedStepComputation )
2104 {
2105 G4Exception("G4Navigator::RecheckDistanceToCurrentBoundary()",
2106 "GeomNav0001", FatalException,
2107 "Method must be called after ComputeStep(), before call to LocateMethod.");
2108 }
2109 #endif
2110
2111 EInside locatedDaug; // = kUndefined;
2112 G4double daughterStep= DBL_MAX;
2113 G4double daughterSafety= DBL_MAX;
2114
2115 if( fEnteredDaughter )
2116 {
2117 if( motherLogical->CharacteriseDaughters() ==kReplica ) { return false; }
2118
2119 // Track arrived at boundary of a daughter volume at
2120 // the last call of ComputeStep().
2121 // In case the proposed displaced point is inside this daughter,
2122 // it must backtrack at least to the entry point.
2123 // NOTE: No check is made against other daughter volumes. It is
2124 // assumed that the proposed displacement is small enough that
2125 // this is not needed.
2126
2127 // Must check boundary of current daughter
2128 //
2129 G4VPhysicalVolume *candPhysical= fBlockedPhysicalVolume;
2130 G4LogicalVolume *candLogical= candPhysical->GetLogicalVolume();
2131 G4VSolid *candSolid= candLogical->GetSolid();
2132
2133 G4AffineTransform nextLevelTrf(candPhysical->GetRotation(),
2134 candPhysical->GetTranslation());
2135
2136 G4ThreeVector dgPosition= nextLevelTrf.TransformPoint(localPosition);
2137 G4ThreeVector dgDirection= nextLevelTrf.TransformAxis(localDirection);
2138 locatedDaug = candSolid->Inside(dgPosition);
2139
2140 if( locatedDaug == kInside )
2141 {
2142 // Reverse direction - and find first exit. ( Is it valid?)
2143 // Must backtrack
2144 G4double distanceBackOut =
2145 candSolid->DistanceToOut(dgPosition,
2146 - dgDirection, // Reverse direction
2147 true, &validExitNormal, &exitNormal);
2148 daughterStep= - distanceBackOut;
2149 // No check is made whether the particle could have arrived at
2150 // at this location without encountering another volume or
2151 // a different psurface of the current volume
2152 if( prNewSafety != nullptr )
2153 {
2154 daughterSafety= candSolid->DistanceToOut(dgPosition);
2155 }
2156 }
2157 else
2158 {
2159 if( locatedDaug == kOutside )
2160 {
2161 // See whether it still intersects it
2162 //
2163 daughterStep= candSolid->DistanceToIn(dgPosition,
2164 dgDirection);
2165 if( prNewSafety != nullptr )
2166 {
2167 daughterSafety= candSolid->DistanceToIn(dgPosition);
2168 }
2169 }
2170 else
2171 {
2172 // The point remains on the surface of candidate solid
2173 //
2174 daughterStep= 0.0;
2175 daughterSafety= 0.0;
2176 }
2177 }
2178
2179 // If trial point is in daughter (or on its surface) we have the
2180 // answer, the rest is not relevant
2181 //
2182 if( locatedDaug != kOutside )
2183 {
2184 *prDistance= daughterStep;
2185 if( prNewSafety != nullptr ) { *prNewSafety= daughterSafety; }
2186 return true;
2187 }
2188 // If ever extended, so that some type of mother cut daughter,
2189 // this would change
2190 }
2191
2192 G4VSolid *motherSolid= motherLogical->GetSolid();
2193
2194 G4double motherStep= DBL_MAX, motherSafety= DBL_MAX;
2195
2196 // Check distance to boundary of mother
2197 //
2198 if ( fHistory.GetTopVolumeType()!=kReplica )
2199 {
2200 EInside locatedMoth = motherSolid->Inside(localPosition);
2201
2202 if( locatedMoth == kInside )
2203 {
2204 motherSafety= motherSolid->DistanceToOut(localPosition);
2205 if( ProposedMove >= motherSafety )
2206 {
2207 motherStep= motherSolid->DistanceToOut(localPosition,
2208 localDirection,
2209 true, &validExitNormal, &exitNormal);
2210 }
2211 else
2212 {
2213 motherStep= ProposedMove;
2214 }
2215 }
2216 else if( locatedMoth == kOutside)
2217 {
2218 motherSafety= motherSolid->DistanceToIn(localPosition);
2219 if( ProposedMove >= motherSafety )
2220 {
2221 motherStep= - motherSolid->DistanceToIn(localPosition,
2222 -localDirection);
2223 }
2224 }
2225 else
2226 {
2227 motherSafety= 0.0;
2228 *prDistance= 0.0; // On surface - no move // motherStep;
2229 if( prNewSafety != nullptr ) { *prNewSafety= motherSafety; }
2230 return false;
2231 }
2232 }
2233 else
2234 {
2235 return false;
2236 }
2237
2238 *prDistance= std::min( motherStep, daughterStep );
2239 if( prNewSafety != nullptr )
2240 {
2241 *prNewSafety= std::min( motherSafety, daughterSafety );
2242 }
2243 return true;
2244 }
2245
2246
2247 // ********************************************************************
2248 // CreateTouchableHistoryHandle
2249 // ********************************************************************
2250 //
2251 G4TouchableHandle G4ITNavigator2::CreateTouchableHistoryHandle() const
2252 {
2253 CheckNavigatorStateIsValid();
2254 return G4TouchableHandle( CreateTouchableHistory() );
2255 }
2256
2257 // ********************************************************************
2258 // PrintState
2259 // ********************************************************************
2260 //
2261 void G4ITNavigator2::PrintState() const
2262 {
2263 CheckNavigatorStateIsValid();
2264 G4long oldcoutPrec = G4cout.precision(4);
2265 if( fVerbose >= 4 )
2266 {
2267 G4cout << "The current state of G4Navigator is: " << G4endl;
2268 G4cout << " ValidExitNormal= " << fValidExitNormal // << G4endl
2269 << " ExitNormal = " << fExitNormal // << G4endl
2270 << " Exiting = " << fExiting // << G4endl
2271 << " Entering = " << fEntering // << G4endl
2272 << " BlockedPhysicalVolume= " ;
2273 if (fBlockedPhysicalVolume==nullptr)
2274 G4cout << "None";
2275 else
2276 G4cout << fBlockedPhysicalVolume->GetName();
2277 G4cout << G4endl
2278 << " BlockedReplicaNo = " << fBlockedReplicaNo // << G4endl
2279 << " LastStepWasZero = " << fLastStepWasZero // << G4endl
2280 << G4endl;
2281 }
2282 if( ( 1 < fVerbose) && (fVerbose < 4) )
2283 {
2284 G4cout << G4endl; // Make sure to line up
2285 G4cout << std::setw(30) << " ExitNormal " << " "
2286 << std::setw( 5) << " Valid " << " "
2287 << std::setw( 9) << " Exiting " << " "
2288 << std::setw( 9) << " Entering" << " "
2289 << std::setw(15) << " Blocked:Volume " << " "
2290 << std::setw( 9) << " ReplicaNo" << " "
2291 << std::setw( 8) << " LastStepZero " << " "
2292 << G4endl;
2293 G4cout << "( " << std::setw(7) << fExitNormal.x()
2294 << ", " << std::setw(7) << fExitNormal.y()
2295 << ", " << std::setw(7) << fExitNormal.z() << " ) "
2296 << std::setw( 5) << fValidExitNormal << " "
2297 << std::setw( 9) << fExiting << " "
2298 << std::setw( 9) << fEntering << " ";
2299 if ( fBlockedPhysicalVolume==nullptr )
2300 {
2301 G4cout << std::setw(15) << "None";
2302 }
2303 else
2304 {
2305 G4cout << std::setw(15)<< fBlockedPhysicalVolume->GetName();
2306 }
2307 G4cout << std::setw( 9) << fBlockedReplicaNo << " "
2308 << std::setw( 8) << fLastStepWasZero << " "
2309 << G4endl;
2310 }
2311 if( fVerbose > 2 )
2312 {
2313 G4cout.precision(8);
2314 G4cout << " Current Localpoint = " << fLastLocatedPointLocal << G4endl;
2315 G4cout << " PreviousSftOrigin = " << fPreviousSftOrigin << G4endl;
2316 G4cout << " PreviousSafety = " << fPreviousSafety << G4endl;
2317 }
2318 G4cout.precision(oldcoutPrec);
2319 }
2320
2321 // ********************************************************************
2322 // ComputeStepLog
2323 // ********************************************************************
2324 //
2325 void G4ITNavigator2::ComputeStepLog(const G4ThreeVector& pGlobalpoint,
2326 G4double moveLenSq) const
2327 {
2328 CheckNavigatorStateIsValid();
2329 // The following checks only make sense if the move is larger
2330 // than the tolerance.
2331
2332 static const G4double fAccuracyForWarning = kCarTolerance,
2333 fAccuracyForException = 1000*kCarTolerance;
2334
2335 G4ThreeVector OriginalGlobalpoint = fHistory.GetTopTransform().Inverse().
2336 TransformPoint(fLastLocatedPointLocal);
2337
2338 G4double shiftOriginSafSq = (fPreviousSftOrigin-pGlobalpoint).mag2();
2339
2340 // Check that the starting point of this step is
2341 // within the isotropic safety sphere of the last point
2342 // to a accuracy/precision given by fAccuracyForWarning.
2343 // If so give warning.
2344 // If it fails by more than fAccuracyForException exit with error.
2345 //
2346 if( shiftOriginSafSq >= sqr(fPreviousSafety) )
2347 {
2348 G4double shiftOrigin = std::sqrt(shiftOriginSafSq);
2349 G4double diffShiftSaf = shiftOrigin - fPreviousSafety;
2350
2351 if( diffShiftSaf > fAccuracyForWarning )
2352 {
2353 G4long oldcoutPrec= G4cout.precision(8);
2354 G4long oldcerrPrec= G4cerr.precision(10);
2355 std::ostringstream message, suggestion;
2356 message << "Accuracy error or slightly inaccurate position shift."
2357 << G4endl
2358 << " The Step's starting point has moved "
2359 << std::sqrt(moveLenSq)/mm << " mm " << G4endl
2360 << " since the last call to a Locate method." << G4endl
2361 << " This has resulted in moving "
2362 << shiftOrigin/mm << " mm "
2363 << " from the last point at which the safety "
2364 << " was calculated " << G4endl
2365 << " which is more than the computed safety= "
2366 << fPreviousSafety/mm << " mm at that point." << G4endl
2367 << " This difference is "
2368 << diffShiftSaf/mm << " mm." << G4endl
2369 << " The tolerated accuracy is "
2370 << fAccuracyForException/mm << " mm.";
2371
2372 suggestion << " ";
2373 static G4ThreadLocal G4int warnNow = 0;
2374 if( ((++warnNow % 100) == 1) )
2375 {
2376 message << G4endl
2377 << " This problem can be due to either " << G4endl
2378 << " - a process that has proposed a displacement"
2379 << " larger than the current safety , or" << G4endl
2380 << " - inaccuracy in the computation of the safety";
2381 suggestion << "We suggest that you " << G4endl
2382 << " - find i) what particle is being tracked, and "
2383 << " ii) through what part of your geometry " << G4endl
2384 << " for example by re-running this event with "
2385 << G4endl
2386 << " /tracking/verbose 1 " << G4endl
2387 << " - check which processes you declare for"
2388 << " this particle (and look at non-standard ones)"
2389 << G4endl
2390 << " - in case, create a detailed logfile"
2391 << " of this event using:" << G4endl
2392 << " /tracking/verbose 6 ";
2393 }
2394 G4Exception("G4ITNavigator2::ComputeStep()",
2395 "GeomNav1002", JustWarning,
2396 message, G4String(suggestion.str()));
2397 G4cout.precision(oldcoutPrec);
2398 G4cerr.precision(oldcerrPrec);
2399 }
2400 #ifdef G4DEBUG_NAVIGATION
2401 else
2402 {
2403 G4cerr << "WARNING - G4ITNavigator2::ComputeStep()" << G4endl
2404 << " The Step's starting point has moved "
2405 << std::sqrt(moveLenSq) << "," << G4endl
2406 << " which has taken it to the limit of"
2407 << " the current safety. " << G4endl;
2408 }
2409 #endif
2410 }
2411 G4double safetyPlus = fPreviousSafety + fAccuracyForException;
2412 if ( shiftOriginSafSq > sqr(safetyPlus) )
2413 {
2414 std::ostringstream message;
2415 message << "May lead to a crash or unreliable results." << G4endl
2416 << " Position has shifted considerably without"
2417 << " notifying the navigator !" << G4endl
2418 << " Tolerated safety: " << safetyPlus << G4endl
2419 << " Computed shift : " << shiftOriginSafSq;
2420 G4Exception("G4ITNavigator2::ComputeStep()", "GeomNav1002",
2421 JustWarning, message);
2422 }
2423 }
2424
2425 // ********************************************************************
2426 // Operator <<
2427 // ********************************************************************
2428 //
2429 std::ostream& operator << (std::ostream &os,const G4ITNavigator2 &n)
2430 {
2431 // Old version did only the following:
2432 // os << "Current History: " << G4endl << n.fHistory;
2433 // Old behaviour is recovered for fVerbose = 0
2434
2435 // Adapted from G4ITNavigator2::PrintState() const
2436
2437 G4long oldcoutPrec = os.precision(4);
2438 if( n.fVerbose >= 4 )
2439 {
2440 os << "The current state of G4ITNavigator2 is: " << G4endl;
2441 os << " ValidExitNormal= " << n.fValidExitNormal << G4endl
2442 << " ExitNormal = " << n.fExitNormal << G4endl
2443 << " Exiting = " << n.fExiting << G4endl
2444 << " Entering = " << n.fEntering << G4endl
2445 << " BlockedPhysicalVolume= " ;
2446
2447 if (n.fBlockedPhysicalVolume==nullptr)
2448 {
2449 os << "None";
2450 }
2451 else
2452 {
2453 os << n.fBlockedPhysicalVolume->GetName();
2454 }
2455
2456 os << G4endl
2457 << " BlockedReplicaNo = " << n.fBlockedReplicaNo << G4endl
2458 << " LastStepWasZero = " << n.fLastStepWasZero << G4endl
2459 << G4endl;
2460 }
2461 if( ( 1 < n.fVerbose) && (n.fVerbose < 4) )
2462 {
2463 os << G4endl; // Make sure to line up
2464 os << std::setw(30) << " ExitNormal " << " "
2465 << std::setw( 5) << " Valid " << " "
2466 << std::setw( 9) << " Exiting " << " "
2467 << std::setw( 9) << " Entering" << " "
2468 << std::setw(15) << " Blocked:Volume " << " "
2469 << std::setw( 9) << " ReplicaNo" << " "
2470 << std::setw( 8) << " LastStepZero " << " "
2471 << G4endl;
2472 os << "( " << std::setw(7) << n.fExitNormal.x()
2473 << ", " << std::setw(7) << n.fExitNormal.y()
2474 << ", " << std::setw(7) << n.fExitNormal.z() << " ) "
2475 << std::setw( 5) << n.fValidExitNormal << " "
2476 << std::setw( 9) << n.fExiting << " "
2477 << std::setw( 9) << n.fEntering << " ";
2478
2479 if ( n.fBlockedPhysicalVolume==nullptr )
2480 { os << std::setw(15) << "None"; }
2481 else
2482 { os << std::setw(15)<< n.fBlockedPhysicalVolume->GetName(); }
2483
2484 os << std::setw( 9) << n.fBlockedReplicaNo << " "
2485 << std::setw( 8) << n.fLastStepWasZero << " "
2486 << G4endl;
2487 }
2488 if( n.fVerbose > 2 )
2489 {
2490 os.precision(8);
2491 os << " Current Localpoint = " << n.fLastLocatedPointLocal << G4endl;
2492 os << " PreviousSftOrigin = " << n.fPreviousSftOrigin << G4endl;
2493 os << " PreviousSafety = " << n.fPreviousSafety << G4endl;
2494 }
2495 if( n.fVerbose > 3 || n.fVerbose == 0 )
2496 {
2497 os << "Current History: " << G4endl << n.fHistory;
2498 }
2499
2500 os.precision(oldcoutPrec);
2501 return os;
2502 }
2503
2504 //------------------------------------------------------------------------------
2505
2506 EInside
2507 G4ITNavigator2::InsideCurrentVolume(const G4ThreeVector& globalPoint) const
2508 {
2509 const G4AffineTransform& transform = GetGlobalToLocalTransform();
2510 G4ThreeVector localPoint(transform.TransformPoint(globalPoint));
2511
2512 G4VSolid* solid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid();
2513
2514 return solid->Inside(localPoint);
2515 }
2516
2517 //------------------------------------------------------------------------------
2518
2519 void G4ITNavigator2::GetRandomInCurrentVolume(G4ThreeVector& _rndmPoint) const
2520 {
2521 const G4AffineTransform& local2Global = GetLocalToGlobalTransform();
2522 G4VSolid* solid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid();
2523
2524 G4VoxelLimits voxelLimits; // Defaults to "infinite" limits.
2525 G4double vmin, vmax;
2526 G4AffineTransform dummy;
2527 std::vector<std::vector<G4double> > fExtend;
2528
2529 solid->CalculateExtent(kXAxis,voxelLimits,dummy,vmin,vmax);
2530 fExtend[kXAxis][BoundingBox::kMin] = vmin;
2531 fExtend[kXAxis][BoundingBox::kMax] = vmax;
2532
2533 solid->CalculateExtent(kYAxis,voxelLimits,dummy,vmin,vmax);
2534 fExtend[kYAxis][BoundingBox::kMin] = vmin;
2535 fExtend[kYAxis][BoundingBox::kMax] = vmax;
2536
2537 solid->CalculateExtent(kZAxis,voxelLimits,dummy,vmin,vmax);
2538 fExtend[kZAxis][BoundingBox::kMin] = vmin;
2539 fExtend[kZAxis][BoundingBox::kMax] = vmax;
2540
2541 G4ThreeVector rndmPos;
2542
2543 while(true)
2544 {
2545 for(G4int i = 0 ; i < 3 ; ++i)
2546 {
2547 G4double min = fExtend[i][BoundingBox::kMin];
2548 G4double max = fExtend[i][BoundingBox::kMax];
2549 rndmPos[i] = G4UniformRand()*(max-min)+min;
2550 }
2551
2552 if(solid->Inside(rndmPos) == kInside) break;
2553 }
2554
2555 _rndmPoint = local2Global.TransformPoint(rndmPos);
2556 }
2557
2558