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