Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PhantomParameterisation.cc
Go to the documentation of this file.
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 101807 2016-11-30 13:42:28Z gunter $
28 // GEANT4 tag $ Name:$
29 //
30 // class G4PhantomParameterisation implementation
31 //
32 // May 2007 Pedro Arce, first version
33 //
34 // --------------------------------------------------------------------
35 
37 
38 #include "globals.hh"
39 #include "G4VSolid.hh"
40 #include "G4VPhysicalVolume.hh"
41 #include "G4LogicalVolume.hh"
43 #include "G4GeometryTolerance.hh"
44 
45 //------------------------------------------------------------------
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 {
54 }
55 
56 
57 //------------------------------------------------------------------
59 {
60 }
61 
62 
63 //------------------------------------------------------------------
66 {
67  fContainerSolid = pMotherPhysical->GetLogicalVolume()->GetSolid();
71 
72  // CheckVoxelsFillContainer();
73 }
74 
75 //------------------------------------------------------------------
77 BuildContainerSolid( G4VSolid *pMotherSolid )
78 {
79  fContainerSolid = pMotherSolid;
83 
84  // CheckVoxelsFillContainer();
85 }
86 
87 
88 //------------------------------------------------------------------
90 ComputeTransformation(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 //------------------------------------------------------------------
102 GetTranslation(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 //------------------------------------------------------------------
121 ComputeSolid(const G4int, G4VPhysicalVolume *pPhysicalVol)
122 {
123  return pPhysicalVol->GetLogicalVolume()->GetSolid();
124 }
125 
126 
127 //------------------------------------------------------------------
130 {
131  CheckCopyNo( copyNo );
132  size_t matIndex = GetMaterialIndex(copyNo);
133 
134  return fMaterials[ matIndex ];
135 }
136 
137 
138 //------------------------------------------------------------------
140 GetMaterialIndex( size_t copyNo ) const
141 {
142  CheckCopyNo( copyNo );
143 
144  if( !fMaterialIndices ) { return 0; }
145  return *(fMaterialIndices+copyNo);
146 }
147 
148 
149 //------------------------------------------------------------------
151 GetMaterialIndex( 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 //------------------------------------------------------------------
159 G4Material*
160 G4PhantomParameterisation::GetMaterial( size_t nx, size_t ny, size_t nz) const
161 {
162  return fMaterials[GetMaterialIndex(nx,ny,nz)];
163 }
164 
165 //------------------------------------------------------------------
167 {
168  return fMaterials[GetMaterialIndex(copyNo)];
169 }
170 
171 //------------------------------------------------------------------
172 void G4PhantomParameterisation::
173 ComputeVoxelIndices(const G4int copyNo, size_t& nx,
174  size_t& ny, size_t& nz ) const
175 {
176  CheckCopyNo( copyNo );
177  nx = size_t(copyNo%fNoVoxelX);
178  ny = size_t( (copyNo/fNoVoxelX)%fNoVoxelY );
179  nz = size_t(copyNo/fNoVoxelXY);
180 }
181 
182 
183 //------------------------------------------------------------------
186 {
187  G4double toleranceForWarning = 0.25*kCarTolerance;
188 
189  // Any bigger value than 0.25*kCarTolerance will give a warning in
190  // G4NormalNavigation::ComputeStep(), because the Inverse of a container
191  // translation that is Z+epsilon gives -Z+epsilon (and the maximum tolerance
192  // in G4Box::Inside is 0.5*kCarTolerance
193  //
194  G4double toleranceForError = 1.*kCarTolerance;
195 
196  // Any bigger value than kCarTolerance will give an error in GetReplicaNo()
197  //
198  if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForError
199  || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForError
200  || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForError )
201  {
202  std::ostringstream message;
203  message << "Voxels do not fully fill the container: "
205  << " DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
206  << " DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
207  << " DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
208  << " Maximum difference is: " << toleranceForError;
209  G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
210  "GeomNav0002", FatalException, message);
211 
212  }
213  else if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForWarning
214  || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForWarning
215  || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForWarning )
216  {
217  std::ostringstream message;
218  message << "Voxels do not fully fill the container: "
220  << " DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
221  << " DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
222  << " DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
223  << " Maximum difference is: " << toleranceForWarning;
224  G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
225  "GeomNav1002", JustWarning, message);
226  }
227 }
228 
229 
230 //------------------------------------------------------------------
232 GetReplicaNo( const G4ThreeVector& localPoint, const G4ThreeVector& localDir )
233 {
234 
235  // Check first that point is really inside voxels
236  //
237  if( fContainerSolid->Inside( localPoint ) == kOutside )
238  {
239  std::ostringstream message;
240  message << "Point outside voxels!" << G4endl
241  << " localPoint - " << localPoint
242  << " - is outside container solid: "
243  << fContainerSolid->GetName() << G4endl;
244  G4Exception("G4PhantomParameterisation::GetReplicaNo()", "GeomNav0003",
245  JustWarning, message);
246  G4cerr << " DIFFERENCE WITH PHANTOM WALLS X: "
247  << std::fabs(localPoint.x()) - fContainerWallX
248  << " Y: " << std::fabs(localPoint.y()) - fContainerWallY
249  << " Z: " << std::fabs(localPoint.z()) - fContainerWallZ << G4endl;
250  }
251 
252  // Check the voxel numbers corresponding to localPoint
253  // When a particle is on a surface, it may be between -kCarTolerance and
254  // +kCartolerance. By a simple distance as:
255  // G4int nx = G4int( (localPoint.x()+)/fVoxelHalfX/2.);
256  // those between -kCartolerance and 0 will be placed on voxel N-1 and those
257  // between 0 and kCarTolerance on voxel N.
258  // To avoid precision problems place the tracks that are on the surface on
259  // voxel N-1 if they have negative direction and on voxel N if they have
260  // positive direction.
261  // Add +kCarTolerance so that they are first placed on voxel N, and then
262  // if the direction is negative substract 1
263 
264  G4double fx = (localPoint.x()+fContainerWallX)/(fVoxelHalfX*2.);
265  G4int nx = G4int(fx);
266 
267  G4double fy = (localPoint.y()+fContainerWallY)/(fVoxelHalfY*2.);
268  G4int ny = G4int(fy);
269 
270  G4double fz = (localPoint.z()+fContainerWallZ)/(fVoxelHalfZ*2.);
271  G4int nz = G4int(fz);
272 
273  // If it is on the surface side, check the direction: if direction is
274  // negative place it in the previous voxel (if direction is positive it is
275  // already in the next voxel).
276  // Correct also cases where n = -1 or n = fNoVoxel. It is always traced to be
277  // due to multiple scattering: track is entering a voxel but multiple
278  // scattering changes the angle towards outside
279  //
280  if( fx - nx < kCarTolerance*fContainerWallX )
281  {
282  if( localDir.x() < 0 )
283  {
284  if( nx != 0 )
285  {
286  nx -= 1;
287  }
288  }
289  else
290  {
291  if( nx == G4int(fNoVoxelX) )
292  {
293  nx -= 1;
294  }
295  }
296  }
297  if( fy - ny < kCarTolerance*fContainerWallY )
298  {
299  if( localDir.y() < 0 )
300  {
301  if( ny != 0 )
302  {
303  ny -= 1;
304  }
305  }
306  else
307  {
308  if( ny == G4int(fNoVoxelY) )
309  {
310  ny -= 1;
311  }
312  }
313  }
314  if( fz - nz < kCarTolerance*fContainerWallZ )
315  {
316  if( localDir.z() < 0 )
317  {
318  if( nz != 0 )
319  {
320  nz -= 1;
321  }
322  }
323  else
324  {
325  if( nz == G4int(fNoVoxelZ) )
326  {
327  nz -= 1;
328  }
329  }
330  }
331 
332  G4int copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
333 
334  // Check if there are still errors
335  //
336  G4bool isOK = true;
337  if( nx < 0 )
338  {
339  nx = 0;
340  isOK = false;
341  }
342  else if( nx >= G4int(fNoVoxelX) )
343  {
344  nx = fNoVoxelX-1;
345  isOK = false;
346  }
347  if( ny < 0 )
348  {
349  ny = 0;
350  isOK = false;
351  }
352  else if( ny >= G4int(fNoVoxelY) )
353  {
354  ny = fNoVoxelY-1;
355  isOK = false;
356  }
357  if( nz < 0 )
358  {
359  nz = 0;
360  isOK = false;
361  }
362  else if( nz >= G4int(fNoVoxelZ) )
363  {
364  nz = fNoVoxelZ-1;
365  isOK = false;
366  }
367  if( !isOK )
368  {
369  std::ostringstream message;
370  message << "Corrected the copy number! It was negative or too big" << G4endl
371  << " LocalPoint: " << localPoint << G4endl
372  << " LocalDir: " << localDir << G4endl
373  << " Voxel container size: " << fContainerWallX
374  << " " << fContainerWallY << " " << fContainerWallZ << G4endl
375  << " LocalPoint - wall: "
376  << localPoint.x()-fContainerWallX << " "
377  << localPoint.y()-fContainerWallY << " "
378  << localPoint.z()-fContainerWallZ;
379  G4Exception("G4PhantomParameterisation::GetReplicaNo()",
380  "GeomNav1002", JustWarning, message);
381  copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
382  }
383 
384  // CheckCopyNo( copyNo ); // not needed, just for debugging code
385  // G4cout << " COPYNO " << copyNo << " " << nx << " " << ny << " " << nz
386  // << G4endl; //GDEB
387 
388  return copyNo;
389 }
390 
391 
392 //------------------------------------------------------------------
393 void G4PhantomParameterisation::CheckCopyNo( const G4int copyNo ) const
394 {
395  if( copyNo < 0 || copyNo >= G4int(fNoVoxel) )
396  {
397  std::ostringstream message;
398  message << "Copy number is negative or too big!" << G4endl
399  << " Copy number: " << copyNo << G4endl
400  << " Total number of voxels: " << fNoVoxel;
401  G4Exception("G4PhantomParameterisation::CheckCopyNo()",
402  "GeomNav0002", FatalErrorInArgument, message);
403  }
404 }
G4String GetName() const
double x() const
G4double GetSurfaceTolerance() const
G4VSolid * GetSolid() const
virtual void ComputeTransformation(const G4int, G4VPhysicalVolume *) const
virtual G4Material * ComputeMaterial(const G4int repNo, G4VPhysicalVolume *currentVol, const G4VTouchable *parentTouch=0)
size_t GetMaterialIndex(size_t nx, size_t ny, size_t nz) const
int G4int
Definition: G4Types.hh:78
virtual G4int GetReplicaNo(const G4ThreeVector &localPoint, const G4ThreeVector &localDir)
void BuildContainerSolid(G4VPhysicalVolume *pPhysicalVol)
double z() const
virtual EInside Inside(const G4ThreeVector &p) const =0
bool G4bool
Definition: G4Types.hh:79
virtual G4VSolid * ComputeSolid(const G4int, G4VPhysicalVolume *)
void SetTranslation(const G4ThreeVector &v)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void CheckVoxelsFillContainer(G4double contX, G4double contY, G4double contZ) const
G4Material * GetMaterial(size_t nx, size_t ny, size_t nz) const
G4LogicalVolume * GetLogicalVolume() const
double y() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
std::vector< G4Material * > fMaterials
static G4GeometryTolerance * GetInstance()
G4GLOB_DLL std::ostream G4cerr
G4ThreeVector GetTranslation(const G4int copyNo) const