| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4SafetyCalculator class Implementation
27 //
28 // Author: John Apostolakis, CERN - February 2
29 // -------------------------------------------
30
31 #include "G4SafetyCalculator.hh"
32
33 #include "G4Navigator.hh"
34 #include "G4GeometryTolerance.hh"
35
36 G4SafetyCalculator::
37 G4SafetyCalculator( const G4Navigator& navigat
38 const G4NavigationHistory&
39 : fNavigator(navigator), fNavHistory(navHist
40 {
41 fkCarTolerance = G4GeometryTolerance::GetIns
42 }
43
44 // *******************************************
45 // ComputeSafety
46 //
47 // It assumes that it will be
48 // i) called at the Point in the same volume
49 // ComputeStep.
50 // ii) after (or at the end of) ComputeStep OR
51 // *******************************************
52 //
53 G4double G4SafetyCalculator::
54 SafetyInCurrentVolume( const G4ThreeVector& pG
55 G4VPhysicalVolume
56 const G4double pMaxLeng
57 G4bool /* verbose
58 {
59 G4double safety = 0.0;
60 G4ThreeVector stepEndPoint = fNavigator.GetL
61
62 G4ThreeVector localPoint = ComputeLocalPoint
63
64 G4double distEndpointSq = (pGlobalpoint-step
65 G4bool stayedOnEndpoint = distEndpointSq < s
66 G4bool endpointOnSurface = fNavigator.Ente
67 || fNavigator.Exit
68
69 G4VPhysicalVolume* motherPhysical = fNavHist
70 if( motherPhysical != physicalVolume )
71 {
72 std::ostringstream msg;
73 msg << " Current (navigation) phys-volume
74 << " name= " << motherPhysical->GetNa
75 << " Request made for phys-volume
76 << " name= " << physicalVolume->GetNa
77 G4Exception("G4SafetyCalculator::SafetyIn
78 FatalException, msg,
79 "This method must be called o
80 }
81
82 if( !(endpointOnSurface && stayedOnEndpoint)
83 {
84 G4LogicalVolume* motherLogical = motherPhy
85 G4SmartVoxelHeader* pVoxelHeader = motherL
86
87 // Pseudo-relocate to this point (updates
88 //
89 QuickLocateWithinVolume( localPoint, mothe
90 //*********************
91
92 // switch(CharacteriseDaughters(motherLogi
93 auto dtype= CharacteriseDaughters(motherLo
94 switch(dtype)
95 {
96 case kNormal:
97 if ( pVoxelHeader )
98 {
99 // New way: best safety
100 safety = fVoxelSafety.ComputeSafety(
101
102 }
103 else
104 {
105 safety=fnormalNav.ComputeSafety(loca
106 }
107 break;
108 case kParameterised:
109 if( GetDaughtersRegularStructureId(mot
110 {
111 safety=fparamNav.ComputeSafety(local
112 }
113 else // Regular structure
114 {
115 safety=fregularNav.ComputeSafety(loc
116 }
117 break;
118 case kReplica:
119 safety = freplicaNav.ComputeSafety(pGl
120 fNa
121 break;
122 case kExternal:
123 safety = fpExternalNav->ComputeSafety(
124
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
141 // is updated in order to start the next st
142 // -> no check is performed whether pGlobalpoi
143 // original volume (this must be the case).
144 //
145 // Note: a direction could be added to the arg
146 // optional checking (via the old code b
147 // [ This would be done only in verbose
148 //
149 // Adapted simplied from G4Navigator::LocateGl
150 // *******************************************
151 //
152 void G4SafetyCalculator::
153 QuickLocateWithinVolume( const G4ThreeVector&
154 G4VPhysicalVolu
155 {
156 // For the case of Voxel (or Parameterised)
157 // sub-navigator must be messaged to update
158
159 // Update the state of the Sub Navigators
160 // - in particular any voxel information the
161 //
162 G4LogicalVolume* motherLogical = motherP
163 G4SmartVoxelHeader* pVoxelHeader = motherL
164
165 switch( CharacteriseDaughters(motherLogical)
166 {
167 case kNormal:
168 if ( pVoxelHeader )
169 {
170 fvoxelNav.VoxelLocate( pVoxelHeader, p
171 }
172 break;
173 case kParameterised:
174 if( GetDaughtersRegularStructureId(mothe
175 {
176 // Resets state & returns voxel node
177 //
178 fparamNav.ParamVoxelLocate( pVoxelHead
179 }
180 break;
181 case kReplica:
182 // Nothing to do
183 break;
184 case kExternal:
185 fpExternalNav->RelocateWithinVolume( mot
186 poi
187 break;
188 }
189 }
190
191 // *******************************************
192 // Accessor for custom external navigation.
193 // *******************************************
194 G4VExternalNavigation* G4SafetyCalculator::Get
195 {
196 return fpExternalNav;
197 }
198
199 // *******************************************
200 // Modifier for custom external navigation.
201 // *******************************************
202 void G4SafetyCalculator::SetExternalNavigation
203 {
204 fpExternalNav = eNav;
205 }
206
207 // *******************************************
208 // CompareSafetyValues
209 // *******************************************
210 void G4SafetyCalculator::
211 CompareSafetyValues( G4double oldSafety,
212 G4double newValue,
213 G4VPhysicalVolume* mother
214 const G4ThreeVector& globalPoin
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(old
228 if( std::fabs( newValue - oldSafety) > repor
229 {
230 G4ExceptionSeverity severity= FatalExcepti
231 std::ostringstream msg;
232 G4double diff= (newValue-oldSafety);
233 G4double relativeDiff= diff / oldSafetyPlu
234
235 msg << " New (G4SafetyCalculator) value *d
236 << " in physical volume '" << motherPh
237 << "copy-no = " << motherPhysical->Get
238 if( enteredDaughterVol ) { msg << " ( Jus
239 if( exitedMotherVol ) { msg << " ( Jus
240 msg << G4endl;
241 msg << " Safeties: old= " << std::setpr
242 << " trial " << newValue
243 << " new-old= " << std::setprecision(
244
245 if( std::fabs(diff) < errorThreshold * ( s
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 = " << keepSt
257 msg << " Location - Global coordinates:
258 << " volume= '" << motherPhysical->
259 << " copy-no= " << motherPhysical->G
260 msg << " Argument maxLength= " << maxLen
261
262 std::size_t depth= fNavHistory.GetDepth(
263 msg << " Navigation History: depth = " <
264 for( G4int i=1; i<(G4int)depth; ++i )
265 {
266 msg << " d= " << i << " " << std:
267 << fNavHistory.GetVolume(i)->GetN
268 << " copyNo= " << fNavHistory.Ge
269 msg << G4endl;
270 }
271 }
272
273 #ifdef G4DEBUG_NAVIGATION
274 G4double redo= SafetyInCurrentVolume(globa
275 maxLe
276 msg << " Redoing estimator: value = " << s
277 << " diff/last= " << std::setprecisio
278 << " diff/old= " << redo - oldSafety
279 #endif
280
281 G4Exception("G4SafetyCalculator::CompareSa
282 severity, msg);
283 }
284 }
285