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