| 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 G4PartialPhantomParameterisation impl
27 //
28 // May 2007 Pedro Arce (CIEMAT), first version
29 // -------------------------------------------
30
31 #include "G4PartialPhantomParameterisation.hh"
32
33 #include "globals.hh"
34 #include "G4Material.hh"
35 #include "G4VSolid.hh"
36 #include "G4VPhysicalVolume.hh"
37 #include "G4LogicalVolume.hh"
38 #include "G4VVolumeMaterialScanner.hh"
39 #include "G4GeometryTolerance.hh"
40
41 #include <list>
42
43 //--------------------------------------------
44 void G4PartialPhantomParameterisation::
45 ComputeTransformation( const G4int copyNo, G4V
46 {
47 // Voxels cannot be rotated, return translat
48 //
49 G4ThreeVector trans = GetTranslation( copyNo
50 physVol->SetTranslation( trans );
51 }
52
53
54 //--------------------------------------------
55 G4ThreeVector G4PartialPhantomParameterisation
56 GetTranslation(const G4int copyNo ) const
57 {
58 CheckCopyNo( copyNo );
59
60 std::size_t nx, ny, nz;
61 ComputeVoxelIndices( copyNo, nx, ny, nz );
62
63 G4ThreeVector trans( (2*nx+1)*fVoxelHalfX -
64 (2*ny+1)*fVoxelHalfY -
65 (2*nz+1)*fVoxelHalfZ -
66 return trans;
67 }
68
69
70 //--------------------------------------------
71 G4Material* G4PartialPhantomParameterisation::
72 ComputeMaterial( const G4int copyNo, G4VPhysic
73 {
74 CheckCopyNo( copyNo );
75 auto matIndex = GetMaterialIndex(copyNo);
76
77 return fMaterials[ matIndex ];
78 }
79
80
81 //--------------------------------------------
82 size_t G4PartialPhantomParameterisation::
83 GetMaterialIndex( std::size_t copyNo ) const
84 {
85 CheckCopyNo( copyNo );
86
87 if( fMaterialIndices == nullptr ) { return 0
88
89 return *(fMaterialIndices+copyNo);
90 }
91
92
93 //--------------------------------------------
94 size_t G4PartialPhantomParameterisation::
95 GetMaterialIndex( std::size_t nx, std::size_t
96 {
97 std::size_t copyNo = nx + fNoVoxelsX*ny + fN
98 return GetMaterialIndex( copyNo );
99 }
100
101
102 //--------------------------------------------
103 G4Material* G4PartialPhantomParameterisation::
104 GetMaterial( std::size_t nx, std::size_t ny, s
105 {
106 return fMaterials[GetMaterialIndex(nx,ny,nz)
107 }
108
109
110 //--------------------------------------------
111 G4Material* G4PartialPhantomParameterisation::
112 GetMaterial( std::size_t copyNo ) const
113 {
114 return fMaterials[GetMaterialIndex(copyNo)];
115 }
116
117
118 //--------------------------------------------
119 void G4PartialPhantomParameterisation::
120 ComputeVoxelIndices(const G4int copyNo, std::s
121 std::size_t& ny, std
122 {
123 CheckCopyNo( copyNo );
124
125 auto ite = fFilledIDs.lower_bound(copyNo);
126 G4long dist = std::distance( fFilledIDs.cbeg
127 nz = std::size_t( dist/fNoVoxelsY );
128 ny = std::size_t( dist%fNoVoxelsY );
129
130 G4int ifmin = (*ite).second;
131 G4int nvoxXprev;
132 if( dist != 0 )
133 {
134 ite--;
135 nvoxXprev = (*ite).first;
136 }
137 else
138 {
139 nvoxXprev = -1;
140 }
141
142 nx = ifmin+copyNo-nvoxXprev-1;
143 }
144
145
146 //--------------------------------------------
147 G4int G4PartialPhantomParameterisation::
148 GetReplicaNo( const G4ThreeVector& localPoint,
149 {
150 // Check the voxel numbers corresponding to
151 // When a particle is on a surface, it may b
152 // +kCartolerance. By a simple distance as:
153 // G4int nx = G4int( (localPoint.x()+)/fVo
154 // those between -kCartolerance and 0 will b
155 // between 0 and kCarTolerance on voxel N.
156 // To avoid precision problems place the tra
157 // voxel N-1 if they have negative direction
158 // positive direction.
159 // Add +kCarTolerance so that they are first
160 // if the direction is negative substract 1
161
162 G4double fx = (localPoint.x()+fContainerWall
163 auto nx = G4int(fx);
164
165 G4double fy = (localPoint.y()+fContainerWall
166 auto ny = G4int(fy);
167
168 G4double fz = (localPoint.z()+fContainerWall
169 auto nz = G4int(fz);
170
171 // If it is on the surface side, check the d
172 // negative place it on the previous voxel (
173 // already in the next voxel...).
174 // Correct also cases where n = -1 or n = fN
175 // due to multiple scattering: track is ente
176 // scattering changes the angle towards outs
177 //
178 if( fx - nx < kCarTolerance/fVoxelHalfX )
179 {
180 if( localDir.x() < 0 )
181 {
182 if( nx != 0 )
183 {
184 nx -= 1;
185 }
186 }
187 else
188 {
189 if( nx == G4int(fNoVoxelsX) )
190 {
191 nx -= 1;
192 }
193 }
194 }
195 if( fy - ny < kCarTolerance/fVoxelHalfY )
196 {
197 if( localDir.y() < 0 )
198 {
199 if( ny != 0 )
200 {
201 ny -= 1;
202 }
203 }
204 else
205 {
206 if( ny == G4int(fNoVoxelsY) )
207 {
208 ny -= 1;
209 }
210 }
211 }
212 if( fz - nz < kCarTolerance/fVoxelHalfZ )
213 {
214 if( localDir.z() < 0 )
215 {
216 if( nz != 0 )
217 {
218 nz -= 1;
219 }
220 }
221 else
222 {
223 if( nz == G4int(fNoVoxelsZ) )
224 {
225 nz -= 1;
226 }
227 }
228 }
229
230 // Check if there are still errors
231 //
232 G4bool isOK = true;
233 if( nx < 0 )
234 {
235 nx = 0;
236 isOK = false;
237 }
238 else if( nx >= G4int(fNoVoxelsX) )
239 {
240 nx = G4int(fNoVoxelsX)-1;
241 isOK = false;
242 }
243 if( ny < 0 )
244 {
245 ny = 0;
246 isOK = false;
247 }
248 else if( ny >= G4int(fNoVoxelsY) )
249 {
250 ny = G4int(fNoVoxelsY)-1;
251 isOK = false;
252 }
253 if( nz < 0 )
254 {
255 nz = 0;
256 isOK = false;
257 }
258 else if( nz >= G4int(fNoVoxelsZ) )
259 {
260 nz = G4int(fNoVoxelsZ)-1;
261 isOK = false;
262 }
263 if( !isOK )
264 {
265 std::ostringstream message;
266 message << "Corrected the copy number! It
267 << G4endl
268 << " LocalPoint: " << loc
269 << " LocalDir: " << local
270 << " Voxel container size
271 << " " << fContainerWallY << " " <
272 << " LocalPoint - wall: "
273 << localPoint.x()-fContainerWallX
274 << localPoint.y()-fContainerWallY
275 << localPoint.z()-fContainerWallZ;
276 G4Exception("G4PartialPhantomParameterisat
277 "GeomNav1002", JustWarning, me
278 }
279
280 auto nyz = G4int(nz*fNoVoxelsY+ny);
281 auto ite = fFilledIDs.cbegin();
282 /*
283 for( ite = fFilledIDs.cbegin(); ite != fFill
284 {
285 G4cout << " G4PartialPhantomParameterisati
286 << (*ite).first << " , " << (*ite).
287 }
288 */
289
290 advance(ite,nyz);
291 auto iteant = ite; iteant--;
292 G4int copyNo = (*iteant).first + 1 + ( nx -
293 /*
294 G4cout << " G4PartialPhantomParameterisation
295 << copyNo << " nyz " << nyz << " (*
296 << (*iteant).first << " (*ite).secon
297
298 G4cout << " G4PartialPhantomParameterisation
299 << " nx " << nx << " ny " << ny << "
300 << " localPoint " << localPoint << "
301 */
302 return copyNo;
303 }
304
305
306 //--------------------------------------------
307 void G4PartialPhantomParameterisation::CheckCo
308 {
309 if( copyNo < 0 || copyNo >= G4int(fNoVoxels)
310 {
311 std::ostringstream message;
312 message << "Copy number is negative or too
313 << " Copy number: " << copy
314 << " Total number of voxels
315 G4Exception("G4PartialPhantomParameterisat
316 "GeomNav0002", FatalErrorInArg
317 }
318 }
319
320
321 //--------------------------------------------
322 void G4PartialPhantomParameterisation::BuildCo
323 {
324 fContainerWallX = fNoVoxelsX * fVoxelHalfX;
325 fContainerWallY = fNoVoxelsY * fVoxelHalfY;
326 fContainerWallZ = fNoVoxelsZ * fVoxelHalfZ;
327 }
328