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