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

Last change on this file since 1191 was 1058, checked in by garnier, 17 years ago

file release beta
  • Property svn:executable set to *
File size: 13.0 KB
RevLine 
[831]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//
[1058]27// $Id: G4PhantomParameterisation.cc,v 1.6 2009/05/20 08:27:10 gcosmo Exp $
[831]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{
48 // Initialise data
49 //
50 fMaterialIndices = 0;
51 fContainerWallX = 0.;
52 fContainerWallY = 0.;
53 fContainerWallZ = 0.;
54
55 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
56
57 bSkipEqualMaterials = 1;
58}
59
60
61//------------------------------------------------------------------
62G4PhantomParameterisation::~G4PhantomParameterisation()
63{
64}
65
66
67//------------------------------------------------------------------
68void G4PhantomParameterisation::
69BuildContainerSolid( G4VPhysicalVolume *pMotherPhysical )
70{
71 fContainerSolid = pMotherPhysical->GetLogicalVolume()->GetSolid();
72 fContainerWallX = fNoVoxelX * fVoxelHalfX;
73 fContainerWallY = fNoVoxelY * fVoxelHalfY;
74 fContainerWallZ = fNoVoxelZ * fVoxelHalfZ;
75
76 // CheckVoxelsFillContainer();
77}
78
[1058]79//------------------------------------------------------------------
80void G4PhantomParameterisation::
81BuildContainerSolid( G4VSolid *pMotherSolid )
82{
83 fContainerSolid = pMotherSolid;
84 fContainerWallX = fNoVoxelX * fVoxelHalfX;
85 fContainerWallY = fNoVoxelY * fVoxelHalfY;
86 fContainerWallZ = fNoVoxelZ * fVoxelHalfZ;
[831]87
[1058]88 // CheckVoxelsFillContainer();
89}
90
91
[831]92//------------------------------------------------------------------
93void G4PhantomParameterisation::
94ComputeTransformation(const G4int copyNo, G4VPhysicalVolume *physVol ) const
95{
96 // Voxels cannot be rotated, return translation
97 //
98 G4ThreeVector trans = GetTranslation( copyNo );
99
100 physVol->SetTranslation( trans );
101}
102
103
104//------------------------------------------------------------------
105G4ThreeVector G4PhantomParameterisation::
106GetTranslation(const G4int copyNo ) const
107{
108 CheckCopyNo( copyNo );
109
110 size_t nx;
111 size_t ny;
112 size_t nz;
113
114 ComputeVoxelIndices( copyNo, nx, ny, nz );
115
116 G4ThreeVector trans( (2*nx+1)*fVoxelHalfX - fContainerWallX,
117 (2*ny+1)*fVoxelHalfY - fContainerWallY,
118 (2*nz+1)*fVoxelHalfZ - fContainerWallZ);
119 return trans;
120}
121
122
123//------------------------------------------------------------------
124G4VSolid* G4PhantomParameterisation::
125ComputeSolid(const G4int, G4VPhysicalVolume *pPhysicalVol)
126{
127 return pPhysicalVol->GetLogicalVolume()->GetSolid();
128}
129
130
131//------------------------------------------------------------------
132G4Material* G4PhantomParameterisation::
133ComputeMaterial(const G4int copyNo, G4VPhysicalVolume *, const G4VTouchable *)
134{
135 CheckCopyNo( copyNo );
136 size_t matIndex = GetMaterialIndex(copyNo);
137
138 return fMaterials[ matIndex ];
139}
140
141
142//------------------------------------------------------------------
143size_t G4PhantomParameterisation::
144GetMaterialIndex( size_t copyNo ) const
145{
146 CheckCopyNo( copyNo );
147
[1058]148 if( !fMaterialIndices ) { return 0; }
[831]149 return *(fMaterialIndices+copyNo);
150}
151
152
153//------------------------------------------------------------------
154size_t G4PhantomParameterisation::
155GetMaterialIndex( size_t nx, size_t ny, size_t nz ) const
156{
157 size_t copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
158 return GetMaterialIndex( copyNo );
159}
160
161
162//------------------------------------------------------------------
163G4Material* G4PhantomParameterisation::GetMaterial( size_t nx, size_t ny, size_t nz) const
164{
165 return fMaterials[GetMaterialIndex(nx,ny,nz)];
166}
167
168//------------------------------------------------------------------
169G4Material* G4PhantomParameterisation::GetMaterial( size_t copyNo ) const
170{
171 return fMaterials[GetMaterialIndex(copyNo)];
172}
173
174//------------------------------------------------------------------
175void G4PhantomParameterisation::
176ComputeVoxelIndices(const G4int copyNo, size_t& nx,
177 size_t& ny, size_t& nz ) const
178{
179 CheckCopyNo( copyNo );
180 nx = size_t(copyNo%fNoVoxelX);
181 ny = size_t( (copyNo/fNoVoxelX)%fNoVoxelY );
182 nz = size_t(copyNo/fNoVoxelXY);
183}
184
185
186//------------------------------------------------------------------
187void G4PhantomParameterisation::
188CheckVoxelsFillContainer( G4double contX, G4double contY, G4double contZ ) const
189{
190 G4double toleranceForWarning = 0.25*kCarTolerance;
191
192 // Any bigger value than 0.25*kCarTolerance will give a warning in
193 // G4NormalNavigation::ComputeStep(), because the Inverse of a container
194 // translation that is Z+epsilon gives -Z+epsilon (and the maximum tolerance
195 // in G4Box::Inside is 0.5*kCarTolerance
196 //
197 G4double toleranceForError = 1.*kCarTolerance;
198
199 // Any bigger value than kCarTolerance will give an error in GetReplicaNo()
200 //
201 if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForError
202 || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForError
203 || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForError )
204 {
205 G4cerr << "ERROR - G4PhantomParameterisation::CheckVoxelsFillContainer()"
206 << G4endl
207 << " Voxels do not fully fill the container: "
208 << fContainerSolid->GetName() << G4endl
209 << " DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
210 << " DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
211 << " DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
212 << " Maximum difference is: " << toleranceForError << G4endl;
213 G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
214 "InvalidSetup", FatalException,
215 "Voxels do not fully fill the container!");
216
217 }
218 else if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForWarning
219 || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForWarning
220 || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForWarning )
221 {
222 G4cerr << "WARNING - G4PhantomParameterisation::CheckVoxelsFillContainer()"
223 << G4endl
224 << " Voxels do not fully fill the container: "
225 << fContainerSolid->GetName() << G4endl
226 << " DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
227 << " DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
228 << " DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
229 << " Maximum difference is: " << toleranceForWarning
230 << G4endl;
231 G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
232 "InvalidSetup", JustWarning,
233 "Voxels do not fully fill the container!");
234 }
235}
236
237
238//------------------------------------------------------------------
239G4int G4PhantomParameterisation::
240GetReplicaNo( const G4ThreeVector& localPoint, const G4ThreeVector& localDir )
241{
242
243 // Check first that point is really inside voxels
244 //
245 if( fContainerSolid->Inside( localPoint ) == kOutside )
246 {
247 G4cerr << "ERROR - G4PhantomParameterisation::GetReplicaNo()" << G4endl
248 << " localPoint - " << localPoint
249 << " - is outside container solid: "
250 << fContainerSolid->GetName() << G4endl;
251 G4Exception("G4PhantomParameterisation::GetReplicaNo()", "InvalidSetup",
252 FatalErrorInArgument, "Point outside voxels!");
253 }
254
255 // Check the voxel numbers corresponding to localPoint
256 // When a particle is on a surface, it may be between -kCarTolerance and
257 // +kCartolerance. By a simple distance as:
258 // G4int nx = G4int( (localPoint.x()+)/fVoxelHalfX/2.);
259 // those between -kCartolerance and 0 will be placed on voxel N-1 and those
260 // between 0 and kCarTolerance on voxel N.
261 // To avoid precision problems place the tracks that are on the surface on
262 // voxel N-1 if they have negative direction and on voxel N if they have
263 // positive direction.
264 // Add +kCarTolerance so that they are first placed on voxel N, and then
265 // if the direction is negative substract 1
266
267 G4double fx = (localPoint.x()+fContainerWallX+kCarTolerance)/(fVoxelHalfX*2.);
268 G4int nx = G4int(fx);
269
270 G4double fy = (localPoint.y()+fContainerWallY+kCarTolerance)/(fVoxelHalfY*2.);
271 G4int ny = G4int(fy);
272
273 G4double fz = (localPoint.z()+fContainerWallZ+kCarTolerance)/(fVoxelHalfZ*2.);
274 G4int nz = G4int(fz);
275
276 // If it is on the surface side, check the direction: if direction is
277 // negative place it on the previous voxel (if direction is positive it is
278 // already in the next voxel...).
279 // Correct also cases where n = -1 or n = fNoVoxel. It is always traced to be
280 // due to multiple scattering: track is entering a voxel but multiple
281 // scattering changes the angle towards outside
282 //
283 if( fx - nx < kCarTolerance/fVoxelHalfX )
284 {
285 if( localDir.x() < 0 )
286 {
287 if( nx != 0 )
288 {
289 nx -= 1;
290 }
291 }
292 else
293 {
294 if( nx == G4int(fNoVoxelX) )
295 {
296 nx -= 1;
297 }
298 }
299 }
300 if( fy - ny < kCarTolerance/fVoxelHalfY )
301 {
302 if( localDir.y() < 0 )
303 {
304 if( ny != 0 )
305 {
306 ny -= 1;
307 }
308 }
309 else
310 {
311 if( ny == G4int(fNoVoxelY) )
312 {
313 ny -= 1;
314 }
315 }
316 }
317 if( fz - nz < kCarTolerance/fVoxelHalfZ )
318 {
319 if( localDir.z() < 0 )
320 {
321 if( nz != 0 )
322 {
323 nz -= 1;
324 }
325 }
326 else
327 {
328 if( nz == G4int(fNoVoxelZ) )
329 {
330 nz -= 1;
331 }
332 }
333 }
334
335 G4int copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
336
337 // Check if there are still errors
338 //
339 G4bool isOK = true;
340 if( nx < 0 )
341 {
342 nx = 0;
343 isOK = false;
344 }
345 else if( nx >= G4int(fNoVoxelX) )
346 {
347 nx = fNoVoxelX-1;
348 isOK = false;
349 }
350 if( ny < 0 )
351 {
352 ny = 0;
353 isOK = false;
354 }
355 else if( ny >= G4int(fNoVoxelY) )
356 {
357 ny = fNoVoxelY-1;
358 isOK = false;
359 }
360 if( nz < 0 )
361 {
362 nz = 0;
363 isOK = false;
364 }
365 else if( nz >= G4int(fNoVoxelZ) )
366 {
367 nz = fNoVoxelZ-1;
368 isOK = false;
369 }
370 if( !isOK )
371 {
372 G4cerr << "WARNING - G4PhantomParameterisation::GetReplicaNo()" << G4endl
373 << " LocalPoint: " << localPoint << G4endl
374 << " LocalDir: " << localDir << G4endl
375 << " Voxel container size: " << fContainerWallX
376 << " " << fContainerWallY << " " << fContainerWallZ << G4endl
377 << " LocalPoint - wall: "
378 << localPoint.x()-fContainerWallX << " "
379 << localPoint.y()-fContainerWallY << " "
380 << localPoint.z()-fContainerWallZ << G4endl;
381 G4Exception("G4PhantomParameterisation::GetReplicaNo()",
382 "Wrong-copy-number", JustWarning,
383 "Corrected the copy number! It was negative or too big");
384 copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
385 }
386
387 // CheckCopyNo( copyNo ); // not needed, just for debugging code
388
389 return copyNo;
390}
391
392
393//------------------------------------------------------------------
394void G4PhantomParameterisation::CheckCopyNo( const G4int copyNo ) const
395{
396 if( copyNo < 0 || copyNo >= G4int(fNoVoxel) )
397 {
398 G4cerr << "ERROR - G4PhantomParameterisation::CheckCopyNo()" << G4endl
399 << " Copy number: " << copyNo << G4endl
400 << " Total number of voxels: " << fNoVoxel << G4endl;
401 G4Exception("G4PhantomParameterisation::CheckCopyNo()",
402 "Wrong-copy-number", FatalErrorInArgument,
403 "Copy number is negative or too big!");
404 }
405}
Note: See TracBrowser for help on using the repository browser.