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