| 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 <<
986 Step = fvoxelNav.ComputeStep(fLastLo 985 Step = fvoxelNav.ComputeStep(fLastLocatedPointLocal,
987 localDi 986 localDirection,
988 pCurren 987 pCurrentProposedStepLength,
989 pNewSaf 988 pNewSafety,
990 fHistor 989 fHistory,
991 fValidE 990 fValidExitNormal,
992 fExitNo 991 fExitNormal,
993 fExitin 992 fExiting,
994 fEnteri 993 fEntering,
995 &fBlock 994 &fBlockedPhysicalVolume,
996 fBlocke 995 fBlockedReplicaNo);
997 996
998 } 997 }
999 else 998 else
1000 { 999 {
1001 if( motherPhysical->GetRegularStruc 1000 if( motherPhysical->GetRegularStructureId() == 0 )
1002 { 1001 {
1003 Step = fnormalNav.ComputeStep(fLa 1002 Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
1004 loc 1003 localDirection,
1005 pCu 1004 pCurrentProposedStepLength,
1006 pNe 1005 pNewSafety,
1007 fHi 1006 fHistory,
1008 fVa 1007 fValidExitNormal,
1009 fEx 1008 fExitNormal,
1010 fEx 1009 fExiting,
1011 fEn 1010 fEntering,
1012 &fB 1011 &fBlockedPhysicalVolume,
1013 fBl 1012 fBlockedReplicaNo);
1014 } 1013 }
1015 else // Regular (non-voxelised) st 1014 else // Regular (non-voxelised) structure
1016 { 1015 {
1017 LocateGlobalPointAndSetup( pGloba 1016 LocateGlobalPointAndSetup( pGlobalpoint, &pDirection, true, true );
1018 fLastTriedStepComputation= true; 1017 fLastTriedStepComputation= true; // Ensure that this is set again !!
1019 // 1018 //
1020 // if physical process limits the 1019 // if physical process limits the step, the voxel will not be the
1021 // one given by ComputeStepSkippi 1020 // one given by ComputeStepSkippingEqualMaterials() and the local
1022 // point will be wrongly calculat 1021 // point will be wrongly calculated.
1023 1022
1024 // There is a problem: when msc l 1023 // There is a problem: when msc limits the step and the point is
1025 // assigned wrongly to phantom in 1024 // assigned wrongly to phantom in previous step (while it is out
1026 // of the container volume). Then 1025 // of the container volume). Then LocateGlobalPointAndSetup() has
1027 // reset the history topvolume to 1026 // reset the history topvolume to world.
1028 // 1027 //
1029 if(fHistory.GetTopVolume()->GetRe 1028 if(fHistory.GetTopVolume()->GetRegularStructureId() == 0 )
1030 { 1029 {
1031 G4Exception("G4ITNavigator2::Co 1030 G4Exception("G4ITNavigator2::ComputeStep()",
1032 "GeomNav1001", Just 1031 "GeomNav1001", JustWarning,
1033 "Point is relocated in voxels 1032 "Point is relocated in voxels, while it should be outside!");
1034 Step = fnormalNav.ComputeStep(f 1033 Step = fnormalNav.ComputeStep(fLastLocatedPointLocal,
1035 l 1034 localDirection,
1036 p 1035 pCurrentProposedStepLength,
1037 p 1036 pNewSafety,
1038 f 1037 fHistory,
1039 f 1038 fValidExitNormal,
1040 f 1039 fExitNormal,
1041 f 1040 fExiting,
1042 f 1041 fEntering,
1043 & 1042 &fBlockedPhysicalVolume,
1044 f 1043 fBlockedReplicaNo);
1045 } 1044 }
1046 else 1045 else
1047 { 1046 {
1048 Step = fregularNav. 1047 Step = fregularNav.
1049 ComputeStepSkippingEqualMa 1048 ComputeStepSkippingEqualMaterials(fLastLocatedPointLocal,
1050 1049 localDirection,
1051 1050 pCurrentProposedStepLength,
1052 1051 pNewSafety,
1053 1052 fHistory,
1054 1053 fValidExitNormal,
1055 1054 fExitNormal,
1056 1055 fExiting,
1057 1056 fEntering,
1058 1057 &fBlockedPhysicalVolume,
1059 1058 fBlockedReplicaNo,
1060 1059 motherPhysical);
1061 } 1060 }
1062 } 1061 }
1063 } 1062 }
1064 break; 1063 break;
1065 case kParameterised: 1064 case kParameterised:
1066 if( GetDaughtersRegularStructureId(mo 1065 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
1067 { 1066 {
1068 Step = fparamNav.ComputeStep(fLastL 1067 Step = fparamNav.ComputeStep(fLastLocatedPointLocal,
1069 localD 1068 localDirection,
1070 pCurre 1069 pCurrentProposedStepLength,
1071 pNewSa 1070 pNewSafety,
1072 fHisto 1071 fHistory,
1073 fValid 1072 fValidExitNormal,
1074 fExitN 1073 fExitNormal,
1075 fExiti 1074 fExiting,
1076 fEnter 1075 fEntering,
1077 &fBloc 1076 &fBlockedPhysicalVolume,
1078 fBlock 1077 fBlockedReplicaNo);
1079 } 1078 }
1080 else // Regular structure 1079 else // Regular structure
1081 { 1080 {
1082 Step = fregularNav.ComputeStep(fLas 1081 Step = fregularNav.ComputeStep(fLastLocatedPointLocal,
1083 loca 1082 localDirection,
1084 pCur 1083 pCurrentProposedStepLength,
1085 pNew 1084 pNewSafety,
1086 fHis 1085 fHistory,
1087 fVal 1086 fValidExitNormal,
1088 fExi 1087 fExitNormal,
1089 fExi 1088 fExiting,
1090 fEnt 1089 fEntering,
1091 &fBl 1090 &fBlockedPhysicalVolume,
1092 fBlo 1091 fBlockedReplicaNo);
1093 } 1092 }
1094 break; 1093 break;
1095 case kReplica: 1094 case kReplica:
1096 G4Exception("G4ITNavigator2::ComputeS 1095 G4Exception("G4ITNavigator2::ComputeStep()", "GeomNav0001",
1097 FatalException, "Not appl 1096 FatalException, "Not applicable for replicated volumes.");
1098 break; 1097 break;
1099 case kExternal: <<
1100 G4Exception("G4ITNavigator2::ComputeS <<
1101 FatalException, "Not appl <<
1102 break; <<
1103 } 1098 }
1104 } 1099 }
1105 else 1100 else
1106 { 1101 {
1107 // In the case of a replica, it must hand 1102 // In the case of a replica, it must handle the exiting
1108 // edge/corner problem by itself 1103 // edge/corner problem by itself
1109 // 1104 //
1110 G4bool exitingReplica = fExitedMother; 1105 G4bool exitingReplica = fExitedMother;
1111 G4bool calculatedExitNormal; 1106 G4bool calculatedExitNormal;
1112 Step = freplicaNav.ComputeStep(pGlobalpoi 1107 Step = freplicaNav.ComputeStep(pGlobalpoint,
1113 pDirection 1108 pDirection,
1114 fLastLocat 1109 fLastLocatedPointLocal,
1115 localDirec 1110 localDirection,
1116 pCurrentPr 1111 pCurrentProposedStepLength,
1117 pNewSafety 1112 pNewSafety,
1118 fHistory, 1113 fHistory,
1119 fValidExit 1114 fValidExitNormal,
1120 calculated 1115 calculatedExitNormal,
1121 fExitNorma 1116 fExitNormal,
1122 exitingRep 1117 exitingReplica,
1123 fEntering, 1118 fEntering,
1124 &fBlockedP 1119 &fBlockedPhysicalVolume,
1125 fBlockedRe 1120 fBlockedReplicaNo);
1126 fExiting= exitingReplica; 1121 fExiting= exitingReplica;
1127 fCalculatedExitNormal= calculatedExitNorm 1122 fCalculatedExitNormal= calculatedExitNormal;
1128 } 1123 }
1129 1124
1130 1125
1131 // G4cout << " !!!! Step = " << Step << G4en 1126 // G4cout << " !!!! Step = " << Step << G4endl;
1132 1127
1133 // Remember last safety origin & value. 1128 // Remember last safety origin & value.
1134 // 1129 //
1135 fPreviousSftOrigin = pGlobalpoint; 1130 fPreviousSftOrigin = pGlobalpoint;
1136 fPreviousSafety = pNewSafety; 1131 fPreviousSafety = pNewSafety;
1137 1132
1138 // Count zero steps - one can occur due to 1133 // Count zero steps - one can occur due to changing momentum at a boundary
1139 // - one, two (or a few) c 1134 // - one, two (or a few) can occur at common edges between
1140 // volumes 1135 // volumes
1141 // - more than two is like 1136 // - more than two is likely a problem in the geometry
1142 // description or the Na 1137 // description or the Navigation
1143 1138
1144 // Rule of thumb: likely at an Edge if two 1139 // Rule of thumb: likely at an Edge if two consecutive steps are zero,
1145 // because at least two cand 1140 // because at least two candidate volumes must have been
1146 // checked 1141 // checked
1147 // 1142 //
1148 fLocatedOnEdge = fLastStepWasZero && (Ste 1143 fLocatedOnEdge = fLastStepWasZero && (Step==0.0);
1149 fLastStepWasZero = (Step==0.0); 1144 fLastStepWasZero = (Step==0.0);
1150 if (fPushed) { fPushed = fLastStepWasZero; 1145 if (fPushed) { fPushed = fLastStepWasZero; }
1151 1146
1152 // Handle large number of consecutive zero 1147 // Handle large number of consecutive zero steps
1153 // 1148 //
1154 if ( fLastStepWasZero ) 1149 if ( fLastStepWasZero )
1155 { 1150 {
1156 fNumberZeroSteps++; 1151 fNumberZeroSteps++;
1157 #ifdef G4DEBUG_NAVIGATION 1152 #ifdef G4DEBUG_NAVIGATION
1158 if( fNumberZeroSteps > 1 ) 1153 if( fNumberZeroSteps > 1 )
1159 { 1154 {
1160 G4cout << "G4ITNavigator2::ComputeStep 1155 G4cout << "G4ITNavigator2::ComputeStep(): another zero step, # "
1161 << fNumberZeroSteps 1156 << fNumberZeroSteps
1162 << " at " << pGlobalpoint 1157 << " at " << pGlobalpoint
1163 << " in volume " << motherPhysi 1158 << " in volume " << motherPhysical->GetName()
>> 1159 << " nav-comp-step calls # " << sNavCScalls
1164 << G4endl; 1160 << G4endl;
1165 } 1161 }
1166 #endif 1162 #endif
1167 if( fNumberZeroSteps > fActionThreshold_N 1163 if( fNumberZeroSteps > fActionThreshold_NoZeroSteps-1 )
1168 { 1164 {
1169 // Act to recover this stuck track. Pu 1165 // Act to recover this stuck track. Pushing it along direction
1170 // 1166 //
1171 Step += 100*kCarTolerance; 1167 Step += 100*kCarTolerance;
1172 #ifdef G4VERBOSE 1168 #ifdef G4VERBOSE
1173 if ((!fPushed) && (fWarnPush)) 1169 if ((!fPushed) && (fWarnPush))
1174 { 1170 {
1175 std::ostringstream message; 1171 std::ostringstream message;
1176 message << "Track stuck or not movin 1172 message << "Track stuck or not moving." << G4endl
1177 << " Track stuck, n 1173 << " Track stuck, not moving for "
1178 << fNumberZeroSteps << " ste 1174 << fNumberZeroSteps << " steps" << G4endl
1179 << " in volume -" < 1175 << " in volume -" << motherPhysical->GetName()
1180 << "- at point " << pGlobalp 1176 << "- at point " << pGlobalpoint << G4endl
1181 << " direction: " < 1177 << " direction: " << pDirection << "." << G4endl
1182 << " Potential geom 1178 << " Potential geometry or navigation problem !"
1183 << G4endl 1179 << G4endl
1184 << " Trying pushing 1180 << " Trying pushing it of " << Step << " mm ...";
1185 G4Exception("G4ITNavigator2::Compute 1181 G4Exception("G4ITNavigator2::ComputeStep()", "GeomNav1002",
1186 JustWarning, message, "P 1182 JustWarning, message, "Potential overlap in geometry!");
1187 } 1183 }
1188 #endif 1184 #endif
1189 fPushed = true; 1185 fPushed = true;
1190 } 1186 }
1191 if( fNumberZeroSteps > fAbandonThreshold_ 1187 if( fNumberZeroSteps > fAbandonThreshold_NoZeroSteps-1 )
1192 { 1188 {
1193 // Must kill this stuck track 1189 // Must kill this stuck track
1194 // 1190 //
1195 std::ostringstream message; 1191 std::ostringstream message;
1196 message << "Stuck Track: potential geom 1192 message << "Stuck Track: potential geometry or navigation problem."
1197 << G4endl 1193 << G4endl
1198 << " Track stuck, not mo 1194 << " Track stuck, not moving for "
1199 << fNumberZeroSteps << " steps" 1195 << fNumberZeroSteps << " steps" << G4endl
1200 << " in volume -" << mot 1196 << " in volume -" << motherPhysical->GetName()
1201 << "- at point " << pGlobalpoin 1197 << "- at point " << pGlobalpoint << G4endl
1202 << " direction: " << pDi 1198 << " direction: " << pDirection << ".";
1203 motherPhysical->CheckOverlaps(5000, 0.0 << 1199 motherPhysical->CheckOverlaps(5000, false);
1204 G4Exception("G4ITNavigator2::ComputeSte 1200 G4Exception("G4ITNavigator2::ComputeStep()", "GeomNav0003",
1205 EventMustBeAborted, message 1201 EventMustBeAborted, message);
1206 } 1202 }
1207 } 1203 }
1208 else 1204 else
1209 { 1205 {
1210 if (!fPushed) fNumberZeroSteps = 0; 1206 if (!fPushed) fNumberZeroSteps = 0;
1211 } 1207 }
1212 1208
1213 fEnteredDaughter = fEntering; // I expect 1209 fEnteredDaughter = fEntering; // I expect to enter a volume in this Step
1214 fExitedMother = fExiting; 1210 fExitedMother = fExiting;
1215 1211
1216 fStepEndPoint = pGlobalpoint 1212 fStepEndPoint = pGlobalpoint
1217 + std::min(Step,pCurrentPropo 1213 + std::min(Step,pCurrentProposedStepLength) * pDirection;
1218 fLastStepEndPointLocal = fLastLocatedPointL 1214 fLastStepEndPointLocal = fLastLocatedPointLocal + Step * localDirection;
1219 1215
1220 if( fExiting ) 1216 if( fExiting )
1221 { 1217 {
1222 #ifdef G4DEBUG_NAVIGATION 1218 #ifdef G4DEBUG_NAVIGATION
1223 if( fVerbose > 2 ) 1219 if( fVerbose > 2 )
1224 { 1220 {
1225 G4cout << " At G4Nav CompStep End - if( 1221 G4cout << " At G4Nav CompStep End - if(exiting) - fExiting= " << fExiting
1226 << " fValidExitNormal = " << fVa 1222 << " fValidExitNormal = " << fValidExitNormal << G4endl;
1227 G4cout << " fExitNormal= " << fExitNorm 1223 G4cout << " fExitNormal= " << fExitNormal << G4endl;
1228 } 1224 }
1229 #endif 1225 #endif
1230 1226
1231 if(fValidExitNormal || fCalculatedExitNor 1227 if(fValidExitNormal || fCalculatedExitNormal)
1232 { 1228 {
1233 if ( fHistory.GetTopVolumeType()!=kRep 1229 if ( fHistory.GetTopVolumeType()!=kReplica )
1234 { 1230 {
1235 // Convention: fExitNormal is in the 1231 // Convention: fExitNormal is in the 'grand-mother' coordinate system
1236 // 1232 //
1237 fGrandMotherExitNormal= fExitNormal; 1233 fGrandMotherExitNormal= fExitNormal;
1238 fCalculatedExitNormal= true; 1234 fCalculatedExitNormal= true;
1239 } 1235 }
1240 else 1236 else
1241 { 1237 {
1242 fGrandMotherExitNormal = fExitNormal; 1238 fGrandMotherExitNormal = fExitNormal;
1243 } 1239 }
1244 } 1240 }
1245 else 1241 else
1246 { 1242 {
1247 // We must calculate the normal anyway 1243 // We must calculate the normal anyway (in order to have it if requested)
1248 // 1244 //
1249 G4ThreeVector finalLocalPoint = 1245 G4ThreeVector finalLocalPoint =
1250 fLastLocatedPointLocal + localDir 1246 fLastLocatedPointLocal + localDirection*Step;
1251 1247
1252 if ( fHistory.GetTopVolumeType()!=kRep 1248 if ( fHistory.GetTopVolumeType()!=kReplica )
1253 { 1249 {
1254 // Find normal in the 'mother' coordi 1250 // Find normal in the 'mother' coordinate system
1255 // 1251 //
1256 G4ThreeVector exitNormalMotherFrame= 1252 G4ThreeVector exitNormalMotherFrame=
1257 motherLogical->GetSolid()->Surface 1253 motherLogical->GetSolid()->SurfaceNormal(finalLocalPoint);
1258 1254
1259 // Transform it to the 'grand-mother' 1255 // Transform it to the 'grand-mother' coordinate system
1260 // 1256 //
1261 const G4RotationMatrix* mRot = mother 1257 const G4RotationMatrix* mRot = motherPhysical->GetRotation();
1262 if( mRot != nullptr ) << 1258 if( mRot )
1263 { 1259 {
1264 fChangedGrandMotherRefFrame= true; 1260 fChangedGrandMotherRefFrame= true;
1265 fGrandMotherExitNormal = (*mRot).in 1261 fGrandMotherExitNormal = (*mRot).inverse() * exitNormalMotherFrame;
1266 } 1262 }
1267 else 1263 else
1268 { 1264 {
1269 fGrandMotherExitNormal = exitNormal 1265 fGrandMotherExitNormal = exitNormalMotherFrame;
1270 } 1266 }
1271 1267
1272 // Do not set fValidExitNormal -- thi 1268 // Do not set fValidExitNormal -- this signifies
1273 // that the solid is convex! 1269 // that the solid is convex!
1274 // 1270 //
1275 fCalculatedExitNormal= true; 1271 fCalculatedExitNormal= true;
1276 } 1272 }
1277 else 1273 else
1278 { 1274 {
1279 fCalculatedExitNormal = false; 1275 fCalculatedExitNormal = false;
1280 // 1276 //
1281 // Nothing can be done at this stage 1277 // Nothing can be done at this stage currently - to solve this
1282 // Replica Navigation must have calcu 1278 // Replica Navigation must have calculated the normal for this case
1283 // already. 1279 // already.
1284 // Cases: mother is not convex, and e 1280 // Cases: mother is not convex, and exit is at previous replica level
1285 1281
1286 #ifdef G4DEBUG_NAVIGATION 1282 #ifdef G4DEBUG_NAVIGATION
1287 G4ExceptionDescription desc; 1283 G4ExceptionDescription desc;
1288 1284
1289 desc << "Problem in ComputeStep: Rep 1285 desc << "Problem in ComputeStep: Replica Navigation did not provide"
1290 << " valid exit Normal. " << G4e 1286 << " valid exit Normal. " << G4endl;
1291 desc << " Do not know how calculate i 1287 desc << " Do not know how calculate it in this case." << G4endl;
1292 desc << " Location = " << finalLo 1288 desc << " Location = " << finalLocalPoint << G4endl;
1293 desc << " Volume name = " << motherP 1289 desc << " Volume name = " << motherPhysical->GetName()
1294 << " copy/replica No = " << mot 1290 << " copy/replica No = " << motherPhysical->GetCopyNo() << G4endl;
1295 G4Exception("G4ITNavigator2::ComputeS 1291 G4Exception("G4ITNavigator2::ComputeStep()", "GeomNav0003",
1296 JustWarning, desc, "Norma 1292 JustWarning, desc, "Normal not available for exiting.");
1297 #endif 1293 #endif
1298 } 1294 }
1299 } 1295 }
1300 1296
1301 // Now transform it to the global referen 1297 // Now transform it to the global reference frame !!
1302 // 1298 //
1303 if( fValidExitNormal || fCalculatedExitNo 1299 if( fValidExitNormal || fCalculatedExitNormal )
1304 { 1300 {
1305 auto depth= (G4int)fHistory.GetDepth() << 1301 G4int depth= fHistory.GetDepth();
1306 if( depth > 0 ) 1302 if( depth > 0 )
1307 { 1303 {
1308 G4AffineTransform GrandMotherToGlobal 1304 G4AffineTransform GrandMotherToGlobalTransf =
1309 fHistory.GetTransform(depth-1).Inve 1305 fHistory.GetTransform(depth-1).Inverse();
1310 fExitNormalGlobalFrame = 1306 fExitNormalGlobalFrame =
1311 GrandMotherToGlobalTransf.Transform 1307 GrandMotherToGlobalTransf.TransformAxis( fGrandMotherExitNormal );
1312 } 1308 }
1313 else 1309 else
1314 { 1310 {
1315 fExitNormalGlobalFrame= fGrandMotherE 1311 fExitNormalGlobalFrame= fGrandMotherExitNormal;
1316 } 1312 }
1317 } 1313 }
1318 else 1314 else
1319 { 1315 {
1320 fExitNormalGlobalFrame= G4ThreeVector( 1316 fExitNormalGlobalFrame= G4ThreeVector( 0., 0., 0.);
1321 } 1317 }
1322 } 1318 }
1323 1319
1324 if( (Step == pCurrentProposedStepLength) && 1320 if( (Step == pCurrentProposedStepLength) && (!fExiting) && (!fEntering) )
1325 { 1321 {
1326 // This if Step is not really limited by 1322 // This if Step is not really limited by the geometry.
1327 // The Navigator is obliged to return "in 1323 // The Navigator is obliged to return "infinity"
1328 // 1324 //
1329 Step = kInfinity; 1325 Step = kInfinity;
1330 } 1326 }
1331 1327
1332 #ifdef G4VERBOSE 1328 #ifdef G4VERBOSE
1333 if( fVerbose > 1 ) 1329 if( fVerbose > 1 )
1334 { 1330 {
1335 if( fVerbose >= 4 ) 1331 if( fVerbose >= 4 )
1336 { 1332 {
1337 G4cout << " ----- Upon exiting :" << 1333 G4cout << " ----- Upon exiting :" << G4endl;
1338 PrintState(); 1334 PrintState();
1339 } 1335 }
1340 G4cout << " Returned step= " << Step; 1336 G4cout << " Returned step= " << Step;
1341 if( fVerbose > 5 ) G4cout << G4endl; 1337 if( fVerbose > 5 ) G4cout << G4endl;
1342 if( Step == kInfinity ) 1338 if( Step == kInfinity )
1343 { 1339 {
1344 G4cout << " Requested step= " << pCurr 1340 G4cout << " Requested step= " << pCurrentProposedStepLength ;
1345 if( fVerbose > 5) G4cout << G4endl; 1341 if( fVerbose > 5) G4cout << G4endl;
1346 } 1342 }
1347 G4cout << " Safety = " << pNewSafety << 1343 G4cout << " Safety = " << pNewSafety << G4endl;
1348 } 1344 }
1349 #endif 1345 #endif
1350 1346
1351 return Step; 1347 return Step;
1352 } 1348 }
1353 1349
1354 // ****************************************** 1350 // ********************************************************************
1355 // CheckNextStep 1351 // CheckNextStep
1356 // 1352 //
1357 // Compute the step without altering the navi 1353 // Compute the step without altering the navigator state
1358 // ****************************************** 1354 // ********************************************************************
1359 // 1355 //
1360 G4double G4ITNavigator2::CheckNextStep( const 1356 G4double G4ITNavigator2::CheckNextStep( const G4ThreeVector& pGlobalpoint,
1361 const 1357 const G4ThreeVector& pDirection,
1362 const 1358 const G4double pCurrentProposedStepLength,
1363 1359 G4double& pNewSafety)
1364 { 1360 {
1365 CheckNavigatorStateIsValid(); 1361 CheckNavigatorStateIsValid();
1366 G4double step; 1362 G4double step;
1367 1363
1368 // Save the state, for this parasitic call 1364 // Save the state, for this parasitic call
1369 // 1365 //
1370 //SetSavedState(); 1366 //SetSavedState();
1371 // G4SaveNavigatorState savedState(fpNavigat 1367 // G4SaveNavigatorState savedState(fpNavigatorState);
1372 G4NavigatorState savedState(*fpNavigatorSta 1368 G4NavigatorState savedState(*fpNavigatorState);
1373 1369
1374 step = ComputeStep ( pGlobalpoint, 1370 step = ComputeStep ( pGlobalpoint,
1375 pDirection, 1371 pDirection,
1376 pCurrentProposedStepLe 1372 pCurrentProposedStepLength,
1377 pNewSafety ); 1373 pNewSafety );
1378 1374
1379 // If a parasitic call, then attempt to res 1375 // If a parasitic call, then attempt to restore the key parts of the state
1380 // 1376 //
1381 *fpNavigatorState = savedState; 1377 *fpNavigatorState = savedState;
1382 //RestoreSavedState(); 1378 //RestoreSavedState();
1383 // NOTE: the state of the current subnaviga 1379 // NOTE: the state of the current subnavigator is NOT restored.
1384 // ***> TODO: restore subnavigator state 1380 // ***> TODO: restore subnavigator state
1385 // if( last_located) Need 1381 // if( last_located) Need Position of last location
1386 // if( last_computed step) Need 1382 // if( last_computed step) Need Endposition of last step
1387 1383
1388 1384
1389 return step; 1385 return step;
1390 } 1386 }
1391 1387
1392 // ****************************************** 1388 // ********************************************************************
1393 // ResetState 1389 // ResetState
1394 // 1390 //
1395 // Resets stack and minimum of navigator stat 1391 // Resets stack and minimum of navigator state `machine'
1396 // ****************************************** 1392 // ********************************************************************
1397 // 1393 //
1398 void G4ITNavigator2::ResetState() 1394 void G4ITNavigator2::ResetState()
1399 { 1395 {
1400 fWasLimitedByGeometry = false; 1396 fWasLimitedByGeometry = false;
1401 fEntering = false; 1397 fEntering = false;
1402 fExiting = false; 1398 fExiting = false;
1403 fLocatedOnEdge = false; 1399 fLocatedOnEdge = false;
1404 fLastStepWasZero = false; 1400 fLastStepWasZero = false;
1405 fEnteredDaughter = false; 1401 fEnteredDaughter = false;
1406 fExitedMother = false; 1402 fExitedMother = false;
1407 fPushed = false; 1403 fPushed = false;
1408 1404
1409 fValidExitNormal = false; 1405 fValidExitNormal = false;
1410 fChangedGrandMotherRefFrame= false; 1406 fChangedGrandMotherRefFrame= false;
1411 fCalculatedExitNormal = false; 1407 fCalculatedExitNormal = false;
1412 1408
1413 fExitNormal = G4ThreeVector(0,0, 1409 fExitNormal = G4ThreeVector(0,0,0);
1414 fGrandMotherExitNormal = G4ThreeVector(0,0, 1410 fGrandMotherExitNormal = G4ThreeVector(0,0,0);
1415 fExitNormalGlobalFrame = G4ThreeVector(0,0, 1411 fExitNormalGlobalFrame = G4ThreeVector(0,0,0);
1416 1412
1417 fPreviousSftOrigin = G4ThreeVector(0,0, 1413 fPreviousSftOrigin = G4ThreeVector(0,0,0);
1418 fPreviousSafety = 0.0; 1414 fPreviousSafety = 0.0;
1419 1415
1420 fNumberZeroSteps = 0; 1416 fNumberZeroSteps = 0;
1421 1417
1422 fBlockedPhysicalVolume = nullptr; << 1418 fBlockedPhysicalVolume = 0;
1423 fBlockedReplicaNo = -1; 1419 fBlockedReplicaNo = -1;
1424 1420
1425 fLastLocatedPointLocal = G4ThreeVector( kIn 1421 fLastLocatedPointLocal = G4ThreeVector( kInfinity, -kInfinity, 0.0 );
1426 fLocatedOutsideWorld = false; 1422 fLocatedOutsideWorld = false;
1427 } 1423 }
1428 1424
1429 // ****************************************** 1425 // ********************************************************************
1430 // SetupHierarchy 1426 // SetupHierarchy
1431 // 1427 //
1432 // Renavigates & resets hierarchy described b 1428 // Renavigates & resets hierarchy described by current history
1433 // o Reset volumes 1429 // o Reset volumes
1434 // o Recompute transforms and/or solids of re 1430 // o Recompute transforms and/or solids of replicated/parameterised volumes
1435 // ****************************************** 1431 // ********************************************************************
1436 // 1432 //
1437 void G4ITNavigator2::SetupHierarchy() 1433 void G4ITNavigator2::SetupHierarchy()
1438 { 1434 {
1439 G4int i; 1435 G4int i;
1440 const auto cdepth = (G4int)fHistory.GetDep << 1436 const G4int cdepth = fHistory.GetDepth();
1441 G4VPhysicalVolume *current; 1437 G4VPhysicalVolume *current;
1442 G4VSolid *pSolid; 1438 G4VSolid *pSolid;
1443 G4VPVParameterisation *pParam; 1439 G4VPVParameterisation *pParam;
1444 1440
1445 for ( i=1; i<=cdepth; i++ ) 1441 for ( i=1; i<=cdepth; i++ )
1446 { 1442 {
1447 current = fHistory.GetVolume(i); 1443 current = fHistory.GetVolume(i);
1448 switch ( fHistory.GetVolumeType(i) ) 1444 switch ( fHistory.GetVolumeType(i) )
1449 { 1445 {
1450 case kNormal: 1446 case kNormal:
1451 break; 1447 break;
1452 case kReplica: 1448 case kReplica:
1453 freplicaNav.ComputeTransformation(fHi 1449 freplicaNav.ComputeTransformation(fHistory.GetReplicaNo(i), current);
1454 break; 1450 break;
1455 case kParameterised: 1451 case kParameterised:
1456 { <<
1457 G4int replicaNo; 1452 G4int replicaNo;
1458 pParam = current->GetParameterisation 1453 pParam = current->GetParameterisation();
1459 replicaNo = fHistory.GetReplicaNo(i); 1454 replicaNo = fHistory.GetReplicaNo(i);
1460 pSolid = pParam->ComputeSolid(replica 1455 pSolid = pParam->ComputeSolid(replicaNo, current);
1461 1456
1462 // Set up dimensions & transform in s 1457 // Set up dimensions & transform in solid/physical volume
1463 // 1458 //
1464 pSolid->ComputeDimensions(pParam, rep 1459 pSolid->ComputeDimensions(pParam, replicaNo, current);
1465 pParam->ComputeTransformation(replica 1460 pParam->ComputeTransformation(replicaNo, current);
1466 1461
1467 G4TouchableHistory *pTouchable= nullp << 1462 G4TouchableHistory *pTouchable= 0;
1468 if( pParam->IsNested() ) 1463 if( pParam->IsNested() )
1469 { 1464 {
1470 pTouchable= new G4TouchableHistory 1465 pTouchable= new G4TouchableHistory( fHistory );
1471 pTouchable->MoveUpHistory(); // Mo 1466 pTouchable->MoveUpHistory(); // Move up to the parent level
1472 // Adequate only if Nested at th 1467 // Adequate only if Nested at the Branch level (last)
1473 // To extend to other cases: 1468 // To extend to other cases:
1474 // pTouchable->MoveUpHistory(cdept 1469 // pTouchable->MoveUpHistory(cdepth-i-1);
1475 // Move to the parent level of * 1470 // Move to the parent level of *Current* level
1476 // Could replace this line and c 1471 // Could replace this line and constructor with a revised
1477 // c-tor for History(levels to d 1472 // c-tor for History(levels to drop)
1478 } 1473 }
1479 // Set up the correct solid and mater 1474 // Set up the correct solid and material in Logical Volume
1480 // 1475 //
1481 G4LogicalVolume *pLogical = current-> 1476 G4LogicalVolume *pLogical = current->GetLogicalVolume();
1482 pLogical->SetSolid( pSolid ); 1477 pLogical->SetSolid( pSolid );
1483 pLogical->UpdateMaterial( pParam -> 1478 pLogical->UpdateMaterial( pParam ->
1484 ComputeMaterial(replicaNo, current, 1479 ComputeMaterial(replicaNo, current, pTouchable) );
1485 delete pTouchable; 1480 delete pTouchable;
1486 } << 1481 break;
1487 break; <<
1488 <<
1489 case kExternal: <<
1490 G4Exception("G4ITNavigator2::SetupHie <<
1491 "GeomNav0001", FatalExcep <<
1492 "Not applicable for exter <<
1493 break; <<
1494 <<
1495 } 1482 }
1496 } 1483 }
1497 } 1484 }
1498 1485
1499 // ****************************************** 1486 // ********************************************************************
1500 // GetLocalExitNormal 1487 // GetLocalExitNormal
1501 // 1488 //
1502 // Obtains the Normal vector to a surface (in 1489 // Obtains the Normal vector to a surface (in local coordinates)
1503 // pointing out of previous volume and into c 1490 // pointing out of previous volume and into current volume
1504 // ****************************************** 1491 // ********************************************************************
1505 // 1492 //
1506 G4ThreeVector G4ITNavigator2::GetLocalExitNor 1493 G4ThreeVector G4ITNavigator2::GetLocalExitNormal( G4bool* valid )
1507 { 1494 {
1508 CheckNavigatorStateIsValid(); 1495 CheckNavigatorStateIsValid();
1509 G4ThreeVector ExitNormal(0.,0.,0.); 1496 G4ThreeVector ExitNormal(0.,0.,0.);
1510 G4VSolid *currentSolid=nullptr; << 1497 G4VSolid *currentSolid=0;
1511 G4LogicalVolume *candidateLogical; 1498 G4LogicalVolume *candidateLogical;
1512 if ( fLastTriedStepComputation ) 1499 if ( fLastTriedStepComputation )
1513 { 1500 {
1514 // use fLastLocatedPointLocal and next ca 1501 // use fLastLocatedPointLocal and next candidate volume
1515 // 1502 //
1516 G4ThreeVector nextSolidExitNormal(0.,0.,0 1503 G4ThreeVector nextSolidExitNormal(0.,0.,0.);
1517 1504
1518 if( fEntering && (fBlockedPhysicalVolume! << 1505 if( fEntering && (fBlockedPhysicalVolume!=0) )
1519 { 1506 {
1520 candidateLogical= fBlockedPhysicalVolum 1507 candidateLogical= fBlockedPhysicalVolume->GetLogicalVolume();
1521 if( candidateLogical != nullptr ) << 1508 if( candidateLogical )
1522 { 1509 {
1523 // fLastStepEndPointLocal is in the c 1510 // fLastStepEndPointLocal is in the coordinates of the mother
1524 // we need it in the daughter's coord 1511 // we need it in the daughter's coordinate system.
1525 1512
1526 // The following code should also wor 1513 // The following code should also work in case of Replica
1527 { 1514 {
1528 // First transform fLastLocatedPoin 1515 // First transform fLastLocatedPointLocal to the new daughter
1529 // coordinates 1516 // coordinates
1530 // 1517 //
1531 G4AffineTransform MotherToDaughterT 1518 G4AffineTransform MotherToDaughterTransform=
1532 GetMotherToDaughterTransform( fBl 1519 GetMotherToDaughterTransform( fBlockedPhysicalVolume,
1533 fBl 1520 fBlockedReplicaNo,
1534 Vol 1521 VolumeType(fBlockedPhysicalVolume) );
1535 G4ThreeVector daughterPointOwnLocal 1522 G4ThreeVector daughterPointOwnLocal=
1536 MotherToDaughterTransform.Transfo 1523 MotherToDaughterTransform.TransformPoint( fLastStepEndPointLocal );
1537 1524
1538 // OK if it is a parameterised volu 1525 // OK if it is a parameterised volume
1539 // 1526 //
1540 EInside inSideIt; 1527 EInside inSideIt;
1541 G4bool onSurface; 1528 G4bool onSurface;
1542 G4double safety= -1.0; 1529 G4double safety= -1.0;
1543 currentSolid= candidateLogical->Get 1530 currentSolid= candidateLogical->GetSolid();
1544 inSideIt = currentSolid->Inside(d 1531 inSideIt = currentSolid->Inside(daughterPointOwnLocal);
1545 onSurface = (inSideIt == kSurface) 1532 onSurface = (inSideIt == kSurface);
1546 if( ! onSurface ) 1533 if( ! onSurface )
1547 { 1534 {
1548 if( inSideIt == kOutside ) 1535 if( inSideIt == kOutside )
1549 { 1536 {
1550 safety = (currentSolid->Distanc 1537 safety = (currentSolid->DistanceToIn(daughterPointOwnLocal));
1551 onSurface = safety < 100.0 * kC 1538 onSurface = safety < 100.0 * kCarTolerance;
1552 } 1539 }
1553 else if (inSideIt == kInside ) 1540 else if (inSideIt == kInside )
1554 { 1541 {
1555 safety = (currentSolid->Distanc 1542 safety = (currentSolid->DistanceToOut(daughterPointOwnLocal));
1556 onSurface = safety < 100.0 * kC 1543 onSurface = safety < 100.0 * kCarTolerance;
1557 } 1544 }
1558 } 1545 }
1559 1546
1560 if( onSurface ) 1547 if( onSurface )
1561 { 1548 {
1562 nextSolidExitNormal = 1549 nextSolidExitNormal =
1563 currentSolid->SurfaceNormal(dau 1550 currentSolid->SurfaceNormal(daughterPointOwnLocal);
1564 1551
1565 // Entering the solid ==> opposit 1552 // Entering the solid ==> opposite
1566 // 1553 //
1567 ExitNormal = -nextSolidExitNormal 1554 ExitNormal = -nextSolidExitNormal;
1568 fCalculatedExitNormal= true; 1555 fCalculatedExitNormal= true;
1569 } 1556 }
1570 else 1557 else
1571 { 1558 {
1572 #ifdef G4VERBOSE 1559 #ifdef G4VERBOSE
1573 if(( fVerbose == 1 ) && ( fCheck 1560 if(( fVerbose == 1 ) && ( fCheck ))
1574 { 1561 {
1575 std::ostringstream message; 1562 std::ostringstream message;
1576 message << "Point not on surfac 1563 message << "Point not on surface ! " << G4endl
1577 << " Point = 1564 << " Point = "
1578 << daughterPointOwnLoca 1565 << daughterPointOwnLocal << G4endl
1579 << " Physical volume = 1566 << " Physical volume = "
1580 << fBlockedPhysicalVolu 1567 << fBlockedPhysicalVolume->GetName() << G4endl
1581 << " Logical volume = 1568 << " Logical volume = "
1582 << candidateLogical->Ge 1569 << candidateLogical->GetName() << G4endl
1583 << " Solid = 1570 << " Solid = " << currentSolid->GetName()
1584 << " Type = 1571 << " Type = "
1585 << currentSolid->GetEnt 1572 << currentSolid->GetEntityType() << G4endl
1586 << *currentSolid << G4e 1573 << *currentSolid << G4endl;
1587 if( inSideIt == kOutside ) 1574 if( inSideIt == kOutside )
1588 { 1575 {
1589 message << "Point is Outside. 1576 message << "Point is Outside. " << G4endl
1590 << " Safety (from ou 1577 << " Safety (from outside) = " << safety << G4endl;
1591 } 1578 }
1592 else // if( inSideIt == kInside 1579 else // if( inSideIt == kInside )
1593 { 1580 {
1594 message << "Point is Inside. 1581 message << "Point is Inside. " << G4endl
1595 << " Safety (from in 1582 << " Safety (from inside) = " << safety << G4endl;
1596 } 1583 }
1597 G4Exception("G4ITNavigator2::Ge 1584 G4Exception("G4ITNavigator2::GetLocalExitNormal()", "GeomNav1001",
1598 JustWarning, messag 1585 JustWarning, message);
1599 } 1586 }
1600 #endif 1587 #endif
1601 } 1588 }
1602 *valid = onSurface; // was =tru 1589 *valid = onSurface; // was =true;
1603 } 1590 }
1604 } 1591 }
1605 } 1592 }
1606 else if ( fExiting ) 1593 else if ( fExiting )
1607 { 1594 {
1608 ExitNormal = fGrandMotherExitNormal; 1595 ExitNormal = fGrandMotherExitNormal;
1609 *valid = true; 1596 *valid = true;
1610 fCalculatedExitNormal= true; // Should 1597 fCalculatedExitNormal= true; // Should be true already
1611 } 1598 }
1612 else // i.e. ( fBlockedPhysicalVolume = 1599 else // i.e. ( fBlockedPhysicalVolume == 0 )
1613 { 1600 {
1614 *valid = false; 1601 *valid = false;
1615 G4Exception("G4ITNavigator2::GetLocalEx 1602 G4Exception("G4ITNavigator2::GetLocalExitNormal()",
1616 "GeomNav0003", JustWarning, 1603 "GeomNav0003", JustWarning,
1617 "Incorrect call to GetLocal 1604 "Incorrect call to GetLocalSurfaceNormal." );
1618 } 1605 }
1619 } 1606 }
1620 else // ( ! fLastTriedStepComputation ) ie 1607 else // ( ! fLastTriedStepComputation ) ie. last call was to Locate
1621 { 1608 {
1622 if ( EnteredDaughterVolume() ) 1609 if ( EnteredDaughterVolume() )
1623 { 1610 {
1624 G4VSolid* daughterSolid =fHistory.GetTo 1611 G4VSolid* daughterSolid =fHistory.GetTopVolume()->GetLogicalVolume()
1625 1612 ->GetSolid();
1626 ExitNormal= -(daughterSolid->SurfaceNor 1613 ExitNormal= -(daughterSolid->SurfaceNormal(fLastLocatedPointLocal));
1627 if( std::fabs(ExitNormal.mag2()-1.0 ) > 1614 if( std::fabs(ExitNormal.mag2()-1.0 ) > CLHEP::perMillion )
1628 { 1615 {
1629 G4ExceptionDescription desc; 1616 G4ExceptionDescription desc;
1630 desc << " Parameters of solid: " << * 1617 desc << " Parameters of solid: " << *daughterSolid
1631 << " Point for surface = " << fL 1618 << " Point for surface = " << fLastLocatedPointLocal << std::endl;
1632 G4Exception("G4ITNavigator2::GetLocal 1619 G4Exception("G4ITNavigator2::GetLocalExitNormal()",
1633 "GeomNav0003", FatalExcep 1620 "GeomNav0003", FatalException, desc,
1634 "Surface Normal returned 1621 "Surface Normal returned by Solid is not a Unit Vector." );
1635 } 1622 }
1636 fCalculatedExitNormal= true; 1623 fCalculatedExitNormal= true;
1637 *valid = true; 1624 *valid = true;
1638 } 1625 }
1639 else 1626 else
1640 { 1627 {
1641 if( fExitedMother ) 1628 if( fExitedMother )
1642 { 1629 {
1643 ExitNormal = fGrandMotherExitNormal; 1630 ExitNormal = fGrandMotherExitNormal;
1644 *valid = true; 1631 *valid = true;
1645 fCalculatedExitNormal= true; 1632 fCalculatedExitNormal= true;
1646 } 1633 }
1647 else // We are not at a boundary. Exit 1634 else // We are not at a boundary. ExitNormal remains (0,0,0)
1648 { 1635 {
1649 *valid = false; 1636 *valid = false;
1650 fCalculatedExitNormal= false; 1637 fCalculatedExitNormal= false;
1651 G4ExceptionDescription message; 1638 G4ExceptionDescription message;
1652 message << "Function called when *NOT 1639 message << "Function called when *NOT* at a Boundary." << G4endl;
1653 G4Exception("G4ITNavigator2::GetLocal 1640 G4Exception("G4ITNavigator2::GetLocalExitNormal()",
1654 "GeomNav0003", JustWarnin 1641 "GeomNav0003", JustWarning, message);
1655 } 1642 }
1656 } 1643 }
1657 } 1644 }
1658 return ExitNormal; 1645 return ExitNormal;
1659 } 1646 }
1660 1647
1661 // ****************************************** 1648 // ********************************************************************
1662 // GetMotherToDaughterTransform 1649 // GetMotherToDaughterTransform
1663 // 1650 //
1664 // Obtains the mother to daughter affine tran 1651 // Obtains the mother to daughter affine transformation
1665 // ****************************************** 1652 // ********************************************************************
1666 // 1653 //
1667 G4AffineTransform 1654 G4AffineTransform
1668 G4ITNavigator2::GetMotherToDaughterTransform( 1655 G4ITNavigator2::GetMotherToDaughterTransform( G4VPhysicalVolume *pEnteringPhysVol, // not Const
1669 G4 1656 G4int enteringReplicaNo,
1670 EV 1657 EVolume enteringVolumeType )
1671 { 1658 {
1672 CheckNavigatorStateIsValid(); 1659 CheckNavigatorStateIsValid();
1673 switch (enteringVolumeType) 1660 switch (enteringVolumeType)
1674 { 1661 {
1675 case kNormal: // Nothing is needed to pr 1662 case kNormal: // Nothing is needed to prepare the transformation
1676 break; // It is stored already in 1663 break; // It is stored already in the physical volume (placement)
1677 case kReplica: // Sets the transform in t 1664 case kReplica: // Sets the transform in the Replica - tbc
1678 G4Exception("G4ITNavigator2::GetMotherT 1665 G4Exception("G4ITNavigator2::GetMotherToDaughterTransform()",
1679 "GeomNav0001", FatalExcepti 1666 "GeomNav0001", FatalException,
1680 "Method NOT Implemented yet 1667 "Method NOT Implemented yet for replica volumes.");
1681 break; 1668 break;
1682 case kParameterised: 1669 case kParameterised:
1683 if( pEnteringPhysVol->GetRegularStructu 1670 if( pEnteringPhysVol->GetRegularStructureId() == 0 )
1684 { 1671 {
1685 G4VPVParameterisation *pParam = 1672 G4VPVParameterisation *pParam =
1686 pEnteringPhysVol->GetParameterisati 1673 pEnteringPhysVol->GetParameterisation();
1687 G4VSolid* pSolid = 1674 G4VSolid* pSolid =
1688 pParam->ComputeSolid(enteringReplic 1675 pParam->ComputeSolid(enteringReplicaNo, pEnteringPhysVol);
1689 pSolid->ComputeDimensions(pParam, ent 1676 pSolid->ComputeDimensions(pParam, enteringReplicaNo, pEnteringPhysVol);
1690 1677
1691 // Sets the transform in the Paramete 1678 // Sets the transform in the Parameterisation
1692 // 1679 //
1693 pParam->ComputeTransformation(enterin 1680 pParam->ComputeTransformation(enteringReplicaNo, pEnteringPhysVol);
1694 1681
1695 // Set the correct solid and material 1682 // Set the correct solid and material in Logical Volume
1696 // 1683 //
1697 G4LogicalVolume* pLogical = pEntering 1684 G4LogicalVolume* pLogical = pEnteringPhysVol->GetLogicalVolume();
1698 pLogical->SetSolid( pSolid ); 1685 pLogical->SetSolid( pSolid );
1699 } 1686 }
1700 break; 1687 break;
1701 case kExternal: <<
1702 G4Exception("G4ITNavigator2::GetMothe <<
1703 "GeomNav0001", FatalExcep <<
1704 "Not applicable for exter <<
1705 break; <<
1706 } 1688 }
1707 return G4AffineTransform(pEnteringPhysVol-> 1689 return G4AffineTransform(pEnteringPhysVol->GetRotation(),
1708 pEnteringPhysVol-> 1690 pEnteringPhysVol->GetTranslation()).Invert();
1709 } 1691 }
1710 1692
1711 // ****************************************** 1693 // ********************************************************************
1712 // GetLocalExitNormalAndCheck 1694 // GetLocalExitNormalAndCheck
1713 // 1695 //
1714 // Obtains the Normal vector to a surface (in 1696 // Obtains the Normal vector to a surface (in local coordinates)
1715 // pointing out of previous volume and into c 1697 // pointing out of previous volume and into current volume, and
1716 // checks the current point against expected 1698 // checks the current point against expected 'local' value.
1717 // ****************************************** 1699 // ********************************************************************
1718 // 1700 //
1719 G4ThreeVector G4ITNavigator2:: 1701 G4ThreeVector G4ITNavigator2::
1720 GetLocalExitNormalAndCheck( 1702 GetLocalExitNormalAndCheck(
1721 #ifdef G4DEBUG_NAVIGATION 1703 #ifdef G4DEBUG_NAVIGATION
1722 const G4ThreeVecto 1704 const G4ThreeVector& ExpectedBoundaryPointGlobal,
1723 #else 1705 #else
1724 const G4ThreeVecto 1706 const G4ThreeVector&,
1725 #endif 1707 #endif
1726 G4bool* 1708 G4bool* pValid)
1727 { 1709 {
1728 CheckNavigatorStateIsValid(); 1710 CheckNavigatorStateIsValid();
1729 #ifdef G4DEBUG_NAVIGATION 1711 #ifdef G4DEBUG_NAVIGATION
1730 // Check Current point against expected 'lo 1712 // Check Current point against expected 'local' value
1731 // 1713 //
1732 if ( fLastTriedStepComputation ) 1714 if ( fLastTriedStepComputation )
1733 { 1715 {
1734 G4ThreeVector ExpectedBoundaryPointLocal; 1716 G4ThreeVector ExpectedBoundaryPointLocal;
1735 1717
1736 const G4AffineTransform& GlobalToLocal= G 1718 const G4AffineTransform& GlobalToLocal= GetGlobalToLocalTransform();
1737 ExpectedBoundaryPointLocal = 1719 ExpectedBoundaryPointLocal =
1738 GlobalToLocal.TransformPoint( ExpectedB 1720 GlobalToLocal.TransformPoint( ExpectedBoundaryPointGlobal );
1739 1721
1740 // Add here: Comparison against expected 1722 // Add here: Comparison against expected position,
1741 // i.e. the endpoint of Comput 1723 // i.e. the endpoint of ComputeStep
1742 } 1724 }
1743 #endif 1725 #endif
1744 1726
1745 return GetLocalExitNormal( pValid); 1727 return GetLocalExitNormal( pValid);
1746 } 1728 }
1747 1729
1748 // ****************************************** 1730 // ********************************************************************
1749 // GetGlobalExitNormal 1731 // GetGlobalExitNormal
1750 // 1732 //
1751 // Obtains the Normal vector to a surface (in 1733 // Obtains the Normal vector to a surface (in global coordinates)
1752 // pointing out of previous volume and into c 1734 // pointing out of previous volume and into current volume
1753 // ****************************************** 1735 // ********************************************************************
1754 // 1736 //
1755 G4ThreeVector 1737 G4ThreeVector
1756 G4ITNavigator2::GetGlobalExitNormal(const G4T 1738 G4ITNavigator2::GetGlobalExitNormal(const G4ThreeVector& IntersectPointGlobal,
1757 G4bool 1739 G4bool* pNormalCalculated)
1758 { 1740 {
1759 CheckNavigatorStateIsValid(); 1741 CheckNavigatorStateIsValid();
1760 G4bool validNormal; 1742 G4bool validNormal;
1761 G4ThreeVector localNormal, globalNormal; 1743 G4ThreeVector localNormal, globalNormal;
1762 G4bool usingStored; 1744 G4bool usingStored;
1763 1745
1764 usingStored= 1746 usingStored=
1765 fCalculatedExitNormal && 1747 fCalculatedExitNormal &&
1766 ( ( fLastTriedStepComputation && fExiti 1748 ( ( fLastTriedStepComputation && fExiting) // Just calculated it
1767 || 1749 || // No locate in between
1768 (!fLastTriedStepComputation 1750 (!fLastTriedStepComputation
1769 && (IntersectPointGlobal-fStepEndPo 1751 && (IntersectPointGlobal-fStepEndPoint).mag2()
1770 < (10.0*kCarTolerance*kCarTolera 1752 < (10.0*kCarTolerance*kCarTolerance)
1771 ) // Calculated it 'just' before & t 1753 ) // Calculated it 'just' before & then called locate
1772 // but it did not move position 1754 // but it did not move position
1773 ); 1755 );
1774 1756
1775 if( usingStored ) 1757 if( usingStored )
1776 { 1758 {
1777 // This was computed in ComputeStep -- an 1759 // This was computed in ComputeStep -- and only on arrival at boundary
1778 // 1760 //
1779 globalNormal = fExitNormalGlobalFrame; 1761 globalNormal = fExitNormalGlobalFrame;
1780 G4double normMag2 = globalNormal.mag2(); 1762 G4double normMag2 = globalNormal.mag2();
1781 if( std::fabs ( normMag2 - 1.0 ) < perMil 1763 if( std::fabs ( normMag2 - 1.0 ) < perMillion ) // Value is good
1782 { 1764 {
1783 *pNormalCalculated = true; // ComputeS 1765 *pNormalCalculated = true; // ComputeStep always computes it if Exiting
1784 // (fExitin 1766 // (fExiting==true)
1785 } 1767 }
1786 else 1768 else
1787 { 1769 {
1788 G4ExceptionDescription message; 1770 G4ExceptionDescription message;
1789 1771
1790 message << " ERROR> Expected normal-gl 1772 message << " ERROR> Expected normal-global-frame to valid (unit vector) "
1791 << " - but |normal| = " << s 1773 << " - but |normal| = " << std::sqrt(normMag2)
1792 << " - and |normal|^ = " << 1774 << " - and |normal|^ = " << normMag2
1793 << " which differs from 1.0 by 1775 << " which differs from 1.0 by " << normMag2 - 1.0 << G4endl
1794 << " n = " << fExitNormalGlo 1776 << " n = " << fExitNormalGlobalFrame << G4endl;
1795 message << "========================== 1777 message << "============================================================"
1796 << G4endl; 1778 << G4endl;
1797 G4int oldVerbose = fVerbose; 1779 G4int oldVerbose = fVerbose;
1798 fVerbose=4; 1780 fVerbose=4;
1799 message << " State of Navigator: " < 1781 message << " State of Navigator: " << G4endl;
1800 message << *this << G4endl; 1782 message << *this << G4endl;
1801 fVerbose = oldVerbose; 1783 fVerbose = oldVerbose;
1802 message << "========================== 1784 message << "============================================================"
1803 << G4endl; 1785 << G4endl;
1804 1786
1805 G4Exception("G4ITNavigator2::GetGlobal 1787 G4Exception("G4ITNavigator2::GetGlobalExitNormal()",
1806 "GeomNav0003",JustWarning, 1788 "GeomNav0003",JustWarning, message,
1807 "Value obtained from stored glo 1789 "Value obtained from stored global-normal is not a unit vector.");
1808 1790
1809 // (Re)Compute it now -- as either it 1791 // (Re)Compute it now -- as either it was not computed, or it is wrong.
1810 // 1792 //
1811 localNormal = GetLocalExitNormalAndChe 1793 localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,
1812 1794 &validNormal);
1813 *pNormalCalculated = fCalculatedExitNo 1795 *pNormalCalculated = fCalculatedExitNormal;
1814 1796
1815 G4AffineTransform localToGlobal = GetL 1797 G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
1816 globalNormal = localToGlobal.Transform 1798 globalNormal = localToGlobal.TransformAxis( localNormal );
1817 } 1799 }
1818 } 1800 }
1819 else 1801 else
1820 { 1802 {
1821 localNormal = GetLocalExitNormalAndCheck( 1803 localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,&validNormal);
1822 *pNormalCalculated = fCalculatedExitNorma 1804 *pNormalCalculated = fCalculatedExitNormal;
1823 1805
1824 #ifdef G4DEBUG_NAVIGATION 1806 #ifdef G4DEBUG_NAVIGATION
1825 usingStored= false; 1807 usingStored= false;
1826 1808
1827 if( (!validNormal) && !fCalculatedExitNor 1809 if( (!validNormal) && !fCalculatedExitNormal)
1828 { 1810 {
1829 G4ExceptionDescription edN; 1811 G4ExceptionDescription edN;
1830 edN << " Calculated = " << fCalculated 1812 edN << " Calculated = " << fCalculatedExitNormal << G4endl;
1831 edN << " Entering= " << fEntering << 1813 edN << " Entering= " << fEntering << G4endl;
1832 G4int oldVerbose= this->GetVerboseLevel 1814 G4int oldVerbose= this->GetVerboseLevel();
1833 this->SetVerboseLevel(4); 1815 this->SetVerboseLevel(4);
1834 edN << " State of Navigator: " << G4e 1816 edN << " State of Navigator: " << G4endl;
1835 edN << *this << G4endl; 1817 edN << *this << G4endl;
1836 this->SetVerboseLevel( oldVerbose ); 1818 this->SetVerboseLevel( oldVerbose );
1837 1819
1838 G4Exception("G4ITNavigator2::GetGlobalE 1820 G4Exception("G4ITNavigator2::GetGlobalExitNormal()",
1839 "GeomNav0003", JustWarning, 1821 "GeomNav0003", JustWarning, edN,
1840 "LocalExitNormalAndCheck() 1822 "LocalExitNormalAndCheck() did not calculate Normal.");
1841 } 1823 }
1842 #endif 1824 #endif
1843 1825
1844 G4double localMag2= localNormal.mag2(); 1826 G4double localMag2= localNormal.mag2();
1845 if( validNormal && (std::fabs(localMag2- 1827 if( validNormal && (std::fabs(localMag2-1.0)) > CLHEP::perMillion )
1846 { 1828 {
1847 G4ExceptionDescription edN; 1829 G4ExceptionDescription edN;
1848 1830
1849 edN << "G4ITNavigator2::GetGlobalExitN 1831 edN << "G4ITNavigator2::GetGlobalExitNormal: "
1850 << " Using Local Normal - from ca 1832 << " Using Local Normal - from call to GetLocalExitNormalAndCheck. "
1851 << G4endl 1833 << G4endl
1852 << " Local Exit Normal : " << " 1834 << " Local Exit Normal : " << " || = " << std::sqrt(localMag2)
1853 << " vec = " << localNormal << G4e 1835 << " vec = " << localNormal << G4endl
1854 << " Global Exit Normal : " << " 1836 << " Global Exit Normal : " << " || = " << globalNormal.mag()
1855 << " vec = " << globalNormal << G4 1837 << " vec = " << globalNormal << G4endl;
1856 edN << " Calculated It = " << fC 1838 edN << " Calculated It = " << fCalculatedExitNormal << G4endl;
1857 1839
1858 G4Exception("G4ITNavigator2::GetGlobal 1840 G4Exception("G4ITNavigator2::GetGlobalExitNormal()",
1859 "GeomNav0003",JustWarning, 1841 "GeomNav0003",JustWarning, edN,
1860 "Value obtained from new l 1842 "Value obtained from new local *solid* is incorrect.");
1861 localNormal = localNormal.unit(); // S 1843 localNormal = localNormal.unit(); // Should we correct it ??
1862 } 1844 }
1863 G4AffineTransform localToGlobal = GetLoc 1845 G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
1864 globalNormal = localToGlobal.TransformAx 1846 globalNormal = localToGlobal.TransformAxis( localNormal );
1865 } 1847 }
1866 1848
1867 #ifdef G4DEBUG_NAVIGATION 1849 #ifdef G4DEBUG_NAVIGATION
1868 if( usingStored ) 1850 if( usingStored )
1869 { 1851 {
1870 G4ThreeVector globalNormAgn; 1852 G4ThreeVector globalNormAgn;
1871 1853
1872 localNormal= GetLocalExitNormalAndCheck(I 1854 localNormal= GetLocalExitNormalAndCheck(IntersectPointGlobal, &validNormal);
1873 1855
1874 G4AffineTransform localToGlobal = GetLoca 1856 G4AffineTransform localToGlobal = GetLocalToGlobalTransform();
1875 globalNormAgn = localToGlobal.TransformAx 1857 globalNormAgn = localToGlobal.TransformAxis( localNormal );
1876 1858
1877 // Check the value computed against fExit 1859 // Check the value computed against fExitNormalGlobalFrame
1878 G4ThreeVector diffNorm = globalNormAgn - 1860 G4ThreeVector diffNorm = globalNormAgn - fExitNormalGlobalFrame;
1879 if( diffNorm.mag2() > perMillion*CLHEP::p 1861 if( diffNorm.mag2() > perMillion*CLHEP::perMillion)
1880 { 1862 {
1881 G4ExceptionDescription edDfn; 1863 G4ExceptionDescription edDfn;
1882 edDfn << "Found difference in normals i 1864 edDfn << "Found difference in normals in case of exiting mother "
1883 << "- when Get is called after Co 1865 << "- when Get is called after ComputingStep " << G4endl;
1884 edDfn << " Magnitude of diff = " 1866 edDfn << " Magnitude of diff = " << diffNorm.mag() << G4endl;
1885 edDfn << " Normal stored (Global) 1867 edDfn << " Normal stored (Global) = " << fExitNormalGlobalFrame
1886 << G4endl; 1868 << G4endl;
1887 edDfn << " Global Computed from Local 1869 edDfn << " Global Computed from Local = " << globalNormAgn << G4endl;
1888 G4Exception("G4ITNavigator::GetGlobalEx 1870 G4Exception("G4ITNavigator::GetGlobalExitNormal()", "GeomNav0003",
1889 JustWarning, edDfn); 1871 JustWarning, edDfn);
1890 } 1872 }
1891 } 1873 }
1892 #endif 1874 #endif
1893 1875
1894 return globalNormal; 1876 return globalNormal;
1895 } 1877 }
1896 1878
1897 // To make the new Voxel Safety the default, 1879 // To make the new Voxel Safety the default, uncomment the next line
1898 #define G4NEW_SAFETY 1 1880 #define G4NEW_SAFETY 1
1899 1881
1900 // ****************************************** 1882 // ********************************************************************
1901 // ComputeSafety 1883 // ComputeSafety
1902 // 1884 //
1903 // It assumes that it will be 1885 // It assumes that it will be
1904 // i) called at the Point in the same volume 1886 // i) called at the Point in the same volume as the EndPoint of the
1905 // ComputeStep. 1887 // ComputeStep.
1906 // ii) after (or at the end of) ComputeStep O 1888 // ii) after (or at the end of) ComputeStep OR after the relocation.
1907 // ****************************************** 1889 // ********************************************************************
1908 // 1890 //
1909 G4double G4ITNavigator2::ComputeSafety( const 1891 G4double G4ITNavigator2::ComputeSafety( const G4ThreeVector &pGlobalpoint,
1910 const G4 1892 const G4double pMaxLength,
1911 const G4 1893 const G4bool keepState)
1912 { 1894 {
1913 CheckNavigatorStateIsValid(); 1895 CheckNavigatorStateIsValid();
1914 G4double newSafety = 0.0; 1896 G4double newSafety = 0.0;
1915 1897
1916 #ifdef G4DEBUG_NAVIGATION 1898 #ifdef G4DEBUG_NAVIGATION
1917 G4int oldcoutPrec = G4cout.precision(8); 1899 G4int oldcoutPrec = G4cout.precision(8);
1918 if( fVerbose > 0 ) 1900 if( fVerbose > 0 )
1919 { 1901 {
1920 G4cout << "*** G4ITNavigator2::ComputeSaf 1902 G4cout << "*** G4ITNavigator2::ComputeSafety: ***" << G4endl
1921 << " Called at point: " << pGlo 1903 << " Called at point: " << pGlobalpoint << G4endl;
1922 1904
1923 G4VPhysicalVolume *motherPhysical = fHis 1905 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1924 G4cout << " Volume = " << motherPhysic 1906 G4cout << " Volume = " << motherPhysical->GetName()
1925 << " - Maximum length = " << pMaxL 1907 << " - Maximum length = " << pMaxLength << G4endl;
1926 if( fVerbose >= 4 ) 1908 if( fVerbose >= 4 )
1927 { 1909 {
1928 G4cout << " ----- Upon entering Com 1910 G4cout << " ----- Upon entering Compute Safety:" << G4endl;
1929 PrintState(); 1911 PrintState();
1930 } 1912 }
1931 } 1913 }
1932 #endif 1914 #endif
1933 1915
1934 G4SaveNavigatorState* savedState(nullptr); << 1916 G4SaveNavigatorState* savedState(0);
1935 1917
1936 G4double distEndpointSq = (pGlobalpoint-fSt 1918 G4double distEndpointSq = (pGlobalpoint-fStepEndPoint).mag2();
1937 G4bool stayedOnEndpoint = distEndpointSq 1919 G4bool stayedOnEndpoint = distEndpointSq < kCarTolerance*kCarTolerance;
1938 G4bool endpointOnSurface = fEnteredDaught 1920 G4bool endpointOnSurface = fEnteredDaughter || fExitedMother;
1939 1921
1940 if( endpointOnSurface && stayedOnEndpoint ) 1922 if( endpointOnSurface && stayedOnEndpoint )
1941 { 1923 {
1942 #ifdef G4DEBUG_NAVIGATION 1924 #ifdef G4DEBUG_NAVIGATION
1943 if( fVerbose >= 2 ) 1925 if( fVerbose >= 2 )
1944 { 1926 {
1945 G4cout << " G4Navigator::ComputeSa 1927 G4cout << " G4Navigator::ComputeSafety() finds that point - "
1946 << pGlobalpoint << " - is on surface 1928 << pGlobalpoint << " - is on surface " << G4endl;
1947 if( fEnteredDaughter ) { G4cout << " 1929 if( fEnteredDaughter ) { G4cout << " entered new daughter volume"; }
1948 if( fExitedMother ) { G4cout << " 1930 if( fExitedMother ) { G4cout << " and exited previous volume."; }
1949 G4cout << G4endl; 1931 G4cout << G4endl;
1950 G4cout << " EndPoint was = " << fStep 1932 G4cout << " EndPoint was = " << fStepEndPoint << G4endl;
1951 } 1933 }
1952 #endif 1934 #endif
1953 newSafety = 0.0; 1935 newSafety = 0.0;
1954 // return newSafety; 1936 // return newSafety;
1955 } 1937 }
1956 else // if( !(endpointOnSurface && stayedOn 1938 else // if( !(endpointOnSurface && stayedOnEndpoint) )
1957 { 1939 {
1958 1940
1959 if (keepState) 1941 if (keepState)
1960 { 1942 {
1961 // SetSavedState(); 1943 // SetSavedState();
1962 savedState = new G4SaveNavigatorState(fpN 1944 savedState = new G4SaveNavigatorState(fpNavigatorState);
1963 } 1945 }
1964 1946
1965 // Pseudo-relocate to this point (updates 1947 // Pseudo-relocate to this point (updates voxel information only)
1966 // 1948 //
1967 LocateGlobalPointWithinVolume( pGlobalpoi 1949 LocateGlobalPointWithinVolume( pGlobalpoint );
1968 // --->> DANGER: Side effects on sub-na 1950 // --->> DANGER: Side effects on sub-navigator voxel information <<---
1969 // Could be replaced again by 'gr 1951 // Could be replaced again by 'granular' calls to sub-navigator
1970 // locates (similar side-effects, 1952 // locates (similar side-effects, but faster.
1971 // Solutions: 1953 // Solutions:
1972 // 1) Re-locate (to where?) 1954 // 1) Re-locate (to where?)
1973 // 2) Insure that the methods us 1955 // 2) Insure that the methods using (G4ComputeStep?)
1974 // does a relocation (if info 1956 // does a relocation (if information is disturbed only ?)
1975 1957
1976 #ifdef G4DEBUG_NAVIGATION 1958 #ifdef G4DEBUG_NAVIGATION
1977 if( fVerbose >= 2 ) 1959 if( fVerbose >= 2 )
1978 { 1960 {
1979 G4cout << " G4ITNavigator2::ComputeSaf 1961 G4cout << " G4ITNavigator2::ComputeSafety() relocates-in-volume to point: "
1980 << pGlobalpoint << G4endl; 1962 << pGlobalpoint << G4endl;
1981 } 1963 }
1982 #endif 1964 #endif
1983 G4VPhysicalVolume *motherPhysical = fHist 1965 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
1984 G4LogicalVolume *motherLogical = motherPh 1966 G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
1985 G4SmartVoxelHeader* pVoxelHeader = mother 1967 G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader();
1986 G4ThreeVector localPoint = ComputeLocalPo 1968 G4ThreeVector localPoint = ComputeLocalPoint(pGlobalpoint);
1987 1969
1988 if ( fHistory.GetTopVolumeType()!=kReplic 1970 if ( fHistory.GetTopVolumeType()!=kReplica )
1989 { 1971 {
1990 switch(CharacteriseDaughters(motherLogi 1972 switch(CharacteriseDaughters(motherLogical))
1991 { 1973 {
1992 case kNormal: 1974 case kNormal:
1993 if ( pVoxelHeader != nullptr ) << 1975 if ( pVoxelHeader )
1994 { 1976 {
1995 #ifdef G4NEW_SAFETY 1977 #ifdef G4NEW_SAFETY
1996 G4double safetyTwo = fpVoxelSafet 1978 G4double safetyTwo = fpVoxelSafety->ComputeSafety(localPoint,
1997 *m 1979 *motherPhysical, pMaxLength);
1998 newSafety= safetyTwo; // Faster 1980 newSafety= safetyTwo; // Faster and best available
1999 #else 1981 #else
2000 G4double safetyOldVoxel; 1982 G4double safetyOldVoxel;
2001 LocateGlobalPointWithinVolume(pGl <<
2002 safetyOldVoxel = 1983 safetyOldVoxel =
2003 fvoxelNav.ComputeSafety(localPo 1984 fvoxelNav.ComputeSafety(localPoint,fHistory,pMaxLength);
2004 newSafety= safetyOldVoxel; 1985 newSafety= safetyOldVoxel;
2005 #endif 1986 #endif
2006 } 1987 }
2007 else 1988 else
2008 { 1989 {
2009 newSafety=fnormalNav.ComputeSafet 1990 newSafety=fnormalNav.ComputeSafety(localPoint,fHistory,pMaxLength);
2010 } 1991 }
2011 break; 1992 break;
2012 case kParameterised: 1993 case kParameterised:
2013 if( GetDaughtersRegularStructureId( 1994 if( GetDaughtersRegularStructureId(motherLogical) != 1 )
2014 { 1995 {
2015 newSafety=fparamNav.ComputeSafety 1996 newSafety=fparamNav.ComputeSafety(localPoint,fHistory,pMaxLength);
2016 } 1997 }
2017 else // Regular structure 1998 else // Regular structure
2018 { 1999 {
2019 newSafety=fregularNav.ComputeSafe 2000 newSafety=fregularNav.ComputeSafety(localPoint,fHistory,pMaxLength);
2020 } 2001 }
2021 break; 2002 break;
2022 case kReplica: 2003 case kReplica:
2023 G4Exception("G4ITNavigator2::Comput 2004 G4Exception("G4ITNavigator2::ComputeSafety()", "GeomNav0001",
2024 FatalException, "Not ap 2005 FatalException, "Not applicable for replicated volumes.");
2025 break; 2006 break;
2026 case kExternal: <<
2027 G4Exception("G4ITNavigator2::Comput <<
2028 "GeomNav0001", FatalExc <<
2029 "Not applicable for ext <<
2030 break; <<
2031 } 2007 }
2032 } 2008 }
2033 else 2009 else
2034 { 2010 {
2035 newSafety = freplicaNav.ComputeSafety(p 2011 newSafety = freplicaNav.ComputeSafety(pGlobalpoint, localPoint,
2036 f 2012 fHistory, pMaxLength);
2037 } 2013 }
2038 2014
2039 if (keepState) 2015 if (keepState)
2040 { 2016 {
2041 *fpNavigatorState = *savedState; 2017 *fpNavigatorState = *savedState;
2042 delete savedState; 2018 delete savedState;
2043 // RestoreSavedState(); 2019 // RestoreSavedState();
2044 // This now overwrites the values of the 2020 // This now overwrites the values of the Safety 'sphere' (correction)
2045 } 2021 }
2046 2022
2047 // Remember last safety origin & value 2023 // Remember last safety origin & value
2048 // 2024 //
2049 // We overwrite the Safety 'sphere' - kee 2025 // We overwrite the Safety 'sphere' - keeping old behaviour
2050 fPreviousSftOrigin = pGlobalpoint; 2026 fPreviousSftOrigin = pGlobalpoint;
2051 fPreviousSafety = newSafety; 2027 fPreviousSafety = newSafety;
2052 } 2028 }
2053 2029
2054 #ifdef G4DEBUG_NAVIGATION 2030 #ifdef G4DEBUG_NAVIGATION
2055 if( fVerbose > 1 ) 2031 if( fVerbose > 1 )
2056 { 2032 {
2057 G4cout << " ---- Exiting ComputeSafety 2033 G4cout << " ---- Exiting ComputeSafety " << G4endl;
2058 if( fVerbose > 2 ) { PrintState(); } 2034 if( fVerbose > 2 ) { PrintState(); }
2059 G4cout << " Returned value of Safety = 2035 G4cout << " Returned value of Safety = " << newSafety << G4endl;
2060 } 2036 }
2061 G4cout.precision(oldcoutPrec); 2037 G4cout.precision(oldcoutPrec);
2062 #endif 2038 #endif
2063 2039
2064 return newSafety; 2040 return newSafety;
2065 } 2041 }
2066 2042
2067 2043
2068 // ****************************************** 2044 // ********************************************************************
2069 // RecheckDistanceToCurrentBoundary 2045 // RecheckDistanceToCurrentBoundary
2070 // 2046 //
2071 // Trial method for checking potential displa 2047 // Trial method for checking potential displacement for MS
2072 // Check position aDisplacedGlobalpoint, to s 2048 // Check position aDisplacedGlobalpoint, to see whether it is in the
2073 // current volume (mother). 2049 // current volume (mother).
2074 // If in mother, check distance to boundary a 2050 // If in mother, check distance to boundary along aNewDirection.
2075 // If in entering daughter, check distance ba 2051 // If in entering daughter, check distance back to boundary.
2076 // NOTE: 2052 // NOTE:
2077 // Can be called only after ComputeStep() is 2053 // Can be called only after ComputeStep() is called - before ReLocation
2078 // Deals only with current volume (and potent 2054 // Deals only with current volume (and potentially entered)
2079 // Caveats 2055 // Caveats
2080 // First VERSION: Does not consider daughter 2056 // First VERSION: Does not consider daughter volumes if it remained in mother
2081 // neither for checking whether it has 2057 // neither for checking whether it has exited current (mother) volume,
2082 // nor for determining distance to exit 2058 // nor for determining distance to exit this (mother) volume.
2083 // ****************************************** 2059 // ********************************************************************
2084 // 2060 //
2085 G4bool G4ITNavigator2::RecheckDistanceToCurre 2061 G4bool G4ITNavigator2::RecheckDistanceToCurrentBoundary(
2086 const G4ThreeVector &aDi 2062 const G4ThreeVector &aDisplacedGlobalPoint,
2087 const G4ThreeVector &aNe 2063 const G4ThreeVector &aNewDirection,
2088 const G4double ProposedM 2064 const G4double ProposedMove,
2089 G4double *prDistance, 2065 G4double *prDistance,
2090 G4double *prNewSafety) c 2066 G4double *prNewSafety) const
2091 { 2067 {
2092 G4ThreeVector localPosition = ComputeLocal 2068 G4ThreeVector localPosition = ComputeLocalPoint(aDisplacedGlobalPoint);
2093 G4ThreeVector localDirection = ComputeLocal 2069 G4ThreeVector localDirection = ComputeLocalAxis(aNewDirection);
2094 // G4double Step = kInfinity; 2070 // G4double Step = kInfinity;
2095 2071
2096 G4bool validExitNormal; 2072 G4bool validExitNormal;
2097 G4ThreeVector exitNormal; 2073 G4ThreeVector exitNormal;
2098 // Check against mother solid 2074 // Check against mother solid
2099 G4VPhysicalVolume *motherPhysical = fHisto 2075 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume();
2100 G4LogicalVolume *motherLogical = motherPhys 2076 G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume();
2101 2077
2102 #ifdef CHECK_ORDER_OF_METHODS 2078 #ifdef CHECK_ORDER_OF_METHODS
2103 if( ! fLastTriedStepComputation ) 2079 if( ! fLastTriedStepComputation )
2104 { 2080 {
2105 G4Exception("G4Navigator::RecheckDistanc 2081 G4Exception("G4Navigator::RecheckDistanceToCurrentBoundary()",
2106 "GeomNav0001", FatalExceptio 2082 "GeomNav0001", FatalException,
2107 "Method must be called after ComputeStep 2083 "Method must be called after ComputeStep(), before call to LocateMethod.");
2108 } 2084 }
2109 #endif 2085 #endif
2110 2086
2111 EInside locatedDaug; // = kUndefined; 2087 EInside locatedDaug; // = kUndefined;
2112 G4double daughterStep= DBL_MAX; 2088 G4double daughterStep= DBL_MAX;
2113 G4double daughterSafety= DBL_MAX; 2089 G4double daughterSafety= DBL_MAX;
2114 2090
2115 if( fEnteredDaughter ) 2091 if( fEnteredDaughter )
2116 { 2092 {
2117 if( motherLogical->CharacteriseDaughters 2093 if( motherLogical->CharacteriseDaughters() ==kReplica ) { return false; }
2118 2094
2119 // Track arrived at boundary of a daught 2095 // Track arrived at boundary of a daughter volume at
2120 // the last call of ComputeStep(). 2096 // the last call of ComputeStep().
2121 // In case the proposed displaced point 2097 // In case the proposed displaced point is inside this daughter,
2122 // it must backtrack at least to the e 2098 // it must backtrack at least to the entry point.
2123 // NOTE: No check is made against other 2099 // NOTE: No check is made against other daughter volumes. It is
2124 // assumed that the proposed displacem 2100 // assumed that the proposed displacement is small enough that
2125 // this is not needed. 2101 // this is not needed.
2126 2102
2127 // Must check boundary of current daught 2103 // Must check boundary of current daughter
2128 // 2104 //
2129 G4VPhysicalVolume *candPhysical= fBlocke 2105 G4VPhysicalVolume *candPhysical= fBlockedPhysicalVolume;
2130 G4LogicalVolume *candLogical= candPhysic 2106 G4LogicalVolume *candLogical= candPhysical->GetLogicalVolume();
2131 G4VSolid *candSolid= candLogica 2107 G4VSolid *candSolid= candLogical->GetSolid();
2132 2108
2133 G4AffineTransform nextLevelTrf(candPhysi 2109 G4AffineTransform nextLevelTrf(candPhysical->GetRotation(),
2134 candPhysi 2110 candPhysical->GetTranslation());
2135 2111
2136 G4ThreeVector dgPosition= nextLevelTrf. 2112 G4ThreeVector dgPosition= nextLevelTrf.TransformPoint(localPosition);
2137 G4ThreeVector dgDirection= nextLevelTrf. 2113 G4ThreeVector dgDirection= nextLevelTrf.TransformAxis(localDirection);
2138 locatedDaug = candSolid->Inside(dgPositi 2114 locatedDaug = candSolid->Inside(dgPosition);
2139 2115
2140 if( locatedDaug == kInside ) 2116 if( locatedDaug == kInside )
2141 { 2117 {
2142 // Reverse direction - and find first 2118 // Reverse direction - and find first exit. ( Is it valid?)
2143 // Must backtrack 2119 // Must backtrack
2144 G4double distanceBackOut = 2120 G4double distanceBackOut =
2145 candSolid->DistanceToOut(dgPosition 2121 candSolid->DistanceToOut(dgPosition,
2146 - dgDirect 2122 - dgDirection, // Reverse direction
2147 true, &val 2123 true, &validExitNormal, &exitNormal);
2148 daughterStep= - distanceBackOut; 2124 daughterStep= - distanceBackOut;
2149 // No check is made whether the par 2125 // No check is made whether the particle could have arrived at
2150 // at this location without encount 2126 // at this location without encountering another volume or
2151 // a different psurface of the curr 2127 // a different psurface of the current volume
2152 if( prNewSafety != nullptr ) << 2128 if( prNewSafety )
2153 { 2129 {
2154 daughterSafety= candSolid->Distanc 2130 daughterSafety= candSolid->DistanceToOut(dgPosition);
2155 } 2131 }
2156 } 2132 }
2157 else 2133 else
2158 { 2134 {
2159 if( locatedDaug == kOutside ) 2135 if( locatedDaug == kOutside )
2160 { 2136 {
2161 // See whether it still intersect 2137 // See whether it still intersects it
2162 // 2138 //
2163 daughterStep= candSolid->Distanc 2139 daughterStep= candSolid->DistanceToIn(dgPosition,
2164 2140 dgDirection);
2165 if( prNewSafety != nullptr ) << 2141 if( prNewSafety )
2166 { 2142 {
2167 daughterSafety= candSolid->Dist 2143 daughterSafety= candSolid->DistanceToIn(dgPosition);
2168 } 2144 }
2169 } 2145 }
2170 else 2146 else
2171 { 2147 {
2172 // The point remains on the surfac 2148 // The point remains on the surface of candidate solid
2173 // 2149 //
2174 daughterStep= 0.0; 2150 daughterStep= 0.0;
2175 daughterSafety= 0.0; 2151 daughterSafety= 0.0;
2176 } 2152 }
2177 } 2153 }
2178 2154
2179 // If trial point is in daughter (or on 2155 // If trial point is in daughter (or on its surface) we have the
2180 // answer, the rest is not relevant 2156 // answer, the rest is not relevant
2181 // 2157 //
2182 if( locatedDaug != kOutside ) 2158 if( locatedDaug != kOutside )
2183 { 2159 {
2184 *prDistance= daughterStep; 2160 *prDistance= daughterStep;
2185 if( prNewSafety != nullptr ) { *prNe << 2161 if( prNewSafety ) { *prNewSafety= daughterSafety; }
2186 return true; 2162 return true;
2187 } 2163 }
2188 // If ever extended, so that some type o 2164 // If ever extended, so that some type of mother cut daughter,
2189 // this would change 2165 // this would change
2190 } 2166 }
2191 2167
2192 G4VSolid *motherSolid= motherLogical->GetSo 2168 G4VSolid *motherSolid= motherLogical->GetSolid();
2193 2169
2194 G4double motherStep= DBL_MAX, motherSafety= 2170 G4double motherStep= DBL_MAX, motherSafety= DBL_MAX;
2195 2171
2196 // Check distance to boundary of mother 2172 // Check distance to boundary of mother
2197 // 2173 //
2198 if ( fHistory.GetTopVolumeType()!=kReplica 2174 if ( fHistory.GetTopVolumeType()!=kReplica )
2199 { 2175 {
2200 EInside locatedMoth = motherSolid->Insid 2176 EInside locatedMoth = motherSolid->Inside(localPosition);
2201 2177
2202 if( locatedMoth == kInside ) 2178 if( locatedMoth == kInside )
2203 { 2179 {
2204 motherSafety= motherSolid->DistanceTo 2180 motherSafety= motherSolid->DistanceToOut(localPosition);
2205 if( ProposedMove >= motherSafety ) 2181 if( ProposedMove >= motherSafety )
2206 { 2182 {
2207 motherStep= motherSolid->DistanceT 2183 motherStep= motherSolid->DistanceToOut(localPosition,
2208 2184 localDirection,
2209 true, &validExit 2185 true, &validExitNormal, &exitNormal);
2210 } 2186 }
2211 else 2187 else
2212 { 2188 {
2213 motherStep= ProposedMove; 2189 motherStep= ProposedMove;
2214 } 2190 }
2215 } 2191 }
2216 else if( locatedMoth == kOutside) 2192 else if( locatedMoth == kOutside)
2217 { 2193 {
2218 motherSafety= motherSolid->DistanceTo 2194 motherSafety= motherSolid->DistanceToIn(localPosition);
2219 if( ProposedMove >= motherSafety ) 2195 if( ProposedMove >= motherSafety )
2220 { 2196 {
2221 motherStep= - motherSolid->Distan 2197 motherStep= - motherSolid->DistanceToIn(localPosition,
2222 2198 -localDirection);
2223 } 2199 }
2224 } 2200 }
2225 else 2201 else
2226 { 2202 {
2227 motherSafety= 0.0; 2203 motherSafety= 0.0;
2228 *prDistance= 0.0; // On surface - n 2204 *prDistance= 0.0; // On surface - no move // motherStep;
2229 if( prNewSafety != nullptr ) { *prNe << 2205 if( prNewSafety ) { *prNewSafety= motherSafety; }
2230 return false; 2206 return false;
2231 } 2207 }
2232 } 2208 }
2233 else 2209 else
2234 { 2210 {
2235 return false; 2211 return false;
2236 } 2212 }
2237 2213
2238 *prDistance= std::min( motherStep, daughte 2214 *prDistance= std::min( motherStep, daughterStep );
2239 if( prNewSafety != nullptr ) << 2215 if( prNewSafety )
2240 { 2216 {
2241 *prNewSafety= std::min( motherSafety, da 2217 *prNewSafety= std::min( motherSafety, daughterSafety );
2242 } 2218 }
2243 return true; 2219 return true;
2244 } 2220 }
2245 2221
2246 2222
2247 // ****************************************** 2223 // ********************************************************************
2248 // CreateTouchableHistoryHandle 2224 // CreateTouchableHistoryHandle
2249 // ****************************************** 2225 // ********************************************************************
2250 // 2226 //
2251 G4TouchableHandle G4ITNavigator2::CreateTouch << 2227 G4TouchableHistoryHandle G4ITNavigator2::CreateTouchableHistoryHandle() const
2252 { 2228 {
2253 CheckNavigatorStateIsValid(); 2229 CheckNavigatorStateIsValid();
2254 return G4TouchableHandle( CreateTouchableHi << 2230 return G4TouchableHistoryHandle( CreateTouchableHistory() );
2255 } 2231 }
2256 2232
2257 // ****************************************** 2233 // ********************************************************************
2258 // PrintState 2234 // PrintState
2259 // ****************************************** 2235 // ********************************************************************
2260 // 2236 //
2261 void G4ITNavigator2::PrintState() const 2237 void G4ITNavigator2::PrintState() const
2262 { 2238 {
2263 CheckNavigatorStateIsValid(); 2239 CheckNavigatorStateIsValid();
2264 G4long oldcoutPrec = G4cout.precision(4); << 2240 G4int oldcoutPrec = G4cout.precision(4);
2265 if( fVerbose >= 4 ) 2241 if( fVerbose >= 4 )
2266 { 2242 {
2267 G4cout << "The current state of G4Navigat 2243 G4cout << "The current state of G4Navigator is: " << G4endl;
2268 G4cout << " ValidExitNormal= " << fValid 2244 G4cout << " ValidExitNormal= " << fValidExitNormal // << G4endl
2269 << " ExitNormal = " << fExitN 2245 << " ExitNormal = " << fExitNormal // << G4endl
2270 << " Exiting = " << fExiti 2246 << " Exiting = " << fExiting // << G4endl
2271 << " Entering = " << fEnter 2247 << " Entering = " << fEntering // << G4endl
2272 << " BlockedPhysicalVolume= " ; 2248 << " BlockedPhysicalVolume= " ;
2273 if (fBlockedPhysicalVolume==nullptr) << 2249 if (fBlockedPhysicalVolume==0)
2274 G4cout << "None"; 2250 G4cout << "None";
2275 else 2251 else
2276 G4cout << fBlockedPhysicalVolume->GetNa 2252 G4cout << fBlockedPhysicalVolume->GetName();
2277 G4cout << G4endl 2253 G4cout << G4endl
2278 << " BlockedReplicaNo = " << 2254 << " BlockedReplicaNo = " << fBlockedReplicaNo // << G4endl
2279 << " LastStepWasZero = " << 2255 << " LastStepWasZero = " << fLastStepWasZero // << G4endl
2280 << G4endl; 2256 << G4endl;
2281 } 2257 }
2282 if( ( 1 < fVerbose) && (fVerbose < 4) ) 2258 if( ( 1 < fVerbose) && (fVerbose < 4) )
2283 { 2259 {
2284 G4cout << G4endl; // Make sure to line up 2260 G4cout << G4endl; // Make sure to line up
2285 G4cout << std::setw(30) << " ExitNormal " 2261 G4cout << std::setw(30) << " ExitNormal " << " "
2286 << std::setw( 5) << " Valid " 2262 << std::setw( 5) << " Valid " << " "
2287 << std::setw( 9) << " Exiting " 2263 << std::setw( 9) << " Exiting " << " "
2288 << std::setw( 9) << " Entering" 2264 << std::setw( 9) << " Entering" << " "
2289 << std::setw(15) << " Blocked:Volu 2265 << std::setw(15) << " Blocked:Volume " << " "
2290 << std::setw( 9) << " ReplicaNo" 2266 << std::setw( 9) << " ReplicaNo" << " "
2291 << std::setw( 8) << " LastStepZero 2267 << std::setw( 8) << " LastStepZero " << " "
2292 << G4endl; 2268 << G4endl;
2293 G4cout << "( " << std::setw(7) << fExitNo 2269 G4cout << "( " << std::setw(7) << fExitNormal.x()
2294 << ", " << std::setw(7) << fExitNo 2270 << ", " << std::setw(7) << fExitNormal.y()
2295 << ", " << std::setw(7) << fExitNo 2271 << ", " << std::setw(7) << fExitNormal.z() << " ) "
2296 << std::setw( 5) << fValidExitNor 2272 << std::setw( 5) << fValidExitNormal << " "
2297 << std::setw( 9) << fExiting 2273 << std::setw( 9) << fExiting << " "
2298 << std::setw( 9) << fEntering 2274 << std::setw( 9) << fEntering << " ";
2299 if ( fBlockedPhysicalVolume==nullptr ) << 2275 if ( fBlockedPhysicalVolume==0 )
2300 { <<
2301 G4cout << std::setw(15) << "None"; 2276 G4cout << std::setw(15) << "None";
2302 } <<
2303 else 2277 else
2304 { <<
2305 G4cout << std::setw(15)<< fBlockedPhysi 2278 G4cout << std::setw(15)<< fBlockedPhysicalVolume->GetName();
2306 } << 2279 G4cout << std::setw( 9) << fBlockedReplicaNo << " "
2307 G4cout << std::setw( 9) << fBlockedRepli << 2280 << std::setw( 8) << fLastStepWasZero << " "
2308 << std::setw( 8) << fLastStepWasZ << 2281 << G4endl;
2309 << G4endl; <<
2310 } 2282 }
2311 if( fVerbose > 2 ) 2283 if( fVerbose > 2 )
2312 { 2284 {
2313 G4cout.precision(8); 2285 G4cout.precision(8);
2314 G4cout << " Current Localpoint = " << fLa 2286 G4cout << " Current Localpoint = " << fLastLocatedPointLocal << G4endl;
2315 G4cout << " PreviousSftOrigin = " << fPr 2287 G4cout << " PreviousSftOrigin = " << fPreviousSftOrigin << G4endl;
2316 G4cout << " PreviousSafety = " << fPr 2288 G4cout << " PreviousSafety = " << fPreviousSafety << G4endl;
2317 } 2289 }
2318 G4cout.precision(oldcoutPrec); 2290 G4cout.precision(oldcoutPrec);
2319 } 2291 }
2320 2292
2321 // ****************************************** 2293 // ********************************************************************
2322 // ComputeStepLog 2294 // ComputeStepLog
2323 // ****************************************** 2295 // ********************************************************************
2324 // 2296 //
2325 void G4ITNavigator2::ComputeStepLog(const G4T 2297 void G4ITNavigator2::ComputeStepLog(const G4ThreeVector& pGlobalpoint,
2326 G4doub 2298 G4double moveLenSq) const
2327 { 2299 {
2328 CheckNavigatorStateIsValid(); 2300 CheckNavigatorStateIsValid();
2329 // The following checks only make sense if 2301 // The following checks only make sense if the move is larger
2330 // than the tolerance. 2302 // than the tolerance.
2331 2303
2332 static const G4double fAccuracyForWarning 2304 static const G4double fAccuracyForWarning = kCarTolerance,
2333 fAccuracyForException 2305 fAccuracyForException = 1000*kCarTolerance;
2334 2306
2335 G4ThreeVector OriginalGlobalpoint = fHistor 2307 G4ThreeVector OriginalGlobalpoint = fHistory.GetTopTransform().Inverse().
2336 Transfo 2308 TransformPoint(fLastLocatedPointLocal);
2337 2309
2338 G4double shiftOriginSafSq = (fPreviousSftOr 2310 G4double shiftOriginSafSq = (fPreviousSftOrigin-pGlobalpoint).mag2();
2339 2311
2340 // Check that the starting point of this st 2312 // Check that the starting point of this step is
2341 // within the isotropic safety sphere of th 2313 // within the isotropic safety sphere of the last point
2342 // to a accuracy/precision given by fAccur 2314 // to a accuracy/precision given by fAccuracyForWarning.
2343 // If so give warning. 2315 // If so give warning.
2344 // If it fails by more than fAccuracyForE 2316 // If it fails by more than fAccuracyForException exit with error.
2345 // 2317 //
2346 if( shiftOriginSafSq >= sqr(fPreviousSafety 2318 if( shiftOriginSafSq >= sqr(fPreviousSafety) )
2347 { 2319 {
2348 G4double shiftOrigin = std::sqrt(shiftOri 2320 G4double shiftOrigin = std::sqrt(shiftOriginSafSq);
2349 G4double diffShiftSaf = shiftOrigin - fPr 2321 G4double diffShiftSaf = shiftOrigin - fPreviousSafety;
2350 2322
2351 if( diffShiftSaf > fAccuracyForWarning ) 2323 if( diffShiftSaf > fAccuracyForWarning )
2352 { 2324 {
2353 G4long oldcoutPrec= G4cout.precision(8) << 2325 G4int oldcoutPrec= G4cout.precision(8);
2354 G4long oldcerrPrec= G4cerr.precision(10 << 2326 G4int oldcerrPrec= G4cerr.precision(10);
2355 std::ostringstream message, suggestion; 2327 std::ostringstream message, suggestion;
2356 message << "Accuracy error or slightly 2328 message << "Accuracy error or slightly inaccurate position shift."
2357 << G4endl 2329 << G4endl
2358 << " The Step's starting po 2330 << " The Step's starting point has moved "
2359 << std::sqrt(moveLenSq)/mm << " 2331 << std::sqrt(moveLenSq)/mm << " mm " << G4endl
2360 << " since the last call to 2332 << " since the last call to a Locate method." << G4endl
2361 << " This has resulted in m 2333 << " This has resulted in moving "
2362 << shiftOrigin/mm << " mm " 2334 << shiftOrigin/mm << " mm "
2363 << " from the last point at whi 2335 << " from the last point at which the safety "
2364 << " was calculated " << G4 2336 << " was calculated " << G4endl
2365 << " which is more than the 2337 << " which is more than the computed safety= "
2366 << fPreviousSafety/mm << " mm 2338 << fPreviousSafety/mm << " mm at that point." << G4endl
2367 << " This difference is " 2339 << " This difference is "
2368 << diffShiftSaf/mm << " mm." << 2340 << diffShiftSaf/mm << " mm." << G4endl
2369 << " The tolerated accuracy 2341 << " The tolerated accuracy is "
2370 << fAccuracyForException/mm << 2342 << fAccuracyForException/mm << " mm.";
2371 2343
2372 suggestion << " "; 2344 suggestion << " ";
2373 static G4ThreadLocal G4int warnNow = 0; 2345 static G4ThreadLocal G4int warnNow = 0;
2374 if( ((++warnNow % 100) == 1) ) 2346 if( ((++warnNow % 100) == 1) )
2375 { 2347 {
2376 message << G4endl 2348 message << G4endl
2377 << " This problem can be due 2349 << " This problem can be due to either " << G4endl
2378 << " - a process that has p 2350 << " - a process that has proposed a displacement"
2379 << " larger than the current s 2351 << " larger than the current safety , or" << G4endl
2380 << " - inaccuracy in the co 2352 << " - inaccuracy in the computation of the safety";
2381 suggestion << "We suggest that you " 2353 suggestion << "We suggest that you " << G4endl
2382 << " - find i) what part 2354 << " - find i) what particle is being tracked, and "
2383 << " ii) through what part 2355 << " ii) through what part of your geometry " << G4endl
2384 << " for example by r 2356 << " for example by re-running this event with "
2385 << G4endl 2357 << G4endl
2386 << " /tracking/ver 2358 << " /tracking/verbose 1 " << G4endl
2387 << " - check which proc 2359 << " - check which processes you declare for"
2388 << " this particle (and lo 2360 << " this particle (and look at non-standard ones)"
2389 << G4endl 2361 << G4endl
2390 << " - in case, create a 2362 << " - in case, create a detailed logfile"
2391 << " of this event using:" 2363 << " of this event using:" << G4endl
2392 << " /tracking/ver 2364 << " /tracking/verbose 6 ";
2393 } 2365 }
2394 G4Exception("G4ITNavigator2::ComputeSte 2366 G4Exception("G4ITNavigator2::ComputeStep()",
2395 "GeomNav1002", JustWarning, 2367 "GeomNav1002", JustWarning,
2396 message, G4String(suggestio 2368 message, G4String(suggestion.str()));
2397 G4cout.precision(oldcoutPrec); 2369 G4cout.precision(oldcoutPrec);
2398 G4cerr.precision(oldcerrPrec); 2370 G4cerr.precision(oldcerrPrec);
2399 } 2371 }
2400 #ifdef G4DEBUG_NAVIGATION 2372 #ifdef G4DEBUG_NAVIGATION
2401 else 2373 else
2402 { 2374 {
2403 G4cerr << "WARNING - G4ITNavigator2::Co 2375 G4cerr << "WARNING - G4ITNavigator2::ComputeStep()" << G4endl
2404 << " The Step's startin 2376 << " The Step's starting point has moved "
2405 << std::sqrt(moveLenSq) << "," < 2377 << std::sqrt(moveLenSq) << "," << G4endl
2406 << " which has taken it 2378 << " which has taken it to the limit of"
2407 << " the current safety. " << G4 2379 << " the current safety. " << G4endl;
2408 } 2380 }
2409 #endif 2381 #endif
2410 } 2382 }
2411 G4double safetyPlus = fPreviousSafety + fAc 2383 G4double safetyPlus = fPreviousSafety + fAccuracyForException;
2412 if ( shiftOriginSafSq > sqr(safetyPlus) ) 2384 if ( shiftOriginSafSq > sqr(safetyPlus) )
2413 { 2385 {
2414 std::ostringstream message; 2386 std::ostringstream message;
2415 message << "May lead to a crash or unreli 2387 message << "May lead to a crash or unreliable results." << G4endl
2416 << " Position has shifted 2388 << " Position has shifted considerably without"
2417 << " notifying the navigator !" < 2389 << " notifying the navigator !" << G4endl
2418 << " Tolerated safety: " < 2390 << " Tolerated safety: " << safetyPlus << G4endl
2419 << " Computed shift : " < 2391 << " Computed shift : " << shiftOriginSafSq;
2420 G4Exception("G4ITNavigator2::ComputeStep( 2392 G4Exception("G4ITNavigator2::ComputeStep()", "GeomNav1002",
2421 JustWarning, message); 2393 JustWarning, message);
2422 } 2394 }
2423 } 2395 }
2424 2396
2425 // ****************************************** 2397 // ********************************************************************
2426 // Operator << 2398 // Operator <<
2427 // ****************************************** 2399 // ********************************************************************
2428 // 2400 //
2429 std::ostream& operator << (std::ostream &os,c 2401 std::ostream& operator << (std::ostream &os,const G4ITNavigator2 &n)
2430 { 2402 {
2431 // Old version did only the following: 2403 // Old version did only the following:
2432 // os << "Current History: " << G4endl << n 2404 // os << "Current History: " << G4endl << n.fHistory;
2433 // Old behaviour is recovered for fVerbose 2405 // Old behaviour is recovered for fVerbose = 0
2434 2406
2435 // Adapted from G4ITNavigator2::PrintState( 2407 // Adapted from G4ITNavigator2::PrintState() const
2436 2408
2437 G4long oldcoutPrec = os.precision(4); << 2409 G4int oldcoutPrec = os.precision(4);
2438 if( n.fVerbose >= 4 ) 2410 if( n.fVerbose >= 4 )
2439 { 2411 {
2440 os << "The current state of G4ITNavigator 2412 os << "The current state of G4ITNavigator2 is: " << G4endl;
2441 os << " ValidExitNormal= " << n.fValidEx 2413 os << " ValidExitNormal= " << n.fValidExitNormal << G4endl
2442 << " ExitNormal = " << n.fExitNormal 2414 << " ExitNormal = " << n.fExitNormal << G4endl
2443 << " Exiting = " << n.fExiting 2415 << " Exiting = " << n.fExiting << G4endl
2444 << " Entering = " << n.fEntering 2416 << " Entering = " << n.fEntering << G4endl
2445 << " BlockedPhysicalVolume= " ; 2417 << " BlockedPhysicalVolume= " ;
2446 << 2418 if (n.fBlockedPhysicalVolume==0)
2447 if (n.fBlockedPhysicalVolume==nullptr) <<
2448 { <<
2449 os << "None"; 2419 os << "None";
2450 } <<
2451 else 2420 else
2452 { <<
2453 os << n.fBlockedPhysicalVolume->GetName 2421 os << n.fBlockedPhysicalVolume->GetName();
2454 } <<
2455 <<
2456 os << G4endl 2422 os << G4endl
2457 << " BlockedReplicaNo = " << n.fBlo 2423 << " BlockedReplicaNo = " << n.fBlockedReplicaNo << G4endl
2458 << " LastStepWasZero = " << n.fLa 2424 << " LastStepWasZero = " << n.fLastStepWasZero << G4endl
2459 << G4endl; 2425 << G4endl;
2460 } 2426 }
2461 if( ( 1 < n.fVerbose) && (n.fVerbose < 4) ) 2427 if( ( 1 < n.fVerbose) && (n.fVerbose < 4) )
2462 { 2428 {
2463 os << G4endl; // Make sure to line up 2429 os << G4endl; // Make sure to line up
2464 os << std::setw(30) << " ExitNormal " << 2430 os << std::setw(30) << " ExitNormal " << " "
2465 << std::setw( 5) << " Valid " << " 2431 << std::setw( 5) << " Valid " << " "
2466 << std::setw( 9) << " Exiting " << " 2432 << std::setw( 9) << " Exiting " << " "
2467 << std::setw( 9) << " Entering" << " 2433 << std::setw( 9) << " Entering" << " "
2468 << std::setw(15) << " Blocked:Volume " < 2434 << std::setw(15) << " Blocked:Volume " << " "
2469 << std::setw( 9) << " ReplicaNo" < 2435 << std::setw( 9) << " ReplicaNo" << " "
2470 << std::setw( 8) << " LastStepZero " < 2436 << std::setw( 8) << " LastStepZero " << " "
2471 << G4endl; 2437 << G4endl;
2472 os << "( " << std::setw(7) << n.fExitNorm 2438 os << "( " << std::setw(7) << n.fExitNormal.x()
2473 << ", " << std::setw(7) << n.fExitNormal. 2439 << ", " << std::setw(7) << n.fExitNormal.y()
2474 << ", " << std::setw(7) << n.fExitNormal. 2440 << ", " << std::setw(7) << n.fExitNormal.z() << " ) "
2475 << std::setw( 5) << n.fValidExitNormal 2441 << std::setw( 5) << n.fValidExitNormal << " "
2476 << std::setw( 9) << n.fExiting 2442 << std::setw( 9) << n.fExiting << " "
2477 << std::setw( 9) << n.fEntering 2443 << std::setw( 9) << n.fEntering << " ";
2478 << 2444 if ( n.fBlockedPhysicalVolume==0 )
2479 if ( n.fBlockedPhysicalVolume==nullptr ) << 2445 { os << std::setw(15) << "None"; }
2480 { os << std::setw(15) << "None"; } <<
2481 else 2446 else
2482 { os << std::setw(15)<< n.fBlockedPhysica << 2447 { os << std::setw(15)<< n.fBlockedPhysicalVolume->GetName(); }
2483 <<
2484 os << std::setw( 9) << n.fBlockedReplica 2448 os << std::setw( 9) << n.fBlockedReplicaNo << " "
2485 << std::setw( 8) << n.fLastStepWasZero 2449 << std::setw( 8) << n.fLastStepWasZero << " "
2486 << G4endl; 2450 << G4endl;
2487 } 2451 }
2488 if( n.fVerbose > 2 ) 2452 if( n.fVerbose > 2 )
2489 { 2453 {
2490 os.precision(8); 2454 os.precision(8);
2491 os << " Current Localpoint = " << n.fLast 2455 os << " Current Localpoint = " << n.fLastLocatedPointLocal << G4endl;
2492 os << " PreviousSftOrigin = " << n.fPrev 2456 os << " PreviousSftOrigin = " << n.fPreviousSftOrigin << G4endl;
2493 os << " PreviousSafety = " << n.fPrev 2457 os << " PreviousSafety = " << n.fPreviousSafety << G4endl;
2494 } 2458 }
2495 if( n.fVerbose > 3 || n.fVerbose == 0 ) 2459 if( n.fVerbose > 3 || n.fVerbose == 0 )
2496 { 2460 {
2497 os << "Current History: " << G4endl << n. 2461 os << "Current History: " << G4endl << n.fHistory;
2498 } 2462 }
2499 2463
2500 os.precision(oldcoutPrec); 2464 os.precision(oldcoutPrec);
2501 return os; 2465 return os;
2502 } 2466 }
2503 2467
2504 //------------------------------------------- 2468 //------------------------------------------------------------------------------
2505 2469
2506 EInside 2470 EInside
2507 G4ITNavigator2::InsideCurrentVolume(const G4T 2471 G4ITNavigator2::InsideCurrentVolume(const G4ThreeVector& globalPoint) const
2508 { 2472 {
2509 const G4AffineTransform& transform = GetGlo 2473 const G4AffineTransform& transform = GetGlobalToLocalTransform();
2510 G4ThreeVector localPoint(transform.Transfor 2474 G4ThreeVector localPoint(transform.TransformPoint(globalPoint));
2511 2475
2512 G4VSolid* solid = fHistory.GetTopVolume()-> 2476 G4VSolid* solid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid();
2513 2477
2514 return solid->Inside(localPoint); 2478 return solid->Inside(localPoint);
2515 } 2479 }
2516 2480
2517 //------------------------------------------- 2481 //------------------------------------------------------------------------------
2518 2482
2519 void G4ITNavigator2::GetRandomInCurrentVolume 2483 void G4ITNavigator2::GetRandomInCurrentVolume(G4ThreeVector& _rndmPoint) const
2520 { 2484 {
2521 const G4AffineTransform& local2Global = Get 2485 const G4AffineTransform& local2Global = GetLocalToGlobalTransform();
2522 G4VSolid* solid = fHistory.GetTopVolume()-> 2486 G4VSolid* solid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid();
2523 2487
2524 G4VoxelLimits voxelLimits; // Defaults to 2488 G4VoxelLimits voxelLimits; // Defaults to "infinite" limits.
2525 G4double vmin, vmax; 2489 G4double vmin, vmax;
2526 G4AffineTransform dummy; 2490 G4AffineTransform dummy;
2527 std::vector<std::vector<G4double> > fExtend << 2491 std::vector<std::vector<double> > fExtend;
2528 2492
2529 solid->CalculateExtent(kXAxis,voxelLimits,d 2493 solid->CalculateExtent(kXAxis,voxelLimits,dummy,vmin,vmax);
2530 fExtend[kXAxis][BoundingBox::kMin] = vmin; 2494 fExtend[kXAxis][BoundingBox::kMin] = vmin;
2531 fExtend[kXAxis][BoundingBox::kMax] = vmax; 2495 fExtend[kXAxis][BoundingBox::kMax] = vmax;
2532 2496
2533 solid->CalculateExtent(kYAxis,voxelLimits,d 2497 solid->CalculateExtent(kYAxis,voxelLimits,dummy,vmin,vmax);
2534 fExtend[kYAxis][BoundingBox::kMin] = vmin; 2498 fExtend[kYAxis][BoundingBox::kMin] = vmin;
2535 fExtend[kYAxis][BoundingBox::kMax] = vmax; 2499 fExtend[kYAxis][BoundingBox::kMax] = vmax;
2536 2500
2537 solid->CalculateExtent(kZAxis,voxelLimits,d 2501 solid->CalculateExtent(kZAxis,voxelLimits,dummy,vmin,vmax);
2538 fExtend[kZAxis][BoundingBox::kMin] = vmin; 2502 fExtend[kZAxis][BoundingBox::kMin] = vmin;
2539 fExtend[kZAxis][BoundingBox::kMax] = vmax; 2503 fExtend[kZAxis][BoundingBox::kMax] = vmax;
2540 2504
2541 G4ThreeVector rndmPos; 2505 G4ThreeVector rndmPos;
2542 2506
2543 while(true) << 2507 while(1)
2544 { 2508 {
2545 for(G4int i = 0 ; i < 3 ; ++i) << 2509 for(size_t i = 0 ; i < 3 ; ++i)
2546 { 2510 {
2547 G4double min = fExtend[i][BoundingBox:: << 2511 double min = fExtend[i][BoundingBox::kMin];
2548 G4double max = fExtend[i][BoundingBox:: << 2512 double max = fExtend[i][BoundingBox::kMax];
2549 rndmPos[i] = G4UniformRand()*(max-min)+ 2513 rndmPos[i] = G4UniformRand()*(max-min)+min;
2550 } 2514 }
2551 2515
2552 if(solid->Inside(rndmPos) == kInside) bre 2516 if(solid->Inside(rndmPos) == kInside) break;
2553 } 2517 }
2554 2518
2555 _rndmPoint = local2Global.TransformPoint(rn 2519 _rndmPoint = local2Global.TransformPoint(rndmPos);
2556 } 2520 }
2557 2521
2558 2522