| 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 //
>> 26 //
>> 27 // $Id: G4Navigator.cc,v 1.37 2007/10/18 14:18:36 gcosmo Exp $
>> 28 // GEANT4 tag $ Name: $
25 // 29 //
26 // G4Navigator class Implementation << 30 // class G4Navigator Implementation
27 // 31 //
28 // Original author: Paul Kent, July 95/96 32 // Original author: Paul Kent, July 95/96
29 // Responsible 1996-present: John Apostolakis, << 33 //
30 // Additional revisions by: Pedro Arce, Vladim <<
31 // ------------------------------------------- 34 // --------------------------------------------------------------------
32 35
33 #include <iomanip> <<
34 <<
35 #include "G4Navigator.hh" 36 #include "G4Navigator.hh"
36 #include "G4ios.hh" 37 #include "G4ios.hh"
37 #include "G4SystemOfUnits.hh" << 38 #include <iomanip>
>> 39
38 #include "G4GeometryTolerance.hh" 40 #include "G4GeometryTolerance.hh"
39 #include "G4VPhysicalVolume.hh" 41 #include "G4VPhysicalVolume.hh"
40 42
41 #include "G4VoxelSafety.hh" <<
42 #include "G4SafetyCalculator.hh" <<
43 <<
44 // Constant determining how precise normals sh <<
45 // vectors). If exceeded, warnings will be iss <<
46 // Can be CLHEP::perMillion (its old default) <<
47 // <<
48 static const G4double kToleranceNormalCheck = <<
49 <<
50 // ******************************************* 43 // ********************************************************************
51 // Constructor 44 // Constructor
52 // ******************************************* 45 // ********************************************************************
53 // 46 //
54 G4Navigator::G4Navigator() 47 G4Navigator::G4Navigator()
>> 48 : fWasLimitedByGeometry(false), fVerbose(0),
>> 49 fTopPhysical(0), fCheck(false), fPushed(false)
55 { 50 {
>> 51 fActive= false;
56 ResetStackAndState(); 52 ResetStackAndState();
57 // Initialises also all <<
58 // - exit / entry flags <<
59 // - flags & variables for exit normals <<
60 // - zero step counters <<
61 // - blocked volume <<
62 53
63 if( fVerbose > 2 ) << 54 fActionThreshold_NoZeroSteps = 10;
64 { << 55 fAbandonThreshold_NoZeroSteps = 25;
65 G4cout << " G4Navigator parameters: Action <<
66 << fActionThreshold_NoZeroSteps <<
67 << " Abandon Threshold (No Zero St <<
68 << fAbandonThreshold_NoZeroSteps << <<
69 } <<
70 kCarTolerance = G4GeometryTolerance::GetInst <<
71 fMinStep = 0.05*kCarTolerance; <<
72 fSqTol = sqr(kCarTolerance); <<
73 56
>> 57 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
74 fregularNav.SetNormalNavigation( &fnormalNav 58 fregularNav.SetNormalNavigation( &fnormalNav );
75 59
76 fStepEndPoint = G4ThreeVector( kInfinity, kI 60 fStepEndPoint = G4ThreeVector( kInfinity, kInfinity, kInfinity );
77 fLastStepEndPointLocal = G4ThreeVector( kInf <<
78 <<
79 fpVoxelSafety = new G4VoxelSafety(); <<
80 fpvoxelNav = new G4VoxelNavigation(); <<
81 fpSafetyCalculator = new G4SafetyCalculator( <<
82 fpSafetyCalculator->SetExternalNavigation(fp <<
83 } 61 }
84 62
85 // ******************************************* 63 // ********************************************************************
86 // Destructor 64 // Destructor
87 // ******************************************* 65 // ********************************************************************
88 // 66 //
89 G4Navigator::~G4Navigator() 67 G4Navigator::~G4Navigator()
90 { << 68 {;}
91 delete fpVoxelSafety; <<
92 delete fpExternalNav; <<
93 delete fpvoxelNav; <<
94 delete fpSafetyCalculator; <<
95 } <<
96 69
97 // ******************************************* 70 // ********************************************************************
98 // ResetHierarchyAndLocate 71 // ResetHierarchyAndLocate
99 // ******************************************* 72 // ********************************************************************
100 // 73 //
101 G4VPhysicalVolume* 74 G4VPhysicalVolume*
102 G4Navigator::ResetHierarchyAndLocate(const G4T << 75 G4Navigator::ResetHierarchyAndLocate(const G4ThreeVector &p,
103 const G4T << 76 const G4ThreeVector &direction,
104 const G4T << 77 const G4TouchableHistory &h)
105 { 78 {
106 ResetState(); 79 ResetState();
107 fHistory = *h.GetHistory(); 80 fHistory = *h.GetHistory();
108 SetupHierarchy(); 81 SetupHierarchy();
109 fLastTriedStepComputation = false; // Redun <<
110 return LocateGlobalPointAndSetup(p, &directi 82 return LocateGlobalPointAndSetup(p, &direction, true, false);
111 } 83 }
112 84
113 // ******************************************* 85 // ********************************************************************
114 // LocateGlobalPointAndSetup 86 // LocateGlobalPointAndSetup
115 // 87 //
116 // Locate the point in the hierarchy return 0 88 // Locate the point in the hierarchy return 0 if outside
117 // The direction is required 89 // The direction is required
118 // - if on an edge shared by more than two 90 // - if on an edge shared by more than two surfaces
119 // (to resolve likely looping in tracking 91 // (to resolve likely looping in tracking)
120 // - at initial location of a particle 92 // - at initial location of a particle
121 // (to resolve potential ambiguity at bou 93 // (to resolve potential ambiguity at boundary)
122 // 94 //
123 // Flags on exit: (comments to be completed) 95 // Flags on exit: (comments to be completed)
124 // fEntering - True if entering `daugh 96 // fEntering - True if entering `daughter' volume (or replica)
125 // whether daughter of las 97 // whether daughter of last mother directly
126 // or daughter of that vol 98 // or daughter of that volume's ancestor.
127 // fExiting - True if exited 'mother' <<
128 // (always ? - how about i <<
129 // ******************************************* 99 // ********************************************************************
130 // 100 //
131 G4VPhysicalVolume* 101 G4VPhysicalVolume*
132 G4Navigator::LocateGlobalPointAndSetup( const 102 G4Navigator::LocateGlobalPointAndSetup( const G4ThreeVector& globalPoint,
133 const 103 const G4ThreeVector* pGlobalDirection,
134 const 104 const G4bool relativeSearch,
135 const 105 const G4bool ignoreDirection )
136 { 106 {
137 G4bool notKnownContained = true, noResult; << 107 G4bool notKnownContained=true, noResult;
138 G4VPhysicalVolume *targetPhysical; 108 G4VPhysicalVolume *targetPhysical;
139 G4LogicalVolume *targetLogical; 109 G4LogicalVolume *targetLogical;
140 G4VSolid *targetSolid = nullptr; << 110 G4VSolid *targetSolid=0;
141 G4ThreeVector localPoint, globalDirection; 111 G4ThreeVector localPoint, globalDirection;
142 EInside insideCode; 112 EInside insideCode;
143 << 113
144 G4bool considerDirection = (pGlobalDirection << 114 G4bool considerDirection = (!ignoreDirection) || fLocatedOnEdge;
145 << 115
146 fLastTriedStepComputation = false; << 116 if( considerDirection && pGlobalDirection != 0 )
147 fChangedGrandMotherRefFrame = false; // For <<
148 <<
149 if( considerDirection ) <<
150 { 117 {
151 globalDirection=*pGlobalDirection; 118 globalDirection=*pGlobalDirection;
152 } 119 }
153 120
>> 121 #ifdef G4DEBUG_NAVIGATION
>> 122 if( fVerbose > 2 )
>> 123 {
>> 124 G4cout << "Upon entering LocateGlobalPointAndSetup():" << G4endl;
>> 125 G4cout << " History = " << G4endl << fHistory << G4endl << G4endl;
>> 126 }
>> 127 #endif
>> 128
154 #ifdef G4VERBOSE 129 #ifdef G4VERBOSE
>> 130 G4int oldcoutPrec = G4cout.precision(8);
155 if( fVerbose > 2 ) 131 if( fVerbose > 2 )
156 { 132 {
157 G4long oldcoutPrec = G4cout.precision(8); <<
158 G4cout << "*** G4Navigator::LocateGlobalPo 133 G4cout << "*** G4Navigator::LocateGlobalPointAndSetup: ***" << G4endl;
159 G4cout << " Called with arguments: " << 134 G4cout << " Called with arguments: " << G4endl
160 << " Globalpoint = " << glob << 135 << " Globalpoint = " << globalPoint << G4endl
161 << " RelativeSearch = " << r << 136 << " RelativeSearch = " << relativeSearch << G4endl;
162 if( fVerbose >= 4 ) << 137 if( fVerbose == 4 )
163 { 138 {
164 G4cout << " ----- Upon entering:" << 139 G4cout << " ----- Upon entering:" << G4endl;
165 PrintState(); 140 PrintState();
166 } 141 }
167 G4cout.precision(oldcoutPrec); <<
168 } 142 }
169 #endif 143 #endif
170 144
171 G4int noLevelsExited = 0; <<
172 <<
173 if ( !relativeSearch ) 145 if ( !relativeSearch )
174 { 146 {
175 ResetStackAndState(); 147 ResetStackAndState();
176 } 148 }
177 else 149 else
178 { 150 {
179 if ( fWasLimitedByGeometry ) 151 if ( fWasLimitedByGeometry )
180 { 152 {
181 fWasLimitedByGeometry = false; 153 fWasLimitedByGeometry = false;
182 fEnteredDaughter = fEntering; // Remem 154 fEnteredDaughter = fEntering; // Remember
183 fExitedMother = fExiting; // Remem 155 fExitedMother = fExiting; // Remember
184 if ( fExiting ) 156 if ( fExiting )
185 { 157 {
186 ++noLevelsExited; // count this first << 158 if ( fHistory.GetDepth() )
187 <<
188 if ( fHistory.GetDepth() != 0 ) <<
189 { 159 {
190 fBlockedPhysicalVolume = fHistory.Ge 160 fBlockedPhysicalVolume = fHistory.GetTopVolume();
191 fBlockedReplicaNo = fHistory.GetTopR 161 fBlockedReplicaNo = fHistory.GetTopReplicaNo();
192 fHistory.BackLevel(); 162 fHistory.BackLevel();
193 } 163 }
194 else 164 else
195 { 165 {
196 fLastLocatedPointLocal = localPoint; 166 fLastLocatedPointLocal = localPoint;
197 fLocatedOutsideWorld = true; 167 fLocatedOutsideWorld = true;
198 fBlockedPhysicalVolume = nullptr; << 168 return 0; // Have exited world volume
199 fBlockedReplicaNo = -1; <<
200 fEntering = false; // No <<
201 fEnteredDaughter = false; <<
202 fExitedMother = true; // ?? <<
203 <<
204 return nullptr; // Have ex <<
205 } 169 }
206 // A fix for the case where a volume i 170 // A fix for the case where a volume is "entered" at an edge
207 // and a coincident surface exists out 171 // and a coincident surface exists outside it.
208 // - This stops it from exiting furth 172 // - This stops it from exiting further volumes and cycling
209 // - However ReplicaNavigator treats 173 // - However ReplicaNavigator treats this case itself
210 // 174 //
211 // assert( fBlockedPhysicalVolume!=0 ) <<
212 <<
213 // Expect to be on edge => on surface <<
214 // <<
215 if ( fLocatedOnEdge && (VolumeType(fBl 175 if ( fLocatedOnEdge && (VolumeType(fBlockedPhysicalVolume)!=kReplica ))
216 { 176 {
217 fExiting = false; << 177 fExiting= false;
218 // Consider effect on Exit Normal !? <<
219 } 178 }
220 } 179 }
221 else 180 else
222 if ( fEntering ) 181 if ( fEntering )
223 { 182 {
224 switch (VolumeType(fBlockedPhysicalV 183 switch (VolumeType(fBlockedPhysicalVolume))
225 { 184 {
226 case kNormal: 185 case kNormal:
227 fHistory.NewLevel(fBlockedPhysic 186 fHistory.NewLevel(fBlockedPhysicalVolume, kNormal,
228 fBlockedPhysic 187 fBlockedPhysicalVolume->GetCopyNo());
229 break; 188 break;
230 case kReplica: 189 case kReplica:
231 freplicaNav.ComputeTransformatio 190 freplicaNav.ComputeTransformation(fBlockedReplicaNo,
232 191 fBlockedPhysicalVolume);
233 fHistory.NewLevel(fBlockedPhysic 192 fHistory.NewLevel(fBlockedPhysicalVolume, kReplica,
234 fBlockedReplic 193 fBlockedReplicaNo);
235 fBlockedPhysicalVolume->SetCopyN 194 fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
236 break; 195 break;
237 case kParameterised: 196 case kParameterised:
238 if( fBlockedPhysicalVolume->GetR << 197 if( fBlockedPhysicalVolume->GetRegularStructureId() != 1 )
239 { 198 {
240 G4VSolid *pSolid; 199 G4VSolid *pSolid;
241 G4VPVParameterisation *pParam; 200 G4VPVParameterisation *pParam;
242 G4TouchableHistory parentTouch 201 G4TouchableHistory parentTouchable( fHistory );
243 pParam = fBlockedPhysicalVolum 202 pParam = fBlockedPhysicalVolume->GetParameterisation();
244 pSolid = pParam->ComputeSolid( 203 pSolid = pParam->ComputeSolid(fBlockedReplicaNo,
245 204 fBlockedPhysicalVolume);
246 pSolid->ComputeDimensions(pPar 205 pSolid->ComputeDimensions(pParam, fBlockedReplicaNo,
247 fBlo 206 fBlockedPhysicalVolume);
248 pParam->ComputeTransformation( 207 pParam->ComputeTransformation(fBlockedReplicaNo,
249 208 fBlockedPhysicalVolume);
250 fHistory.NewLevel(fBlockedPhys 209 fHistory.NewLevel(fBlockedPhysicalVolume, kParameterised,
251 fBlockedRepl 210 fBlockedReplicaNo);
252 fBlockedPhysicalVolume->SetCop 211 fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
253 // 212 //
254 // Set the correct solid and m 213 // Set the correct solid and material in Logical Volume
255 // 214 //
256 G4LogicalVolume *pLogical; 215 G4LogicalVolume *pLogical;
257 pLogical = fBlockedPhysicalVol 216 pLogical = fBlockedPhysicalVolume->GetLogicalVolume();
258 pLogical->SetSolid( pSolid ); 217 pLogical->SetSolid( pSolid );
259 pLogical->UpdateMaterial(pPara 218 pLogical->UpdateMaterial(pParam ->
260 ComputeMaterial(fBlockedRepl 219 ComputeMaterial(fBlockedReplicaNo,
261 fBlockedPhys 220 fBlockedPhysicalVolume,
262 &parentTouch 221 &parentTouchable));
263 } 222 }
264 break; 223 break;
265 case kExternal: <<
266 G4Exception("G4Navigator::Locate <<
267 "GeomNav0001", Fatal <<
268 "Extra levels not ap <<
269 break; <<
270 } 224 }
271 fEntering = false; 225 fEntering = false;
272 fBlockedPhysicalVolume = nullptr; << 226 fBlockedPhysicalVolume = 0;
273 localPoint = fHistory.GetTopTransfor 227 localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
274 notKnownContained = false; 228 notKnownContained = false;
275 } 229 }
276 } 230 }
277 else 231 else
278 { 232 {
279 fBlockedPhysicalVolume = nullptr; << 233 fBlockedPhysicalVolume = 0;
280 fEntering = false; 234 fEntering = false;
281 fEnteredDaughter = false; // Full Step 235 fEnteredDaughter = false; // Full Step was not taken, did not enter
282 fExiting = false; 236 fExiting = false;
283 fExitedMother = false; // Full Step 237 fExitedMother = false; // Full Step was not taken, did not exit
284 } 238 }
285 } 239 }
286 // 240 //
287 // Search from top of history up through geo 241 // Search from top of history up through geometry until
288 // containing volume found: 242 // containing volume found:
289 // If on 243 // If on
290 // o OUTSIDE - Back up level, not/no longer 244 // o OUTSIDE - Back up level, not/no longer exiting volumes
291 // o SURFACE and EXITING - Back up level, se 245 // o SURFACE and EXITING - Back up level, setting new blocking no.s
292 // else 246 // else
293 // o containing volume found 247 // o containing volume found
294 // 248 //
295 << 249 while (notKnownContained)
296 while (notKnownContained) // Loop checking, <<
297 { 250 {
298 EVolume topVolumeType = fHistory.GetTopVol << 251 if ( fHistory.GetTopVolumeType()!=kReplica )
299 if (topVolumeType!=kReplica && topVolumeTy <<
300 { 252 {
301 targetSolid = fHistory.GetTopVolume()->G 253 targetSolid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid();
302 localPoint = fHistory.GetTopTransform(). 254 localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
303 insideCode = targetSolid->Inside(localPo 255 insideCode = targetSolid->Inside(localPoint);
304 #ifdef G4VERBOSE 256 #ifdef G4VERBOSE
305 if(( fVerbose == 1 ) && ( fCheck )) 257 if(( fVerbose == 1 ) && ( fCheck ))
306 { 258 {
307 G4String solidResponse = "-kInside-"; << 259 G4String solidResponse = "-kInside-";
308 if (insideCode == kOutside) << 260 if (insideCode == kOutside)
309 { << 261 solidResponse = "-kOutside-";
310 solidResponse = "-kOutside-"; << 262 else if (insideCode == kSurface)
311 } << 263 solidResponse = "-kSurface-";
312 else if (insideCode == kSurface) << 264 G4cout << "*** G4Navigator::LocateGlobalPointAndSetup(): ***" << G4endl
313 { << 265 << " Invoked Inside() for solid: " << targetSolid->GetName()
314 solidResponse = "-kSurface-"; << 266 << ". Solid replied: " << solidResponse << G4endl
315 } << 267 << " For local point p: " << localPoint << G4endl;
316 G4cout << "*** G4Navigator::LocateGlob <<
317 << " Invoked Inside() for so <<
318 << ". Solid replied: " << solid <<
319 << " For local point p: " << <<
320 } 268 }
321 #endif 269 #endif
322 } 270 }
323 else 271 else
324 { 272 {
325 if( topVolumeType == kReplica ) << 273 insideCode = freplicaNav.BackLocate(fHistory, globalPoint, localPoint,
326 { << 274 fExiting, notKnownContained);
327 insideCode = freplicaNav.BackLocate( << 275 // !CARE! if notKnownContained returns false then the point is within
328 << 276 // the containing placement volume of the replica(s). If insidecode
329 // !CARE! if notKnownContained retur << 277 // will result in the history being backed up one level, then the
330 // the containing placement volume o << 278 // local point returned is the point in the system of this new level
331 // will result in the history being <<
332 // local point returned is the point <<
333 } <<
334 else <<
335 { <<
336 targetSolid = fHistory.GetTopVolume( <<
337 localPoint = fHistory.GetTopTransfor <<
338 G4ThreeVector localDirection = <<
339 fHistory.GetTopTransform().Transf <<
340 insideCode = fpExternalNav->Inside(t <<
341 } <<
342 } 279 }
343 << 280 if ( insideCode==kOutside )
344 // Point is inside current volume, break o <<
345 if ( insideCode == kInside ) { break; } <<
346 <<
347 // Point is outside current volume, move u <<
348 if ( insideCode == kOutside ) <<
349 { 281 {
350 ++noLevelsExited; << 282 if ( fHistory.GetDepth() )
351 <<
352 // Exiting world volume <<
353 if ( fHistory.GetDepth() == 0 ) <<
354 { 283 {
355 fLocatedOutsideWorld = true; << 284 fBlockedPhysicalVolume = fHistory.GetTopVolume();
356 fLastLocatedPointLocal = localPoint; << 285 fBlockedReplicaNo = fHistory.GetTopReplicaNo();
357 return nullptr; << 286 fHistory.BackLevel();
>> 287 fExiting = false;
358 } 288 }
359 << 289 else
360 fBlockedPhysicalVolume = fHistory.GetTop <<
361 fBlockedReplicaNo = fHistory.GetTopRepli <<
362 fHistory.BackLevel(); <<
363 fExiting = false; <<
364 <<
365 if( noLevelsExited > 1 ) <<
366 { 290 {
367 // The first transformation was done b << 291 fLastLocatedPointLocal = localPoint;
368 // << 292 fLocatedOutsideWorld = true;
369 if(const auto *mRot = fBlockedPhysical << 293 return 0; // Have exited world volume
370 { <<
371 fGrandMotherExitNormal *= (*mRot).in <<
372 fChangedGrandMotherRefFrame = true; <<
373 } <<
374 } 294 }
375 continue; <<
376 } 295 }
377 << 296 else
378 // Point is on the surface of a volume << 297 if ( insideCode==kSurface )
379 G4bool isExiting = fExiting; <<
380 if( (!fExiting) && considerDirection ) <<
381 { <<
382 // Figure out whether we are exiting thi <<
383 // by using the direction <<
384 // <<
385 G4bool directionExiting = false; <<
386 G4ThreeVector localDirection = <<
387 fHistory.GetTopTransform().TransformAx <<
388 <<
389 // Make sure localPoint in correct refer <<
390 // ( Was it already correct ? How ? <<
391 // <<
392 localPoint= fHistory.GetTopTransform().T <<
393 if ( fHistory.GetTopVolumeType() != kRep <<
394 { 298 {
395 G4ThreeVector normal = targetSolid->Su << 299 G4bool isExiting = fExiting;
396 directionExiting = normal.dot(localDir << 300 if( (!fExiting)&&considerDirection )
397 isExiting = isExiting || directionExit << 301 {
>> 302 // Figure out whether we are exiting this level's volume
>> 303 // by using the direction
>> 304 //
>> 305 G4bool directionExiting = false;
>> 306 G4ThreeVector localDirection =
>> 307 fHistory.GetTopTransform().TransformAxis(globalDirection);
>> 308 if ( fHistory.GetTopVolumeType()!=kReplica )
>> 309 {
>> 310 G4ThreeVector normal = targetSolid->SurfaceNormal(localPoint);
>> 311 directionExiting = normal.dot(localDirection) > 0.0;
>> 312 isExiting = isExiting || directionExiting;
>> 313 }
>> 314 }
>> 315 if( isExiting )
>> 316 {
>> 317 if ( fHistory.GetDepth() )
>> 318 {
>> 319 fBlockedPhysicalVolume = fHistory.GetTopVolume();
>> 320 fBlockedReplicaNo = fHistory.GetTopReplicaNo();
>> 321 fHistory.BackLevel();
>> 322 //
>> 323 // Still on surface but exited volume not necessarily convex
>> 324 //
>> 325 fValidExitNormal = false;
>> 326 }
>> 327 else
>> 328 {
>> 329 fLastLocatedPointLocal = localPoint;
>> 330 fLocatedOutsideWorld = true;
>> 331 return 0; // Have exited world volume
>> 332 }
>> 333 }
>> 334 else
>> 335 {
>> 336 notKnownContained=false;
>> 337 }
398 } 338 }
399 } << 339 else
400 <<
401 // Point is on a surface, but no longer ex <<
402 if ( !isExiting ) { break; } <<
403 <<
404 ++noLevelsExited; <<
405 <<
406 // Point is on the outer surface, leaving <<
407 if ( fHistory.GetDepth() == 0 ) <<
408 { <<
409 fLocatedOutsideWorld = true; <<
410 fLastLocatedPointLocal = localPoint; <<
411 return nullptr; <<
412 } <<
413 <<
414 // Point is still on a surface, but exited <<
415 fValidExitNormal = false; <<
416 fBlockedPhysicalVolume = fHistory.GetTopVo <<
417 fBlockedReplicaNo = fHistory.GetTopReplica <<
418 fHistory.BackLevel(); <<
419 <<
420 if( noLevelsExited > 1 ) <<
421 { <<
422 // The first transformation was done by <<
423 // <<
424 const G4RotationMatrix* mRot = <<
425 fBlockedPhysicalVolume->GetRotation(); <<
426 if( mRot != nullptr ) <<
427 { 340 {
428 fGrandMotherExitNormal *= (*mRot).inve << 341 notKnownContained=false;
429 fChangedGrandMotherRefFrame = true; <<
430 } 342 }
431 } <<
432 } // END while (notKnownContained) 343 } // END while (notKnownContained)
433 // 344 //
434 // Search downwards until deepest containing 345 // Search downwards until deepest containing volume found,
435 // blocking fBlockedPhysicalVolume/BlockedRe 346 // blocking fBlockedPhysicalVolume/BlockedReplicaNum
436 // 347 //
437 // 3 Cases: 348 // 3 Cases:
438 // 349 //
439 // o Parameterised daughters 350 // o Parameterised daughters
440 // =>Must be one G4PVParameterised daughte 351 // =>Must be one G4PVParameterised daughter & voxels
441 // o Positioned daughters & voxels 352 // o Positioned daughters & voxels
442 // o Positioned daughters & no voxels 353 // o Positioned daughters & no voxels
443 354
444 noResult = true; // noResult should be rena << 355 noResult = true; // noResult should be renamed to
445 // something like enteredL 356 // something like enteredLevel, as that is its meaning.
446 do 357 do
447 { 358 {
448 // Determine `type' of current mother volu 359 // Determine `type' of current mother volume
449 // 360 //
450 targetPhysical = fHistory.GetTopVolume(); 361 targetPhysical = fHistory.GetTopVolume();
451 if (targetPhysical == nullptr) { break; } <<
452 targetLogical = targetPhysical->GetLogical 362 targetLogical = targetPhysical->GetLogicalVolume();
453 switch( CharacteriseDaughters(targetLogica 363 switch( CharacteriseDaughters(targetLogical) )
454 { 364 {
455 case kNormal: 365 case kNormal:
456 if ( targetLogical->GetVoxelHeader() ! << 366 if ( targetLogical->GetVoxelHeader() ) // use optimised navigation
457 { 367 {
458 noResult = GetVoxelNavigator().Level << 368 noResult = fvoxelNav.LevelLocate(fHistory,
459 fBl 369 fBlockedPhysicalVolume,
460 fBl 370 fBlockedReplicaNo,
461 glo 371 globalPoint,
462 pGl 372 pGlobalDirection,
463 con 373 considerDirection,
464 loc 374 localPoint);
465 } 375 }
466 else // do not u 376 else // do not use optimised navigation
467 { 377 {
468 noResult = fnormalNav.LevelLocate(fH 378 noResult = fnormalNav.LevelLocate(fHistory,
469 fB 379 fBlockedPhysicalVolume,
470 fB 380 fBlockedReplicaNo,
471 gl 381 globalPoint,
472 pG 382 pGlobalDirection,
473 co 383 considerDirection,
474 lo 384 localPoint);
475 } 385 }
476 break; 386 break;
477 case kReplica: 387 case kReplica:
478 noResult = freplicaNav.LevelLocate(fHi 388 noResult = freplicaNav.LevelLocate(fHistory,
479 fBl 389 fBlockedPhysicalVolume,
480 fBl 390 fBlockedReplicaNo,
481 glo 391 globalPoint,
482 pGl 392 pGlobalDirection,
483 con 393 considerDirection,
484 loc 394 localPoint);
485 break; 395 break;
486 case kParameterised: 396 case kParameterised:
487 if( GetDaughtersRegularStructureId(tar 397 if( GetDaughtersRegularStructureId(targetLogical) != 1 )
488 { 398 {
489 noResult = fparamNav.LevelLocate(fHi 399 noResult = fparamNav.LevelLocate(fHistory,
490 fBl 400 fBlockedPhysicalVolume,
491 fBl 401 fBlockedReplicaNo,
492 glo 402 globalPoint,
493 pGl 403 pGlobalDirection,
494 con 404 considerDirection,
495 loc 405 localPoint);
496 } 406 }
497 else // Regular structure 407 else // Regular structure
498 { 408 {
499 noResult = fregularNav.LevelLocate(f 409 noResult = fregularNav.LevelLocate(fHistory,
500 f 410 fBlockedPhysicalVolume,
501 f 411 fBlockedReplicaNo,
502 g 412 globalPoint,
503 p 413 pGlobalDirection,
504 c 414 considerDirection,
505 l 415 localPoint);
506 } 416 }
507 break; 417 break;
508 case kExternal: <<
509 noResult = fpExternalNav->LevelLocate( <<
510 <<
511 <<
512 <<
513 <<
514 <<
515 <<
516 break; <<
517 } 418 }
518 419
519 // LevelLocate returns true if it finds a 420 // LevelLocate returns true if it finds a daughter volume
520 // in which globalPoint is inside (or on t 421 // in which globalPoint is inside (or on the surface).
521 422
522 if ( noResult ) 423 if ( noResult )
523 { 424 {
524 // Entering a daughter after ascending 425 // Entering a daughter after ascending
525 // 426 //
526 // The blocked volume is no longer valid 427 // The blocked volume is no longer valid - it was for another level
527 // 428 //
528 fBlockedPhysicalVolume = nullptr; << 429 fBlockedPhysicalVolume = 0;
529 fBlockedReplicaNo = -1; 430 fBlockedReplicaNo = -1;
530 431
531 // fEntering should be false -- else blo 432 // fEntering should be false -- else blockedVolume is assumed good.
532 // fEnteredDaughter is used for ExitNorm 433 // fEnteredDaughter is used for ExitNormal
533 // 434 //
534 fEntering = false; 435 fEntering = false;
535 fEnteredDaughter = true; 436 fEnteredDaughter = true;
536 <<
537 if( fExitedMother ) <<
538 { <<
539 G4VPhysicalVolume* enteredPhysical = f <<
540 const G4RotationMatrix* mRot = entered <<
541 if( mRot != nullptr ) <<
542 { <<
543 // Go deeper, i.e. move 'down' in th <<
544 // Apply direct rotation, not invers <<
545 // <<
546 fGrandMotherExitNormal *= (*mRot); <<
547 fChangedGrandMotherRefFrame= true; <<
548 } <<
549 } <<
550 <<
551 #ifdef G4DEBUG_NAVIGATION 437 #ifdef G4DEBUG_NAVIGATION
552 if( fVerbose > 2 ) 438 if( fVerbose > 2 )
553 { 439 {
554 G4VPhysicalVolume* enteredPhysical = 440 G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume();
555 G4cout << "*** G4Navigator::LocateGlo << 441 G4cout << "*** G4Navigator::LocateGlobalPointAndSetup() ***" << G4endl;
556 G4cout << " Entering volume: " << 442 G4cout << " Entering volume: " << enteredPhysical->GetName()
557 << G4endl; 443 << G4endl;
558 } 444 }
559 #endif 445 #endif
560 } 446 }
561 } while (noResult); // Loop checking, 07.10 << 447 } while (noResult);
562 448
563 fLastLocatedPointLocal = localPoint; 449 fLastLocatedPointLocal = localPoint;
564 450
565 #ifdef G4VERBOSE 451 #ifdef G4VERBOSE
566 if( fVerbose >= 4 ) << 452 if( fVerbose == 4 )
567 { 453 {
568 G4long oldcoutPrec = G4cout.precision(8); << 454 G4cout.precision(6);
569 G4String curPhysVol_Name("None"); 455 G4String curPhysVol_Name("None");
570 if (targetPhysical != nullptr) { curPhysV << 456 if (targetPhysical!=0)
>> 457 {
>> 458 curPhysVol_Name = targetPhysical->GetName();
>> 459 }
571 G4cout << " Return value = new volume = 460 G4cout << " Return value = new volume = " << curPhysVol_Name << G4endl;
572 G4cout << " ----- Upon exiting:" << G4e 461 G4cout << " ----- Upon exiting:" << G4endl;
573 PrintState(); 462 PrintState();
574 if( fVerbose >= 5 ) << 463 #ifdef G4DEBUG_NAVIGATION
575 { << 464 G4cout << "Upon exiting LocateGlobalPointAndSetup():" << G4endl;
576 G4cout << "Upon exiting LocateGlobalPoin << 465 G4cout << " History = " << G4endl << fHistory << G4endl << G4endl;
577 G4cout << " History = " << G4endl << << 466 #endif
578 } <<
579 G4cout.precision(oldcoutPrec); <<
580 } 467 }
>> 468 G4cout.precision(oldcoutPrec);
581 #endif 469 #endif
582 470
583 fLocatedOutsideWorld = false; << 471 fLocatedOutsideWorld= false;
584 472
585 return targetPhysical; 473 return targetPhysical;
586 } 474 }
587 475
588 // ******************************************* 476 // ********************************************************************
589 // LocateGlobalPointWithinVolume 477 // LocateGlobalPointWithinVolume
590 // 478 //
591 // -> the state information of this Navigator 479 // -> the state information of this Navigator and its subNavigators
592 // is updated in order to start the next st 480 // is updated in order to start the next step at pGlobalpoint
593 // -> no check is performed whether pGlobalpoi 481 // -> no check is performed whether pGlobalpoint is inside the
594 // original volume (this must be the case). 482 // original volume (this must be the case).
595 // 483 //
596 // Note: a direction could be added to the arg 484 // Note: a direction could be added to the arguments, to aid in future
597 // optional checking (via the old code b 485 // optional checking (via the old code below, flagged by OLD_LOCATE).
598 // [ This would be done only in verbose 486 // [ This would be done only in verbose mode ]
599 // ******************************************* 487 // ********************************************************************
600 // 488 //
601 void 489 void
602 G4Navigator::LocateGlobalPointWithinVolume(con 490 G4Navigator::LocateGlobalPointWithinVolume(const G4ThreeVector& pGlobalpoint)
603 { << 491 {
>> 492 fLastLocatedPointLocal = ComputeLocalPoint(pGlobalpoint);
>> 493
604 #ifdef G4DEBUG_NAVIGATION 494 #ifdef G4DEBUG_NAVIGATION
605 assert( !fWasLimitedByGeometry ); << 495 if( fVerbose > 2 )
606 // Check: Either step was not limited by a << 496 {
607 // else the full step is no longer << 497 G4cout << "Entering LocateGlobalWithinVolume(): History = " << G4endl;
>> 498 G4cout << fHistory << G4endl;
>> 499 }
608 #endif 500 #endif
609 <<
610 fLastLocatedPointLocal = ComputeLocalPoint( <<
611 fLastTriedStepComputation = false; <<
612 fChangedGrandMotherRefFrame = false; // F <<
613 501
614 // For the case of Voxel (or Parameterised) 502 // For the case of Voxel (or Parameterised) volume the respective
615 // Navigator must be messaged to update its 503 // Navigator must be messaged to update its voxel information etc
616 504
617 // Update the state of the Sub Navigators 505 // Update the state of the Sub Navigators
618 // - in particular any voxel information th 506 // - in particular any voxel information they store/cache
619 // 507 //
620 G4VPhysicalVolume* motherPhysical = fHisto 508 G4VPhysicalVolume* motherPhysical = fHistory.GetTopVolume();
621 G4LogicalVolume* motherLogical = mother 509 G4LogicalVolume* motherLogical = motherPhysical->GetLogicalVolume();
>> 510 G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
622 511
623 switch( CharacteriseDaughters(motherLogical << 512 if ( fHistory.GetTopVolumeType()!=kReplica )
624 { 513 {
>> 514 switch( CharacteriseDaughters(motherLogical) )
>> 515 {
625 case kNormal: 516 case kNormal:
626 GetVoxelNavigator().RelocateWithinVol << 517 if ( pVoxelHeader )
>> 518 {
>> 519 fvoxelNav.VoxelLocate( pVoxelHeader, fLastLocatedPointLocal );
>> 520 }
627 break; 521 break;
628 case kParameterised: 522 case kParameterised:
629 fparamNav.RelocateWithinVolume( mothe << 523 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
>> 524 {
>> 525 // Resets state & returns voxel node
>> 526 //
>> 527 fparamNav.ParamVoxelLocate( pVoxelHeader, fLastLocatedPointLocal );
>> 528 }
630 break; 529 break;
631 case kReplica: 530 case kReplica:
632 // Nothing to do << 531 G4Exception("G4Navigator::LocateGlobalPointWithinVolume()",
633 break; << 532 "NotApplicable", FatalException,
634 case kExternal: << 533 "Not applicable for replicated volumes.");
635 fpExternalNav->RelocateWithinVolume( <<
636 <<
637 break; 534 break;
>> 535 }
638 } 536 }
639 537
640 // Reset the state variables 538 // Reset the state variables
641 // - which would have been affected 539 // - which would have been affected
642 // by the 'equivalent' call to LocateGl 540 // by the 'equivalent' call to LocateGlobalPointAndSetup
643 // - who's values have been invalidated b 541 // - who's values have been invalidated by the 'move'.
644 // 542 //
645 fBlockedPhysicalVolume = nullptr; << 543 fBlockedPhysicalVolume = 0;
646 fBlockedReplicaNo = -1; 544 fBlockedReplicaNo = -1;
647 fEntering = false; 545 fEntering = false;
648 fEnteredDaughter = false; // Boundary not 546 fEnteredDaughter = false; // Boundary not encountered, did not enter
649 fExiting = false; 547 fExiting = false;
650 fExitedMother = false; // Boundary not 548 fExitedMother = false; // Boundary not encountered, did not exit
651 } 549 }
652 550
653 // ******************************************* 551 // ********************************************************************
654 // SetSavedState 552 // SetSavedState
655 // 553 //
656 // Save the state, in case this is a parasitic 554 // Save the state, in case this is a parasitic call
657 // Save fValidExitNormal, fExitNormal, fExitin 555 // Save fValidExitNormal, fExitNormal, fExiting, fEntering,
658 // fBlockedPhysicalVolume, fBlockedReplic 556 // fBlockedPhysicalVolume, fBlockedReplicaNo, fLastStepWasZero;
659 // ******************************************* 557 // ********************************************************************
660 // 558 //
661 void G4Navigator::SetSavedState() 559 void G4Navigator::SetSavedState()
662 { 560 {
663 // Note: the state of dependent objects is n << 561 // fSaveExitNormal = fExitNormal;
664 // ( This means that the full state is cha <<
665 // SetSavedState() and RestoreSavedState <<
666 <<
667 fSaveState.sExitNormal = fExitNormal; 562 fSaveState.sExitNormal = fExitNormal;
668 fSaveState.sValidExitNormal = fValidExitNorm 563 fSaveState.sValidExitNormal = fValidExitNormal;
669 fSaveState.sExiting = fExiting; 564 fSaveState.sExiting = fExiting;
670 fSaveState.sEntering = fEntering; 565 fSaveState.sEntering = fEntering;
671 566
672 fSaveState.spBlockedPhysicalVolume = fBlocke 567 fSaveState.spBlockedPhysicalVolume = fBlockedPhysicalVolume;
673 fSaveState.sBlockedReplicaNo = fBlockedRepli << 568 fSaveState.sBlockedReplicaNo = fBlockedReplicaNo,
674 <<
675 fSaveState.sLastStepWasZero = static_cast<G4 <<
676 <<
677 fSaveState.sLocatedOutsideWorld = fLocatedOu <<
678 fSaveState.sLastLocatedPointLocal = fLastLoc <<
679 fSaveState.sEnteredDaughter = fEnteredDaught <<
680 fSaveState.sExitedMother = fExitedMother; <<
681 fSaveState.sWasLimitedByGeometry = fWasLimit <<
682 569
683 // Even the safety sphere - if you want to c << 570 fSaveState.sLastStepWasZero = fLastStepWasZero;
684 // <<
685 fSaveState.sPreviousSftOrigin = fPreviousSft <<
686 fSaveState.sPreviousSafety = fPreviousSafety <<
687 } 571 }
688 572
689 // ******************************************* 573 // ********************************************************************
690 // RestoreSavedState 574 // RestoreSavedState
691 // 575 //
692 // Restore the state (in Compute Step), in cas 576 // Restore the state (in Compute Step), in case this is a parasitic call
693 // ******************************************* 577 // ********************************************************************
694 // 578 //
695 void G4Navigator::RestoreSavedState() 579 void G4Navigator::RestoreSavedState()
696 { 580 {
697 fExitNormal = fSaveState.sExitNormal; 581 fExitNormal = fSaveState.sExitNormal;
698 fValidExitNormal = fSaveState.sValidExitNorm 582 fValidExitNormal = fSaveState.sValidExitNormal;
699 fExiting = fSaveState.sExiting; 583 fExiting = fSaveState.sExiting;
700 fEntering = fSaveState.sEntering; 584 fEntering = fSaveState.sEntering;
701 585
702 fBlockedPhysicalVolume = fSaveState.spBlocke 586 fBlockedPhysicalVolume = fSaveState.spBlockedPhysicalVolume;
703 fBlockedReplicaNo = fSaveState.sBlockedRepli << 587 fBlockedReplicaNo = fSaveState.sBlockedReplicaNo,
704 588
705 fLastStepWasZero = (fSaveState.sLastStepWasZ << 589 fLastStepWasZero = fSaveState.sLastStepWasZero;
706 <<
707 fLocatedOutsideWorld = fSaveState.sLocatedOu <<
708 fLastLocatedPointLocal = fSaveState.sLastLoc <<
709 fEnteredDaughter = fSaveState.sEnteredDaught <<
710 fExitedMother = fSaveState.sExitedMother; <<
711 fWasLimitedByGeometry = fSaveState.sWasLimit <<
712 <<
713 // The 'expected' behaviour is to restore th <<
714 fPreviousSftOrigin = fSaveState.sPreviousSft <<
715 fPreviousSafety = fSaveState.sPreviousSafety <<
716 } 590 }
717 591
718 // ******************************************* 592 // ********************************************************************
719 // ComputeStep 593 // ComputeStep
720 // 594 //
721 // Computes the next geometric Step: intersect 595 // Computes the next geometric Step: intersections with current
722 // mother and `daughter' volumes. 596 // mother and `daughter' volumes.
723 // 597 //
724 // NOTE: 598 // NOTE:
725 // 599 //
726 // Flags on entry: 600 // Flags on entry:
727 // -------------- 601 // --------------
728 // fValidExitNormal - Normal of exited volume 602 // fValidExitNormal - Normal of exited volume is valid (convex, not a
729 // coincident boundary) 603 // coincident boundary)
730 // fExitNormal - Surface normal of exite 604 // fExitNormal - Surface normal of exited volume
731 // fExiting - True if have exited sol 605 // fExiting - True if have exited solid
732 // 606 //
733 // fBlockedPhysicalVolume - Ptr to exited volu 607 // fBlockedPhysicalVolume - Ptr to exited volume (or 0)
734 // fBlockedReplicaNo - Replication no of exite 608 // fBlockedReplicaNo - Replication no of exited volume
735 // fLastStepWasZero - True if last Step size << 609 // fLastStepWasZero - True if last Step size was zero.
736 // 610 //
737 // Flags on exit: 611 // Flags on exit:
738 // ------------- 612 // -------------
739 // fValidExitNormal - True if surface normal 613 // fValidExitNormal - True if surface normal of exited volume is valid
740 // fExitNormal - Surface normal of exite 614 // fExitNormal - Surface normal of exited volume rotated to mothers
741 // reference system 615 // reference system
742 // fExiting - True if exiting mother 616 // fExiting - True if exiting mother
743 // fEntering - True if entering `daugh 617 // fEntering - True if entering `daughter' volume (or replica)
744 // fBlockedPhysicalVolume - Ptr to candidate ( 618 // fBlockedPhysicalVolume - Ptr to candidate (entered) volume
745 // fBlockedReplicaNo - Replication no of candi 619 // fBlockedReplicaNo - Replication no of candidate (entered) volume
746 // fLastStepWasZero - True if this Step size << 620 // fLastStepWasZero - True if this Step size was zero.
747 // ******************************************* 621 // ********************************************************************
748 // 622 //
749 G4double G4Navigator::ComputeStep( const G4Thr << 623 G4double G4Navigator::ComputeStep( const G4ThreeVector &pGlobalpoint,
750 const G4Thr << 624 const G4ThreeVector &pDirection,
751 const G4dou 625 const G4double pCurrentProposedStepLength,
752 G4dou << 626 G4double &pNewSafety)
753 { 627 {
754 #ifdef G4DEBUG_NAVIGATION <<
755 static G4ThreadLocal G4int sNavCScalls = 0; <<
756 ++sNavCScalls; <<
757 #endif <<
758 <<
759 G4ThreeVector localDirection = ComputeLocalA 628 G4ThreeVector localDirection = ComputeLocalAxis(pDirection);
760 G4double Step = kInfinity; << 629 G4double Step = DBL_MAX;
761 G4VPhysicalVolume *motherPhysical = fHistor 630 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
762 G4LogicalVolume *motherLogical = motherPhysi 631 G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
763 632
764 // All state relating to exiting normals mus << 633 static G4int sNavCScalls=0;
765 // << 634 sNavCScalls++;
766 fExitNormalGlobalFrame = G4ThreeVector( 0., <<
767 // Reset value - to erase its memory <<
768 fChangedGrandMotherRefFrame = false; <<
769 // Reset - used for local exit normal <<
770 fGrandMotherExitNormal = G4ThreeVector( 0., <<
771 fCalculatedExitNormal = false; <<
772 // Reset for new step <<
773 635
774 #ifdef G4VERBOSE 636 #ifdef G4VERBOSE
>> 637 G4int oldcoutPrec= G4cout.precision(8);
>> 638 G4int oldcerrPrec= G4cerr.precision(10);
775 if( fVerbose > 0 ) 639 if( fVerbose > 0 )
776 { 640 {
777 G4cout << "*** G4Navigator::ComputeStep: * 641 G4cout << "*** G4Navigator::ComputeStep: ***" << G4endl;
778 G4cout << " Volume = " << motherPhysica 642 G4cout << " Volume = " << motherPhysical->GetName()
779 << " - Proposed step length = " << 643 << " - Proposed step length = " << pCurrentProposedStepLength
780 << G4endl; 644 << G4endl;
781 #ifdef G4DEBUG_NAVIGATION 645 #ifdef G4DEBUG_NAVIGATION
782 if( fVerbose >= 2 ) << 646 if( fVerbose >= 4 )
783 { 647 {
784 G4cout << " Called with the arguments: 648 G4cout << " Called with the arguments: " << G4endl
785 << " Globalpoint = " << std::set 649 << " Globalpoint = " << std::setw(25) << pGlobalpoint << G4endl
786 << " Direction = " << std::set 650 << " Direction = " << std::setw(25) << pDirection << G4endl;
787 if( fVerbose >= 4 ) << 651 G4cout << " ---- Upon entering :" << G4endl;
788 { << 652 PrintState();
789 G4cout << " ---- Upon entering : Stat <<
790 PrintState(); <<
791 } <<
792 } 653 }
793 #endif 654 #endif
794 } 655 }
>> 656
>> 657 static G4double fAccuracyForWarning = kCarTolerance,
>> 658 fAccuracyForException = 1000*kCarTolerance;
795 #endif 659 #endif
796 660
797 G4ThreeVector newLocalPoint = ComputeLocalPo 661 G4ThreeVector newLocalPoint = ComputeLocalPoint(pGlobalpoint);
798 <<
799 if( newLocalPoint != fLastLocatedPointLocal 662 if( newLocalPoint != fLastLocatedPointLocal )
800 { 663 {
801 // Check whether the relocation is within 664 // Check whether the relocation is within safety
802 // 665 //
803 G4ThreeVector oldLocalPoint = fLastLocated 666 G4ThreeVector oldLocalPoint = fLastLocatedPointLocal;
804 G4double moveLenSq = (newLocalPoint-oldLoc 667 G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2();
805 668
806 if ( moveLenSq >= fSqTol ) << 669 if ( moveLenSq >= kCarTolerance*kCarTolerance )
807 { 670 {
808 #ifdef G4VERBOSE 671 #ifdef G4VERBOSE
809 ComputeStepLog(pGlobalpoint, moveLenSq); << 672 // The following checks only make sense if the move is larger
>> 673 // than the tolerance.
>> 674 //
>> 675 G4ThreeVector OriginalGlobalpoint =
>> 676 fHistory.GetTopTransform().Inverse().
>> 677 TransformPoint(fLastLocatedPointLocal);
>> 678
>> 679 G4double shiftOriginSafSq = (fPreviousSftOrigin-pGlobalpoint).mag2();
>> 680
>> 681 // Check that the starting point of this step is
>> 682 // within the isotropic safety sphere of the last point
>> 683 // to a accuracy/precision given by fAccuracyForWarning.
>> 684 // If so give warning.
>> 685 // If it fails by more than fAccuracyForException exit with error.
>> 686 //
>> 687 if( shiftOriginSafSq >= sqr(fPreviousSafety) )
>> 688 {
>> 689 G4double shiftOrigin = std::sqrt(shiftOriginSafSq);
>> 690 G4double diffShiftSaf = shiftOrigin - fPreviousSafety;
>> 691
>> 692 if( diffShiftSaf > fAccuracyForWarning )
>> 693 {
>> 694 G4Exception("G4Navigator::ComputeStep()",
>> 695 "UnexpectedPositionShift", JustWarning,
>> 696 "Accuracy ERROR or slightly inaccurate position shift.");
>> 697 G4cerr << " The Step's starting point has moved "
>> 698 << std::sqrt(moveLenSq)/mm << " mm " << G4endl
>> 699 << " since the last call to a Locate method." << G4endl;
>> 700 G4cerr << " This has resulted in moving "
>> 701 << shiftOrigin/mm << " mm "
>> 702 << " from the last point at which the safety "
>> 703 << " was calculated " << G4endl;
>> 704 G4cerr << " which is more than the computed safety= "
>> 705 << fPreviousSafety/mm << " mm at that point." << G4endl;
>> 706 G4cerr << " This difference is "
>> 707 << diffShiftSaf/mm << " mm." << G4endl
>> 708 << " The tolerated accuracy is "
>> 709 << fAccuracyForException/mm << " mm." << G4endl;
>> 710
>> 711 static G4int warnNow = 0;
>> 712 if( ((++warnNow % 100) == 1) )
>> 713 {
>> 714 G4cerr << " This problem can be due to either " << G4endl;
>> 715 G4cerr << " - a process that has proposed a displacement"
>> 716 << " larger than the current safety , or" << G4endl;
>> 717 G4cerr << " - inaccuracy in the computation of the safety"
>> 718 << G4endl;
>> 719 G4cerr << " We suggest that you " << G4endl
>> 720 << " - find i) what particle is being tracked, and "
>> 721 << " ii) through what part of your geometry " << G4endl
>> 722 << " for example by re-running this event with "
>> 723 << G4endl
>> 724 << " /tracking/verbose 1 " << G4endl
>> 725 << " - check which processes you declare for"
>> 726 << " this particle (and look at non-standard ones)"
>> 727 << G4endl
>> 728 << " - in case, create a detailed logfile"
>> 729 << " of this event using:" << G4endl
>> 730 << " /tracking/verbose 6 "
>> 731 << G4endl;
>> 732 }
>> 733 }
>> 734 #ifdef G4DEBUG_NAVIGATION
>> 735 else
>> 736 {
>> 737 G4cerr << "WARNING - G4Navigator::ComputeStep()" << G4endl
>> 738 << " The Step's starting point has moved "
>> 739 << std::sqrt(moveLenSq) << "," << G4endl
>> 740 << " which has taken it to the limit of"
>> 741 << " the current safety. " << G4endl;
>> 742 }
810 #endif 743 #endif
>> 744 }
>> 745 G4double safetyPlus = fPreviousSafety + fAccuracyForException;
>> 746 if ( shiftOriginSafSq > sqr(safetyPlus) )
>> 747 {
>> 748 G4cerr << "ERROR - G4Navigator::ComputeStep()" << G4endl
>> 749 << " Position has shifted considerably without"
>> 750 << " notifying the navigator !" << G4endl
>> 751 << " Tolerated safety: " << safetyPlus << G4endl
>> 752 << " Computed shift : " << shiftOriginSafSq << G4endl;
>> 753 G4Exception("G4Navigator::ComputeStep()",
>> 754 "SignificantPositionShift", JustWarning,
>> 755 "May lead to a crash or unreliable results.");
>> 756 }
>> 757 #endif // end G4VERBOSE
>> 758
811 // Relocate the point within the same vo 759 // Relocate the point within the same volume
812 // 760 //
813 LocateGlobalPointWithinVolume( pGlobalpo 761 LocateGlobalPointWithinVolume( pGlobalpoint );
814 } 762 }
815 } 763 }
816 if ( fHistory.GetTopVolumeType()!=kReplica ) 764 if ( fHistory.GetTopVolumeType()!=kReplica )
817 { 765 {
818 switch( CharacteriseDaughters(motherLogica 766 switch( CharacteriseDaughters(motherLogical) )
819 { 767 {
820 case kNormal: 768 case kNormal:
821 if ( motherLogical->GetVoxelHeader() ! << 769 if ( motherLogical->GetVoxelHeader() )
822 { 770 {
823 Step = GetVoxelNavigator().ComputeSt << 771 Step = fvoxelNav.ComputeStep(fLastLocatedPointLocal,
824 localDi 772 localDirection,
825 pCurren 773 pCurrentProposedStepLength,
826 pNewSaf 774 pNewSafety,
827 fHistor 775 fHistory,
828 fValidE 776 fValidExitNormal,
829 fExitNo 777 fExitNormal,
830 fExitin 778 fExiting,
831 fEnteri 779 fEntering,
832 &fBlock 780 &fBlockedPhysicalVolume,
833 fBlocke 781 fBlockedReplicaNo);
834 782
835 } 783 }
836 else 784 else
837 { 785 {
838 if( motherPhysical->GetRegularStruct << 786 if( motherPhysical->GetRegularStructureId() != 1 )
839 { 787 {
840 Step = fnormalNav.ComputeStep(fLas 788 Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
841 loca 789 localDirection,
842 pCur 790 pCurrentProposedStepLength,
843 pNew 791 pNewSafety,
844 fHis 792 fHistory,
845 fVal 793 fValidExitNormal,
846 fExi 794 fExitNormal,
847 fExi 795 fExiting,
848 fEnt 796 fEntering,
849 &fBl 797 &fBlockedPhysicalVolume,
850 fBlo 798 fBlockedReplicaNo);
851 } 799 }
852 else // Regular (non-voxelised) str 800 else // Regular (non-voxelised) structure
853 { 801 {
854 LocateGlobalPointAndSetup( pGlobal 802 LocateGlobalPointAndSetup( pGlobalpoint, &pDirection, true, true );
855 // 803 //
856 // if physical process limits the 804 // if physical process limits the step, the voxel will not be the
857 // one given by ComputeStepSkippin 805 // one given by ComputeStepSkippingEqualMaterials() and the local
858 // point will be wrongly calculate 806 // point will be wrongly calculated.
859 807
860 // There is a problem: when msc li 808 // There is a problem: when msc limits the step and the point is
861 // assigned wrongly to phantom in 809 // assigned wrongly to phantom in previous step (while it is out
862 // of the container volume). Then 810 // of the container volume). Then LocateGlobalPointAndSetup() has
863 // reset the history topvolume to 811 // reset the history topvolume to world.
864 // 812 //
865 if(fHistory.GetTopVolume()->GetReg << 813 if(fHistory.GetTopVolume()->GetRegularStructureId() != 1 )
866 { 814 {
867 G4Exception("G4Navigator::Comput 815 G4Exception("G4Navigator::ComputeStep()",
868 "GeomNav1001", JustW << 816 "Bad-location-of-point", JustWarning,
869 "Point is relocated in voxels, 817 "Point is relocated in voxels, while it should be outside!");
870 Step = fnormalNav.ComputeStep(fL 818 Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
871 lo 819 localDirection,
872 pC 820 pCurrentProposedStepLength,
873 pN 821 pNewSafety,
874 fH 822 fHistory,
875 fV 823 fValidExitNormal,
876 fE 824 fExitNormal,
877 fE 825 fExiting,
878 fE 826 fEntering,
879 &f 827 &fBlockedPhysicalVolume,
880 fB 828 fBlockedReplicaNo);
881 } 829 }
882 else 830 else
883 { 831 {
884 Step = fregularNav. 832 Step = fregularNav.
885 ComputeStepSkippingEqualMat 833 ComputeStepSkippingEqualMaterials(fLastLocatedPointLocal,
886 834 localDirection,
887 835 pCurrentProposedStepLength,
888 836 pNewSafety,
889 837 fHistory,
890 838 fValidExitNormal,
891 839 fExitNormal,
892 840 fExiting,
893 841 fEntering,
894 842 &fBlockedPhysicalVolume,
895 843 fBlockedReplicaNo,
896 844 motherPhysical);
897 } 845 }
898 } 846 }
899 } 847 }
900 break; 848 break;
901 case kParameterised: 849 case kParameterised:
902 if( GetDaughtersRegularStructureId(mot 850 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
903 { 851 {
904 Step = fparamNav.ComputeStep(fLastLo 852 Step = fparamNav.ComputeStep(fLastLocatedPointLocal,
905 localDi 853 localDirection,
906 pCurren 854 pCurrentProposedStepLength,
907 pNewSaf 855 pNewSafety,
908 fHistor 856 fHistory,
909 fValidE 857 fValidExitNormal,
910 fExitNo 858 fExitNormal,
911 fExitin 859 fExiting,
912 fEnteri 860 fEntering,
913 &fBlock 861 &fBlockedPhysicalVolume,
914 fBlocke 862 fBlockedReplicaNo);
915 } 863 }
916 else // Regular structure 864 else // Regular structure
917 { 865 {
918 Step = fregularNav.ComputeStep(fLast 866 Step = fregularNav.ComputeStep(fLastLocatedPointLocal,
919 local 867 localDirection,
920 pCurr 868 pCurrentProposedStepLength,
921 pNewS 869 pNewSafety,
922 fHist 870 fHistory,
923 fVali 871 fValidExitNormal,
924 fExit 872 fExitNormal,
925 fExit 873 fExiting,
926 fEnte 874 fEntering,
927 &fBlo 875 &fBlockedPhysicalVolume,
928 fBloc 876 fBlockedReplicaNo);
929 } 877 }
930 break; 878 break;
931 case kReplica: 879 case kReplica:
932 G4Exception("G4Navigator::ComputeStep( << 880 G4Exception("G4Navigator::ComputeStep()", "NotApplicable",
933 FatalException, "Not appli 881 FatalException, "Not applicable for replicated volumes.");
934 break; 882 break;
935 case kExternal: <<
936 Step = fpExternalNav->ComputeStep(fLas <<
937 loca <<
938 pCur <<
939 pNew <<
940 fHis <<
941 fVal <<
942 fExi <<
943 fExi <<
944 fEnt <<
945 &fBl <<
946 fBlo <<
947 break; <<
948 } 883 }
949 } 884 }
950 else 885 else
951 { 886 {
952 // In the case of a replica, it must handl 887 // In the case of a replica, it must handle the exiting
953 // edge/corner problem by itself 888 // edge/corner problem by itself
954 // 889 //
955 fExiting = fExitedMother; << 890 G4bool exitingReplica = fExitedMother;
956 Step = freplicaNav.ComputeStep(pGlobalpoin 891 Step = freplicaNav.ComputeStep(pGlobalpoint,
957 pDirection, 892 pDirection,
958 fLastLocate 893 fLastLocatedPointLocal,
959 localDirect 894 localDirection,
960 pCurrentPro 895 pCurrentProposedStepLength,
961 pNewSafety, 896 pNewSafety,
962 fHistory, 897 fHistory,
963 fValidExitN 898 fValidExitNormal,
964 fCalculated <<
965 fExitNormal 899 fExitNormal,
966 fExiting, << 900 exitingReplica,
967 fEntering, 901 fEntering,
968 &fBlockedPh 902 &fBlockedPhysicalVolume,
969 fBlockedRep 903 fBlockedReplicaNo);
>> 904 fExiting= exitingReplica; // still ok to set it ??
970 } 905 }
971 906
972 // Remember last safety origin & value. 907 // Remember last safety origin & value.
973 // 908 //
974 fPreviousSftOrigin = pGlobalpoint; 909 fPreviousSftOrigin = pGlobalpoint;
975 fPreviousSafety = pNewSafety; 910 fPreviousSafety = pNewSafety;
976 911
977 // Count zero steps - one can occur due to c 912 // Count zero steps - one can occur due to changing momentum at a boundary
978 // - one, two (or a few) ca 913 // - one, two (or a few) can occur at common edges between
979 // volumes 914 // volumes
980 // - more than two is likel 915 // - more than two is likely a problem in the geometry
981 // description or the Nav 916 // description or the Navigation
982 917
983 // Rule of thumb: likely at an Edge if two c 918 // Rule of thumb: likely at an Edge if two consecutive steps are zero,
984 // because at least two candi 919 // because at least two candidate volumes must have been
985 // checked 920 // checked
986 // 921 //
987 fLocatedOnEdge = fLastStepWasZero && (Step 922 fLocatedOnEdge = fLastStepWasZero && (Step==0.0);
988 fLastStepWasZero = (Step<fMinStep); << 923 fLastStepWasZero = (Step==0.0);
989 if (fPushed) { fPushed = fLastStepWasZero; << 924 if (fPushed) fPushed = fLastStepWasZero;
990 925
991 // Handle large number of consecutive zero s 926 // Handle large number of consecutive zero steps
992 // 927 //
993 if ( fLastStepWasZero ) 928 if ( fLastStepWasZero )
994 { 929 {
995 ++fNumberZeroSteps; << 930 fNumberZeroSteps++;
996 <<
997 G4bool act = fNumberZeroSteps >= fActi <<
998 G4bool actAndReport = false; <<
999 G4bool abandon = fNumberZeroSteps >= fAban <<
1000 G4bool inform = false; <<
1001 #ifdef G4VERBOSE <<
1002 actAndReport = act && (!fPushed) && fWarn <<
1003 #endif <<
1004 #ifdef G4DEBUG_NAVIGATION 931 #ifdef G4DEBUG_NAVIGATION
1005 inform = fNumberZeroSteps > 1; << 932 if( fNumberZeroSteps > 1 )
1006 #endif <<
1007 <<
1008 if ( act || inform ) <<
1009 { 933 {
1010 if( act && !abandon ) << 934 G4cout << "G4Navigator::ComputeStep(): another zero step, # "
1011 { << 935 << fNumberZeroSteps
1012 // Act to recover this stuck track. P << 936 << " at " << pGlobalpoint
1013 // << 937 << " in volume " << motherPhysical->GetName()
1014 Step += 100*kCarTolerance; << 938 << " nav-comp-step calls # " << sNavCScalls
1015 fPushed = true; << 939 << G4endl;
1016 } << 940 }
1017 <<
1018 if( actAndReport || abandon || inform ) <<
1019 { <<
1020 std::ostringstream message; <<
1021 <<
1022 message.precision(16); <<
1023 message << "Stuck Track: potential ge <<
1024 << G4endl; <<
1025 message << " Track stuck, not moving <<
1026 << fNumberZeroSteps << " step <<
1027 << " Current phys volume: ' <<
1028 << "'" << G4endl <<
1029 << " - at position : " << p <<
1030 << " in direction: " << p <<
1031 << " (local position: " << <<
1032 << " (local direction: " < <<
1033 << " Previous phys volume: ' <<
1034 << ( fLastMotherPhys != nullp <<
1035 << "'" << G4endl << G4endl; <<
1036 if( actAndReport || abandon ) <<
1037 { <<
1038 message << " Likely geometry over <<
1039 << G4endl; <<
1040 } <<
1041 if( abandon ) // i.e. fNumberZeroStep <<
1042 { <<
1043 // Must kill this stuck track <<
1044 #ifdef G4VERBOSE <<
1045 if ( fWarnPush ) { CheckOverlapsIte <<
1046 #endif 941 #endif
1047 message << " Track *abandoned* due << 942 if( fNumberZeroSteps > fActionThreshold_NoZeroSteps-1 )
1048 << " Event aborted. " << G4 << 943 {
1049 G4Exception("G4Navigator::ComputeSt << 944 // Act to recover this stuck track. Pushing it along direction
1050 EventMustBeAborted, mes << 945 //
1051 } << 946 Step += 0.9*kCarTolerance;
1052 else <<
1053 { <<
1054 #ifdef G4VERBOSE 947 #ifdef G4VERBOSE
1055 if ( actAndReport ) // (!fPushed = << 948 if (!fPushed)
1056 { << 949 {
1057 message << " *** Trying to get << 950 G4cerr << "WARNING - G4Navigator::ComputeStep()" << G4endl
1058 << " - expanding step to << 951 << " Track stuck, not moving for "
1059 << " Potential ove << 952 << fNumberZeroSteps << " steps" << G4endl
1060 G4Exception("G4Navigator::Comput << 953 << " in volume -" << motherPhysical->GetName()
1061 JustWarning, message << 954 << "- at point " << pGlobalpoint << G4endl
1062 } << 955 << " direction: " << pDirection << "." << G4endl
1063 #endif << 956 << " Potential geometry or navigation problem !"
1064 #ifdef G4DEBUG_NAVIGATION << 957 << G4endl
1065 else << 958 << " Trying pushing it of " << Step << " mm ..."
1066 { << 959 << G4endl;
1067 if( fNumberZeroSteps > 1 ) << 960 }
1068 { <<
1069 message << ", nav-comp-step ca <<
1070 << ", Step= " << Step <<
1071 G4cout << message.str(); <<
1072 } <<
1073 } <<
1074 #endif 961 #endif
1075 } // end of else if ( abandon ) << 962 fPushed = true;
1076 } // end of if( actAndReport || abandon << 963 }
1077 } // end of if ( act || inform ) << 964 if( fNumberZeroSteps > fAbandonThreshold_NoZeroSteps-1 )
>> 965 {
>> 966 // Must kill this stuck track
>> 967 //
>> 968 G4cerr << "ERROR - G4Navigator::ComputeStep()" << G4endl
>> 969 << " Track stuck, not moving for "
>> 970 << fNumberZeroSteps << " steps" << G4endl
>> 971 << " in volume -" << motherPhysical->GetName()
>> 972 << "- at point " << pGlobalpoint << G4endl
>> 973 << " direction: " << pDirection << "." << G4endl;
>> 974 motherPhysical->CheckOverlaps(5000, false);
>> 975 G4Exception("G4Navigator::ComputeStep()",
>> 976 "StuckTrack", EventMustBeAborted,
>> 977 "Stuck Track: potential geometry or navigation problem.");
>> 978 }
1078 } 979 }
1079 else 980 else
1080 { 981 {
1081 if (!fPushed) { fNumberZeroSteps = 0; } << 982 if (!fPushed) fNumberZeroSteps = 0;
1082 } 983 }
1083 fLastMotherPhys = motherPhysical; <<
1084 984
1085 fEnteredDaughter = fEntering; // I expect 985 fEnteredDaughter = fEntering; // I expect to enter a volume in this Step
1086 fExitedMother = fExiting; 986 fExitedMother = fExiting;
1087 987
1088 fStepEndPoint = pGlobalpoint <<
1089 + std::min(Step,pCurrentPropo <<
1090 fLastStepEndPointLocal = fLastLocatedPointL <<
1091 <<
1092 if( fExiting ) 988 if( fExiting )
1093 { 989 {
1094 #ifdef G4DEBUG_NAVIGATION 990 #ifdef G4DEBUG_NAVIGATION
1095 if( fVerbose > 2 ) 991 if( fVerbose > 2 )
1096 { 992 {
1097 G4cout << " At G4Nav CompStep End - if( << 993 G4cout << " At G4Nav CompStep End - if(exiting) - fExiting= " << fExiting
1098 << " fValidExitNormal = " << fVa 994 << " fValidExitNormal = " << fValidExitNormal << G4endl;
1099 G4cout << " fExitNormal= " << fExitNorm 995 G4cout << " fExitNormal= " << fExitNormal << G4endl;
1100 } 996 }
1101 #endif 997 #endif
1102 998
1103 if ( fValidExitNormal || fCalculatedExitN << 999 if(fValidExitNormal)
1104 { 1000 {
1105 // Convention: fExitNormal is in the 'g 1001 // Convention: fExitNormal is in the 'grand-mother' coordinate system
1106 fGrandMotherExitNormal = fExitNormal; << 1002 //
>> 1003 fGrandMotherExitNormal= fExitNormal;
1107 } 1004 }
1108 else 1005 else
1109 { 1006 {
1110 // We must calculate the normal anyway 1007 // We must calculate the normal anyway (in order to have it if requested)
1111 // 1008 //
1112 G4ThreeVector finalLocalPoint = fLastLo << 1009 G4ThreeVector finalLocalPoint =
1113 + localDi << 1010 fLastLocatedPointLocal + localDirection*Step;
1114 <<
1115 if ( fHistory.GetTopVolumeType() != kR <<
1116 { <<
1117 // Find normal in the 'mother' coordi <<
1118 // <<
1119 G4ThreeVector exitNormalMotherFrame= <<
1120 motherLogical->GetSolid()->Surface <<
1121 <<
1122 // Transform it to the 'grand-mother' <<
1123 // <<
1124 const G4RotationMatrix* mRot = mother <<
1125 if( mRot != nullptr ) <<
1126 { <<
1127 fChangedGrandMotherRefFrame = true; <<
1128 fGrandMotherExitNormal = (*mRot).in <<
1129 } <<
1130 else <<
1131 { <<
1132 fGrandMotherExitNormal = exitNormal <<
1133 } <<
1134 <<
1135 // Do not set fValidExitNormal -- thi <<
1136 // that the solid is convex! <<
1137 } <<
1138 else <<
1139 { <<
1140 fCalculatedExitNormal = false; <<
1141 // <<
1142 // Nothing can be done at this stage <<
1143 // Replica Navigation must have calcu <<
1144 // already. <<
1145 // Cases: mother is not convex, and e <<
1146 <<
1147 #ifdef G4DEBUG_NAVIGATION <<
1148 G4ExceptionDescription desc; <<
1149 <<
1150 desc << "Problem in ComputeStep: Rep <<
1151 << " valid exit Normal. " << G4e <<
1152 desc << " Do not know how calculate i <<
1153 desc << " Location = " << finalLo <<
1154 desc << " Volume name = " << motherP <<
1155 << " copy/replica No = " << mot <<
1156 G4Exception("G4Navigator::ComputeStep <<
1157 JustWarning, desc, "Norma <<
1158 #endif <<
1159 } <<
1160 } <<
1161 1011
1162 if ( fHistory.GetTopVolumeType() != kRepl << 1012 // Now fGrandMotherExitNormal is in the 'grand-mother' coordinate system
1163 { << 1013 //
1164 fCalculatedExitNormal = true; << 1014 fGrandMotherExitNormal =
1165 } << 1015 motherLogical->GetSolid()->SurfaceNormal(finalLocalPoint);
1166 1016
1167 // Now transform it to the global referen << 1017 const G4RotationMatrix* mRot = motherPhysical->GetRotation();
1168 // << 1018 if( mRot )
1169 if( fValidExitNormal || fCalculatedExitNo << 1019 {
1170 { << 1020 fGrandMotherExitNormal *= (*mRot);
1171 auto depth = (G4int)fHistory.GetDepth( <<
1172 if( depth > 0 ) <<
1173 { <<
1174 fExitNormalGlobalFrame = fHistory.Get <<
1175 .InverseTrans <<
1176 } <<
1177 else <<
1178 { <<
1179 fExitNormalGlobalFrame = fGrandMother <<
1180 } 1021 }
1181 } 1022 }
1182 else <<
1183 { <<
1184 fExitNormalGlobalFrame = G4ThreeVector( <<
1185 } <<
1186 } 1023 }
>> 1024 fStepEndPoint= pGlobalpoint+Step*pDirection;
1187 1025
1188 if( (Step == pCurrentProposedStepLength) && 1026 if( (Step == pCurrentProposedStepLength) && (!fExiting) && (!fEntering) )
1189 { 1027 {
1190 // This if Step is not really limited by 1028 // This if Step is not really limited by the geometry.
1191 // The Navigator is obliged to return "in 1029 // The Navigator is obliged to return "infinity"
1192 // 1030 //
1193 Step = kInfinity; 1031 Step = kInfinity;
1194 } 1032 }
1195 1033
1196 #ifdef G4VERBOSE 1034 #ifdef G4VERBOSE
1197 if( fVerbose > 1 ) 1035 if( fVerbose > 1 )
1198 { 1036 {
1199 if( fVerbose >= 4 ) 1037 if( fVerbose >= 4 )
1200 { 1038 {
1201 G4cout << " ----- Upon exiting :" << 1039 G4cout << " ----- Upon exiting :" << G4endl;
1202 PrintState(); 1040 PrintState();
1203 } 1041 }
1204 G4cout << " Returned step= " << Step; << 1042 G4cout <<" Returned step = " << Step << G4endl;
1205 if( fVerbose > 5 ) { G4cout << G4endl; } <<
1206 if( Step == kInfinity ) 1043 if( Step == kInfinity )
1207 { 1044 {
1208 G4cout << " Requested step= " << pCurr << 1045 G4cout << " Original proposed step = "
1209 if( fVerbose > 5) { G4cout << G4endl; << 1046 << pCurrentProposedStepLength << G4endl;
1210 } 1047 }
1211 G4cout << " Safety = " << pNewSafety << << 1048 G4cout << " Safety = " << pNewSafety << G4endl;
1212 } 1049 }
>> 1050 G4cout.precision(oldcoutPrec);
>> 1051 G4cerr.precision(oldcerrPrec);
1213 #endif 1052 #endif
1214 1053
1215 fLastTriedStepComputation = true; <<
1216 <<
1217 return Step; 1054 return Step;
1218 } 1055 }
1219 1056
1220 // ****************************************** 1057 // ********************************************************************
1221 // CheckNextStep 1058 // CheckNextStep
1222 // 1059 //
1223 // Compute the step without altering the navi 1060 // Compute the step without altering the navigator state
1224 // ****************************************** 1061 // ********************************************************************
1225 // 1062 //
1226 G4double G4Navigator::CheckNextStep( const G4 1063 G4double G4Navigator::CheckNextStep( const G4ThreeVector& pGlobalpoint,
1227 const G4 1064 const G4ThreeVector& pDirection,
1228 const G4 1065 const G4double pCurrentProposedStepLength,
1229 G4 1066 G4double& pNewSafety)
1230 { 1067 {
1231 G4double step; 1068 G4double step;
1232 1069
1233 // Save the state, for this parasitic call 1070 // Save the state, for this parasitic call
1234 // 1071 //
1235 SetSavedState(); 1072 SetSavedState();
1236 1073
1237 step = ComputeStep ( pGlobalpoint, 1074 step = ComputeStep ( pGlobalpoint,
1238 pDirection, 1075 pDirection,
1239 pCurrentProposedStepLe 1076 pCurrentProposedStepLength,
1240 pNewSafety ); 1077 pNewSafety );
1241 1078
1242 // It is a parasitic call, so attempt to re << 1079 // If a parasitic call, then attempt to restore the key parts of the state
1243 // 1080 //
1244 RestoreSavedState(); 1081 RestoreSavedState();
1245 // NOTE: the state of the current subnaviga << 1082
1246 // ***> TODO: restore subnavigator state <<
1247 // if( last_located) Need <<
1248 // if( last_computed step) Need <<
1249 <<
1250 return step; 1083 return step;
1251 } 1084 }
1252 1085
1253 // ****************************************** 1086 // ********************************************************************
1254 // ResetState 1087 // ResetState
1255 // 1088 //
1256 // Resets stack and minimum of navigator stat 1089 // Resets stack and minimum of navigator state `machine'
1257 // ****************************************** 1090 // ********************************************************************
1258 // 1091 //
1259 void G4Navigator::ResetState() 1092 void G4Navigator::ResetState()
1260 { 1093 {
1261 fWasLimitedByGeometry = false; 1094 fWasLimitedByGeometry = false;
1262 fEntering = false; 1095 fEntering = false;
1263 fExiting = false; 1096 fExiting = false;
1264 fLocatedOnEdge = false; 1097 fLocatedOnEdge = false;
1265 fLastStepWasZero = false; 1098 fLastStepWasZero = false;
1266 fEnteredDaughter = false; 1099 fEnteredDaughter = false;
1267 fExitedMother = false; 1100 fExitedMother = false;
1268 fPushed = false; 1101 fPushed = false;
1269 1102
1270 fValidExitNormal = false; 1103 fValidExitNormal = false;
1271 fChangedGrandMotherRefFrame = false; <<
1272 fCalculatedExitNormal = false; <<
1273 <<
1274 fExitNormal = G4ThreeVector(0,0, 1104 fExitNormal = G4ThreeVector(0,0,0);
1275 fGrandMotherExitNormal = G4ThreeVector(0,0, <<
1276 fExitNormalGlobalFrame = G4ThreeVector(0,0, <<
1277 1105
1278 fPreviousSftOrigin = G4ThreeVector(0,0, 1106 fPreviousSftOrigin = G4ThreeVector(0,0,0);
1279 fPreviousSafety = 0.0; 1107 fPreviousSafety = 0.0;
1280 1108
1281 fNumberZeroSteps = 0; 1109 fNumberZeroSteps = 0;
1282 << 1110
1283 fBlockedPhysicalVolume = nullptr; << 1111 fBlockedPhysicalVolume = 0;
1284 fBlockedReplicaNo = -1; 1112 fBlockedReplicaNo = -1;
1285 1113
1286 fLastLocatedPointLocal = G4ThreeVector( kIn << 1114 fLastLocatedPointLocal = G4ThreeVector( DBL_MAX, -DBL_MAX, 0.0 );
1287 fLocatedOutsideWorld = false; 1115 fLocatedOutsideWorld = false;
1288 <<
1289 fLastMotherPhys = nullptr; <<
1290 } 1116 }
1291 1117
1292 // ****************************************** 1118 // ********************************************************************
1293 // SetupHierarchy 1119 // SetupHierarchy
1294 // 1120 //
1295 // Renavigates & resets hierarchy described b 1121 // Renavigates & resets hierarchy described by current history
1296 // o Reset volumes 1122 // o Reset volumes
1297 // o Recompute transforms and/or solids of re 1123 // o Recompute transforms and/or solids of replicated/parameterised volumes
1298 // ****************************************** 1124 // ********************************************************************
1299 // 1125 //
1300 void G4Navigator::SetupHierarchy() 1126 void G4Navigator::SetupHierarchy()
1301 { 1127 {
1302 const auto depth = (G4int)fHistory.GetDept << 1128 G4int i;
1303 for ( auto i = 1; i <= depth; ++i ) << 1129 const G4int cdepth = fHistory.GetDepth();
>> 1130 G4VPhysicalVolume *mother, *current;
>> 1131 G4VSolid *pSolid;
>> 1132 G4VPVParameterisation *pParam;
>> 1133
>> 1134 mother = fHistory.GetVolume(0);
>> 1135 for ( i=1; i<=cdepth; i++ )
1304 { 1136 {
>> 1137 current = fHistory.GetVolume(i);
1305 switch ( fHistory.GetVolumeType(i) ) 1138 switch ( fHistory.GetVolumeType(i) )
1306 { 1139 {
1307 case kNormal: 1140 case kNormal:
1308 case kExternal: <<
1309 break; 1141 break;
1310 case kReplica: 1142 case kReplica:
1311 freplicaNav.ComputeTransformation(fHi << 1143 freplicaNav.ComputeTransformation(fHistory.GetReplicaNo(i), current);
1312 break; 1144 break;
1313 case kParameterised: 1145 case kParameterised:
1314 G4VPhysicalVolume* current = fHistory << 1146 G4int replicaNo;
1315 G4int replicaNo = fHistory.GetReplica << 1147 pParam = current->GetParameterisation();
1316 G4VPVParameterisation* pParam = curre << 1148 replicaNo = fHistory.GetReplicaNo(i);
1317 G4VSolid* pSolid = pParam->ComputeSol << 1149 pSolid = pParam->ComputeSolid(replicaNo, current);
1318 1150
1319 // Set up dimensions & transform in s 1151 // Set up dimensions & transform in solid/physical volume
1320 // 1152 //
1321 pSolid->ComputeDimensions(pParam, rep 1153 pSolid->ComputeDimensions(pParam, replicaNo, current);
1322 pParam->ComputeTransformation(replica 1154 pParam->ComputeTransformation(replicaNo, current);
1323 1155
1324 G4TouchableHistory* pTouchable = null << 1156 G4TouchableHistory touchable( fHistory );
1325 if( pParam->IsNested() ) << 1157 touchable.MoveUpHistory(); // move up to the parent level
1326 { << 1158
1327 pTouchable= new G4TouchableHistory( <<
1328 pTouchable->MoveUpHistory(); // Mov <<
1329 // Adequate only if Nested at the <<
1330 // To extend to other cases: <<
1331 // pTouchable->MoveUpHistory(cdep <<
1332 // Move to the parent level of *C <<
1333 // Could replace this line and co <<
1334 // c-tor for History(levels to dr <<
1335 } <<
1336 // Set up the correct solid and mater 1159 // Set up the correct solid and material in Logical Volume
1337 // 1160 //
1338 G4LogicalVolume* pLogical = current-> << 1161 G4LogicalVolume *pLogical = current->GetLogicalVolume();
1339 pLogical->SetSolid( pSolid ); 1162 pLogical->SetSolid( pSolid );
1340 pLogical->UpdateMaterial( pParam -> 1163 pLogical->UpdateMaterial( pParam ->
1341 ComputeMaterial(replicaNo, current, << 1164 ComputeMaterial(replicaNo, current, &touchable) );
1342 delete pTouchable; <<
1343 break; 1165 break;
1344 } 1166 }
>> 1167 mother = current;
1345 } 1168 }
1346 } 1169 }
1347 1170
1348 // ****************************************** 1171 // ********************************************************************
1349 // GetLocalExitNormal 1172 // GetLocalExitNormal
1350 // 1173 //
1351 // Obtains the Normal vector to a surface (in 1174 // Obtains the Normal vector to a surface (in local coordinates)
1352 // pointing out of previous volume and into c 1175 // pointing out of previous volume and into current volume
1353 // ****************************************** 1176 // ********************************************************************
1354 // 1177 //
1355 G4ThreeVector G4Navigator::GetLocalExitNormal 1178 G4ThreeVector G4Navigator::GetLocalExitNormal( G4bool* valid )
1356 { 1179 {
1357 G4ThreeVector ExitNormal(0.,0.,0.); << 1180 G4ThreeVector ExitNormal(0.,0.,0.);
1358 G4VSolid* currentSolid = nullptr; <<
1359 G4LogicalVolume* candidateLogical; <<
1360 1181
1361 if ( fLastTriedStepComputation ) << 1182 if ( EnteredDaughterVolume() )
1362 { 1183 {
1363 // use fLastLocatedPointLocal and next ca << 1184 ExitNormal= -(fHistory.GetTopVolume()->GetLogicalVolume()->
1364 // << 1185 GetSolid()->SurfaceNormal(fLastLocatedPointLocal));
1365 G4ThreeVector nextSolidExitNormal(0.,0.,0 << 1186 *valid = true;
1366 <<
1367 if( fEntering && (fBlockedPhysicalVolume! <<
1368 { <<
1369 candidateLogical = fBlockedPhysicalVolu <<
1370 if( candidateLogical != nullptr ) <<
1371 { <<
1372 // fLastStepEndPointLocal is in the c <<
1373 // we need it in the daughter's coord <<
1374 <<
1375 // The following code should also wor <<
1376 { <<
1377 // First transform fLastLocatedPoin <<
1378 // coordinates <<
1379 // <<
1380 G4AffineTransform MotherToDaughterT <<
1381 GetMotherToDaughterTransform( fBl <<
1382 fBl <<
1383 Vol <<
1384 G4ThreeVector daughterPointOwnLocal <<
1385 MotherToDaughterTransform.Transfo <<
1386 <<
1387 // OK if it is a parameterised volu <<
1388 // <<
1389 EInside inSideIt; <<
1390 G4bool onSurface; <<
1391 G4double safety = -1.0; <<
1392 currentSolid = candidateLogical->Ge <<
1393 inSideIt = currentSolid->Inside(dau <<
1394 onSurface = (inSideIt == kSurface); <<
1395 if( !onSurface ) <<
1396 { <<
1397 if( inSideIt == kOutside ) <<
1398 { <<
1399 safety = (currentSolid->Distanc <<
1400 onSurface = safety < 100.0 * kC <<
1401 } <<
1402 else if (inSideIt == kInside ) <<
1403 { <<
1404 safety = (currentSolid->Distanc <<
1405 onSurface = safety < 100.0 * kC <<
1406 } <<
1407 } <<
1408 <<
1409 if( onSurface ) <<
1410 { <<
1411 nextSolidExitNormal = <<
1412 currentSolid->SurfaceNormal(dau <<
1413 <<
1414 // Entering the solid ==> opposit <<
1415 // <<
1416 // First flip ( ExitNormal = -nex <<
1417 // and then rotate the the norma <<
1418 ExitNormal = MotherToDaughterTran <<
1419 .InverseTransformAxis <<
1420 fCalculatedExitNormal = true; <<
1421 } <<
1422 else <<
1423 { <<
1424 #ifdef G4VERBOSE <<
1425 if(( fVerbose == 1 ) && ( fCheck <<
1426 { <<
1427 std::ostringstream message; <<
1428 message << "Point not on surfac <<
1429 << " Point = <<
1430 << daughterPointOwnLoca <<
1431 << " Physical volume = <<
1432 << fBlockedPhysicalVolu <<
1433 << " Logical volume = <<
1434 << candidateLogical->Ge <<
1435 << " Solid = <<
1436 << " Type = <<
1437 << currentSolid->GetEnt <<
1438 << *currentSolid << G4e <<
1439 if( inSideIt == kOutside ) <<
1440 { <<
1441 message << "Point is Outside. <<
1442 << " Safety (from ou <<
1443 } <<
1444 else // if( inSideIt == kInside <<
1445 { <<
1446 message << "Point is Inside. <<
1447 << " Safety (from in <<
1448 } <<
1449 G4Exception("G4Navigator::GetLo <<
1450 JustWarning, messag <<
1451 } <<
1452 #endif <<
1453 } <<
1454 *valid = onSurface; // was =tru <<
1455 } <<
1456 } <<
1457 } <<
1458 else if ( fExiting ) <<
1459 { <<
1460 ExitNormal = fGrandMotherExitNormal; <<
1461 *valid = true; <<
1462 fCalculatedExitNormal = true; // Shoul <<
1463 } <<
1464 else // i.e. ( fBlockedPhysicalVolume = <<
1465 { <<
1466 *valid = false; <<
1467 G4Exception("G4Navigator::GetLocalExitN <<
1468 "GeomNav0003", JustWarning, <<
1469 "Incorrect call to GetLocal <<
1470 } <<
1471 } 1187 }
1472 else // ( ! fLastTriedStepComputation ) i. << 1188 else
1473 { 1189 {
1474 if ( EnteredDaughterVolume() ) << 1190 if( fExitedMother )
1475 { 1191 {
1476 G4VSolid* daughterSolid = fHistory.GetT << 1192 ExitNormal = fGrandMotherExitNormal;
1477 <<
1478 ExitNormal = -(daughterSolid->SurfaceNo <<
1479 if( std::fabs(ExitNormal.mag2()-1.0 ) > <<
1480 { <<
1481 G4ExceptionDescription desc; <<
1482 desc << " Parameters of solid: " << * <<
1483 << " Point for surface = " << fL <<
1484 G4Exception("G4Navigator::GetLocalExi <<
1485 "GeomNav0003", FatalExcep <<
1486 "Surface Normal returned <<
1487 } <<
1488 fCalculatedExitNormal = true; <<
1489 *valid = true; 1193 *valid = true;
1490 } 1194 }
1491 else 1195 else
1492 { 1196 {
1493 if( fExitedMother ) << 1197 // We are not at a boundary.
1494 { << 1198 // ExitNormal remains (0,0,0)
1495 ExitNormal = fGrandMotherExitNormal; << 1199 //
1496 *valid = true; << 1200 *valid = false;
1497 fCalculatedExitNormal = true; <<
1498 } <<
1499 else // We are not at a boundary. Exit <<
1500 { <<
1501 *valid = false; <<
1502 fCalculatedExitNormal = false; <<
1503 G4ExceptionDescription message; <<
1504 message << "Function called when *NOT <<
1505 message << "Exit Normal not calculate <<
1506 G4Exception("G4Navigator::GetLocalExi <<
1507 "GeomNav0003", JustWarnin <<
1508 } <<
1509 } 1201 }
1510 } 1202 }
1511 return ExitNormal; 1203 return ExitNormal;
1512 } 1204 }
1513 1205
1514 // ****************************************** 1206 // ********************************************************************
1515 // GetMotherToDaughterTransform << 1207 // ComputeSafety
1516 // 1208 //
1517 // Obtains the mother to daughter affine tran << 1209 // It assumes that it will be
>> 1210 // i) called at the Point in the same volume as the EndPoint of the
>> 1211 // ComputeStep.
>> 1212 // ii) after (or at the end of) ComputeStep OR after the relocation.
1518 // ****************************************** 1213 // ********************************************************************
1519 // 1214 //
1520 G4AffineTransform << 1215 G4double G4Navigator::ComputeSafety( const G4ThreeVector &pGlobalpoint,
1521 G4Navigator::GetMotherToDaughterTransform( G4 << 1216 const G4double pMaxLength)
1522 G4 <<
1523 EV <<
1524 { 1217 {
1525 switch (enteringVolumeType) << 1218 G4double newSafety = 0.0;
1526 { <<
1527 case kNormal: // Nothing is needed to pr <<
1528 break; // It is stored already in <<
1529 case kReplica: // Sets the transform in t <<
1530 G4Exception("G4Navigator::GetMotherToDa <<
1531 "GeomNav0001", FatalExcepti <<
1532 "Method NOT Implemented yet <<
1533 break; <<
1534 case kParameterised: <<
1535 if( pEnteringPhysVol->GetRegularStructu <<
1536 { <<
1537 G4VPVParameterisation *pParam = <<
1538 pEnteringPhysVol->GetParameterisati <<
1539 G4VSolid* pSolid = <<
1540 pParam->ComputeSolid(enteringReplic <<
1541 pSolid->ComputeDimensions(pParam, ent <<
1542 <<
1543 // Sets the transform in the Paramete <<
1544 // <<
1545 pParam->ComputeTransformation(enterin <<
1546 <<
1547 // Set the correct solid and material <<
1548 // <<
1549 G4LogicalVolume* pLogical = pEntering <<
1550 pLogical->SetSolid( pSolid ); <<
1551 } <<
1552 break; <<
1553 case kExternal: <<
1554 // Expect that nothing is needed to pre <<
1555 // It is stored already in the physical <<
1556 break; <<
1557 } <<
1558 return G4AffineTransform(pEnteringPhysVol-> <<
1559 pEnteringPhysVol-> <<
1560 } <<
1561 1219
1562 <<
1563 // ****************************************** <<
1564 // GetLocalExitNormalAndCheck <<
1565 // <<
1566 // Obtains the Normal vector to a surface (in <<
1567 // pointing out of previous volume and into c <<
1568 // checks the current point against expected <<
1569 // ****************************************** <<
1570 // <<
1571 G4ThreeVector <<
1572 G4Navigator::GetLocalExitNormalAndCheck( <<
1573 #ifdef G4DEBUG_NAVIGATION <<
1574 const G4ThreeVecto <<
1575 #else <<
1576 const G4ThreeVecto <<
1577 #endif <<
1578 G4bool* pVal <<
1579 { <<
1580 #ifdef G4DEBUG_NAVIGATION 1220 #ifdef G4DEBUG_NAVIGATION
1581 // Check Current point against expected 'lo << 1221 G4int oldcoutPrec = G4cout.precision(8);
1582 // << 1222 if( fVerbose > 0 )
1583 if ( fLastTriedStepComputation ) <<
1584 { 1223 {
1585 G4ThreeVector ExpectedBoundaryPointLocal; << 1224 G4cout << "*** G4Navigator::ComputeSafety: ***" << G4endl
1586 << 1225 << " Called at point: " << pGlobalpoint << G4endl;
1587 const G4AffineTransform& GlobalToLocal = <<
1588 ExpectedBoundaryPointLocal = <<
1589 GlobalToLocal.TransformPoint( ExpectedB <<
1590 1226
1591 // Add here: Comparison against expected << 1227 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1592 // i.e. the endpoint of Comput << 1228 G4cout << " Volume = " << motherPhysical->GetName()
>> 1229 << " - Maximum length = " << pMaxLength << G4endl;
>> 1230 if( fVerbose >= 4 )
>> 1231 {
>> 1232 G4cout << " ----- Upon entering Compute Safety:" << G4endl;
>> 1233 PrintState();
>> 1234 }
1593 } 1235 }
1594 #endif 1236 #endif
1595 <<
1596 return GetLocalExitNormal( pValid ); <<
1597 } <<
1598 1237
1599 // ****************************************** << 1238 G4double distEndpointSq = (pGlobalpoint-fStepEndPoint).mag2();
1600 // GetGlobalExitNormal << 1239 G4bool stayedOnEndpoint = distEndpointSq < kCarTolerance*kCarTolerance;
1601 // << 1240 G4bool endpointOnSurface = fEnteredDaughter || fExitedMother;
1602 // Obtains the Normal vector to a surface (in << 1241
1603 // pointing out of previous volume and into c << 1242 if( !(endpointOnSurface && stayedOnEndpoint) )
1604 // ****************************************** <<
1605 // <<
1606 G4ThreeVector <<
1607 G4Navigator::GetGlobalExitNormal(const G4Thre <<
1608 G4bool <<
1609 { <<
1610 G4bool validNormal; <<
1611 G4ThreeVector localNormal, globalNormal; <<
1612 <<
1613 G4bool usingStored = fCalculatedExitNormal <<
1614 ( fLastTriedStepComputation && fExitin <<
1615 || <<
1616 ( !fLastTriedStepComputation <<
1617 && (IntersectPointGlobal-fStepEndPo <<
1618 // Calculated it 'just' before & t <<
1619 // but it did not move position <<
1620 <<
1621 if( usingStored ) <<
1622 { 1243 {
1623 // This was computed in last call to Comp << 1244 // Pseudo-relocate to this point (updates voxel information only)
1624 // and only if it arrived at boundary <<
1625 // 1245 //
1626 globalNormal = fExitNormalGlobalFrame; << 1246 LocateGlobalPointWithinVolume( pGlobalpoint );
1627 G4double normMag2 = globalNormal.mag2(); << 1247 // --->> Danger: Side effects on sub-navigator voxel information <<---
1628 if( std::fabs ( normMag2 - 1.0 ) < perTho << 1248 // Could be replaced again by 'granular' calls to sub-navigator
>> 1249 // locates (similar side-effects, but faster.
>> 1250 // Solutions:
>> 1251 // 1) Re-locate (to where?)
>> 1252 // 2) Insure that the methods using (G4ComputeStep?)
>> 1253 // does a relocation (if information is disturbed only ?)
>> 1254
>> 1255 #ifdef G4DEBUG_NAVIGATION
>> 1256 if( fVerbose >= 2 )
1629 { 1257 {
1630 *pNormalCalculated = true; // ComputeS << 1258 G4cout << " G4Navigator::ComputeSafety() relocates-in-volume to point: "
1631 // (fExitin << 1259 << pGlobalpoint << G4endl;
>> 1260 }
>> 1261 #endif
>> 1262 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
>> 1263 G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
>> 1264 G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
>> 1265 G4ThreeVector localPoint = ComputeLocalPoint(pGlobalpoint);
>> 1266
>> 1267 if ( fHistory.GetTopVolumeType()!=kReplica )
>> 1268 {
>> 1269 switch(CharacteriseDaughters(motherLogical))
>> 1270 {
>> 1271 case kNormal:
>> 1272 if ( pVoxelHeader )
>> 1273 {
>> 1274 newSafety=fvoxelNav.ComputeSafety(localPoint,fHistory,pMaxLength);
>> 1275 }
>> 1276 else
>> 1277 {
>> 1278 newSafety=fnormalNav.ComputeSafety(localPoint,fHistory,pMaxLength);
>> 1279 }
>> 1280 break;
>> 1281 case kParameterised:
>> 1282 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
>> 1283 {
>> 1284 newSafety = fparamNav.ComputeSafety(localPoint,fHistory,pMaxLength);
>> 1285 }
>> 1286 else // Regular structure
>> 1287 {
>> 1288 newSafety = fregularNav.ComputeSafety(localPoint,fHistory,pMaxLength);
>> 1289 }
>> 1290 break;
>> 1291 case kReplica:
>> 1292 G4Exception("G4Navigator::ComputeSafety()", "NotApplicable",
>> 1293 FatalException, "Not applicable for replicated volumes.");
>> 1294 break;
>> 1295 }
1632 } 1296 }
1633 else 1297 else
1634 { 1298 {
1635 G4ExceptionDescription message; << 1299 newSafety = freplicaNav.ComputeSafety(pGlobalpoint, localPoint,
1636 message.precision(10); << 1300 fHistory, pMaxLength);
1637 message << " WARNING> Expected normal- <<
1638 << " i.e. a unit vector!" << G <<
1639 << " - but |normal| = " << <<
1640 << " - and |normal|^2 = " << <<
1641 << " which differs from 1.0 by <<
1642 << " n = " << fExitNormalGlo <<
1643 << " Global point: " << Inters <<
1644 << " Volume: " << fHistory.Get <<
1645 #ifdef G4VERBOSE <<
1646 G4LogicalVolume* candLog = fHistory.Ge <<
1647 if ( candLog != nullptr ) <<
1648 { <<
1649 message << " Solid: " << candLog->Ge <<
1650 << ", Type: " << candLog->Ge <<
1651 << *candLog->GetSolid() << G <<
1652 } <<
1653 #endif <<
1654 message << "========================== <<
1655 << G4endl; <<
1656 G4int oldVerbose = fVerbose; <<
1657 fVerbose = 4; <<
1658 message << " State of Navigator: " < <<
1659 message << *this << G4endl; <<
1660 fVerbose = oldVerbose; <<
1661 message << "========================== <<
1662 << G4endl; <<
1663 <<
1664 G4Exception("G4Navigator::GetGlobalExi <<
1665 "GeomNav0003",JustWarning, <<
1666 "Value obtained from stored glo <<
1667 <<
1668 // (Re)Compute it now -- as either it <<
1669 // <<
1670 localNormal = GetLocalExitNormalAndChe <<
1671 <<
1672 *pNormalCalculated = fCalculatedExitNo <<
1673 globalNormal = fHistory.GetTopTransfor <<
1674 .InverseTransformAxis(lo <<
1675 } 1301 }
1676 } 1302 }
1677 else << 1303 else // if( endpointOnSurface && stayedOnEndpoint )
1678 { 1304 {
1679 localNormal = GetLocalExitNormalAndCheck( <<
1680 *pNormalCalculated = fCalculatedExitNorma <<
1681 <<
1682 #ifdef G4DEBUG_NAVIGATION 1305 #ifdef G4DEBUG_NAVIGATION
1683 usingStored = false; << 1306 if( fVerbose >= 2 )
1684 <<
1685 if( (!validNormal) && !fCalculatedExitNor <<
1686 { 1307 {
1687 G4ExceptionDescription edN; << 1308 G4cout << " G4Navigator::ComputeSafety() finds that point - "
1688 edN << " Calculated = " << fCalculated << 1309 << pGlobalpoint << " - is on surface " << G4endl;
1689 edN << " Entering= " << fEntering << << 1310 if( fEnteredDaughter ) { G4cout << " entered new daughter volume"; }
1690 G4int oldVerbose = this->GetVerboseLeve << 1311 if( fExitedMother ) { G4cout << " and exited previous volume."; }
1691 this->SetVerboseLevel(4); << 1312 G4cout << G4endl;
1692 edN << " State of Navigator: " << G4e << 1313 G4cout << " EndPoint was = " << fStepEndPoint << G4endl;
1693 edN << *this << G4endl; << 1314 }
1694 this->SetVerboseLevel( oldVerbose ); <<
1695 <<
1696 G4Exception("G4Navigator::GetGlobalExit <<
1697 "GeomNav0003", JustWarning, <<
1698 "LocalExitNormalAndCheck() <<
1699 } <<
1700 #endif <<
1701 <<
1702 G4double localMag2 = localNormal.mag2(); <<
1703 if( validNormal && (std::fabs(localMag2- <<
1704 { <<
1705 G4ExceptionDescription edN; <<
1706 edN.precision(10); <<
1707 edN << "G4Navigator::GetGlobalExitNorm <<
1708 << " Using Local Normal - from ca <<
1709 << G4endl <<
1710 << " Local Exit Normal : " << " <<
1711 << " vec = " << localNormal << G4e <<
1712 << " Global Exit Normal : " << " <<
1713 << " vec = " << globalNormal << G4 <<
1714 << " Global point: " << Intersect <<
1715 edN << " Calculated It = " << fC <<
1716 << " Volume: " << fHistory.GetTop <<
1717 #ifdef G4VERBOSE <<
1718 G4LogicalVolume* candLog = fHistory.Ge <<
1719 if ( candLog != nullptr ) <<
1720 { <<
1721 edN << " Solid: " << candLog->GetSo <<
1722 << ", Type: " << candLog->GetSol <<
1723 << *candLog->GetSolid(); <<
1724 } <<
1725 #endif 1315 #endif
1726 G4Exception("G4Navigator::GetGlobalExi << 1316 newSafety = 0.0;
1727 "GeomNav0003",JustWarning, <<
1728 "Value obtained from new l <<
1729 localNormal = localNormal.unit(); // S <<
1730 } <<
1731 globalNormal = fHistory.GetTopTransform( <<
1732 .InverseTransformAxis(loca <<
1733 } 1317 }
1734 1318
1735 #ifdef G4DEBUG_NAVIGATION << 1319 // Remember last safety origin & value
1736 if( usingStored ) <<
1737 { <<
1738 G4ThreeVector globalNormAgn; <<
1739 <<
1740 localNormal = GetLocalExitNormalAndCheck( <<
1741 <<
1742 globalNormAgn = fHistory.GetTopTransform( <<
1743 .InverseTransformAxis(loca <<
1744 <<
1745 // Check the value computed against fExit <<
1746 G4ThreeVector diffNorm = globalNormAgn - <<
1747 if( diffNorm.mag2() > kToleranceNormalChe <<
1748 { <<
1749 G4ExceptionDescription edDfn; <<
1750 edDfn << "Found difference in normals i <<
1751 << "- when Get is called after Co <<
1752 edDfn << " Magnitude of diff = " <<
1753 edDfn << " Normal stored (Global) <<
1754 << G4endl; <<
1755 edDfn << " Global Computed from Local <<
1756 G4Exception("G4Navigator::GetGlobalExit <<
1757 JustWarning, edDfn); <<
1758 } <<
1759 } <<
1760 #endif <<
1761 <<
1762 // Synchronise stored global exit normal as <<
1763 // 1320 //
1764 fExitNormalGlobalFrame = globalNormal; << 1321 fPreviousSftOrigin = pGlobalpoint;
1765 << 1322 fPreviousSafety = newSafety;
1766 return globalNormal; <<
1767 } <<
1768 <<
1769 // ****************************************** <<
1770 // ComputeSafety <<
1771 // <<
1772 // It assumes that it will be <<
1773 // i) called at the Point in the same volume <<
1774 // ComputeStep. <<
1775 // ii) after (or at the end of) ComputeStep O <<
1776 // ****************************************** <<
1777 // <<
1778 G4double G4Navigator::ComputeSafety( const G4 <<
1779 const G4 <<
1780 const G4 <<
1781 { <<
1782 G4VPhysicalVolume *motherPhysical = fHisto <<
1783 G4double safety = 0.0; <<
1784 <<
1785 G4double distEndpointSq = (pGlobalpoint-fSt <<
1786 G4bool stayedOnEndpoint = distEndpointSq < <<
1787 G4bool endpointOnSurface = fEnteredDaughter <<
1788 1323
1789 G4bool onSurface = endpointOnSurface && sta << 1324 #ifdef G4DEBUG_NAVIGATION
1790 if( ! onSurface ) << 1325 if( fVerbose > 1 )
1791 { 1326 {
1792 safety= fpSafetyCalculator->SafetyInCurre << 1327 G4cout << " ---- Exiting ComputeSafety " << G4endl;
1793 // offload to G4SafetyCalculator - avoids << 1328 if( fVerbose > 2 ) { PrintState(); }
1794 << 1329 G4cout << " Returned value of Safety = " << newSafety << G4endl;
1795 // Remember last safety origin & value <<
1796 // <<
1797 fPreviousSftOrigin = pGlobalpoint; <<
1798 fPreviousSafety = safety; <<
1799 // We overwrite the Safety 'sphere' - kee <<
1800 } 1330 }
>> 1331 G4cout.precision(oldcoutPrec);
>> 1332 #endif
1801 1333
1802 return safety; << 1334 return newSafety;
1803 } 1335 }
1804 1336
1805 // ****************************************** 1337 // ********************************************************************
1806 // CreateTouchableHistoryHandle 1338 // CreateTouchableHistoryHandle
1807 // ****************************************** 1339 // ********************************************************************
1808 // 1340 //
1809 G4TouchableHandle G4Navigator::CreateTouchabl << 1341 G4TouchableHistoryHandle G4Navigator::CreateTouchableHistoryHandle() const
1810 { 1342 {
1811 return G4TouchableHandle( CreateTouchableHi << 1343 return G4TouchableHistoryHandle( CreateTouchableHistory() );
1812 } 1344 }
1813 1345
1814 // ****************************************** 1346 // ********************************************************************
1815 // PrintState 1347 // PrintState
1816 // ****************************************** 1348 // ********************************************************************
1817 // 1349 //
1818 void G4Navigator::PrintState() const 1350 void G4Navigator::PrintState() const
1819 { 1351 {
1820 G4long oldcoutPrec = G4cout.precision(4); << 1352 G4int oldcoutPrec = G4cout.precision(4);
1821 if( fVerbose >= 4 ) << 1353 if( fVerbose == 4 )
1822 { 1354 {
1823 G4cout << "The current state of G4Navigat 1355 G4cout << "The current state of G4Navigator is: " << G4endl;
1824 G4cout << " ValidExitNormal= " << fValid << 1356 G4cout << " ValidExitNormal= " << fValidExitNormal << G4endl
1825 << " ExitNormal = " << fExitN << 1357 << " ExitNormal = " << fExitNormal << G4endl
1826 << " Exiting = " << fExiti << 1358 << " Exiting = " << fExiting << G4endl
1827 << " Entering = " << fEnter << 1359 << " Entering = " << fEntering << G4endl
1828 << " BlockedPhysicalVolume= " ; 1360 << " BlockedPhysicalVolume= " ;
1829 if (fBlockedPhysicalVolume==nullptr) << 1361 if (fBlockedPhysicalVolume==0)
1830 { <<
1831 G4cout << "None"; 1362 G4cout << "None";
1832 } <<
1833 else 1363 else
1834 { <<
1835 G4cout << fBlockedPhysicalVolume->GetNa 1364 G4cout << fBlockedPhysicalVolume->GetName();
1836 } <<
1837 G4cout << G4endl 1365 G4cout << G4endl
1838 << " BlockedReplicaNo = " << << 1366 << " BlockedReplicaNo = " << fBlockedReplicaNo << G4endl
1839 << " LastStepWasZero = " << << 1367 << " LastStepWasZero = " << fLastStepWasZero << G4endl
1840 << G4endl; 1368 << G4endl;
1841 } 1369 }
1842 if( ( 1 < fVerbose) && (fVerbose < 4) ) 1370 if( ( 1 < fVerbose) && (fVerbose < 4) )
1843 { 1371 {
1844 G4cout << G4endl; // Make sure to line up << 1372 G4cout << std::setw(30) << " ExitNormal " << " "
1845 G4cout << std::setw(30) << " ExitNormal " <<
1846 << std::setw( 5) << " Valid " 1373 << std::setw( 5) << " Valid " << " "
1847 << std::setw( 9) << " Exiting " 1374 << std::setw( 9) << " Exiting " << " "
1848 << std::setw( 9) << " Entering" 1375 << std::setw( 9) << " Entering" << " "
1849 << std::setw(15) << " Blocked:Volu 1376 << std::setw(15) << " Blocked:Volume " << " "
1850 << std::setw( 9) << " ReplicaNo" 1377 << std::setw( 9) << " ReplicaNo" << " "
1851 << std::setw( 8) << " LastStepZero 1378 << std::setw( 8) << " LastStepZero " << " "
1852 << G4endl; 1379 << G4endl;
1853 G4cout << "( " << std::setw(7) << fExitNo 1380 G4cout << "( " << std::setw(7) << fExitNormal.x()
1854 << ", " << std::setw(7) << fExitNo 1381 << ", " << std::setw(7) << fExitNormal.y()
1855 << ", " << std::setw(7) << fExitNo 1382 << ", " << std::setw(7) << fExitNormal.z() << " ) "
1856 << std::setw( 5) << fValidExitNor 1383 << std::setw( 5) << fValidExitNormal << " "
1857 << std::setw( 9) << fExiting 1384 << std::setw( 9) << fExiting << " "
1858 << std::setw( 9) << fEntering 1385 << std::setw( 9) << fEntering << " ";
1859 if ( fBlockedPhysicalVolume == nullptr ) << 1386 if ( fBlockedPhysicalVolume==0 )
1860 { G4cout << std::setw(15) << "None"; } << 1387 G4cout << std::setw(15) << "None";
1861 else 1388 else
1862 { G4cout << std::setw(15)<< fBlockedPhysi << 1389 G4cout << std::setw(15)<< fBlockedPhysicalVolume->GetName();
1863 G4cout << std::setw( 9) << fBlockedRepli << 1390 G4cout << std::setw( 9) << fBlockedReplicaNo << " "
1864 << std::setw( 8) << fLastStepWasZ << 1391 << std::setw( 8) << fLastStepWasZero << " "
1865 << G4endl; << 1392 << G4endl;
1866 } 1393 }
1867 if( fVerbose > 2 ) 1394 if( fVerbose > 2 )
1868 { 1395 {
1869 G4cout.precision(8); 1396 G4cout.precision(8);
1870 G4cout << " Current Localpoint = " << fLa 1397 G4cout << " Current Localpoint = " << fLastLocatedPointLocal << G4endl;
1871 G4cout << " PreviousSftOrigin = " << fPr 1398 G4cout << " PreviousSftOrigin = " << fPreviousSftOrigin << G4endl;
1872 G4cout << " PreviousSafety = " << fPr 1399 G4cout << " PreviousSafety = " << fPreviousSafety << G4endl;
1873 } 1400 }
1874 G4cout.precision(oldcoutPrec); 1401 G4cout.precision(oldcoutPrec);
1875 } 1402 }
1876 1403
1877 // ****************************************** 1404 // ********************************************************************
1878 // ComputeStepLog <<
1879 // ****************************************** <<
1880 // <<
1881 void G4Navigator::ComputeStepLog(const G4Thre <<
1882 G4doub <<
1883 { <<
1884 // The following checks only make sense if <<
1885 // than the tolerance. <<
1886 <<
1887 const G4double fAccuracyForWarning = kCar <<
1888 fAccuracyForException = 1000 <<
1889 <<
1890 G4ThreeVector OriginalGlobalpoint = fHistor <<
1891 InverseTransfo <<
1892 <<
1893 G4double shiftOriginSafSq = (fPreviousSftOr <<
1894 <<
1895 // Check that the starting point of this st <<
1896 // within the isotropic safety sphere of th <<
1897 // to a accuracy/precision given by fAccur <<
1898 // If so give warning. <<
1899 // If it fails by more than fAccuracyForE <<
1900 // <<
1901 if( shiftOriginSafSq >= sqr(fPreviousSafety <<
1902 { <<
1903 G4double shiftOrigin = std::sqrt(shiftOri <<
1904 G4double diffShiftSaf = shiftOrigin - fPr <<
1905 <<
1906 if( diffShiftSaf > fAccuracyForWarning ) <<
1907 { <<
1908 G4long oldcoutPrec = G4cout.precision(8 <<
1909 G4long oldcerrPrec = G4cerr.precision(1 <<
1910 std::ostringstream message, suggestion; <<
1911 message << "Accuracy error or slightly <<
1912 << G4endl <<
1913 << " The Step's starting po <<
1914 << std::sqrt(moveLenSq)/mm << " <<
1915 << " since the last call to <<
1916 << " This has resulted in m <<
1917 << shiftOrigin/mm << " mm " <<
1918 << " from the last point at whi <<
1919 << " was calculated " << G4 <<
1920 << " which is more than the <<
1921 << fPreviousSafety/mm << " mm <<
1922 << " This difference is " <<
1923 << diffShiftSaf/mm << " mm." << <<
1924 << " The tolerated accuracy <<
1925 << fAccuracyForException/mm << <<
1926 <<
1927 suggestion << " "; <<
1928 static G4ThreadLocal G4int warnNow = 0; <<
1929 if( ((++warnNow % 100) == 1) ) <<
1930 { <<
1931 message << G4endl <<
1932 << " This problem can be due <<
1933 << " - a process that has p <<
1934 << " larger than the current s <<
1935 << " - inaccuracy in the co <<
1936 suggestion << "We suggest that you " <<
1937 << " - find i) what part <<
1938 << " ii) through what part <<
1939 << " for example by r <<
1940 << G4endl <<
1941 << " /tracking/ver <<
1942 << " - check which proc <<
1943 << " this particle (and lo <<
1944 << G4endl <<
1945 << " - in case, create a <<
1946 << " of this event using:" <<
1947 << " /tracking/ver <<
1948 } <<
1949 G4Exception("G4Navigator::ComputeStep() <<
1950 "GeomNav1002", JustWarning, <<
1951 message, G4String(suggestio <<
1952 G4cout.precision(oldcoutPrec); <<
1953 G4cerr.precision(oldcerrPrec); <<
1954 } <<
1955 #ifdef G4DEBUG_NAVIGATION <<
1956 else <<
1957 { <<
1958 G4cerr << "WARNING - G4Navigator::Compu <<
1959 << " The Step's startin <<
1960 << std::sqrt(moveLenSq) << "," < <<
1961 << " which has taken it <<
1962 << " the current safety. " << G4 <<
1963 } <<
1964 #endif <<
1965 } <<
1966 G4double safetyPlus = fPreviousSafety + fAc <<
1967 if ( shiftOriginSafSq > sqr(safetyPlus) ) <<
1968 { <<
1969 std::ostringstream message; <<
1970 message << "May lead to a crash or unreli <<
1971 << " Position has shifted <<
1972 << " notifying the navigator !" < <<
1973 << " Tolerated safety: " < <<
1974 << " Computed shift : " < <<
1975 G4Exception("G4Navigator::ComputeStep()", <<
1976 JustWarning, message); <<
1977 } <<
1978 } <<
1979 <<
1980 // ****************************************** <<
1981 // CheckOverlapsIterative <<
1982 // ****************************************** <<
1983 // <<
1984 G4bool G4Navigator::CheckOverlapsIterative(G4 <<
1985 { <<
1986 // Check and report overlaps <<
1987 // <<
1988 G4bool foundOverlap = false; <<
1989 G4int nPoints = 300000, ntrials = 9, numO <<
1990 G4double trialLength = 1.0 * CLHEP::centim <<
1991 while ( ntrials-- > 0 && !foundOverlap ) <<
1992 { <<
1993 if ( fVerbose > 1 ) <<
1994 { <<
1995 G4cout << " ** Running overlap checks <<
1996 << vol->GetName() <<
1997 << " with length = " << trialLe <<
1998 } <<
1999 foundOverlap = vol->CheckOverlaps(nPoints <<
2000 fVerbos <<
2001 trialLength *= 0.1; <<
2002 if ( trialLength <= 1.0e-5 ) { numOverlap <<
2003 } <<
2004 return foundOverlap; <<
2005 } <<
2006 <<
2007 // ****************************************** <<
2008 // Operator << 1405 // Operator <<
2009 // ****************************************** 1406 // ********************************************************************
2010 // 1407 //
2011 std::ostream& operator << (std::ostream &os,c 1408 std::ostream& operator << (std::ostream &os,const G4Navigator &n)
2012 { 1409 {
2013 // Old version did only the following: << 1410 os << "Current History: " << G4endl << n.fHistory;
2014 // os << "Current History: " << G4endl << n <<
2015 // Old behaviour is recovered for fVerbose <<
2016 <<
2017 // Adapted from G4Navigator::PrintState() c <<
2018 <<
2019 G4long oldcoutPrec = os.precision(4); <<
2020 if( n.fVerbose >= 4 ) <<
2021 { <<
2022 os << "The current state of G4Navigator i <<
2023 os << " ValidExitNormal= " << n.fValidEx <<
2024 << " ExitNormal = " << n.fExitNormal <<
2025 << " Exiting = " << n.fExiting <<
2026 << " Entering = " << n.fEntering <<
2027 << " BlockedPhysicalVolume= " ; <<
2028 if (n.fBlockedPhysicalVolume==nullptr) <<
2029 { <<
2030 os << "None"; <<
2031 } <<
2032 else <<
2033 { <<
2034 os << n.fBlockedPhysicalVolume->GetName <<
2035 } <<
2036 os << G4endl <<
2037 << " BlockedReplicaNo = " << n.fBlo <<
2038 << " LastStepWasZero = " << n.fLa <<
2039 << G4endl; <<
2040 } <<
2041 if( ( 1 < n.fVerbose) && (n.fVerbose < 4) ) <<
2042 { <<
2043 os << G4endl; // Make sure to line up <<
2044 os << std::setw(30) << " ExitNormal " << <<
2045 << std::setw( 5) << " Valid " << " <<
2046 << std::setw( 9) << " Exiting " << " <<
2047 << std::setw( 9) << " Entering" << " <<
2048 << std::setw(15) << " Blocked:Volume " < <<
2049 << std::setw( 9) << " ReplicaNo" < <<
2050 << std::setw( 8) << " LastStepZero " < <<
2051 << G4endl; <<
2052 os << "( " << std::setw(7) << n.fExitNorm <<
2053 << ", " << std::setw(7) << n.fExitNormal. <<
2054 << ", " << std::setw(7) << n.fExitNormal. <<
2055 << std::setw( 5) << n.fValidExitNormal <<
2056 << std::setw( 9) << n.fExiting <<
2057 << std::setw( 9) << n.fEntering <<
2058 if ( n.fBlockedPhysicalVolume==nullptr ) <<
2059 { os << std::setw(15) << "None"; } <<
2060 else <<
2061 { os << std::setw(15)<< n.fBlockedPhysi <<
2062 os << std::setw( 9) << n.fBlockedReplica <<
2063 << std::setw( 8) << n.fLastStepWasZero <<
2064 << G4endl; <<
2065 } <<
2066 if( n.fVerbose > 2 ) <<
2067 { <<
2068 os.precision(8); <<
2069 os << " Current Localpoint = " << n.fLast <<
2070 os << " PreviousSftOrigin = " << n.fPrev <<
2071 os << " PreviousSafety = " << n.fPrev <<
2072 } <<
2073 if( n.fVerbose > 3 || n.fVerbose == 0 ) <<
2074 { <<
2075 os << "Current History: " << G4endl << n. <<
2076 } <<
2077 <<
2078 os.precision(oldcoutPrec); <<
2079 return os; 1411 return os;
2080 } <<
2081 <<
2082 // ****************************************** <<
2083 // SetVoxelNavigation: alternative navigator <<
2084 // ****************************************** <<
2085 // <<
2086 void G4Navigator::SetVoxelNavigation(G4VoxelN <<
2087 { <<
2088 delete fpvoxelNav; <<
2089 fpvoxelNav = voxelNav; <<
2090 } <<
2091 <<
2092 // ****************************************** <<
2093 // InformLastStep: derived navigators can inf <<
2094 // used to update fLastStepWa <<
2095 // ****************************************** <<
2096 void G4Navigator::InformLastStep(G4double la <<
2097 G4bool exit <<
2098 { <<
2099 G4bool zeroStep = ( lastStep == 0.0 ); <<
2100 fLocatedOnEdge = fLastStepWasZero && zero <<
2101 fLastStepWasZero = zeroStep; <<
2102 <<
2103 fExiting = exitsMotherVol; <<
2104 fEntering = entersDaughtVol; <<
2105 } <<
2106 <<
2107 // ****************************************** <<
2108 // SetExternalNavigation <<
2109 // ****************************************** <<
2110 // <<
2111 void G4Navigator::SetExternalNavigation(G4VEx <<
2112 { <<
2113 fpExternalNav = externalNav; <<
2114 fpSafetyCalculator->SetExternalNavigation(e <<
2115 } 1412 }
2116 1413