Geant4 Cross Reference

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


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 //
                                                   >>  27 // $Id: G4PhantomParameterisation.cc 66356 2012-12-18 09:02:32Z gcosmo $
                                                   >>  28 // GEANT4 tag $ Name:$
                                                   >>  29 //
 26 // class G4PhantomParameterisation implementat     30 // class G4PhantomParameterisation implementation
 27 //                                                 31 //
 28 // May 2007 Pedro Arce, first version          <<  32 // May 2007 Pedro Arce,   first version
 29 //                                                 33 //
 30 // -------------------------------------------     34 // --------------------------------------------------------------------
 31                                                    35 
 32 #include "G4PhantomParameterisation.hh"            36 #include "G4PhantomParameterisation.hh"
 33                                                    37 
 34 #include "globals.hh"                              38 #include "globals.hh"
 35 #include "G4VSolid.hh"                             39 #include "G4VSolid.hh"
 36 #include "G4VPhysicalVolume.hh"                    40 #include "G4VPhysicalVolume.hh"
 37 #include "G4LogicalVolume.hh"                      41 #include "G4LogicalVolume.hh"
 38 #include "G4VVolumeMaterialScanner.hh"             42 #include "G4VVolumeMaterialScanner.hh"
 39 #include "G4GeometryTolerance.hh"                  43 #include "G4GeometryTolerance.hh"
 40                                                    44 
 41 //--------------------------------------------     45 //------------------------------------------------------------------
 42 G4PhantomParameterisation::G4PhantomParameteri     46 G4PhantomParameterisation::G4PhantomParameterisation()
                                                   >>  47   : fVoxelHalfX(0.), fVoxelHalfY(0.), fVoxelHalfZ(0.),
                                                   >>  48     fNoVoxelX(0), fNoVoxelY(0), fNoVoxelZ(0), fNoVoxelXY(0), fNoVoxel(0),
                                                   >>  49     fMaterialIndices(0), fContainerSolid(0),
                                                   >>  50     fContainerWallX(0.), fContainerWallY(0.), fContainerWallZ(0.),
                                                   >>  51     bSkipEqualMaterials(true)
 43 {                                                  52 {
 44   kCarTolerance = G4GeometryTolerance::GetInst     53   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
 45 }                                                  54 }
 46                                                    55 
 47                                                <<  56 
 48 //--------------------------------------------     57 //------------------------------------------------------------------
 49 G4PhantomParameterisation::~G4PhantomParameter <<  58 G4PhantomParameterisation::~G4PhantomParameterisation()
                                                   >>  59 {
                                                   >>  60 }
 50                                                    61 
 51                                                    62 
 52 //--------------------------------------------     63 //------------------------------------------------------------------
 53 void G4PhantomParameterisation::                   64 void G4PhantomParameterisation::
 54 BuildContainerSolid( G4VPhysicalVolume* pMothe <<  65 BuildContainerSolid( G4VPhysicalVolume *pMotherPhysical )
 55 {                                                  66 {
 56   fContainerSolid = pMotherPhysical->GetLogica     67   fContainerSolid = pMotherPhysical->GetLogicalVolume()->GetSolid();
 57   fContainerWallX = fNoVoxelsX * fVoxelHalfX;  <<  68   fContainerWallX = fNoVoxelX * fVoxelHalfX;
 58   fContainerWallY = fNoVoxelsY * fVoxelHalfY;  <<  69   fContainerWallY = fNoVoxelY * fVoxelHalfY;
 59   fContainerWallZ = fNoVoxelsZ * fVoxelHalfZ;  <<  70   fContainerWallZ = fNoVoxelZ * fVoxelHalfZ;
 60                                                    71 
 61   // CheckVoxelsFillContainer();                   72   // CheckVoxelsFillContainer();
 62 }                                                  73 }
 63                                                    74 
 64                                                << 
 65 //--------------------------------------------     75 //------------------------------------------------------------------
 66 void G4PhantomParameterisation::                   76 void G4PhantomParameterisation::
 67 BuildContainerSolid( G4VSolid* pMotherSolid )  <<  77 BuildContainerSolid( G4VSolid *pMotherSolid )
 68 {                                                  78 {
 69   fContainerSolid = pMotherSolid;                  79   fContainerSolid = pMotherSolid;
 70   fContainerWallX = fNoVoxelsX * fVoxelHalfX;  <<  80   fContainerWallX = fNoVoxelX * fVoxelHalfX;
 71   fContainerWallY = fNoVoxelsY * fVoxelHalfY;  <<  81   fContainerWallY = fNoVoxelY * fVoxelHalfY;
 72   fContainerWallZ = fNoVoxelsZ * fVoxelHalfZ;  <<  82   fContainerWallZ = fNoVoxelZ * fVoxelHalfZ;
 73                                                    83 
 74   // CheckVoxelsFillContainer();                   84   // CheckVoxelsFillContainer();
 75 }                                                  85 }
 76                                                    86 
 77                                                    87 
 78 //--------------------------------------------     88 //------------------------------------------------------------------
 79 void G4PhantomParameterisation::                   89 void G4PhantomParameterisation::
 80 ComputeTransformation(const G4int copyNo, G4VP <<  90 ComputeTransformation(const G4int copyNo, G4VPhysicalVolume *physVol ) const
 81 {                                                  91 {
 82   // Voxels cannot be rotated, return translat     92   // Voxels cannot be rotated, return translation
 83   //                                               93   //
 84   G4ThreeVector trans = GetTranslation( copyNo     94   G4ThreeVector trans = GetTranslation( copyNo );
 85                                                    95 
 86   physVol->SetTranslation( trans );                96   physVol->SetTranslation( trans );
 87 }                                                  97 }
 88                                                    98 
 89                                                    99 
 90 //--------------------------------------------    100 //------------------------------------------------------------------
 91 G4ThreeVector G4PhantomParameterisation::         101 G4ThreeVector G4PhantomParameterisation::
 92 GetTranslation(const G4int copyNo ) const         102 GetTranslation(const G4int copyNo ) const
 93 {                                                 103 {
 94   CheckCopyNo( copyNo );                          104   CheckCopyNo( copyNo );
 95                                                   105 
 96   std::size_t nx;                              << 106   size_t nx;
 97   std::size_t ny;                              << 107   size_t ny;
 98   std::size_t nz;                              << 108   size_t nz;
 99                                                   109 
100   ComputeVoxelIndices( copyNo, nx, ny, nz );      110   ComputeVoxelIndices( copyNo, nx, ny, nz );
101                                                   111 
102   G4ThreeVector trans( (2*nx+1)*fVoxelHalfX -     112   G4ThreeVector trans( (2*nx+1)*fVoxelHalfX - fContainerWallX,
103                        (2*ny+1)*fVoxelHalfY -     113                        (2*ny+1)*fVoxelHalfY - fContainerWallY,
104                        (2*nz+1)*fVoxelHalfZ -     114                        (2*nz+1)*fVoxelHalfZ - fContainerWallZ);
105   return trans;                                   115   return trans;
106 }                                                 116 }
107                                                   117 
108                                                   118 
109 //--------------------------------------------    119 //------------------------------------------------------------------
110 G4VSolid* G4PhantomParameterisation::             120 G4VSolid* G4PhantomParameterisation::
111 ComputeSolid(const G4int, G4VPhysicalVolume* p << 121 ComputeSolid(const G4int, G4VPhysicalVolume *pPhysicalVol) 
112 {                                                 122 {
113   return pPhysicalVol->GetLogicalVolume()->Get    123   return pPhysicalVol->GetLogicalVolume()->GetSolid();
114 }                                                 124 }
115                                                   125        
116                                                   126 
117 //--------------------------------------------    127 //------------------------------------------------------------------
118 G4Material* G4PhantomParameterisation::           128 G4Material* G4PhantomParameterisation::
119 ComputeMaterial(const G4int copyNo, G4VPhysica    129 ComputeMaterial(const G4int copyNo, G4VPhysicalVolume *, const G4VTouchable *) 
120 {                                                 130 { 
121   CheckCopyNo( copyNo );                          131   CheckCopyNo( copyNo );
122   std::size_t matIndex = GetMaterialIndex(copy << 132   size_t matIndex = GetMaterialIndex(copyNo);
123                                                   133 
124   return fMaterials[ matIndex ];                  134   return fMaterials[ matIndex ];
125 }                                                 135 }
126                                                   136 
127                                                   137 
128 //--------------------------------------------    138 //------------------------------------------------------------------
129 std::size_t G4PhantomParameterisation::        << 139 size_t G4PhantomParameterisation::
130 GetMaterialIndex( std::size_t copyNo ) const   << 140 GetMaterialIndex( size_t copyNo ) const
131 {                                                 141 {
132   CheckCopyNo( copyNo );                          142   CheckCopyNo( copyNo );
133                                                   143 
134   if( fMaterialIndices == nullptr ) { return 0 << 144   if( !fMaterialIndices ) { return 0; }
135   return *(fMaterialIndices+copyNo);              145   return *(fMaterialIndices+copyNo);
136 }                                                 146 }
137                                                   147 
138                                                   148 
139 //--------------------------------------------    149 //------------------------------------------------------------------
140 std::size_t G4PhantomParameterisation::        << 150 size_t G4PhantomParameterisation::
141 GetMaterialIndex( std::size_t nx, std::size_t  << 151 GetMaterialIndex( size_t nx, size_t ny, size_t nz ) const
142 {                                                 152 {
143   std::size_t copyNo = nx + fNoVoxelsX*ny + fN << 153   size_t copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
144   return GetMaterialIndex( copyNo );              154   return GetMaterialIndex( copyNo );
145 }                                                 155 }
146                                                   156 
147                                                   157 
148 //--------------------------------------------    158 //------------------------------------------------------------------
149 G4Material*                                    << 159 G4Material* G4PhantomParameterisation::GetMaterial( size_t nx, size_t ny, size_t nz) const
150 G4PhantomParameterisation::GetMaterial( std::s << 
151 {                                                 160 {
152   return fMaterials[GetMaterialIndex(nx,ny,nz)    161   return fMaterials[GetMaterialIndex(nx,ny,nz)];
153 }                                                 162 }
154                                                   163 
155                                                << 
156 //--------------------------------------------    164 //------------------------------------------------------------------
157 G4Material* G4PhantomParameterisation::GetMate << 165 G4Material* G4PhantomParameterisation::GetMaterial( size_t copyNo ) const
158 {                                                 166 {
159   return fMaterials[GetMaterialIndex(copyNo)];    167   return fMaterials[GetMaterialIndex(copyNo)];
160 }                                                 168 }
161                                                   169 
162                                                << 
163 //--------------------------------------------    170 //------------------------------------------------------------------
164 void G4PhantomParameterisation::                  171 void G4PhantomParameterisation::
165 ComputeVoxelIndices(const G4int copyNo, std::s << 172 ComputeVoxelIndices(const G4int copyNo, size_t& nx,
166                           std::size_t& ny, std << 173                             size_t& ny, size_t& nz ) const
167 {                                                 174 {
168   CheckCopyNo( copyNo );                          175   CheckCopyNo( copyNo );
169   nx = std::size_t(copyNo%fNoVoxelsX);         << 176   nx = size_t(copyNo%fNoVoxelX);
170   ny = std::size_t( (copyNo/fNoVoxelsX)%fNoVox << 177   ny = size_t( (copyNo/fNoVoxelX)%fNoVoxelY );
171   nz = std::size_t(copyNo/fNoVoxelsXY);        << 178   nz = size_t(copyNo/fNoVoxelXY);
172 }                                                 179 }
173                                                   180 
174                                                   181 
175 //--------------------------------------------    182 //------------------------------------------------------------------
176 void G4PhantomParameterisation::                  183 void G4PhantomParameterisation::
177 CheckVoxelsFillContainer( G4double contX, G4do    184 CheckVoxelsFillContainer( G4double contX, G4double contY, G4double contZ ) const
178 {                                                 185 {
179   G4double toleranceForWarning = 0.25*kCarTole    186   G4double toleranceForWarning = 0.25*kCarTolerance;
180                                                   187 
181   // Any bigger value than 0.25*kCarTolerance     188   // Any bigger value than 0.25*kCarTolerance will give a warning in
182   // G4NormalNavigation::ComputeStep(), becaus    189   // G4NormalNavigation::ComputeStep(), because the Inverse of a container
183   // translation that is Z+epsilon gives -Z+ep    190   // translation that is Z+epsilon gives -Z+epsilon (and the maximum tolerance
184   // in G4Box::Inside is 0.5*kCarTolerance        191   // in G4Box::Inside is 0.5*kCarTolerance
185   //                                              192   //
186   G4double toleranceForError = 1.*kCarToleranc    193   G4double toleranceForError = 1.*kCarTolerance;  
187                                                   194 
188   // Any bigger value than kCarTolerance will     195   // Any bigger value than kCarTolerance will give an error in GetReplicaNo()
189   //                                              196   //
190   if( std::fabs(contX-fNoVoxelsX*fVoxelHalfX)  << 197   if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForError
191    || std::fabs(contY-fNoVoxelsY*fVoxelHalfY)  << 198    || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForError  
192    || std::fabs(contZ-fNoVoxelsZ*fVoxelHalfZ)  << 199    || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForError )
193   {                                               200   {
194     std::ostringstream message;                   201     std::ostringstream message;
195     message << "Voxels do not fully fill the c    202     message << "Voxels do not fully fill the container: "
196             << fContainerSolid->GetName() << G    203             << fContainerSolid->GetName() << G4endl
197             << "        DiffX= " << contX-fNoV << 204             << "        DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
198             << "        DiffY= " << contY-fNoV << 205             << "        DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
199             << "        DiffZ= " << contZ-fNoV << 206             << "        DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
200             << "        Maximum difference is:    207             << "        Maximum difference is: " << toleranceForError;
201     G4Exception("G4PhantomParameterisation::Ch    208     G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
202                 "GeomNav0002", FatalException,    209                 "GeomNav0002", FatalException, message);
203                                                   210 
204   }                                               211   }
205   else if( std::fabs(contX-fNoVoxelsX*fVoxelHa << 212   else if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForWarning
206         || std::fabs(contY-fNoVoxelsY*fVoxelHa << 213         || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForWarning  
207         || std::fabs(contZ-fNoVoxelsZ*fVoxelHa << 214         || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForWarning )
208   {                                               215   {
209     std::ostringstream message;                   216     std::ostringstream message;
210     message << "Voxels do not fully fill the c    217     message << "Voxels do not fully fill the container: "
211             << fContainerSolid->GetName() << G    218             << fContainerSolid->GetName() << G4endl
212             << "          DiffX= " << contX-fN << 219             << "          DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
213             << "          DiffY= " << contY-fN << 220             << "          DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
214             << "          DiffZ= " << contZ-fN << 221             << "          DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
215             << "          Maximum difference i    222             << "          Maximum difference is: " << toleranceForWarning;
216     G4Exception("G4PhantomParameterisation::Ch    223     G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
217                 "GeomNav1002", JustWarning, me    224                 "GeomNav1002", JustWarning, message);
218   }                                               225   }
219 }                                                 226 }
220                                                   227   
221                                                   228  
222 //--------------------------------------------    229 //------------------------------------------------------------------
223 G4int G4PhantomParameterisation::                 230 G4int G4PhantomParameterisation::
224 GetReplicaNo( const G4ThreeVector& localPoint,    231 GetReplicaNo( const G4ThreeVector& localPoint, const G4ThreeVector& localDir )
225 {                                                 232 {
226                                                   233 
227   // Check first that point is really inside v    234   // Check first that point is really inside voxels
228   //                                              235   //
229   if( fContainerSolid->Inside( localPoint ) ==    236   if( fContainerSolid->Inside( localPoint ) == kOutside )
230   {                                               237   {
231     if( std::fabs(localPoint.x()) - fContainer << 238     std::ostringstream message;
232   && std::fabs(localPoint.y()) - fContainerWal << 239     message << "Point outside voxels!" << G4endl
233   && std::fabs(localPoint.z()) - fContainerWal << 240             << "        localPoint - " << localPoint
234     {                                          << 241             << " - is outside container solid: "
235       std::ostringstream message;              << 242             << fContainerSolid->GetName() << G4endl;
236       message << "Point outside voxels!" << G4 << 243     G4Exception("G4PhantomParameterisation::GetReplicaNo()", "GeomNav0003",
237         << "        localPoint - " << localPoi << 244                 FatalErrorInArgument, message);
238         << " - is outside container solid: "   << 
239         << fContainerSolid->GetName() << G4end << 
240         << "DIFFERENCE WITH PHANTOM WALLS X: " << 
241         << std::fabs(localPoint.x()) - fContai << 
242         << " Y: " << std::fabs(localPoint.y()) << 
243         << " Z: " << std::fabs(localPoint.z()) << 
244       G4Exception("G4PhantomParameterisation:: << 
245       FatalErrorInArgument, message);          << 
246     }                                          << 
247   }                                               245   }
248                                                   246   
249   // Check the voxel numbers corresponding to     247   // Check the voxel numbers corresponding to localPoint
250   // When a particle is on a surface, it may b    248   // When a particle is on a surface, it may be between -kCarTolerance and
251   // +kCartolerance. By a simple distance as:     249   // +kCartolerance. By a simple distance as:
252   //   G4int nx = G4int( (localPoint.x()+)/fVo    250   //   G4int nx = G4int( (localPoint.x()+)/fVoxelHalfX/2.);
253   // those between -kCartolerance and 0 will b    251   // those between -kCartolerance and 0 will be placed on voxel N-1 and those
254   // between 0 and kCarTolerance on voxel N.      252   // between 0 and kCarTolerance on voxel N.
255   // To avoid precision problems place the tra    253   // To avoid precision problems place the tracks that are on the surface on
256   // voxel N-1 if they have negative direction    254   // voxel N-1 if they have negative direction and on voxel N if they have
257   // positive direction.                          255   // positive direction.
258   // Add +kCarTolerance so that they are first    256   // Add +kCarTolerance so that they are first placed on voxel N, and then
259   // if the direction is negative substract 1     257   // if the direction is negative substract 1
260                                                   258 
261   G4double fx = (localPoint.x()+fContainerWall    259   G4double fx = (localPoint.x()+fContainerWallX+kCarTolerance)/(fVoxelHalfX*2.);
262   auto  nx = G4int(fx);                        << 260   G4int nx = G4int(fx);
263                                                   261 
264   G4double fy = (localPoint.y()+fContainerWall << 262   G4double fy = (localPoint.y()+fContainerWallY+kCarTolerance)/(fVoxelHalfY*2.);
265   auto  ny = G4int(fy);                        << 263   G4int ny = G4int(fy);
266                                                   264 
267   G4double fz = (localPoint.z()+fContainerWall    265   G4double fz = (localPoint.z()+fContainerWallZ+kCarTolerance)/(fVoxelHalfZ*2.);
268   auto  nz = G4int(fz);                        << 266   G4int nz = G4int(fz);
269                                                   267 
270   // If it is on the surface side, check the d    268   // If it is on the surface side, check the direction: if direction is
271   // negative place it in the previous voxel ( << 269   // negative place it on the previous voxel (if direction is positive it is
272   // already in the next voxel).               << 270   // already in the next voxel...). 
273   // Correct also cases where n = -1 or n = fN << 271   // Correct also cases where n = -1 or n = fNoVoxel. It is always traced to be
274   // due to multiple scattering: track is ente    272   // due to multiple scattering: track is entering a voxel but multiple
275   // scattering changes the angle towards outs    273   // scattering changes the angle towards outside
276   //                                              274   //
277   if( fx - nx < kCarTolerance*fVoxelHalfX )    << 275   if( fx - nx < kCarTolerance/fVoxelHalfX )
278   {                                               276   {
279     if( localDir.x() < 0 )                        277     if( localDir.x() < 0 )
280     {                                             278     {
281       if( nx != 0 )                               279       if( nx != 0 )
282       {                                           280       {
283         nx -= 1;                                  281         nx -= 1;
284       }                                           282       }
285     }                                             283     }
286     else                                          284     else
287     {                                             285     {
288       if( nx == G4int(fNoVoxelsX) )            << 286       if( nx == G4int(fNoVoxelX) )  
289       {                                           287       {
290         nx -= 1;                                  288         nx -= 1;       
291       }                                           289       }
292     }                                             290     }
293   }                                               291   }
294   if( fy - ny < kCarTolerance*fVoxelHalfY )    << 292   if( fy - ny < kCarTolerance/fVoxelHalfY )
295   {                                               293   {
296     if( localDir.y() < 0 )                        294     if( localDir.y() < 0 )
297     {                                             295     {
298       if( ny != 0 )                               296       if( ny != 0 )
299       {                                           297       {
300         ny -= 1;                                  298         ny -= 1;
301       }                                           299       }
302     }                                             300     }
303     else                                          301     else
304     {                                             302     {
305       if( ny == G4int(fNoVoxelsY) )            << 303       if( ny == G4int(fNoVoxelY) )  
306       {                                           304       {
307         ny -= 1;                                  305         ny -= 1;       
308       }                                           306       }
309     }                                             307     }
310   }                                               308   }
311   if( fz - nz < kCarTolerance*fVoxelHalfZ )    << 309   if( fz - nz < kCarTolerance/fVoxelHalfZ )
312   {                                               310   {
313     if( localDir.z() < 0 )                        311     if( localDir.z() < 0 )
314     {                                             312     {
315       if( nz != 0 )                               313       if( nz != 0 )
316       {                                           314       {
317         nz -= 1;                                  315         nz -= 1;
318       }                                           316       }
319     }                                             317     }
320     else                                          318     else
321     {                                             319     {
322       if( nz == G4int(fNoVoxelsZ) )            << 320       if( nz == G4int(fNoVoxelZ) )  
323       {                                           321       {
324         nz -= 1;                                  322         nz -= 1;       
325       }                                           323       }
326     }                                             324     }
327   }                                               325   }
328                                                   326   
329   auto  copyNo = G4int(nx + fNoVoxelsX*ny + fN << 327   G4int copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
330                                                   328 
331   // Check if there are still errors              329   // Check if there are still errors 
332   //                                              330   //
333   G4bool isOK = true;                             331   G4bool isOK = true;
334   if( nx < 0 )                                    332   if( nx < 0 )
335   {                                               333   {
336     nx = 0;                                       334     nx = 0;
337     isOK = false;                                 335     isOK = false;
338   }                                               336   }
339   else if( nx >= G4int(fNoVoxelsX) )           << 337   else if( nx >= G4int(fNoVoxelX) )
340   {                                               338   {
341     nx = G4int(fNoVoxelsX)-1;                  << 339     nx = fNoVoxelX-1;
342     isOK = false;                                 340     isOK = false;
343   }                                               341   }
344   if( ny < 0 )                                    342   if( ny < 0 )
345   {                                               343   {
346     ny = 0;                                       344     ny = 0;
347     isOK = false;                                 345     isOK = false;
348   }                                               346   }
349   else if( ny >= G4int(fNoVoxelsY) )           << 347   else if( ny >= G4int(fNoVoxelY) )
350   {                                               348   {
351     ny = G4int(fNoVoxelsY)-1;                  << 349     ny = fNoVoxelY-1;
352     isOK = false;                                 350     isOK = false;
353   }                                               351   }
354   if( nz < 0 )                                    352   if( nz < 0 )
355   {                                               353   {
356     nz = 0;                                       354     nz = 0;
357     isOK = false;                                 355     isOK = false;
358   }                                               356   }
359   else if( nz >= G4int(fNoVoxelsZ) )           << 357   else if( nz >= G4int(fNoVoxelZ) )
360   {                                               358   {
361     nz = G4int(fNoVoxelsZ)-1;                  << 359     nz = fNoVoxelZ-1;
362     isOK = false;                                 360     isOK = false;
363   }                                               361   }
364   if( !isOK )                                     362   if( !isOK )
365   {                                               363   {
366     if( std::fabs(localPoint.x()-fContainerWal << 364     std::ostringstream message;
367   std::fabs(localPoint.y()-fContainerWallY) >  << 365     message << "Corrected the copy number! It was negative or too big" << G4endl
368   std::fabs(localPoint.z()-fContainerWallZ) >  << 366             << "          LocalPoint: " << localPoint << G4endl
369       std::ostringstream message;              << 367             << "          LocalDir: " << localDir << G4endl
370       message << "Corrected the copy number! I << 368             << "          Voxel container size: " << fContainerWallX
371         << "          LocalPoint: " << localPo << 369             << " " << fContainerWallY << " " << fContainerWallZ << G4endl
372         << "          LocalDir: " << localDir  << 370             << "          LocalPoint - wall: "
373         << "          Voxel container size: "  << 371             << localPoint.x()-fContainerWallX << " "
374         << " " << fContainerWallY << " " << fC << 372             << localPoint.y()-fContainerWallY << " "
375         << "          LocalPoint - wall: "     << 373             << localPoint.z()-fContainerWallZ;
376         << localPoint.x()-fContainerWallX << " << 374     G4Exception("G4PhantomParameterisation::GetReplicaNo()",
377         << localPoint.y()-fContainerWallY << " << 375                 "GeomNav1002", JustWarning, message);
378         << localPoint.z()-fContainerWallZ;     << 376     copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
379       G4Exception("G4PhantomParameterisation:: << 
380       "GeomNav1002", JustWarning, message);    << 
381     }                                          << 
382                                                << 
383     copyNo = G4int(nx + fNoVoxelsX*ny + fNoVox << 
384   }                                               377   }
385                                                   378 
                                                   >> 379   // CheckCopyNo( copyNo ); // not needed, just for debugging code
                                                   >> 380 
386   return copyNo;                                  381   return copyNo;
387 }                                                 382 }
388                                                   383 
389                                                   384 
390 //--------------------------------------------    385 //------------------------------------------------------------------
391 void G4PhantomParameterisation::CheckCopyNo( c << 386 void G4PhantomParameterisation::CheckCopyNo( const G4int copyNo ) const
392 {                                                 387 { 
393   if( copyNo < 0 || copyNo >= G4int(fNoVoxels) << 388   if( copyNo < 0 || copyNo >= G4int(fNoVoxel) )
394   {                                               389   {
395     std::ostringstream message;                   390     std::ostringstream message;
396     message << "Copy number is negative or too    391     message << "Copy number is negative or too big!" << G4endl
397             << "        Copy number: " << copy    392             << "        Copy number: " << copyNo << G4endl
398             << "        Total number of voxels << 393             << "        Total number of voxels: " << fNoVoxel;
399     G4Exception("G4PhantomParameterisation::Ch    394     G4Exception("G4PhantomParameterisation::CheckCopyNo()",
400                 "GeomNav0002", FatalErrorInArg    395                 "GeomNav0002", FatalErrorInArgument, message);
401   }                                               396   }
402 }                                                 397 }
403                                                   398