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