Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/navigation/src/G4VoxelSafety.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/navigation/src/G4VoxelSafety.cc (Version 11.3.0) and /geometry/navigation/src/G4VoxelSafety.cc (Version 11.1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27 //  Author:  John Apostolakis                      27 //  Author:  John Apostolakis
 28 //  First version:  31 May 2010                    28 //  First version:  31 May 2010
 29 //                                                 29 // 
 30 // -------------------------------------------     30 // --------------------------------------------------------------------
 31 #include "G4VoxelSafety.hh"                        31 #include "G4VoxelSafety.hh"
 32                                                    32 
 33 #include "G4GeometryTolerance.hh"                  33 #include "G4GeometryTolerance.hh"
 34                                                    34 
 35 #include "G4SmartVoxelProxy.hh"                    35 #include "G4SmartVoxelProxy.hh"
 36 #include "G4SmartVoxelNode.hh"                     36 #include "G4SmartVoxelNode.hh"
 37 #include "G4SmartVoxelHeader.hh"                   37 #include "G4SmartVoxelHeader.hh"
 38                                                    38 
 39 // *******************************************     39 // ********************************************************************
 40 // Constructor                                     40 // Constructor
 41 //     - copied from G4VoxelNavigation  (1st v     41 //     - copied from G4VoxelNavigation  (1st version)
 42 // *******************************************     42 // ********************************************************************
 43 //                                                 43 //
 44 G4VoxelSafety::G4VoxelSafety()                     44 G4VoxelSafety::G4VoxelSafety()
 45   : fVoxelAxisStack(kNavigatorVoxelStackMax,kX <<  45   : fBlockList(),
                                                   >>  46     fVoxelAxisStack(kNavigatorVoxelStackMax,kXAxis),
 46     fVoxelNoSlicesStack(kNavigatorVoxelStackMa     47     fVoxelNoSlicesStack(kNavigatorVoxelStackMax,0),
 47     fVoxelSliceWidthStack(kNavigatorVoxelStack     48     fVoxelSliceWidthStack(kNavigatorVoxelStackMax,0.),
 48     fVoxelNodeNoStack(kNavigatorVoxelStackMax,     49     fVoxelNodeNoStack(kNavigatorVoxelStackMax,0),
 49     fVoxelHeaderStack(kNavigatorVoxelStackMax,     50     fVoxelHeaderStack(kNavigatorVoxelStackMax,(G4SmartVoxelHeader*)nullptr)
 50 {                                                  51 {
 51   kCarTolerance = G4GeometryTolerance::GetInst     52   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
 52 }                                                  53 }
 53                                                    54 
 54 // *******************************************     55 // ********************************************************************
 55 // Destructor                                      56 // Destructor
 56 // *******************************************     57 // ********************************************************************
 57 //                                                 58 //
 58 G4VoxelSafety::~G4VoxelSafety() = default;     <<  59 G4VoxelSafety::~G4VoxelSafety()
                                                   >>  60 {
                                                   >>  61 }
 59                                                    62 
 60 // *******************************************     63 // ********************************************************************
 61 // ComputeSafety                                   64 // ComputeSafety
 62 //                                                 65 //
 63 // Calculates the isotropic distance to the ne     66 // Calculates the isotropic distance to the nearest boundary from the
 64 // specified point in the local coordinate sys     67 // specified point in the local coordinate system. 
 65 // The localpoint utilised must be within the      68 // The localpoint utilised must be within the current volume.
 66 // *******************************************     69 // ********************************************************************
 67 //                                                 70 //
 68 G4double                                           71 G4double
 69 G4VoxelSafety::ComputeSafety(const G4ThreeVect     72 G4VoxelSafety::ComputeSafety(const G4ThreeVector& localPoint,
 70                              const G4VPhysical     73                              const G4VPhysicalVolume& currentPhysical,
 71                                    G4double ma     74                                    G4double maxLength)
 72 {                                                  75 {
 73   G4LogicalVolume *motherLogical;                  76   G4LogicalVolume *motherLogical;
 74   G4VSolid *motherSolid;                           77   G4VSolid *motherSolid;
 75   G4SmartVoxelHeader *motherVoxelHeader;           78   G4SmartVoxelHeader *motherVoxelHeader;
 76   G4double motherSafety, ourSafety;                79   G4double motherSafety, ourSafety;
 77   G4int localNoDaughters;                          80   G4int localNoDaughters;
 78   G4double daughterSafety;                         81   G4double daughterSafety;
 79                                                    82 
 80   motherLogical = currentPhysical.GetLogicalVo     83   motherLogical = currentPhysical.GetLogicalVolume();
 81   fpMotherLogical= motherLogical;   // For use     84   fpMotherLogical= motherLogical;   // For use by the other methods
 82   motherSolid = motherLogical->GetSolid();         85   motherSolid = motherLogical->GetSolid();
 83   motherVoxelHeader = motherLogical->GetVoxelH     86   motherVoxelHeader = motherLogical->GetVoxelHeader();
 84                                                    87 
 85 #ifdef G4VERBOSE                                   88 #ifdef G4VERBOSE  
 86   if( fVerbose > 0 )                               89   if( fVerbose > 0 )
 87   {                                                90   { 
 88     G4cout << "*** G4VoxelSafety::ComputeSafet     91     G4cout << "*** G4VoxelSafety::ComputeSafety(): ***" << G4endl; 
 89   }                                                92   }
 90 #endif                                             93 #endif
 91                                                    94 
 92   // Check that point is inside mother volume      95   // Check that point is inside mother volume
 93   //                                               96   //
 94   EInside  insideMother = motherSolid->Inside(     97   EInside  insideMother = motherSolid->Inside(localPoint); 
 95   if( insideMother != kInside  )                   98   if( insideMother != kInside  )
 96   {                                                99   { 
 97 #ifdef G4DEBUG_NAVIGATION                         100 #ifdef G4DEBUG_NAVIGATION
 98     if( insideMother == kOutside )                101     if( insideMother == kOutside )
 99     {                                             102     {
100       std::ostringstream message;                 103       std::ostringstream message;
101       message << "Safety method called for loc    104       message << "Safety method called for location outside current Volume."
102               << G4endl                           105               << G4endl
103               << "Location for safety is Outsi    106               << "Location for safety is Outside this volume. " << G4endl
104               << "The approximate distance to     107               << "The approximate distance to the solid "
105               << "(safety from outside) is: "     108               << "(safety from outside) is: " 
106               << motherSolid->DistanceToIn( lo    109               << motherSolid->DistanceToIn( localPoint ) << G4endl;
107       message << "  Problem occurred with phys    110       message << "  Problem occurred with physical volume: "
108               << " Name: " << currentPhysical.    111               << " Name: " << currentPhysical.GetName()
109               << " Copy No: " << currentPhysic    112               << " Copy No: " << currentPhysical.GetCopyNo() << G4endl
110               << "    Local Point = " << local    113               << "    Local Point = " << localPoint << G4endl;
111       message << "  Description of solid: " <<    114       message << "  Description of solid: " << G4endl
112               << *motherSolid << G4endl;          115               << *motherSolid << G4endl;
113       G4Exception("G4VoxelSafety::ComputeSafet    116       G4Exception("G4VoxelSafety::ComputeSafety()", "GeomNav0003",
114                   FatalException, message);       117                   FatalException, message);
115     }                                             118     }
116 #endif                                            119 #endif
117     return 0.0;                                   120     return 0.0;
118   }                                               121   }   
119                                                   122 
120   //  First limit:  mother safety - distance t    123   //  First limit:  mother safety - distance to outer boundaries
121   //                                              124   //
122   motherSafety = motherSolid->DistanceToOut(lo    125   motherSafety = motherSolid->DistanceToOut(localPoint);
123   ourSafety = motherSafety;                 //    126   ourSafety = motherSafety;                 // Working isotropic safety
124                                                   127    
125 #ifdef G4VERBOSE                                  128 #ifdef G4VERBOSE
126   if(( fCheck ) ) // && ( fVerbose == 1 ))        129   if(( fCheck ) ) // && ( fVerbose == 1 ))
127   {                                               130   {
128     G4cout << "    Invoked DistanceToOut(p) fo    131     G4cout << "    Invoked DistanceToOut(p) for mother solid: "
129            << motherSolid->GetName()              132            << motherSolid->GetName()
130            << ". Solid replied: " << motherSaf    133            << ". Solid replied: " << motherSafety << G4endl
131            << "    For local point p: " << loc    134            << "    For local point p: " << localPoint
132            << ", to be considered as 'mother s    135            << ", to be considered as 'mother safety'." << G4endl;
133   }                                               136   }
134 #endif                                            137 #endif
135   localNoDaughters = (G4int)motherLogical->Get    138   localNoDaughters = (G4int)motherLogical->GetNoDaughters();
136                                                   139 
137   fBlockList.Enlarge(localNoDaughters);           140   fBlockList.Enlarge(localNoDaughters);
138   fBlockList.Reset();                             141   fBlockList.Reset();
139                                                   142 
140   fVoxelDepth = -1;  // Resets the depth -- mu    143   fVoxelDepth = -1;  // Resets the depth -- must be done for next method
141   daughterSafety= SafetyForVoxelHeader(motherV    144   daughterSafety= SafetyForVoxelHeader(motherVoxelHeader, localPoint, maxLength,
142                                        current    145                                        currentPhysical, 0.0, ourSafety);
143   ourSafety= std::min( motherSafety, daughterS    146   ourSafety= std::min( motherSafety, daughterSafety ); 
144                                                   147 
145   return ourSafety;                               148   return ourSafety;
146 }                                                 149 }
147                                                   150 
148 // *******************************************    151 // ********************************************************************
149 // SafetyForVoxelNode                             152 // SafetyForVoxelNode
150 //                                                153 //
151 // Calculate the safety for volumes included i    154 // Calculate the safety for volumes included in current Voxel Node
152 // *******************************************    155 // ********************************************************************
153 //                                                156 // 
154 G4double                                          157 G4double
155 G4VoxelSafety::SafetyForVoxelNode( const G4Sma    158 G4VoxelSafety::SafetyForVoxelNode( const G4SmartVoxelNode* curVoxelNode,
156                                    const G4Thr    159                                    const G4ThreeVector& localPoint )  
157 {                                                 160 {
158    G4double ourSafety = DBL_MAX;                  161    G4double ourSafety = DBL_MAX;
159                                                   162 
160    G4long curNoVolumes, contentNo;                163    G4long curNoVolumes, contentNo;
161    G4int sampleNo;                                164    G4int sampleNo;
162    G4VPhysicalVolume* samplePhysical;             165    G4VPhysicalVolume* samplePhysical;
163                                                   166 
164    G4double      sampleSafety = 0.0;              167    G4double      sampleSafety = 0.0; 
165    G4ThreeVector samplePoint;                     168    G4ThreeVector samplePoint;
166    G4VSolid*     ptrSolid = nullptr;              169    G4VSolid*     ptrSolid = nullptr;
167                                                   170 
168    curNoVolumes = curVoxelNode->GetNoContained    171    curNoVolumes = curVoxelNode->GetNoContained();
169                                                   172 
170    for ( contentNo=curNoVolumes-1; contentNo>=    173    for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
171    {                                              174    {
172       sampleNo = curVoxelNode->GetVolume((G4in    175       sampleNo = curVoxelNode->GetVolume((G4int)contentNo);
173       if ( !fBlockList.IsBlocked(sampleNo) )      176       if ( !fBlockList.IsBlocked(sampleNo) ) 
174       {                                           177       { 
175         fBlockList.BlockVolume(sampleNo);         178         fBlockList.BlockVolume(sampleNo);
176                                                   179 
177         samplePhysical = fpMotherLogical->GetD    180         samplePhysical = fpMotherLogical->GetDaughter(sampleNo);
178         G4AffineTransform sampleTf(samplePhysi    181         G4AffineTransform sampleTf(samplePhysical->GetRotation(),
179                                    samplePhysi    182                                    samplePhysical->GetTranslation());
180         sampleTf.Invert();                        183         sampleTf.Invert();
181         samplePoint = sampleTf.TransformPoint(    184         samplePoint = sampleTf.TransformPoint(localPoint);
182         ptrSolid = samplePhysical->GetLogicalV    185         ptrSolid = samplePhysical->GetLogicalVolume()->GetSolid();
183                                                   186 
184         sampleSafety = ptrSolid->DistanceToIn(    187         sampleSafety = ptrSolid->DistanceToIn(samplePoint);
185         ourSafety = std::min( sampleSafety, ou    188         ourSafety = std::min( sampleSafety, ourSafety ); 
186 #ifdef G4VERBOSE                                  189 #ifdef G4VERBOSE
187         if(( fCheck ) && ( fVerbose == 1 ))       190         if(( fCheck ) && ( fVerbose == 1 ))
188         {                                         191         {
189           G4cout << "*** G4VoxelSafety::Safety    192           G4cout << "*** G4VoxelSafety::SafetyForVoxelNode(): ***" << G4endl
190                  << "    Invoked DistanceToIn(    193                  << "    Invoked DistanceToIn(p) for daughter solid: "
191                  << ptrSolid->GetName()           194                  << ptrSolid->GetName()
192                  << ". Solid replied: " << sam    195                  << ". Solid replied: " << sampleSafety << G4endl
193                  << "    For local point p: "     196                  << "    For local point p: " << samplePoint
194                  << ", to be considered as 'da    197                  << ", to be considered as 'daughter safety'." << G4endl;
195         }                                         198         }
196 #endif                                            199 #endif
197       }                                           200       }
198    }  // end for contentNo                        201    }  // end for contentNo
199                                                   202 
200    return ourSafety;                              203    return ourSafety; 
201 }                                                 204 }
202                                                   205 
203 // *******************************************    206 // ********************************************************************
204 // SafetyForVoxelHeader                           207 // SafetyForVoxelHeader
205 //                                                208 //
206 // Cycles through levels of headers to process    209 // Cycles through levels of headers to process each node level
207 // Obtained by modifying VoxelLocate (to cycle    210 // Obtained by modifying VoxelLocate (to cycle through Node Headers)
208 // *******************************************    211 // *********************************************************************
209 //                                                212 //
210 G4double                                          213 G4double
211 G4VoxelSafety::SafetyForVoxelHeader( const G4S    214 G4VoxelSafety::SafetyForVoxelHeader( const G4SmartVoxelHeader* pHeader,
212                                      const G4T    215                                      const G4ThreeVector& localPoint,
213                                            G4d    216                                            G4double maxLength,
214                                      const G4V    217                                      const G4VPhysicalVolume&  currentPhysical, //Debug
215                                            G4d    218                                            G4double distUpperDepth_Sq,
216                                            G4d    219                                            G4double previousMinSafety
217                                    )              220                                    )
218 {                                                 221 {
219   const G4SmartVoxelHeader* const targetVoxelH    222   const G4SmartVoxelHeader* const targetVoxelHeader = pHeader;
220   G4SmartVoxelNode* targetVoxelNode = nullptr;    223   G4SmartVoxelNode* targetVoxelNode = nullptr;
221                                                   224 
222   const G4SmartVoxelProxy* sampleProxy;           225   const G4SmartVoxelProxy* sampleProxy;
223   EAxis    targetHeaderAxis;                      226   EAxis    targetHeaderAxis;
224   G4double targetHeaderMin, targetHeaderMax, t    227   G4double targetHeaderMin, targetHeaderMax, targetHeaderNodeWidth;
225   G4int    targetHeaderNoSlices;                  228   G4int    targetHeaderNoSlices;
226   G4int    targetNodeNo;                          229   G4int    targetNodeNo;
227                                                   230 
228   G4double minSafety = previousMinSafety;         231   G4double minSafety = previousMinSafety;
229   G4double ourSafety = DBL_MAX;                   232   G4double ourSafety = DBL_MAX;
230   unsigned int checkedNum= 0;                     233   unsigned int checkedNum= 0;
231                                                   234   
232   ++fVoxelDepth;                                  235   ++fVoxelDepth;
233   // fVoxelDepth  set by ComputeSafety or prev    236   // fVoxelDepth  set by ComputeSafety or previous level call
234                                                   237 
235   targetHeaderAxis =      targetVoxelHeader->G    238   targetHeaderAxis =      targetVoxelHeader->GetAxis();
236   targetHeaderNoSlices =  (G4int)targetVoxelHe    239   targetHeaderNoSlices =  (G4int)targetVoxelHeader->GetNoSlices();
237   targetHeaderMin =       targetVoxelHeader->G    240   targetHeaderMin =       targetVoxelHeader->GetMinExtent();
238   targetHeaderMax =       targetVoxelHeader->G    241   targetHeaderMax =       targetVoxelHeader->GetMaxExtent();
239                                                   242 
240   targetHeaderNodeWidth = (targetHeaderMax-tar    243   targetHeaderNodeWidth = (targetHeaderMax-targetHeaderMin)
241                           / targetHeaderNoSlic    244                           / targetHeaderNoSlices;
242                                                   245 
243   G4double localCrd = localPoint(targetHeaderA    246   G4double localCrd = localPoint(targetHeaderAxis);
244                                                   247 
245   const auto  candNodeNo = G4int( (localCrd-ta << 248   const G4int candNodeNo = G4int( (localCrd-targetHeaderMin)
246                                  / targetHeade    249                                  / targetHeaderNodeWidth );
247   // Ensure that it is between 0 and targetHea    250   // Ensure that it is between 0 and targetHeader->GetMaxExtent() - 1
248                                                   251 
249 #ifdef G4DEBUG_VOXELISATION                       252 #ifdef G4DEBUG_VOXELISATION  
250   if( candNodeNo < 0 || candNodeNo > targetHea    253   if( candNodeNo < 0 || candNodeNo > targetHeaderNoSlices-1 )
251   {                                               254   {
252      G4ExceptionDescription ed;                   255      G4ExceptionDescription ed;
253      ed << " Potential ERROR."                    256      ed << " Potential ERROR."
254         << " Point is outside range of Voxel i    257         << " Point is outside range of Voxel in current coordinate" << G4endl;
255      ed << "  Node number of point " << localP    258      ed << "  Node number of point " << localPoint
256         << "is outside the range. " << G4endl;    259         << "is outside the range. " << G4endl;
257      ed << "    Voxel node Num= " << candNodeN    260      ed << "    Voxel node Num= " << candNodeNo << " versus  minimum= " << 0
258         << " and maximum= " << targetHeaderNoS    261         << " and maximum= " << targetHeaderNoSlices-1 << G4endl;
259      ed << "    Axis = " << targetHeaderAxis      262      ed << "    Axis = " << targetHeaderAxis
260         << "  No of slices = " << targetHeader    263         << "  No of slices = " << targetHeaderNoSlices << G4endl;
261      ed << "    Local coord = " << localCrd       264      ed << "    Local coord = " << localCrd
262         << "  Voxel Min = " << targetHeaderMin    265         << "  Voxel Min = " << targetHeaderMin
263         <<            " Max = " << targetHeade    266         <<            " Max = " << targetHeaderMax << G4endl;
264      G4LogicalVolume *pLogical= currentPhysica    267      G4LogicalVolume *pLogical= currentPhysical.GetLogicalVolume();
265      ed << "  Current volume (physical) = " <<    268      ed << "  Current volume (physical) = " << currentPhysical.GetName()
266         << "   (logical) = " << pLogical->GetN    269         << "   (logical) = " << pLogical->GetName()     << G4endl;
267      G4VSolid* pSolid= pLogical->GetSolid();      270      G4VSolid* pSolid= pLogical->GetSolid();
268      ed << "  Solid type = " << pSolid->GetEnt    271      ed << "  Solid type = " << pSolid->GetEntityType() << G4endl;
269      ed << *pSolid << G4endl;                     272      ed << *pSolid << G4endl;
270      G4Exception("G4VoxelSafety::SafetyForVoxe    273      G4Exception("G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav1003",
271                  JustWarning, ed,                 274                  JustWarning, ed,
272                  "Point is outside range of Vo    275                  "Point is outside range of Voxel in current coordinate");
273   }                                               276   }
274 #endif                                            277 #endif
275                                                   278   
276   const G4int pointNodeNo =                       279   const G4int pointNodeNo =
277               std::max( 0, std::min( candNodeN    280               std::max( 0, std::min( candNodeNo, targetHeaderNoSlices-1 ) );
278                                                   281   
279 #ifdef G4VERBOSE                                  282 #ifdef G4VERBOSE  
280   if( fVerbose > 2 )                              283   if( fVerbose > 2 )
281   {                                               284   { 
282     G4cout << G4endl;                             285     G4cout << G4endl;
283     G4cout << "**** G4VoxelSafety::SafetyForVo    286     G4cout << "**** G4VoxelSafety::SafetyForVoxelHeader  " << G4endl;
284     G4cout << "  Called at Depth     = " << fV    287     G4cout << "  Called at Depth     = " << fVoxelDepth;
285     G4cout << " Calculated pointNodeNo= " << p    288     G4cout << " Calculated pointNodeNo= " << pointNodeNo
286            << "  from position= " <<  localPoi    289            << "  from position= " <<  localPoint(targetHeaderAxis)
287            << "  min= "    << targetHeaderMin     290            << "  min= "    << targetHeaderMin
288            << "  max= "    << targetHeaderMax     291            << "  max= "    << targetHeaderMax
289            << "  width= "  << targetHeaderNode    292            << "  width= "  << targetHeaderNodeWidth 
290            << "  no-slices= " << targetHeaderN    293            << "  no-slices= " << targetHeaderNoSlices
291            << "  axis=  "  << targetHeaderAxis    294            << "  axis=  "  << targetHeaderAxis   << G4endl;
292   }                                               295   }
293   else if (fVerbose == 1)                         296   else if (fVerbose == 1)
294   {                                               297   {
295     G4cout << " VoxelSafety: Depth  = " << fVo    298     G4cout << " VoxelSafety: Depth  = " << fVoxelDepth 
296            << " Number of Slices = " << target    299            << " Number of Slices = " << targetHeaderNoSlices
297            << " Header (address) = " << target    300            << " Header (address) = " << targetVoxelHeader  << G4endl;
298   }                                               301   }
299 #endif                                            302 #endif
300                                                   303 
301   // Stack info for stepping                      304   // Stack info for stepping
302   //                                              305   //
303   fVoxelAxisStack[fVoxelDepth] = targetHeaderA    306   fVoxelAxisStack[fVoxelDepth] = targetHeaderAxis;
304   fVoxelNoSlicesStack[fVoxelDepth] = targetHea    307   fVoxelNoSlicesStack[fVoxelDepth] = targetHeaderNoSlices;
305   fVoxelSliceWidthStack[fVoxelDepth] = targetH    308   fVoxelSliceWidthStack[fVoxelDepth] = targetHeaderNodeWidth;
306                                                   309 
307   fVoxelHeaderStack[fVoxelDepth] = pHeader;       310   fVoxelHeaderStack[fVoxelDepth] = pHeader;
308                                                   311 
309   G4int   trialUp = -1,     trialDown = -1;       312   G4int   trialUp = -1,     trialDown = -1;
310   G4double distUp = DBL_MAX, distDown = DBL_MA    313   G4double distUp = DBL_MAX, distDown = DBL_MAX;
311                                                   314 
312   // Using Equivalent voxels - this is pre-ini    315   // Using Equivalent voxels - this is pre-initialisation only
313   //                                              316   //
314   G4int nextUp =   pointNodeNo+1;                 317   G4int nextUp =   pointNodeNo+1;
315   G4int nextDown = pointNodeNo-1;                 318   G4int nextDown = pointNodeNo-1;
316                                                   319 
317   G4int    nextNodeNo = pointNodeNo;              320   G4int    nextNodeNo = pointNodeNo;
318   G4double distAxis;  // Distance in current A    321   G4double distAxis;  // Distance in current Axis
319   distAxis = 0.0;  // Starting in node contain    322   distAxis = 0.0;  // Starting in node containing local Coordinate
320                                                   323 
321   G4bool nextIsInside = false;                    324   G4bool nextIsInside = false;
322                                                   325 
323   G4double distMaxInterest= std::min( previous    326   G4double distMaxInterest= std::min( previousMinSafety, maxLength);
324     // We will not look beyond this distance.     327     // We will not look beyond this distance.
325     // This distance will be updated to reflec    328     // This distance will be updated to reflect the
326     // max ( minSafety, maxLength ) at each st    329     // max ( minSafety, maxLength ) at each step
327                                                   330 
328   targetNodeNo = pointNodeNo;                     331   targetNodeNo = pointNodeNo;
329   do                                              332   do
330   {                                               333   {
331      G4double nodeSafety = DBL_MAX, headerSafe    334      G4double nodeSafety = DBL_MAX, headerSafety = DBL_MAX;
332      fVoxelNodeNoStack[fVoxelDepth] = targetNo    335      fVoxelNodeNoStack[fVoxelDepth] = targetNodeNo;
333                                                   336     
334      ++checkedNum;                                337      ++checkedNum; 
335                                                   338     
336      sampleProxy = targetVoxelHeader->GetSlice    339      sampleProxy = targetVoxelHeader->GetSlice(targetNodeNo);
337                                                   340 
338 #ifdef G4DEBUG_NAVIGATION                         341 #ifdef G4DEBUG_NAVIGATION
339      assert( sampleProxy != 0);                   342      assert( sampleProxy != 0);
340      if( fVerbose > 2 )                           343      if( fVerbose > 2 ) 
341      {                                            344      {
342        G4cout << " -Checking node " << targetN    345        G4cout << " -Checking node " << targetNodeNo
343               << " is proxy with address " <<     346               << " is proxy with address " << sampleProxy << G4endl;
344      }                                            347      }
345 #endif                                            348 #endif 
346                                                   349 
347      if ( sampleProxy == nullptr )             << 350      if ( sampleProxy == 0 )
348      {                                            351      {
349        G4ExceptionDescription ed;                 352        G4ExceptionDescription ed;
350        ed << " Problem for node number= " << t    353        ed << " Problem for node number= " << targetNodeNo
351           << "    Number of slides = " << targ    354           << "    Number of slides = " << targetHeaderNoSlices
352           << G4endl;                              355           << G4endl;
353        G4Exception( "G4VoxelSafety::SafetyForV    356        G4Exception( "G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav0003",
354                     FatalException, ed,           357                     FatalException, ed,
355                     "Problem sampleProxy is Ze    358                     "Problem sampleProxy is Zero. Failure in loop.");
356      }                                            359      }
357      else if ( sampleProxy->IsNode() )            360      else if ( sampleProxy->IsNode() )
358      {                                            361      {
359        targetVoxelNode = sampleProxy->GetNode(    362        targetVoxelNode = sampleProxy->GetNode();
360                                                   363 
361        // Deal with the node here [ i.e. the l    364        // Deal with the node here [ i.e. the last level ]
362        //                                         365        //
363        nodeSafety= SafetyForVoxelNode( targetV    366        nodeSafety= SafetyForVoxelNode( targetVoxelNode, localPoint);
364 #ifdef G4DEBUG_NAVIGATION                         367 #ifdef G4DEBUG_NAVIGATION
365        if( fVerbose > 2 )                         368        if( fVerbose > 2 )
366        {                                          369        {
367           G4cout << " -- It is a Node ";          370           G4cout << " -- It is a Node ";
368           G4cout << "    its safety= " << node    371           G4cout << "    its safety= " << nodeSafety
369                  << " our level Saf = " << our    372                  << " our level Saf = " << ourSafety
370                  << "  MinSaf= " << minSafety     373                  << "  MinSaf= " << minSafety << G4endl;
371        }                                          374        }
372 #endif                                            375 #endif
373         ourSafety= std::min( ourSafety, nodeSa    376         ourSafety= std::min( ourSafety, nodeSafety );
374                                                   377        
375         trialUp = targetVoxelNode->GetMaxEquiv    378         trialUp = targetVoxelNode->GetMaxEquivalentSliceNo()+1;
376         trialDown = targetVoxelNode->GetMinEqu    379         trialDown = targetVoxelNode->GetMinEquivalentSliceNo()-1;
377      }                                            380      }
378      else                                         381      else  
379      {                                            382      {
380         const G4SmartVoxelHeader* pNewVoxelHea    383         const G4SmartVoxelHeader* pNewVoxelHeader = sampleProxy->GetHeader();
381                                                   384 
382         G4double distCombined_Sq;                 385         G4double distCombined_Sq;
383         distCombined_Sq = distUpperDepth_Sq +     386         distCombined_Sq = distUpperDepth_Sq + distAxis*distAxis;
384                                                   387 
385 #ifdef G4DEBUG_NAVIGATION                         388 #ifdef G4DEBUG_NAVIGATION
386         if( fVerbose > 2 )                        389         if( fVerbose > 2 )
387         {                                         390         {
388           G4double distCombined= std::sqrt( di    391           G4double distCombined= std::sqrt( distCombined_Sq );
389           G4double distUpperDepth= std::sqrt (    392           G4double distUpperDepth= std::sqrt ( distUpperDepth_Sq );
390           G4cout << " -- It is a Header " << G    393           G4cout << " -- It is a Header " << G4endl;
391           G4cout << "  Recurse to deal with ne    394           G4cout << "  Recurse to deal with next level, fVoxelDepth= " 
392                  << fVoxelDepth+1 << G4endl;      395                  << fVoxelDepth+1 << G4endl;
393           G4cout << "  Distances:  Upper= " <<    396           G4cout << "  Distances:  Upper= " << distUpperDepth
394                  << " Axis= " << distAxis         397                  << " Axis= " << distAxis
395                  << " Combined= " << distCombi    398                  << " Combined= " << distCombined << G4endl;
396         }                                         399         }
397 #endif                                            400 #endif
398                                                   401         
399         // Recurse to deal with lower levels      402         // Recurse to deal with lower levels
400         //                                        403         //
401         headerSafety= SafetyForVoxelHeader( pN    404         headerSafety= SafetyForVoxelHeader( pNewVoxelHeader, localPoint,
402                                             ma    405                                             maxLength, currentPhysical,
403                                             di    406                                             distCombined_Sq, minSafety);
404         ourSafety = std::min( ourSafety, heade    407         ourSafety = std::min( ourSafety, headerSafety );
405                                                   408        
406 #ifdef G4DEBUG_NAVIGATION                         409 #ifdef G4DEBUG_NAVIGATION
407         if( fVerbose > 2 )                        410         if( fVerbose > 2 )
408         {                                         411         { 
409           G4cout << "      >> Header safety =     412           G4cout << "      >> Header safety = " << headerSafety
410                  << " our level Saf = " << our    413                  << " our level Saf = " << ourSafety << G4endl;
411         }                                         414         }
412 #endif                                            415 #endif
413         trialUp =   pNewVoxelHeader->GetMaxEqu    416         trialUp =   pNewVoxelHeader->GetMaxEquivalentSliceNo()+1;
414         trialDown = pNewVoxelHeader->GetMinEqu    417         trialDown = pNewVoxelHeader->GetMinEquivalentSliceNo()-1;
415      }                                            418      }
416      minSafety = std::min( minSafety, ourSafet    419      minSafety = std::min( minSafety, ourSafety );
417                                                   420 
418      // Find next closest Voxel                   421      // Find next closest Voxel
419      //    - first try: by simple subtraction     422      //    - first try: by simple subtraction
420      //    - later:  using distance  (TODO - t    423      //    - later:  using distance  (TODO - tbc)
421      //                                           424      //
422      if( targetNodeNo >= pointNodeNo )            425      if( targetNodeNo >= pointNodeNo )
423      {                                            426      {
424         nextUp = trialUp;                         427         nextUp = trialUp;
425         // distUp   = std::max( targetHeaderMa    428         // distUp   = std::max( targetHeaderMax-localCrd,  0.0 );
426         G4double lowerEdgeOfNext = targetHeade    429         G4double lowerEdgeOfNext = targetHeaderMin
427                                  + nextUp * ta    430                                  + nextUp * targetHeaderNodeWidth;
428         distUp = lowerEdgeOfNext-localCrd ;       431         distUp = lowerEdgeOfNext-localCrd ;
429         if( distUp < 0.0 )                        432         if( distUp < 0.0 )
430         {                                         433         {
431            distUp = DBL_MAX;  // On the wrong     434            distUp = DBL_MAX;  // On the wrong side - must not be considered
432         }                                         435         }
433 #ifdef G4DEBUG_NAVIGATION                         436 #ifdef G4DEBUG_NAVIGATION
434        if( fVerbose > 2 )                         437        if( fVerbose > 2 )
435        {                                          438        {
436          G4cout << " > Updated nextUp= " << ne    439          G4cout << " > Updated nextUp= " << nextUp << G4endl;
437        }                                          440        }
438 #endif                                            441 #endif
439      }                                            442      }
440                                                   443      
441      if( targetNodeNo <= pointNodeNo )            444      if( targetNodeNo <= pointNodeNo ) 
442      {                                            445      {
443          nextDown = trialDown;                    446          nextDown = trialDown;
444          // distDown = std::max( localCrd-targ    447          // distDown = std::max( localCrd-targetHeaderMin,  0.0);
445          G4double upperEdgeOfNext = targetHead    448          G4double upperEdgeOfNext = targetHeaderMin
446                                   + (nextDown+    449                                   + (nextDown+1) * targetHeaderNodeWidth;
447          distDown = localCrd-upperEdgeOfNext;     450          distDown = localCrd-upperEdgeOfNext;
448          if( distDown < 0.0 )                     451          if( distDown < 0.0 )
449          {                                        452          {
450             distDown= DBL_MAX; // On the wrong    453             distDown= DBL_MAX; // On the wrong side - must not be considered 
451          }                                        454          }
452 #ifdef G4DEBUG_NAVIGATION                         455 #ifdef G4DEBUG_NAVIGATION
453        if( fVerbose > 2 )                         456        if( fVerbose > 2 )
454        {                                          457        {
455          G4cout << " > Updated nextDown= " <<     458          G4cout << " > Updated nextDown= " << nextDown << G4endl;
456        }                                          459        }
457 #endif                                            460 #endif
458      }                                            461      }
459                                                   462 
460 #ifdef G4DEBUG_NAVIGATION                         463 #ifdef G4DEBUG_NAVIGATION
461      if( fVerbose > 2 )                           464      if( fVerbose > 2 )
462      {                                            465      {
463        G4cout << " Node= " << pointNodeNo         466        G4cout << " Node= " << pointNodeNo
464               << " Up:   next= " << nextUp  <<    467               << " Up:   next= " << nextUp  << " d# "
465               << nextUp - pointNodeNo             468               << nextUp - pointNodeNo
466               << " trialUp=  " << trialUp << "    469               << " trialUp=  " << trialUp << " d# "
467               << trialUp - pointNodeNo            470               << trialUp - pointNodeNo
468               << " Down: next= " << nextDown <    471               << " Down: next= " << nextDown << " d# "
469               << targetNodeNo - nextDown          472               << targetNodeNo - nextDown
470               << " trialDwn= " << trialDown <<    473               << " trialDwn= " << trialDown << " d# "
471               << targetNodeNo - trialDown         474               << targetNodeNo - trialDown
472               << " condition (next is Inside)=    475               << " condition (next is Inside)= " << nextIsInside
473               << G4endl;                          476               << G4endl;
474      }                                            477      }
475 #endif                                            478 #endif     
476                                                   479      
477      G4bool UpIsClosest;                          480      G4bool UpIsClosest;
478      UpIsClosest = distUp < distDown;             481      UpIsClosest = distUp < distDown;
479                                                   482      
480      if( (nextUp < targetHeaderNoSlices)          483      if( (nextUp < targetHeaderNoSlices)
481          && (UpIsClosest || (nextDown < 0)) )     484          && (UpIsClosest || (nextDown < 0)) )
482      {                                            485      {
483        nextNodeNo = nextUp;                       486        nextNodeNo = nextUp;
484        distAxis = distUp;                         487        distAxis = distUp;
485        ++nextUp; // Default                       488        ++nextUp; // Default
486 #ifdef G4VERBOSE                                  489 #ifdef G4VERBOSE
487        if( fVerbose > 2 )                         490        if( fVerbose > 2 )
488        {                                          491        {
489           G4cout << " > Chose Up.   Depth= " <    492           G4cout << " > Chose Up.   Depth= " << fVoxelDepth
490                  << "   Nodes: next= " << next    493                  << "   Nodes: next= " << nextNodeNo
491                  << " new nextUp=   " << nextU    494                  << " new nextUp=   " << nextUp
492                  << " Dist = " << distAxis <<     495                  << " Dist = " << distAxis << G4endl;
493        }                                          496        }
494 #endif                                            497 #endif
495      }                                            498      }
496      else                                         499      else
497      {                                            500      {
498        nextNodeNo = nextDown;                     501        nextNodeNo = nextDown;
499        distAxis = distDown;                       502        distAxis = distDown;
500        --nextDown; // A default value             503        --nextDown; // A default value
501 #ifdef G4VERBOSE                                  504 #ifdef G4VERBOSE
502        if( fVerbose > 2 )                         505        if( fVerbose > 2 )
503        {                                          506        {
504          G4cout << " > Chose Down.  Depth= " <    507          G4cout << " > Chose Down.  Depth= " << fVoxelDepth
505                 << "  Nodes: next= " << nextNo    508                 << "  Nodes: next= " << nextNodeNo
506                 << " new nextDown=  " << nextD    509                 << " new nextDown=  " << nextDown
507                 << " Dist = " << distAxis << G    510                 << " Dist = " << distAxis << G4endl;
508        }                                          511        }
509 #endif                                            512 #endif
510      }                                            513      }
511                                                   514 
512      nextIsInside = (nextNodeNo >= 0) && (next    515      nextIsInside = (nextNodeNo >= 0) && (nextNodeNo < targetHeaderNoSlices);
513      if( nextIsInside )                           516      if( nextIsInside )
514      {                                            517      {
515        targetNodeNo= nextNodeNo;                  518        targetNodeNo= nextNodeNo;
516                                                   519 
517 #ifdef G4DEBUG_NAVIGATION                         520 #ifdef G4DEBUG_NAVIGATION
518        assert( targetVoxelHeader->GetSlice(nex    521        assert( targetVoxelHeader->GetSlice(nextNodeNo) != 0 );
519        G4bool bContinue = (distAxis<minSafety)    522        G4bool bContinue = (distAxis<minSafety);
520        if( !bContinue )                           523        if( !bContinue )
521        {                                          524        {
522          if( fVerbose > 2 )                       525          if( fVerbose > 2 )
523          {                                        526          {
524            G4cout << " Can skip remaining at d    527            G4cout << " Can skip remaining at depth " << targetHeaderAxis
525                   << " >>  distAxis= " << dist    528                   << " >>  distAxis= " << distAxis
526                   << " minSafety= " << minSafe    529                   << " minSafety= " << minSafety << G4endl;
527          }                                        530          }
528        }                                          531        }
529 #endif                                            532 #endif
530      }                                            533      }
531      else                                         534      else
532      {                                            535      {
533 #ifdef G4DEBUG_NAVIGATION                         536 #ifdef G4DEBUG_NAVIGATION
534        if( fVerbose > 2)                          537        if( fVerbose > 2)
535        {                                          538        {
536          G4cout << " VoxSaf> depth= " << fVoxe    539          G4cout << " VoxSaf> depth= " << fVoxelDepth << G4endl;
537          G4cout << " VoxSaf> No more candidate    540          G4cout << " VoxSaf> No more candidates:  nodeDown= " << nextDown
538                 << " nodeUp= " << nextUp          541                 << " nodeUp= " << nextUp
539                 << " lastSlice= " << targetHea    542                 << " lastSlice= " << targetHeaderNoSlices << G4endl;
540        }                                          543        }
541 #endif                                            544 #endif
542      }                                            545      }
543                                                   546 
544      // This calculation can be 'hauled'-up to    547      // This calculation can be 'hauled'-up to where minSafety is calculated
545      //                                           548      //
546      distMaxInterest = std::min( minSafety, ma    549      distMaxInterest = std::min( minSafety, maxLength ); 
547                                                   550 
548   } while ( nextIsInside && ( distAxis*distAxi    551   } while ( nextIsInside && ( distAxis*distAxis + distUpperDepth_Sq
549                             < distMaxInterest*    552                             < distMaxInterest*distMaxInterest ) ); 
550                                                   553 
551 #ifdef G4VERBOSE                                  554 #ifdef G4VERBOSE
552   if( fVerbose > 0 )                              555   if( fVerbose > 0 )
553   {                                               556   { 
554     G4cout << " Ended for targetNodeNo -- chec    557     G4cout << " Ended for targetNodeNo -- checked " << checkedNum << " out of "
555            << targetHeaderNoSlices << " slices    558            << targetHeaderNoSlices << " slices." << G4endl;
556     G4cout << " ===== Returning from SafetyFor    559     G4cout << " ===== Returning from SafetyForVoxelHeader " 
557            << "  Depth= " << fVoxelDepth << G4    560            << "  Depth= " << fVoxelDepth << G4endl
558            << G4endl;                             561            << G4endl;  
559   }                                               562   }
560 #endif                                            563 #endif
561                                                   564 
562   // Go back one level                            565   // Go back one level
563   fVoxelDepth--;                                  566   fVoxelDepth--; 
564                                                   567   
565   return ourSafety;                               568   return ourSafety;
566 }                                                 569 }
567                                                   570