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

Last change on this file since 836 was 831, checked in by garnier, 16 years ago

import all except CVS
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:  $
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.