source: trunk/source/geometry/navigation/src/G4VoxelNavigation.cc@ 892

Last change on this file since 892 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 26.4 KB
Line 
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: G4VoxelNavigation.cc,v 1.7 2007/05/11 13:43:59 gcosmo Exp $
28// GEANT4 tag $Name: HEAD $
29//
30//
31// class G4VoxelNavigation Implementation
32//
33// Author: P.Kent, 1996
34//
35// --------------------------------------------------------------------
36
37#include "G4VoxelNavigation.hh"
38#include "G4GeometryTolerance.hh"
39
40// ********************************************************************
41// Constructor
42// ********************************************************************
43//
44G4VoxelNavigation::G4VoxelNavigation()
45 : fVoxelDepth(-1),
46 fVoxelAxisStack(kNavigatorVoxelStackMax,kXAxis),
47 fVoxelNoSlicesStack(kNavigatorVoxelStackMax,0),
48 fVoxelSliceWidthStack(kNavigatorVoxelStackMax,0.),
49 fVoxelNodeNoStack(kNavigatorVoxelStackMax,0),
50 fVoxelHeaderStack(kNavigatorVoxelStackMax,(G4SmartVoxelHeader*)0),
51 fVoxelNode(0),
52 fCheck(false),
53 fVerbose(0)
54{
55 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
56}
57
58// ********************************************************************
59// Destructor
60// ********************************************************************
61//
62G4VoxelNavigation::~G4VoxelNavigation()
63{
64#ifdef G4DEBUG_NAVIGATION
65 G4cout << "G4VoxelNavigation::~G4VoxelNavigation() called." << G4endl;
66#endif
67}
68
69// ********************************************************************
70// ComputeStep
71// ********************************************************************
72//
73G4double
74G4VoxelNavigation::ComputeStep( const G4ThreeVector& localPoint,
75 const G4ThreeVector& localDirection,
76 const G4double currentProposedStepLength,
77 G4double& newSafety,
78 G4NavigationHistory& history,
79 G4bool& validExitNormal,
80 G4ThreeVector& exitNormal,
81 G4bool& exiting,
82 G4bool& entering,
83 G4VPhysicalVolume *(*pBlockedPhysical),
84 G4int& blockedReplicaNo )
85{
86 G4VPhysicalVolume *motherPhysical, *samplePhysical, *blockedExitedVol=0;
87 G4LogicalVolume *motherLogical;
88 G4VSolid *motherSolid;
89 G4ThreeVector sampleDirection;
90 G4double ourStep=currentProposedStepLength, motherSafety, ourSafety;
91 G4int localNoDaughters, sampleNo;
92
93 G4bool initialNode, noStep;
94 G4SmartVoxelNode *curVoxelNode;
95 G4int curNoVolumes, contentNo;
96 G4double voxelSafety;
97
98 motherPhysical = history.GetTopVolume();
99 motherLogical = motherPhysical->GetLogicalVolume();
100 motherSolid = motherLogical->GetSolid();
101
102 //
103 // Compute mother safety
104 //
105
106 motherSafety = motherSolid->DistanceToOut(localPoint);
107 ourSafety = motherSafety; // Working isotropic safety
108
109#ifdef G4VERBOSE
110 if ( fCheck )
111 {
112 if(fVerbose == 1 )
113 {
114 G4cout << "*** G4VoxelNavigation::ComputeStep(): ***" << G4endl
115 << " Invoked DistanceToOut(p) for mother solid: "
116 << motherSolid->GetName()
117 << ". Solid replied: " << motherSafety << G4endl
118 << " For local point p: " << localPoint
119 << ", to be considered as 'mother safety'." << G4endl;
120 }
121 if( motherSafety < 0.0 )
122 {
123 G4cout << "ERROR - G4VoxelNavigation::ComputeStep()" << G4endl
124 << " Current solid " << motherSolid->GetName()
125 << " gave negative safety: " << motherSafety << G4endl
126 << " for the current (local) point " << localPoint
127 << G4endl;
128 motherSolid->DumpInfo();
129 G4Exception("G4VoxelNavigation::ComputeStep()",
130 "NegativeSafetyMotherVol", FatalException,
131 "Negative Safety In Voxel Navigation !" );
132 }
133 if( motherSolid->Inside(localPoint)==kOutside )
134 {
135 G4cout << "WARNING - G4VoxelNavigation::ComputeStep()" << G4endl
136 << " Point " << localPoint
137 << " is outside current volume " << motherPhysical->GetName()
138 << G4endl;
139 G4double estDistToSolid= motherSolid->DistanceToIn(localPoint);
140 G4cout << " Estimated isotropic distance to solid (distToIn)= "
141 << estDistToSolid << G4endl;
142 if( estDistToSolid > 100.0 * kCarTolerance )
143 {
144 motherSolid->DumpInfo();
145 G4Exception("G4VoxelNavigation::ComputeStep()",
146 "FarOutsideCurrentVolume", FatalException,
147 "Point is far outside Current Volume !");
148 }
149 else
150 G4Exception("G4VoxelNavigation::ComputeStep()", "OutsideCurrentVolume",
151 JustWarning, "Point is a little outside Current Volume.");
152 }
153 }
154#endif
155
156 //
157 // Compute daughter safeties & intersections
158 //
159
160 // Exiting normal optimisation
161 //
162 if ( exiting && validExitNormal )
163 {
164 if ( localDirection.dot(exitNormal)>=kMinExitingNormalCosine )
165 {
166 // Block exited daughter volume
167 //
168 blockedExitedVol = *pBlockedPhysical;
169 ourSafety = 0;
170 }
171 }
172 exiting = false;
173 entering = false;
174
175 localNoDaughters = motherLogical->GetNoDaughters();
176
177 fBList.Enlarge(localNoDaughters);
178 fBList.Reset();
179
180 initialNode = true;
181 noStep = true;
182
183 while (noStep)
184 {
185 curVoxelNode = fVoxelNode;
186 curNoVolumes = curVoxelNode->GetNoContained();
187 for (contentNo=curNoVolumes-1; contentNo>=0; contentNo--)
188 {
189 sampleNo = curVoxelNode->GetVolume(contentNo);
190 if ( !fBList.IsBlocked(sampleNo) )
191 {
192 fBList.BlockVolume(sampleNo);
193 samplePhysical = motherLogical->GetDaughter(sampleNo);
194 if ( samplePhysical!=blockedExitedVol )
195 {
196 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
197 samplePhysical->GetTranslation());
198 sampleTf.Invert();
199 const G4ThreeVector samplePoint =
200 sampleTf.TransformPoint(localPoint);
201 const G4VSolid *sampleSolid =
202 samplePhysical->GetLogicalVolume()->GetSolid();
203 const G4double sampleSafety =
204 sampleSolid->DistanceToIn(samplePoint);
205#ifdef G4VERBOSE
206 if(( fCheck ) && ( fVerbose == 1 ))
207 {
208 G4cout << "*** G4VoxelNavigation::ComputeStep(): ***" << G4endl
209 << " Invoked DistanceToIn(p) for daughter solid: "
210 << sampleSolid->GetName()
211 << ". Solid replied: " << sampleSafety << G4endl
212 << " For local point p: " << samplePoint
213 << ", to be considered as 'daughter safety'." << G4endl;
214 }
215#endif
216 if ( sampleSafety<ourSafety )
217 {
218 ourSafety = sampleSafety;
219 }
220 if ( sampleSafety<=ourStep )
221 {
222 sampleDirection = sampleTf.TransformAxis(localDirection);
223 G4double sampleStep =
224 sampleSolid->DistanceToIn(samplePoint, sampleDirection);
225#ifdef G4VERBOSE
226 if(( fCheck ) && ( fVerbose == 1 ))
227 {
228 G4cout << "*** G4VoxelNavigation::ComputeStep(): ***" << G4endl
229 << " Invoked DistanceToIn(p,v) for daughter solid: "
230 << sampleSolid->GetName()
231 << ". Solid replied: " << sampleStep << G4endl
232 << " For local point p: " << samplePoint << G4endl
233 << " Direction v: " << sampleDirection
234 << ", to be considered as 'daughter step'." << G4endl;
235 }
236#endif
237 if ( sampleStep<=ourStep )
238 {
239 ourStep = sampleStep;
240 entering = true;
241 exiting = false;
242 *pBlockedPhysical = samplePhysical;
243 blockedReplicaNo = -1;
244#ifdef G4VERBOSE
245 // Check to see that the resulting point is indeed in/on volume.
246 // This check could eventually be made only for successful
247 // candidate.
248
249 if ( ( fCheck ) && ( sampleStep < kInfinity ) )
250 {
251 G4ThreeVector intersectionPoint;
252 intersectionPoint= samplePoint + sampleStep * sampleDirection;
253 EInside insideIntPt= sampleSolid->Inside(intersectionPoint);
254 G4String solidResponse = "-kInside-";
255 if (insideIntPt == kOutside)
256 solidResponse = "-kOutside-";
257 else if (insideIntPt == kSurface)
258 solidResponse = "-kSurface-";
259 if( fVerbose == 1 )
260 {
261 G4cout << "*** G4VoxelNavigation::ComputeStep(): ***"<<G4endl
262 << " Invoked Inside() for solid: "
263 << sampleSolid->GetName()
264 << ". Solid replied: " << solidResponse << G4endl
265 << " For point p: " << intersectionPoint
266 << ", considered as 'intersection' point." << G4endl;
267 }
268 if( insideIntPt != kSurface )
269 {
270 G4int oldcoutPrec = G4cout.precision(16);
271 G4cout << "WARNING - G4VoxelNavigation::ComputeStep()"
272 << G4endl
273 << " Inaccurate solid DistanceToIn"
274 << " for solid " << sampleSolid->GetName() << G4endl;
275 G4cout << " Solid gave DistanceToIn = "
276 << sampleStep << " yet returns " << solidResponse
277 << " for this point !" << G4endl;
278 G4cout << " Point = " << intersectionPoint << G4endl;
279 if ( insideIntPt != kInside )
280 G4cout << " DistanceToIn(p) = "
281 << sampleSolid->DistanceToIn(intersectionPoint)
282 << G4endl;
283 if ( insideIntPt != kOutside )
284 G4cout << " DistanceToOut(p) = "
285 << sampleSolid->DistanceToOut(intersectionPoint)
286 << G4endl;
287 G4Exception("G4VoxelNavigation::ComputeStep()",
288 "InaccurateDistanceToIn", JustWarning,
289 "Navigator gets conflicting response from Solid.");
290 G4cout.precision(oldcoutPrec);
291 }
292 }
293#endif
294 }
295 }
296 }
297 }
298 }
299 if (initialNode)
300 {
301 initialNode = false;
302 voxelSafety = ComputeVoxelSafety(localPoint);
303 if ( voxelSafety<ourSafety )
304 {
305 ourSafety = voxelSafety;
306 }
307 if ( currentProposedStepLength<ourSafety )
308 {
309 // Guaranteed physics limited
310 //
311 noStep = false;
312 entering = false;
313 exiting = false;
314 *pBlockedPhysical = 0;
315 ourStep = kInfinity;
316 }
317 else
318 {
319 //
320 // Compute mother intersection if required
321 //
322 if ( motherSafety<=ourStep )
323 {
324 G4double motherStep =
325 motherSolid->DistanceToOut(localPoint,
326 localDirection,
327 true, &validExitNormal, &exitNormal);
328#ifdef G4VERBOSE
329 if ( fCheck )
330 {
331 if(fVerbose == 1)
332 {
333 G4cout << "*** G4VoxelNavigation::ComputeStep(): ***" << G4endl
334 << " Invoked DistanceToOut(p,v,...) for mother solid: "
335 << motherSolid->GetName()
336 << ". Solid replied: " << motherStep << G4endl
337 << " For local point p: " << localPoint << G4endl
338 << " Direction v: " << localDirection
339 << ", to be considered as 'mother step'." << G4endl;
340 }
341 if( ( motherStep < 0.0 ) || ( motherStep >= kInfinity) )
342 {
343 G4int oldPrOut= G4cout.precision(16);
344 G4int oldPrErr= G4cerr.precision(16);
345 G4cerr << "ERROR - G4VoxelNavigation::ComputeStep()" << G4endl
346 << " Problem in Navigation" << G4endl
347 << " Point (local coordinates): "
348 << localPoint << G4endl
349 << " Local Direction: " << localDirection << G4endl
350 << " Solid: " << motherSolid->GetName() << G4endl;
351 motherSolid->DumpInfo();
352 G4Exception("G4VoxelNavigation::ComputeStep()",
353 "PointOutsideCurrentVolume", FatalException,
354 "Current point is outside the current solid !");
355 G4cout.precision(oldPrOut);
356 G4cerr.precision(oldPrErr);
357 }
358 }
359#endif
360 if ( motherStep<=ourStep )
361 {
362 ourStep = motherStep;
363 exiting = true;
364 entering = false;
365 if ( validExitNormal )
366 {
367 const G4RotationMatrix *rot = motherPhysical->GetRotation();
368 if (rot)
369 {
370 exitNormal *= rot->inverse();
371 }
372 }
373 }
374 else
375 {
376 validExitNormal = false;
377 }
378 }
379 }
380 newSafety = ourSafety;
381 }
382 if (noStep)
383 {
384 noStep = LocateNextVoxel(localPoint, localDirection, ourStep);
385 }
386 } // end -while (noStep)- loop
387
388 return ourStep;
389}
390
391// ********************************************************************
392// ComputeVoxelSafety
393//
394// Computes safety from specified point to voxel boundaries
395// using already located point
396// o collected boundaries for most derived level
397// o adjacent boundaries for previous levels
398// ********************************************************************
399//
400G4double
401G4VoxelNavigation::ComputeVoxelSafety(const G4ThreeVector& localPoint) const
402{
403 G4SmartVoxelHeader *curHeader;
404 G4double voxelSafety, curNodeWidth;
405 G4double curNodeOffset, minCurCommonDelta, maxCurCommonDelta;
406 G4int minCurNodeNoDelta, maxCurNodeNoDelta;
407 G4int localVoxelDepth, curNodeNo;
408 EAxis curHeaderAxis;
409
410 localVoxelDepth = fVoxelDepth;
411
412 curHeader = fVoxelHeaderStack[localVoxelDepth];
413 curHeaderAxis = fVoxelAxisStack[localVoxelDepth];
414 curNodeNo = fVoxelNodeNoStack[localVoxelDepth];
415 curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth];
416
417 // Compute linear intersection distance to boundaries of max/min
418 // to collected nodes at current level
419 //
420 curNodeOffset = curNodeNo*curNodeWidth;
421 maxCurNodeNoDelta = fVoxelNode->GetMaxEquivalentSliceNo()-curNodeNo;
422 minCurNodeNoDelta = curNodeNo-fVoxelNode->GetMinEquivalentSliceNo();
423 minCurCommonDelta = localPoint(curHeaderAxis)
424 - curHeader->GetMinExtent() - curNodeOffset;
425 maxCurCommonDelta = curNodeWidth-minCurCommonDelta;
426
427 if ( minCurNodeNoDelta<maxCurNodeNoDelta )
428 {
429 voxelSafety = minCurNodeNoDelta*curNodeWidth;
430 voxelSafety += minCurCommonDelta;
431 }
432 else if (maxCurNodeNoDelta < minCurNodeNoDelta)
433 {
434 voxelSafety = maxCurNodeNoDelta*curNodeWidth;
435 voxelSafety += maxCurCommonDelta;
436 }
437 else // (maxCurNodeNoDelta == minCurNodeNoDelta)
438 {
439 voxelSafety = minCurNodeNoDelta*curNodeWidth;
440 voxelSafety += std::min(minCurCommonDelta,maxCurCommonDelta);
441 }
442
443 // Compute isotropic safety to boundaries of previous levels
444 // [NOT to collected boundaries]
445 //
446 while ( (localVoxelDepth>0) && (voxelSafety>0) )
447 {
448 localVoxelDepth--;
449 curHeader = fVoxelHeaderStack[localVoxelDepth];
450 curHeaderAxis = fVoxelAxisStack[localVoxelDepth];
451 curNodeNo = fVoxelNodeNoStack[localVoxelDepth];
452 curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth];
453 curNodeOffset = curNodeNo*curNodeWidth;
454 minCurCommonDelta = localPoint(curHeaderAxis)
455 - curHeader->GetMinExtent() - curNodeOffset;
456 maxCurCommonDelta = curNodeWidth-minCurCommonDelta;
457
458 if ( minCurCommonDelta<voxelSafety )
459 {
460 voxelSafety = minCurCommonDelta;
461 }
462 if ( maxCurCommonDelta<voxelSafety )
463 {
464 voxelSafety = maxCurCommonDelta;
465 }
466 }
467 if ( voxelSafety<0 )
468 {
469 voxelSafety = 0;
470 }
471
472 return voxelSafety;
473}
474
475// ********************************************************************
476// LocateNextVoxel
477//
478// Finds the next voxel from the current voxel and point
479// in the specified direction
480//
481// Returns false if all voxels considered
482// [current Step ends inside same voxel or leaves all voxels]
483// true otherwise
484// [the information on the next voxel is put into the set of
485// fVoxel* variables & "stacks"]
486// ********************************************************************
487//
488G4bool
489G4VoxelNavigation::LocateNextVoxel(const G4ThreeVector& localPoint,
490 const G4ThreeVector& localDirection,
491 const G4double currentStep)
492{
493 G4SmartVoxelHeader *workHeader=0, *newHeader=0;
494 G4SmartVoxelProxy *newProxy=0;
495 G4SmartVoxelNode *newVoxelNode=0;
496 G4ThreeVector targetPoint, voxelPoint;
497 G4double workNodeWidth, workMinExtent, workCoord;
498 G4double minVal, maxVal, newDistance=0.;
499 G4double newHeaderMin, newHeaderNodeWidth;
500 G4int depth=0, newDepth=0, workNodeNo=0, newNodeNo=0, newHeaderNoSlices=0;
501 EAxis workHeaderAxis, newHeaderAxis;
502 G4bool isNewVoxel=false;
503
504 G4double currentDistance = currentStep;
505
506 // Determine if end of Step within current voxel
507 //
508 for (depth=0; depth<fVoxelDepth; depth++)
509 {
510 targetPoint = localPoint+localDirection*currentDistance;
511 newDistance = currentDistance;
512 workHeader = fVoxelHeaderStack[depth];
513 workHeaderAxis = fVoxelAxisStack[depth];
514 workNodeNo = fVoxelNodeNoStack[depth];
515 workNodeWidth = fVoxelSliceWidthStack[depth];
516 workMinExtent = workHeader->GetMinExtent();
517 workCoord = targetPoint(workHeaderAxis);
518 minVal = workMinExtent+workNodeNo*workNodeWidth;
519
520 if ( minVal<=workCoord+kCarTolerance*0.5 )
521 {
522 maxVal = minVal+workNodeWidth;
523 if ( maxVal<=workCoord-kCarTolerance*0.5 )
524 {
525 // Must consider next voxel
526 //
527 newNodeNo = workNodeNo+1;
528 newHeader = workHeader;
529 newDistance = (maxVal-localPoint(workHeaderAxis))
530 / localDirection(workHeaderAxis);
531 isNewVoxel = true;
532 newDepth = depth;
533 }
534 }
535 else
536 {
537 newNodeNo = workNodeNo-1;
538 newHeader = workHeader;
539 newDistance = (minVal-localPoint(workHeaderAxis))
540 / localDirection(workHeaderAxis);
541 isNewVoxel = true;
542 newDepth = depth;
543 }
544 currentDistance = newDistance;
545 }
546 targetPoint = localPoint+localDirection*currentDistance;
547
548 // Check if end of Step within collected boundaries of current voxel
549 //
550 depth = fVoxelDepth;
551 {
552 workHeader = fVoxelHeaderStack[depth];
553 workHeaderAxis = fVoxelAxisStack[depth];
554 workNodeNo = fVoxelNodeNoStack[depth];
555 workNodeWidth = fVoxelSliceWidthStack[depth];
556 workMinExtent = workHeader->GetMinExtent();
557 workCoord = targetPoint(workHeaderAxis);
558 minVal = workMinExtent+fVoxelNode->GetMinEquivalentSliceNo()*workNodeWidth;
559
560 if ( minVal<=workCoord+kCarTolerance*0.5 )
561 {
562 maxVal = workMinExtent+(fVoxelNode->GetMaxEquivalentSliceNo()+1)
563 *workNodeWidth;
564 if ( maxVal<=workCoord-kCarTolerance*0.5 )
565 {
566 newNodeNo = fVoxelNode->GetMaxEquivalentSliceNo()+1;
567 newHeader = workHeader;
568 newDistance = (maxVal-localPoint(workHeaderAxis))
569 / localDirection(workHeaderAxis);
570 isNewVoxel = true;
571 newDepth = depth;
572 }
573 }
574 else
575 {
576 newNodeNo = fVoxelNode->GetMinEquivalentSliceNo()-1;
577 newHeader = workHeader;
578 newDistance = (minVal-localPoint(workHeaderAxis))
579 / localDirection(workHeaderAxis);
580 isNewVoxel = true;
581 newDepth = depth;
582 }
583 currentDistance = newDistance;
584 }
585 if (isNewVoxel)
586 {
587 // Compute new voxel & adjust voxel stack
588 //
589 // newNodeNo=Candidate node no at
590 // newDepth =refinement depth of crossed voxel boundary
591 // newHeader=Header for crossed voxel
592 // newDistance=distance to crossed voxel boundary (along the track)
593 //
594 if ( (newNodeNo<0) || (newNodeNo>=newHeader->GetNoSlices()))
595 {
596 // Leaving mother volume
597 //
598 isNewVoxel = false;
599 }
600 else
601 {
602 // Compute intersection point on the least refined
603 // voxel boundary that is hit
604 //
605 voxelPoint = localPoint+localDirection*newDistance;
606 fVoxelNodeNoStack[newDepth] = newNodeNo;
607 fVoxelDepth = newDepth;
608 newVoxelNode = 0;
609 while ( !newVoxelNode )
610 {
611 newProxy = newHeader->GetSlice(newNodeNo);
612 if (newProxy->IsNode())
613 {
614 newVoxelNode = newProxy->GetNode();
615 }
616 else
617 {
618 fVoxelDepth++;
619 newHeader = newProxy->GetHeader();
620 newHeaderAxis = newHeader->GetAxis();
621 newHeaderNoSlices = newHeader->GetNoSlices();
622 newHeaderMin = newHeader->GetMinExtent();
623 newHeaderNodeWidth = (newHeader->GetMaxExtent()-newHeaderMin)
624 / newHeaderNoSlices;
625 newNodeNo = G4int( (voxelPoint(newHeaderAxis)-newHeaderMin)
626 / newHeaderNodeWidth );
627 // Rounding protection
628 //
629 if ( newNodeNo<0 )
630 {
631 newNodeNo=0;
632 }
633 else if ( newNodeNo>=newHeaderNoSlices )
634 {
635 newNodeNo = newHeaderNoSlices-1;
636 }
637 // Stack info for stepping
638 //
639 fVoxelAxisStack[fVoxelDepth] = newHeaderAxis;
640 fVoxelNoSlicesStack[fVoxelDepth] = newHeaderNoSlices;
641 fVoxelSliceWidthStack[fVoxelDepth] = newHeaderNodeWidth;
642 fVoxelNodeNoStack[fVoxelDepth] = newNodeNo;
643 fVoxelHeaderStack[fVoxelDepth] = newHeader;
644 }
645 }
646 fVoxelNode = newVoxelNode;
647 }
648 }
649 return isNewVoxel;
650}
651
652// ********************************************************************
653// ComputeSafety
654//
655// Calculates the isotropic distance to the nearest boundary from the
656// specified point in the local coordinate system.
657// The localpoint utilised must be within the current volume.
658// ********************************************************************
659//
660G4double
661G4VoxelNavigation::ComputeSafety(const G4ThreeVector& localPoint,
662 const G4NavigationHistory& history,
663 const G4double )
664{
665 G4VPhysicalVolume *motherPhysical, *samplePhysical;
666 G4LogicalVolume *motherLogical;
667 G4VSolid *motherSolid;
668 G4double motherSafety, ourSafety;
669 G4int localNoDaughters, sampleNo;
670 G4SmartVoxelNode *curVoxelNode;
671 G4int curNoVolumes, contentNo;
672 G4double voxelSafety;
673
674 motherPhysical = history.GetTopVolume();
675 motherLogical = motherPhysical->GetLogicalVolume();
676 motherSolid = motherLogical->GetSolid();
677
678 //
679 // Compute mother safety
680 //
681
682 motherSafety = motherSolid->DistanceToOut(localPoint);
683 ourSafety = motherSafety; // Working isotropic safety
684
685#ifdef G4VERBOSE
686 if(( fCheck ) && ( fVerbose == 1 ))
687 {
688 G4cout << "*** G4VoxelNavigation::ComputeSafety(): ***" << G4endl
689 << " Invoked DistanceToOut(p) for mother solid: "
690 << motherSolid->GetName()
691 << ". Solid replied: " << motherSafety << G4endl
692 << " For local point p: " << localPoint
693 << ", to be considered as 'mother safety'." << G4endl;
694 }
695#endif
696 //
697 // Compute daughter safeties
698 //
699
700 localNoDaughters = motherLogical->GetNoDaughters();
701
702 // Look only inside the current Voxel only (in the first version).
703 //
704 curVoxelNode = fVoxelNode;
705 curNoVolumes = curVoxelNode->GetNoContained();
706
707 for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
708 {
709 sampleNo = curVoxelNode->GetVolume(contentNo);
710 samplePhysical = motherLogical->GetDaughter(sampleNo);
711
712 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
713 samplePhysical->GetTranslation());
714 sampleTf.Invert();
715 const G4ThreeVector samplePoint =
716 sampleTf.TransformPoint(localPoint);
717 const G4VSolid *sampleSolid =
718 samplePhysical->GetLogicalVolume()->GetSolid();
719 G4double sampleSafety = sampleSolid->DistanceToIn(samplePoint);
720 if ( sampleSafety<ourSafety )
721 {
722 ourSafety = sampleSafety;
723 }
724#ifdef G4VERBOSE
725 if(( fCheck ) && ( fVerbose == 1 ))
726 {
727 G4cout << "*** G4VoxelNavigation::ComputeSafety(): ***" << G4endl
728 << " Invoked DistanceToIn(p) for daughter solid: "
729 << sampleSolid->GetName()
730 << ". Solid replied: " << sampleSafety << G4endl
731 << " For local point p: " << samplePoint
732 << ", to be considered as 'daughter safety'." << G4endl;
733 }
734#endif
735 }
736 voxelSafety = ComputeVoxelSafety(localPoint);
737 if ( voxelSafety<ourSafety )
738 {
739 ourSafety = voxelSafety;
740 }
741 return ourSafety;
742}
Note: See TracBrowser for help on using the repository browser.