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