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 9.1.p2)


  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,v 1.4 2008/01/22 15:02:36 gcosmo Exp $
                                                   >>  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()
 43 {                                                  47 {
                                                   >>  48   // Initialise data
                                                   >>  49   //
                                                   >>  50   fMaterialIndices = 0;
                                                   >>  51   fContainerWallX = 0.;
                                                   >>  52   fContainerWallY = 0.;
                                                   >>  53   fContainerWallZ = 0.;
                                                   >>  54 
 44   kCarTolerance = G4GeometryTolerance::GetInst     55   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
 45 }                                              << 
 46                                                    56 
 47                                                <<  57   bSkipEqualMaterials = 1;
 48 //-------------------------------------------- <<  58 }
 49 G4PhantomParameterisation::~G4PhantomParameter << 
 50                                                    59 
 51                                                    60 
 52 //--------------------------------------------     61 //------------------------------------------------------------------
 53 void G4PhantomParameterisation::               <<  62 G4PhantomParameterisation::~G4PhantomParameterisation()
 54 BuildContainerSolid( G4VPhysicalVolume* pMothe << 
 55 {                                                  63 {
 56   fContainerSolid = pMotherPhysical->GetLogica << 
 57   fContainerWallX = fNoVoxelsX * fVoxelHalfX;  << 
 58   fContainerWallY = fNoVoxelsY * fVoxelHalfY;  << 
 59   fContainerWallZ = fNoVoxelsZ * fVoxelHalfZ;  << 
 60                                                << 
 61   // CheckVoxelsFillContainer();               << 
 62 }                                                  64 }
 63                                                    65 
 64                                                    66 
 65 //--------------------------------------------     67 //------------------------------------------------------------------
 66 void G4PhantomParameterisation::                   68 void G4PhantomParameterisation::
 67 BuildContainerSolid( G4VSolid* pMotherSolid )  <<  69 BuildContainerSolid( G4VPhysicalVolume *pMotherPhysical )
 68 {                                                  70 {
 69   fContainerSolid = pMotherSolid;              <<  71   fContainerSolid = pMotherPhysical->GetLogicalVolume()->GetSolid();
 70   fContainerWallX = fNoVoxelsX * fVoxelHalfX;  <<  72   fContainerWallX = fNoVoxelX * fVoxelHalfX;
 71   fContainerWallY = fNoVoxelsY * fVoxelHalfY;  <<  73   fContainerWallY = fNoVoxelY * fVoxelHalfY;
 72   fContainerWallZ = fNoVoxelsZ * fVoxelHalfZ;  <<  74   fContainerWallZ = fNoVoxelZ * fVoxelHalfZ;
 73                                                    75 
 74   // CheckVoxelsFillContainer();                   76   // CheckVoxelsFillContainer();
 75 }                                                  77 }
 76                                                    78 
 77                                                    79 
 78 //--------------------------------------------     80 //------------------------------------------------------------------
 79 void G4PhantomParameterisation::                   81 void G4PhantomParameterisation::
 80 ComputeTransformation(const G4int copyNo, G4VP <<  82 ComputeTransformation(const G4int copyNo, G4VPhysicalVolume *physVol ) const
 81 {                                                  83 {
 82   // Voxels cannot be rotated, return translat     84   // Voxels cannot be rotated, return translation
 83   //                                               85   //
 84   G4ThreeVector trans = GetTranslation( copyNo     86   G4ThreeVector trans = GetTranslation( copyNo );
 85                                                    87 
 86   physVol->SetTranslation( trans );                88   physVol->SetTranslation( trans );
 87 }                                                  89 }
 88                                                    90 
 89                                                    91 
 90 //--------------------------------------------     92 //------------------------------------------------------------------
 91 G4ThreeVector G4PhantomParameterisation::          93 G4ThreeVector G4PhantomParameterisation::
 92 GetTranslation(const G4int copyNo ) const          94 GetTranslation(const G4int copyNo ) const
 93 {                                                  95 {
 94   CheckCopyNo( copyNo );                           96   CheckCopyNo( copyNo );
 95                                                    97 
 96   std::size_t nx;                              <<  98   size_t nx;
 97   std::size_t ny;                              <<  99   size_t ny;
 98   std::size_t nz;                              << 100   size_t nz;
 99                                                   101 
100   ComputeVoxelIndices( copyNo, nx, ny, nz );      102   ComputeVoxelIndices( copyNo, nx, ny, nz );
101                                                   103 
102   G4ThreeVector trans( (2*nx+1)*fVoxelHalfX -     104   G4ThreeVector trans( (2*nx+1)*fVoxelHalfX - fContainerWallX,
103                        (2*ny+1)*fVoxelHalfY -     105                        (2*ny+1)*fVoxelHalfY - fContainerWallY,
104                        (2*nz+1)*fVoxelHalfZ -     106                        (2*nz+1)*fVoxelHalfZ - fContainerWallZ);
105   return trans;                                   107   return trans;
106 }                                                 108 }
107                                                   109 
108                                                   110 
109 //--------------------------------------------    111 //------------------------------------------------------------------
110 G4VSolid* G4PhantomParameterisation::             112 G4VSolid* G4PhantomParameterisation::
111 ComputeSolid(const G4int, G4VPhysicalVolume* p << 113 ComputeSolid(const G4int, G4VPhysicalVolume *pPhysicalVol) 
112 {                                                 114 {
113   return pPhysicalVol->GetLogicalVolume()->Get    115   return pPhysicalVol->GetLogicalVolume()->GetSolid();
114 }                                                 116 }
115                                                   117        
116                                                   118 
117 //--------------------------------------------    119 //------------------------------------------------------------------
118 G4Material* G4PhantomParameterisation::           120 G4Material* G4PhantomParameterisation::
119 ComputeMaterial(const G4int copyNo, G4VPhysica    121 ComputeMaterial(const G4int copyNo, G4VPhysicalVolume *, const G4VTouchable *) 
120 {                                                 122 { 
121   CheckCopyNo( copyNo );                          123   CheckCopyNo( copyNo );
122   std::size_t matIndex = GetMaterialIndex(copy << 124   size_t matIndex = GetMaterialIndex(copyNo);
123                                                   125 
124   return fMaterials[ matIndex ];                  126   return fMaterials[ matIndex ];
125 }                                                 127 }
126                                                   128 
127                                                   129 
128 //--------------------------------------------    130 //------------------------------------------------------------------
129 std::size_t G4PhantomParameterisation::        << 131 size_t G4PhantomParameterisation::
130 GetMaterialIndex( std::size_t copyNo ) const   << 132 GetMaterialIndex( size_t copyNo ) const
131 {                                                 133 {
132   CheckCopyNo( copyNo );                          134   CheckCopyNo( copyNo );
133                                                   135 
134   if( fMaterialIndices == nullptr ) { return 0 << 
135   return *(fMaterialIndices+copyNo);              136   return *(fMaterialIndices+copyNo);
136 }                                                 137 }
137                                                   138 
138                                                   139 
139 //--------------------------------------------    140 //------------------------------------------------------------------
140 std::size_t G4PhantomParameterisation::        << 141 size_t G4PhantomParameterisation::
141 GetMaterialIndex( std::size_t nx, std::size_t  << 142 GetMaterialIndex( size_t nx, size_t ny, size_t nz ) const
142 {                                                 143 {
143   std::size_t copyNo = nx + fNoVoxelsX*ny + fN << 144   size_t copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
144   return GetMaterialIndex( copyNo );              145   return GetMaterialIndex( copyNo );
145 }                                                 146 }
146                                                   147 
147                                                   148 
148 //--------------------------------------------    149 //------------------------------------------------------------------
149 G4Material*                                    << 150 G4Material* G4PhantomParameterisation::GetMaterial( size_t nx, size_t ny, size_t nz) const
150 G4PhantomParameterisation::GetMaterial( std::s << 
151 {                                                 151 {
152   return fMaterials[GetMaterialIndex(nx,ny,nz)    152   return fMaterials[GetMaterialIndex(nx,ny,nz)];
153 }                                                 153 }
154                                                   154 
155                                                << 
156 //--------------------------------------------    155 //------------------------------------------------------------------
157 G4Material* G4PhantomParameterisation::GetMate << 156 G4Material* G4PhantomParameterisation::GetMaterial( size_t copyNo ) const
158 {                                                 157 {
159   return fMaterials[GetMaterialIndex(copyNo)];    158   return fMaterials[GetMaterialIndex(copyNo)];
160 }                                                 159 }
161                                                   160 
162                                                << 
163 //--------------------------------------------    161 //------------------------------------------------------------------
164 void G4PhantomParameterisation::                  162 void G4PhantomParameterisation::
165 ComputeVoxelIndices(const G4int copyNo, std::s << 163 ComputeVoxelIndices(const G4int copyNo, size_t& nx,
166                           std::size_t& ny, std << 164                             size_t& ny, size_t& nz ) const
167 {                                                 165 {
168   CheckCopyNo( copyNo );                          166   CheckCopyNo( copyNo );
169   nx = std::size_t(copyNo%fNoVoxelsX);         << 167   nx = size_t(copyNo%fNoVoxelX);
170   ny = std::size_t( (copyNo/fNoVoxelsX)%fNoVox << 168   ny = size_t( (copyNo/fNoVoxelX)%fNoVoxelY );
171   nz = std::size_t(copyNo/fNoVoxelsXY);        << 169   nz = size_t(copyNo/fNoVoxelXY);
172 }                                                 170 }
173                                                   171 
174                                                   172 
175 //--------------------------------------------    173 //------------------------------------------------------------------
176 void G4PhantomParameterisation::                  174 void G4PhantomParameterisation::
177 CheckVoxelsFillContainer( G4double contX, G4do    175 CheckVoxelsFillContainer( G4double contX, G4double contY, G4double contZ ) const
178 {                                                 176 {
179   G4double toleranceForWarning = 0.25*kCarTole    177   G4double toleranceForWarning = 0.25*kCarTolerance;
180                                                   178 
181   // Any bigger value than 0.25*kCarTolerance     179   // Any bigger value than 0.25*kCarTolerance will give a warning in
182   // G4NormalNavigation::ComputeStep(), becaus    180   // G4NormalNavigation::ComputeStep(), because the Inverse of a container
183   // translation that is Z+epsilon gives -Z+ep    181   // translation that is Z+epsilon gives -Z+epsilon (and the maximum tolerance
184   // in G4Box::Inside is 0.5*kCarTolerance        182   // in G4Box::Inside is 0.5*kCarTolerance
185   //                                              183   //
186   G4double toleranceForError = 1.*kCarToleranc    184   G4double toleranceForError = 1.*kCarTolerance;  
187                                                   185 
188   // Any bigger value than kCarTolerance will     186   // Any bigger value than kCarTolerance will give an error in GetReplicaNo()
189   //                                              187   //
190   if( std::fabs(contX-fNoVoxelsX*fVoxelHalfX)  << 188   if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForError
191    || std::fabs(contY-fNoVoxelsY*fVoxelHalfY)  << 189    || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForError  
192    || std::fabs(contZ-fNoVoxelsZ*fVoxelHalfZ)  << 190    || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForError )
193   {                                            << 191   {
194     std::ostringstream message;                << 192     G4cerr << "ERROR - G4PhantomParameterisation::CheckVoxelsFillContainer()"
195     message << "Voxels do not fully fill the c << 193            << G4endl
196             << fContainerSolid->GetName() << G << 194            << "        Voxels do not fully fill the container: "
197             << "        DiffX= " << contX-fNoV << 195            << fContainerSolid->GetName() << G4endl
198             << "        DiffY= " << contY-fNoV << 196            << "        DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
199             << "        DiffZ= " << contZ-fNoV << 197            << "        DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
200             << "        Maximum difference is: << 198            << "        DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
                                                   >> 199            << "        Maximum difference is: " << toleranceForError << G4endl;
201     G4Exception("G4PhantomParameterisation::Ch    200     G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
202                 "GeomNav0002", FatalException, << 201                 "InvalidSetup", FatalException,
                                                   >> 202                 "Voxels do not fully fill the container!");
203                                                   203 
204   }                                               204   }
205   else if( std::fabs(contX-fNoVoxelsX*fVoxelHa << 205   else if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForWarning
206         || std::fabs(contY-fNoVoxelsY*fVoxelHa << 206         || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForWarning  
207         || std::fabs(contZ-fNoVoxelsZ*fVoxelHa << 207         || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForWarning )
208   {                                            << 208   {
209     std::ostringstream message;                << 209     G4cerr << "WARNING - G4PhantomParameterisation::CheckVoxelsFillContainer()"
210     message << "Voxels do not fully fill the c << 210            << G4endl
211             << fContainerSolid->GetName() << G << 211            << "          Voxels do not fully fill the container: "
212             << "          DiffX= " << contX-fN << 212            << fContainerSolid->GetName() << G4endl
213             << "          DiffY= " << contY-fN << 213            << "          DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
214             << "          DiffZ= " << contZ-fN << 214            << "          DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
215             << "          Maximum difference i << 215            << "          DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
                                                   >> 216            << "          Maximum difference is: " << toleranceForWarning
                                                   >> 217            << G4endl;
216     G4Exception("G4PhantomParameterisation::Ch    218     G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
217                 "GeomNav1002", JustWarning, me << 219                 "InvalidSetup", JustWarning,
                                                   >> 220                 "Voxels do not fully fill the container!");
218   }                                               221   }
219 }                                                 222 }
220                                                   223   
221                                                   224  
222 //--------------------------------------------    225 //------------------------------------------------------------------
223 G4int G4PhantomParameterisation::                 226 G4int G4PhantomParameterisation::
224 GetReplicaNo( const G4ThreeVector& localPoint,    227 GetReplicaNo( const G4ThreeVector& localPoint, const G4ThreeVector& localDir )
225 {                                                 228 {
226                                                   229 
227   // Check first that point is really inside v    230   // Check first that point is really inside voxels
228   //                                              231   //
229   if( fContainerSolid->Inside( localPoint ) ==    232   if( fContainerSolid->Inside( localPoint ) == kOutside )
230   {                                               233   {
231     if( std::fabs(localPoint.x()) - fContainer << 234     G4cerr << "ERROR - G4PhantomParameterisation::GetReplicaNo()" << G4endl
232   && std::fabs(localPoint.y()) - fContainerWal << 235            << "        localPoint - " << localPoint
233   && std::fabs(localPoint.z()) - fContainerWal << 236            << " - is outside container solid: "
234     {                                          << 237            << fContainerSolid->GetName() << G4endl;
235       std::ostringstream message;              << 238     G4Exception("G4PhantomParameterisation::GetReplicaNo()", "InvalidSetup",
236       message << "Point outside voxels!" << G4 << 239                 FatalErrorInArgument, "Point outside voxels!");
237         << "        localPoint - " << localPoi << 
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   }                                               240   }
248                                                   241   
249   // Check the voxel numbers corresponding to     242   // Check the voxel numbers corresponding to localPoint
250   // When a particle is on a surface, it may b    243   // When a particle is on a surface, it may be between -kCarTolerance and
251   // +kCartolerance. By a simple distance as:     244   // +kCartolerance. By a simple distance as:
252   //   G4int nx = G4int( (localPoint.x()+)/fVo    245   //   G4int nx = G4int( (localPoint.x()+)/fVoxelHalfX/2.);
253   // those between -kCartolerance and 0 will b    246   // those between -kCartolerance and 0 will be placed on voxel N-1 and those
254   // between 0 and kCarTolerance on voxel N.      247   // between 0 and kCarTolerance on voxel N.
255   // To avoid precision problems place the tra    248   // To avoid precision problems place the tracks that are on the surface on
256   // voxel N-1 if they have negative direction    249   // voxel N-1 if they have negative direction and on voxel N if they have
257   // positive direction.                          250   // positive direction.
258   // Add +kCarTolerance so that they are first    251   // Add +kCarTolerance so that they are first placed on voxel N, and then
259   // if the direction is negative substract 1     252   // if the direction is negative substract 1
260                                                   253 
261   G4double fx = (localPoint.x()+fContainerWall    254   G4double fx = (localPoint.x()+fContainerWallX+kCarTolerance)/(fVoxelHalfX*2.);
262   auto  nx = G4int(fx);                        << 255   G4int nx = G4int(fx);
263                                                   256 
264   G4double fy = (localPoint.y()+fContainerWall << 257   G4double fy = (localPoint.y()+fContainerWallY+kCarTolerance)/(fVoxelHalfY*2.);
265   auto  ny = G4int(fy);                        << 258   G4int ny = G4int(fy);
266                                                   259 
267   G4double fz = (localPoint.z()+fContainerWall    260   G4double fz = (localPoint.z()+fContainerWallZ+kCarTolerance)/(fVoxelHalfZ*2.);
268   auto  nz = G4int(fz);                        << 261   G4int nz = G4int(fz);
269                                                   262 
270   // If it is on the surface side, check the d    263   // If it is on the surface side, check the direction: if direction is
271   // negative place it in the previous voxel ( << 264   // negative place it on the previous voxel (if direction is positive it is
272   // already in the next voxel).               << 265   // already in the next voxel...). 
273   // Correct also cases where n = -1 or n = fN << 266   // Correct also cases where n = -1 or n = fNoVoxel. It is always traced to be
274   // due to multiple scattering: track is ente    267   // due to multiple scattering: track is entering a voxel but multiple
275   // scattering changes the angle towards outs    268   // scattering changes the angle towards outside
276   //                                              269   //
277   if( fx - nx < kCarTolerance*fVoxelHalfX )    << 270   if( fx - nx < kCarTolerance/fVoxelHalfX )
278   {                                               271   {
279     if( localDir.x() < 0 )                        272     if( localDir.x() < 0 )
280     {                                             273     {
281       if( nx != 0 )                               274       if( nx != 0 )
282       {                                           275       {
283         nx -= 1;                                  276         nx -= 1;
284       }                                           277       }
285     }                                             278     }
286     else                                          279     else
287     {                                             280     {
288       if( nx == G4int(fNoVoxelsX) )            << 281       if( nx == G4int(fNoVoxelX) )  
289       {                                           282       {
290         nx -= 1;                                  283         nx -= 1;       
291       }                                           284       }
292     }                                             285     }
293   }                                               286   }
294   if( fy - ny < kCarTolerance*fVoxelHalfY )    << 287   if( fy - ny < kCarTolerance/fVoxelHalfY )
295   {                                               288   {
296     if( localDir.y() < 0 )                        289     if( localDir.y() < 0 )
297     {                                             290     {
298       if( ny != 0 )                               291       if( ny != 0 )
299       {                                           292       {
300         ny -= 1;                                  293         ny -= 1;
301       }                                           294       }
302     }                                             295     }
303     else                                          296     else
304     {                                             297     {
305       if( ny == G4int(fNoVoxelsY) )            << 298       if( ny == G4int(fNoVoxelY) )  
306       {                                           299       {
307         ny -= 1;                                  300         ny -= 1;       
308       }                                           301       }
309     }                                             302     }
310   }                                               303   }
311   if( fz - nz < kCarTolerance*fVoxelHalfZ )    << 304   if( fz - nz < kCarTolerance/fVoxelHalfZ )
312   {                                               305   {
313     if( localDir.z() < 0 )                        306     if( localDir.z() < 0 )
314     {                                             307     {
315       if( nz != 0 )                               308       if( nz != 0 )
316       {                                           309       {
317         nz -= 1;                                  310         nz -= 1;
318       }                                           311       }
319     }                                             312     }
320     else                                          313     else
321     {                                             314     {
322       if( nz == G4int(fNoVoxelsZ) )            << 315       if( nz == G4int(fNoVoxelZ) )  
323       {                                           316       {
324         nz -= 1;                                  317         nz -= 1;       
325       }                                           318       }
326     }                                             319     }
327   }                                               320   }
328                                                   321   
329   auto  copyNo = G4int(nx + fNoVoxelsX*ny + fN << 322   G4int copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
330                                                   323 
331   // Check if there are still errors              324   // Check if there are still errors 
332   //                                              325   //
333   G4bool isOK = true;                             326   G4bool isOK = true;
334   if( nx < 0 )                                    327   if( nx < 0 )
335   {                                               328   {
336     nx = 0;                                       329     nx = 0;
337     isOK = false;                                 330     isOK = false;
338   }                                               331   }
339   else if( nx >= G4int(fNoVoxelsX) )           << 332   else if( nx >= G4int(fNoVoxelX) )
340   {                                               333   {
341     nx = G4int(fNoVoxelsX)-1;                  << 334     nx = fNoVoxelX-1;
342     isOK = false;                                 335     isOK = false;
343   }                                               336   }
344   if( ny < 0 )                                    337   if( ny < 0 )
345   {                                               338   {
346     ny = 0;                                       339     ny = 0;
347     isOK = false;                                 340     isOK = false;
348   }                                               341   }
349   else if( ny >= G4int(fNoVoxelsY) )           << 342   else if( ny >= G4int(fNoVoxelY) )
350   {                                               343   {
351     ny = G4int(fNoVoxelsY)-1;                  << 344     ny = fNoVoxelY-1;
352     isOK = false;                                 345     isOK = false;
353   }                                               346   }
354   if( nz < 0 )                                    347   if( nz < 0 )
355   {                                               348   {
356     nz = 0;                                       349     nz = 0;
357     isOK = false;                                 350     isOK = false;
358   }                                               351   }
359   else if( nz >= G4int(fNoVoxelsZ) )           << 352   else if( nz >= G4int(fNoVoxelZ) )
360   {                                               353   {
361     nz = G4int(fNoVoxelsZ)-1;                  << 354     nz = fNoVoxelZ-1;
362     isOK = false;                                 355     isOK = false;
363   }                                               356   }
364   if( !isOK )                                     357   if( !isOK )
365   {                                               358   {
366     if( std::fabs(localPoint.x()-fContainerWal << 359     G4cerr << "WARNING - G4PhantomParameterisation::GetReplicaNo()" << G4endl
367   std::fabs(localPoint.y()-fContainerWallY) >  << 360            << "          LocalPoint: " << localPoint << G4endl
368   std::fabs(localPoint.z()-fContainerWallZ) >  << 361            << "          LocalDir: " << localDir << G4endl
369       std::ostringstream message;              << 362            << "          Voxel container size: " << fContainerWallX
370       message << "Corrected the copy number! I << 363            << " " << fContainerWallY << " " << fContainerWallZ << G4endl
371         << "          LocalPoint: " << localPo << 364            << "          LocalPoint - wall: "
372         << "          LocalDir: " << localDir  << 365            << localPoint.x()-fContainerWallX << " "
373         << "          Voxel container size: "  << 366            << localPoint.y()-fContainerWallY << " "
374         << " " << fContainerWallY << " " << fC << 367            << localPoint.z()-fContainerWallZ << G4endl;
375         << "          LocalPoint - wall: "     << 368     G4Exception("G4PhantomParameterisation::GetReplicaNo()",
376         << localPoint.x()-fContainerWallX << " << 369                 "Wrong-copy-number", JustWarning,
377         << localPoint.y()-fContainerWallY << " << 370                 "Corrected the copy number! It was negative or too big");
378         << localPoint.z()-fContainerWallZ;     << 371     copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
379       G4Exception("G4PhantomParameterisation:: << 
380       "GeomNav1002", JustWarning, message);    << 
381     }                                          << 
382                                                << 
383     copyNo = G4int(nx + fNoVoxelsX*ny + fNoVox << 
384   }                                               372   }
385                                                   373 
                                                   >> 374   // CheckCopyNo( copyNo ); // not needed, just for debugging code
                                                   >> 375 
386   return copyNo;                                  376   return copyNo;
387 }                                                 377 }
388                                                   378 
389                                                   379 
390 //--------------------------------------------    380 //------------------------------------------------------------------
391 void G4PhantomParameterisation::CheckCopyNo( c << 381 void G4PhantomParameterisation::CheckCopyNo( const G4int copyNo ) const
392 {                                                 382 { 
393   if( copyNo < 0 || copyNo >= G4int(fNoVoxels) << 383   if( copyNo < 0 || copyNo >= G4int(fNoVoxel) )
394   {                                               384   {
395     std::ostringstream message;                << 385     G4cerr << "ERROR - G4PhantomParameterisation::CheckCopyNo()" << G4endl
396     message << "Copy number is negative or too << 386            << "        Copy number: " << copyNo << G4endl
397             << "        Copy number: " << copy << 387            << "        Total number of voxels: " << fNoVoxel << G4endl;
398             << "        Total number of voxels << 
399     G4Exception("G4PhantomParameterisation::Ch    388     G4Exception("G4PhantomParameterisation::CheckCopyNo()",
400                 "GeomNav0002", FatalErrorInArg << 389                 "Wrong-copy-number", FatalErrorInArgument,
                                                   >> 390                 "Copy number is negative or too big!");
401   }                                               391   }
402 }                                                 392 }
403                                                   393