| Geant4 Cross Reference |
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 102290 2017-01-20 11:19:44Z 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*
150 G4PhantomParameterisation::GetMaterial( std::s << 160 G4PhantomParameterisation::GetMaterial( size_t nx, size_t ny, size_t nz) const
151 { 161 {
152 return fMaterials[GetMaterialIndex(nx,ny,nz) 162 return fMaterials[GetMaterialIndex(nx,ny,nz)];
153 } 163 }
154 164
155 <<
156 //-------------------------------------------- 165 //------------------------------------------------------------------
157 G4Material* G4PhantomParameterisation::GetMate << 166 G4Material* G4PhantomParameterisation::GetMaterial( size_t copyNo ) const
158 { 167 {
159 return fMaterials[GetMaterialIndex(copyNo)]; 168 return fMaterials[GetMaterialIndex(copyNo)];
160 } 169 }
161 170
162 <<
163 //-------------------------------------------- 171 //------------------------------------------------------------------
164 void G4PhantomParameterisation:: 172 void G4PhantomParameterisation::
165 ComputeVoxelIndices(const G4int copyNo, std::s << 173 ComputeVoxelIndices(const G4int copyNo, size_t& nx,
166 std::size_t& ny, std << 174 size_t& ny, size_t& nz ) const
167 { 175 {
168 CheckCopyNo( copyNo ); 176 CheckCopyNo( copyNo );
169 nx = std::size_t(copyNo%fNoVoxelsX); << 177 nx = size_t(copyNo%fNoVoxelX);
170 ny = std::size_t( (copyNo/fNoVoxelsX)%fNoVox << 178 ny = size_t( (copyNo/fNoVoxelX)%fNoVoxelY );
171 nz = std::size_t(copyNo/fNoVoxelsXY); << 179 nz = size_t(copyNo/fNoVoxelXY);
172 } 180 }
173 181
174 182
175 //-------------------------------------------- 183 //------------------------------------------------------------------
176 void G4PhantomParameterisation:: 184 void G4PhantomParameterisation::
177 CheckVoxelsFillContainer( G4double contX, G4do 185 CheckVoxelsFillContainer( G4double contX, G4double contY, G4double contZ ) const
178 { 186 {
179 G4double toleranceForWarning = 0.25*kCarTole 187 G4double toleranceForWarning = 0.25*kCarTolerance;
180 188
181 // Any bigger value than 0.25*kCarTolerance 189 // Any bigger value than 0.25*kCarTolerance will give a warning in
182 // G4NormalNavigation::ComputeStep(), becaus 190 // G4NormalNavigation::ComputeStep(), because the Inverse of a container
183 // translation that is Z+epsilon gives -Z+ep 191 // translation that is Z+epsilon gives -Z+epsilon (and the maximum tolerance
184 // in G4Box::Inside is 0.5*kCarTolerance 192 // in G4Box::Inside is 0.5*kCarTolerance
185 // 193 //
186 G4double toleranceForError = 1.*kCarToleranc 194 G4double toleranceForError = 1.*kCarTolerance;
187 195
188 // Any bigger value than kCarTolerance will 196 // Any bigger value than kCarTolerance will give an error in GetReplicaNo()
189 // 197 //
190 if( std::fabs(contX-fNoVoxelsX*fVoxelHalfX) << 198 if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForError
191 || std::fabs(contY-fNoVoxelsY*fVoxelHalfY) << 199 || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForError
192 || std::fabs(contZ-fNoVoxelsZ*fVoxelHalfZ) << 200 || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForError )
193 { 201 {
194 std::ostringstream message; 202 std::ostringstream message;
195 message << "Voxels do not fully fill the c 203 message << "Voxels do not fully fill the container: "
196 << fContainerSolid->GetName() << G 204 << fContainerSolid->GetName() << G4endl
197 << " DiffX= " << contX-fNoV << 205 << " DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
198 << " DiffY= " << contY-fNoV << 206 << " DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
199 << " DiffZ= " << contZ-fNoV << 207 << " DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
200 << " Maximum difference is: 208 << " Maximum difference is: " << toleranceForError;
201 G4Exception("G4PhantomParameterisation::Ch 209 G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
202 "GeomNav0002", FatalException, 210 "GeomNav0002", FatalException, message);
203 211
204 } 212 }
205 else if( std::fabs(contX-fNoVoxelsX*fVoxelHa << 213 else if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForWarning
206 || std::fabs(contY-fNoVoxelsY*fVoxelHa << 214 || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForWarning
207 || std::fabs(contZ-fNoVoxelsZ*fVoxelHa << 215 || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForWarning )
208 { 216 {
209 std::ostringstream message; 217 std::ostringstream message;
210 message << "Voxels do not fully fill the c 218 message << "Voxels do not fully fill the container: "
211 << fContainerSolid->GetName() << G 219 << fContainerSolid->GetName() << G4endl
212 << " DiffX= " << contX-fN << 220 << " DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
213 << " DiffY= " << contY-fN << 221 << " DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
214 << " DiffZ= " << contZ-fN << 222 << " DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
215 << " Maximum difference i 223 << " Maximum difference is: " << toleranceForWarning;
216 G4Exception("G4PhantomParameterisation::Ch 224 G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
217 "GeomNav1002", JustWarning, me 225 "GeomNav1002", JustWarning, message);
218 } 226 }
219 } 227 }
220 228
221 229
222 //-------------------------------------------- 230 //------------------------------------------------------------------
223 G4int G4PhantomParameterisation:: 231 G4int G4PhantomParameterisation::
224 GetReplicaNo( const G4ThreeVector& localPoint, 232 GetReplicaNo( const G4ThreeVector& localPoint, const G4ThreeVector& localDir )
225 { 233 {
226 234
227 // Check first that point is really inside v 235 // Check first that point is really inside voxels
228 // 236 //
229 if( fContainerSolid->Inside( localPoint ) == 237 if( fContainerSolid->Inside( localPoint ) == kOutside )
230 { 238 {
231 if( std::fabs(localPoint.x()) - fContainer << 239 std::ostringstream message;
232 && std::fabs(localPoint.y()) - fContainerWal << 240 message << "Point outside voxels!" << G4endl
233 && std::fabs(localPoint.z()) - fContainerWal << 241 << " localPoint - " << localPoint
234 { << 242 << " - is outside container solid: "
235 std::ostringstream message; << 243 << fContainerSolid->GetName() << G4endl;
236 message << "Point outside voxels!" << G4 << 244 G4Exception("G4PhantomParameterisation::GetReplicaNo()", "GeomNav0003",
237 << " localPoint - " << localPoi << 245 JustWarning, message);
238 << " - is outside container solid: " << 246 G4cerr << " DIFFERENCE WITH PHANTOM WALLS X: "
239 << fContainerSolid->GetName() << G4end << 247 << fabs(localPoint.x()) - fContainerWallX
240 << "DIFFERENCE WITH PHANTOM WALLS X: " << 248 << " Y: " << fabs(localPoint.y()) - fContainerWallY
241 << std::fabs(localPoint.x()) - fContai << 249 << " Z: " << fabs(localPoint.z()) - fContainerWallZ << G4endl;
242 << " Y: " << std::fabs(localPoint.y()) <<
243 << " Z: " << std::fabs(localPoint.z()) <<
244 G4Exception("G4PhantomParameterisation:: <<
245 FatalErrorInArgument, message); <<
246 } <<
247 } 250 }
248 251
249 // Check the voxel numbers corresponding to 252 // Check the voxel numbers corresponding to localPoint
250 // When a particle is on a surface, it may b 253 // When a particle is on a surface, it may be between -kCarTolerance and
251 // +kCartolerance. By a simple distance as: 254 // +kCartolerance. By a simple distance as:
252 // G4int nx = G4int( (localPoint.x()+)/fVo 255 // G4int nx = G4int( (localPoint.x()+)/fVoxelHalfX/2.);
253 // those between -kCartolerance and 0 will b 256 // those between -kCartolerance and 0 will be placed on voxel N-1 and those
254 // between 0 and kCarTolerance on voxel N. 257 // between 0 and kCarTolerance on voxel N.
255 // To avoid precision problems place the tra 258 // To avoid precision problems place the tracks that are on the surface on
256 // voxel N-1 if they have negative direction 259 // voxel N-1 if they have negative direction and on voxel N if they have
257 // positive direction. 260 // positive direction.
258 // Add +kCarTolerance so that they are first 261 // Add +kCarTolerance so that they are first placed on voxel N, and then
259 // if the direction is negative substract 1 262 // if the direction is negative substract 1
260 263
261 G4double fx = (localPoint.x()+fContainerWall << 264 G4double fx = (localPoint.x()+fContainerWallX)/(fVoxelHalfX*2.);
262 auto nx = G4int(fx); << 265 G4int nx = G4int(fx);
263 266
264 G4double fy = (localPoint.y()+fContainerWall << 267 G4double fy = (localPoint.y()+fContainerWallY)/(fVoxelHalfY*2.);
265 auto ny = G4int(fy); << 268 G4int ny = G4int(fy);
266 269
267 G4double fz = (localPoint.z()+fContainerWall << 270 G4double fz = (localPoint.z()+fContainerWallZ)/(fVoxelHalfZ*2.);
268 auto nz = G4int(fz); << 271 G4int nz = G4int(fz);
269 272
270 // If it is on the surface side, check the d 273 // If it is on the surface side, check the direction: if direction is
271 // negative place it in the previous voxel ( 274 // negative place it in the previous voxel (if direction is positive it is
272 // already in the next voxel). 275 // already in the next voxel).
273 // Correct also cases where n = -1 or n = fN << 276 // Correct also cases where n = -1 or n = fNoVoxel. It is always traced to be
274 // due to multiple scattering: track is ente 277 // due to multiple scattering: track is entering a voxel but multiple
275 // scattering changes the angle towards outs 278 // scattering changes the angle towards outside
276 // 279 //
277 if( fx - nx < kCarTolerance*fVoxelHalfX ) << 280 if( fx - nx < kCarTolerance*fContainerWallX )
278 { 281 {
279 if( localDir.x() < 0 ) 282 if( localDir.x() < 0 )
280 { 283 {
281 if( nx != 0 ) 284 if( nx != 0 )
282 { 285 {
283 nx -= 1; 286 nx -= 1;
284 } 287 }
285 } 288 }
286 else 289 else
287 { 290 {
288 if( nx == G4int(fNoVoxelsX) ) << 291 if( nx == G4int(fNoVoxelX) )
289 { 292 {
290 nx -= 1; 293 nx -= 1;
291 } 294 }
292 } 295 }
293 } 296 }
294 if( fy - ny < kCarTolerance*fVoxelHalfY ) << 297 if( fy - ny < kCarTolerance*fContainerWallY )
295 { 298 {
296 if( localDir.y() < 0 ) 299 if( localDir.y() < 0 )
297 { 300 {
298 if( ny != 0 ) 301 if( ny != 0 )
299 { 302 {
300 ny -= 1; 303 ny -= 1;
301 } 304 }
302 } 305 }
303 else 306 else
304 { 307 {
305 if( ny == G4int(fNoVoxelsY) ) << 308 if( ny == G4int(fNoVoxelY) )
306 { 309 {
307 ny -= 1; 310 ny -= 1;
308 } 311 }
309 } 312 }
310 } 313 }
311 if( fz - nz < kCarTolerance*fVoxelHalfZ ) << 314 if( fz - nz < kCarTolerance*fContainerWallZ )
312 { 315 {
313 if( localDir.z() < 0 ) 316 if( localDir.z() < 0 )
314 { 317 {
315 if( nz != 0 ) 318 if( nz != 0 )
316 { 319 {
317 nz -= 1; 320 nz -= 1;
318 } 321 }
319 } 322 }
320 else 323 else
321 { 324 {
322 if( nz == G4int(fNoVoxelsZ) ) << 325 if( nz == G4int(fNoVoxelZ) )
323 { 326 {
324 nz -= 1; 327 nz -= 1;
325 } 328 }
326 } 329 }
327 } 330 }
328 331
329 auto copyNo = G4int(nx + fNoVoxelsX*ny + fN << 332 G4int copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
330 333
331 // Check if there are still errors 334 // Check if there are still errors
332 // 335 //
333 G4bool isOK = true; 336 G4bool isOK = true;
334 if( nx < 0 ) 337 if( nx < 0 )
335 { 338 {
336 nx = 0; 339 nx = 0;
337 isOK = false; 340 isOK = false;
338 } 341 }
339 else if( nx >= G4int(fNoVoxelsX) ) << 342 else if( nx >= G4int(fNoVoxelX) )
340 { 343 {
341 nx = G4int(fNoVoxelsX)-1; << 344 nx = fNoVoxelX-1;
342 isOK = false; 345 isOK = false;
343 } 346 }
344 if( ny < 0 ) 347 if( ny < 0 )
345 { 348 {
346 ny = 0; 349 ny = 0;
347 isOK = false; 350 isOK = false;
348 } 351 }
349 else if( ny >= G4int(fNoVoxelsY) ) << 352 else if( ny >= G4int(fNoVoxelY) )
350 { 353 {
351 ny = G4int(fNoVoxelsY)-1; << 354 ny = fNoVoxelY-1;
352 isOK = false; 355 isOK = false;
353 } 356 }
354 if( nz < 0 ) 357 if( nz < 0 )
355 { 358 {
356 nz = 0; 359 nz = 0;
357 isOK = false; 360 isOK = false;
358 } 361 }
359 else if( nz >= G4int(fNoVoxelsZ) ) << 362 else if( nz >= G4int(fNoVoxelZ) )
360 { 363 {
361 nz = G4int(fNoVoxelsZ)-1; << 364 nz = fNoVoxelZ-1;
362 isOK = false; 365 isOK = false;
363 } 366 }
364 if( !isOK ) 367 if( !isOK )
365 { 368 {
366 if( std::fabs(localPoint.x()-fContainerWal << 369 std::ostringstream message;
367 std::fabs(localPoint.y()-fContainerWallY) > << 370 message << "Corrected the copy number! It was negative or too big" << G4endl
368 std::fabs(localPoint.z()-fContainerWallZ) > << 371 << " LocalPoint: " << localPoint << G4endl
369 std::ostringstream message; << 372 << " LocalDir: " << localDir << G4endl
370 message << "Corrected the copy number! I << 373 << " Voxel container size: " << fContainerWallX
371 << " LocalPoint: " << localPo << 374 << " " << fContainerWallY << " " << fContainerWallZ << G4endl
372 << " LocalDir: " << localDir << 375 << " LocalPoint - wall: "
373 << " Voxel container size: " << 376 << localPoint.x()-fContainerWallX << " "
374 << " " << fContainerWallY << " " << fC << 377 << localPoint.y()-fContainerWallY << " "
375 << " LocalPoint - wall: " << 378 << localPoint.z()-fContainerWallZ;
376 << localPoint.x()-fContainerWallX << " << 379 G4Exception("G4PhantomParameterisation::GetReplicaNo()",
377 << localPoint.y()-fContainerWallY << " << 380 "GeomNav1002", JustWarning, message);
378 << localPoint.z()-fContainerWallZ; << 381 copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
379 G4Exception("G4PhantomParameterisation:: <<
380 "GeomNav1002", JustWarning, message); <<
381 } <<
382 <<
383 copyNo = G4int(nx + fNoVoxelsX*ny + fNoVox <<
384 } 382 }
385 383
>> 384 // CheckCopyNo( copyNo ); // not needed, just for debugging code
>> 385 // G4cout << " COPYNO " << copyNo << " " << nx << " " << ny << " " << nz
>> 386 // << G4endl; //GDEB
>> 387
386 return copyNo; 388 return copyNo;
387 } 389 }
388 390
389 391
390 //-------------------------------------------- 392 //------------------------------------------------------------------
391 void G4PhantomParameterisation::CheckCopyNo( c << 393 void G4PhantomParameterisation::CheckCopyNo( const G4int copyNo ) const
392 { 394 {
393 if( copyNo < 0 || copyNo >= G4int(fNoVoxels) << 395 if( copyNo < 0 || copyNo >= G4int(fNoVoxel) )
394 { 396 {
395 std::ostringstream message; 397 std::ostringstream message;
396 message << "Copy number is negative or too 398 message << "Copy number is negative or too big!" << G4endl
397 << " Copy number: " << copy 399 << " Copy number: " << copyNo << G4endl
398 << " Total number of voxels << 400 << " Total number of voxels: " << fNoVoxel;
399 G4Exception("G4PhantomParameterisation::Ch 401 G4Exception("G4PhantomParameterisation::CheckCopyNo()",
400 "GeomNav0002", FatalErrorInArg 402 "GeomNav0002", FatalErrorInArgument, message);
401 } 403 }
402 } 404 }
403 405