Geant4 Cross Reference

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


  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 //
                                                   >>  27 // $Id: G4RegularNavigation.cc,v 1.9 2009/01/27 09:31:29 gcosmo Exp $
                                                   >>  28 // GEANT4 tag $ Name:$
                                                   >>  29 //
 26 // class G4RegularNavigation implementation        30 // class G4RegularNavigation implementation
 27 //                                                 31 //
 28 // Author: Pedro Arce, May 2007                    32 // Author: Pedro Arce, May 2007
 29 //                                                 33 //
 30 // -------------------------------------------     34 // --------------------------------------------------------------------
 31                                                    35 
 32 #include "G4RegularNavigation.hh"                  36 #include "G4RegularNavigation.hh"
 33 #include "G4TouchableHistory.hh"                   37 #include "G4TouchableHistory.hh"
 34 #include "G4PhantomParameterisation.hh"            38 #include "G4PhantomParameterisation.hh"
 35 #include "G4Material.hh"                           39 #include "G4Material.hh"
 36 #include "G4NormalNavigation.hh"                   40 #include "G4NormalNavigation.hh"
 37 #include "G4Navigator.hh"                          41 #include "G4Navigator.hh"
 38 #include "G4GeometryTolerance.hh"                  42 #include "G4GeometryTolerance.hh"
 39 #include "G4RegularNavigationHelper.hh"            43 #include "G4RegularNavigationHelper.hh"
 40                                                    44 
 41 //--------------------------------------------     45 //------------------------------------------------------------------
 42 G4RegularNavigation::G4RegularNavigation()         46 G4RegularNavigation::G4RegularNavigation()
                                                   >>  47   : fVerbose(1), fCheck(true)
 43 {                                                  48 {
 44   kCarTolerance = G4GeometryTolerance::GetInst     49   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
 45   fMinStep = 101*kCarTolerance;                << 
 46 }                                                  50 }
 47                                                    51 
 48                                                    52 
 49 //--------------------------------------------     53 //------------------------------------------------------------------
 50 G4RegularNavigation::~G4RegularNavigation() =  <<  54 G4RegularNavigation::~G4RegularNavigation()
                                                   >>  55 {
                                                   >>  56 }
 51                                                    57 
 52                                                    58 
 53 //--------------------------------------------     59 //------------------------------------------------------------------
 54 G4double G4RegularNavigation::                     60 G4double G4RegularNavigation::
 55                     ComputeStep(const G4ThreeV     61                     ComputeStep(const G4ThreeVector& localPoint,
 56                                 const G4ThreeV     62                                 const G4ThreeVector& localDirection,
 57                                 const G4double     63                                 const G4double currentProposedStepLength,
 58                                       G4double     64                                       G4double& newSafety,
 59                                       G4Naviga     65                                       G4NavigationHistory& history,
 60                                       G4bool&      66                                       G4bool& validExitNormal,
 61                                       G4ThreeV     67                                       G4ThreeVector& exitNormal,
 62                                       G4bool&      68                                       G4bool& exiting,
 63                                       G4bool&      69                                       G4bool& entering,
 64                                       G4VPhysi     70                                       G4VPhysicalVolume *(*pBlockedPhysical),
 65                                       G4int& b     71                                       G4int& blockedReplicaNo)
 66 {                                                  72 {
 67   // This method is never called because to be     73   // This method is never called because to be called the daughter has to be
 68   // a regular structure. This would only happ     74   // a regular structure. This would only happen if the track is in the mother
 69   // of voxels volume. But the voxels fill com     75   // of voxels volume. But the voxels fill completely their mother, so when a
 70   // track enters the mother it automatically      76   // track enters the mother it automatically enters a voxel. Only precision
 71   // problems would make this method to be cal     77   // problems would make this method to be called
 72                                                    78 
 73   G4ThreeVector globalPoint =                      79   G4ThreeVector globalPoint =
 74     history.GetTopTransform().InverseTransform <<  80          history.GetTopTransform().Inverse().TransformPoint(localPoint);
 75   G4ThreeVector globalDirection =                  81   G4ThreeVector globalDirection =
 76          history.GetTopTransform().InverseTran <<  82          history.GetTopTransform().Inverse().TransformAxis(localDirection);
 77                                                    83 
 78   G4ThreeVector localPoint2 = localPoint; // t     84   G4ThreeVector localPoint2 = localPoint; // take away constantness
 79                                                    85 
 80   LevelLocate( history, *pBlockedPhysical, blo     86   LevelLocate( history, *pBlockedPhysical, blockedReplicaNo, 
 81                globalPoint, &globalDirection,      87                globalPoint, &globalDirection, true, localPoint2 );
 82                                                    88 
 83                                                    89 
 84   // Get in which voxel it is                      90   // Get in which voxel it is
 85   //                                               91   //
 86   G4VPhysicalVolume *motherPhysical, *daughter     92   G4VPhysicalVolume *motherPhysical, *daughterPhysical;
 87   G4LogicalVolume *motherLogical;                  93   G4LogicalVolume *motherLogical;
 88   motherPhysical = history.GetTopVolume();         94   motherPhysical = history.GetTopVolume();
 89   motherLogical = motherPhysical->GetLogicalVo     95   motherLogical = motherPhysical->GetLogicalVolume();
 90   daughterPhysical = motherLogical->GetDaughte     96   daughterPhysical = motherLogical->GetDaughter(0);
 91                                                    97 
 92   auto daughterParam =                         <<  98   G4PhantomParameterisation * daughterParam =
 93        (G4PhantomParameterisation*)(daughterPh <<  99         (G4PhantomParameterisation*)(daughterPhysical->GetParameterisation());
 94   G4int copyNo = daughterParam ->GetReplicaNo(    100   G4int copyNo = daughterParam ->GetReplicaNo(localPoint,localDirection);
 95                                                   101 
 96   G4ThreeVector voxelTranslation = daughterPar    102   G4ThreeVector voxelTranslation = daughterParam->GetTranslation( copyNo );
 97   G4ThreeVector daughterPoint = localPoint - v    103   G4ThreeVector daughterPoint = localPoint - voxelTranslation;
 98                                                   104 
 99                                                   105 
100   // Compute step in voxel                        106   // Compute step in voxel
101   //                                              107   //
102   return fnormalNav->ComputeStep(daughterPoint    108   return fnormalNav->ComputeStep(daughterPoint,
103                                  localDirectio    109                                  localDirection,
104                                  currentPropos    110                                  currentProposedStepLength,
105                                  newSafety,       111                                  newSafety,
106                                  history,         112                                  history,
107                                  validExitNorm    113                                  validExitNormal,
108                                  exitNormal,      114                                  exitNormal,
109                                  exiting,         115                                  exiting,
110                                  entering,        116                                  entering,
111                                  pBlockedPhysi    117                                  pBlockedPhysical,
112                                  blockedReplic    118                                  blockedReplicaNo);
113 }                                                 119 }
114                                                   120 
115                                                   121 
116 //--------------------------------------------    122 //------------------------------------------------------------------
117 G4double G4RegularNavigation::ComputeStepSkipp    123 G4double G4RegularNavigation::ComputeStepSkippingEqualMaterials(
118                                       G4ThreeV << 124                                 G4ThreeVector localPoint,
119                                 const G4ThreeV    125                                 const G4ThreeVector& localDirection,
120                                 const G4double    126                                 const G4double currentProposedStepLength,
121                                 G4double& newS    127                                 G4double& newSafety,
122                                 G4NavigationHi    128                                 G4NavigationHistory& history,
123                                 G4bool& validE    129                                 G4bool& validExitNormal,
124                                 G4ThreeVector&    130                                 G4ThreeVector& exitNormal,
125                                 G4bool& exitin    131                                 G4bool& exiting,
126                                 G4bool& enteri    132                                 G4bool& entering,
127                                 G4VPhysicalVol    133                                 G4VPhysicalVolume *(*pBlockedPhysical),
128                                 G4int& blocked    134                                 G4int& blockedReplicaNo,
129                                 G4VPhysicalVol    135                                 G4VPhysicalVolume* pCurrentPhysical)
130 {                                                 136 {
131   G4RegularNavigationHelper::Instance()->Clear << 137   G4RegularNavigationHelper::ClearStepLengths();
132                                                   138 
133   auto param =                                 << 139   G4PhantomParameterisation *param =
134     (G4PhantomParameterisation*)(pCurrentPhysi    140     (G4PhantomParameterisation*)(pCurrentPhysical->GetParameterisation());
135                                                   141 
136   if( !param->SkipEqualMaterials() )              142   if( !param->SkipEqualMaterials() )
137   {                                               143   {
138     return fnormalNav->ComputeStep(localPoint,    144     return fnormalNav->ComputeStep(localPoint,
139                                    localDirect    145                                    localDirection,
140                                    currentProp    146                                    currentProposedStepLength,
141                                    newSafety,     147                                    newSafety,
142                                    history,       148                                    history,
143                                    validExitNo    149                                    validExitNormal,
144                                    exitNormal,    150                                    exitNormal,
145                                    exiting,       151                                    exiting,
146                                    entering,      152                                    entering,
147                                    pBlockedPhy    153                                    pBlockedPhysical,
148                                    blockedRepl    154                                    blockedReplicaNo);
149   }                                               155   }
150                                                   156 
151                                                   157 
152   G4double ourStep = 0.;                          158   G4double ourStep = 0.;
153                                                   159 
154   // To get replica No: transform local point     160   // To get replica No: transform local point to the reference system of the
155   // param container volume                       161   // param container volume
156   //                                              162   //
157   auto  ide = (G4int)history.GetDepth();       << 163   G4int ide = history.GetDepth();
158   G4ThreeVector containerPoint = history.GetTr << 164   G4ThreeVector containerPoint = history.GetTransform(ide).Inverse().TransformPoint(localPoint);
159                                  .InverseTrans << 
160                                                   165 
161   // Point in global frame                        166   // Point in global frame
162   //                                              167   //
163   containerPoint = history.GetTransform(ide).I << 168   containerPoint = history.GetTransform(ide).Inverse().TransformPoint(localPoint);
164                                                   169 
165   // Point in voxel parent volume frame           170   // Point in voxel parent volume frame
166   //                                              171   //
167   containerPoint = history.GetTransform(ide-1)    172   containerPoint = history.GetTransform(ide-1).TransformPoint(containerPoint);
168                                                   173 
169   // Store previous voxel translation to move     174   // Store previous voxel translation to move localPoint by the difference
170   // with the new one                             175   // with the new one
171   //                                              176   //
172   G4ThreeVector prevVoxelTranslation = contain    177   G4ThreeVector prevVoxelTranslation = containerPoint - localPoint;
173                                                   178 
174   // Do not use the expression below: There ar    179   // Do not use the expression below: There are cases where the
175   // fLastLocatedPointLocal does not give the     180   // fLastLocatedPointLocal does not give the correct answer
176   // (particle reaching a wall and bounced bac    181   // (particle reaching a wall and bounced back, particle travelling through
177   // the wall that is deviated in an step, ...    182   // the wall that is deviated in an step, ...; these are pathological cases
178   // that give wrong answers in G4PhantomParam    183   // that give wrong answers in G4PhantomParameterisation::GetReplicaNo()
179   //                                              184   //
180   // G4ThreeVector prevVoxelTranslation = para    185   // G4ThreeVector prevVoxelTranslation = param->GetTranslation( copyNo );
181                                                   186 
182   G4int copyNo = param->GetReplicaNo(container    187   G4int copyNo = param->GetReplicaNo(containerPoint,localDirection);
183                                                   188 
184   G4Material* currentMate = param->ComputeMate << 189   G4Material* currentMate = param->ComputeMaterial( copyNo, 0, 0 );
185   G4VSolid* voxelBox = pCurrentPhysical->GetLo    190   G4VSolid* voxelBox = pCurrentPhysical->GetLogicalVolume()->GetSolid();
186                                                   191 
187   G4VSolid* containerSolid = param->GetContain    192   G4VSolid* containerSolid = param->GetContainerSolid();
188   G4Material* nextMate;                           193   G4Material* nextMate;
                                                   >> 194   G4bool bLocatedOnEdge = false;
189   G4bool bFirstStep = true;                       195   G4bool bFirstStep = true;
190   G4double newStep;                               196   G4double newStep;
191   G4double totalNewStep = 0.;                     197   G4double totalNewStep = 0.;
192                                                   198 
193   // Loop while same material is found            199   // Loop while same material is found 
194   //                                              200   //
195   //                                           << 201   for( ;; )
196   fNumberZeroSteps = 0;                        << 
197   for( G4int ii = 0; ii < fNoStepsAllowed+1; + << 
198   {                                               202   {
199     if( ii == fNoStepsAllowed ) {              << 
200       // Must kill this stuck track            << 
201       //                                       << 
202       G4ThreeVector pGlobalpoint = history.Get << 
203                                    .InverseTra << 
204       std::ostringstream message;              << 
205       message << "G4RegularNavigation::Compute << 
206         << "Stuck Track: potential geometry or << 
207         << G4endl                              << 
208         << "        Track stuck, moving for mo << 
209         << ii << " steps" << G4endl            << 
210         << "- at point " << pGlobalpoint << G4 << 
211         << "        local direction: " << loca << 
212       G4Exception("G4RegularNavigation::Comput << 
213       "GeomRegNav1001",                        << 
214       EventMustBeAborted,                      << 
215       message);                                << 
216     }                                          << 
217     newStep = voxelBox->DistanceToOut( localPo    203     newStep = voxelBox->DistanceToOut( localPoint, localDirection );
218     fLastStepWasZero = (newStep<fMinStep);     << 204 
219     if( fLastStepWasZero )                     << 
220     {                                          << 
221       ++fNumberZeroSteps;                      << 
222 #ifdef G4DEBUG_NAVIGATION                      << 
223       if( fNumberZeroSteps > 1 )               << 
224       {                                        << 
225         G4ThreeVector pGlobalpoint = history.G << 
226                                      .InverseT << 
227       std::ostringstream message;              << 
228       message.precision(16);                   << 
229       message << "G4RegularNavigation::Compute << 
230             << fNumberZeroSteps                << 
231             << ", at " << pGlobalpoint         << 
232             << ", nav-comp-step calls # " << i << 
233             << ", Step= " << newStep;          << 
234             G4Exception("G4RegularNavigation:: << 
235                         "GeomRegNav1002", Just << 
236                         "Potential overlap in  << 
237       }                                        << 
238 #endif                                         << 
239       if( fNumberZeroSteps > fActionThreshold_ << 
240       {                                        << 
241         // Act to recover this stuck track. Pu << 
242         //                                     << 
243         newStep = std::min(101*kCarTolerance*s << 
244 #ifdef G4DEBUG_NAVIGATION                      << 
245         G4ThreeVector pGlobalpoint = history.G << 
246                                        .Invers << 
247         std::ostringstream message;            << 
248         message.precision(16);                 << 
249         message << "Track stuck or not moving. << 
250                       << "          Track stuc << 
251                       << fNumberZeroSteps << " << 
252                       << "- at point " << pGlo << 
253                       << " (local point " << l << 
254                       << "        local direct << 
255                       << "          Potential  << 
256                       << G4endl                << 
257                       << "          Trying pus << 
258               G4Exception("G4RegularNavigation << 
259                           "GeomRegNav1003", Ju << 
260                           "Potential overlap i << 
261 #endif                                         << 
262       }                                        << 
263       if( fNumberZeroSteps > fAbandonThreshold << 
264       {                                        << 
265         // Must kill this stuck track          << 
266         //                                     << 
267         G4ThreeVector pGlobalpoint = history.G << 
268                                           .Inv << 
269         std::ostringstream message;            << 
270         message << "G4RegularNavigation::Compu << 
271           << "Stuck Track: potential geometry  << 
272           << G4endl                            << 
273           << "        Track stuck, not moving  << 
274           << fNumberZeroSteps << " steps" << G << 
275           << "- at point " << pGlobalpoint <<  << 
276           << "        local direction: " << lo << 
277         G4Exception("G4RegularNavigation::Comp << 
278               "GeomRegNav1004",                << 
279               EventMustBeAborted,              << 
280               message);                        << 
281       }                                        << 
282     }                                          << 
283     else                                       << 
284     {                                          << 
285       // reset the zero step counter when a no << 
286       fNumberZeroSteps = 0;                    << 
287     }                                          << 
288     if( (bFirstStep) && (newStep < currentProp    205     if( (bFirstStep) && (newStep < currentProposedStepLength) )
289     {                                             206     {
290       exiting  = true;                            207       exiting  = true;
291     }                                             208     }
292     bFirstStep = false;                           209     bFirstStep = false;
293                                                   210  
294     newStep += kCarTolerance;   // Avoid preci    211     newStep += kCarTolerance;   // Avoid precision problems
295     ourStep += newStep;                           212     ourStep += newStep;
296     totalNewStep += newStep;                      213     totalNewStep += newStep;
297                                                   214 
298     // Physical process is limiting the step,     215     // Physical process is limiting the step, don't continue
299     //                                            216     //
300     if(std::fabs(totalNewStep-currentProposedS    217     if(std::fabs(totalNewStep-currentProposedStepLength) < kCarTolerance)
301     {                                             218     { 
302       return currentProposedStepLength;           219       return currentProposedStepLength;
303     }                                             220     }
304     if(totalNewStep > currentProposedStepLengt    221     if(totalNewStep > currentProposedStepLength) 
305     {                                             222     { 
306       G4RegularNavigationHelper::Instance()->  << 223       G4RegularNavigationHelper::
307         AddStepLength(copyNo, newStep-totalNew    224         AddStepLength(copyNo, newStep-totalNewStep+currentProposedStepLength);
308       return currentProposedStepLength;           225       return currentProposedStepLength;
309     }                                             226     }
310     G4RegularNavigationHelper::Instance()->Add << 227     else
                                                   >> 228     {
                                                   >> 229       G4RegularNavigationHelper::AddStepLength( copyNo, newStep );
                                                   >> 230     }
                                                   >> 231 
311                                                   232 
312     // Move container point until wall of voxe    233     // Move container point until wall of voxel
313     //                                            234     //
314     containerPoint += newStep*localDirection;     235     containerPoint += newStep*localDirection;
315     if( containerSolid->Inside( containerPoint    236     if( containerSolid->Inside( containerPoint ) != kInside )
316     {                                             237     {
                                                   >> 238       bLocatedOnEdge = true; // NEEDED??
317       break;                                      239       break;
318     }                                             240     }
319                                                   241 
320     // Get copyNo and translation of new voxel    242     // Get copyNo and translation of new voxel
321     //                                            243     //
322     copyNo = param->GetReplicaNo(containerPoin << 244     copyNo = param->GetReplicaNo(containerPoint,localDirection);
323     G4ThreeVector voxelTranslation = param->Ge    245     G4ThreeVector voxelTranslation = param->GetTranslation( copyNo );
324                                                   246 
                                                   >> 247     //    G4cout << " copyNo " << copyNo << " = " << pCurrentPhysical->GetCopyNo() << G4endl;
325     // Move local point until wall of voxel an    248     // Move local point until wall of voxel and then put it in the new voxel
326     // local coordinates                          249     // local coordinates
327     //                                            250     //
328     localPoint += newStep*localDirection;         251     localPoint += newStep*localDirection;
329     localPoint += prevVoxelTranslation - voxel    252     localPoint += prevVoxelTranslation - voxelTranslation;
330                                                   253 
331     prevVoxelTranslation = voxelTranslation;      254     prevVoxelTranslation = voxelTranslation;
332                                                   255 
333     // Check if material of next voxel is the     256     // Check if material of next voxel is the same as that of the current voxel
334     nextMate = param->ComputeMaterial( copyNo, << 257     nextMate = param->ComputeMaterial( copyNo, 0, 0 );
335                                                   258 
336     if( currentMate != nextMate ) { break; }      259     if( currentMate != nextMate ) { break; }
337   }                                               260   }
338                                                   261 
339   return ourStep;                                 262   return ourStep;
340 }                                                 263 }
341                                                   264 
342                                                   265 
343 //--------------------------------------------    266 //------------------------------------------------------------------
344 G4double                                          267 G4double
345 G4RegularNavigation::ComputeSafety(const G4Thr    268 G4RegularNavigation::ComputeSafety(const G4ThreeVector& localPoint,
346                                    const G4Nav    269                                    const G4NavigationHistory& history,
347                                    const G4dou    270                                    const G4double pMaxLength)
348 {                                                 271 {
349   // This method is never called because to be    272   // This method is never called because to be called the daughter has to be a
350   // regular structure. This would only happen    273   // regular structure. This would only happen if the track is in the mother of
351   // voxels volume. But the voxels fill comple    274   // voxels volume. But the voxels fill completely their mother, so when a
352   // track enters the mother it automatically     275   // track enters the mother it automatically enters a voxel. Only precision
353   // problems would make this method to be cal    276   // problems would make this method to be called
354                                                   277 
355   // Compute step in voxel                        278   // Compute step in voxel
356   //                                              279   //
357   return fnormalNav->ComputeSafety(localPoint,    280   return fnormalNav->ComputeSafety(localPoint,
358                                    history,       281                                    history,
359                                    pMaxLength     282                                    pMaxLength );
360 }                                                 283 }
361                                                   284 
362                                                   285 
363 //--------------------------------------------    286 //------------------------------------------------------------------
364 G4bool                                            287 G4bool
365 G4RegularNavigation::LevelLocate( G4Navigation    288 G4RegularNavigation::LevelLocate( G4NavigationHistory& history,
366                                   const G4VPhy    289                                   const G4VPhysicalVolume* ,
367                                   const G4int     290                                   const G4int ,
368                                   const G4Thre    291                                   const G4ThreeVector& globalPoint,
369                                   const G4Thre    292                                   const G4ThreeVector* globalDirection,
370                                   const G4bool << 293                                   const G4bool pLocatedOnEdge, 
371                                   G4ThreeVecto    294                                   G4ThreeVector& localPoint )
372 {                                                 295 {
                                                   >> 296   G4SmartVoxelHeader *motherVoxelHeader;
373   G4VPhysicalVolume *motherPhysical, *pPhysica    297   G4VPhysicalVolume *motherPhysical, *pPhysical;
374   G4PhantomParameterisation *pParam;              298   G4PhantomParameterisation *pParam;
375   G4LogicalVolume *motherLogical;                 299   G4LogicalVolume *motherLogical;
                                                   >> 300   G4VSolid *pSolid;
376   G4ThreeVector localDir;                         301   G4ThreeVector localDir;
377   G4int replicaNo;                                302   G4int replicaNo;
378                                                   303   
379   motherPhysical = history.GetTopVolume();        304   motherPhysical = history.GetTopVolume();
380   motherLogical = motherPhysical->GetLogicalVo    305   motherLogical = motherPhysical->GetLogicalVolume();
                                                   >> 306   motherVoxelHeader = motherLogical->GetVoxelHeader();
381                                                   307   
382   pPhysical = motherLogical->GetDaughter(0);      308   pPhysical = motherLogical->GetDaughter(0);
383   pParam = (G4PhantomParameterisation*)(pPhysi    309   pParam = (G4PhantomParameterisation*)(pPhysical->GetParameterisation());
384                                                   310   
                                                   >> 311   pSolid = pParam->GetContainerSolid();
                                                   >> 312   
385   // Save parent history in touchable history     313   // Save parent history in touchable history
386   // ... for use as parent t-h in ComputeMater    314   // ... for use as parent t-h in ComputeMaterial method of param
387   //                                              315   //
388   G4TouchableHistory parentTouchable( history     316   G4TouchableHistory parentTouchable( history ); 
389                                                   317   
390   // Get local direction                          318   // Get local direction
391   //                                              319   //
392   if( globalDirection != nullptr )             << 320   if( globalDirection )
393   {                                               321   {
394     localDir = history.GetTopTransform().Trans    322     localDir = history.GetTopTransform().TransformAxis(*globalDirection);
395   }                                               323   }
396   else                                            324   else
397   {                                               325   {
398     localDir = G4ThreeVector(0.,0.,0.);           326     localDir = G4ThreeVector(0.,0.,0.);
399   }                                               327   }
400                                                   328 
401   // Enter this daughter                       << 329   // Check that track is not on the surface and check that track is not
                                                   >> 330   // exiting the voxel parent volume
402   //                                              331   //
403   replicaNo = pParam->GetReplicaNo( localPoint << 332   if ( !G4AuxiliaryNavServices::CheckPointOnSurface(pSolid, localPoint,
404                                                << 333           globalDirection, history.GetTopTransform(), pLocatedOnEdge) 
405   if( replicaNo < 0 || replicaNo >= G4int(pPar << 334      || G4AuxiliaryNavServices::CheckPointExiting(pSolid, localPoint,
                                                   >> 335           globalDirection, history.GetTopTransform() ) )
406   {                                               336   {
407     return false;                              << 
408   }                                               337   }
                                                   >> 338   else
                                                   >> 339   {
                                                   >> 340     // Enter this daughter
                                                   >> 341     //
                                                   >> 342     replicaNo = pParam->GetReplicaNo( localPoint, localDir );
409                                                   343 
410   // Set the correct copy number in physical   << 344     if( replicaNo < 0 || replicaNo >= G4int(pParam->GetNoVoxel()) )
411   //                                           << 345     {
412   pPhysical->SetCopyNo(replicaNo);             << 346       return false;
413   pParam->ComputeTransformation(replicaNo,pPhy << 347     }
                                                   >> 348 
                                                   >> 349     // Set the correct copy number in physical
                                                   >> 350     //
                                                   >> 351     pPhysical->SetCopyNo(replicaNo);
                                                   >> 352     pParam->ComputeTransformation(replicaNo,pPhysical);
414                                                   353 
415   history.NewLevel(pPhysical, kParameterised,  << 354     history.NewLevel(pPhysical, kParameterised, replicaNo );
416   localPoint = history.GetTopTransform().Trans << 355     localPoint = history.GetTopTransform().TransformPoint(globalPoint);
417                                                   356 
418   // Set the correct solid and material in Log << 357     // Set the correct solid and material in Logical Volume
419   //                                           << 358     //
420   G4LogicalVolume *pLogical = pPhysical->GetLo << 359     G4LogicalVolume *pLogical = pPhysical->GetLogicalVolume();
421                                                   360       
422   pLogical->UpdateMaterial(pParam->ComputeMate << 361     pLogical->UpdateMaterial(pParam->ComputeMaterial(replicaNo,
423                            pPhysical, &parentT << 362                              pPhysical, &parentTouchable) );
424   return true;                                 << 363     return true;
                                                   >> 364   }
                                                   >> 365 
                                                   >> 366   return false;
425 }                                                 367 }
426                                                   368