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