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