source: trunk/source/geometry/navigation/src/G4PhantomParameterisation.cc@ 1351

Last change on this file since 1351 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
  • Property svn:executable set to *
File size: 13.1 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4PhantomParameterisation.cc,v 1.7 2010/07/13 15:59:42 gcosmo Exp $
28// GEANT4 tag $ Name:$
29//
30// class G4PhantomParameterisation implementation
31//
32// May 2007 Pedro Arce, first version
33//
34// --------------------------------------------------------------------
35
36#include "G4PhantomParameterisation.hh"
37
38#include "globals.hh"
39#include "G4VSolid.hh"
40#include "G4VPhysicalVolume.hh"
41#include "G4LogicalVolume.hh"
42#include "G4VVolumeMaterialScanner.hh"
43#include "G4GeometryTolerance.hh"
44
45//------------------------------------------------------------------
46G4PhantomParameterisation::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)
52{
53 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
54}
55
56
57//------------------------------------------------------------------
58G4PhantomParameterisation::~G4PhantomParameterisation()
59{
60}
61
62
63//------------------------------------------------------------------
64void G4PhantomParameterisation::
65BuildContainerSolid( G4VPhysicalVolume *pMotherPhysical )
66{
67 fContainerSolid = pMotherPhysical->GetLogicalVolume()->GetSolid();
68 fContainerWallX = fNoVoxelX * fVoxelHalfX;
69 fContainerWallY = fNoVoxelY * fVoxelHalfY;
70 fContainerWallZ = fNoVoxelZ * fVoxelHalfZ;
71
72 // CheckVoxelsFillContainer();
73}
74
75//------------------------------------------------------------------
76void G4PhantomParameterisation::
77BuildContainerSolid( G4VSolid *pMotherSolid )
78{
79 fContainerSolid = pMotherSolid;
80 fContainerWallX = fNoVoxelX * fVoxelHalfX;
81 fContainerWallY = fNoVoxelY * fVoxelHalfY;
82 fContainerWallZ = fNoVoxelZ * fVoxelHalfZ;
83
84 // CheckVoxelsFillContainer();
85}
86
87
88//------------------------------------------------------------------
89void G4PhantomParameterisation::
90ComputeTransformation(const G4int copyNo, G4VPhysicalVolume *physVol ) const
91{
92 // Voxels cannot be rotated, return translation
93 //
94 G4ThreeVector trans = GetTranslation( copyNo );
95
96 physVol->SetTranslation( trans );
97}
98
99
100//------------------------------------------------------------------
101G4ThreeVector G4PhantomParameterisation::
102GetTranslation(const G4int copyNo ) const
103{
104 CheckCopyNo( copyNo );
105
106 size_t nx;
107 size_t ny;
108 size_t nz;
109
110 ComputeVoxelIndices( copyNo, nx, ny, nz );
111
112 G4ThreeVector trans( (2*nx+1)*fVoxelHalfX - fContainerWallX,
113 (2*ny+1)*fVoxelHalfY - fContainerWallY,
114 (2*nz+1)*fVoxelHalfZ - fContainerWallZ);
115 return trans;
116}
117
118
119//------------------------------------------------------------------
120G4VSolid* G4PhantomParameterisation::
121ComputeSolid(const G4int, G4VPhysicalVolume *pPhysicalVol)
122{
123 return pPhysicalVol->GetLogicalVolume()->GetSolid();
124}
125
126
127//------------------------------------------------------------------
128G4Material* G4PhantomParameterisation::
129ComputeMaterial(const G4int copyNo, G4VPhysicalVolume *, const G4VTouchable *)
130{
131 CheckCopyNo( copyNo );
132 size_t matIndex = GetMaterialIndex(copyNo);
133
134 return fMaterials[ matIndex ];
135}
136
137
138//------------------------------------------------------------------
139size_t G4PhantomParameterisation::
140GetMaterialIndex( size_t copyNo ) const
141{
142 CheckCopyNo( copyNo );
143
144 if( !fMaterialIndices ) { return 0; }
145 return *(fMaterialIndices+copyNo);
146}
147
148
149//------------------------------------------------------------------
150size_t G4PhantomParameterisation::
151GetMaterialIndex( size_t nx, size_t ny, size_t nz ) const
152{
153 size_t copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
154 return GetMaterialIndex( copyNo );
155}
156
157
158//------------------------------------------------------------------
159G4Material* G4PhantomParameterisation::GetMaterial( size_t nx, size_t ny, size_t nz) const
160{
161 return fMaterials[GetMaterialIndex(nx,ny,nz)];
162}
163
164//------------------------------------------------------------------
165G4Material* G4PhantomParameterisation::GetMaterial( size_t copyNo ) const
166{
167 return fMaterials[GetMaterialIndex(copyNo)];
168}
169
170//------------------------------------------------------------------
171void G4PhantomParameterisation::
172ComputeVoxelIndices(const G4int copyNo, size_t& nx,
173 size_t& ny, size_t& nz ) const
174{
175 CheckCopyNo( copyNo );
176 nx = size_t(copyNo%fNoVoxelX);
177 ny = size_t( (copyNo/fNoVoxelX)%fNoVoxelY );
178 nz = size_t(copyNo/fNoVoxelXY);
179}
180
181
182//------------------------------------------------------------------
183void G4PhantomParameterisation::
184CheckVoxelsFillContainer( G4double contX, G4double contY, G4double contZ ) const
185{
186 G4double toleranceForWarning = 0.25*kCarTolerance;
187
188 // Any bigger value than 0.25*kCarTolerance will give a warning in
189 // G4NormalNavigation::ComputeStep(), because the Inverse of a container
190 // translation that is Z+epsilon gives -Z+epsilon (and the maximum tolerance
191 // in G4Box::Inside is 0.5*kCarTolerance
192 //
193 G4double toleranceForError = 1.*kCarTolerance;
194
195 // Any bigger value than kCarTolerance will give an error in GetReplicaNo()
196 //
197 if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForError
198 || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForError
199 || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForError )
200 {
201 G4cerr << "ERROR - G4PhantomParameterisation::CheckVoxelsFillContainer()"
202 << G4endl
203 << " Voxels do not fully fill the container: "
204 << fContainerSolid->GetName() << G4endl
205 << " DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
206 << " DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
207 << " DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
208 << " Maximum difference is: " << toleranceForError << G4endl;
209 G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
210 "InvalidSetup", FatalException,
211 "Voxels do not fully fill the container!");
212
213 }
214 else if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForWarning
215 || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForWarning
216 || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForWarning )
217 {
218 G4cerr << "WARNING - G4PhantomParameterisation::CheckVoxelsFillContainer()"
219 << G4endl
220 << " Voxels do not fully fill the container: "
221 << fContainerSolid->GetName() << G4endl
222 << " DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
223 << " DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
224 << " DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
225 << " Maximum difference is: " << toleranceForWarning
226 << G4endl;
227 G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
228 "InvalidSetup", JustWarning,
229 "Voxels do not fully fill the container!");
230 }
231}
232
233
234//------------------------------------------------------------------
235G4int G4PhantomParameterisation::
236GetReplicaNo( const G4ThreeVector& localPoint, const G4ThreeVector& localDir )
237{
238
239 // Check first that point is really inside voxels
240 //
241 if( fContainerSolid->Inside( localPoint ) == kOutside )
242 {
243 G4cerr << "ERROR - G4PhantomParameterisation::GetReplicaNo()" << G4endl
244 << " localPoint - " << localPoint
245 << " - is outside container solid: "
246 << fContainerSolid->GetName() << G4endl;
247 G4Exception("G4PhantomParameterisation::GetReplicaNo()", "InvalidSetup",
248 FatalErrorInArgument, "Point outside voxels!");
249 }
250
251 // Check the voxel numbers corresponding to localPoint
252 // When a particle is on a surface, it may be between -kCarTolerance and
253 // +kCartolerance. By a simple distance as:
254 // G4int nx = G4int( (localPoint.x()+)/fVoxelHalfX/2.);
255 // those between -kCartolerance and 0 will be placed on voxel N-1 and those
256 // between 0 and kCarTolerance on voxel N.
257 // To avoid precision problems place the tracks that are on the surface on
258 // voxel N-1 if they have negative direction and on voxel N if they have
259 // positive direction.
260 // Add +kCarTolerance so that they are first placed on voxel N, and then
261 // if the direction is negative substract 1
262
263 G4double fx = (localPoint.x()+fContainerWallX+kCarTolerance)/(fVoxelHalfX*2.);
264 G4int nx = G4int(fx);
265
266 G4double fy = (localPoint.y()+fContainerWallY+kCarTolerance)/(fVoxelHalfY*2.);
267 G4int ny = G4int(fy);
268
269 G4double fz = (localPoint.z()+fContainerWallZ+kCarTolerance)/(fVoxelHalfZ*2.);
270 G4int nz = G4int(fz);
271
272 // If it is on the surface side, check the direction: if direction is
273 // negative place it on the previous voxel (if direction is positive it is
274 // already in the next voxel...).
275 // Correct also cases where n = -1 or n = fNoVoxel. It is always traced to be
276 // due to multiple scattering: track is entering a voxel but multiple
277 // scattering changes the angle towards outside
278 //
279 if( fx - nx < kCarTolerance/fVoxelHalfX )
280 {
281 if( localDir.x() < 0 )
282 {
283 if( nx != 0 )
284 {
285 nx -= 1;
286 }
287 }
288 else
289 {
290 if( nx == G4int(fNoVoxelX) )
291 {
292 nx -= 1;
293 }
294 }
295 }
296 if( fy - ny < kCarTolerance/fVoxelHalfY )
297 {
298 if( localDir.y() < 0 )
299 {
300 if( ny != 0 )
301 {
302 ny -= 1;
303 }
304 }
305 else
306 {
307 if( ny == G4int(fNoVoxelY) )
308 {
309 ny -= 1;
310 }
311 }
312 }
313 if( fz - nz < kCarTolerance/fVoxelHalfZ )
314 {
315 if( localDir.z() < 0 )
316 {
317 if( nz != 0 )
318 {
319 nz -= 1;
320 }
321 }
322 else
323 {
324 if( nz == G4int(fNoVoxelZ) )
325 {
326 nz -= 1;
327 }
328 }
329 }
330
331 G4int copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
332
333 // Check if there are still errors
334 //
335 G4bool isOK = true;
336 if( nx < 0 )
337 {
338 nx = 0;
339 isOK = false;
340 }
341 else if( nx >= G4int(fNoVoxelX) )
342 {
343 nx = fNoVoxelX-1;
344 isOK = false;
345 }
346 if( ny < 0 )
347 {
348 ny = 0;
349 isOK = false;
350 }
351 else if( ny >= G4int(fNoVoxelY) )
352 {
353 ny = fNoVoxelY-1;
354 isOK = false;
355 }
356 if( nz < 0 )
357 {
358 nz = 0;
359 isOK = false;
360 }
361 else if( nz >= G4int(fNoVoxelZ) )
362 {
363 nz = fNoVoxelZ-1;
364 isOK = false;
365 }
366 if( !isOK )
367 {
368 G4cerr << "WARNING - G4PhantomParameterisation::GetReplicaNo()" << G4endl
369 << " LocalPoint: " << localPoint << G4endl
370 << " LocalDir: " << localDir << G4endl
371 << " Voxel container size: " << fContainerWallX
372 << " " << fContainerWallY << " " << fContainerWallZ << G4endl
373 << " LocalPoint - wall: "
374 << localPoint.x()-fContainerWallX << " "
375 << localPoint.y()-fContainerWallY << " "
376 << localPoint.z()-fContainerWallZ << G4endl;
377 G4Exception("G4PhantomParameterisation::GetReplicaNo()",
378 "Wrong-copy-number", JustWarning,
379 "Corrected the copy number! It was negative or too big");
380 copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
381 }
382
383 // CheckCopyNo( copyNo ); // not needed, just for debugging code
384
385 return copyNo;
386}
387
388
389//------------------------------------------------------------------
390void G4PhantomParameterisation::CheckCopyNo( const G4int copyNo ) const
391{
392 if( copyNo < 0 || copyNo >= G4int(fNoVoxel) )
393 {
394 G4cerr << "ERROR - G4PhantomParameterisation::CheckCopyNo()" << G4endl
395 << " Copy number: " << copyNo << G4endl
396 << " Total number of voxels: " << fNoVoxel << G4endl;
397 G4Exception("G4PhantomParameterisation::CheckCopyNo()",
398 "Wrong-copy-number", FatalErrorInArgument,
399 "Copy number is negative or too big!");
400 }
401}
Note: See TracBrowser for help on using the repository browser.