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