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