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