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 10.7.p1)


  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 = 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    G4int curNoVolumes, contentNo, sampleNo;
161    G4int sampleNo;                             << 
162    G4VPhysicalVolume* samplePhysical;             164    G4VPhysicalVolume* samplePhysical;
163                                                   165 
164    G4double      sampleSafety = 0.0;              166    G4double      sampleSafety = 0.0; 
165    G4ThreeVector samplePoint;                     167    G4ThreeVector samplePoint;
166    G4VSolid*     ptrSolid = nullptr;              168    G4VSolid*     ptrSolid = nullptr;
167                                                   169 
168    curNoVolumes = curVoxelNode->GetNoContained    170    curNoVolumes = curVoxelNode->GetNoContained();
169                                                   171 
170    for ( contentNo=curNoVolumes-1; contentNo>=    172    for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
171    {                                              173    {
172       sampleNo = curVoxelNode->GetVolume((G4in << 174       sampleNo = curVoxelNode->GetVolume(contentNo);
173       if ( !fBlockList.IsBlocked(sampleNo) )      175       if ( !fBlockList.IsBlocked(sampleNo) ) 
174       {                                           176       { 
175         fBlockList.BlockVolume(sampleNo);         177         fBlockList.BlockVolume(sampleNo);
176                                                   178 
177         samplePhysical = fpMotherLogical->GetD    179         samplePhysical = fpMotherLogical->GetDaughter(sampleNo);
178         G4AffineTransform sampleTf(samplePhysi    180         G4AffineTransform sampleTf(samplePhysical->GetRotation(),
179                                    samplePhysi    181                                    samplePhysical->GetTranslation());
180         sampleTf.Invert();                        182         sampleTf.Invert();
181         samplePoint = sampleTf.TransformPoint(    183         samplePoint = sampleTf.TransformPoint(localPoint);
182         ptrSolid = samplePhysical->GetLogicalV    184         ptrSolid = samplePhysical->GetLogicalVolume()->GetSolid();
183                                                   185 
184         sampleSafety = ptrSolid->DistanceToIn(    186         sampleSafety = ptrSolid->DistanceToIn(samplePoint);
185         ourSafety = std::min( sampleSafety, ou    187         ourSafety = std::min( sampleSafety, ourSafety ); 
186 #ifdef G4VERBOSE                                  188 #ifdef G4VERBOSE
187         if(( fCheck ) && ( fVerbose == 1 ))       189         if(( fCheck ) && ( fVerbose == 1 ))
188         {                                         190         {
189           G4cout << "*** G4VoxelSafety::Safety    191           G4cout << "*** G4VoxelSafety::SafetyForVoxelNode(): ***" << G4endl
190                  << "    Invoked DistanceToIn(    192                  << "    Invoked DistanceToIn(p) for daughter solid: "
191                  << ptrSolid->GetName()           193                  << ptrSolid->GetName()
192                  << ". Solid replied: " << sam    194                  << ". Solid replied: " << sampleSafety << G4endl
193                  << "    For local point p: "     195                  << "    For local point p: " << samplePoint
194                  << ", to be considered as 'da    196                  << ", to be considered as 'daughter safety'." << G4endl;
195         }                                         197         }
196 #endif                                            198 #endif
197       }                                           199       }
198    }  // end for contentNo                        200    }  // end for contentNo
199                                                   201 
200    return ourSafety;                              202    return ourSafety; 
201 }                                                 203 }
202                                                   204 
203 // *******************************************    205 // ********************************************************************
204 // SafetyForVoxelHeader                           206 // SafetyForVoxelHeader
205 //                                                207 //
206 // Cycles through levels of headers to process    208 // Cycles through levels of headers to process each node level
207 // Obtained by modifying VoxelLocate (to cycle    209 // Obtained by modifying VoxelLocate (to cycle through Node Headers)
208 // *******************************************    210 // *********************************************************************
209 //                                                211 //
210 G4double                                          212 G4double
211 G4VoxelSafety::SafetyForVoxelHeader( const G4S    213 G4VoxelSafety::SafetyForVoxelHeader( const G4SmartVoxelHeader* pHeader,
212                                      const G4T    214                                      const G4ThreeVector& localPoint,
213                                            G4d    215                                            G4double maxLength,
214                                      const G4V    216                                      const G4VPhysicalVolume&  currentPhysical, //Debug
215                                            G4d    217                                            G4double distUpperDepth_Sq,
216                                            G4d    218                                            G4double previousMinSafety
217                                    )              219                                    )
218 {                                                 220 {
219   const G4SmartVoxelHeader* const targetVoxelH    221   const G4SmartVoxelHeader* const targetVoxelHeader = pHeader;
220   G4SmartVoxelNode* targetVoxelNode = nullptr;    222   G4SmartVoxelNode* targetVoxelNode = nullptr;
221                                                   223 
222   const G4SmartVoxelProxy* sampleProxy;           224   const G4SmartVoxelProxy* sampleProxy;
223   EAxis    targetHeaderAxis;                      225   EAxis    targetHeaderAxis;
224   G4double targetHeaderMin, targetHeaderMax, t    226   G4double targetHeaderMin, targetHeaderMax, targetHeaderNodeWidth;
225   G4int    targetHeaderNoSlices;                  227   G4int    targetHeaderNoSlices;
226   G4int    targetNodeNo;                          228   G4int    targetNodeNo;
227                                                   229 
228   G4double minSafety = previousMinSafety;         230   G4double minSafety = previousMinSafety;
229   G4double ourSafety = DBL_MAX;                   231   G4double ourSafety = DBL_MAX;
230   unsigned int checkedNum= 0;                     232   unsigned int checkedNum= 0;
231                                                   233   
232   ++fVoxelDepth;                                  234   ++fVoxelDepth;
233   // fVoxelDepth  set by ComputeSafety or prev    235   // fVoxelDepth  set by ComputeSafety or previous level call
234                                                   236 
235   targetHeaderAxis =      targetVoxelHeader->G    237   targetHeaderAxis =      targetVoxelHeader->GetAxis();
236   targetHeaderNoSlices =  (G4int)targetVoxelHe << 238   targetHeaderNoSlices =  targetVoxelHeader->GetNoSlices();
237   targetHeaderMin =       targetVoxelHeader->G    239   targetHeaderMin =       targetVoxelHeader->GetMinExtent();
238   targetHeaderMax =       targetVoxelHeader->G    240   targetHeaderMax =       targetVoxelHeader->GetMaxExtent();
239                                                   241 
240   targetHeaderNodeWidth = (targetHeaderMax-tar    242   targetHeaderNodeWidth = (targetHeaderMax-targetHeaderMin)
241                           / targetHeaderNoSlic    243                           / targetHeaderNoSlices;
242                                                   244 
243   G4double localCrd = localPoint(targetHeaderA    245   G4double localCrd = localPoint(targetHeaderAxis);
244                                                   246 
245   const auto  candNodeNo = G4int( (localCrd-ta << 247   const G4int candNodeNo = G4int( (localCrd-targetHeaderMin)
246                                  / targetHeade    248                                  / targetHeaderNodeWidth );
247   // Ensure that it is between 0 and targetHea    249   // Ensure that it is between 0 and targetHeader->GetMaxExtent() - 1
248                                                   250 
249 #ifdef G4DEBUG_VOXELISATION                       251 #ifdef G4DEBUG_VOXELISATION  
250   if( candNodeNo < 0 || candNodeNo > targetHea    252   if( candNodeNo < 0 || candNodeNo > targetHeaderNoSlices-1 )
251   {                                               253   {
252      G4ExceptionDescription ed;                   254      G4ExceptionDescription ed;
253      ed << " Potential ERROR."                    255      ed << " Potential ERROR."
254         << " Point is outside range of Voxel i    256         << " Point is outside range of Voxel in current coordinate" << G4endl;
255      ed << "  Node number of point " << localP    257      ed << "  Node number of point " << localPoint
256         << "is outside the range. " << G4endl;    258         << "is outside the range. " << G4endl;
257      ed << "    Voxel node Num= " << candNodeN    259      ed << "    Voxel node Num= " << candNodeNo << " versus  minimum= " << 0
258         << " and maximum= " << targetHeaderNoS    260         << " and maximum= " << targetHeaderNoSlices-1 << G4endl;
259      ed << "    Axis = " << targetHeaderAxis      261      ed << "    Axis = " << targetHeaderAxis
260         << "  No of slices = " << targetHeader    262         << "  No of slices = " << targetHeaderNoSlices << G4endl;
261      ed << "    Local coord = " << localCrd       263      ed << "    Local coord = " << localCrd
262         << "  Voxel Min = " << targetHeaderMin    264         << "  Voxel Min = " << targetHeaderMin
263         <<            " Max = " << targetHeade    265         <<            " Max = " << targetHeaderMax << G4endl;
264      G4LogicalVolume *pLogical= currentPhysica    266      G4LogicalVolume *pLogical= currentPhysical.GetLogicalVolume();
265      ed << "  Current volume (physical) = " <<    267      ed << "  Current volume (physical) = " << currentPhysical.GetName()
266         << "   (logical) = " << pLogical->GetN    268         << "   (logical) = " << pLogical->GetName()     << G4endl;
267      G4VSolid* pSolid= pLogical->GetSolid();      269      G4VSolid* pSolid= pLogical->GetSolid();
268      ed << "  Solid type = " << pSolid->GetEnt    270      ed << "  Solid type = " << pSolid->GetEntityType() << G4endl;
269      ed << *pSolid << G4endl;                     271      ed << *pSolid << G4endl;
270      G4Exception("G4VoxelSafety::SafetyForVoxe    272      G4Exception("G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav1003",
271                  JustWarning, ed,                 273                  JustWarning, ed,
272                  "Point is outside range of Vo    274                  "Point is outside range of Voxel in current coordinate");
273   }                                               275   }
274 #endif                                            276 #endif
275                                                   277   
276   const G4int pointNodeNo =                       278   const G4int pointNodeNo =
277               std::max( 0, std::min( candNodeN    279               std::max( 0, std::min( candNodeNo, targetHeaderNoSlices-1 ) );
278                                                   280   
279 #ifdef G4VERBOSE                                  281 #ifdef G4VERBOSE  
280   if( fVerbose > 2 )                              282   if( fVerbose > 2 )
281   {                                               283   { 
282     G4cout << G4endl;                             284     G4cout << G4endl;
283     G4cout << "**** G4VoxelSafety::SafetyForVo    285     G4cout << "**** G4VoxelSafety::SafetyForVoxelHeader  " << G4endl;
284     G4cout << "  Called at Depth     = " << fV    286     G4cout << "  Called at Depth     = " << fVoxelDepth;
285     G4cout << " Calculated pointNodeNo= " << p    287     G4cout << " Calculated pointNodeNo= " << pointNodeNo
286            << "  from position= " <<  localPoi    288            << "  from position= " <<  localPoint(targetHeaderAxis)
287            << "  min= "    << targetHeaderMin     289            << "  min= "    << targetHeaderMin
288            << "  max= "    << targetHeaderMax     290            << "  max= "    << targetHeaderMax
289            << "  width= "  << targetHeaderNode    291            << "  width= "  << targetHeaderNodeWidth 
290            << "  no-slices= " << targetHeaderN    292            << "  no-slices= " << targetHeaderNoSlices
291            << "  axis=  "  << targetHeaderAxis    293            << "  axis=  "  << targetHeaderAxis   << G4endl;
292   }                                               294   }
293   else if (fVerbose == 1)                         295   else if (fVerbose == 1)
294   {                                               296   {
295     G4cout << " VoxelSafety: Depth  = " << fVo    297     G4cout << " VoxelSafety: Depth  = " << fVoxelDepth 
296            << " Number of Slices = " << target    298            << " Number of Slices = " << targetHeaderNoSlices
297            << " Header (address) = " << target    299            << " Header (address) = " << targetVoxelHeader  << G4endl;
298   }                                               300   }
299 #endif                                            301 #endif
300                                                   302 
301   // Stack info for stepping                      303   // Stack info for stepping
302   //                                              304   //
303   fVoxelAxisStack[fVoxelDepth] = targetHeaderA    305   fVoxelAxisStack[fVoxelDepth] = targetHeaderAxis;
304   fVoxelNoSlicesStack[fVoxelDepth] = targetHea    306   fVoxelNoSlicesStack[fVoxelDepth] = targetHeaderNoSlices;
305   fVoxelSliceWidthStack[fVoxelDepth] = targetH    307   fVoxelSliceWidthStack[fVoxelDepth] = targetHeaderNodeWidth;
306                                                   308 
307   fVoxelHeaderStack[fVoxelDepth] = pHeader;       309   fVoxelHeaderStack[fVoxelDepth] = pHeader;
308                                                   310 
309   G4int   trialUp = -1,     trialDown = -1;       311   G4int   trialUp = -1,     trialDown = -1;
310   G4double distUp = DBL_MAX, distDown = DBL_MA    312   G4double distUp = DBL_MAX, distDown = DBL_MAX;
311                                                   313 
312   // Using Equivalent voxels - this is pre-ini    314   // Using Equivalent voxels - this is pre-initialisation only
313   //                                              315   //
314   G4int nextUp =   pointNodeNo+1;                 316   G4int nextUp =   pointNodeNo+1;
315   G4int nextDown = pointNodeNo-1;                 317   G4int nextDown = pointNodeNo-1;
316                                                   318 
317   G4int    nextNodeNo = pointNodeNo;              319   G4int    nextNodeNo = pointNodeNo;
318   G4double distAxis;  // Distance in current A    320   G4double distAxis;  // Distance in current Axis
319   distAxis = 0.0;  // Starting in node contain    321   distAxis = 0.0;  // Starting in node containing local Coordinate
320                                                   322 
321   G4bool nextIsInside = false;                    323   G4bool nextIsInside = false;
322                                                   324 
323   G4double distMaxInterest= std::min( previous    325   G4double distMaxInterest= std::min( previousMinSafety, maxLength);
324     // We will not look beyond this distance.     326     // We will not look beyond this distance.
325     // This distance will be updated to reflec    327     // This distance will be updated to reflect the
326     // max ( minSafety, maxLength ) at each st    328     // max ( minSafety, maxLength ) at each step
327                                                   329 
328   targetNodeNo = pointNodeNo;                     330   targetNodeNo = pointNodeNo;
329   do                                              331   do
330   {                                               332   {
331      G4double nodeSafety = DBL_MAX, headerSafe    333      G4double nodeSafety = DBL_MAX, headerSafety = DBL_MAX;
332      fVoxelNodeNoStack[fVoxelDepth] = targetNo    334      fVoxelNodeNoStack[fVoxelDepth] = targetNodeNo;
333                                                   335     
334      ++checkedNum;                                336      ++checkedNum; 
335                                                   337     
336      sampleProxy = targetVoxelHeader->GetSlice    338      sampleProxy = targetVoxelHeader->GetSlice(targetNodeNo);
337                                                   339 
338 #ifdef G4DEBUG_NAVIGATION                         340 #ifdef G4DEBUG_NAVIGATION
339      assert( sampleProxy != 0);                   341      assert( sampleProxy != 0);
340      if( fVerbose > 2 )                           342      if( fVerbose > 2 ) 
341      {                                            343      {
342        G4cout << " -Checking node " << targetN    344        G4cout << " -Checking node " << targetNodeNo
343               << " is proxy with address " <<     345               << " is proxy with address " << sampleProxy << G4endl;
344      }                                            346      }
345 #endif                                            347 #endif 
346                                                   348 
347      if ( sampleProxy == nullptr )             << 349      if ( sampleProxy == 0 )
348      {                                            350      {
349        G4ExceptionDescription ed;                 351        G4ExceptionDescription ed;
350        ed << " Problem for node number= " << t    352        ed << " Problem for node number= " << targetNodeNo
351           << "    Number of slides = " << targ    353           << "    Number of slides = " << targetHeaderNoSlices
352           << G4endl;                              354           << G4endl;
353        G4Exception( "G4VoxelSafety::SafetyForV    355        G4Exception( "G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav0003",
354                     FatalException, ed,           356                     FatalException, ed,
355                     "Problem sampleProxy is Ze    357                     "Problem sampleProxy is Zero. Failure in loop.");
356      }                                            358      }
357      else if ( sampleProxy->IsNode() )            359      else if ( sampleProxy->IsNode() )
358      {                                            360      {
359        targetVoxelNode = sampleProxy->GetNode(    361        targetVoxelNode = sampleProxy->GetNode();
360                                                   362 
361        // Deal with the node here [ i.e. the l    363        // Deal with the node here [ i.e. the last level ]
362        //                                         364        //
363        nodeSafety= SafetyForVoxelNode( targetV    365        nodeSafety= SafetyForVoxelNode( targetVoxelNode, localPoint);
364 #ifdef G4DEBUG_NAVIGATION                         366 #ifdef G4DEBUG_NAVIGATION
365        if( fVerbose > 2 )                         367        if( fVerbose > 2 )
366        {                                          368        {
367           G4cout << " -- It is a Node ";          369           G4cout << " -- It is a Node ";
368           G4cout << "    its safety= " << node    370           G4cout << "    its safety= " << nodeSafety
369                  << " our level Saf = " << our    371                  << " our level Saf = " << ourSafety
370                  << "  MinSaf= " << minSafety     372                  << "  MinSaf= " << minSafety << G4endl;
371        }                                          373        }
372 #endif                                            374 #endif
373         ourSafety= std::min( ourSafety, nodeSa    375         ourSafety= std::min( ourSafety, nodeSafety );
374                                                   376        
375         trialUp = targetVoxelNode->GetMaxEquiv    377         trialUp = targetVoxelNode->GetMaxEquivalentSliceNo()+1;
376         trialDown = targetVoxelNode->GetMinEqu    378         trialDown = targetVoxelNode->GetMinEquivalentSliceNo()-1;
377      }                                            379      }
378      else                                         380      else  
379      {                                            381      {
380         const G4SmartVoxelHeader* pNewVoxelHea    382         const G4SmartVoxelHeader* pNewVoxelHeader = sampleProxy->GetHeader();
381                                                   383 
382         G4double distCombined_Sq;                 384         G4double distCombined_Sq;
383         distCombined_Sq = distUpperDepth_Sq +     385         distCombined_Sq = distUpperDepth_Sq + distAxis*distAxis;
384                                                   386 
385 #ifdef G4DEBUG_NAVIGATION                         387 #ifdef G4DEBUG_NAVIGATION
386         if( fVerbose > 2 )                        388         if( fVerbose > 2 )
387         {                                         389         {
388           G4double distCombined= std::sqrt( di    390           G4double distCombined= std::sqrt( distCombined_Sq );
389           G4double distUpperDepth= std::sqrt (    391           G4double distUpperDepth= std::sqrt ( distUpperDepth_Sq );
390           G4cout << " -- It is a Header " << G    392           G4cout << " -- It is a Header " << G4endl;
391           G4cout << "  Recurse to deal with ne    393           G4cout << "  Recurse to deal with next level, fVoxelDepth= " 
392                  << fVoxelDepth+1 << G4endl;      394                  << fVoxelDepth+1 << G4endl;
393           G4cout << "  Distances:  Upper= " <<    395           G4cout << "  Distances:  Upper= " << distUpperDepth
394                  << " Axis= " << distAxis         396                  << " Axis= " << distAxis
395                  << " Combined= " << distCombi    397                  << " Combined= " << distCombined << G4endl;
396         }                                         398         }
397 #endif                                            399 #endif
398                                                   400         
399         // Recurse to deal with lower levels      401         // Recurse to deal with lower levels
400         //                                        402         //
401         headerSafety= SafetyForVoxelHeader( pN    403         headerSafety= SafetyForVoxelHeader( pNewVoxelHeader, localPoint,
402                                             ma    404                                             maxLength, currentPhysical,
403                                             di    405                                             distCombined_Sq, minSafety);
404         ourSafety = std::min( ourSafety, heade    406         ourSafety = std::min( ourSafety, headerSafety );
405                                                   407        
406 #ifdef G4DEBUG_NAVIGATION                         408 #ifdef G4DEBUG_NAVIGATION
407         if( fVerbose > 2 )                        409         if( fVerbose > 2 )
408         {                                         410         { 
409           G4cout << "      >> Header safety =     411           G4cout << "      >> Header safety = " << headerSafety
410                  << " our level Saf = " << our    412                  << " our level Saf = " << ourSafety << G4endl;
411         }                                         413         }
412 #endif                                            414 #endif
413         trialUp =   pNewVoxelHeader->GetMaxEqu    415         trialUp =   pNewVoxelHeader->GetMaxEquivalentSliceNo()+1;
414         trialDown = pNewVoxelHeader->GetMinEqu    416         trialDown = pNewVoxelHeader->GetMinEquivalentSliceNo()-1;
415      }                                            417      }
416      minSafety = std::min( minSafety, ourSafet    418      minSafety = std::min( minSafety, ourSafety );
417                                                   419 
418      // Find next closest Voxel                   420      // Find next closest Voxel
419      //    - first try: by simple subtraction     421      //    - first try: by simple subtraction
420      //    - later:  using distance  (TODO - t    422      //    - later:  using distance  (TODO - tbc)
421      //                                           423      //
422      if( targetNodeNo >= pointNodeNo )            424      if( targetNodeNo >= pointNodeNo )
423      {                                            425      {
424         nextUp = trialUp;                         426         nextUp = trialUp;
425         // distUp   = std::max( targetHeaderMa    427         // distUp   = std::max( targetHeaderMax-localCrd,  0.0 );
426         G4double lowerEdgeOfNext = targetHeade    428         G4double lowerEdgeOfNext = targetHeaderMin
427                                  + nextUp * ta    429                                  + nextUp * targetHeaderNodeWidth;
428         distUp = lowerEdgeOfNext-localCrd ;       430         distUp = lowerEdgeOfNext-localCrd ;
429         if( distUp < 0.0 )                        431         if( distUp < 0.0 )
430         {                                         432         {
431            distUp = DBL_MAX;  // On the wrong     433            distUp = DBL_MAX;  // On the wrong side - must not be considered
432         }                                         434         }
433 #ifdef G4DEBUG_NAVIGATION                         435 #ifdef G4DEBUG_NAVIGATION
434        if( fVerbose > 2 )                         436        if( fVerbose > 2 )
435        {                                          437        {
436          G4cout << " > Updated nextUp= " << ne    438          G4cout << " > Updated nextUp= " << nextUp << G4endl;
437        }                                          439        }
438 #endif                                            440 #endif
439      }                                            441      }
440                                                   442      
441      if( targetNodeNo <= pointNodeNo )            443      if( targetNodeNo <= pointNodeNo ) 
442      {                                            444      {
443          nextDown = trialDown;                    445          nextDown = trialDown;
444          // distDown = std::max( localCrd-targ    446          // distDown = std::max( localCrd-targetHeaderMin,  0.0);
445          G4double upperEdgeOfNext = targetHead    447          G4double upperEdgeOfNext = targetHeaderMin
446                                   + (nextDown+    448                                   + (nextDown+1) * targetHeaderNodeWidth;
447          distDown = localCrd-upperEdgeOfNext;     449          distDown = localCrd-upperEdgeOfNext;
448          if( distDown < 0.0 )                     450          if( distDown < 0.0 )
449          {                                        451          {
450             distDown= DBL_MAX; // On the wrong    452             distDown= DBL_MAX; // On the wrong side - must not be considered 
451          }                                        453          }
452 #ifdef G4DEBUG_NAVIGATION                         454 #ifdef G4DEBUG_NAVIGATION
453        if( fVerbose > 2 )                         455        if( fVerbose > 2 )
454        {                                          456        {
455          G4cout << " > Updated nextDown= " <<     457          G4cout << " > Updated nextDown= " << nextDown << G4endl;
456        }                                          458        }
457 #endif                                            459 #endif
458      }                                            460      }
459                                                   461 
460 #ifdef G4DEBUG_NAVIGATION                         462 #ifdef G4DEBUG_NAVIGATION
461      if( fVerbose > 2 )                           463      if( fVerbose > 2 )
462      {                                            464      {
463        G4cout << " Node= " << pointNodeNo         465        G4cout << " Node= " << pointNodeNo
464               << " Up:   next= " << nextUp  <<    466               << " Up:   next= " << nextUp  << " d# "
465               << nextUp - pointNodeNo             467               << nextUp - pointNodeNo
466               << " trialUp=  " << trialUp << "    468               << " trialUp=  " << trialUp << " d# "
467               << trialUp - pointNodeNo            469               << trialUp - pointNodeNo
468               << " Down: next= " << nextDown <    470               << " Down: next= " << nextDown << " d# "
469               << targetNodeNo - nextDown          471               << targetNodeNo - nextDown
470               << " trialDwn= " << trialDown <<    472               << " trialDwn= " << trialDown << " d# "
471               << targetNodeNo - trialDown         473               << targetNodeNo - trialDown
472               << " condition (next is Inside)=    474               << " condition (next is Inside)= " << nextIsInside
473               << G4endl;                          475               << G4endl;
474      }                                            476      }
475 #endif                                            477 #endif     
476                                                   478      
477      G4bool UpIsClosest;                          479      G4bool UpIsClosest;
478      UpIsClosest = distUp < distDown;             480      UpIsClosest = distUp < distDown;
479                                                   481      
480      if( (nextUp < targetHeaderNoSlices)          482      if( (nextUp < targetHeaderNoSlices)
481          && (UpIsClosest || (nextDown < 0)) )     483          && (UpIsClosest || (nextDown < 0)) )
482      {                                            484      {
483        nextNodeNo = nextUp;                       485        nextNodeNo = nextUp;
484        distAxis = distUp;                         486        distAxis = distUp;
485        ++nextUp; // Default                       487        ++nextUp; // Default
486 #ifdef G4VERBOSE                                  488 #ifdef G4VERBOSE
487        if( fVerbose > 2 )                         489        if( fVerbose > 2 )
488        {                                          490        {
489           G4cout << " > Chose Up.   Depth= " <    491           G4cout << " > Chose Up.   Depth= " << fVoxelDepth
490                  << "   Nodes: next= " << next    492                  << "   Nodes: next= " << nextNodeNo
491                  << " new nextUp=   " << nextU    493                  << " new nextUp=   " << nextUp
492                  << " Dist = " << distAxis <<     494                  << " Dist = " << distAxis << G4endl;
493        }                                          495        }
494 #endif                                            496 #endif
495      }                                            497      }
496      else                                         498      else
497      {                                            499      {
498        nextNodeNo = nextDown;                     500        nextNodeNo = nextDown;
499        distAxis = distDown;                       501        distAxis = distDown;
500        --nextDown; // A default value             502        --nextDown; // A default value
501 #ifdef G4VERBOSE                                  503 #ifdef G4VERBOSE
502        if( fVerbose > 2 )                         504        if( fVerbose > 2 )
503        {                                          505        {
504          G4cout << " > Chose Down.  Depth= " <    506          G4cout << " > Chose Down.  Depth= " << fVoxelDepth
505                 << "  Nodes: next= " << nextNo    507                 << "  Nodes: next= " << nextNodeNo
506                 << " new nextDown=  " << nextD    508                 << " new nextDown=  " << nextDown
507                 << " Dist = " << distAxis << G    509                 << " Dist = " << distAxis << G4endl;
508        }                                          510        }
509 #endif                                            511 #endif
510      }                                            512      }
511                                                   513 
512      nextIsInside = (nextNodeNo >= 0) && (next    514      nextIsInside = (nextNodeNo >= 0) && (nextNodeNo < targetHeaderNoSlices);
513      if( nextIsInside )                           515      if( nextIsInside )
514      {                                            516      {
515        targetNodeNo= nextNodeNo;                  517        targetNodeNo= nextNodeNo;
516                                                   518 
517 #ifdef G4DEBUG_NAVIGATION                         519 #ifdef G4DEBUG_NAVIGATION
518        assert( targetVoxelHeader->GetSlice(nex    520        assert( targetVoxelHeader->GetSlice(nextNodeNo) != 0 );
519        G4bool bContinue = (distAxis<minSafety)    521        G4bool bContinue = (distAxis<minSafety);
520        if( !bContinue )                           522        if( !bContinue )
521        {                                          523        {
522          if( fVerbose > 2 )                       524          if( fVerbose > 2 )
523          {                                        525          {
524            G4cout << " Can skip remaining at d    526            G4cout << " Can skip remaining at depth " << targetHeaderAxis
525                   << " >>  distAxis= " << dist    527                   << " >>  distAxis= " << distAxis
526                   << " minSafety= " << minSafe    528                   << " minSafety= " << minSafety << G4endl;
527          }                                        529          }
528        }                                          530        }
529 #endif                                            531 #endif
530      }                                            532      }
531      else                                         533      else
532      {                                            534      {
533 #ifdef G4DEBUG_NAVIGATION                         535 #ifdef G4DEBUG_NAVIGATION
534        if( fVerbose > 2)                          536        if( fVerbose > 2)
535        {                                          537        {
536          G4cout << " VoxSaf> depth= " << fVoxe    538          G4cout << " VoxSaf> depth= " << fVoxelDepth << G4endl;
537          G4cout << " VoxSaf> No more candidate    539          G4cout << " VoxSaf> No more candidates:  nodeDown= " << nextDown
538                 << " nodeUp= " << nextUp          540                 << " nodeUp= " << nextUp
539                 << " lastSlice= " << targetHea    541                 << " lastSlice= " << targetHeaderNoSlices << G4endl;
540        }                                          542        }
541 #endif                                            543 #endif
542      }                                            544      }
543                                                   545 
544      // This calculation can be 'hauled'-up to    546      // This calculation can be 'hauled'-up to where minSafety is calculated
545      //                                           547      //
546      distMaxInterest = std::min( minSafety, ma    548      distMaxInterest = std::min( minSafety, maxLength ); 
547                                                   549 
548   } while ( nextIsInside && ( distAxis*distAxi    550   } while ( nextIsInside && ( distAxis*distAxis + distUpperDepth_Sq
549                             < distMaxInterest*    551                             < distMaxInterest*distMaxInterest ) ); 
550                                                   552 
551 #ifdef G4VERBOSE                                  553 #ifdef G4VERBOSE
552   if( fVerbose > 0 )                              554   if( fVerbose > 0 )
553   {                                               555   { 
554     G4cout << " Ended for targetNodeNo -- chec    556     G4cout << " Ended for targetNodeNo -- checked " << checkedNum << " out of "
555            << targetHeaderNoSlices << " slices    557            << targetHeaderNoSlices << " slices." << G4endl;
556     G4cout << " ===== Returning from SafetyFor    558     G4cout << " ===== Returning from SafetyForVoxelHeader " 
557            << "  Depth= " << fVoxelDepth << G4    559            << "  Depth= " << fVoxelDepth << G4endl
558            << G4endl;                             560            << G4endl;  
559   }                                               561   }
560 #endif                                            562 #endif
561                                                   563 
562   // Go back one level                            564   // Go back one level
563   fVoxelDepth--;                                  565   fVoxelDepth--; 
564                                                   566   
565   return ourSafety;                               567   return ourSafety;
566 }                                                 568 }
567                                                   569