| 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 // G4SafetyCalculator class Implementation
27 //
28 // Author: John Apostolakis, CERN - February 2023
29 // --------------------------------------------------------------------
30
31 #include "G4SafetyCalculator.hh"
32
33 #include "G4Navigator.hh"
34 #include "G4GeometryTolerance.hh"
35
36 G4SafetyCalculator::
37 G4SafetyCalculator( const G4Navigator& navigator,
38 const G4NavigationHistory& navHistory )
39 : fNavigator(navigator), fNavHistory(navHistory)
40 {
41 fkCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
42 }
43
44 // ********************************************************************
45 // ComputeSafety
46 //
47 // It assumes that it will be
48 // i) called at the Point in the same volume as the EndPoint of the
49 // ComputeStep.
50 // ii) after (or at the end of) ComputeStep OR after the relocation.
51 // ********************************************************************
52 //
53 G4double G4SafetyCalculator::
54 SafetyInCurrentVolume( const G4ThreeVector& pGlobalpoint,
55 G4VPhysicalVolume* physicalVolume,
56 const G4double pMaxLength,
57 G4bool /* verbose */ )
58 {
59 G4double safety = 0.0;
60 G4ThreeVector stepEndPoint = fNavigator.GetLastStepEndPoint();
61
62 G4ThreeVector localPoint = ComputeLocalPoint(pGlobalpoint);
63
64 G4double distEndpointSq = (pGlobalpoint-stepEndPoint).mag2();
65 G4bool stayedOnEndpoint = distEndpointSq < sqr(fkCarTolerance);
66 G4bool endpointOnSurface = fNavigator.EnteredDaughterVolume()
67 || fNavigator.ExitedMotherVolume();
68
69 G4VPhysicalVolume* motherPhysical = fNavHistory.GetTopVolume();
70 if( motherPhysical != physicalVolume )
71 {
72 std::ostringstream msg;
73 msg << " Current (navigation) phys-volume: " << motherPhysical
74 << " name= " << motherPhysical->GetName() << G4endl
75 << " Request made for phys-volume: " << physicalVolume
76 << " name= " << physicalVolume->GetName() << G4endl;
77 G4Exception("G4SafetyCalculator::SafetyInCurrentVolume", "GeomNav0001",
78 FatalException, msg,
79 "This method must be called only in the Current volume.");
80 }
81
82 if( !(endpointOnSurface && stayedOnEndpoint) )
83 {
84 G4LogicalVolume* motherLogical = motherPhysical->GetLogicalVolume();
85 G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
86
87 // Pseudo-relocate to this point (updates voxel information only)
88 //
89 QuickLocateWithinVolume( localPoint, motherPhysical );
90 //*********************
91
92 // switch(CharacteriseDaughters(motherLogical))
93 auto dtype= CharacteriseDaughters(motherLogical);
94 switch(dtype)
95 {
96 case kNormal:
97 if ( pVoxelHeader )
98 {
99 // New way: best safety
100 safety = fVoxelSafety.ComputeSafety(localPoint,
101 *motherPhysical, pMaxLength);
102 }
103 else
104 {
105 safety=fnormalNav.ComputeSafety(localPoint,fNavHistory,pMaxLength);
106 }
107 break;
108 case kParameterised:
109 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
110 {
111 safety=fparamNav.ComputeSafety(localPoint,fNavHistory,pMaxLength);
112 }
113 else // Regular structure
114 {
115 safety=fregularNav.ComputeSafety(localPoint,fNavHistory,pMaxLength);
116 }
117 break;
118 case kReplica:
119 safety = freplicaNav.ComputeSafety(pGlobalpoint, localPoint,
120 fNavHistory, pMaxLength);
121 break;
122 case kExternal:
123 safety = fpExternalNav->ComputeSafety(localPoint, fNavHistory,
124 pMaxLength);
125 break;
126 }
127
128 // Remember last safety origin & value
129 //
130 fPreviousSftOrigin = pGlobalpoint;
131 fPreviousSafety = safety;
132 }
133
134 return safety;
135 }
136
137 // ********************************************************************
138 // QuickLocateWithinVolume
139 //
140 // -> the state information of this Navigator and its subNavigators
141 // is updated in order to start the next step at pGlobalpoint
142 // -> no check is performed whether pGlobalpoint is inside the
143 // original volume (this must be the case).
144 //
145 // Note: a direction could be added to the arguments, to aid in future
146 // optional checking (via the old code below, flagged by OLD_LOCATE).
147 // [ This would be done only in verbose mode ]
148 //
149 // Adapted simplied from G4Navigator::LocateGlobalPointWithinVolume()
150 // ********************************************************************
151 //
152 void G4SafetyCalculator::
153 QuickLocateWithinVolume( const G4ThreeVector& pointLocal,
154 G4VPhysicalVolume* motherPhysical )
155 {
156 // For the case of Voxel (or Parameterised) volume the respective
157 // sub-navigator must be messaged to update its voxel information etc
158
159 // Update the state of the Sub Navigators
160 // - in particular any voxel information they store/cache
161 //
162 G4LogicalVolume* motherLogical = motherPhysical->GetLogicalVolume();
163 G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
164
165 switch( CharacteriseDaughters(motherLogical) )
166 {
167 case kNormal:
168 if ( pVoxelHeader )
169 {
170 fvoxelNav.VoxelLocate( pVoxelHeader, pointLocal );
171 }
172 break;
173 case kParameterised:
174 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
175 {
176 // Resets state & returns voxel node
177 //
178 fparamNav.ParamVoxelLocate( pVoxelHeader, pointLocal );
179 }
180 break;
181 case kReplica:
182 // Nothing to do
183 break;
184 case kExternal:
185 fpExternalNav->RelocateWithinVolume( motherPhysical,
186 pointLocal );
187 break;
188 }
189 }
190
191 // ********************************************************************
192 // Accessor for custom external navigation.
193 // ********************************************************************
194 G4VExternalNavigation* G4SafetyCalculator::GetExternalNavigation() const
195 {
196 return fpExternalNav;
197 }
198
199 // ********************************************************************
200 // Modifier for custom external navigation.
201 // ********************************************************************
202 void G4SafetyCalculator::SetExternalNavigation(G4VExternalNavigation* eNav)
203 {
204 fpExternalNav = eNav;
205 }
206
207 // ********************************************************************
208 // CompareSafetyValues
209 // ********************************************************************
210 void G4SafetyCalculator::
211 CompareSafetyValues( G4double oldSafety,
212 G4double newValue,
213 G4VPhysicalVolume* motherPhysical,
214 const G4ThreeVector& globalPoint,
215 G4bool keepState,
216 G4double maxLength,
217 G4bool enteredDaughterVol,
218 G4bool exitedMotherVol )
219 {
220 constexpr G4double reportThreshold= 3.0e-14;
221 // At least warn if rel-error exceeds it
222 constexpr G4double errorThreshold= 1.0e-08;
223 // Fatal if relative error is larger
224 constexpr G4double epsilonLen= 1.0e-20;
225 // Baseline minimal value for divisor
226
227 const G4double oldSafetyPlus = std::fabs(oldSafety)+epsilonLen;
228 if( std::fabs( newValue - oldSafety) > reportThreshold * oldSafetyPlus )
229 {
230 G4ExceptionSeverity severity= FatalException;
231 std::ostringstream msg;
232 G4double diff= (newValue-oldSafety);
233 G4double relativeDiff= diff / oldSafetyPlus;
234
235 msg << " New (G4SafetyCalculator) value *disagrees* by relative diff " << relativeDiff
236 << " in physical volume '" << motherPhysical->GetName() << "' "
237 << "copy-no = " << motherPhysical->GetCopyNo();
238 if( enteredDaughterVol ) { msg << " ( Just Entered new daughter volume. ) "; }
239 if( exitedMotherVol ) { msg << " ( Just Exited previous volume. ) "; }
240 msg << G4endl;
241 msg << " Safeties: old= " << std::setprecision(12) << oldSafety
242 << " trial " << newValue
243 << " new-old= " << std::setprecision(7) << diff << G4endl;
244
245 if( std::fabs(diff) < errorThreshold * ( std::fabs(oldSafety)+1.0e-20 ) )
246 {
247 msg << " (tiny difference) ";
248 severity= JustWarning;
249 }
250 else
251 {
252 msg << " (real difference) ";
253 severity= FatalException;
254
255 // Extra information -- for big errors
256 msg << " NOTE: keepState = " << keepState << G4endl;
257 msg << " Location - Global coordinates: " << globalPoint
258 << " volume= '" << motherPhysical->GetName() << "'"
259 << " copy-no= " << motherPhysical->GetCopyNo() << G4endl;
260 msg << " Argument maxLength= " << maxLength << G4endl;
261
262 std::size_t depth= fNavHistory.GetDepth();
263 msg << " Navigation History: depth = " << depth << G4endl;
264 for( G4int i=1; i<(G4int)depth; ++i )
265 {
266 msg << " d= " << i << " " << std::setw(32)
267 << fNavHistory.GetVolume(i)->GetName()
268 << " copyNo= " << fNavHistory.GetReplicaNo(i);
269 msg << G4endl;
270 }
271 }
272
273 #ifdef G4DEBUG_NAVIGATION
274 G4double redo= SafetyInCurrentVolume(globalPoint, motherPhysical,
275 maxLength, true);
276 msg << " Redoing estimator: value = " << std::setprecision(16) << redo
277 << " diff/last= " << std::setprecision(7) << redo - newValue
278 << " diff/old= " << redo - oldSafety << G4endl;
279 #endif
280
281 G4Exception("G4SafetyCalculator::CompareSafetyValues()", "GeomNav1007",
282 severity, msg);
283 }
284 }
285