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

Last change on this file since 1058 was 921, checked in by garnier, 15 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.