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