| 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.46 2010-11-15 14:03:27 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), fWarnPush(true)
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
155 if( fVerbose > 2 ) 130 if( fVerbose > 2 )
156 { 131 {
157 G4long oldcoutPrec = G4cout.precision(8); << 132 G4int 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); 142 G4cout.precision(oldcoutPrec);
168 } 143 }
169 #endif 144 #endif
170 145
171 G4int noLevelsExited = 0; <<
172 <<
173 if ( !relativeSearch ) 146 if ( !relativeSearch )
174 { 147 {
175 ResetStackAndState(); 148 ResetStackAndState();
176 } 149 }
177 else 150 else
178 { 151 {
179 if ( fWasLimitedByGeometry ) 152 if ( fWasLimitedByGeometry )
180 { 153 {
181 fWasLimitedByGeometry = false; 154 fWasLimitedByGeometry = false;
182 fEnteredDaughter = fEntering; // Remem 155 fEnteredDaughter = fEntering; // Remember
183 fExitedMother = fExiting; // Remem 156 fExitedMother = fExiting; // Remember
184 if ( fExiting ) 157 if ( fExiting )
185 { 158 {
186 ++noLevelsExited; // count this first << 159 if ( fHistory.GetDepth() )
187 <<
188 if ( fHistory.GetDepth() != 0 ) <<
189 { 160 {
190 fBlockedPhysicalVolume = fHistory.Ge 161 fBlockedPhysicalVolume = fHistory.GetTopVolume();
191 fBlockedReplicaNo = fHistory.GetTopR 162 fBlockedReplicaNo = fHistory.GetTopReplicaNo();
192 fHistory.BackLevel(); 163 fHistory.BackLevel();
193 } 164 }
194 else 165 else
195 { 166 {
196 fLastLocatedPointLocal = localPoint; 167 fLastLocatedPointLocal = localPoint;
197 fLocatedOutsideWorld = true; 168 fLocatedOutsideWorld = true;
198 fBlockedPhysicalVolume = nullptr; << 169 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 } 170 }
206 // A fix for the case where a volume i 171 // A fix for the case where a volume is "entered" at an edge
207 // and a coincident surface exists out 172 // and a coincident surface exists outside it.
208 // - This stops it from exiting furth 173 // - This stops it from exiting further volumes and cycling
209 // - However ReplicaNavigator treats 174 // - However ReplicaNavigator treats this case itself
210 // 175 //
211 // assert( fBlockedPhysicalVolume!=0 ) <<
212 <<
213 // Expect to be on edge => on surface <<
214 // <<
215 if ( fLocatedOnEdge && (VolumeType(fBl 176 if ( fLocatedOnEdge && (VolumeType(fBlockedPhysicalVolume)!=kReplica ))
216 { 177 {
217 fExiting = false; << 178 fExiting= false;
218 // Consider effect on Exit Normal !? <<
219 } 179 }
220 } 180 }
221 else 181 else
222 if ( fEntering ) 182 if ( fEntering )
223 { 183 {
224 switch (VolumeType(fBlockedPhysicalV 184 switch (VolumeType(fBlockedPhysicalVolume))
225 { 185 {
226 case kNormal: 186 case kNormal:
227 fHistory.NewLevel(fBlockedPhysic 187 fHistory.NewLevel(fBlockedPhysicalVolume, kNormal,
228 fBlockedPhysic 188 fBlockedPhysicalVolume->GetCopyNo());
229 break; 189 break;
230 case kReplica: 190 case kReplica:
231 freplicaNav.ComputeTransformatio 191 freplicaNav.ComputeTransformation(fBlockedReplicaNo,
232 192 fBlockedPhysicalVolume);
233 fHistory.NewLevel(fBlockedPhysic 193 fHistory.NewLevel(fBlockedPhysicalVolume, kReplica,
234 fBlockedReplic 194 fBlockedReplicaNo);
235 fBlockedPhysicalVolume->SetCopyN 195 fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
236 break; 196 break;
237 case kParameterised: 197 case kParameterised:
238 if( fBlockedPhysicalVolume->GetR 198 if( fBlockedPhysicalVolume->GetRegularStructureId() == 0 )
239 { 199 {
240 G4VSolid *pSolid; 200 G4VSolid *pSolid;
241 G4VPVParameterisation *pParam; 201 G4VPVParameterisation *pParam;
242 G4TouchableHistory parentTouch 202 G4TouchableHistory parentTouchable( fHistory );
243 pParam = fBlockedPhysicalVolum 203 pParam = fBlockedPhysicalVolume->GetParameterisation();
244 pSolid = pParam->ComputeSolid( 204 pSolid = pParam->ComputeSolid(fBlockedReplicaNo,
245 205 fBlockedPhysicalVolume);
246 pSolid->ComputeDimensions(pPar 206 pSolid->ComputeDimensions(pParam, fBlockedReplicaNo,
247 fBlo 207 fBlockedPhysicalVolume);
248 pParam->ComputeTransformation( 208 pParam->ComputeTransformation(fBlockedReplicaNo,
249 209 fBlockedPhysicalVolume);
250 fHistory.NewLevel(fBlockedPhys 210 fHistory.NewLevel(fBlockedPhysicalVolume, kParameterised,
251 fBlockedRepl 211 fBlockedReplicaNo);
252 fBlockedPhysicalVolume->SetCop 212 fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo);
253 // 213 //
254 // Set the correct solid and m 214 // Set the correct solid and material in Logical Volume
255 // 215 //
256 G4LogicalVolume *pLogical; 216 G4LogicalVolume *pLogical;
257 pLogical = fBlockedPhysicalVol 217 pLogical = fBlockedPhysicalVolume->GetLogicalVolume();
258 pLogical->SetSolid( pSolid ); 218 pLogical->SetSolid( pSolid );
259 pLogical->UpdateMaterial(pPara 219 pLogical->UpdateMaterial(pParam ->
260 ComputeMaterial(fBlockedRepl 220 ComputeMaterial(fBlockedReplicaNo,
261 fBlockedPhys 221 fBlockedPhysicalVolume,
262 &parentTouch 222 &parentTouchable));
263 } 223 }
264 break; 224 break;
265 case kExternal: <<
266 G4Exception("G4Navigator::Locate <<
267 "GeomNav0001", Fatal <<
268 "Extra levels not ap <<
269 break; <<
270 } 225 }
271 fEntering = false; 226 fEntering = false;
272 fBlockedPhysicalVolume = nullptr; << 227 fBlockedPhysicalVolume = 0;
273 localPoint = fHistory.GetTopTransfor 228 localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
274 notKnownContained = false; 229 notKnownContained = false;
275 } 230 }
276 } 231 }
277 else 232 else
278 { 233 {
279 fBlockedPhysicalVolume = nullptr; << 234 fBlockedPhysicalVolume = 0;
280 fEntering = false; 235 fEntering = false;
281 fEnteredDaughter = false; // Full Step 236 fEnteredDaughter = false; // Full Step was not taken, did not enter
282 fExiting = false; 237 fExiting = false;
283 fExitedMother = false; // Full Step 238 fExitedMother = false; // Full Step was not taken, did not exit
284 } 239 }
285 } 240 }
286 // 241 //
287 // Search from top of history up through geo 242 // Search from top of history up through geometry until
288 // containing volume found: 243 // containing volume found:
289 // If on 244 // If on
290 // o OUTSIDE - Back up level, not/no longer 245 // o OUTSIDE - Back up level, not/no longer exiting volumes
291 // o SURFACE and EXITING - Back up level, se 246 // o SURFACE and EXITING - Back up level, setting new blocking no.s
292 // else 247 // else
293 // o containing volume found 248 // o containing volume found
294 // 249 //
295 << 250 while (notKnownContained)
296 while (notKnownContained) // Loop checking, <<
297 { 251 {
298 EVolume topVolumeType = fHistory.GetTopVol << 252 if ( fHistory.GetTopVolumeType()!=kReplica )
299 if (topVolumeType!=kReplica && topVolumeTy <<
300 { 253 {
301 targetSolid = fHistory.GetTopVolume()->G 254 targetSolid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid();
302 localPoint = fHistory.GetTopTransform(). 255 localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint);
303 insideCode = targetSolid->Inside(localPo 256 insideCode = targetSolid->Inside(localPoint);
304 #ifdef G4VERBOSE 257 #ifdef G4VERBOSE
305 if(( fVerbose == 1 ) && ( fCheck )) 258 if(( fVerbose == 1 ) && ( fCheck ))
306 { 259 {
307 G4String solidResponse = "-kInside-"; << 260 G4String solidResponse = "-kInside-";
308 if (insideCode == kOutside) << 261 if (insideCode == kOutside)
309 { << 262 solidResponse = "-kOutside-";
310 solidResponse = "-kOutside-"; << 263 else if (insideCode == kSurface)
311 } << 264 solidResponse = "-kSurface-";
312 else if (insideCode == kSurface) << 265 G4cout << "*** G4Navigator::LocateGlobalPointAndSetup(): ***" << G4endl
313 { << 266 << " Invoked Inside() for solid: " << targetSolid->GetName()
314 solidResponse = "-kSurface-"; << 267 << ". Solid replied: " << solidResponse << G4endl
315 } << 268 << " 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 } 269 }
321 #endif 270 #endif
322 } 271 }
323 else 272 else
324 { 273 {
325 if( topVolumeType == kReplica ) << 274 insideCode = freplicaNav.BackLocate(fHistory, globalPoint, localPoint,
326 { << 275 fExiting, notKnownContained);
327 insideCode = freplicaNav.BackLocate( << 276 // !CARE! if notKnownContained returns false then the point is within
328 << 277 // the containing placement volume of the replica(s). If insidecode
329 // !CARE! if notKnownContained retur << 278 // will result in the history being backed up one level, then the
330 // the containing placement volume o << 279 // 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 } 280 }
343 << 281 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 { 282 {
350 ++noLevelsExited; << 283 if ( fHistory.GetDepth() )
351 <<
352 // Exiting world volume <<
353 if ( fHistory.GetDepth() == 0 ) <<
354 { 284 {
355 fLocatedOutsideWorld = true; << 285 fBlockedPhysicalVolume = fHistory.GetTopVolume();
356 fLastLocatedPointLocal = localPoint; << 286 fBlockedReplicaNo = fHistory.GetTopReplicaNo();
357 return nullptr; << 287 fHistory.BackLevel();
>> 288 fExiting = false;
358 } 289 }
359 << 290 else
360 fBlockedPhysicalVolume = fHistory.GetTop <<
361 fBlockedReplicaNo = fHistory.GetTopRepli <<
362 fHistory.BackLevel(); <<
363 fExiting = false; <<
364 <<
365 if( noLevelsExited > 1 ) <<
366 { 291 {
367 // The first transformation was done b << 292 fLastLocatedPointLocal = localPoint;
368 // << 293 fLocatedOutsideWorld = true;
369 if(const auto *mRot = fBlockedPhysical << 294 return 0; // Have exited world volume
370 { <<
371 fGrandMotherExitNormal *= (*mRot).in <<
372 fChangedGrandMotherRefFrame = true; <<
373 } <<
374 } 295 }
375 continue; <<
376 } 296 }
377 << 297 else
378 // Point is on the surface of a volume << 298 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 { 299 {
395 G4ThreeVector normal = targetSolid->Su << 300 G4bool isExiting = fExiting;
396 directionExiting = normal.dot(localDir << 301 if( (!fExiting)&&considerDirection )
397 isExiting = isExiting || directionExit << 302 {
>> 303 // Figure out whether we are exiting this level's volume
>> 304 // by using the direction
>> 305 //
>> 306 G4bool directionExiting = false;
>> 307 G4ThreeVector localDirection =
>> 308 fHistory.GetTopTransform().TransformAxis(globalDirection);
>> 309 if ( fHistory.GetTopVolumeType()!=kReplica )
>> 310 {
>> 311 G4ThreeVector normal = targetSolid->SurfaceNormal(localPoint);
>> 312 directionExiting = normal.dot(localDirection) > 0.0;
>> 313 isExiting = isExiting || directionExiting;
>> 314 }
>> 315 }
>> 316 if( isExiting )
>> 317 {
>> 318 if ( fHistory.GetDepth() )
>> 319 {
>> 320 fBlockedPhysicalVolume = fHistory.GetTopVolume();
>> 321 fBlockedReplicaNo = fHistory.GetTopReplicaNo();
>> 322 fHistory.BackLevel();
>> 323 //
>> 324 // Still on surface but exited volume not necessarily convex
>> 325 //
>> 326 fValidExitNormal = false;
>> 327 }
>> 328 else
>> 329 {
>> 330 fLastLocatedPointLocal = localPoint;
>> 331 fLocatedOutsideWorld = true;
>> 332 return 0; // Have exited world volume
>> 333 }
>> 334 }
>> 335 else
>> 336 {
>> 337 notKnownContained=false;
>> 338 }
398 } 339 }
399 } << 340 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 { 341 {
428 fGrandMotherExitNormal *= (*mRot).inve << 342 notKnownContained=false;
429 fChangedGrandMotherRefFrame = true; <<
430 } 343 }
431 } <<
432 } // END while (notKnownContained) 344 } // END while (notKnownContained)
433 // 345 //
434 // Search downwards until deepest containing 346 // Search downwards until deepest containing volume found,
435 // blocking fBlockedPhysicalVolume/BlockedRe 347 // blocking fBlockedPhysicalVolume/BlockedReplicaNum
436 // 348 //
437 // 3 Cases: 349 // 3 Cases:
438 // 350 //
439 // o Parameterised daughters 351 // o Parameterised daughters
440 // =>Must be one G4PVParameterised daughte 352 // =>Must be one G4PVParameterised daughter & voxels
441 // o Positioned daughters & voxels 353 // o Positioned daughters & voxels
442 // o Positioned daughters & no voxels 354 // o Positioned daughters & no voxels
443 355
444 noResult = true; // noResult should be rena << 356 noResult = true; // noResult should be renamed to
445 // something like enteredL 357 // something like enteredLevel, as that is its meaning.
446 do 358 do
447 { 359 {
448 // Determine `type' of current mother volu 360 // Determine `type' of current mother volume
449 // 361 //
450 targetPhysical = fHistory.GetTopVolume(); 362 targetPhysical = fHistory.GetTopVolume();
451 if (targetPhysical == nullptr) { break; } << 363 if (!targetPhysical) { break; }
452 targetLogical = targetPhysical->GetLogical 364 targetLogical = targetPhysical->GetLogicalVolume();
453 switch( CharacteriseDaughters(targetLogica 365 switch( CharacteriseDaughters(targetLogical) )
454 { 366 {
455 case kNormal: 367 case kNormal:
456 if ( targetLogical->GetVoxelHeader() ! << 368 if ( targetLogical->GetVoxelHeader() ) // use optimised navigation
457 { 369 {
458 noResult = GetVoxelNavigator().Level << 370 noResult = fvoxelNav.LevelLocate(fHistory,
459 fBl 371 fBlockedPhysicalVolume,
460 fBl 372 fBlockedReplicaNo,
461 glo 373 globalPoint,
462 pGl 374 pGlobalDirection,
463 con 375 considerDirection,
464 loc 376 localPoint);
465 } 377 }
466 else // do not u 378 else // do not use optimised navigation
467 { 379 {
468 noResult = fnormalNav.LevelLocate(fH 380 noResult = fnormalNav.LevelLocate(fHistory,
469 fB 381 fBlockedPhysicalVolume,
470 fB 382 fBlockedReplicaNo,
471 gl 383 globalPoint,
472 pG 384 pGlobalDirection,
473 co 385 considerDirection,
474 lo 386 localPoint);
475 } 387 }
476 break; 388 break;
477 case kReplica: 389 case kReplica:
478 noResult = freplicaNav.LevelLocate(fHi 390 noResult = freplicaNav.LevelLocate(fHistory,
479 fBl 391 fBlockedPhysicalVolume,
480 fBl 392 fBlockedReplicaNo,
481 glo 393 globalPoint,
482 pGl 394 pGlobalDirection,
483 con 395 considerDirection,
484 loc 396 localPoint);
485 break; 397 break;
486 case kParameterised: 398 case kParameterised:
487 if( GetDaughtersRegularStructureId(tar 399 if( GetDaughtersRegularStructureId(targetLogical) != 1 )
488 { 400 {
489 noResult = fparamNav.LevelLocate(fHi 401 noResult = fparamNav.LevelLocate(fHistory,
490 fBl 402 fBlockedPhysicalVolume,
491 fBl 403 fBlockedReplicaNo,
492 glo 404 globalPoint,
493 pGl 405 pGlobalDirection,
494 con 406 considerDirection,
495 loc 407 localPoint);
496 } 408 }
497 else // Regular structure 409 else // Regular structure
498 { 410 {
499 noResult = fregularNav.LevelLocate(f 411 noResult = fregularNav.LevelLocate(fHistory,
500 f 412 fBlockedPhysicalVolume,
501 f 413 fBlockedReplicaNo,
502 g 414 globalPoint,
503 p 415 pGlobalDirection,
504 c 416 considerDirection,
505 l 417 localPoint);
506 } 418 }
507 break; 419 break;
508 case kExternal: <<
509 noResult = fpExternalNav->LevelLocate( <<
510 <<
511 <<
512 <<
513 <<
514 <<
515 <<
516 break; <<
517 } 420 }
518 421
519 // LevelLocate returns true if it finds a 422 // LevelLocate returns true if it finds a daughter volume
520 // in which globalPoint is inside (or on t 423 // in which globalPoint is inside (or on the surface).
521 424
522 if ( noResult ) 425 if ( noResult )
523 { 426 {
524 // Entering a daughter after ascending 427 // Entering a daughter after ascending
525 // 428 //
526 // The blocked volume is no longer valid 429 // The blocked volume is no longer valid - it was for another level
527 // 430 //
528 fBlockedPhysicalVolume = nullptr; << 431 fBlockedPhysicalVolume = 0;
529 fBlockedReplicaNo = -1; 432 fBlockedReplicaNo = -1;
530 433
531 // fEntering should be false -- else blo 434 // fEntering should be false -- else blockedVolume is assumed good.
532 // fEnteredDaughter is used for ExitNorm 435 // fEnteredDaughter is used for ExitNormal
533 // 436 //
534 fEntering = false; 437 fEntering = false;
535 fEnteredDaughter = true; 438 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 439 #ifdef G4DEBUG_NAVIGATION
552 if( fVerbose > 2 ) 440 if( fVerbose > 2 )
553 { 441 {
554 G4VPhysicalVolume* enteredPhysical = 442 G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume();
555 G4cout << "*** G4Navigator::LocateGlo << 443 G4cout << "*** G4Navigator::LocateGlobalPointAndSetup() ***" << G4endl;
556 G4cout << " Entering volume: " << 444 G4cout << " Entering volume: " << enteredPhysical->GetName()
557 << G4endl; 445 << G4endl;
558 } 446 }
559 #endif 447 #endif
560 } 448 }
561 } while (noResult); // Loop checking, 07.10 << 449 } while (noResult);
562 450
563 fLastLocatedPointLocal = localPoint; 451 fLastLocatedPointLocal = localPoint;
564 452
565 #ifdef G4VERBOSE 453 #ifdef G4VERBOSE
566 if( fVerbose >= 4 ) << 454 if( fVerbose == 4 )
567 { 455 {
568 G4long oldcoutPrec = G4cout.precision(8); << 456 G4int oldcoutPrec = G4cout.precision(8);
569 G4String curPhysVol_Name("None"); 457 G4String curPhysVol_Name("None");
570 if (targetPhysical != nullptr) { curPhysV << 458 if (targetPhysical) { curPhysVol_Name = targetPhysical->GetName(); }
571 G4cout << " Return value = new volume = 459 G4cout << " Return value = new volume = " << curPhysVol_Name << G4endl;
572 G4cout << " ----- Upon exiting:" << G4e 460 G4cout << " ----- Upon exiting:" << G4endl;
573 PrintState(); 461 PrintState();
574 if( fVerbose >= 5 ) << 462 #ifdef G4DEBUG_NAVIGATION
575 { << 463 G4cout << "Upon exiting LocateGlobalPointAndSetup():" << G4endl;
576 G4cout << "Upon exiting LocateGlobalPoin << 464 G4cout << " History = " << G4endl << fHistory << G4endl << G4endl;
577 G4cout << " History = " << G4endl << << 465 #endif
578 } <<
579 G4cout.precision(oldcoutPrec); 466 G4cout.precision(oldcoutPrec);
580 } 467 }
581 #endif 468 #endif
582 469
583 fLocatedOutsideWorld = false; << 470 fLocatedOutsideWorld= false;
584 471
585 return targetPhysical; 472 return targetPhysical;
586 } 473 }
587 474
588 // ******************************************* 475 // ********************************************************************
589 // LocateGlobalPointWithinVolume 476 // LocateGlobalPointWithinVolume
590 // 477 //
591 // -> the state information of this Navigator 478 // -> the state information of this Navigator and its subNavigators
592 // is updated in order to start the next st 479 // is updated in order to start the next step at pGlobalpoint
593 // -> no check is performed whether pGlobalpoi 480 // -> no check is performed whether pGlobalpoint is inside the
594 // original volume (this must be the case). 481 // original volume (this must be the case).
595 // 482 //
596 // Note: a direction could be added to the arg 483 // Note: a direction could be added to the arguments, to aid in future
597 // optional checking (via the old code b 484 // optional checking (via the old code below, flagged by OLD_LOCATE).
598 // [ This would be done only in verbose 485 // [ This would be done only in verbose mode ]
599 // ******************************************* 486 // ********************************************************************
600 // 487 //
601 void 488 void
602 G4Navigator::LocateGlobalPointWithinVolume(con 489 G4Navigator::LocateGlobalPointWithinVolume(const G4ThreeVector& pGlobalpoint)
603 { << 490 {
>> 491 fLastLocatedPointLocal = ComputeLocalPoint(pGlobalpoint);
>> 492
604 #ifdef G4DEBUG_NAVIGATION 493 #ifdef G4DEBUG_NAVIGATION
605 assert( !fWasLimitedByGeometry ); << 494 if( fVerbose > 2 )
606 // Check: Either step was not limited by a << 495 {
607 // else the full step is no longer << 496 G4cout << "Entering LocateGlobalWithinVolume(): History = " << G4endl;
>> 497 G4cout << fHistory << G4endl;
>> 498 }
608 #endif 499 #endif
609 <<
610 fLastLocatedPointLocal = ComputeLocalPoint( <<
611 fLastTriedStepComputation = false; <<
612 fChangedGrandMotherRefFrame = false; // F <<
613 500
614 // For the case of Voxel (or Parameterised) 501 // For the case of Voxel (or Parameterised) volume the respective
615 // Navigator must be messaged to update its 502 // Navigator must be messaged to update its voxel information etc
616 503
617 // Update the state of the Sub Navigators 504 // Update the state of the Sub Navigators
618 // - in particular any voxel information th 505 // - in particular any voxel information they store/cache
619 // 506 //
620 G4VPhysicalVolume* motherPhysical = fHisto 507 G4VPhysicalVolume* motherPhysical = fHistory.GetTopVolume();
621 G4LogicalVolume* motherLogical = mother 508 G4LogicalVolume* motherLogical = motherPhysical->GetLogicalVolume();
>> 509 G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
622 510
623 switch( CharacteriseDaughters(motherLogical << 511 if ( fHistory.GetTopVolumeType()!=kReplica )
624 { 512 {
>> 513 switch( CharacteriseDaughters(motherLogical) )
>> 514 {
625 case kNormal: 515 case kNormal:
626 GetVoxelNavigator().RelocateWithinVol << 516 if ( pVoxelHeader )
>> 517 {
>> 518 fvoxelNav.VoxelLocate( pVoxelHeader, fLastLocatedPointLocal );
>> 519 }
627 break; 520 break;
628 case kParameterised: 521 case kParameterised:
629 fparamNav.RelocateWithinVolume( mothe << 522 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
>> 523 {
>> 524 // Resets state & returns voxel node
>> 525 //
>> 526 fparamNav.ParamVoxelLocate( pVoxelHeader, fLastLocatedPointLocal );
>> 527 }
630 break; 528 break;
631 case kReplica: 529 case kReplica:
632 // Nothing to do << 530 G4Exception("G4Navigator::LocateGlobalPointWithinVolume()",
633 break; << 531 "NotApplicable", FatalException,
634 case kExternal: << 532 "Not applicable for replicated volumes.");
635 fpExternalNav->RelocateWithinVolume( <<
636 <<
637 break; 533 break;
>> 534 }
638 } 535 }
639 536
640 // Reset the state variables 537 // Reset the state variables
641 // - which would have been affected 538 // - which would have been affected
642 // by the 'equivalent' call to LocateGl 539 // by the 'equivalent' call to LocateGlobalPointAndSetup
643 // - who's values have been invalidated b 540 // - who's values have been invalidated by the 'move'.
644 // 541 //
645 fBlockedPhysicalVolume = nullptr; << 542 fBlockedPhysicalVolume = 0;
646 fBlockedReplicaNo = -1; 543 fBlockedReplicaNo = -1;
647 fEntering = false; 544 fEntering = false;
648 fEnteredDaughter = false; // Boundary not 545 fEnteredDaughter = false; // Boundary not encountered, did not enter
649 fExiting = false; 546 fExiting = false;
650 fExitedMother = false; // Boundary not 547 fExitedMother = false; // Boundary not encountered, did not exit
651 } 548 }
652 549
653 // ******************************************* 550 // ********************************************************************
654 // SetSavedState 551 // SetSavedState
655 // 552 //
656 // Save the state, in case this is a parasitic 553 // Save the state, in case this is a parasitic call
657 // Save fValidExitNormal, fExitNormal, fExitin 554 // Save fValidExitNormal, fExitNormal, fExiting, fEntering,
658 // fBlockedPhysicalVolume, fBlockedReplic 555 // fBlockedPhysicalVolume, fBlockedReplicaNo, fLastStepWasZero;
659 // ******************************************* 556 // ********************************************************************
660 // 557 //
661 void G4Navigator::SetSavedState() 558 void G4Navigator::SetSavedState()
662 { 559 {
663 // Note: the state of dependent objects is n << 560 // fSaveExitNormal = fExitNormal;
664 // ( This means that the full state is cha <<
665 // SetSavedState() and RestoreSavedState <<
666 <<
667 fSaveState.sExitNormal = fExitNormal; 561 fSaveState.sExitNormal = fExitNormal;
668 fSaveState.sValidExitNormal = fValidExitNorm 562 fSaveState.sValidExitNormal = fValidExitNormal;
669 fSaveState.sExiting = fExiting; 563 fSaveState.sExiting = fExiting;
670 fSaveState.sEntering = fEntering; 564 fSaveState.sEntering = fEntering;
671 565
672 fSaveState.spBlockedPhysicalVolume = fBlocke 566 fSaveState.spBlockedPhysicalVolume = fBlockedPhysicalVolume;
673 fSaveState.sBlockedReplicaNo = fBlockedRepli << 567 fSaveState.sBlockedReplicaNo = fBlockedReplicaNo,
674 568
675 fSaveState.sLastStepWasZero = static_cast<G4 << 569 fSaveState.sLastStepWasZero = fLastStepWasZero;
676 <<
677 fSaveState.sLocatedOutsideWorld = fLocatedOu <<
678 fSaveState.sLastLocatedPointLocal = fLastLoc <<
679 fSaveState.sEnteredDaughter = fEnteredDaught <<
680 fSaveState.sExitedMother = fExitedMother; <<
681 fSaveState.sWasLimitedByGeometry = fWasLimit <<
682 <<
683 // Even the safety sphere - if you want to c <<
684 // <<
685 fSaveState.sPreviousSftOrigin = fPreviousSft <<
686 fSaveState.sPreviousSafety = fPreviousSafety <<
687 } 570 }
688 571
689 // ******************************************* 572 // ********************************************************************
690 // RestoreSavedState 573 // RestoreSavedState
691 // 574 //
692 // Restore the state (in Compute Step), in cas 575 // Restore the state (in Compute Step), in case this is a parasitic call
693 // ******************************************* 576 // ********************************************************************
694 // 577 //
695 void G4Navigator::RestoreSavedState() 578 void G4Navigator::RestoreSavedState()
696 { 579 {
697 fExitNormal = fSaveState.sExitNormal; 580 fExitNormal = fSaveState.sExitNormal;
698 fValidExitNormal = fSaveState.sValidExitNorm 581 fValidExitNormal = fSaveState.sValidExitNormal;
699 fExiting = fSaveState.sExiting; 582 fExiting = fSaveState.sExiting;
700 fEntering = fSaveState.sEntering; 583 fEntering = fSaveState.sEntering;
701 584
702 fBlockedPhysicalVolume = fSaveState.spBlocke 585 fBlockedPhysicalVolume = fSaveState.spBlockedPhysicalVolume;
703 fBlockedReplicaNo = fSaveState.sBlockedRepli << 586 fBlockedReplicaNo = fSaveState.sBlockedReplicaNo,
704 587
705 fLastStepWasZero = (fSaveState.sLastStepWasZ << 588 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 } 589 }
717 590
718 // ******************************************* 591 // ********************************************************************
719 // ComputeStep 592 // ComputeStep
720 // 593 //
721 // Computes the next geometric Step: intersect 594 // Computes the next geometric Step: intersections with current
722 // mother and `daughter' volumes. 595 // mother and `daughter' volumes.
723 // 596 //
724 // NOTE: 597 // NOTE:
725 // 598 //
726 // Flags on entry: 599 // Flags on entry:
727 // -------------- 600 // --------------
728 // fValidExitNormal - Normal of exited volume 601 // fValidExitNormal - Normal of exited volume is valid (convex, not a
729 // coincident boundary) 602 // coincident boundary)
730 // fExitNormal - Surface normal of exite 603 // fExitNormal - Surface normal of exited volume
731 // fExiting - True if have exited sol 604 // fExiting - True if have exited solid
732 // 605 //
733 // fBlockedPhysicalVolume - Ptr to exited volu 606 // fBlockedPhysicalVolume - Ptr to exited volume (or 0)
734 // fBlockedReplicaNo - Replication no of exite 607 // fBlockedReplicaNo - Replication no of exited volume
735 // fLastStepWasZero - True if last Step size << 608 // fLastStepWasZero - True if last Step size was zero.
736 // 609 //
737 // Flags on exit: 610 // Flags on exit:
738 // ------------- 611 // -------------
739 // fValidExitNormal - True if surface normal 612 // fValidExitNormal - True if surface normal of exited volume is valid
740 // fExitNormal - Surface normal of exite 613 // fExitNormal - Surface normal of exited volume rotated to mothers
741 // reference system 614 // reference system
742 // fExiting - True if exiting mother 615 // fExiting - True if exiting mother
743 // fEntering - True if entering `daugh 616 // fEntering - True if entering `daughter' volume (or replica)
744 // fBlockedPhysicalVolume - Ptr to candidate ( 617 // fBlockedPhysicalVolume - Ptr to candidate (entered) volume
745 // fBlockedReplicaNo - Replication no of candi 618 // fBlockedReplicaNo - Replication no of candidate (entered) volume
746 // fLastStepWasZero - True if this Step size << 619 // fLastStepWasZero - True if this Step size was zero.
747 // ******************************************* 620 // ********************************************************************
748 // 621 //
749 G4double G4Navigator::ComputeStep( const G4Thr << 622 G4double G4Navigator::ComputeStep( const G4ThreeVector &pGlobalpoint,
750 const G4Thr << 623 const G4ThreeVector &pDirection,
751 const G4dou 624 const G4double pCurrentProposedStepLength,
752 G4dou << 625 G4double &pNewSafety)
753 { 626 {
754 #ifdef G4DEBUG_NAVIGATION <<
755 static G4ThreadLocal G4int sNavCScalls = 0; <<
756 ++sNavCScalls; <<
757 #endif <<
758 <<
759 G4ThreeVector localDirection = ComputeLocalA 627 G4ThreeVector localDirection = ComputeLocalAxis(pDirection);
760 G4double Step = kInfinity; 628 G4double Step = kInfinity;
761 G4VPhysicalVolume *motherPhysical = fHistor 629 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
762 G4LogicalVolume *motherLogical = motherPhysi 630 G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
763 631
764 // All state relating to exiting normals mus << 632 static G4int sNavCScalls=0;
765 // << 633 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 634
774 #ifdef G4VERBOSE 635 #ifdef G4VERBOSE
775 if( fVerbose > 0 ) 636 if( fVerbose > 0 )
776 { 637 {
777 G4cout << "*** G4Navigator::ComputeStep: * 638 G4cout << "*** G4Navigator::ComputeStep: ***" << G4endl;
778 G4cout << " Volume = " << motherPhysica 639 G4cout << " Volume = " << motherPhysical->GetName()
779 << " - Proposed step length = " << 640 << " - Proposed step length = " << pCurrentProposedStepLength
780 << G4endl; 641 << G4endl;
781 #ifdef G4DEBUG_NAVIGATION 642 #ifdef G4DEBUG_NAVIGATION
782 if( fVerbose >= 2 ) << 643 if( fVerbose >= 4 )
783 { 644 {
784 G4cout << " Called with the arguments: 645 G4cout << " Called with the arguments: " << G4endl
785 << " Globalpoint = " << std::set 646 << " Globalpoint = " << std::setw(25) << pGlobalpoint << G4endl
786 << " Direction = " << std::set 647 << " Direction = " << std::setw(25) << pDirection << G4endl;
787 if( fVerbose >= 4 ) << 648 G4cout << " ---- Upon entering :" << G4endl;
788 { << 649 PrintState();
789 G4cout << " ---- Upon entering : Stat <<
790 PrintState(); <<
791 } <<
792 } 650 }
793 #endif 651 #endif
794 } 652 }
795 #endif 653 #endif
796 654
797 G4ThreeVector newLocalPoint = ComputeLocalPo 655 G4ThreeVector newLocalPoint = ComputeLocalPoint(pGlobalpoint);
798 <<
799 if( newLocalPoint != fLastLocatedPointLocal 656 if( newLocalPoint != fLastLocatedPointLocal )
800 { 657 {
801 // Check whether the relocation is within 658 // Check whether the relocation is within safety
802 // 659 //
803 G4ThreeVector oldLocalPoint = fLastLocated 660 G4ThreeVector oldLocalPoint = fLastLocatedPointLocal;
804 G4double moveLenSq = (newLocalPoint-oldLoc 661 G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2();
805 662
806 if ( moveLenSq >= fSqTol ) << 663 if ( moveLenSq >= kCarTolerance*kCarTolerance )
807 { 664 {
808 #ifdef G4VERBOSE 665 #ifdef G4VERBOSE
809 ComputeStepLog(pGlobalpoint, moveLenSq); 666 ComputeStepLog(pGlobalpoint, moveLenSq);
810 #endif 667 #endif
811 // Relocate the point within the same vo 668 // Relocate the point within the same volume
812 // 669 //
813 LocateGlobalPointWithinVolume( pGlobalpo 670 LocateGlobalPointWithinVolume( pGlobalpoint );
814 } 671 }
815 } 672 }
816 if ( fHistory.GetTopVolumeType()!=kReplica ) 673 if ( fHistory.GetTopVolumeType()!=kReplica )
817 { 674 {
818 switch( CharacteriseDaughters(motherLogica 675 switch( CharacteriseDaughters(motherLogical) )
819 { 676 {
820 case kNormal: 677 case kNormal:
821 if ( motherLogical->GetVoxelHeader() ! << 678 if ( motherLogical->GetVoxelHeader() )
822 { 679 {
823 Step = GetVoxelNavigator().ComputeSt << 680 Step = fvoxelNav.ComputeStep(fLastLocatedPointLocal,
824 localDi 681 localDirection,
825 pCurren 682 pCurrentProposedStepLength,
826 pNewSaf 683 pNewSafety,
827 fHistor 684 fHistory,
828 fValidE 685 fValidExitNormal,
829 fExitNo 686 fExitNormal,
830 fExitin 687 fExiting,
831 fEnteri 688 fEntering,
832 &fBlock 689 &fBlockedPhysicalVolume,
833 fBlocke 690 fBlockedReplicaNo);
834 691
835 } 692 }
836 else 693 else
837 { 694 {
838 if( motherPhysical->GetRegularStruct 695 if( motherPhysical->GetRegularStructureId() == 0 )
839 { 696 {
840 Step = fnormalNav.ComputeStep(fLas 697 Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
841 loca 698 localDirection,
842 pCur 699 pCurrentProposedStepLength,
843 pNew 700 pNewSafety,
844 fHis 701 fHistory,
845 fVal 702 fValidExitNormal,
846 fExi 703 fExitNormal,
847 fExi 704 fExiting,
848 fEnt 705 fEntering,
849 &fBl 706 &fBlockedPhysicalVolume,
850 fBlo 707 fBlockedReplicaNo);
851 } 708 }
852 else // Regular (non-voxelised) str 709 else // Regular (non-voxelised) structure
853 { 710 {
854 LocateGlobalPointAndSetup( pGlobal 711 LocateGlobalPointAndSetup( pGlobalpoint, &pDirection, true, true );
855 // 712 //
856 // if physical process limits the 713 // if physical process limits the step, the voxel will not be the
857 // one given by ComputeStepSkippin 714 // one given by ComputeStepSkippingEqualMaterials() and the local
858 // point will be wrongly calculate 715 // point will be wrongly calculated.
859 716
860 // There is a problem: when msc li 717 // There is a problem: when msc limits the step and the point is
861 // assigned wrongly to phantom in 718 // assigned wrongly to phantom in previous step (while it is out
862 // of the container volume). Then 719 // of the container volume). Then LocateGlobalPointAndSetup() has
863 // reset the history topvolume to 720 // reset the history topvolume to world.
864 // 721 //
865 if(fHistory.GetTopVolume()->GetReg 722 if(fHistory.GetTopVolume()->GetRegularStructureId() == 0 )
866 { 723 {
867 G4Exception("G4Navigator::Comput 724 G4Exception("G4Navigator::ComputeStep()",
868 "GeomNav1001", JustW << 725 "Bad-location-of-point", JustWarning,
869 "Point is relocated in voxels, 726 "Point is relocated in voxels, while it should be outside!");
870 Step = fnormalNav.ComputeStep(fL 727 Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
871 lo 728 localDirection,
872 pC 729 pCurrentProposedStepLength,
873 pN 730 pNewSafety,
874 fH 731 fHistory,
875 fV 732 fValidExitNormal,
876 fE 733 fExitNormal,
877 fE 734 fExiting,
878 fE 735 fEntering,
879 &f 736 &fBlockedPhysicalVolume,
880 fB 737 fBlockedReplicaNo);
881 } 738 }
882 else 739 else
883 { 740 {
884 Step = fregularNav. 741 Step = fregularNav.
885 ComputeStepSkippingEqualMat 742 ComputeStepSkippingEqualMaterials(fLastLocatedPointLocal,
886 743 localDirection,
887 744 pCurrentProposedStepLength,
888 745 pNewSafety,
889 746 fHistory,
890 747 fValidExitNormal,
891 748 fExitNormal,
892 749 fExiting,
893 750 fEntering,
894 751 &fBlockedPhysicalVolume,
895 752 fBlockedReplicaNo,
896 753 motherPhysical);
897 } 754 }
898 } 755 }
899 } 756 }
900 break; 757 break;
901 case kParameterised: 758 case kParameterised:
902 if( GetDaughtersRegularStructureId(mot 759 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
903 { 760 {
904 Step = fparamNav.ComputeStep(fLastLo 761 Step = fparamNav.ComputeStep(fLastLocatedPointLocal,
905 localDi 762 localDirection,
906 pCurren 763 pCurrentProposedStepLength,
907 pNewSaf 764 pNewSafety,
908 fHistor 765 fHistory,
909 fValidE 766 fValidExitNormal,
910 fExitNo 767 fExitNormal,
911 fExitin 768 fExiting,
912 fEnteri 769 fEntering,
913 &fBlock 770 &fBlockedPhysicalVolume,
914 fBlocke 771 fBlockedReplicaNo);
915 } 772 }
916 else // Regular structure 773 else // Regular structure
917 { 774 {
918 Step = fregularNav.ComputeStep(fLast 775 Step = fregularNav.ComputeStep(fLastLocatedPointLocal,
919 local 776 localDirection,
920 pCurr 777 pCurrentProposedStepLength,
921 pNewS 778 pNewSafety,
922 fHist 779 fHistory,
923 fVali 780 fValidExitNormal,
924 fExit 781 fExitNormal,
925 fExit 782 fExiting,
926 fEnte 783 fEntering,
927 &fBlo 784 &fBlockedPhysicalVolume,
928 fBloc 785 fBlockedReplicaNo);
929 } 786 }
930 break; 787 break;
931 case kReplica: 788 case kReplica:
932 G4Exception("G4Navigator::ComputeStep( << 789 G4Exception("G4Navigator::ComputeStep()", "NotApplicable",
933 FatalException, "Not appli 790 FatalException, "Not applicable for replicated volumes.");
934 break; 791 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 } 792 }
949 } 793 }
950 else 794 else
951 { 795 {
952 // In the case of a replica, it must handl 796 // In the case of a replica, it must handle the exiting
953 // edge/corner problem by itself 797 // edge/corner problem by itself
954 // 798 //
955 fExiting = fExitedMother; << 799 G4bool exitingReplica = fExitedMother;
956 Step = freplicaNav.ComputeStep(pGlobalpoin 800 Step = freplicaNav.ComputeStep(pGlobalpoint,
957 pDirection, 801 pDirection,
958 fLastLocate 802 fLastLocatedPointLocal,
959 localDirect 803 localDirection,
960 pCurrentPro 804 pCurrentProposedStepLength,
961 pNewSafety, 805 pNewSafety,
962 fHistory, 806 fHistory,
963 fValidExitN 807 fValidExitNormal,
964 fCalculated <<
965 fExitNormal 808 fExitNormal,
966 fExiting, << 809 exitingReplica,
967 fEntering, 810 fEntering,
968 &fBlockedPh 811 &fBlockedPhysicalVolume,
969 fBlockedRep 812 fBlockedReplicaNo);
>> 813 fExiting= exitingReplica; // still ok to set it ??
970 } 814 }
971 815
972 // Remember last safety origin & value. 816 // Remember last safety origin & value.
973 // 817 //
974 fPreviousSftOrigin = pGlobalpoint; 818 fPreviousSftOrigin = pGlobalpoint;
975 fPreviousSafety = pNewSafety; 819 fPreviousSafety = pNewSafety;
976 820
977 // Count zero steps - one can occur due to c 821 // Count zero steps - one can occur due to changing momentum at a boundary
978 // - one, two (or a few) ca 822 // - one, two (or a few) can occur at common edges between
979 // volumes 823 // volumes
980 // - more than two is likel 824 // - more than two is likely a problem in the geometry
981 // description or the Nav 825 // description or the Navigation
982 826
983 // Rule of thumb: likely at an Edge if two c 827 // Rule of thumb: likely at an Edge if two consecutive steps are zero,
984 // because at least two candi 828 // because at least two candidate volumes must have been
985 // checked 829 // checked
986 // 830 //
987 fLocatedOnEdge = fLastStepWasZero && (Step 831 fLocatedOnEdge = fLastStepWasZero && (Step==0.0);
988 fLastStepWasZero = (Step<fMinStep); << 832 fLastStepWasZero = (Step==0.0);
989 if (fPushed) { fPushed = fLastStepWasZero; << 833 if (fPushed) fPushed = fLastStepWasZero;
990 834
991 // Handle large number of consecutive zero s 835 // Handle large number of consecutive zero steps
992 // 836 //
993 if ( fLastStepWasZero ) 837 if ( fLastStepWasZero )
994 { 838 {
995 ++fNumberZeroSteps; << 839 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 840 #ifdef G4DEBUG_NAVIGATION
1005 inform = fNumberZeroSteps > 1; << 841 if( fNumberZeroSteps > 1 )
1006 #endif <<
1007 <<
1008 if ( act || inform ) <<
1009 { 842 {
1010 if( act && !abandon ) << 843 G4cout << "G4Navigator::ComputeStep(): another zero step, # "
1011 { << 844 << fNumberZeroSteps
1012 // Act to recover this stuck track. P << 845 << " at " << pGlobalpoint
1013 // << 846 << " in volume " << motherPhysical->GetName()
1014 Step += 100*kCarTolerance; << 847 << " nav-comp-step calls # " << sNavCScalls
1015 fPushed = true; << 848 << G4endl;
1016 } << 849 }
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 850 #endif
1047 message << " Track *abandoned* due << 851 if( fNumberZeroSteps > fActionThreshold_NoZeroSteps-1 )
1048 << " Event aborted. " << G4 << 852 {
1049 G4Exception("G4Navigator::ComputeSt << 853 // Act to recover this stuck track. Pushing it along direction
1050 EventMustBeAborted, mes << 854 //
1051 } << 855 Step += 100*kCarTolerance;
1052 else <<
1053 { <<
1054 #ifdef G4VERBOSE 856 #ifdef G4VERBOSE
1055 if ( actAndReport ) // (!fPushed = << 857 if ((!fPushed) && (fWarnPush))
1056 { << 858 {
1057 message << " *** Trying to get << 859 G4cerr << "WARNING - G4Navigator::ComputeStep()" << G4endl
1058 << " - expanding step to << 860 << " Track stuck, not moving for "
1059 << " Potential ove << 861 << fNumberZeroSteps << " steps" << G4endl
1060 G4Exception("G4Navigator::Comput << 862 << " in volume -" << motherPhysical->GetName()
1061 JustWarning, message << 863 << "- at point " << pGlobalpoint << G4endl
1062 } << 864 << " direction: " << pDirection << "." << G4endl
1063 #endif << 865 << " Potential geometry or navigation problem !"
1064 #ifdef G4DEBUG_NAVIGATION << 866 << G4endl
1065 else << 867 << " Trying pushing it of " << Step << " mm ..."
1066 { << 868 << G4endl;
1067 if( fNumberZeroSteps > 1 ) << 869 }
1068 { <<
1069 message << ", nav-comp-step ca <<
1070 << ", Step= " << Step <<
1071 G4cout << message.str(); <<
1072 } <<
1073 } <<
1074 #endif 870 #endif
1075 } // end of else if ( abandon ) << 871 fPushed = true;
1076 } // end of if( actAndReport || abandon << 872 }
1077 } // end of if ( act || inform ) << 873 if( fNumberZeroSteps > fAbandonThreshold_NoZeroSteps-1 )
>> 874 {
>> 875 // Must kill this stuck track
>> 876 //
>> 877 G4cerr << "ERROR - G4Navigator::ComputeStep()" << G4endl
>> 878 << " Track stuck, not moving for "
>> 879 << fNumberZeroSteps << " steps" << G4endl
>> 880 << " in volume -" << motherPhysical->GetName()
>> 881 << "- at point " << pGlobalpoint << G4endl
>> 882 << " direction: " << pDirection << "." << G4endl;
>> 883 motherPhysical->CheckOverlaps(5000, false);
>> 884 G4Exception("G4Navigator::ComputeStep()",
>> 885 "StuckTrack", EventMustBeAborted,
>> 886 "Stuck Track: potential geometry or navigation problem.");
>> 887 }
1078 } 888 }
1079 else 889 else
1080 { 890 {
1081 if (!fPushed) { fNumberZeroSteps = 0; } << 891 if (!fPushed) fNumberZeroSteps = 0;
1082 } 892 }
1083 fLastMotherPhys = motherPhysical; <<
1084 893
1085 fEnteredDaughter = fEntering; // I expect 894 fEnteredDaughter = fEntering; // I expect to enter a volume in this Step
1086 fExitedMother = fExiting; 895 fExitedMother = fExiting;
1087 896
1088 fStepEndPoint = pGlobalpoint <<
1089 + std::min(Step,pCurrentPropo <<
1090 fLastStepEndPointLocal = fLastLocatedPointL <<
1091 <<
1092 if( fExiting ) 897 if( fExiting )
1093 { 898 {
1094 #ifdef G4DEBUG_NAVIGATION 899 #ifdef G4DEBUG_NAVIGATION
1095 if( fVerbose > 2 ) 900 if( fVerbose > 2 )
1096 { 901 {
1097 G4cout << " At G4Nav CompStep End - if( << 902 G4cout << " At G4Nav CompStep End - if(exiting) - fExiting= " << fExiting
1098 << " fValidExitNormal = " << fVa 903 << " fValidExitNormal = " << fValidExitNormal << G4endl;
1099 G4cout << " fExitNormal= " << fExitNorm 904 G4cout << " fExitNormal= " << fExitNormal << G4endl;
1100 } 905 }
1101 #endif 906 #endif
1102 907
1103 if ( fValidExitNormal || fCalculatedExitN << 908 if(fValidExitNormal)
1104 { 909 {
1105 // Convention: fExitNormal is in the 'g 910 // Convention: fExitNormal is in the 'grand-mother' coordinate system
1106 fGrandMotherExitNormal = fExitNormal; << 911 //
>> 912 fGrandMotherExitNormal= fExitNormal;
1107 } 913 }
1108 else 914 else
1109 { 915 {
1110 // We must calculate the normal anyway 916 // We must calculate the normal anyway (in order to have it if requested)
1111 // 917 //
1112 G4ThreeVector finalLocalPoint = fLastLo << 918 G4ThreeVector finalLocalPoint =
1113 + localDi << 919 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 920
1135 // Do not set fValidExitNormal -- thi << 921 // Now fGrandMotherExitNormal is in the 'grand-mother' coordinate system
1136 // that the solid is convex! << 922 //
1137 } << 923 fGrandMotherExitNormal =
1138 else << 924 motherLogical->GetSolid()->SurfaceNormal(finalLocalPoint);
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 <<
1162 if ( fHistory.GetTopVolumeType() != kRepl <<
1163 { <<
1164 fCalculatedExitNormal = true; <<
1165 } <<
1166 925
1167 // Now transform it to the global referen << 926 const G4RotationMatrix* mRot = motherPhysical->GetRotation();
1168 // << 927 if( mRot )
1169 if( fValidExitNormal || fCalculatedExitNo << 928 {
1170 { << 929 fGrandMotherExitNormal *= (*mRot).inverse();
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 } 930 }
1181 } 931 }
1182 else <<
1183 { <<
1184 fExitNormalGlobalFrame = G4ThreeVector( <<
1185 } <<
1186 } 932 }
>> 933 fStepEndPoint= pGlobalpoint+Step*pDirection;
1187 934
1188 if( (Step == pCurrentProposedStepLength) && 935 if( (Step == pCurrentProposedStepLength) && (!fExiting) && (!fEntering) )
1189 { 936 {
1190 // This if Step is not really limited by 937 // This if Step is not really limited by the geometry.
1191 // The Navigator is obliged to return "in 938 // The Navigator is obliged to return "infinity"
1192 // 939 //
1193 Step = kInfinity; 940 Step = kInfinity;
1194 } 941 }
1195 942
1196 #ifdef G4VERBOSE 943 #ifdef G4VERBOSE
1197 if( fVerbose > 1 ) 944 if( fVerbose > 1 )
1198 { 945 {
1199 if( fVerbose >= 4 ) 946 if( fVerbose >= 4 )
1200 { 947 {
1201 G4cout << " ----- Upon exiting :" << 948 G4cout << " ----- Upon exiting :" << G4endl;
1202 PrintState(); 949 PrintState();
1203 } 950 }
1204 G4cout << " Returned step= " << Step; << 951 G4cout <<" Returned step = " << Step << G4endl;
1205 if( fVerbose > 5 ) { G4cout << G4endl; } <<
1206 if( Step == kInfinity ) 952 if( Step == kInfinity )
1207 { 953 {
1208 G4cout << " Requested step= " << pCurr << 954 G4cout << " Original proposed step = "
1209 if( fVerbose > 5) { G4cout << G4endl; << 955 << pCurrentProposedStepLength << G4endl;
1210 } 956 }
1211 G4cout << " Safety = " << pNewSafety << << 957 G4cout << " Safety = " << pNewSafety << G4endl;
1212 } 958 }
1213 #endif 959 #endif
1214 960
1215 fLastTriedStepComputation = true; <<
1216 <<
1217 return Step; 961 return Step;
1218 } 962 }
1219 963
1220 // ****************************************** 964 // ********************************************************************
1221 // CheckNextStep 965 // CheckNextStep
1222 // 966 //
1223 // Compute the step without altering the navi 967 // Compute the step without altering the navigator state
1224 // ****************************************** 968 // ********************************************************************
1225 // 969 //
1226 G4double G4Navigator::CheckNextStep( const G4 970 G4double G4Navigator::CheckNextStep( const G4ThreeVector& pGlobalpoint,
1227 const G4 971 const G4ThreeVector& pDirection,
1228 const G4 972 const G4double pCurrentProposedStepLength,
1229 G4 973 G4double& pNewSafety)
1230 { 974 {
1231 G4double step; 975 G4double step;
1232 976
1233 // Save the state, for this parasitic call 977 // Save the state, for this parasitic call
1234 // 978 //
1235 SetSavedState(); 979 SetSavedState();
1236 980
1237 step = ComputeStep ( pGlobalpoint, 981 step = ComputeStep ( pGlobalpoint,
1238 pDirection, 982 pDirection,
1239 pCurrentProposedStepLe 983 pCurrentProposedStepLength,
1240 pNewSafety ); 984 pNewSafety );
1241 985
1242 // It is a parasitic call, so attempt to re << 986 // If a parasitic call, then attempt to restore the key parts of the state
1243 // 987 //
1244 RestoreSavedState(); 988 RestoreSavedState();
1245 // NOTE: the state of the current subnaviga << 989
1246 // ***> TODO: restore subnavigator state <<
1247 // if( last_located) Need <<
1248 // if( last_computed step) Need <<
1249 <<
1250 return step; 990 return step;
1251 } 991 }
1252 992
1253 // ****************************************** 993 // ********************************************************************
1254 // ResetState 994 // ResetState
1255 // 995 //
1256 // Resets stack and minimum of navigator stat 996 // Resets stack and minimum of navigator state `machine'
1257 // ****************************************** 997 // ********************************************************************
1258 // 998 //
1259 void G4Navigator::ResetState() 999 void G4Navigator::ResetState()
1260 { 1000 {
1261 fWasLimitedByGeometry = false; 1001 fWasLimitedByGeometry = false;
1262 fEntering = false; 1002 fEntering = false;
1263 fExiting = false; 1003 fExiting = false;
1264 fLocatedOnEdge = false; 1004 fLocatedOnEdge = false;
1265 fLastStepWasZero = false; 1005 fLastStepWasZero = false;
1266 fEnteredDaughter = false; 1006 fEnteredDaughter = false;
1267 fExitedMother = false; 1007 fExitedMother = false;
1268 fPushed = false; 1008 fPushed = false;
1269 1009
1270 fValidExitNormal = false; 1010 fValidExitNormal = false;
1271 fChangedGrandMotherRefFrame = false; <<
1272 fCalculatedExitNormal = false; <<
1273 <<
1274 fExitNormal = G4ThreeVector(0,0, 1011 fExitNormal = G4ThreeVector(0,0,0);
1275 fGrandMotherExitNormal = G4ThreeVector(0,0, <<
1276 fExitNormalGlobalFrame = G4ThreeVector(0,0, <<
1277 1012
1278 fPreviousSftOrigin = G4ThreeVector(0,0, 1013 fPreviousSftOrigin = G4ThreeVector(0,0,0);
1279 fPreviousSafety = 0.0; 1014 fPreviousSafety = 0.0;
1280 1015
1281 fNumberZeroSteps = 0; 1016 fNumberZeroSteps = 0;
1282 << 1017
1283 fBlockedPhysicalVolume = nullptr; << 1018 fBlockedPhysicalVolume = 0;
1284 fBlockedReplicaNo = -1; 1019 fBlockedReplicaNo = -1;
1285 1020
1286 fLastLocatedPointLocal = G4ThreeVector( kIn 1021 fLastLocatedPointLocal = G4ThreeVector( kInfinity, -kInfinity, 0.0 );
1287 fLocatedOutsideWorld = false; 1022 fLocatedOutsideWorld = false;
1288 <<
1289 fLastMotherPhys = nullptr; <<
1290 } 1023 }
1291 1024
1292 // ****************************************** 1025 // ********************************************************************
1293 // SetupHierarchy 1026 // SetupHierarchy
1294 // 1027 //
1295 // Renavigates & resets hierarchy described b 1028 // Renavigates & resets hierarchy described by current history
1296 // o Reset volumes 1029 // o Reset volumes
1297 // o Recompute transforms and/or solids of re 1030 // o Recompute transforms and/or solids of replicated/parameterised volumes
1298 // ****************************************** 1031 // ********************************************************************
1299 // 1032 //
1300 void G4Navigator::SetupHierarchy() 1033 void G4Navigator::SetupHierarchy()
1301 { 1034 {
1302 const auto depth = (G4int)fHistory.GetDept << 1035 G4int i;
1303 for ( auto i = 1; i <= depth; ++i ) << 1036 const G4int cdepth = fHistory.GetDepth();
>> 1037 G4VPhysicalVolume *current;
>> 1038 G4VSolid *pSolid;
>> 1039 G4VPVParameterisation *pParam;
>> 1040
>> 1041 for ( i=1; i<=cdepth; i++ )
1304 { 1042 {
>> 1043 current = fHistory.GetVolume(i);
1305 switch ( fHistory.GetVolumeType(i) ) 1044 switch ( fHistory.GetVolumeType(i) )
1306 { 1045 {
1307 case kNormal: 1046 case kNormal:
1308 case kExternal: <<
1309 break; 1047 break;
1310 case kReplica: 1048 case kReplica:
1311 freplicaNav.ComputeTransformation(fHi << 1049 freplicaNav.ComputeTransformation(fHistory.GetReplicaNo(i), current);
1312 break; 1050 break;
1313 case kParameterised: 1051 case kParameterised:
1314 G4VPhysicalVolume* current = fHistory << 1052 G4int replicaNo;
1315 G4int replicaNo = fHistory.GetReplica << 1053 pParam = current->GetParameterisation();
1316 G4VPVParameterisation* pParam = curre << 1054 replicaNo = fHistory.GetReplicaNo(i);
1317 G4VSolid* pSolid = pParam->ComputeSol << 1055 pSolid = pParam->ComputeSolid(replicaNo, current);
1318 1056
1319 // Set up dimensions & transform in s 1057 // Set up dimensions & transform in solid/physical volume
1320 // 1058 //
1321 pSolid->ComputeDimensions(pParam, rep 1059 pSolid->ComputeDimensions(pParam, replicaNo, current);
1322 pParam->ComputeTransformation(replica 1060 pParam->ComputeTransformation(replicaNo, current);
1323 1061
1324 G4TouchableHistory* pTouchable = null << 1062 G4TouchableHistory touchable( fHistory );
1325 if( pParam->IsNested() ) << 1063 touchable.MoveUpHistory(); // move up to the parent level
1326 { << 1064
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 1065 // Set up the correct solid and material in Logical Volume
1337 // 1066 //
1338 G4LogicalVolume* pLogical = current-> << 1067 G4LogicalVolume *pLogical = current->GetLogicalVolume();
1339 pLogical->SetSolid( pSolid ); 1068 pLogical->SetSolid( pSolid );
1340 pLogical->UpdateMaterial( pParam -> 1069 pLogical->UpdateMaterial( pParam ->
1341 ComputeMaterial(replicaNo, current, << 1070 ComputeMaterial(replicaNo, current, &touchable) );
1342 delete pTouchable; <<
1343 break; 1071 break;
1344 } 1072 }
1345 } 1073 }
1346 } 1074 }
1347 1075
1348 // ****************************************** 1076 // ********************************************************************
1349 // GetLocalExitNormal 1077 // GetLocalExitNormal
1350 // 1078 //
1351 // Obtains the Normal vector to a surface (in 1079 // Obtains the Normal vector to a surface (in local coordinates)
1352 // pointing out of previous volume and into c 1080 // pointing out of previous volume and into current volume
1353 // ****************************************** 1081 // ********************************************************************
1354 // 1082 //
1355 G4ThreeVector G4Navigator::GetLocalExitNormal 1083 G4ThreeVector G4Navigator::GetLocalExitNormal( G4bool* valid )
1356 { 1084 {
1357 G4ThreeVector ExitNormal(0.,0.,0.); << 1085 G4ThreeVector ExitNormal(0.,0.,0.);
1358 G4VSolid* currentSolid = nullptr; <<
1359 G4LogicalVolume* candidateLogical; <<
1360 1086
1361 if ( fLastTriedStepComputation ) << 1087 if ( EnteredDaughterVolume() )
1362 { 1088 {
1363 // use fLastLocatedPointLocal and next ca << 1089 ExitNormal= -(fHistory.GetTopVolume()->GetLogicalVolume()->
1364 // << 1090 GetSolid()->SurfaceNormal(fLastLocatedPointLocal));
1365 G4ThreeVector nextSolidExitNormal(0.,0.,0 << 1091 *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 } 1092 }
1472 else // ( ! fLastTriedStepComputation ) i. << 1093 else
1473 { 1094 {
1474 if ( EnteredDaughterVolume() ) << 1095 if( fExitedMother )
1475 { 1096 {
1476 G4VSolid* daughterSolid = fHistory.GetT << 1097 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; 1098 *valid = true;
1490 } 1099 }
1491 else 1100 else
1492 { 1101 {
1493 if( fExitedMother ) << 1102 // We are not at a boundary.
1494 { << 1103 // ExitNormal remains (0,0,0)
1495 ExitNormal = fGrandMotherExitNormal; << 1104 //
1496 *valid = true; << 1105 *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 } 1106 }
1510 } 1107 }
1511 return ExitNormal; 1108 return ExitNormal;
1512 } 1109 }
1513 1110
1514 // ****************************************** 1111 // ********************************************************************
1515 // GetMotherToDaughterTransform << 1112 // ComputeSafety
1516 // 1113 //
1517 // Obtains the mother to daughter affine tran << 1114 // It assumes that it will be
>> 1115 // i) called at the Point in the same volume as the EndPoint of the
>> 1116 // ComputeStep.
>> 1117 // ii) after (or at the end of) ComputeStep OR after the relocation.
1518 // ****************************************** 1118 // ********************************************************************
1519 // 1119 //
1520 G4AffineTransform << 1120 G4double G4Navigator::ComputeSafety( const G4ThreeVector &pGlobalpoint,
1521 G4Navigator::GetMotherToDaughterTransform( G4 << 1121 const G4double pMaxLength,
1522 G4 << 1122 const G4bool keepState)
1523 EV <<
1524 { 1123 {
1525 switch (enteringVolumeType) << 1124 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 <<
1562 1125
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 1126 #ifdef G4DEBUG_NAVIGATION
1574 const G4ThreeVecto << 1127 G4int oldcoutPrec = G4cout.precision(8);
1575 #else << 1128 if( fVerbose > 0 )
1576 const G4ThreeVecto <<
1577 #endif <<
1578 G4bool* pVal <<
1579 { <<
1580 #ifdef G4DEBUG_NAVIGATION <<
1581 // Check Current point against expected 'lo <<
1582 // <<
1583 if ( fLastTriedStepComputation ) <<
1584 { 1129 {
1585 G4ThreeVector ExpectedBoundaryPointLocal; << 1130 G4cout << "*** G4Navigator::ComputeSafety: ***" << G4endl
>> 1131 << " Called at point: " << pGlobalpoint << G4endl;
1586 1132
1587 const G4AffineTransform& GlobalToLocal = << 1133 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1588 ExpectedBoundaryPointLocal = << 1134 G4cout << " Volume = " << motherPhysical->GetName()
1589 GlobalToLocal.TransformPoint( ExpectedB << 1135 << " - Maximum length = " << pMaxLength << G4endl;
1590 << 1136 if( fVerbose >= 4 )
1591 // Add here: Comparison against expected << 1137 {
1592 // i.e. the endpoint of Comput << 1138 G4cout << " ----- Upon entering Compute Safety:" << G4endl;
>> 1139 PrintState();
>> 1140 }
1593 } 1141 }
1594 #endif 1142 #endif
1595 <<
1596 return GetLocalExitNormal( pValid ); <<
1597 } <<
1598 1143
1599 // ****************************************** << 1144 if (keepState) { SetSavedState(); }
1600 // GetGlobalExitNormal << 1145
1601 // << 1146 G4double distEndpointSq = (pGlobalpoint-fStepEndPoint).mag2();
1602 // Obtains the Normal vector to a surface (in << 1147 G4bool stayedOnEndpoint = distEndpointSq < kCarTolerance*kCarTolerance;
1603 // pointing out of previous volume and into c << 1148 G4bool endpointOnSurface = fEnteredDaughter || fExitedMother;
1604 // ****************************************** << 1149
1605 // << 1150 if( !(endpointOnSurface && stayedOnEndpoint) )
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 { 1151 {
1623 // This was computed in last call to Comp << 1152 // Pseudo-relocate to this point (updates voxel information only)
1624 // and only if it arrived at boundary <<
1625 // 1153 //
1626 globalNormal = fExitNormalGlobalFrame; << 1154 LocateGlobalPointWithinVolume( pGlobalpoint );
1627 G4double normMag2 = globalNormal.mag2(); << 1155 // --->> Danger: Side effects on sub-navigator voxel information <<---
1628 if( std::fabs ( normMag2 - 1.0 ) < perTho << 1156 // Could be replaced again by 'granular' calls to sub-navigator
>> 1157 // locates (similar side-effects, but faster.
>> 1158 // Solutions:
>> 1159 // 1) Re-locate (to where?)
>> 1160 // 2) Insure that the methods using (G4ComputeStep?)
>> 1161 // does a relocation (if information is disturbed only ?)
>> 1162
>> 1163 #ifdef G4DEBUG_NAVIGATION
>> 1164 if( fVerbose >= 2 )
1629 { 1165 {
1630 *pNormalCalculated = true; // ComputeS << 1166 G4cout << " G4Navigator::ComputeSafety() relocates-in-volume to point: "
1631 // (fExitin << 1167 << pGlobalpoint << G4endl;
>> 1168 }
>> 1169 #endif
>> 1170 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
>> 1171 G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
>> 1172 G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
>> 1173 G4ThreeVector localPoint = ComputeLocalPoint(pGlobalpoint);
>> 1174
>> 1175 if ( fHistory.GetTopVolumeType()!=kReplica )
>> 1176 {
>> 1177 switch(CharacteriseDaughters(motherLogical))
>> 1178 {
>> 1179 case kNormal:
>> 1180 if ( pVoxelHeader )
>> 1181 {
>> 1182 newSafety=fvoxelNav.ComputeSafety(localPoint,fHistory,pMaxLength);
>> 1183 }
>> 1184 else
>> 1185 {
>> 1186 newSafety=fnormalNav.ComputeSafety(localPoint,fHistory,pMaxLength);
>> 1187 }
>> 1188 break;
>> 1189 case kParameterised:
>> 1190 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
>> 1191 {
>> 1192 newSafety = fparamNav.ComputeSafety(localPoint,fHistory,pMaxLength);
>> 1193 }
>> 1194 else // Regular structure
>> 1195 {
>> 1196 newSafety = fregularNav.ComputeSafety(localPoint,fHistory,pMaxLength);
>> 1197 }
>> 1198 break;
>> 1199 case kReplica:
>> 1200 G4Exception("G4Navigator::ComputeSafety()", "NotApplicable",
>> 1201 FatalException, "Not applicable for replicated volumes.");
>> 1202 break;
>> 1203 }
1632 } 1204 }
1633 else 1205 else
1634 { 1206 {
1635 G4ExceptionDescription message; << 1207 newSafety = freplicaNav.ComputeSafety(pGlobalpoint, localPoint,
1636 message.precision(10); << 1208 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 } 1209 }
1676 } 1210 }
1677 else << 1211 else // if( endpointOnSurface && stayedOnEndpoint )
1678 { 1212 {
1679 localNormal = GetLocalExitNormalAndCheck( <<
1680 *pNormalCalculated = fCalculatedExitNorma <<
1681 <<
1682 #ifdef G4DEBUG_NAVIGATION 1213 #ifdef G4DEBUG_NAVIGATION
1683 usingStored = false; << 1214 if( fVerbose >= 2 )
1684 <<
1685 if( (!validNormal) && !fCalculatedExitNor <<
1686 { 1215 {
1687 G4ExceptionDescription edN; << 1216 G4cout << " G4Navigator::ComputeSafety() finds that point - "
1688 edN << " Calculated = " << fCalculated << 1217 << pGlobalpoint << " - is on surface " << G4endl;
1689 edN << " Entering= " << fEntering << << 1218 if( fEnteredDaughter ) { G4cout << " entered new daughter volume"; }
1690 G4int oldVerbose = this->GetVerboseLeve << 1219 if( fExitedMother ) { G4cout << " and exited previous volume."; }
1691 this->SetVerboseLevel(4); << 1220 G4cout << G4endl;
1692 edN << " State of Navigator: " << G4e << 1221 G4cout << " EndPoint was = " << fStepEndPoint << G4endl;
1693 edN << *this << G4endl; << 1222 }
1694 this->SetVerboseLevel( oldVerbose ); <<
1695 <<
1696 G4Exception("G4Navigator::GetGlobalExit <<
1697 "GeomNav0003", JustWarning, <<
1698 "LocalExitNormalAndCheck() <<
1699 } <<
1700 #endif 1223 #endif
1701 << 1224 newSafety = 0.0;
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 <<
1726 G4Exception("G4Navigator::GetGlobalExi <<
1727 "GeomNav0003",JustWarning, <<
1728 "Value obtained from new l <<
1729 localNormal = localNormal.unit(); // S <<
1730 } <<
1731 globalNormal = fHistory.GetTopTransform( <<
1732 .InverseTransformAxis(loca <<
1733 } 1225 }
1734 1226
1735 #ifdef G4DEBUG_NAVIGATION << 1227 // 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 // 1228 //
1764 fExitNormalGlobalFrame = globalNormal; << 1229 fPreviousSftOrigin = pGlobalpoint;
1765 << 1230 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 1231
1785 G4double distEndpointSq = (pGlobalpoint-fSt << 1232 if (keepState) { RestoreSavedState(); }
1786 G4bool stayedOnEndpoint = distEndpointSq < <<
1787 G4bool endpointOnSurface = fEnteredDaughter <<
1788 1233
1789 G4bool onSurface = endpointOnSurface && sta << 1234 #ifdef G4DEBUG_NAVIGATION
1790 if( ! onSurface ) << 1235 if( fVerbose > 1 )
1791 { 1236 {
1792 safety= fpSafetyCalculator->SafetyInCurre << 1237 G4cout << " ---- Exiting ComputeSafety " << G4endl;
1793 // offload to G4SafetyCalculator - avoids << 1238 if( fVerbose > 2 ) { PrintState(); }
1794 << 1239 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 } 1240 }
>> 1241 G4cout.precision(oldcoutPrec);
>> 1242 #endif
1801 1243
1802 return safety; << 1244 return newSafety;
1803 } 1245 }
1804 1246
1805 // ****************************************** 1247 // ********************************************************************
1806 // CreateTouchableHistoryHandle 1248 // CreateTouchableHistoryHandle
1807 // ****************************************** 1249 // ********************************************************************
1808 // 1250 //
1809 G4TouchableHandle G4Navigator::CreateTouchabl << 1251 G4TouchableHistoryHandle G4Navigator::CreateTouchableHistoryHandle() const
1810 { 1252 {
1811 return G4TouchableHandle( CreateTouchableHi << 1253 return G4TouchableHistoryHandle( CreateTouchableHistory() );
1812 } 1254 }
1813 1255
1814 // ****************************************** 1256 // ********************************************************************
1815 // PrintState 1257 // PrintState
1816 // ****************************************** 1258 // ********************************************************************
1817 // 1259 //
1818 void G4Navigator::PrintState() const 1260 void G4Navigator::PrintState() const
1819 { 1261 {
1820 G4long oldcoutPrec = G4cout.precision(4); << 1262 G4int oldcoutPrec = G4cout.precision(4);
1821 if( fVerbose >= 4 ) << 1263 if( fVerbose == 4 )
1822 { 1264 {
1823 G4cout << "The current state of G4Navigat 1265 G4cout << "The current state of G4Navigator is: " << G4endl;
1824 G4cout << " ValidExitNormal= " << fValid << 1266 G4cout << " ValidExitNormal= " << fValidExitNormal << G4endl
1825 << " ExitNormal = " << fExitN << 1267 << " ExitNormal = " << fExitNormal << G4endl
1826 << " Exiting = " << fExiti << 1268 << " Exiting = " << fExiting << G4endl
1827 << " Entering = " << fEnter << 1269 << " Entering = " << fEntering << G4endl
1828 << " BlockedPhysicalVolume= " ; 1270 << " BlockedPhysicalVolume= " ;
1829 if (fBlockedPhysicalVolume==nullptr) << 1271 if (fBlockedPhysicalVolume==0)
1830 { <<
1831 G4cout << "None"; 1272 G4cout << "None";
1832 } <<
1833 else 1273 else
1834 { <<
1835 G4cout << fBlockedPhysicalVolume->GetNa 1274 G4cout << fBlockedPhysicalVolume->GetName();
1836 } <<
1837 G4cout << G4endl 1275 G4cout << G4endl
1838 << " BlockedReplicaNo = " << << 1276 << " BlockedReplicaNo = " << fBlockedReplicaNo << G4endl
1839 << " LastStepWasZero = " << << 1277 << " LastStepWasZero = " << fLastStepWasZero << G4endl
1840 << G4endl; 1278 << G4endl;
1841 } 1279 }
1842 if( ( 1 < fVerbose) && (fVerbose < 4) ) 1280 if( ( 1 < fVerbose) && (fVerbose < 4) )
1843 { 1281 {
1844 G4cout << G4endl; // Make sure to line up << 1282 G4cout << std::setw(30) << " ExitNormal " << " "
1845 G4cout << std::setw(30) << " ExitNormal " <<
1846 << std::setw( 5) << " Valid " 1283 << std::setw( 5) << " Valid " << " "
1847 << std::setw( 9) << " Exiting " 1284 << std::setw( 9) << " Exiting " << " "
1848 << std::setw( 9) << " Entering" 1285 << std::setw( 9) << " Entering" << " "
1849 << std::setw(15) << " Blocked:Volu 1286 << std::setw(15) << " Blocked:Volume " << " "
1850 << std::setw( 9) << " ReplicaNo" 1287 << std::setw( 9) << " ReplicaNo" << " "
1851 << std::setw( 8) << " LastStepZero 1288 << std::setw( 8) << " LastStepZero " << " "
1852 << G4endl; 1289 << G4endl;
1853 G4cout << "( " << std::setw(7) << fExitNo 1290 G4cout << "( " << std::setw(7) << fExitNormal.x()
1854 << ", " << std::setw(7) << fExitNo 1291 << ", " << std::setw(7) << fExitNormal.y()
1855 << ", " << std::setw(7) << fExitNo 1292 << ", " << std::setw(7) << fExitNormal.z() << " ) "
1856 << std::setw( 5) << fValidExitNor 1293 << std::setw( 5) << fValidExitNormal << " "
1857 << std::setw( 9) << fExiting 1294 << std::setw( 9) << fExiting << " "
1858 << std::setw( 9) << fEntering 1295 << std::setw( 9) << fEntering << " ";
1859 if ( fBlockedPhysicalVolume == nullptr ) << 1296 if ( fBlockedPhysicalVolume==0 )
1860 { G4cout << std::setw(15) << "None"; } << 1297 G4cout << std::setw(15) << "None";
1861 else 1298 else
1862 { G4cout << std::setw(15)<< fBlockedPhysi << 1299 G4cout << std::setw(15)<< fBlockedPhysicalVolume->GetName();
1863 G4cout << std::setw( 9) << fBlockedRepli << 1300 G4cout << std::setw( 9) << fBlockedReplicaNo << " "
1864 << std::setw( 8) << fLastStepWasZ << 1301 << std::setw( 8) << fLastStepWasZero << " "
1865 << G4endl; << 1302 << G4endl;
1866 } 1303 }
1867 if( fVerbose > 2 ) 1304 if( fVerbose > 2 )
1868 { 1305 {
1869 G4cout.precision(8); 1306 G4cout.precision(8);
1870 G4cout << " Current Localpoint = " << fLa 1307 G4cout << " Current Localpoint = " << fLastLocatedPointLocal << G4endl;
1871 G4cout << " PreviousSftOrigin = " << fPr 1308 G4cout << " PreviousSftOrigin = " << fPreviousSftOrigin << G4endl;
1872 G4cout << " PreviousSafety = " << fPr 1309 G4cout << " PreviousSafety = " << fPreviousSafety << G4endl;
1873 } 1310 }
1874 G4cout.precision(oldcoutPrec); 1311 G4cout.precision(oldcoutPrec);
1875 } 1312 }
1876 1313
1877 // ****************************************** 1314 // ********************************************************************
1878 // ComputeStepLog 1315 // ComputeStepLog
1879 // ****************************************** 1316 // ********************************************************************
1880 // 1317 //
1881 void G4Navigator::ComputeStepLog(const G4Thre 1318 void G4Navigator::ComputeStepLog(const G4ThreeVector& pGlobalpoint,
1882 G4doub 1319 G4double moveLenSq) const
1883 { 1320 {
1884 // The following checks only make sense if 1321 // The following checks only make sense if the move is larger
1885 // than the tolerance. 1322 // than the tolerance.
1886 1323
1887 const G4double fAccuracyForWarning = kCar << 1324 static const G4double fAccuracyForWarning = kCarTolerance,
1888 fAccuracyForException = 1000 << 1325 fAccuracyForException = 1000*kCarTolerance;
1889 1326
1890 G4ThreeVector OriginalGlobalpoint = fHistor << 1327 G4ThreeVector OriginalGlobalpoint = fHistory.GetTopTransform().Inverse().
1891 InverseTransfo << 1328 TransformPoint(fLastLocatedPointLocal);
1892 1329
1893 G4double shiftOriginSafSq = (fPreviousSftOr 1330 G4double shiftOriginSafSq = (fPreviousSftOrigin-pGlobalpoint).mag2();
1894 1331
1895 // Check that the starting point of this st 1332 // Check that the starting point of this step is
1896 // within the isotropic safety sphere of th 1333 // within the isotropic safety sphere of the last point
1897 // to a accuracy/precision given by fAccur 1334 // to a accuracy/precision given by fAccuracyForWarning.
1898 // If so give warning. 1335 // If so give warning.
1899 // If it fails by more than fAccuracyForE 1336 // If it fails by more than fAccuracyForException exit with error.
1900 // 1337 //
1901 if( shiftOriginSafSq >= sqr(fPreviousSafety 1338 if( shiftOriginSafSq >= sqr(fPreviousSafety) )
1902 { 1339 {
1903 G4double shiftOrigin = std::sqrt(shiftOri 1340 G4double shiftOrigin = std::sqrt(shiftOriginSafSq);
1904 G4double diffShiftSaf = shiftOrigin - fPr 1341 G4double diffShiftSaf = shiftOrigin - fPreviousSafety;
1905 1342
1906 if( diffShiftSaf > fAccuracyForWarning ) 1343 if( diffShiftSaf > fAccuracyForWarning )
1907 { 1344 {
1908 G4long oldcoutPrec = G4cout.precision(8 << 1345 G4int oldcoutPrec= G4cout.precision(8);
1909 G4long oldcerrPrec = G4cerr.precision(1 << 1346 G4int oldcerrPrec= G4cerr.precision(10);
1910 std::ostringstream message, suggestion; << 1347 G4Exception("G4Navigator::ComputeStep()",
1911 message << "Accuracy error or slightly << 1348 "UnexpectedPositionShift", JustWarning,
1912 << G4endl << 1349 "Accuracy error or slightly inaccurate position shift.");
1913 << " The Step's starting po << 1350 G4cerr << " The Step's starting point has moved "
1914 << std::sqrt(moveLenSq)/mm << " << 1351 << std::sqrt(moveLenSq)/mm << " mm " << G4endl
1915 << " since the last call to << 1352 << " since the last call to a Locate method." << G4endl;
1916 << " This has resulted in m << 1353 G4cerr << " This has resulted in moving "
1917 << shiftOrigin/mm << " mm " << 1354 << shiftOrigin/mm << " mm "
1918 << " from the last point at whi << 1355 << " from the last point at which the safety "
1919 << " was calculated " << G4 << 1356 << " was calculated " << G4endl;
1920 << " which is more than the << 1357 G4cerr << " which is more than the computed safety= "
1921 << fPreviousSafety/mm << " mm << 1358 << fPreviousSafety/mm << " mm at that point." << G4endl;
1922 << " This difference is " << 1359 G4cerr << " This difference is "
1923 << diffShiftSaf/mm << " mm." << << 1360 << diffShiftSaf/mm << " mm." << G4endl
1924 << " The tolerated accuracy << 1361 << " The tolerated accuracy is "
1925 << fAccuracyForException/mm << << 1362 << fAccuracyForException/mm << " mm." << G4endl;
1926 1363
1927 suggestion << " "; << 1364 static G4int warnNow = 0;
1928 static G4ThreadLocal G4int warnNow = 0; <<
1929 if( ((++warnNow % 100) == 1) ) 1365 if( ((++warnNow % 100) == 1) )
1930 { 1366 {
1931 message << G4endl << 1367 G4cerr << " This problem can be due to either " << G4endl;
1932 << " This problem can be due << 1368 G4cerr << " - a process that has proposed a displacement"
1933 << " - a process that has p << 1369 << " larger than the current safety , or" << G4endl;
1934 << " larger than the current s << 1370 G4cerr << " - inaccuracy in the computation of the safety"
1935 << " - inaccuracy in the co << 1371 << G4endl;
1936 suggestion << "We suggest that you " << 1372 G4cerr << " We suggest that you " << G4endl
1937 << " - find i) what part << 1373 << " - find i) what particle is being tracked, and "
1938 << " ii) through what part << 1374 << " ii) through what part of your geometry " << G4endl
1939 << " for example by r << 1375 << " for example by re-running this event with "
1940 << G4endl << 1376 << G4endl
1941 << " /tracking/ver << 1377 << " /tracking/verbose 1 " << G4endl
1942 << " - check which proc << 1378 << " - check which processes you declare for"
1943 << " this particle (and lo << 1379 << " this particle (and look at non-standard ones)"
1944 << G4endl << 1380 << G4endl
1945 << " - in case, create a << 1381 << " - in case, create a detailed logfile"
1946 << " of this event using:" << 1382 << " of this event using:" << G4endl
1947 << " /tracking/ver << 1383 << " /tracking/verbose 6 "
>> 1384 << G4endl;
1948 } 1385 }
1949 G4Exception("G4Navigator::ComputeStep() <<
1950 "GeomNav1002", JustWarning, <<
1951 message, G4String(suggestio <<
1952 G4cout.precision(oldcoutPrec); 1386 G4cout.precision(oldcoutPrec);
1953 G4cerr.precision(oldcerrPrec); 1387 G4cerr.precision(oldcerrPrec);
1954 } 1388 }
1955 #ifdef G4DEBUG_NAVIGATION 1389 #ifdef G4DEBUG_NAVIGATION
1956 else 1390 else
1957 { 1391 {
1958 G4cerr << "WARNING - G4Navigator::Compu 1392 G4cerr << "WARNING - G4Navigator::ComputeStep()" << G4endl
1959 << " The Step's startin 1393 << " The Step's starting point has moved "
1960 << std::sqrt(moveLenSq) << "," < 1394 << std::sqrt(moveLenSq) << "," << G4endl
1961 << " which has taken it 1395 << " which has taken it to the limit of"
1962 << " the current safety. " << G4 1396 << " the current safety. " << G4endl;
1963 } 1397 }
1964 #endif 1398 #endif
1965 } 1399 }
1966 G4double safetyPlus = fPreviousSafety + fAc 1400 G4double safetyPlus = fPreviousSafety + fAccuracyForException;
1967 if ( shiftOriginSafSq > sqr(safetyPlus) ) 1401 if ( shiftOriginSafSq > sqr(safetyPlus) )
1968 { 1402 {
1969 std::ostringstream message; << 1403 G4cerr << "ERROR - G4Navigator::ComputeStep()" << G4endl
1970 message << "May lead to a crash or unreli << 1404 << " Position has shifted considerably without"
1971 << " Position has shifted << 1405 << " notifying the navigator !" << G4endl
1972 << " notifying the navigator !" < << 1406 << " Tolerated safety: " << safetyPlus << G4endl
1973 << " Tolerated safety: " < << 1407 << " Computed shift : " << shiftOriginSafSq << G4endl;
1974 << " Computed shift : " < << 1408 G4Exception("G4Navigator::ComputeStep()",
1975 G4Exception("G4Navigator::ComputeStep()", << 1409 "SignificantPositionShift", JustWarning,
1976 JustWarning, message); << 1410 "May lead to a crash or unreliable results.");
1977 } 1411 }
1978 } 1412 }
1979 1413
1980 // ****************************************** 1414 // ********************************************************************
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 << 1415 // Operator <<
2009 // ****************************************** 1416 // ********************************************************************
2010 // 1417 //
2011 std::ostream& operator << (std::ostream &os,c 1418 std::ostream& operator << (std::ostream &os,const G4Navigator &n)
2012 { 1419 {
2013 // Old version did only the following: << 1420 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; 1421 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 } 1422 }
2116 1423