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

Last change on this file since 1200 was 921, checked in by garnier, 17 years ago

en test de gl2ps. Problemes de libraries
File size: 28.9 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.9 2008/11/14 18:26:35 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-02-cand-01 $
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 G4double safetyIn= -1, safetyOut= -1; // Set to invalid values
269 G4double newDistIn= -1, newDistOut= -1;
270 if( insideIntPt != kInside )
271 {
272 safetyIn= sampleSolid->DistanceToIn(intersectionPoint);
273 newDistIn= sampleSolid->DistanceToIn(intersectionPoint,
274 sampleDirection);
275 }
276 if( insideIntPt != kOutside )
277 {
278 safetyOut= sampleSolid->DistanceToOut(intersectionPoint);
279 newDistOut= sampleSolid->DistanceToOut(intersectionPoint,
280 sampleDirection);
281 }
282 if( insideIntPt != kSurface )
283 {
284 G4int oldcoutPrec = G4cout.precision(16);
285 G4cout << "WARNING - G4VoxelNavigation::ComputeStep()"
286 << G4endl
287 << " Inaccurate solid DistanceToIn"
288 << " for solid " << sampleSolid->GetName() << G4endl;
289 G4cout << " Solid gave DistanceToIn = "
290 << sampleStep << " yet returns " << solidResponse
291 << " for this point !" << G4endl;
292 G4cout << " Point = " << intersectionPoint << G4endl;
293 G4cout << " Safety values: " << G4endl;
294 if ( insideIntPt != kInside )
295 {
296 G4cout << " DistanceToIn(p) = " << safetyIn
297 << G4endl;
298 }
299 if ( insideIntPt != kOutside )
300 {
301 G4cout << " DistanceToOut(p) = " << safetyOut
302 << G4endl;
303 }
304 G4Exception("G4VoxelNavigation::ComputeStep()",
305 "InaccurateDistanceToIn", JustWarning,
306 "Conflicting response from Solid.");
307 G4cout.precision(oldcoutPrec);
308 }
309 else
310 {
311 // If it is on the surface, *ensure* that either DistanceToIn
312 // or DistanceToOut returns a finite value ( >= Tolerance).
313 //
314 if( std::max( newDistIn, newDistOut ) <= kCarTolerance )
315 {
316 G4cout << "ERROR - G4VoxelNavigation::ComputeStep()"
317 << G4endl
318 << " Identified point for which the solid "
319 << sampleSolid->GetName() << G4endl
320 << " has MAJOR problem: " << G4endl
321 << " --> Both DistanceToIn(p,v) and DistanceToOut(p,v) "
322 << "return Zero, an equivalent value or negative value."
323 << G4endl;
324 G4cout << " Solid: " << sampleSolid << G4endl;
325 G4cout << " Point p= " << intersectionPoint << G4endl;
326 G4cout << " Direction v= " << sampleDirection << G4endl;
327 G4cout << " DistanceToIn(p,v) = " << newDistIn
328 << G4endl;
329 G4cout << " DistanceToOut(p,v,..) = " << newDistOut
330 << G4endl;
331 G4cout << " Safety values: " << G4endl;
332 G4cout << " DistanceToIn(p) = " << safetyIn
333 << G4endl;
334 G4cout << " DistanceToOut(p) = " << safetyOut
335 << G4endl;
336 G4Exception("G4VoxelNavigation::ComputeStep()",
337 "DistanceToInAndOutAreZero", FatalException,
338 "Zero from both Solid DistanceIn and Out(p,v).");
339 }
340 }
341 }
342#endif
343 }
344 }
345 }
346 }
347 }
348 if (initialNode)
349 {
350 initialNode = false;
351 voxelSafety = ComputeVoxelSafety(localPoint);
352 if ( voxelSafety<ourSafety )
353 {
354 ourSafety = voxelSafety;
355 }
356 if ( currentProposedStepLength<ourSafety )
357 {
358 // Guaranteed physics limited
359 //
360 noStep = false;
361 entering = false;
362 exiting = false;
363 *pBlockedPhysical = 0;
364 ourStep = kInfinity;
365 }
366 else
367 {
368 //
369 // Compute mother intersection if required
370 //
371 if ( motherSafety<=ourStep )
372 {
373 G4double motherStep =
374 motherSolid->DistanceToOut(localPoint,
375 localDirection,
376 true, &validExitNormal, &exitNormal);
377#ifdef G4VERBOSE
378 if ( fCheck )
379 {
380 if(fVerbose == 1)
381 {
382 G4cout << "*** G4VoxelNavigation::ComputeStep(): ***" << G4endl
383 << " Invoked DistanceToOut(p,v,...) for mother solid: "
384 << motherSolid->GetName()
385 << ". Solid replied: " << motherStep << G4endl
386 << " For local point p: " << localPoint << G4endl
387 << " Direction v: " << localDirection
388 << ", to be considered as 'mother step'." << G4endl;
389 }
390 if( ( motherStep < 0.0 ) || ( motherStep >= kInfinity) )
391 {
392 G4int oldPrOut= G4cout.precision(16);
393 G4int oldPrErr= G4cerr.precision(16);
394 G4cerr << "ERROR - G4VoxelNavigation::ComputeStep()" << G4endl
395 << " Problem in Navigation" << G4endl
396 << " Point (local coordinates): "
397 << localPoint << G4endl
398 << " Local Direction: " << localDirection << G4endl
399 << " Solid: " << motherSolid->GetName() << G4endl;
400 motherSolid->DumpInfo();
401 G4Exception("G4VoxelNavigation::ComputeStep()",
402 "PointOutsideCurrentVolume", FatalException,
403 "Current point is outside the current solid !");
404 G4cout.precision(oldPrOut);
405 G4cerr.precision(oldPrErr);
406 }
407 }
408#endif
409 if ( motherStep<=ourStep )
410 {
411 ourStep = motherStep;
412 exiting = true;
413 entering = false;
414 if ( validExitNormal )
415 {
416 const G4RotationMatrix *rot = motherPhysical->GetRotation();
417 if (rot)
418 {
419 exitNormal *= rot->inverse();
420 }
421 }
422 }
423 else
424 {
425 validExitNormal = false;
426 }
427 }
428 }
429 newSafety = ourSafety;
430 }
431 if (noStep)
432 {
433 noStep = LocateNextVoxel(localPoint, localDirection, ourStep);
434 }
435 } // end -while (noStep)- loop
436
437 return ourStep;
438}
439
440// ********************************************************************
441// ComputeVoxelSafety
442//
443// Computes safety from specified point to voxel boundaries
444// using already located point
445// o collected boundaries for most derived level
446// o adjacent boundaries for previous levels
447// ********************************************************************
448//
449G4double
450G4VoxelNavigation::ComputeVoxelSafety(const G4ThreeVector& localPoint) const
451{
452 G4SmartVoxelHeader *curHeader;
453 G4double voxelSafety, curNodeWidth;
454 G4double curNodeOffset, minCurCommonDelta, maxCurCommonDelta;
455 G4int minCurNodeNoDelta, maxCurNodeNoDelta;
456 G4int localVoxelDepth, curNodeNo;
457 EAxis curHeaderAxis;
458
459 localVoxelDepth = fVoxelDepth;
460
461 curHeader = fVoxelHeaderStack[localVoxelDepth];
462 curHeaderAxis = fVoxelAxisStack[localVoxelDepth];
463 curNodeNo = fVoxelNodeNoStack[localVoxelDepth];
464 curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth];
465
466 // Compute linear intersection distance to boundaries of max/min
467 // to collected nodes at current level
468 //
469 curNodeOffset = curNodeNo*curNodeWidth;
470 maxCurNodeNoDelta = fVoxelNode->GetMaxEquivalentSliceNo()-curNodeNo;
471 minCurNodeNoDelta = curNodeNo-fVoxelNode->GetMinEquivalentSliceNo();
472 minCurCommonDelta = localPoint(curHeaderAxis)
473 - curHeader->GetMinExtent() - curNodeOffset;
474 maxCurCommonDelta = curNodeWidth-minCurCommonDelta;
475
476 if ( minCurNodeNoDelta<maxCurNodeNoDelta )
477 {
478 voxelSafety = minCurNodeNoDelta*curNodeWidth;
479 voxelSafety += minCurCommonDelta;
480 }
481 else if (maxCurNodeNoDelta < minCurNodeNoDelta)
482 {
483 voxelSafety = maxCurNodeNoDelta*curNodeWidth;
484 voxelSafety += maxCurCommonDelta;
485 }
486 else // (maxCurNodeNoDelta == minCurNodeNoDelta)
487 {
488 voxelSafety = minCurNodeNoDelta*curNodeWidth;
489 voxelSafety += std::min(minCurCommonDelta,maxCurCommonDelta);
490 }
491
492 // Compute isotropic safety to boundaries of previous levels
493 // [NOT to collected boundaries]
494 //
495 while ( (localVoxelDepth>0) && (voxelSafety>0) )
496 {
497 localVoxelDepth--;
498 curHeader = fVoxelHeaderStack[localVoxelDepth];
499 curHeaderAxis = fVoxelAxisStack[localVoxelDepth];
500 curNodeNo = fVoxelNodeNoStack[localVoxelDepth];
501 curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth];
502 curNodeOffset = curNodeNo*curNodeWidth;
503 minCurCommonDelta = localPoint(curHeaderAxis)
504 - curHeader->GetMinExtent() - curNodeOffset;
505 maxCurCommonDelta = curNodeWidth-minCurCommonDelta;
506
507 if ( minCurCommonDelta<voxelSafety )
508 {
509 voxelSafety = minCurCommonDelta;
510 }
511 if ( maxCurCommonDelta<voxelSafety )
512 {
513 voxelSafety = maxCurCommonDelta;
514 }
515 }
516 if ( voxelSafety<0 )
517 {
518 voxelSafety = 0;
519 }
520
521 return voxelSafety;
522}
523
524// ********************************************************************
525// LocateNextVoxel
526//
527// Finds the next voxel from the current voxel and point
528// in the specified direction
529//
530// Returns false if all voxels considered
531// [current Step ends inside same voxel or leaves all voxels]
532// true otherwise
533// [the information on the next voxel is put into the set of
534// fVoxel* variables & "stacks"]
535// ********************************************************************
536//
537G4bool
538G4VoxelNavigation::LocateNextVoxel(const G4ThreeVector& localPoint,
539 const G4ThreeVector& localDirection,
540 const G4double currentStep)
541{
542 G4SmartVoxelHeader *workHeader=0, *newHeader=0;
543 G4SmartVoxelProxy *newProxy=0;
544 G4SmartVoxelNode *newVoxelNode=0;
545 G4ThreeVector targetPoint, voxelPoint;
546 G4double workNodeWidth, workMinExtent, workCoord;
547 G4double minVal, maxVal, newDistance=0.;
548 G4double newHeaderMin, newHeaderNodeWidth;
549 G4int depth=0, newDepth=0, workNodeNo=0, newNodeNo=0, newHeaderNoSlices=0;
550 EAxis workHeaderAxis, newHeaderAxis;
551 G4bool isNewVoxel=false;
552
553 G4double currentDistance = currentStep;
554
555 // Determine if end of Step within current voxel
556 //
557 for (depth=0; depth<fVoxelDepth; depth++)
558 {
559 targetPoint = localPoint+localDirection*currentDistance;
560 newDistance = currentDistance;
561 workHeader = fVoxelHeaderStack[depth];
562 workHeaderAxis = fVoxelAxisStack[depth];
563 workNodeNo = fVoxelNodeNoStack[depth];
564 workNodeWidth = fVoxelSliceWidthStack[depth];
565 workMinExtent = workHeader->GetMinExtent();
566 workCoord = targetPoint(workHeaderAxis);
567 minVal = workMinExtent+workNodeNo*workNodeWidth;
568
569 if ( minVal<=workCoord+kCarTolerance*0.5 )
570 {
571 maxVal = minVal+workNodeWidth;
572 if ( maxVal<=workCoord-kCarTolerance*0.5 )
573 {
574 // Must consider next voxel
575 //
576 newNodeNo = workNodeNo+1;
577 newHeader = workHeader;
578 newDistance = (maxVal-localPoint(workHeaderAxis))
579 / localDirection(workHeaderAxis);
580 isNewVoxel = true;
581 newDepth = depth;
582 }
583 }
584 else
585 {
586 newNodeNo = workNodeNo-1;
587 newHeader = workHeader;
588 newDistance = (minVal-localPoint(workHeaderAxis))
589 / localDirection(workHeaderAxis);
590 isNewVoxel = true;
591 newDepth = depth;
592 }
593 currentDistance = newDistance;
594 }
595 targetPoint = localPoint+localDirection*currentDistance;
596
597 // Check if end of Step within collected boundaries of current voxel
598 //
599 depth = fVoxelDepth;
600 {
601 workHeader = fVoxelHeaderStack[depth];
602 workHeaderAxis = fVoxelAxisStack[depth];
603 workNodeNo = fVoxelNodeNoStack[depth];
604 workNodeWidth = fVoxelSliceWidthStack[depth];
605 workMinExtent = workHeader->GetMinExtent();
606 workCoord = targetPoint(workHeaderAxis);
607 minVal = workMinExtent+fVoxelNode->GetMinEquivalentSliceNo()*workNodeWidth;
608
609 if ( minVal<=workCoord+kCarTolerance*0.5 )
610 {
611 maxVal = workMinExtent+(fVoxelNode->GetMaxEquivalentSliceNo()+1)
612 *workNodeWidth;
613 if ( maxVal<=workCoord-kCarTolerance*0.5 )
614 {
615 newNodeNo = fVoxelNode->GetMaxEquivalentSliceNo()+1;
616 newHeader = workHeader;
617 newDistance = (maxVal-localPoint(workHeaderAxis))
618 / localDirection(workHeaderAxis);
619 isNewVoxel = true;
620 newDepth = depth;
621 }
622 }
623 else
624 {
625 newNodeNo = fVoxelNode->GetMinEquivalentSliceNo()-1;
626 newHeader = workHeader;
627 newDistance = (minVal-localPoint(workHeaderAxis))
628 / localDirection(workHeaderAxis);
629 isNewVoxel = true;
630 newDepth = depth;
631 }
632 currentDistance = newDistance;
633 }
634 if (isNewVoxel)
635 {
636 // Compute new voxel & adjust voxel stack
637 //
638 // newNodeNo=Candidate node no at
639 // newDepth =refinement depth of crossed voxel boundary
640 // newHeader=Header for crossed voxel
641 // newDistance=distance to crossed voxel boundary (along the track)
642 //
643 if ( (newNodeNo<0) || (newNodeNo>=newHeader->GetNoSlices()))
644 {
645 // Leaving mother volume
646 //
647 isNewVoxel = false;
648 }
649 else
650 {
651 // Compute intersection point on the least refined
652 // voxel boundary that is hit
653 //
654 voxelPoint = localPoint+localDirection*newDistance;
655 fVoxelNodeNoStack[newDepth] = newNodeNo;
656 fVoxelDepth = newDepth;
657 newVoxelNode = 0;
658 while ( !newVoxelNode )
659 {
660 newProxy = newHeader->GetSlice(newNodeNo);
661 if (newProxy->IsNode())
662 {
663 newVoxelNode = newProxy->GetNode();
664 }
665 else
666 {
667 fVoxelDepth++;
668 newHeader = newProxy->GetHeader();
669 newHeaderAxis = newHeader->GetAxis();
670 newHeaderNoSlices = newHeader->GetNoSlices();
671 newHeaderMin = newHeader->GetMinExtent();
672 newHeaderNodeWidth = (newHeader->GetMaxExtent()-newHeaderMin)
673 / newHeaderNoSlices;
674 newNodeNo = G4int( (voxelPoint(newHeaderAxis)-newHeaderMin)
675 / newHeaderNodeWidth );
676 // Rounding protection
677 //
678 if ( newNodeNo<0 )
679 {
680 newNodeNo=0;
681 }
682 else if ( newNodeNo>=newHeaderNoSlices )
683 {
684 newNodeNo = newHeaderNoSlices-1;
685 }
686 // Stack info for stepping
687 //
688 fVoxelAxisStack[fVoxelDepth] = newHeaderAxis;
689 fVoxelNoSlicesStack[fVoxelDepth] = newHeaderNoSlices;
690 fVoxelSliceWidthStack[fVoxelDepth] = newHeaderNodeWidth;
691 fVoxelNodeNoStack[fVoxelDepth] = newNodeNo;
692 fVoxelHeaderStack[fVoxelDepth] = newHeader;
693 }
694 }
695 fVoxelNode = newVoxelNode;
696 }
697 }
698 return isNewVoxel;
699}
700
701// ********************************************************************
702// ComputeSafety
703//
704// Calculates the isotropic distance to the nearest boundary from the
705// specified point in the local coordinate system.
706// The localpoint utilised must be within the current volume.
707// ********************************************************************
708//
709G4double
710G4VoxelNavigation::ComputeSafety(const G4ThreeVector& localPoint,
711 const G4NavigationHistory& history,
712 const G4double )
713{
714 G4VPhysicalVolume *motherPhysical, *samplePhysical;
715 G4LogicalVolume *motherLogical;
716 G4VSolid *motherSolid;
717 G4double motherSafety, ourSafety;
718 G4int localNoDaughters, sampleNo;
719 G4SmartVoxelNode *curVoxelNode;
720 G4int curNoVolumes, contentNo;
721 G4double voxelSafety;
722
723 motherPhysical = history.GetTopVolume();
724 motherLogical = motherPhysical->GetLogicalVolume();
725 motherSolid = motherLogical->GetSolid();
726
727 //
728 // Compute mother safety
729 //
730
731 motherSafety = motherSolid->DistanceToOut(localPoint);
732 ourSafety = motherSafety; // Working isotropic safety
733
734#ifdef G4VERBOSE
735 if(( fCheck ) && ( fVerbose == 1 ))
736 {
737 G4cout << "*** G4VoxelNavigation::ComputeSafety(): ***" << G4endl
738 << " Invoked DistanceToOut(p) for mother solid: "
739 << motherSolid->GetName()
740 << ". Solid replied: " << motherSafety << G4endl
741 << " For local point p: " << localPoint
742 << ", to be considered as 'mother safety'." << G4endl;
743 }
744#endif
745 //
746 // Compute daughter safeties
747 //
748
749 localNoDaughters = motherLogical->GetNoDaughters();
750
751 // Look only inside the current Voxel only (in the first version).
752 //
753 curVoxelNode = fVoxelNode;
754 curNoVolumes = curVoxelNode->GetNoContained();
755
756 for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
757 {
758 sampleNo = curVoxelNode->GetVolume(contentNo);
759 samplePhysical = motherLogical->GetDaughter(sampleNo);
760
761 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
762 samplePhysical->GetTranslation());
763 sampleTf.Invert();
764 const G4ThreeVector samplePoint =
765 sampleTf.TransformPoint(localPoint);
766 const G4VSolid *sampleSolid =
767 samplePhysical->GetLogicalVolume()->GetSolid();
768 G4double sampleSafety = sampleSolid->DistanceToIn(samplePoint);
769 if ( sampleSafety<ourSafety )
770 {
771 ourSafety = sampleSafety;
772 }
773#ifdef G4VERBOSE
774 if(( fCheck ) && ( fVerbose == 1 ))
775 {
776 G4cout << "*** G4VoxelNavigation::ComputeSafety(): ***" << G4endl
777 << " Invoked DistanceToIn(p) for daughter solid: "
778 << sampleSolid->GetName()
779 << ". Solid replied: " << sampleSafety << G4endl
780 << " For local point p: " << samplePoint
781 << ", to be considered as 'daughter safety'." << G4endl;
782 }
783#endif
784 }
785 voxelSafety = ComputeVoxelSafety(localPoint);
786 if ( voxelSafety<ourSafety )
787 {
788 ourSafety = voxelSafety;
789 }
790 return ourSafety;
791}
Note: See TracBrowser for help on using the repository browser.