1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
27 | // $Id: G4Navigator.cc,v 1.39 2009/05/08 06:47:32 tnikitin Exp $ |
---|
28 | // GEANT4 tag $ Name: $ |
---|
29 | // |
---|
30 | // class G4Navigator Implementation |
---|
31 | // |
---|
32 | // Original author: Paul Kent, July 95/96 |
---|
33 | // |
---|
34 | // -------------------------------------------------------------------- |
---|
35 | |
---|
36 | #include "G4Navigator.hh" |
---|
37 | #include "G4ios.hh" |
---|
38 | #include <iomanip> |
---|
39 | |
---|
40 | #include "G4GeometryTolerance.hh" |
---|
41 | #include "G4VPhysicalVolume.hh" |
---|
42 | |
---|
43 | // ******************************************************************** |
---|
44 | // Constructor |
---|
45 | // ******************************************************************** |
---|
46 | // |
---|
47 | G4Navigator::G4Navigator() |
---|
48 | : fWasLimitedByGeometry(false), fVerbose(0), |
---|
49 | fTopPhysical(0), fCheck(false), fPushed(false) |
---|
50 | { |
---|
51 | fActive= false; |
---|
52 | ResetStackAndState(); |
---|
53 | |
---|
54 | fActionThreshold_NoZeroSteps = 10; |
---|
55 | fAbandonThreshold_NoZeroSteps = 25; |
---|
56 | |
---|
57 | kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); |
---|
58 | fregularNav.SetNormalNavigation( &fnormalNav ); |
---|
59 | |
---|
60 | fStepEndPoint = G4ThreeVector( kInfinity, kInfinity, kInfinity ); |
---|
61 | } |
---|
62 | |
---|
63 | // ******************************************************************** |
---|
64 | // Destructor |
---|
65 | // ******************************************************************** |
---|
66 | // |
---|
67 | G4Navigator::~G4Navigator() |
---|
68 | {;} |
---|
69 | |
---|
70 | // ******************************************************************** |
---|
71 | // ResetHierarchyAndLocate |
---|
72 | // ******************************************************************** |
---|
73 | // |
---|
74 | G4VPhysicalVolume* |
---|
75 | G4Navigator::ResetHierarchyAndLocate(const G4ThreeVector &p, |
---|
76 | const G4ThreeVector &direction, |
---|
77 | const G4TouchableHistory &h) |
---|
78 | { |
---|
79 | ResetState(); |
---|
80 | fHistory = *h.GetHistory(); |
---|
81 | SetupHierarchy(); |
---|
82 | return LocateGlobalPointAndSetup(p, &direction, true, false); |
---|
83 | } |
---|
84 | |
---|
85 | // ******************************************************************** |
---|
86 | // LocateGlobalPointAndSetup |
---|
87 | // |
---|
88 | // Locate the point in the hierarchy return 0 if outside |
---|
89 | // The direction is required |
---|
90 | // - if on an edge shared by more than two surfaces |
---|
91 | // (to resolve likely looping in tracking) |
---|
92 | // - at initial location of a particle |
---|
93 | // (to resolve potential ambiguity at boundary) |
---|
94 | // |
---|
95 | // Flags on exit: (comments to be completed) |
---|
96 | // fEntering - True if entering `daughter' volume (or replica) |
---|
97 | // whether daughter of last mother directly |
---|
98 | // or daughter of that volume's ancestor. |
---|
99 | // ******************************************************************** |
---|
100 | // |
---|
101 | G4VPhysicalVolume* |
---|
102 | G4Navigator::LocateGlobalPointAndSetup( const G4ThreeVector& globalPoint, |
---|
103 | const G4ThreeVector* pGlobalDirection, |
---|
104 | const G4bool relativeSearch, |
---|
105 | const G4bool ignoreDirection ) |
---|
106 | { |
---|
107 | G4bool notKnownContained=true, noResult; |
---|
108 | G4VPhysicalVolume *targetPhysical; |
---|
109 | G4LogicalVolume *targetLogical; |
---|
110 | G4VSolid *targetSolid=0; |
---|
111 | G4ThreeVector localPoint, globalDirection; |
---|
112 | EInside insideCode; |
---|
113 | |
---|
114 | G4bool considerDirection = (!ignoreDirection) || fLocatedOnEdge; |
---|
115 | |
---|
116 | if( considerDirection && pGlobalDirection != 0 ) |
---|
117 | { |
---|
118 | globalDirection=*pGlobalDirection; |
---|
119 | } |
---|
120 | |
---|
121 | #ifdef G4DEBUG_NAVIGATION |
---|
122 | if( fVerbose > 2 ) |
---|
123 | { |
---|
124 | G4cout << "Upon entering LocateGlobalPointAndSetup():" << G4endl; |
---|
125 | G4cout << " History = " << G4endl << fHistory << G4endl << G4endl; |
---|
126 | } |
---|
127 | #endif |
---|
128 | |
---|
129 | #ifdef G4VERBOSE |
---|
130 | G4int oldcoutPrec = G4cout.precision(8); |
---|
131 | if( fVerbose > 2 ) |
---|
132 | { |
---|
133 | G4cout << "*** G4Navigator::LocateGlobalPointAndSetup: ***" << G4endl; |
---|
134 | G4cout << " Called with arguments: " << G4endl |
---|
135 | << " Globalpoint = " << globalPoint << G4endl |
---|
136 | << " RelativeSearch = " << relativeSearch << G4endl; |
---|
137 | if( fVerbose == 4 ) |
---|
138 | { |
---|
139 | G4cout << " ----- Upon entering:" << G4endl; |
---|
140 | PrintState(); |
---|
141 | } |
---|
142 | } |
---|
143 | #endif |
---|
144 | |
---|
145 | if ( !relativeSearch ) |
---|
146 | { |
---|
147 | ResetStackAndState(); |
---|
148 | } |
---|
149 | else |
---|
150 | { |
---|
151 | if ( fWasLimitedByGeometry ) |
---|
152 | { |
---|
153 | fWasLimitedByGeometry = false; |
---|
154 | fEnteredDaughter = fEntering; // Remember |
---|
155 | fExitedMother = fExiting; // Remember |
---|
156 | if ( fExiting ) |
---|
157 | { |
---|
158 | if ( fHistory.GetDepth() ) |
---|
159 | { |
---|
160 | fBlockedPhysicalVolume = fHistory.GetTopVolume(); |
---|
161 | fBlockedReplicaNo = fHistory.GetTopReplicaNo(); |
---|
162 | fHistory.BackLevel(); |
---|
163 | } |
---|
164 | else |
---|
165 | { |
---|
166 | fLastLocatedPointLocal = localPoint; |
---|
167 | fLocatedOutsideWorld = true; |
---|
168 | return 0; // Have exited world volume |
---|
169 | } |
---|
170 | // A fix for the case where a volume is "entered" at an edge |
---|
171 | // and a coincident surface exists outside it. |
---|
172 | // - This stops it from exiting further volumes and cycling |
---|
173 | // - However ReplicaNavigator treats this case itself |
---|
174 | // |
---|
175 | if ( fLocatedOnEdge && (VolumeType(fBlockedPhysicalVolume)!=kReplica )) |
---|
176 | { |
---|
177 | fExiting= false; |
---|
178 | } |
---|
179 | } |
---|
180 | else |
---|
181 | if ( fEntering ) |
---|
182 | { |
---|
183 | switch (VolumeType(fBlockedPhysicalVolume)) |
---|
184 | { |
---|
185 | case kNormal: |
---|
186 | fHistory.NewLevel(fBlockedPhysicalVolume, kNormal, |
---|
187 | fBlockedPhysicalVolume->GetCopyNo()); |
---|
188 | break; |
---|
189 | case kReplica: |
---|
190 | freplicaNav.ComputeTransformation(fBlockedReplicaNo, |
---|
191 | fBlockedPhysicalVolume); |
---|
192 | fHistory.NewLevel(fBlockedPhysicalVolume, kReplica, |
---|
193 | fBlockedReplicaNo); |
---|
194 | fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo); |
---|
195 | break; |
---|
196 | case kParameterised: |
---|
197 | if( fBlockedPhysicalVolume->GetRegularStructureId() != 1 ) |
---|
198 | { |
---|
199 | G4VSolid *pSolid; |
---|
200 | G4VPVParameterisation *pParam; |
---|
201 | G4TouchableHistory parentTouchable( fHistory ); |
---|
202 | pParam = fBlockedPhysicalVolume->GetParameterisation(); |
---|
203 | pSolid = pParam->ComputeSolid(fBlockedReplicaNo, |
---|
204 | fBlockedPhysicalVolume); |
---|
205 | pSolid->ComputeDimensions(pParam, fBlockedReplicaNo, |
---|
206 | fBlockedPhysicalVolume); |
---|
207 | pParam->ComputeTransformation(fBlockedReplicaNo, |
---|
208 | fBlockedPhysicalVolume); |
---|
209 | fHistory.NewLevel(fBlockedPhysicalVolume, kParameterised, |
---|
210 | fBlockedReplicaNo); |
---|
211 | fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo); |
---|
212 | // |
---|
213 | // Set the correct solid and material in Logical Volume |
---|
214 | // |
---|
215 | G4LogicalVolume *pLogical; |
---|
216 | pLogical = fBlockedPhysicalVolume->GetLogicalVolume(); |
---|
217 | pLogical->SetSolid( pSolid ); |
---|
218 | pLogical->UpdateMaterial(pParam -> |
---|
219 | ComputeMaterial(fBlockedReplicaNo, |
---|
220 | fBlockedPhysicalVolume, |
---|
221 | &parentTouchable)); |
---|
222 | } |
---|
223 | break; |
---|
224 | } |
---|
225 | fEntering = false; |
---|
226 | fBlockedPhysicalVolume = 0; |
---|
227 | localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint); |
---|
228 | notKnownContained = false; |
---|
229 | } |
---|
230 | } |
---|
231 | else |
---|
232 | { |
---|
233 | fBlockedPhysicalVolume = 0; |
---|
234 | fEntering = false; |
---|
235 | fEnteredDaughter = false; // Full Step was not taken, did not enter |
---|
236 | fExiting = false; |
---|
237 | fExitedMother = false; // Full Step was not taken, did not exit |
---|
238 | } |
---|
239 | } |
---|
240 | // |
---|
241 | // Search from top of history up through geometry until |
---|
242 | // containing volume found: |
---|
243 | // If on |
---|
244 | // o OUTSIDE - Back up level, not/no longer exiting volumes |
---|
245 | // o SURFACE and EXITING - Back up level, setting new blocking no.s |
---|
246 | // else |
---|
247 | // o containing volume found |
---|
248 | // |
---|
249 | while (notKnownContained) |
---|
250 | { |
---|
251 | if ( fHistory.GetTopVolumeType()!=kReplica ) |
---|
252 | { |
---|
253 | targetSolid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid(); |
---|
254 | localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint); |
---|
255 | insideCode = targetSolid->Inside(localPoint); |
---|
256 | #ifdef G4VERBOSE |
---|
257 | if(( fVerbose == 1 ) && ( fCheck )) |
---|
258 | { |
---|
259 | G4String solidResponse = "-kInside-"; |
---|
260 | if (insideCode == kOutside) |
---|
261 | solidResponse = "-kOutside-"; |
---|
262 | else if (insideCode == kSurface) |
---|
263 | solidResponse = "-kSurface-"; |
---|
264 | G4cout << "*** G4Navigator::LocateGlobalPointAndSetup(): ***" << G4endl |
---|
265 | << " Invoked Inside() for solid: " << targetSolid->GetName() |
---|
266 | << ". Solid replied: " << solidResponse << G4endl |
---|
267 | << " For local point p: " << localPoint << G4endl; |
---|
268 | } |
---|
269 | #endif |
---|
270 | } |
---|
271 | else |
---|
272 | { |
---|
273 | insideCode = freplicaNav.BackLocate(fHistory, globalPoint, localPoint, |
---|
274 | fExiting, notKnownContained); |
---|
275 | // !CARE! if notKnownContained returns false then the point is within |
---|
276 | // the containing placement volume of the replica(s). If insidecode |
---|
277 | // will result in the history being backed up one level, then the |
---|
278 | // local point returned is the point in the system of this new level |
---|
279 | } |
---|
280 | if ( insideCode==kOutside ) |
---|
281 | { |
---|
282 | if ( fHistory.GetDepth() ) |
---|
283 | { |
---|
284 | fBlockedPhysicalVolume = fHistory.GetTopVolume(); |
---|
285 | fBlockedReplicaNo = fHistory.GetTopReplicaNo(); |
---|
286 | fHistory.BackLevel(); |
---|
287 | fExiting = false; |
---|
288 | } |
---|
289 | else |
---|
290 | { |
---|
291 | fLastLocatedPointLocal = localPoint; |
---|
292 | fLocatedOutsideWorld = true; |
---|
293 | return 0; // Have exited world volume |
---|
294 | } |
---|
295 | } |
---|
296 | else |
---|
297 | if ( insideCode==kSurface ) |
---|
298 | { |
---|
299 | G4bool isExiting = fExiting; |
---|
300 | if( (!fExiting)&&considerDirection ) |
---|
301 | { |
---|
302 | // Figure out whether we are exiting this level's volume |
---|
303 | // by using the direction |
---|
304 | // |
---|
305 | G4bool directionExiting = false; |
---|
306 | G4ThreeVector localDirection = |
---|
307 | fHistory.GetTopTransform().TransformAxis(globalDirection); |
---|
308 | if ( fHistory.GetTopVolumeType()!=kReplica ) |
---|
309 | { |
---|
310 | G4ThreeVector normal = targetSolid->SurfaceNormal(localPoint); |
---|
311 | directionExiting = normal.dot(localDirection) > 0.0; |
---|
312 | isExiting = isExiting || directionExiting; |
---|
313 | } |
---|
314 | } |
---|
315 | if( isExiting ) |
---|
316 | { |
---|
317 | if ( fHistory.GetDepth() ) |
---|
318 | { |
---|
319 | fBlockedPhysicalVolume = fHistory.GetTopVolume(); |
---|
320 | fBlockedReplicaNo = fHistory.GetTopReplicaNo(); |
---|
321 | fHistory.BackLevel(); |
---|
322 | // |
---|
323 | // Still on surface but exited volume not necessarily convex |
---|
324 | // |
---|
325 | fValidExitNormal = false; |
---|
326 | } |
---|
327 | else |
---|
328 | { |
---|
329 | fLastLocatedPointLocal = localPoint; |
---|
330 | fLocatedOutsideWorld = true; |
---|
331 | return 0; // Have exited world volume |
---|
332 | } |
---|
333 | } |
---|
334 | else |
---|
335 | { |
---|
336 | notKnownContained=false; |
---|
337 | } |
---|
338 | } |
---|
339 | else |
---|
340 | { |
---|
341 | notKnownContained=false; |
---|
342 | } |
---|
343 | } // END while (notKnownContained) |
---|
344 | // |
---|
345 | // Search downwards until deepest containing volume found, |
---|
346 | // blocking fBlockedPhysicalVolume/BlockedReplicaNum |
---|
347 | // |
---|
348 | // 3 Cases: |
---|
349 | // |
---|
350 | // o Parameterised daughters |
---|
351 | // =>Must be one G4PVParameterised daughter & voxels |
---|
352 | // o Positioned daughters & voxels |
---|
353 | // o Positioned daughters & no voxels |
---|
354 | |
---|
355 | noResult = true; // noResult should be renamed to |
---|
356 | // something like enteredLevel, as that is its meaning. |
---|
357 | do |
---|
358 | { |
---|
359 | // Determine `type' of current mother volume |
---|
360 | // |
---|
361 | targetPhysical = fHistory.GetTopVolume(); |
---|
362 | targetLogical = targetPhysical->GetLogicalVolume(); |
---|
363 | switch( CharacteriseDaughters(targetLogical) ) |
---|
364 | { |
---|
365 | case kNormal: |
---|
366 | if ( targetLogical->GetVoxelHeader() ) // use optimised navigation |
---|
367 | { |
---|
368 | noResult = fvoxelNav.LevelLocate(fHistory, |
---|
369 | fBlockedPhysicalVolume, |
---|
370 | fBlockedReplicaNo, |
---|
371 | globalPoint, |
---|
372 | pGlobalDirection, |
---|
373 | considerDirection, |
---|
374 | localPoint); |
---|
375 | } |
---|
376 | else // do not use optimised navigation |
---|
377 | { |
---|
378 | noResult = fnormalNav.LevelLocate(fHistory, |
---|
379 | fBlockedPhysicalVolume, |
---|
380 | fBlockedReplicaNo, |
---|
381 | globalPoint, |
---|
382 | pGlobalDirection, |
---|
383 | considerDirection, |
---|
384 | localPoint); |
---|
385 | } |
---|
386 | break; |
---|
387 | case kReplica: |
---|
388 | noResult = freplicaNav.LevelLocate(fHistory, |
---|
389 | fBlockedPhysicalVolume, |
---|
390 | fBlockedReplicaNo, |
---|
391 | globalPoint, |
---|
392 | pGlobalDirection, |
---|
393 | considerDirection, |
---|
394 | localPoint); |
---|
395 | break; |
---|
396 | case kParameterised: |
---|
397 | if( GetDaughtersRegularStructureId(targetLogical) != 1 ) |
---|
398 | { |
---|
399 | noResult = fparamNav.LevelLocate(fHistory, |
---|
400 | fBlockedPhysicalVolume, |
---|
401 | fBlockedReplicaNo, |
---|
402 | globalPoint, |
---|
403 | pGlobalDirection, |
---|
404 | considerDirection, |
---|
405 | localPoint); |
---|
406 | } |
---|
407 | else // Regular structure |
---|
408 | { |
---|
409 | noResult = fregularNav.LevelLocate(fHistory, |
---|
410 | fBlockedPhysicalVolume, |
---|
411 | fBlockedReplicaNo, |
---|
412 | globalPoint, |
---|
413 | pGlobalDirection, |
---|
414 | considerDirection, |
---|
415 | localPoint); |
---|
416 | } |
---|
417 | break; |
---|
418 | } |
---|
419 | |
---|
420 | // LevelLocate returns true if it finds a daughter volume |
---|
421 | // in which globalPoint is inside (or on the surface). |
---|
422 | |
---|
423 | if ( noResult ) |
---|
424 | { |
---|
425 | // Entering a daughter after ascending |
---|
426 | // |
---|
427 | // The blocked volume is no longer valid - it was for another level |
---|
428 | // |
---|
429 | fBlockedPhysicalVolume = 0; |
---|
430 | fBlockedReplicaNo = -1; |
---|
431 | |
---|
432 | // fEntering should be false -- else blockedVolume is assumed good. |
---|
433 | // fEnteredDaughter is used for ExitNormal |
---|
434 | // |
---|
435 | fEntering = false; |
---|
436 | fEnteredDaughter = true; |
---|
437 | #ifdef G4DEBUG_NAVIGATION |
---|
438 | if( fVerbose > 2 ) |
---|
439 | { |
---|
440 | G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume(); |
---|
441 | G4cout << "*** G4Navigator::LocateGlobalPointAndSetup() ***" << G4endl; |
---|
442 | G4cout << " Entering volume: " << enteredPhysical->GetName() |
---|
443 | << G4endl; |
---|
444 | } |
---|
445 | #endif |
---|
446 | } |
---|
447 | } while (noResult); |
---|
448 | |
---|
449 | fLastLocatedPointLocal = localPoint; |
---|
450 | |
---|
451 | #ifdef G4VERBOSE |
---|
452 | if( fVerbose == 4 ) |
---|
453 | { |
---|
454 | G4cout.precision(6); |
---|
455 | G4String curPhysVol_Name("None"); |
---|
456 | if (targetPhysical!=0) |
---|
457 | { |
---|
458 | curPhysVol_Name = targetPhysical->GetName(); |
---|
459 | } |
---|
460 | G4cout << " Return value = new volume = " << curPhysVol_Name << G4endl; |
---|
461 | G4cout << " ----- Upon exiting:" << G4endl; |
---|
462 | PrintState(); |
---|
463 | #ifdef G4DEBUG_NAVIGATION |
---|
464 | G4cout << "Upon exiting LocateGlobalPointAndSetup():" << G4endl; |
---|
465 | G4cout << " History = " << G4endl << fHistory << G4endl << G4endl; |
---|
466 | #endif |
---|
467 | } |
---|
468 | G4cout.precision(oldcoutPrec); |
---|
469 | #endif |
---|
470 | |
---|
471 | fLocatedOutsideWorld= false; |
---|
472 | |
---|
473 | return targetPhysical; |
---|
474 | } |
---|
475 | |
---|
476 | // ******************************************************************** |
---|
477 | // LocateGlobalPointWithinVolume |
---|
478 | // |
---|
479 | // -> the state information of this Navigator and its subNavigators |
---|
480 | // is updated in order to start the next step at pGlobalpoint |
---|
481 | // -> no check is performed whether pGlobalpoint is inside the |
---|
482 | // original volume (this must be the case). |
---|
483 | // |
---|
484 | // Note: a direction could be added to the arguments, to aid in future |
---|
485 | // optional checking (via the old code below, flagged by OLD_LOCATE). |
---|
486 | // [ This would be done only in verbose mode ] |
---|
487 | // ******************************************************************** |
---|
488 | // |
---|
489 | void |
---|
490 | G4Navigator::LocateGlobalPointWithinVolume(const G4ThreeVector& pGlobalpoint) |
---|
491 | { |
---|
492 | fLastLocatedPointLocal = ComputeLocalPoint(pGlobalpoint); |
---|
493 | |
---|
494 | #ifdef G4DEBUG_NAVIGATION |
---|
495 | if( fVerbose > 2 ) |
---|
496 | { |
---|
497 | G4cout << "Entering LocateGlobalWithinVolume(): History = " << G4endl; |
---|
498 | G4cout << fHistory << G4endl; |
---|
499 | } |
---|
500 | #endif |
---|
501 | |
---|
502 | // For the case of Voxel (or Parameterised) volume the respective |
---|
503 | // Navigator must be messaged to update its voxel information etc |
---|
504 | |
---|
505 | // Update the state of the Sub Navigators |
---|
506 | // - in particular any voxel information they store/cache |
---|
507 | // |
---|
508 | G4VPhysicalVolume* motherPhysical = fHistory.GetTopVolume(); |
---|
509 | G4LogicalVolume* motherLogical = motherPhysical->GetLogicalVolume(); |
---|
510 | G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader(); |
---|
511 | |
---|
512 | if ( fHistory.GetTopVolumeType()!=kReplica ) |
---|
513 | { |
---|
514 | switch( CharacteriseDaughters(motherLogical) ) |
---|
515 | { |
---|
516 | case kNormal: |
---|
517 | if ( pVoxelHeader ) |
---|
518 | { |
---|
519 | fvoxelNav.VoxelLocate( pVoxelHeader, fLastLocatedPointLocal ); |
---|
520 | } |
---|
521 | break; |
---|
522 | case kParameterised: |
---|
523 | if( GetDaughtersRegularStructureId(motherLogical) != 1 ) |
---|
524 | { |
---|
525 | // Resets state & returns voxel node |
---|
526 | // |
---|
527 | fparamNav.ParamVoxelLocate( pVoxelHeader, fLastLocatedPointLocal ); |
---|
528 | } |
---|
529 | break; |
---|
530 | case kReplica: |
---|
531 | G4Exception("G4Navigator::LocateGlobalPointWithinVolume()", |
---|
532 | "NotApplicable", FatalException, |
---|
533 | "Not applicable for replicated volumes."); |
---|
534 | break; |
---|
535 | } |
---|
536 | } |
---|
537 | |
---|
538 | // Reset the state variables |
---|
539 | // - which would have been affected |
---|
540 | // by the 'equivalent' call to LocateGlobalPointAndSetup |
---|
541 | // - who's values have been invalidated by the 'move'. |
---|
542 | // |
---|
543 | fBlockedPhysicalVolume = 0; |
---|
544 | fBlockedReplicaNo = -1; |
---|
545 | fEntering = false; |
---|
546 | fEnteredDaughter = false; // Boundary not encountered, did not enter |
---|
547 | fExiting = false; |
---|
548 | fExitedMother = false; // Boundary not encountered, did not exit |
---|
549 | } |
---|
550 | |
---|
551 | // ******************************************************************** |
---|
552 | // SetSavedState |
---|
553 | // |
---|
554 | // Save the state, in case this is a parasitic call |
---|
555 | // Save fValidExitNormal, fExitNormal, fExiting, fEntering, |
---|
556 | // fBlockedPhysicalVolume, fBlockedReplicaNo, fLastStepWasZero; |
---|
557 | // ******************************************************************** |
---|
558 | // |
---|
559 | void G4Navigator::SetSavedState() |
---|
560 | { |
---|
561 | // fSaveExitNormal = fExitNormal; |
---|
562 | fSaveState.sExitNormal = fExitNormal; |
---|
563 | fSaveState.sValidExitNormal = fValidExitNormal; |
---|
564 | fSaveState.sExiting = fExiting; |
---|
565 | fSaveState.sEntering = fEntering; |
---|
566 | |
---|
567 | fSaveState.spBlockedPhysicalVolume = fBlockedPhysicalVolume; |
---|
568 | fSaveState.sBlockedReplicaNo = fBlockedReplicaNo, |
---|
569 | |
---|
570 | fSaveState.sLastStepWasZero = fLastStepWasZero; |
---|
571 | } |
---|
572 | |
---|
573 | // ******************************************************************** |
---|
574 | // RestoreSavedState |
---|
575 | // |
---|
576 | // Restore the state (in Compute Step), in case this is a parasitic call |
---|
577 | // ******************************************************************** |
---|
578 | // |
---|
579 | void G4Navigator::RestoreSavedState() |
---|
580 | { |
---|
581 | fExitNormal = fSaveState.sExitNormal; |
---|
582 | fValidExitNormal = fSaveState.sValidExitNormal; |
---|
583 | fExiting = fSaveState.sExiting; |
---|
584 | fEntering = fSaveState.sEntering; |
---|
585 | |
---|
586 | fBlockedPhysicalVolume = fSaveState.spBlockedPhysicalVolume; |
---|
587 | fBlockedReplicaNo = fSaveState.sBlockedReplicaNo, |
---|
588 | |
---|
589 | fLastStepWasZero = fSaveState.sLastStepWasZero; |
---|
590 | } |
---|
591 | |
---|
592 | // ******************************************************************** |
---|
593 | // ComputeStep |
---|
594 | // |
---|
595 | // Computes the next geometric Step: intersections with current |
---|
596 | // mother and `daughter' volumes. |
---|
597 | // |
---|
598 | // NOTE: |
---|
599 | // |
---|
600 | // Flags on entry: |
---|
601 | // -------------- |
---|
602 | // fValidExitNormal - Normal of exited volume is valid (convex, not a |
---|
603 | // coincident boundary) |
---|
604 | // fExitNormal - Surface normal of exited volume |
---|
605 | // fExiting - True if have exited solid |
---|
606 | // |
---|
607 | // fBlockedPhysicalVolume - Ptr to exited volume (or 0) |
---|
608 | // fBlockedReplicaNo - Replication no of exited volume |
---|
609 | // fLastStepWasZero - True if last Step size was zero. |
---|
610 | // |
---|
611 | // Flags on exit: |
---|
612 | // ------------- |
---|
613 | // fValidExitNormal - True if surface normal of exited volume is valid |
---|
614 | // fExitNormal - Surface normal of exited volume rotated to mothers |
---|
615 | // reference system |
---|
616 | // fExiting - True if exiting mother |
---|
617 | // fEntering - True if entering `daughter' volume (or replica) |
---|
618 | // fBlockedPhysicalVolume - Ptr to candidate (entered) volume |
---|
619 | // fBlockedReplicaNo - Replication no of candidate (entered) volume |
---|
620 | // fLastStepWasZero - True if this Step size was zero. |
---|
621 | // ******************************************************************** |
---|
622 | // |
---|
623 | G4double G4Navigator::ComputeStep( const G4ThreeVector &pGlobalpoint, |
---|
624 | const G4ThreeVector &pDirection, |
---|
625 | const G4double pCurrentProposedStepLength, |
---|
626 | G4double &pNewSafety) |
---|
627 | { |
---|
628 | G4ThreeVector localDirection = ComputeLocalAxis(pDirection); |
---|
629 | G4double Step = DBL_MAX; |
---|
630 | G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume(); |
---|
631 | G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume(); |
---|
632 | |
---|
633 | static G4int sNavCScalls=0; |
---|
634 | sNavCScalls++; |
---|
635 | |
---|
636 | #ifdef G4VERBOSE |
---|
637 | G4int oldcoutPrec= G4cout.precision(8); |
---|
638 | G4int oldcerrPrec= G4cerr.precision(10); |
---|
639 | if( fVerbose > 0 ) |
---|
640 | { |
---|
641 | G4cout << "*** G4Navigator::ComputeStep: ***" << G4endl; |
---|
642 | G4cout << " Volume = " << motherPhysical->GetName() |
---|
643 | << " - Proposed step length = " << pCurrentProposedStepLength |
---|
644 | << G4endl; |
---|
645 | #ifdef G4DEBUG_NAVIGATION |
---|
646 | if( fVerbose >= 4 ) |
---|
647 | { |
---|
648 | G4cout << " Called with the arguments: " << G4endl |
---|
649 | << " Globalpoint = " << std::setw(25) << pGlobalpoint << G4endl |
---|
650 | << " Direction = " << std::setw(25) << pDirection << G4endl; |
---|
651 | G4cout << " ---- Upon entering :" << G4endl; |
---|
652 | PrintState(); |
---|
653 | } |
---|
654 | #endif |
---|
655 | } |
---|
656 | |
---|
657 | static G4double fAccuracyForWarning = kCarTolerance, |
---|
658 | fAccuracyForException = 1000*kCarTolerance; |
---|
659 | #endif |
---|
660 | |
---|
661 | G4ThreeVector newLocalPoint = ComputeLocalPoint(pGlobalpoint); |
---|
662 | if( newLocalPoint != fLastLocatedPointLocal ) |
---|
663 | { |
---|
664 | // Check whether the relocation is within safety |
---|
665 | // |
---|
666 | G4ThreeVector oldLocalPoint = fLastLocatedPointLocal; |
---|
667 | G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2(); |
---|
668 | |
---|
669 | if ( moveLenSq >= kCarTolerance*kCarTolerance ) |
---|
670 | { |
---|
671 | #ifdef G4VERBOSE |
---|
672 | // The following checks only make sense if the move is larger |
---|
673 | // than the tolerance. |
---|
674 | // |
---|
675 | G4ThreeVector OriginalGlobalpoint = |
---|
676 | fHistory.GetTopTransform().Inverse(). |
---|
677 | TransformPoint(fLastLocatedPointLocal); |
---|
678 | |
---|
679 | G4double shiftOriginSafSq = (fPreviousSftOrigin-pGlobalpoint).mag2(); |
---|
680 | |
---|
681 | // Check that the starting point of this step is |
---|
682 | // within the isotropic safety sphere of the last point |
---|
683 | // to a accuracy/precision given by fAccuracyForWarning. |
---|
684 | // If so give warning. |
---|
685 | // If it fails by more than fAccuracyForException exit with error. |
---|
686 | // |
---|
687 | if( shiftOriginSafSq >= sqr(fPreviousSafety) ) |
---|
688 | { |
---|
689 | G4double shiftOrigin = std::sqrt(shiftOriginSafSq); |
---|
690 | G4double diffShiftSaf = shiftOrigin - fPreviousSafety; |
---|
691 | |
---|
692 | if( diffShiftSaf > fAccuracyForWarning ) |
---|
693 | { |
---|
694 | G4Exception("G4Navigator::ComputeStep()", |
---|
695 | "UnexpectedPositionShift", JustWarning, |
---|
696 | "Accuracy ERROR or slightly inaccurate position shift."); |
---|
697 | G4cerr << " The Step's starting point has moved " |
---|
698 | << std::sqrt(moveLenSq)/mm << " mm " << G4endl |
---|
699 | << " since the last call to a Locate method." << G4endl; |
---|
700 | G4cerr << " This has resulted in moving " |
---|
701 | << shiftOrigin/mm << " mm " |
---|
702 | << " from the last point at which the safety " |
---|
703 | << " was calculated " << G4endl; |
---|
704 | G4cerr << " which is more than the computed safety= " |
---|
705 | << fPreviousSafety/mm << " mm at that point." << G4endl; |
---|
706 | G4cerr << " This difference is " |
---|
707 | << diffShiftSaf/mm << " mm." << G4endl |
---|
708 | << " The tolerated accuracy is " |
---|
709 | << fAccuracyForException/mm << " mm." << G4endl; |
---|
710 | |
---|
711 | static G4int warnNow = 0; |
---|
712 | if( ((++warnNow % 100) == 1) ) |
---|
713 | { |
---|
714 | G4cerr << " This problem can be due to either " << G4endl; |
---|
715 | G4cerr << " - a process that has proposed a displacement" |
---|
716 | << " larger than the current safety , or" << G4endl; |
---|
717 | G4cerr << " - inaccuracy in the computation of the safety" |
---|
718 | << G4endl; |
---|
719 | G4cerr << " We suggest that you " << G4endl |
---|
720 | << " - find i) what particle is being tracked, and " |
---|
721 | << " ii) through what part of your geometry " << G4endl |
---|
722 | << " for example by re-running this event with " |
---|
723 | << G4endl |
---|
724 | << " /tracking/verbose 1 " << G4endl |
---|
725 | << " - check which processes you declare for" |
---|
726 | << " this particle (and look at non-standard ones)" |
---|
727 | << G4endl |
---|
728 | << " - in case, create a detailed logfile" |
---|
729 | << " of this event using:" << G4endl |
---|
730 | << " /tracking/verbose 6 " |
---|
731 | << G4endl; |
---|
732 | } |
---|
733 | } |
---|
734 | #ifdef G4DEBUG_NAVIGATION |
---|
735 | else |
---|
736 | { |
---|
737 | G4cerr << "WARNING - G4Navigator::ComputeStep()" << G4endl |
---|
738 | << " The Step's starting point has moved " |
---|
739 | << std::sqrt(moveLenSq) << "," << G4endl |
---|
740 | << " which has taken it to the limit of" |
---|
741 | << " the current safety. " << G4endl; |
---|
742 | } |
---|
743 | #endif |
---|
744 | } |
---|
745 | G4double safetyPlus = fPreviousSafety + fAccuracyForException; |
---|
746 | if ( shiftOriginSafSq > sqr(safetyPlus) ) |
---|
747 | { |
---|
748 | G4cerr << "ERROR - G4Navigator::ComputeStep()" << G4endl |
---|
749 | << " Position has shifted considerably without" |
---|
750 | << " notifying the navigator !" << G4endl |
---|
751 | << " Tolerated safety: " << safetyPlus << G4endl |
---|
752 | << " Computed shift : " << shiftOriginSafSq << G4endl; |
---|
753 | G4Exception("G4Navigator::ComputeStep()", |
---|
754 | "SignificantPositionShift", JustWarning, |
---|
755 | "May lead to a crash or unreliable results."); |
---|
756 | } |
---|
757 | #endif // end G4VERBOSE |
---|
758 | |
---|
759 | // Relocate the point within the same volume |
---|
760 | // |
---|
761 | LocateGlobalPointWithinVolume( pGlobalpoint ); |
---|
762 | } |
---|
763 | } |
---|
764 | if ( fHistory.GetTopVolumeType()!=kReplica ) |
---|
765 | { |
---|
766 | switch( CharacteriseDaughters(motherLogical) ) |
---|
767 | { |
---|
768 | case kNormal: |
---|
769 | if ( motherLogical->GetVoxelHeader() ) |
---|
770 | { |
---|
771 | Step = fvoxelNav.ComputeStep(fLastLocatedPointLocal, |
---|
772 | localDirection, |
---|
773 | pCurrentProposedStepLength, |
---|
774 | pNewSafety, |
---|
775 | fHistory, |
---|
776 | fValidExitNormal, |
---|
777 | fExitNormal, |
---|
778 | fExiting, |
---|
779 | fEntering, |
---|
780 | &fBlockedPhysicalVolume, |
---|
781 | fBlockedReplicaNo); |
---|
782 | |
---|
783 | } |
---|
784 | else |
---|
785 | { |
---|
786 | if( motherPhysical->GetRegularStructureId() != 1 ) |
---|
787 | { |
---|
788 | Step = fnormalNav.ComputeStep(fLastLocatedPointLocal, |
---|
789 | localDirection, |
---|
790 | pCurrentProposedStepLength, |
---|
791 | pNewSafety, |
---|
792 | fHistory, |
---|
793 | fValidExitNormal, |
---|
794 | fExitNormal, |
---|
795 | fExiting, |
---|
796 | fEntering, |
---|
797 | &fBlockedPhysicalVolume, |
---|
798 | fBlockedReplicaNo); |
---|
799 | } |
---|
800 | else // Regular (non-voxelised) structure |
---|
801 | { |
---|
802 | LocateGlobalPointAndSetup( pGlobalpoint, &pDirection, true, true ); |
---|
803 | // |
---|
804 | // if physical process limits the step, the voxel will not be the |
---|
805 | // one given by ComputeStepSkippingEqualMaterials() and the local |
---|
806 | // point will be wrongly calculated. |
---|
807 | |
---|
808 | // There is a problem: when msc limits the step and the point is |
---|
809 | // assigned wrongly to phantom in previous step (while it is out |
---|
810 | // of the container volume). Then LocateGlobalPointAndSetup() has |
---|
811 | // reset the history topvolume to world. |
---|
812 | // |
---|
813 | if(fHistory.GetTopVolume()->GetRegularStructureId() != 1 ) |
---|
814 | { |
---|
815 | G4Exception("G4Navigator::ComputeStep()", |
---|
816 | "Bad-location-of-point", JustWarning, |
---|
817 | "Point is relocated in voxels, while it should be outside!"); |
---|
818 | Step = fnormalNav.ComputeStep(fLastLocatedPointLocal, |
---|
819 | localDirection, |
---|
820 | pCurrentProposedStepLength, |
---|
821 | pNewSafety, |
---|
822 | fHistory, |
---|
823 | fValidExitNormal, |
---|
824 | fExitNormal, |
---|
825 | fExiting, |
---|
826 | fEntering, |
---|
827 | &fBlockedPhysicalVolume, |
---|
828 | fBlockedReplicaNo); |
---|
829 | } |
---|
830 | else |
---|
831 | { |
---|
832 | Step = fregularNav. |
---|
833 | ComputeStepSkippingEqualMaterials(fLastLocatedPointLocal, |
---|
834 | localDirection, |
---|
835 | pCurrentProposedStepLength, |
---|
836 | pNewSafety, |
---|
837 | fHistory, |
---|
838 | fValidExitNormal, |
---|
839 | fExitNormal, |
---|
840 | fExiting, |
---|
841 | fEntering, |
---|
842 | &fBlockedPhysicalVolume, |
---|
843 | fBlockedReplicaNo, |
---|
844 | motherPhysical); |
---|
845 | } |
---|
846 | } |
---|
847 | } |
---|
848 | break; |
---|
849 | case kParameterised: |
---|
850 | if( GetDaughtersRegularStructureId(motherLogical) != 1 ) |
---|
851 | { |
---|
852 | Step = fparamNav.ComputeStep(fLastLocatedPointLocal, |
---|
853 | localDirection, |
---|
854 | pCurrentProposedStepLength, |
---|
855 | pNewSafety, |
---|
856 | fHistory, |
---|
857 | fValidExitNormal, |
---|
858 | fExitNormal, |
---|
859 | fExiting, |
---|
860 | fEntering, |
---|
861 | &fBlockedPhysicalVolume, |
---|
862 | fBlockedReplicaNo); |
---|
863 | } |
---|
864 | else // Regular structure |
---|
865 | { |
---|
866 | Step = fregularNav.ComputeStep(fLastLocatedPointLocal, |
---|
867 | localDirection, |
---|
868 | pCurrentProposedStepLength, |
---|
869 | pNewSafety, |
---|
870 | fHistory, |
---|
871 | fValidExitNormal, |
---|
872 | fExitNormal, |
---|
873 | fExiting, |
---|
874 | fEntering, |
---|
875 | &fBlockedPhysicalVolume, |
---|
876 | fBlockedReplicaNo); |
---|
877 | } |
---|
878 | break; |
---|
879 | case kReplica: |
---|
880 | G4Exception("G4Navigator::ComputeStep()", "NotApplicable", |
---|
881 | FatalException, "Not applicable for replicated volumes."); |
---|
882 | break; |
---|
883 | } |
---|
884 | } |
---|
885 | else |
---|
886 | { |
---|
887 | // In the case of a replica, it must handle the exiting |
---|
888 | // edge/corner problem by itself |
---|
889 | // |
---|
890 | G4bool exitingReplica = fExitedMother; |
---|
891 | Step = freplicaNav.ComputeStep(pGlobalpoint, |
---|
892 | pDirection, |
---|
893 | fLastLocatedPointLocal, |
---|
894 | localDirection, |
---|
895 | pCurrentProposedStepLength, |
---|
896 | pNewSafety, |
---|
897 | fHistory, |
---|
898 | fValidExitNormal, |
---|
899 | fExitNormal, |
---|
900 | exitingReplica, |
---|
901 | fEntering, |
---|
902 | &fBlockedPhysicalVolume, |
---|
903 | fBlockedReplicaNo); |
---|
904 | fExiting= exitingReplica; // still ok to set it ?? |
---|
905 | } |
---|
906 | |
---|
907 | // Remember last safety origin & value. |
---|
908 | // |
---|
909 | fPreviousSftOrigin = pGlobalpoint; |
---|
910 | fPreviousSafety = pNewSafety; |
---|
911 | |
---|
912 | // Count zero steps - one can occur due to changing momentum at a boundary |
---|
913 | // - one, two (or a few) can occur at common edges between |
---|
914 | // volumes |
---|
915 | // - more than two is likely a problem in the geometry |
---|
916 | // description or the Navigation |
---|
917 | |
---|
918 | // Rule of thumb: likely at an Edge if two consecutive steps are zero, |
---|
919 | // because at least two candidate volumes must have been |
---|
920 | // checked |
---|
921 | // |
---|
922 | fLocatedOnEdge = fLastStepWasZero && (Step==0.0); |
---|
923 | fLastStepWasZero = (Step==0.0); |
---|
924 | if (fPushed) fPushed = fLastStepWasZero; |
---|
925 | |
---|
926 | // Handle large number of consecutive zero steps |
---|
927 | // |
---|
928 | if ( fLastStepWasZero ) |
---|
929 | { |
---|
930 | fNumberZeroSteps++; |
---|
931 | #ifdef G4DEBUG_NAVIGATION |
---|
932 | if( fNumberZeroSteps > 1 ) |
---|
933 | { |
---|
934 | G4cout << "G4Navigator::ComputeStep(): another zero step, # " |
---|
935 | << fNumberZeroSteps |
---|
936 | << " at " << pGlobalpoint |
---|
937 | << " in volume " << motherPhysical->GetName() |
---|
938 | << " nav-comp-step calls # " << sNavCScalls |
---|
939 | << G4endl; |
---|
940 | } |
---|
941 | #endif |
---|
942 | if( fNumberZeroSteps > fActionThreshold_NoZeroSteps-1 ) |
---|
943 | { |
---|
944 | // Act to recover this stuck track. Pushing it along direction |
---|
945 | // |
---|
946 | Step += 0.9*kCarTolerance; |
---|
947 | #ifdef G4VERBOSE |
---|
948 | if (!fPushed) |
---|
949 | { |
---|
950 | G4cerr << "WARNING - G4Navigator::ComputeStep()" << G4endl |
---|
951 | << " Track stuck, not moving for " |
---|
952 | << fNumberZeroSteps << " steps" << G4endl |
---|
953 | << " in volume -" << motherPhysical->GetName() |
---|
954 | << "- at point " << pGlobalpoint << G4endl |
---|
955 | << " direction: " << pDirection << "." << G4endl |
---|
956 | << " Potential geometry or navigation problem !" |
---|
957 | << G4endl |
---|
958 | << " Trying pushing it of " << Step << " mm ..." |
---|
959 | << G4endl; |
---|
960 | } |
---|
961 | #endif |
---|
962 | fPushed = true; |
---|
963 | } |
---|
964 | if( fNumberZeroSteps > fAbandonThreshold_NoZeroSteps-1 ) |
---|
965 | { |
---|
966 | // Must kill this stuck track |
---|
967 | // |
---|
968 | G4cerr << "ERROR - G4Navigator::ComputeStep()" << G4endl |
---|
969 | << " Track stuck, not moving for " |
---|
970 | << fNumberZeroSteps << " steps" << G4endl |
---|
971 | << " in volume -" << motherPhysical->GetName() |
---|
972 | << "- at point " << pGlobalpoint << G4endl |
---|
973 | << " direction: " << pDirection << "." << G4endl; |
---|
974 | motherPhysical->CheckOverlaps(5000, false); |
---|
975 | G4Exception("G4Navigator::ComputeStep()", |
---|
976 | "StuckTrack", EventMustBeAborted, |
---|
977 | "Stuck Track: potential geometry or navigation problem."); |
---|
978 | } |
---|
979 | } |
---|
980 | else |
---|
981 | { |
---|
982 | if (!fPushed) fNumberZeroSteps = 0; |
---|
983 | } |
---|
984 | |
---|
985 | fEnteredDaughter = fEntering; // I expect to enter a volume in this Step |
---|
986 | fExitedMother = fExiting; |
---|
987 | |
---|
988 | if( fExiting ) |
---|
989 | { |
---|
990 | #ifdef G4DEBUG_NAVIGATION |
---|
991 | if( fVerbose > 2 ) |
---|
992 | { |
---|
993 | G4cout << " At G4Nav CompStep End - if(exiting) - fExiting= " << fExiting |
---|
994 | << " fValidExitNormal = " << fValidExitNormal << G4endl; |
---|
995 | G4cout << " fExitNormal= " << fExitNormal << G4endl; |
---|
996 | } |
---|
997 | #endif |
---|
998 | |
---|
999 | if(fValidExitNormal) |
---|
1000 | { |
---|
1001 | // Convention: fExitNormal is in the 'grand-mother' coordinate system |
---|
1002 | // |
---|
1003 | fGrandMotherExitNormal= fExitNormal; |
---|
1004 | } |
---|
1005 | else |
---|
1006 | { |
---|
1007 | // We must calculate the normal anyway (in order to have it if requested) |
---|
1008 | // |
---|
1009 | G4ThreeVector finalLocalPoint = |
---|
1010 | fLastLocatedPointLocal + localDirection*Step; |
---|
1011 | |
---|
1012 | // Now fGrandMotherExitNormal is in the 'grand-mother' coordinate system |
---|
1013 | // |
---|
1014 | fGrandMotherExitNormal = |
---|
1015 | motherLogical->GetSolid()->SurfaceNormal(finalLocalPoint); |
---|
1016 | |
---|
1017 | const G4RotationMatrix* mRot = motherPhysical->GetRotation(); |
---|
1018 | if( mRot ) |
---|
1019 | { |
---|
1020 | fGrandMotherExitNormal *= (*mRot).inverse(); |
---|
1021 | } |
---|
1022 | } |
---|
1023 | } |
---|
1024 | fStepEndPoint= pGlobalpoint+Step*pDirection; |
---|
1025 | |
---|
1026 | if( (Step == pCurrentProposedStepLength) && (!fExiting) && (!fEntering) ) |
---|
1027 | { |
---|
1028 | // This if Step is not really limited by the geometry. |
---|
1029 | // The Navigator is obliged to return "infinity" |
---|
1030 | // |
---|
1031 | Step = kInfinity; |
---|
1032 | } |
---|
1033 | |
---|
1034 | #ifdef G4VERBOSE |
---|
1035 | if( fVerbose > 1 ) |
---|
1036 | { |
---|
1037 | if( fVerbose >= 4 ) |
---|
1038 | { |
---|
1039 | G4cout << " ----- Upon exiting :" << G4endl; |
---|
1040 | PrintState(); |
---|
1041 | } |
---|
1042 | G4cout <<" Returned step = " << Step << G4endl; |
---|
1043 | if( Step == kInfinity ) |
---|
1044 | { |
---|
1045 | G4cout << " Original proposed step = " |
---|
1046 | << pCurrentProposedStepLength << G4endl; |
---|
1047 | } |
---|
1048 | G4cout << " Safety = " << pNewSafety << G4endl; |
---|
1049 | } |
---|
1050 | G4cout.precision(oldcoutPrec); |
---|
1051 | G4cerr.precision(oldcerrPrec); |
---|
1052 | #endif |
---|
1053 | |
---|
1054 | return Step; |
---|
1055 | } |
---|
1056 | |
---|
1057 | // ******************************************************************** |
---|
1058 | // CheckNextStep |
---|
1059 | // |
---|
1060 | // Compute the step without altering the navigator state |
---|
1061 | // ******************************************************************** |
---|
1062 | // |
---|
1063 | G4double G4Navigator::CheckNextStep( const G4ThreeVector& pGlobalpoint, |
---|
1064 | const G4ThreeVector& pDirection, |
---|
1065 | const G4double pCurrentProposedStepLength, |
---|
1066 | G4double& pNewSafety) |
---|
1067 | { |
---|
1068 | G4double step; |
---|
1069 | |
---|
1070 | // Save the state, for this parasitic call |
---|
1071 | // |
---|
1072 | SetSavedState(); |
---|
1073 | |
---|
1074 | step = ComputeStep ( pGlobalpoint, |
---|
1075 | pDirection, |
---|
1076 | pCurrentProposedStepLength, |
---|
1077 | pNewSafety ); |
---|
1078 | |
---|
1079 | // If a parasitic call, then attempt to restore the key parts of the state |
---|
1080 | // |
---|
1081 | RestoreSavedState(); |
---|
1082 | |
---|
1083 | return step; |
---|
1084 | } |
---|
1085 | |
---|
1086 | // ******************************************************************** |
---|
1087 | // ResetState |
---|
1088 | // |
---|
1089 | // Resets stack and minimum of navigator state `machine' |
---|
1090 | // ******************************************************************** |
---|
1091 | // |
---|
1092 | void G4Navigator::ResetState() |
---|
1093 | { |
---|
1094 | fWasLimitedByGeometry = false; |
---|
1095 | fEntering = false; |
---|
1096 | fExiting = false; |
---|
1097 | fLocatedOnEdge = false; |
---|
1098 | fLastStepWasZero = false; |
---|
1099 | fEnteredDaughter = false; |
---|
1100 | fExitedMother = false; |
---|
1101 | fPushed = false; |
---|
1102 | |
---|
1103 | fValidExitNormal = false; |
---|
1104 | fExitNormal = G4ThreeVector(0,0,0); |
---|
1105 | |
---|
1106 | fPreviousSftOrigin = G4ThreeVector(0,0,0); |
---|
1107 | fPreviousSafety = 0.0; |
---|
1108 | |
---|
1109 | fNumberZeroSteps = 0; |
---|
1110 | |
---|
1111 | fBlockedPhysicalVolume = 0; |
---|
1112 | fBlockedReplicaNo = -1; |
---|
1113 | |
---|
1114 | fLastLocatedPointLocal = G4ThreeVector( DBL_MAX, -DBL_MAX, 0.0 ); |
---|
1115 | fLocatedOutsideWorld = false; |
---|
1116 | } |
---|
1117 | |
---|
1118 | // ******************************************************************** |
---|
1119 | // SetupHierarchy |
---|
1120 | // |
---|
1121 | // Renavigates & resets hierarchy described by current history |
---|
1122 | // o Reset volumes |
---|
1123 | // o Recompute transforms and/or solids of replicated/parameterised volumes |
---|
1124 | // ******************************************************************** |
---|
1125 | // |
---|
1126 | void G4Navigator::SetupHierarchy() |
---|
1127 | { |
---|
1128 | G4int i; |
---|
1129 | const G4int cdepth = fHistory.GetDepth(); |
---|
1130 | G4VPhysicalVolume *mother, *current; |
---|
1131 | G4VSolid *pSolid; |
---|
1132 | G4VPVParameterisation *pParam; |
---|
1133 | |
---|
1134 | mother = fHistory.GetVolume(0); |
---|
1135 | for ( i=1; i<=cdepth; i++ ) |
---|
1136 | { |
---|
1137 | current = fHistory.GetVolume(i); |
---|
1138 | switch ( fHistory.GetVolumeType(i) ) |
---|
1139 | { |
---|
1140 | case kNormal: |
---|
1141 | break; |
---|
1142 | case kReplica: |
---|
1143 | freplicaNav.ComputeTransformation(fHistory.GetReplicaNo(i), current); |
---|
1144 | break; |
---|
1145 | case kParameterised: |
---|
1146 | G4int replicaNo; |
---|
1147 | pParam = current->GetParameterisation(); |
---|
1148 | replicaNo = fHistory.GetReplicaNo(i); |
---|
1149 | pSolid = pParam->ComputeSolid(replicaNo, current); |
---|
1150 | |
---|
1151 | // Set up dimensions & transform in solid/physical volume |
---|
1152 | // |
---|
1153 | pSolid->ComputeDimensions(pParam, replicaNo, current); |
---|
1154 | pParam->ComputeTransformation(replicaNo, current); |
---|
1155 | |
---|
1156 | G4TouchableHistory touchable( fHistory ); |
---|
1157 | touchable.MoveUpHistory(); // move up to the parent level |
---|
1158 | |
---|
1159 | // Set up the correct solid and material in Logical Volume |
---|
1160 | // |
---|
1161 | G4LogicalVolume *pLogical = current->GetLogicalVolume(); |
---|
1162 | pLogical->SetSolid( pSolid ); |
---|
1163 | pLogical->UpdateMaterial( pParam -> |
---|
1164 | ComputeMaterial(replicaNo, current, &touchable) ); |
---|
1165 | break; |
---|
1166 | } |
---|
1167 | mother = current; |
---|
1168 | } |
---|
1169 | } |
---|
1170 | |
---|
1171 | // ******************************************************************** |
---|
1172 | // GetLocalExitNormal |
---|
1173 | // |
---|
1174 | // Obtains the Normal vector to a surface (in local coordinates) |
---|
1175 | // pointing out of previous volume and into current volume |
---|
1176 | // ******************************************************************** |
---|
1177 | // |
---|
1178 | G4ThreeVector G4Navigator::GetLocalExitNormal( G4bool* valid ) |
---|
1179 | { |
---|
1180 | G4ThreeVector ExitNormal(0.,0.,0.); |
---|
1181 | |
---|
1182 | if ( EnteredDaughterVolume() ) |
---|
1183 | { |
---|
1184 | ExitNormal= -(fHistory.GetTopVolume()->GetLogicalVolume()-> |
---|
1185 | GetSolid()->SurfaceNormal(fLastLocatedPointLocal)); |
---|
1186 | *valid = true; |
---|
1187 | } |
---|
1188 | else |
---|
1189 | { |
---|
1190 | if( fExitedMother ) |
---|
1191 | { |
---|
1192 | ExitNormal = fGrandMotherExitNormal; |
---|
1193 | *valid = true; |
---|
1194 | } |
---|
1195 | else |
---|
1196 | { |
---|
1197 | // We are not at a boundary. |
---|
1198 | // ExitNormal remains (0,0,0) |
---|
1199 | // |
---|
1200 | *valid = false; |
---|
1201 | } |
---|
1202 | } |
---|
1203 | return ExitNormal; |
---|
1204 | } |
---|
1205 | |
---|
1206 | // ******************************************************************** |
---|
1207 | // ComputeSafety |
---|
1208 | // |
---|
1209 | // It assumes that it will be |
---|
1210 | // i) called at the Point in the same volume as the EndPoint of the |
---|
1211 | // ComputeStep. |
---|
1212 | // ii) after (or at the end of) ComputeStep OR after the relocation. |
---|
1213 | // ******************************************************************** |
---|
1214 | // |
---|
1215 | G4double G4Navigator::ComputeSafety( const G4ThreeVector &pGlobalpoint, |
---|
1216 | const G4double pMaxLength, |
---|
1217 | const G4bool keepState) |
---|
1218 | { |
---|
1219 | G4double newSafety = 0.0; |
---|
1220 | |
---|
1221 | #ifdef G4DEBUG_NAVIGATION |
---|
1222 | G4int oldcoutPrec = G4cout.precision(8); |
---|
1223 | if( fVerbose > 0 ) |
---|
1224 | { |
---|
1225 | G4cout << "*** G4Navigator::ComputeSafety: ***" << G4endl |
---|
1226 | << " Called at point: " << pGlobalpoint << G4endl; |
---|
1227 | |
---|
1228 | G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume(); |
---|
1229 | G4cout << " Volume = " << motherPhysical->GetName() |
---|
1230 | << " - Maximum length = " << pMaxLength << G4endl; |
---|
1231 | if( fVerbose >= 4 ) |
---|
1232 | { |
---|
1233 | G4cout << " ----- Upon entering Compute Safety:" << G4endl; |
---|
1234 | PrintState(); |
---|
1235 | } |
---|
1236 | } |
---|
1237 | #endif |
---|
1238 | |
---|
1239 | if (keepState) { SetSavedState(); } |
---|
1240 | |
---|
1241 | G4double distEndpointSq = (pGlobalpoint-fStepEndPoint).mag2(); |
---|
1242 | G4bool stayedOnEndpoint = distEndpointSq < kCarTolerance*kCarTolerance; |
---|
1243 | G4bool endpointOnSurface = fEnteredDaughter || fExitedMother; |
---|
1244 | |
---|
1245 | if( !(endpointOnSurface && stayedOnEndpoint) ) |
---|
1246 | { |
---|
1247 | // Pseudo-relocate to this point (updates voxel information only) |
---|
1248 | // |
---|
1249 | LocateGlobalPointWithinVolume( pGlobalpoint ); |
---|
1250 | // --->> Danger: Side effects on sub-navigator voxel information <<--- |
---|
1251 | // Could be replaced again by 'granular' calls to sub-navigator |
---|
1252 | // locates (similar side-effects, but faster. |
---|
1253 | // Solutions: |
---|
1254 | // 1) Re-locate (to where?) |
---|
1255 | // 2) Insure that the methods using (G4ComputeStep?) |
---|
1256 | // does a relocation (if information is disturbed only ?) |
---|
1257 | |
---|
1258 | #ifdef G4DEBUG_NAVIGATION |
---|
1259 | if( fVerbose >= 2 ) |
---|
1260 | { |
---|
1261 | G4cout << " G4Navigator::ComputeSafety() relocates-in-volume to point: " |
---|
1262 | << pGlobalpoint << G4endl; |
---|
1263 | } |
---|
1264 | #endif |
---|
1265 | G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume(); |
---|
1266 | G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume(); |
---|
1267 | G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader(); |
---|
1268 | G4ThreeVector localPoint = ComputeLocalPoint(pGlobalpoint); |
---|
1269 | |
---|
1270 | if ( fHistory.GetTopVolumeType()!=kReplica ) |
---|
1271 | { |
---|
1272 | switch(CharacteriseDaughters(motherLogical)) |
---|
1273 | { |
---|
1274 | case kNormal: |
---|
1275 | if ( pVoxelHeader ) |
---|
1276 | { |
---|
1277 | newSafety=fvoxelNav.ComputeSafety(localPoint,fHistory,pMaxLength); |
---|
1278 | } |
---|
1279 | else |
---|
1280 | { |
---|
1281 | newSafety=fnormalNav.ComputeSafety(localPoint,fHistory,pMaxLength); |
---|
1282 | } |
---|
1283 | break; |
---|
1284 | case kParameterised: |
---|
1285 | if( GetDaughtersRegularStructureId(motherLogical) != 1 ) |
---|
1286 | { |
---|
1287 | newSafety = fparamNav.ComputeSafety(localPoint,fHistory,pMaxLength); |
---|
1288 | } |
---|
1289 | else // Regular structure |
---|
1290 | { |
---|
1291 | newSafety = fregularNav.ComputeSafety(localPoint,fHistory,pMaxLength); |
---|
1292 | } |
---|
1293 | break; |
---|
1294 | case kReplica: |
---|
1295 | G4Exception("G4Navigator::ComputeSafety()", "NotApplicable", |
---|
1296 | FatalException, "Not applicable for replicated volumes."); |
---|
1297 | break; |
---|
1298 | } |
---|
1299 | } |
---|
1300 | else |
---|
1301 | { |
---|
1302 | newSafety = freplicaNav.ComputeSafety(pGlobalpoint, localPoint, |
---|
1303 | fHistory, pMaxLength); |
---|
1304 | } |
---|
1305 | } |
---|
1306 | else // if( endpointOnSurface && stayedOnEndpoint ) |
---|
1307 | { |
---|
1308 | #ifdef G4DEBUG_NAVIGATION |
---|
1309 | if( fVerbose >= 2 ) |
---|
1310 | { |
---|
1311 | G4cout << " G4Navigator::ComputeSafety() finds that point - " |
---|
1312 | << pGlobalpoint << " - is on surface " << G4endl; |
---|
1313 | if( fEnteredDaughter ) { G4cout << " entered new daughter volume"; } |
---|
1314 | if( fExitedMother ) { G4cout << " and exited previous volume."; } |
---|
1315 | G4cout << G4endl; |
---|
1316 | G4cout << " EndPoint was = " << fStepEndPoint << G4endl; |
---|
1317 | } |
---|
1318 | #endif |
---|
1319 | newSafety = 0.0; |
---|
1320 | } |
---|
1321 | |
---|
1322 | // Remember last safety origin & value |
---|
1323 | // |
---|
1324 | fPreviousSftOrigin = pGlobalpoint; |
---|
1325 | fPreviousSafety = newSafety; |
---|
1326 | |
---|
1327 | if (keepState) { RestoreSavedState(); } |
---|
1328 | |
---|
1329 | #ifdef G4DEBUG_NAVIGATION |
---|
1330 | if( fVerbose > 1 ) |
---|
1331 | { |
---|
1332 | G4cout << " ---- Exiting ComputeSafety " << G4endl; |
---|
1333 | if( fVerbose > 2 ) { PrintState(); } |
---|
1334 | G4cout << " Returned value of Safety = " << newSafety << G4endl; |
---|
1335 | } |
---|
1336 | G4cout.precision(oldcoutPrec); |
---|
1337 | #endif |
---|
1338 | |
---|
1339 | return newSafety; |
---|
1340 | } |
---|
1341 | |
---|
1342 | // ******************************************************************** |
---|
1343 | // CreateTouchableHistoryHandle |
---|
1344 | // ******************************************************************** |
---|
1345 | // |
---|
1346 | G4TouchableHistoryHandle G4Navigator::CreateTouchableHistoryHandle() const |
---|
1347 | { |
---|
1348 | return G4TouchableHistoryHandle( CreateTouchableHistory() ); |
---|
1349 | } |
---|
1350 | |
---|
1351 | // ******************************************************************** |
---|
1352 | // PrintState |
---|
1353 | // ******************************************************************** |
---|
1354 | // |
---|
1355 | void G4Navigator::PrintState() const |
---|
1356 | { |
---|
1357 | G4int oldcoutPrec = G4cout.precision(4); |
---|
1358 | if( fVerbose == 4 ) |
---|
1359 | { |
---|
1360 | G4cout << "The current state of G4Navigator is: " << G4endl; |
---|
1361 | G4cout << " ValidExitNormal= " << fValidExitNormal << G4endl |
---|
1362 | << " ExitNormal = " << fExitNormal << G4endl |
---|
1363 | << " Exiting = " << fExiting << G4endl |
---|
1364 | << " Entering = " << fEntering << G4endl |
---|
1365 | << " BlockedPhysicalVolume= " ; |
---|
1366 | if (fBlockedPhysicalVolume==0) |
---|
1367 | G4cout << "None"; |
---|
1368 | else |
---|
1369 | G4cout << fBlockedPhysicalVolume->GetName(); |
---|
1370 | G4cout << G4endl |
---|
1371 | << " BlockedReplicaNo = " << fBlockedReplicaNo << G4endl |
---|
1372 | << " LastStepWasZero = " << fLastStepWasZero << G4endl |
---|
1373 | << G4endl; |
---|
1374 | } |
---|
1375 | if( ( 1 < fVerbose) && (fVerbose < 4) ) |
---|
1376 | { |
---|
1377 | G4cout << std::setw(30) << " ExitNormal " << " " |
---|
1378 | << std::setw( 5) << " Valid " << " " |
---|
1379 | << std::setw( 9) << " Exiting " << " " |
---|
1380 | << std::setw( 9) << " Entering" << " " |
---|
1381 | << std::setw(15) << " Blocked:Volume " << " " |
---|
1382 | << std::setw( 9) << " ReplicaNo" << " " |
---|
1383 | << std::setw( 8) << " LastStepZero " << " " |
---|
1384 | << G4endl; |
---|
1385 | G4cout << "( " << std::setw(7) << fExitNormal.x() |
---|
1386 | << ", " << std::setw(7) << fExitNormal.y() |
---|
1387 | << ", " << std::setw(7) << fExitNormal.z() << " ) " |
---|
1388 | << std::setw( 5) << fValidExitNormal << " " |
---|
1389 | << std::setw( 9) << fExiting << " " |
---|
1390 | << std::setw( 9) << fEntering << " "; |
---|
1391 | if ( fBlockedPhysicalVolume==0 ) |
---|
1392 | G4cout << std::setw(15) << "None"; |
---|
1393 | else |
---|
1394 | G4cout << std::setw(15)<< fBlockedPhysicalVolume->GetName(); |
---|
1395 | G4cout << std::setw( 9) << fBlockedReplicaNo << " " |
---|
1396 | << std::setw( 8) << fLastStepWasZero << " " |
---|
1397 | << G4endl; |
---|
1398 | } |
---|
1399 | if( fVerbose > 2 ) |
---|
1400 | { |
---|
1401 | G4cout.precision(8); |
---|
1402 | G4cout << " Current Localpoint = " << fLastLocatedPointLocal << G4endl; |
---|
1403 | G4cout << " PreviousSftOrigin = " << fPreviousSftOrigin << G4endl; |
---|
1404 | G4cout << " PreviousSafety = " << fPreviousSafety << G4endl; |
---|
1405 | } |
---|
1406 | G4cout.precision(oldcoutPrec); |
---|
1407 | } |
---|
1408 | |
---|
1409 | // ******************************************************************** |
---|
1410 | // Operator << |
---|
1411 | // ******************************************************************** |
---|
1412 | // |
---|
1413 | std::ostream& operator << (std::ostream &os,const G4Navigator &n) |
---|
1414 | { |
---|
1415 | os << "Current History: " << G4endl << n.fHistory; |
---|
1416 | return os; |
---|
1417 | } |
---|