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