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