| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4SafetyCalculator class Implementation 26 // G4SafetyCalculator class Implementation
27 // 27 //
28 // Author: John Apostolakis, CERN - February 2 28 // Author: John Apostolakis, CERN - February 2023
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "G4SafetyCalculator.hh" 31 #include "G4SafetyCalculator.hh"
32 32
33 #include "G4Navigator.hh" 33 #include "G4Navigator.hh"
34 #include "G4GeometryTolerance.hh" 34 #include "G4GeometryTolerance.hh"
35 35
36 G4SafetyCalculator:: 36 G4SafetyCalculator::
37 G4SafetyCalculator( const G4Navigator& navigat 37 G4SafetyCalculator( const G4Navigator& navigator,
38 const G4NavigationHistory& 38 const G4NavigationHistory& navHistory )
39 : fNavigator(navigator), fNavHistory(navHist 39 : fNavigator(navigator), fNavHistory(navHistory)
40 { 40 {
41 fkCarTolerance = G4GeometryTolerance::GetIns 41 fkCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
42 } 42 }
43 43
44 // ******************************************* 44 // ********************************************************************
45 // ComputeSafety 45 // ComputeSafety
46 // 46 //
47 // It assumes that it will be 47 // It assumes that it will be
48 // i) called at the Point in the same volume 48 // i) called at the Point in the same volume as the EndPoint of the
49 // ComputeStep. 49 // ComputeStep.
50 // ii) after (or at the end of) ComputeStep OR 50 // ii) after (or at the end of) ComputeStep OR after the relocation.
51 // ******************************************* 51 // ********************************************************************
52 // 52 //
53 G4double G4SafetyCalculator:: 53 G4double G4SafetyCalculator::
54 SafetyInCurrentVolume( const G4ThreeVector& pG 54 SafetyInCurrentVolume( const G4ThreeVector& pGlobalpoint,
55 G4VPhysicalVolume 55 G4VPhysicalVolume* physicalVolume,
56 const G4double pMaxLeng 56 const G4double pMaxLength,
57 G4bool /* verbose 57 G4bool /* verbose */ )
58 { 58 {
59 G4double safety = 0.0; 59 G4double safety = 0.0;
60 G4ThreeVector stepEndPoint = fNavigator.GetL 60 G4ThreeVector stepEndPoint = fNavigator.GetLastStepEndPoint();
61 61
62 G4ThreeVector localPoint = ComputeLocalPoint 62 G4ThreeVector localPoint = ComputeLocalPoint(pGlobalpoint);
63 63
64 G4double distEndpointSq = (pGlobalpoint-step 64 G4double distEndpointSq = (pGlobalpoint-stepEndPoint).mag2();
65 G4bool stayedOnEndpoint = distEndpointSq < s 65 G4bool stayedOnEndpoint = distEndpointSq < sqr(fkCarTolerance);
66 G4bool endpointOnSurface = fNavigator.Ente 66 G4bool endpointOnSurface = fNavigator.EnteredDaughterVolume()
67 || fNavigator.Exit 67 || fNavigator.ExitedMotherVolume();
68 68
69 G4VPhysicalVolume* motherPhysical = fNavHist 69 G4VPhysicalVolume* motherPhysical = fNavHistory.GetTopVolume();
70 if( motherPhysical != physicalVolume ) 70 if( motherPhysical != physicalVolume )
71 { 71 {
72 std::ostringstream msg; 72 std::ostringstream msg;
73 msg << " Current (navigation) phys-volume 73 msg << " Current (navigation) phys-volume: " << motherPhysical
74 << " name= " << motherPhysical->GetNa 74 << " name= " << motherPhysical->GetName() << G4endl
75 << " Request made for phys-volume 75 << " Request made for phys-volume: " << physicalVolume
76 << " name= " << physicalVolume->GetNa 76 << " name= " << physicalVolume->GetName() << G4endl;
77 G4Exception("G4SafetyCalculator::SafetyIn 77 G4Exception("G4SafetyCalculator::SafetyInCurrentVolume", "GeomNav0001",
78 FatalException, msg, 78 FatalException, msg,
79 "This method must be called o 79 "This method must be called only in the Current volume.");
80 } 80 }
81 81
82 if( !(endpointOnSurface && stayedOnEndpoint) 82 if( !(endpointOnSurface && stayedOnEndpoint) )
83 { 83 {
84 G4LogicalVolume* motherLogical = motherPhy 84 G4LogicalVolume* motherLogical = motherPhysical->GetLogicalVolume();
85 G4SmartVoxelHeader* pVoxelHeader = motherL 85 G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
86 86
87 // Pseudo-relocate to this point (updates 87 // Pseudo-relocate to this point (updates voxel information only)
88 // 88 //
89 QuickLocateWithinVolume( localPoint, mothe 89 QuickLocateWithinVolume( localPoint, motherPhysical );
90 //********************* 90 //*********************
91 91
92 // switch(CharacteriseDaughters(motherLogi 92 // switch(CharacteriseDaughters(motherLogical))
93 auto dtype= CharacteriseDaughters(motherLo 93 auto dtype= CharacteriseDaughters(motherLogical);
94 switch(dtype) 94 switch(dtype)
95 { 95 {
96 case kNormal: 96 case kNormal:
97 if ( pVoxelHeader ) 97 if ( pVoxelHeader )
98 { 98 {
99 // New way: best safety 99 // New way: best safety
100 safety = fVoxelSafety.ComputeSafety( 100 safety = fVoxelSafety.ComputeSafety(localPoint,
101 101 *motherPhysical, pMaxLength);
102 } 102 }
103 else 103 else
104 { 104 {
105 safety=fnormalNav.ComputeSafety(loca 105 safety=fnormalNav.ComputeSafety(localPoint,fNavHistory,pMaxLength);
106 } 106 }
107 break; 107 break;
108 case kParameterised: 108 case kParameterised:
109 if( GetDaughtersRegularStructureId(mot 109 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
110 { 110 {
111 safety=fparamNav.ComputeSafety(local 111 safety=fparamNav.ComputeSafety(localPoint,fNavHistory,pMaxLength);
112 } 112 }
113 else // Regular structure 113 else // Regular structure
114 { 114 {
115 safety=fregularNav.ComputeSafety(loc 115 safety=fregularNav.ComputeSafety(localPoint,fNavHistory,pMaxLength);
116 } 116 }
117 break; 117 break;
118 case kReplica: 118 case kReplica:
119 safety = freplicaNav.ComputeSafety(pGl 119 safety = freplicaNav.ComputeSafety(pGlobalpoint, localPoint,
120 fNa 120 fNavHistory, pMaxLength);
121 break; 121 break;
122 case kExternal: 122 case kExternal:
123 safety = fpExternalNav->ComputeSafety( 123 safety = fpExternalNav->ComputeSafety(localPoint, fNavHistory,
124 124 pMaxLength);
125 break; 125 break;
126 } 126 }
127 127
128 // Remember last safety origin & value 128 // Remember last safety origin & value
129 // 129 //
130 fPreviousSftOrigin = pGlobalpoint; 130 fPreviousSftOrigin = pGlobalpoint;
131 fPreviousSafety = safety; 131 fPreviousSafety = safety;
132 } 132 }
133 133
134 return safety; 134 return safety;
135 } 135 }
136 136
137 // ******************************************* 137 // ********************************************************************
138 // QuickLocateWithinVolume 138 // QuickLocateWithinVolume
139 // 139 //
140 // -> the state information of this Navigator 140 // -> the state information of this Navigator and its subNavigators
141 // is updated in order to start the next st 141 // is updated in order to start the next step at pGlobalpoint
142 // -> no check is performed whether pGlobalpoi 142 // -> no check is performed whether pGlobalpoint is inside the
143 // original volume (this must be the case). 143 // original volume (this must be the case).
144 // 144 //
145 // Note: a direction could be added to the arg 145 // Note: a direction could be added to the arguments, to aid in future
146 // optional checking (via the old code b 146 // optional checking (via the old code below, flagged by OLD_LOCATE).
147 // [ This would be done only in verbose 147 // [ This would be done only in verbose mode ]
148 // 148 //
149 // Adapted simplied from G4Navigator::LocateGl 149 // Adapted simplied from G4Navigator::LocateGlobalPointWithinVolume()
150 // ******************************************* 150 // ********************************************************************
151 // 151 //
152 void G4SafetyCalculator:: 152 void G4SafetyCalculator::
153 QuickLocateWithinVolume( const G4ThreeVector& 153 QuickLocateWithinVolume( const G4ThreeVector& pointLocal,
154 G4VPhysicalVolu 154 G4VPhysicalVolume* motherPhysical )
155 { 155 {
156 // For the case of Voxel (or Parameterised) 156 // For the case of Voxel (or Parameterised) volume the respective
157 // sub-navigator must be messaged to update 157 // sub-navigator must be messaged to update its voxel information etc
158 158
159 // Update the state of the Sub Navigators 159 // Update the state of the Sub Navigators
160 // - in particular any voxel information the 160 // - in particular any voxel information they store/cache
161 // 161 //
162 G4LogicalVolume* motherLogical = motherP 162 G4LogicalVolume* motherLogical = motherPhysical->GetLogicalVolume();
163 G4SmartVoxelHeader* pVoxelHeader = motherL 163 G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
164 164
165 switch( CharacteriseDaughters(motherLogical) 165 switch( CharacteriseDaughters(motherLogical) )
166 { 166 {
167 case kNormal: 167 case kNormal:
168 if ( pVoxelHeader ) 168 if ( pVoxelHeader )
169 { 169 {
170 fvoxelNav.VoxelLocate( pVoxelHeader, p 170 fvoxelNav.VoxelLocate( pVoxelHeader, pointLocal );
171 } 171 }
172 break; 172 break;
173 case kParameterised: 173 case kParameterised:
174 if( GetDaughtersRegularStructureId(mothe 174 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
175 { 175 {
176 // Resets state & returns voxel node 176 // Resets state & returns voxel node
177 // 177 //
178 fparamNav.ParamVoxelLocate( pVoxelHead 178 fparamNav.ParamVoxelLocate( pVoxelHeader, pointLocal );
179 } 179 }
180 break; 180 break;
181 case kReplica: 181 case kReplica:
182 // Nothing to do 182 // Nothing to do
183 break; 183 break;
184 case kExternal: 184 case kExternal:
185 fpExternalNav->RelocateWithinVolume( mot 185 fpExternalNav->RelocateWithinVolume( motherPhysical,
186 poi 186 pointLocal );
187 break; 187 break;
188 } 188 }
189 } 189 }
190 190
191 // ******************************************* 191 // ********************************************************************
192 // Accessor for custom external navigation. 192 // Accessor for custom external navigation.
193 // ******************************************* 193 // ********************************************************************
194 G4VExternalNavigation* G4SafetyCalculator::Get 194 G4VExternalNavigation* G4SafetyCalculator::GetExternalNavigation() const
195 { 195 {
196 return fpExternalNav; 196 return fpExternalNav;
197 } 197 }
198 198
199 // ******************************************* 199 // ********************************************************************
200 // Modifier for custom external navigation. 200 // Modifier for custom external navigation.
201 // ******************************************* 201 // ********************************************************************
202 void G4SafetyCalculator::SetExternalNavigation 202 void G4SafetyCalculator::SetExternalNavigation(G4VExternalNavigation* eNav)
203 { 203 {
204 fpExternalNav = eNav; 204 fpExternalNav = eNav;
205 } 205 }
206 206
207 // ******************************************* 207 // ********************************************************************
208 // CompareSafetyValues 208 // CompareSafetyValues
209 // ******************************************* 209 // ********************************************************************
210 void G4SafetyCalculator:: 210 void G4SafetyCalculator::
211 CompareSafetyValues( G4double oldSafety, 211 CompareSafetyValues( G4double oldSafety,
212 G4double newValue, 212 G4double newValue,
213 G4VPhysicalVolume* mother 213 G4VPhysicalVolume* motherPhysical,
214 const G4ThreeVector& globalPoin 214 const G4ThreeVector& globalPoint,
215 G4bool keepState, 215 G4bool keepState,
216 G4double maxLength, 216 G4double maxLength,
217 G4bool enteredDaughterVol 217 G4bool enteredDaughterVol,
218 G4bool exitedMotherVol ) 218 G4bool exitedMotherVol )
219 { 219 {
220 constexpr G4double reportThreshold= 3.0e-14; 220 constexpr G4double reportThreshold= 3.0e-14;
221 // At least warn if rel-error exceeds it 221 // At least warn if rel-error exceeds it
222 constexpr G4double errorThreshold= 1.0e-08; 222 constexpr G4double errorThreshold= 1.0e-08;
223 // Fatal if relative error is larger 223 // Fatal if relative error is larger
224 constexpr G4double epsilonLen= 1.0e-20; 224 constexpr G4double epsilonLen= 1.0e-20;
225 // Baseline minimal value for divisor 225 // Baseline minimal value for divisor
226 226
227 const G4double oldSafetyPlus = std::fabs(old 227 const G4double oldSafetyPlus = std::fabs(oldSafety)+epsilonLen;
228 if( std::fabs( newValue - oldSafety) > repor 228 if( std::fabs( newValue - oldSafety) > reportThreshold * oldSafetyPlus )
229 { 229 {
230 G4ExceptionSeverity severity= FatalExcepti 230 G4ExceptionSeverity severity= FatalException;
231 std::ostringstream msg; 231 std::ostringstream msg;
232 G4double diff= (newValue-oldSafety); 232 G4double diff= (newValue-oldSafety);
233 G4double relativeDiff= diff / oldSafetyPlu 233 G4double relativeDiff= diff / oldSafetyPlus;
234 234
235 msg << " New (G4SafetyCalculator) value *d 235 msg << " New (G4SafetyCalculator) value *disagrees* by relative diff " << relativeDiff
236 << " in physical volume '" << motherPh 236 << " in physical volume '" << motherPhysical->GetName() << "' "
237 << "copy-no = " << motherPhysical->Get 237 << "copy-no = " << motherPhysical->GetCopyNo();
238 if( enteredDaughterVol ) { msg << " ( Jus 238 if( enteredDaughterVol ) { msg << " ( Just Entered new daughter volume. ) "; }
239 if( exitedMotherVol ) { msg << " ( Jus 239 if( exitedMotherVol ) { msg << " ( Just Exited previous volume. ) "; }
240 msg << G4endl; 240 msg << G4endl;
241 msg << " Safeties: old= " << std::setpr 241 msg << " Safeties: old= " << std::setprecision(12) << oldSafety
242 << " trial " << newValue 242 << " trial " << newValue
243 << " new-old= " << std::setprecision( 243 << " new-old= " << std::setprecision(7) << diff << G4endl;
244 244
245 if( std::fabs(diff) < errorThreshold * ( s 245 if( std::fabs(diff) < errorThreshold * ( std::fabs(oldSafety)+1.0e-20 ) )
246 { 246 {
247 msg << " (tiny difference) "; 247 msg << " (tiny difference) ";
248 severity= JustWarning; 248 severity= JustWarning;
249 } 249 }
250 else 250 else
251 { 251 {
252 msg << " (real difference) "; 252 msg << " (real difference) ";
253 severity= FatalException; 253 severity= FatalException;
254 254
255 // Extra information -- for big errors 255 // Extra information -- for big errors
256 msg << " NOTE: keepState = " << keepSt 256 msg << " NOTE: keepState = " << keepState << G4endl;
257 msg << " Location - Global coordinates: 257 msg << " Location - Global coordinates: " << globalPoint
258 << " volume= '" << motherPhysical-> 258 << " volume= '" << motherPhysical->GetName() << "'"
259 << " copy-no= " << motherPhysical->G 259 << " copy-no= " << motherPhysical->GetCopyNo() << G4endl;
260 msg << " Argument maxLength= " << maxLen 260 msg << " Argument maxLength= " << maxLength << G4endl;
261 261
262 std::size_t depth= fNavHistory.GetDepth( 262 std::size_t depth= fNavHistory.GetDepth();
263 msg << " Navigation History: depth = " < 263 msg << " Navigation History: depth = " << depth << G4endl;
264 for( G4int i=1; i<(G4int)depth; ++i ) 264 for( G4int i=1; i<(G4int)depth; ++i )
265 { 265 {
266 msg << " d= " << i << " " << std: 266 msg << " d= " << i << " " << std::setw(32)
267 << fNavHistory.GetVolume(i)->GetN 267 << fNavHistory.GetVolume(i)->GetName()
268 << " copyNo= " << fNavHistory.Ge 268 << " copyNo= " << fNavHistory.GetReplicaNo(i);
269 msg << G4endl; 269 msg << G4endl;
270 } 270 }
271 } 271 }
272 272
273 #ifdef G4DEBUG_NAVIGATION 273 #ifdef G4DEBUG_NAVIGATION
274 G4double redo= SafetyInCurrentVolume(globa 274 G4double redo= SafetyInCurrentVolume(globalPoint, motherPhysical,
275 maxLe 275 maxLength, true);
276 msg << " Redoing estimator: value = " << s 276 msg << " Redoing estimator: value = " << std::setprecision(16) << redo
277 << " diff/last= " << std::setprecisio 277 << " diff/last= " << std::setprecision(7) << redo - newValue
278 << " diff/old= " << redo - oldSafety 278 << " diff/old= " << redo - oldSafety << G4endl;
279 #endif 279 #endif
280 280
281 G4Exception("G4SafetyCalculator::CompareSa 281 G4Exception("G4SafetyCalculator::CompareSafetyValues()", "GeomNav1007",
282 severity, msg); 282 severity, msg);
283 } 283 }
284 } 284 }
285 285