Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/navigation/src/G4VoxelNavigation.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/G4VoxelNavigation.cc (Version 11.3.0) and /geometry/navigation/src/G4VoxelNavigation.cc (Version 6.2.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
                                                   >>  23 //
                                                   >>  24 // $Id: G4VoxelNavigation.cc,v 1.4 2004/05/17 13:30:26 gcosmo Exp $
                                                   >>  25 // GEANT4 tag $Name: geant4-06-02 $
                                                   >>  26 //
                                                   >>  27 //
 26 // class G4VoxelNavigation Implementation          28 // class G4VoxelNavigation Implementation
 27 //                                                 29 //
 28 // Author: P.Kent, 1996                            30 // Author: P.Kent, 1996
 29 //                                                 31 //
 30 // -------------------------------------------     32 // --------------------------------------------------------------------
 31 #include "G4VoxelNavigation.hh"                << 
 32 #include "G4GeometryTolerance.hh"              << 
 33 #include "G4VoxelSafety.hh"                    << 
 34                                                << 
 35 #include "G4AuxiliaryNavServices.hh"           << 
 36                                                    33 
 37 #include <cassert>                             <<  34 #include "G4VoxelNavigation.hh"
 38 #include <ostream>                             << 
 39                                                    35 
 40 // *******************************************     36 // ********************************************************************
 41 // Constructor                                     37 // Constructor
 42 // *******************************************     38 // ********************************************************************
 43 //                                                 39 //
 44 G4VoxelNavigation::G4VoxelNavigation()             40 G4VoxelNavigation::G4VoxelNavigation()
 45   : fVoxelAxisStack(kNavigatorVoxelStackMax,kX <<  41   : fVoxelDepth(-1),
                                                   >>  42     fVoxelAxisStack(kNavigatorVoxelStackMax,kXAxis),
 46     fVoxelNoSlicesStack(kNavigatorVoxelStackMa     43     fVoxelNoSlicesStack(kNavigatorVoxelStackMax,0),
 47     fVoxelSliceWidthStack(kNavigatorVoxelStack     44     fVoxelSliceWidthStack(kNavigatorVoxelStackMax,0.),
 48     fVoxelNodeNoStack(kNavigatorVoxelStackMax,     45     fVoxelNodeNoStack(kNavigatorVoxelStackMax,0),
 49     fVoxelHeaderStack(kNavigatorVoxelStackMax, <<  46     fVoxelHeaderStack(kNavigatorVoxelStackMax,(G4SmartVoxelHeader*)0),
                                                   >>  47     fVoxelNode(0),
                                                   >>  48     fCheck(false),
                                                   >>  49     fVerbose(0)
 50 {                                                  50 {
 51   fLogger= new G4NavigationLogger("G4VoxelNavi << 
 52   fpVoxelSafety= new G4VoxelSafety();          << 
 53   fHalfTolerance= 0.5*G4GeometryTolerance::Get << 
 54                                                << 
 55 #ifdef G4DEBUG_NAVIGATION                      << 
 56   SetVerboseLevel(5);   // Reports most about  << 
 57 #endif                                         << 
 58 }                                                  51 }
 59                                                    52 
 60 // *******************************************     53 // ********************************************************************
 61 // Destructor                                      54 // Destructor
 62 // *******************************************     55 // ********************************************************************
 63 //                                                 56 //
 64 G4VoxelNavigation::~G4VoxelNavigation()            57 G4VoxelNavigation::~G4VoxelNavigation()
 65 {                                                  58 {
 66   delete fpVoxelSafety;                        <<  59 #ifdef G4DEBUG_NAVIGATION
 67   delete fLogger;                              <<  60   G4cout << "G4VoxelNavigation::~G4VoxelNavigation() called." << G4endl;
                                                   >>  61 #endif
 68 }                                                  62 }
 69                                                    63 
 70 // ------------------------------------------- << 
 71 // Input:                                      << 
 72 //    exiting:         : last step exited      << 
 73 //    blockedPhysical  : phys volume last exit << 
 74 //    blockedReplicaNo : copy/replica number o << 
 75 // Output:                                     << 
 76 //    entering         : if true, found candid << 
 77 //    blockedPhysical  : candidate phys volume << 
 78 //    blockedReplicaNo : copy/replica number   << 
 79 //    exiting:         : will exit current (mo << 
 80 // In/Out                                      << 
 81 // ------------------------------------------- << 
 82                                                << 
 83 // *******************************************     64 // ********************************************************************
 84 // ComputeStep                                     65 // ComputeStep
 85 // *******************************************     66 // ********************************************************************
 86 //                                                 67 //
 87 G4double                                           68 G4double
 88 G4VoxelNavigation::ComputeStep( const G4ThreeV     69 G4VoxelNavigation::ComputeStep( const G4ThreeVector& localPoint,
 89                                 const G4ThreeV     70                                 const G4ThreeVector& localDirection,
 90                                 const G4double     71                                 const G4double currentProposedStepLength,
 91                                       G4double     72                                       G4double& newSafety,
 92                           /* const */ G4Naviga <<  73                                       G4NavigationHistory& history,
 93                                       G4bool&      74                                       G4bool& validExitNormal,
 94                                       G4ThreeV     75                                       G4ThreeVector& exitNormal,
 95                                       G4bool&      76                                       G4bool& exiting,
 96                                       G4bool&      77                                       G4bool& entering,
 97                                       G4VPhysi <<  78                                       G4VPhysicalVolume *(*pBlockedPhysical),
 98                                       G4int& b     79                                       G4int& blockedReplicaNo )
 99 {                                                  80 {
100   G4VPhysicalVolume *motherPhysical, *samplePh <<  81   G4VPhysicalVolume *motherPhysical, *samplePhysical, *blockedExitedVol=0;
101   G4LogicalVolume *motherLogical;                  82   G4LogicalVolume *motherLogical;
102   G4VSolid *motherSolid;                           83   G4VSolid *motherSolid;
103   G4ThreeVector sampleDirection;                   84   G4ThreeVector sampleDirection;
104   G4double ourStep=currentProposedStepLength,  <<  85   G4double ourStep=currentProposedStepLength, motherSafety, ourSafety;
105   G4double motherSafety, motherStep = DBL_MAX; << 
106   G4int localNoDaughters, sampleNo;                86   G4int localNoDaughters, sampleNo;
107                                                    87 
108   G4bool initialNode, noStep;                      88   G4bool initialNode, noStep;
109   G4SmartVoxelNode *curVoxelNode;                  89   G4SmartVoxelNode *curVoxelNode;
110   G4long curNoVolumes, contentNo;              <<  90   G4int curNoVolumes, contentNo;
111   G4double voxelSafety;                            91   G4double voxelSafety;
112                                                    92 
113   motherPhysical = history.GetTopVolume();         93   motherPhysical = history.GetTopVolume();
114   motherLogical = motherPhysical->GetLogicalVo     94   motherLogical = motherPhysical->GetLogicalVolume();
115   motherSolid = motherLogical->GetSolid();         95   motherSolid = motherLogical->GetSolid();
116                                                    96 
117   //                                               97   //
118   // Compute mother safety                         98   // Compute mother safety
119   //                                               99   //
120                                                   100 
121   motherSafety = motherSolid->DistanceToOut(lo    101   motherSafety = motherSolid->DistanceToOut(localPoint);
122   ourSafety = motherSafety;                 //    102   ourSafety = motherSafety;                 // Working isotropic safety
123                                                   103   
124 #ifdef G4VERBOSE                                  104 #ifdef G4VERBOSE
125   if ( fCheck )                                   105   if ( fCheck )
126   {                                               106   {
127     fLogger->PreComputeStepLog (motherPhysical << 107     if(fVerbose == 1 )
                                                   >> 108     {
                                                   >> 109       G4cout << "*** G4VoxelNavigation::ComputeStep(): ***" << G4endl
                                                   >> 110              << "    Invoked DistanceToOut(p) for mother solid: "
                                                   >> 111              << motherSolid->GetName()
                                                   >> 112              << ". Solid replied: " << motherSafety << G4endl
                                                   >> 113              << "    For local point p: " << localPoint
                                                   >> 114              << ", to be considered as 'mother safety'." << G4endl;
                                                   >> 115     }
                                                   >> 116     if( motherSafety < 0.0 )
                                                   >> 117     {
                                                   >> 118       G4cout << "ERROR - G4VoxelNavigation::ComputeStep()" << G4endl
                                                   >> 119              << "        Current solid " << motherSolid->GetName()
                                                   >> 120              << " gave negative safety: " << motherSafety << G4endl
                                                   >> 121              << "        for the current (local) point " << localPoint
                                                   >> 122              << G4endl;
                                                   >> 123       G4Exception("G4VoxelNavigation::ComputeStep()",
                                                   >> 124                   "NegativeSafetyMotherVol", FatalException,
                                                   >> 125                   "Negative Safety In Voxel Navigation !" ); 
                                                   >> 126     }
                                                   >> 127     if( motherSolid->Inside(localPoint)==kOutside )
                                                   >> 128     { 
                                                   >> 129       G4cout << "WARNING - G4VoxelNavigation::ComputeStep()" << G4endl
                                                   >> 130              << "          Point " << localPoint
                                                   >> 131              << " is outside current volume " << motherPhysical->GetName()
                                                   >> 132              << G4endl;
                                                   >> 133       G4double  estDistToSolid= motherSolid->DistanceToIn(localPoint); 
                                                   >> 134       G4cout << "          Estimated isotropic distance to solid (distToIn)= " 
                                                   >> 135              << estDistToSolid << G4endl;
                                                   >> 136       if( estDistToSolid > 100.0 * kCarTolerance ) 
                                                   >> 137         G4Exception("G4VoxelNavigation::ComputeStep()",
                                                   >> 138                     "FarOutsideCurrentVolume", FatalException,
                                                   >> 139                     "Point is far outside Current Volume !"); 
                                                   >> 140       else
                                                   >> 141         G4Exception("G4VoxelNavigation::ComputeStep()", "OutsideCurrentVolume", 
                                                   >> 142                     JustWarning, "Point is a little outside Current Volume."); 
                                                   >> 143     }
128   }                                               144   }
129 #endif                                            145 #endif
130                                                   146 
131   //                                              147   //
132   // Compute daughter safeties & intersections    148   // Compute daughter safeties & intersections
133   //                                              149   //
134                                                   150 
135   // Exiting normal optimisation                  151   // Exiting normal optimisation
136   //                                              152   //
137   if ( exiting && validExitNormal )               153   if ( exiting && validExitNormal )
138   {                                               154   {
139     if ( localDirection.dot(exitNormal)>=kMinE    155     if ( localDirection.dot(exitNormal)>=kMinExitingNormalCosine )
140     {                                             156     {
141       // Block exited daughter volume             157       // Block exited daughter volume
142       //                                          158       //
143       blockedExitedVol = *pBlockedPhysical;       159       blockedExitedVol = *pBlockedPhysical;
144       ourSafety = 0;                              160       ourSafety = 0;
145     }                                             161     }
146   }                                               162   }
147   exiting = false;                                163   exiting = false;
148   entering = false;                               164   entering = false;
149                                                   165 
150   // For extra checking,  get the distance to  << 166   localNoDaughters = motherLogical->GetNoDaughters();
151   G4bool motherValidExitNormal = false;        << 
152   G4ThreeVector motherExitNormal(0.0, 0.0, 0.0 << 
153                                                << 
154 #ifdef G4VERBOSE                               << 
155   if ( fCheck )                                << 
156   {                                            << 
157     // Compute early -- a) for validity        << 
158     //                  b) to check against an << 
159     motherStep = motherSolid->DistanceToOut(lo << 
160                                             lo << 
161                                             tr << 
162                                            &mo << 
163                                            &mo << 
164   }                                            << 
165 #endif                                         << 
166                                                << 
167   localNoDaughters = (G4int)motherLogical->Get << 
168                                                   167 
169   fBList.Enlarge(localNoDaughters);               168   fBList.Enlarge(localNoDaughters);
170   fBList.Reset();                                 169   fBList.Reset();
171                                                   170 
172   initialNode = true;                             171   initialNode = true;
173   noStep = true;                                  172   noStep = true;
174                                                   173 
175   while (noStep)                                  174   while (noStep)
176   {                                               175   {
177     curVoxelNode = fVoxelNode;                    176     curVoxelNode = fVoxelNode;
178     curNoVolumes = curVoxelNode->GetNoContaine    177     curNoVolumes = curVoxelNode->GetNoContained();
179     for (contentNo=curNoVolumes-1; contentNo>=    178     for (contentNo=curNoVolumes-1; contentNo>=0; contentNo--)
180     {                                             179     {
181       sampleNo = curVoxelNode->GetVolume((G4in << 180       sampleNo = curVoxelNode->GetVolume(contentNo);
182       if ( !fBList.IsBlocked(sampleNo) )          181       if ( !fBList.IsBlocked(sampleNo) )
183       {                                           182       {
184         fBList.BlockVolume(sampleNo);             183         fBList.BlockVolume(sampleNo);
185         samplePhysical = motherLogical->GetDau    184         samplePhysical = motherLogical->GetDaughter(sampleNo);
186         if ( samplePhysical!=blockedExitedVol     185         if ( samplePhysical!=blockedExitedVol )
187         {                                         186         {
188           G4AffineTransform sampleTf(samplePhy    187           G4AffineTransform sampleTf(samplePhysical->GetRotation(),
189                                      samplePhy    188                                      samplePhysical->GetTranslation());
190           sampleTf.Invert();                      189           sampleTf.Invert();
191           const G4ThreeVector samplePoint =       190           const G4ThreeVector samplePoint =
192                      sampleTf.TransformPoint(l    191                      sampleTf.TransformPoint(localPoint);
193           const G4VSolid *sampleSolid     =       192           const G4VSolid *sampleSolid     =
194                      samplePhysical->GetLogica    193                      samplePhysical->GetLogicalVolume()->GetSolid();
195           const G4double sampleSafety     =       194           const G4double sampleSafety     =
196                      sampleSolid->DistanceToIn    195                      sampleSolid->DistanceToIn(samplePoint);
197                                                << 196 #ifdef G4VERBOSE
                                                   >> 197           if(( fCheck ) && ( fVerbose == 1 ))
                                                   >> 198           {
                                                   >> 199             G4cout << "*** G4VoxelNavigation::ComputeStep(): ***" << G4endl
                                                   >> 200                    << "    Invoked DistanceToIn(p) for daughter solid: "
                                                   >> 201                    << sampleSolid->GetName()
                                                   >> 202                    << ". Solid replied: " << sampleSafety << G4endl
                                                   >> 203                    << "    For local point p: " << samplePoint
                                                   >> 204                    << ", to be considered as 'daughter safety'." << G4endl;
                                                   >> 205           }
                                                   >> 206 #endif
198           if ( sampleSafety<ourSafety )           207           if ( sampleSafety<ourSafety )
199           {                                       208           {
200             ourSafety = sampleSafety;             209             ourSafety = sampleSafety;
201           }                                       210           }
202           if ( sampleSafety<=ourStep )            211           if ( sampleSafety<=ourStep )
203           {                                       212           {
204             sampleDirection = sampleTf.Transfo    213             sampleDirection = sampleTf.TransformAxis(localDirection);
205             G4double sampleStep =                 214             G4double sampleStep =
206                      sampleSolid->DistanceToIn    215                      sampleSolid->DistanceToIn(samplePoint, sampleDirection);
207 #ifdef G4VERBOSE                                  216 #ifdef G4VERBOSE
208             if( fCheck )                       << 217             if(( fCheck ) && ( fVerbose == 1 ))
209             {                                     218             {
210               fLogger->PrintDaughterLog(sample << 219               G4cout << "*** G4VoxelNavigation::ComputeStep(): ***" << G4endl
211                                         sample << 220                      << "    Invoked DistanceToIn(p,v) for daughter solid: "
212                                         sample << 221                      << sampleSolid->GetName()
                                                   >> 222                      << ". Solid replied: " << sampleStep << G4endl
                                                   >> 223                      << "    For local point p: " << samplePoint << G4endl
                                                   >> 224                      << "    Direction v: " << sampleDirection
                                                   >> 225                      << ", to be considered as 'daughter step'." << G4endl;
213             }                                     226             }
214 #endif                                            227 #endif
215             if ( sampleStep<=ourStep )            228             if ( sampleStep<=ourStep )
216             {                                     229             {
217               ourStep = sampleStep;               230               ourStep = sampleStep;
218               entering = true;                    231               entering = true;
219               exiting = false;                    232               exiting = false;
220               *pBlockedPhysical = samplePhysic    233               *pBlockedPhysical = samplePhysical;
221               blockedReplicaNo = -1;              234               blockedReplicaNo = -1;
222 #ifdef G4VERBOSE                                  235 #ifdef G4VERBOSE
223               // Check to see that the resulti    236               // Check to see that the resulting point is indeed in/on volume.
224               // This could be done only for s << 237               // This check could eventually be made only for successful
225               if ( fCheck )                    << 238               // candidate.
                                                   >> 239 
                                                   >> 240               if ( ( fCheck ) && ( sampleStep < kInfinity ) )
226               {                                   241               {
227                 fLogger->AlongComputeStepLog ( << 242                 G4ThreeVector intersectionPoint;
228                   sampleDirection, localDirect << 243                 intersectionPoint= samplePoint + sampleStep * sampleDirection;
                                                   >> 244                 EInside insideIntPt= sampleSolid->Inside(intersectionPoint); 
                                                   >> 245                 G4String solidResponse = "-kInside-";
                                                   >> 246                 if (insideIntPt == kOutside)
                                                   >> 247                   solidResponse = "-kOutside-";
                                                   >> 248                 else if (insideIntPt == kSurface)
                                                   >> 249                   solidResponse = "-kSurface-";
                                                   >> 250                 if( fVerbose == 1 )
                                                   >> 251                 {
                                                   >> 252                   G4cout << "*** G4VoxelNavigation::ComputeStep(): ***"<<G4endl
                                                   >> 253                          << "    Invoked Inside() for solid: "
                                                   >> 254                          << sampleSolid->GetName()
                                                   >> 255                          << ". Solid replied: " << solidResponse << G4endl
                                                   >> 256                          << "    For point p: " << intersectionPoint
                                                   >> 257                          << ", considered as 'intersection' point." << G4endl;
                                                   >> 258                 }
                                                   >> 259                 if( insideIntPt != kSurface )
                                                   >> 260                 {
                                                   >> 261                   G4int oldcoutPrec = G4cout.precision(16); 
                                                   >> 262                   G4cout << "WARNING - G4VoxelNavigation::ComputeStep()"
                                                   >> 263                          << G4endl
                                                   >> 264                          << "          Inaccurate solid DistanceToIn"
                                                   >> 265                          << " for solid " << sampleSolid->GetName() << G4endl;
                                                   >> 266                   G4cout << "          Solid gave DistanceToIn = "
                                                   >> 267                          << sampleStep << " yet returns " << solidResponse
                                                   >> 268                          << " for this point !" << G4endl; 
                                                   >> 269                   G4cout << "          Point = " << intersectionPoint << G4endl;
                                                   >> 270                   if ( insideIntPt != kInside )
                                                   >> 271                     G4cout << "        DistanceToIn(p) = " 
                                                   >> 272                            << sampleSolid->DistanceToIn(intersectionPoint)
                                                   >> 273                            << G4endl;
                                                   >> 274                   if ( insideIntPt != kOutside ) 
                                                   >> 275                     G4cout << "        DistanceToOut(p) = " 
                                                   >> 276                            << sampleSolid->DistanceToOut(intersectionPoint)
                                                   >> 277                            << G4endl;
                                                   >> 278                   G4Exception("G4VoxelNavigation::ComputeStep()", 
                                                   >> 279                               "InaccurateDistanceToIn", JustWarning,
                                                   >> 280                               "Navigator gets conflicting response from Solid.");
                                                   >> 281                   G4cout.precision(oldcoutPrec);
                                                   >> 282                 }
229               }                                   283               }
230 #endif                                            284 #endif
231             }                                     285             }
232 #ifdef G4VERBOSE                               << 
233             if ( fCheck && ( sampleStep < kInf << 
234                         && ( sampleStep >= mot << 
235             {                                  << 
236                // The intersection point with  << 
237                // point from the mother volume << 
238                fLogger->CheckDaughterEntryPoin << 
239                                                << 
240                                                << 
241                                                << 
242                                                << 
243             }                                  << 
244 #endif                                         << 
245           }                                    << 
246 #ifdef G4VERBOSE                               << 
247           else // ie if sampleSafety > outStep << 
248           {                                    << 
249             if( fCheck )                       << 
250             {                                  << 
251               fLogger->PrintDaughterLog(sample << 
252                                         sample << 
253                                         G4Thre << 
254             }                                  << 
255           }                                       286           }
256 #endif                                         << 
257         }                                         287         }
258       }                                           288       }
259     }                                             289     }
260     if (initialNode)                              290     if (initialNode)
261     {                                             291     {
262       initialNode = false;                        292       initialNode = false;
263       voxelSafety = ComputeVoxelSafety(localPo    293       voxelSafety = ComputeVoxelSafety(localPoint);
264       if ( voxelSafety<ourSafety )                294       if ( voxelSafety<ourSafety )
265       {                                           295       {
266         ourSafety = voxelSafety;                  296         ourSafety = voxelSafety;
267       }                                           297       }
268       if ( currentProposedStepLength<ourSafety    298       if ( currentProposedStepLength<ourSafety )
269       {                                           299       {
270         // Guaranteed physics limited             300         // Guaranteed physics limited
271         //                                        301         //      
272         noStep = false;                           302         noStep = false;
273         entering = false;                         303         entering = false;
274         exiting = false;                          304         exiting = false;
275         *pBlockedPhysical = nullptr;           << 305         *pBlockedPhysical = 0;
276         ourStep = kInfinity;                      306         ourStep = kInfinity;
277       }                                           307       }
278       else                                        308       else
279       {                                           309       {
280         //                                        310         //
281         // Compute mother intersection if requ    311         // Compute mother intersection if required
282         //                                        312         //
283         if ( motherSafety<=ourStep )              313         if ( motherSafety<=ourStep )
284         {                                         314         {
285           // In case of check mode this is a d << 315           G4double motherStep =
286           motherStep = motherSolid->DistanceTo << 316               motherSolid->DistanceToOut(localPoint,
287                               true, &motherVal << 317                                          localDirection,
                                                   >> 318                                          true, &validExitNormal, &exitNormal);
288 #ifdef G4VERBOSE                                  319 #ifdef G4VERBOSE
289           if ( fCheck )                           320           if ( fCheck )
290           {                                       321           {
291             fLogger->PostComputeStepLog(mother << 322             if(fVerbose == 1)
292                                         mother << 
293             if( motherValidExitNormal )        << 
294             {                                     323             {
295               fLogger->CheckAndReportBadNormal << 324               G4cout << "*** G4VoxelNavigation::ComputeStep(): ***" << G4endl
296                                                << 325                      << "    Invoked DistanceToOut(p,v,...) for mother solid: "
297                                                << 326                      << motherSolid->GetName()
298                                         "From  << 327                      << ". Solid replied: " << motherStep << G4endl
                                                   >> 328                      << "    For local point p: " << localPoint << G4endl
                                                   >> 329                      << "    Direction v: " << localDirection
                                                   >> 330                      << ", to be considered as 'mother step'." << G4endl;
299             }                                     331             }
300           }                                    << 332             if( ( motherStep < 0.0 ) || ( motherStep >= kInfinity) )
301 #endif                                         << 
302           if( (motherStep >= kInfinity) || (mo << 
303           {                                    << 
304 #ifdef G4VERBOSE                               << 
305             if( fCheck ) // Error - indication << 
306             {                                     333             {
307               fLogger->ReportOutsideMother(loc << 334               G4int oldPrOut= G4cout.precision(16); 
308                                            mot << 335               G4int oldPrErr= G4cerr.precision(16);
                                                   >> 336               G4cerr << "ERROR - G4VoxelNavigation::ComputeStep()" << G4endl
                                                   >> 337                      << "        Problem in Navigation"  << G4endl
                                                   >> 338                      << "        Point (local coordinates): "
                                                   >> 339                      << localPoint << G4endl
                                                   >> 340                      << "        Local Direction: " << localDirection << G4endl
                                                   >> 341                      << "        Solid: " << motherSolid->GetName() << G4endl; 
                                                   >> 342               G4Exception("G4VoxelNavigation::ComputeStep()",
                                                   >> 343                           "PointOutsideCurrentVolume", FatalException,
                                                   >> 344                           "Current point is outside the current solid !");
                                                   >> 345               G4cout.precision(oldPrOut);
                                                   >> 346               G4cerr.precision(oldPrErr);
309             }                                     347             }
                                                   >> 348           }
310 #endif                                            349 #endif
311             motherStep = 0.0;                  << 
312             ourStep = 0.0;                     << 
313             exiting = true;                    << 
314             entering = false;                  << 
315                                                << 
316             // validExitNormal= motherValidExi << 
317             // exitNormal= motherExitNormal;   << 
318             // Useful only if the point is ver << 
319             // => but it would need to be rota << 
320             validExitNormal= false;            << 
321                                                << 
322             *pBlockedPhysical = nullptr; // or << 
323             blockedReplicaNo = 0;  // or mothe << 
324                                                << 
325             newSafety = 0.0;                   << 
326             return ourStep;                    << 
327           }                                    << 
328                                                << 
329           if ( motherStep<=ourStep )              350           if ( motherStep<=ourStep )
330           {                                       351           {
331             ourStep = motherStep;                 352             ourStep = motherStep;
332             exiting = true;                       353             exiting = true;
333             entering = false;                     354             entering = false;
334                                                << 
335             // Exit normal: Natural location t << 
336             //                                 << 
337             validExitNormal = motherValidExitN << 
338             exitNormal = motherExitNormal;     << 
339                                                << 
340             if ( validExitNormal )                355             if ( validExitNormal )
341             {                                     356             {
342               const G4RotationMatrix *rot = mo    357               const G4RotationMatrix *rot = motherPhysical->GetRotation();
343               if (rot != nullptr)              << 358               if (rot)
344               {                                   359               {
345                 exitNormal *= rot->inverse();     360                 exitNormal *= rot->inverse();
346 #ifdef G4VERBOSE                               << 
347                 if( fCheck )                   << 
348                 {                              << 
349                   fLogger->CheckAndReportBadNo << 
350                                                << 
351                                                << 
352                                                << 
353                 }                              << 
354 #endif                                         << 
355               }                                   361               }
356             }                                  << 362             }  
357           }                                       363           }
358           else                                    364           else
359           {                                       365           {
360             validExitNormal = false;              366             validExitNormal = false;
361           }                                       367           }
362         }                                         368         }
363       }                                           369       }
364       newSafety = ourSafety;                      370       newSafety = ourSafety;
365     }                                             371     }
366     if (noStep)                                   372     if (noStep)
367     {                                             373     {
368       noStep = LocateNextVoxel(localPoint, loc    374       noStep = LocateNextVoxel(localPoint, localDirection, ourStep);
369     }                                             375     }
370   }  // end -while (noStep)- loop                 376   }  // end -while (noStep)- loop
371                                                   377 
372   return ourStep;                                 378   return ourStep;
373 }                                                 379 }
374                                                   380 
375 // *******************************************    381 // ********************************************************************
376 // ComputeVoxelSafety                             382 // ComputeVoxelSafety
377 //                                                383 //
378 // Computes safety from specified point to vox    384 // Computes safety from specified point to voxel boundaries
379 // using already located point                    385 // using already located point
380 // o collected boundaries for most derived lev    386 // o collected boundaries for most derived level
381 // o adjacent boundaries for previous levels      387 // o adjacent boundaries for previous levels
382 // *******************************************    388 // ********************************************************************
383 //                                                389 //
384 G4double                                          390 G4double
385 G4VoxelNavigation::ComputeVoxelSafety(const G4    391 G4VoxelNavigation::ComputeVoxelSafety(const G4ThreeVector& localPoint) const
386 {                                                 392 {
387   G4SmartVoxelHeader *curHeader;                  393   G4SmartVoxelHeader *curHeader;
388   G4double voxelSafety, curNodeWidth;             394   G4double voxelSafety, curNodeWidth;
389   G4double curNodeOffset, minCurCommonDelta, m    395   G4double curNodeOffset, minCurCommonDelta, maxCurCommonDelta;
390   G4int minCurNodeNoDelta, maxCurNodeNoDelta;     396   G4int minCurNodeNoDelta, maxCurNodeNoDelta;
391   G4int localVoxelDepth, curNodeNo;               397   G4int localVoxelDepth, curNodeNo;
392   EAxis curHeaderAxis;                            398   EAxis curHeaderAxis;
393                                                   399 
394   localVoxelDepth = fVoxelDepth;                  400   localVoxelDepth = fVoxelDepth;
395                                                   401 
396   curHeader = fVoxelHeaderStack[localVoxelDept    402   curHeader = fVoxelHeaderStack[localVoxelDepth];
397   curHeaderAxis = fVoxelAxisStack[localVoxelDe    403   curHeaderAxis = fVoxelAxisStack[localVoxelDepth];
398   curNodeNo = fVoxelNodeNoStack[localVoxelDept    404   curNodeNo = fVoxelNodeNoStack[localVoxelDepth];
399   curNodeWidth = fVoxelSliceWidthStack[localVo    405   curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth];
400                                                   406   
401   // Compute linear intersection distance to b    407   // Compute linear intersection distance to boundaries of max/min
402   // to collected nodes at current level          408   // to collected nodes at current level
403   //                                              409   //
404   curNodeOffset = curNodeNo*curNodeWidth;         410   curNodeOffset = curNodeNo*curNodeWidth;
405   maxCurNodeNoDelta = fVoxelNode->GetMaxEquiva    411   maxCurNodeNoDelta = fVoxelNode->GetMaxEquivalentSliceNo()-curNodeNo;
406   minCurNodeNoDelta = curNodeNo-fVoxelNode->Ge    412   minCurNodeNoDelta = curNodeNo-fVoxelNode->GetMinEquivalentSliceNo();
407   minCurCommonDelta = localPoint(curHeaderAxis    413   minCurCommonDelta = localPoint(curHeaderAxis)
408                       - curHeader->GetMinExten    414                       - curHeader->GetMinExtent() - curNodeOffset;
409   maxCurCommonDelta = curNodeWidth-minCurCommo    415   maxCurCommonDelta = curNodeWidth-minCurCommonDelta;
410                                                   416 
411   if ( minCurNodeNoDelta<maxCurNodeNoDelta )      417   if ( minCurNodeNoDelta<maxCurNodeNoDelta )
412   {                                               418   {
413     voxelSafety = minCurNodeNoDelta*curNodeWid    419     voxelSafety = minCurNodeNoDelta*curNodeWidth;
414     voxelSafety += minCurCommonDelta;             420     voxelSafety += minCurCommonDelta;
415   }                                               421   }
416   else if (maxCurNodeNoDelta < minCurNodeNoDel    422   else if (maxCurNodeNoDelta < minCurNodeNoDelta)
417   {                                            << 423        {
418     voxelSafety = maxCurNodeNoDelta*curNodeWid << 424          voxelSafety = maxCurNodeNoDelta*curNodeWidth;
419     voxelSafety += maxCurCommonDelta;          << 425          voxelSafety += maxCurCommonDelta;
420   }                                            << 426         }
421   else    // (maxCurNodeNoDelta == minCurNodeN << 427         else    // (maxCurNodeNoDelta == minCurNodeNoDelta)
422   {                                            << 428         {
423     voxelSafety = minCurNodeNoDelta*curNodeWid << 429           voxelSafety = minCurNodeNoDelta*curNodeWidth;
424     voxelSafety += std::min(minCurCommonDelta, << 430           voxelSafety += std::min(minCurCommonDelta,maxCurCommonDelta);
425   }                                            << 431         }
426                                                   432 
427   // Compute isotropic safety to boundaries of    433   // Compute isotropic safety to boundaries of previous levels
428   // [NOT to collected boundaries]                434   // [NOT to collected boundaries]
429                                                << 435   //
430   // Loop checking, 07.10.2016, JA             << 
431   while ( (localVoxelDepth>0) && (voxelSafety>    436   while ( (localVoxelDepth>0) && (voxelSafety>0) )
432   {                                               437   {
433     localVoxelDepth--;                            438     localVoxelDepth--;
434     curHeader = fVoxelHeaderStack[localVoxelDe    439     curHeader = fVoxelHeaderStack[localVoxelDepth];
435     curHeaderAxis = fVoxelAxisStack[localVoxel    440     curHeaderAxis = fVoxelAxisStack[localVoxelDepth];
436     curNodeNo = fVoxelNodeNoStack[localVoxelDe    441     curNodeNo = fVoxelNodeNoStack[localVoxelDepth];
437     curNodeWidth = fVoxelSliceWidthStack[local    442     curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth];
438     curNodeOffset = curNodeNo*curNodeWidth;       443     curNodeOffset = curNodeNo*curNodeWidth;
439     minCurCommonDelta = localPoint(curHeaderAx    444     minCurCommonDelta = localPoint(curHeaderAxis)
440                         - curHeader->GetMinExt    445                         - curHeader->GetMinExtent() - curNodeOffset;
441     maxCurCommonDelta = curNodeWidth-minCurCom    446     maxCurCommonDelta = curNodeWidth-minCurCommonDelta;
442                                                   447     
443     if ( minCurCommonDelta<voxelSafety )          448     if ( minCurCommonDelta<voxelSafety )
444     {                                             449     {
445       voxelSafety = minCurCommonDelta;            450       voxelSafety = minCurCommonDelta;
446     }                                             451     }
447     if ( maxCurCommonDelta<voxelSafety )          452     if ( maxCurCommonDelta<voxelSafety )
448     {                                             453     {
449       voxelSafety = maxCurCommonDelta;            454       voxelSafety = maxCurCommonDelta;
450     }                                             455     }
451   }                                               456   }
452   if ( voxelSafety<0 )                            457   if ( voxelSafety<0 )
453   {                                               458   {
454     voxelSafety = 0;                              459     voxelSafety = 0;
455   }                                               460   }
456                                                   461 
457   return voxelSafety;                             462   return voxelSafety;
458 }                                                 463 }
459                                                   464 
460 // *******************************************    465 // ********************************************************************
461 // LocateNextVoxel                                466 // LocateNextVoxel
462 //                                                467 //
463 // Finds the next voxel from the current voxel    468 // Finds the next voxel from the current voxel and point
464 // in the specified direction                     469 // in the specified direction
465 //                                                470 //
466 // Returns false if all voxels considered         471 // Returns false if all voxels considered
467 //              [current Step ends inside same    472 //              [current Step ends inside same voxel or leaves all voxels]
468 //         true  otherwise                        473 //         true  otherwise
469 //              [the information on the next v    474 //              [the information on the next voxel is put into the set of
470 //               fVoxel* variables & "stacks"]    475 //               fVoxel* variables & "stacks"] 
471 // *******************************************    476 // ********************************************************************
472 //                                                477 // 
473 G4bool                                            478 G4bool
474 G4VoxelNavigation::LocateNextVoxel(const G4Thr    479 G4VoxelNavigation::LocateNextVoxel(const G4ThreeVector& localPoint,
475                                    const G4Thr    480                                    const G4ThreeVector& localDirection,
476                                    const G4dou    481                                    const G4double currentStep)
477 {                                                 482 {
478   G4SmartVoxelHeader *workHeader=nullptr, *new << 483   G4SmartVoxelHeader *workHeader=0, *newHeader=0;
479   G4SmartVoxelProxy *newProxy=nullptr;         << 484   G4SmartVoxelProxy *newProxy=0;
480   G4SmartVoxelNode *newVoxelNode=nullptr;      << 485   G4SmartVoxelNode *newVoxelNode=0;
481   G4ThreeVector targetPoint, voxelPoint;          486   G4ThreeVector targetPoint, voxelPoint;
482   G4double workNodeWidth, workMinExtent, workC    487   G4double workNodeWidth, workMinExtent, workCoord;
483   G4double minVal, maxVal, newDistance=0.;        488   G4double minVal, maxVal, newDistance=0.;
484   G4double newHeaderMin, newHeaderNodeWidth;      489   G4double newHeaderMin, newHeaderNodeWidth;
485   G4int depth=0, newDepth=0, workNodeNo=0, new    490   G4int depth=0, newDepth=0, workNodeNo=0, newNodeNo=0, newHeaderNoSlices=0;
486   EAxis workHeaderAxis, newHeaderAxis;            491   EAxis workHeaderAxis, newHeaderAxis;
487   G4bool isNewVoxel = false;                   << 492   G4bool isNewVoxel=false;
488                                                   493   
489   G4double currentDistance = currentStep;         494   G4double currentDistance = currentStep;
490                                                   495 
491   // Determine if end of Step within current v    496   // Determine if end of Step within current voxel
492   //                                              497   //
493   for (depth=0; depth<fVoxelDepth; ++depth)    << 498   for (depth=0; depth<fVoxelDepth; depth++)
494   {                                               499   {
495     targetPoint = localPoint+localDirection*cu    500     targetPoint = localPoint+localDirection*currentDistance;
496     newDistance = currentDistance;                501     newDistance = currentDistance;
497     workHeader = fVoxelHeaderStack[depth];        502     workHeader = fVoxelHeaderStack[depth];
498     workHeaderAxis = fVoxelAxisStack[depth];      503     workHeaderAxis = fVoxelAxisStack[depth];
499     workNodeNo = fVoxelNodeNoStack[depth];        504     workNodeNo = fVoxelNodeNoStack[depth];
500     workNodeWidth = fVoxelSliceWidthStack[dept    505     workNodeWidth = fVoxelSliceWidthStack[depth];
501     workMinExtent = workHeader->GetMinExtent()    506     workMinExtent = workHeader->GetMinExtent();
502     workCoord = targetPoint(workHeaderAxis);      507     workCoord = targetPoint(workHeaderAxis);
503     minVal = workMinExtent+workNodeNo*workNode    508     minVal = workMinExtent+workNodeNo*workNodeWidth;
504                                                   509 
505     if ( minVal<=workCoord+fHalfTolerance )    << 510     if ( minVal<=workCoord+kCarTolerance*0.5 )
506     {                                             511     {
507       maxVal = minVal+workNodeWidth;              512       maxVal = minVal+workNodeWidth;
508       if ( maxVal<=workCoord-fHalfTolerance )  << 513       if ( maxVal<=workCoord-kCarTolerance*0.5 )
509       {                                           514       {
510         // Must consider next voxel               515         // Must consider next voxel
511         //                                        516         //
512         newNodeNo = workNodeNo+1;                 517         newNodeNo = workNodeNo+1;
513         newHeader = workHeader;                   518         newHeader = workHeader;
514         newDistance = (maxVal-localPoint(workH    519         newDistance = (maxVal-localPoint(workHeaderAxis))
515                     / localDirection(workHeade    520                     / localDirection(workHeaderAxis);
516         isNewVoxel = true;                        521         isNewVoxel = true;
517         newDepth = depth;                         522         newDepth = depth;
518       }                                           523       }
519     }                                             524     }
520     else                                          525     else
521     {                                             526     {
522       newNodeNo = workNodeNo-1;                   527       newNodeNo = workNodeNo-1;
523       newHeader = workHeader;                     528       newHeader = workHeader;
524       newDistance = (minVal-localPoint(workHea    529       newDistance = (minVal-localPoint(workHeaderAxis))
525                   / localDirection(workHeaderA    530                   / localDirection(workHeaderAxis);
526       isNewVoxel = true;                          531       isNewVoxel = true;
527       newDepth = depth;                           532       newDepth = depth;
528     }                                             533     }
529     currentDistance = newDistance;                534     currentDistance = newDistance;
530   }                                               535   }
531   targetPoint = localPoint+localDirection*curr    536   targetPoint = localPoint+localDirection*currentDistance;
532                                                   537 
533   // Check if end of Step within collected bou    538   // Check if end of Step within collected boundaries of current voxel
534   //                                              539   //
535   depth = fVoxelDepth;                            540   depth = fVoxelDepth;
536   {                                               541   {
537     workHeader = fVoxelHeaderStack[depth];        542     workHeader = fVoxelHeaderStack[depth];
538     workHeaderAxis = fVoxelAxisStack[depth];      543     workHeaderAxis = fVoxelAxisStack[depth];
539     workNodeNo = fVoxelNodeNoStack[depth];        544     workNodeNo = fVoxelNodeNoStack[depth];
540     workNodeWidth = fVoxelSliceWidthStack[dept    545     workNodeWidth = fVoxelSliceWidthStack[depth];
541     workMinExtent = workHeader->GetMinExtent()    546     workMinExtent = workHeader->GetMinExtent();
542     workCoord = targetPoint(workHeaderAxis);      547     workCoord = targetPoint(workHeaderAxis);
543     minVal = workMinExtent+fVoxelNode->GetMinE    548     minVal = workMinExtent+fVoxelNode->GetMinEquivalentSliceNo()*workNodeWidth;
544                                                   549 
545     if ( minVal<=workCoord+fHalfTolerance )    << 550     if ( minVal<=workCoord+kCarTolerance*0.5 )
546     {                                             551     {
547       maxVal = workMinExtent+(fVoxelNode->GetM    552       maxVal = workMinExtent+(fVoxelNode->GetMaxEquivalentSliceNo()+1)
548                             *workNodeWidth;       553                             *workNodeWidth;
549       if ( maxVal<=workCoord-fHalfTolerance )  << 554       if ( maxVal<=workCoord-kCarTolerance*0.5 )
550       {                                           555       {
551         newNodeNo = fVoxelNode->GetMaxEquivale    556         newNodeNo = fVoxelNode->GetMaxEquivalentSliceNo()+1;
552         newHeader = workHeader;                   557         newHeader = workHeader;
553         newDistance = (maxVal-localPoint(workH    558         newDistance = (maxVal-localPoint(workHeaderAxis))
554                     / localDirection(workHeade    559                     / localDirection(workHeaderAxis);
555         isNewVoxel = true;                        560         isNewVoxel = true;
556         newDepth = depth;                         561         newDepth = depth;
557       }                                           562       }
558     }                                             563     }
559     else                                          564     else
560     {                                             565     {
561       newNodeNo = fVoxelNode->GetMinEquivalent    566       newNodeNo = fVoxelNode->GetMinEquivalentSliceNo()-1;
562       newHeader = workHeader;                     567       newHeader = workHeader;
563       newDistance = (minVal-localPoint(workHea    568       newDistance = (minVal-localPoint(workHeaderAxis))
564                   / localDirection(workHeaderA    569                   / localDirection(workHeaderAxis);
565       isNewVoxel = true;                          570       isNewVoxel = true;
566       newDepth = depth;                           571       newDepth = depth;
567     }                                             572     }
568     currentDistance = newDistance;                573     currentDistance = newDistance;
569   }                                               574   }
570   if (isNewVoxel)                                 575   if (isNewVoxel)
571   {                                               576   {
572     // Compute new voxel & adjust voxel stack     577     // Compute new voxel & adjust voxel stack
573     //                                            578     //
574     // newNodeNo=Candidate node no at             579     // newNodeNo=Candidate node no at 
575     // newDepth =refinement depth of crossed v    580     // newDepth =refinement depth of crossed voxel boundary
576     // newHeader=Header for crossed voxel         581     // newHeader=Header for crossed voxel
577     // newDistance=distance to crossed voxel b    582     // newDistance=distance to crossed voxel boundary (along the track)
578     //                                            583     //
579     if ( (newNodeNo<0) || (newNodeNo>=G4int(ne << 584     if ( (newNodeNo<0) || (newNodeNo>=newHeader->GetNoSlices()))
580     {                                             585     {
581       // Leaving mother volume                    586       // Leaving mother volume
582       //                                          587       //
583       isNewVoxel = false;                         588       isNewVoxel = false;
584     }                                             589     }
585     else                                          590     else
586     {                                             591     {
587       // Compute intersection point on the lea    592       // Compute intersection point on the least refined
588       // voxel boundary that is hit               593       // voxel boundary that is hit
589       //                                          594       //
590       voxelPoint = localPoint+localDirection*n    595       voxelPoint = localPoint+localDirection*newDistance;
591       fVoxelNodeNoStack[newDepth] = newNodeNo;    596       fVoxelNodeNoStack[newDepth] = newNodeNo;
592       fVoxelDepth = newDepth;                     597       fVoxelDepth = newDepth;
593       newVoxelNode = nullptr;                  << 598       newVoxelNode = 0;
594       while ( newVoxelNode == nullptr )        << 599       while ( !newVoxelNode )
595       {                                           600       {
596         newProxy = newHeader->GetSlice(newNode    601         newProxy = newHeader->GetSlice(newNodeNo);
597         if (newProxy->IsNode())                   602         if (newProxy->IsNode())
598         {                                         603         {
599           newVoxelNode = newProxy->GetNode();     604           newVoxelNode = newProxy->GetNode();
600         }                                         605         }
601         else                                      606         else
602         {                                         607         {
603           ++fVoxelDepth;                       << 608           fVoxelDepth++;
604           newHeader = newProxy->GetHeader();      609           newHeader = newProxy->GetHeader();
605           newHeaderAxis = newHeader->GetAxis()    610           newHeaderAxis = newHeader->GetAxis();
606           newHeaderNoSlices = (G4int)newHeader << 611           newHeaderNoSlices = newHeader->GetNoSlices();
607           newHeaderMin = newHeader->GetMinExte    612           newHeaderMin = newHeader->GetMinExtent();
608           newHeaderNodeWidth = (newHeader->Get    613           newHeaderNodeWidth = (newHeader->GetMaxExtent()-newHeaderMin)
609                              / newHeaderNoSlic    614                              / newHeaderNoSlices;
610           newNodeNo = G4int( (voxelPoint(newHe    615           newNodeNo = G4int( (voxelPoint(newHeaderAxis)-newHeaderMin)
611                              / newHeaderNodeWi    616                              / newHeaderNodeWidth );
612           // Rounding protection                  617           // Rounding protection
613           //                                      618           //
614           if ( newNodeNo<0 )                      619           if ( newNodeNo<0 )
615           {                                       620           {
616             newNodeNo=0;                          621             newNodeNo=0;
617           }                                       622           }
618           else if ( newNodeNo>=newHeaderNoSlic    623           else if ( newNodeNo>=newHeaderNoSlices )
619           {                                    << 624                {
620             newNodeNo = newHeaderNoSlices-1;   << 625                  newNodeNo = newHeaderNoSlices-1;
621           }                                    << 626                }
622           // Stack info for stepping              627           // Stack info for stepping
623           //                                      628           //
624           fVoxelAxisStack[fVoxelDepth] = newHe    629           fVoxelAxisStack[fVoxelDepth] = newHeaderAxis;
625           fVoxelNoSlicesStack[fVoxelDepth] = n    630           fVoxelNoSlicesStack[fVoxelDepth] = newHeaderNoSlices;
626           fVoxelSliceWidthStack[fVoxelDepth] =    631           fVoxelSliceWidthStack[fVoxelDepth] = newHeaderNodeWidth;
627           fVoxelNodeNoStack[fVoxelDepth] = new    632           fVoxelNodeNoStack[fVoxelDepth] = newNodeNo;
628           fVoxelHeaderStack[fVoxelDepth] = new    633           fVoxelHeaderStack[fVoxelDepth] = newHeader;
629         }                                         634         }
630       }                                           635       }
631       fVoxelNode = newVoxelNode;                  636       fVoxelNode = newVoxelNode;
632     }                                             637     }
633   }                                               638   }
634   return isNewVoxel;                              639   return isNewVoxel;        
635 }                                                 640 }
636                                                   641 
637 // *******************************************    642 // ********************************************************************
638 // ComputeSafety                                  643 // ComputeSafety
639 //                                                644 //
640 // Calculates the isotropic distance to the ne    645 // Calculates the isotropic distance to the nearest boundary from the
641 // specified point in the local coordinate sys    646 // specified point in the local coordinate system. 
642 // The localpoint utilised must be within the     647 // The localpoint utilised must be within the current volume.
643 // *******************************************    648 // ********************************************************************
644 //                                                649 //
645 G4double                                          650 G4double
646 G4VoxelNavigation::ComputeSafety(const G4Three    651 G4VoxelNavigation::ComputeSafety(const G4ThreeVector& localPoint,
647                                  const G4Navig    652                                  const G4NavigationHistory& history,
648                                  const G4doubl << 653                                  const G4double )
649 {                                                 654 {
650   G4VPhysicalVolume *motherPhysical, *samplePh    655   G4VPhysicalVolume *motherPhysical, *samplePhysical;
651   G4LogicalVolume *motherLogical;                 656   G4LogicalVolume *motherLogical;
652   G4VSolid *motherSolid;                          657   G4VSolid *motherSolid;
653   G4double motherSafety, ourSafety;               658   G4double motherSafety, ourSafety;
654   G4int sampleNo;                              << 659   G4int localNoDaughters, sampleNo;
655   G4SmartVoxelNode *curVoxelNode;                 660   G4SmartVoxelNode *curVoxelNode;
656   G4long curNoVolumes, contentNo;              << 661   G4int curNoVolumes, contentNo;
657   G4double voxelSafety;                           662   G4double voxelSafety;
658                                                   663 
659   motherPhysical = history.GetTopVolume();        664   motherPhysical = history.GetTopVolume();
660   motherLogical = motherPhysical->GetLogicalVo    665   motherLogical = motherPhysical->GetLogicalVolume();
661   motherSolid = motherLogical->GetSolid();        666   motherSolid = motherLogical->GetSolid();
662                                                   667 
663   if( fBestSafety )                            << 
664   {                                            << 
665     return fpVoxelSafety->ComputeSafety( local << 
666   }                                            << 
667                                                << 
668   //                                              668   //
669   // Compute mother safety                        669   // Compute mother safety
670   //                                              670   //
671                                                   671 
672   motherSafety = motherSolid->DistanceToOut(lo    672   motherSafety = motherSolid->DistanceToOut(localPoint);
673   ourSafety = motherSafety;                 //    673   ourSafety = motherSafety;                 // Working isotropic safety
674                                                   674 
675   if( motherSafety == 0.0 )                    << 
676   {                                            << 
677 #ifdef G4DEBUG_NAVIGATION                      << 
678     // Check that point is inside mother volum << 
679     EInside  insideMother = motherSolid->Insid << 
680                                                << 
681     if( insideMother == kOutside )             << 
682     {                                          << 
683       G4ExceptionDescription message;          << 
684       message << "Safety method called for loc << 
685          << "Location for safety is Outside th << 
686          << "The approximate distance to the s << 
687          << "(safety from outside) is: "       << 
688          << motherSolid->DistanceToIn( localPo << 
689       message << "  Problem occurred with phys << 
690          << " Name: " << motherPhysical->GetNa << 
691          << " Copy No: " << motherPhysical->Ge << 
692          << "    Local Point = " << localPoint << 
693       message << "  Description of solid: " << << 
694             << *motherSolid << G4endl;         << 
695       G4Exception("G4VoxelNavigation::ComputeS << 
696                   JustWarning, message);       << 
697     }                                          << 
698                                                << 
699     // Following check is NOT for an issue - i << 
700     //  It is allowed that a solid gives appro << 
701     //                                         << 
702     if( insideMother == kInside ) // && fVerbo << 
703     {                                          << 
704       G4ExceptionDescription messageIn;        << 
705                                                << 
706       messageIn << " Point is Inside, but safe << 
707       messageIn << " Inexact safety for volume << 
708              << "  Solid: Name= " << motherSol << 
709              << "   Type= " << motherSolid->Ge << 
710       messageIn << "  Local point= " << localP << 
711       messageIn << "  Solid parameters: " << G << 
712       G4Exception("G4VoxelNavigation::ComputeS << 
713                   JustWarning, messageIn);     << 
714     }                                          << 
715 #endif                                         << 
716     // if( insideMother != kInside )           << 
717     return 0.0;                                << 
718   }                                            << 
719                                                << 
720 #ifdef G4VERBOSE                                  675 #ifdef G4VERBOSE
721   if( fCheck )                                 << 676   if(( fCheck ) && ( fVerbose == 1 ))
722   {                                               677   {
723     fLogger->ComputeSafetyLog (motherSolid,loc << 678     G4cout << "*** G4VoxelNavigation::ComputeSafety(): ***" << G4endl
                                                   >> 679            << "    Invoked DistanceToOut(p) for mother solid: "
                                                   >> 680            << motherSolid->GetName()
                                                   >> 681            << ". Solid replied: " << motherSafety << G4endl
                                                   >> 682            << "    For local point p: " << localPoint
                                                   >> 683            << ", to be considered as 'mother safety'." << G4endl;
724   }                                               684   }
725 #endif                                            685 #endif
726   //                                              686   //
727   // Compute daughter safeties                 << 687   // Compute daughter safeties 
728   //                                              688   //
729   // Look only inside the current Voxel only ( << 689 
                                                   >> 690   localNoDaughters = motherLogical->GetNoDaughters();
                                                   >> 691 
                                                   >> 692   //  Look only inside the current Voxel only (in the first version).
730   //                                              693   //
731   curVoxelNode = fVoxelNode;                      694   curVoxelNode = fVoxelNode;
732   curNoVolumes = curVoxelNode->GetNoContained(    695   curNoVolumes = curVoxelNode->GetNoContained();
733                                                   696 
734   for ( contentNo=curNoVolumes-1; contentNo>=0    697   for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
735   {                                               698   {
736     sampleNo = curVoxelNode->GetVolume((G4int) << 699     sampleNo = curVoxelNode->GetVolume(contentNo);
737     samplePhysical = motherLogical->GetDaughte    700     samplePhysical = motherLogical->GetDaughter(sampleNo);
738                                                   701 
739     G4AffineTransform sampleTf(samplePhysical-    702     G4AffineTransform sampleTf(samplePhysical->GetRotation(),
740                                samplePhysical-    703                                samplePhysical->GetTranslation());
741     sampleTf.Invert();                            704     sampleTf.Invert();
742     const G4ThreeVector samplePoint = sampleTf << 705     const G4ThreeVector samplePoint =
743     const G4VSolid* sampleSolid= samplePhysica << 706                           sampleTf.TransformPoint(localPoint);
                                                   >> 707     const G4VSolid *sampleSolid     =
                                                   >> 708                           samplePhysical->GetLogicalVolume()->GetSolid();
744     G4double sampleSafety = sampleSolid->Dista    709     G4double sampleSafety = sampleSolid->DistanceToIn(samplePoint);
745     if ( sampleSafety<ourSafety )                 710     if ( sampleSafety<ourSafety )
746     {                                             711     {
747       ourSafety = sampleSafety;                   712       ourSafety = sampleSafety;
748     }                                             713     }
749 #ifdef G4VERBOSE                                  714 #ifdef G4VERBOSE
750     if( fCheck )                               << 715     if(( fCheck ) && ( fVerbose == 1 ))
751     {                                             716     {
752       fLogger->ComputeSafetyLog(sampleSolid, s << 717       G4cout << "*** G4VoxelNavigation::ComputeSafety(): ***" << G4endl
753                                 sampleSafety,  << 718              << "    Invoked DistanceToIn(p) for daughter solid: "
                                                   >> 719              << sampleSolid->GetName()
                                                   >> 720              << ". Solid replied: " << sampleSafety << G4endl
                                                   >> 721              << "    For local point p: " << samplePoint
                                                   >> 722              << ", to be considered as 'daughter safety'." << G4endl;
754     }                                             723     }
755 #endif                                            724 #endif
756   }                                               725   }
757   voxelSafety = ComputeVoxelSafety(localPoint)    726   voxelSafety = ComputeVoxelSafety(localPoint);
758   if ( voxelSafety<ourSafety )                    727   if ( voxelSafety<ourSafety )
759   {                                               728   {
760     ourSafety = voxelSafety;                      729     ourSafety = voxelSafety;
761   }                                               730   }
762   return ourSafety;                               731   return ourSafety;
763 }                                              << 
764                                                << 
765 void G4VoxelNavigation::RelocateWithinVolume(  << 
766                                                << 
767 {                                              << 
768   auto motherLogical = motherPhysical->GetLogi << 
769                                                << 
770   assert(motherLogical != nullptr);            << 
771                                                << 
772   if ( auto pVoxelHeader = motherLogical->GetV << 
773     VoxelLocate( pVoxelHeader, localPoint );   << 
774 }                                              << 
775                                                << 
776 // ******************************************* << 
777 // SetVerboseLevel                             << 
778 // ******************************************* << 
779 //                                             << 
780 void  G4VoxelNavigation::SetVerboseLevel(G4int << 
781 {                                              << 
782   if( fLogger != nullptr ) { fLogger->SetVerbo << 
783   if( fpVoxelSafety != nullptr) { fpVoxelSafet << 
784 }                                                 732 }
785                                                   733