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