Geant4  10.01.p01
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 66356 2012-12-18 09:02:32Z gcosmo $
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* G4PhantomParameterisation::GetMaterial( size_t nx, size_t ny, size_t nz) const
160 {
161  return fMaterials[GetMaterialIndex(nx,ny,nz)];
162 }
163 
164 //------------------------------------------------------------------
166 {
167  return fMaterials[GetMaterialIndex(copyNo)];
168 }
169 
170 //------------------------------------------------------------------
172 ComputeVoxelIndices(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 //------------------------------------------------------------------
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  std::ostringstream message;
202  message << "Voxels do not fully fill the container: "
204  << " DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
205  << " DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
206  << " DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
207  << " Maximum difference is: " << toleranceForError;
208  G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
209  "GeomNav0002", FatalException, message);
210 
211  }
212  else if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForWarning
213  || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForWarning
214  || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForWarning )
215  {
216  std::ostringstream message;
217  message << "Voxels do not fully fill the container: "
219  << " DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
220  << " DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
221  << " DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
222  << " Maximum difference is: " << toleranceForWarning;
223  G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
224  "GeomNav1002", JustWarning, message);
225  }
226 }
227 
228 
229 //------------------------------------------------------------------
231 GetReplicaNo( const G4ThreeVector& localPoint, const G4ThreeVector& localDir )
232 {
233 
234  // Check first that point is really inside voxels
235  //
236  if( fContainerSolid->Inside( localPoint ) == kOutside )
237  {
238  std::ostringstream message;
239  message << "Point outside voxels!" << G4endl
240  << " localPoint - " << localPoint
241  << " - is outside container solid: "
242  << fContainerSolid->GetName() << G4endl;
243  G4Exception("G4PhantomParameterisation::GetReplicaNo()", "GeomNav0003",
244  FatalErrorInArgument, message);
245  }
246 
247  // Check the voxel numbers corresponding to localPoint
248  // When a particle is on a surface, it may be between -kCarTolerance and
249  // +kCartolerance. By a simple distance as:
250  // G4int nx = G4int( (localPoint.x()+)/fVoxelHalfX/2.);
251  // those between -kCartolerance and 0 will be placed on voxel N-1 and those
252  // between 0 and kCarTolerance on voxel N.
253  // To avoid precision problems place the tracks that are on the surface on
254  // voxel N-1 if they have negative direction and on voxel N if they have
255  // positive direction.
256  // Add +kCarTolerance so that they are first placed on voxel N, and then
257  // if the direction is negative substract 1
258 
259  G4double fx = (localPoint.x()+fContainerWallX+kCarTolerance)/(fVoxelHalfX*2.);
260  G4int nx = G4int(fx);
261 
262  G4double fy = (localPoint.y()+fContainerWallY+kCarTolerance)/(fVoxelHalfY*2.);
263  G4int ny = G4int(fy);
264 
265  G4double fz = (localPoint.z()+fContainerWallZ+kCarTolerance)/(fVoxelHalfZ*2.);
266  G4int nz = G4int(fz);
267 
268  // If it is on the surface side, check the direction: if direction is
269  // negative place it on the previous voxel (if direction is positive it is
270  // already in the next voxel...).
271  // Correct also cases where n = -1 or n = fNoVoxel. It is always traced to be
272  // due to multiple scattering: track is entering a voxel but multiple
273  // scattering changes the angle towards outside
274  //
275  if( fx - nx < kCarTolerance/fVoxelHalfX )
276  {
277  if( localDir.x() < 0 )
278  {
279  if( nx != 0 )
280  {
281  nx -= 1;
282  }
283  }
284  else
285  {
286  if( nx == G4int(fNoVoxelX) )
287  {
288  nx -= 1;
289  }
290  }
291  }
292  if( fy - ny < kCarTolerance/fVoxelHalfY )
293  {
294  if( localDir.y() < 0 )
295  {
296  if( ny != 0 )
297  {
298  ny -= 1;
299  }
300  }
301  else
302  {
303  if( ny == G4int(fNoVoxelY) )
304  {
305  ny -= 1;
306  }
307  }
308  }
309  if( fz - nz < kCarTolerance/fVoxelHalfZ )
310  {
311  if( localDir.z() < 0 )
312  {
313  if( nz != 0 )
314  {
315  nz -= 1;
316  }
317  }
318  else
319  {
320  if( nz == G4int(fNoVoxelZ) )
321  {
322  nz -= 1;
323  }
324  }
325  }
326 
327  G4int copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
328 
329  // Check if there are still errors
330  //
331  G4bool isOK = true;
332  if( nx < 0 )
333  {
334  nx = 0;
335  isOK = false;
336  }
337  else if( nx >= G4int(fNoVoxelX) )
338  {
339  nx = fNoVoxelX-1;
340  isOK = false;
341  }
342  if( ny < 0 )
343  {
344  ny = 0;
345  isOK = false;
346  }
347  else if( ny >= G4int(fNoVoxelY) )
348  {
349  ny = fNoVoxelY-1;
350  isOK = false;
351  }
352  if( nz < 0 )
353  {
354  nz = 0;
355  isOK = false;
356  }
357  else if( nz >= G4int(fNoVoxelZ) )
358  {
359  nz = fNoVoxelZ-1;
360  isOK = false;
361  }
362  if( !isOK )
363  {
364  std::ostringstream message;
365  message << "Corrected the copy number! It was negative or too big" << G4endl
366  << " LocalPoint: " << localPoint << G4endl
367  << " LocalDir: " << localDir << G4endl
368  << " Voxel container size: " << fContainerWallX
369  << " " << fContainerWallY << " " << fContainerWallZ << G4endl
370  << " LocalPoint - wall: "
371  << localPoint.x()-fContainerWallX << " "
372  << localPoint.y()-fContainerWallY << " "
373  << localPoint.z()-fContainerWallZ;
374  G4Exception("G4PhantomParameterisation::GetReplicaNo()",
375  "GeomNav1002", JustWarning, message);
376  copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
377  }
378 
379  // CheckCopyNo( copyNo ); // not needed, just for debugging code
380 
381  return copyNo;
382 }
383 
384 
385 //------------------------------------------------------------------
387 {
388  if( copyNo < 0 || copyNo >= G4int(fNoVoxel) )
389  {
390  std::ostringstream message;
391  message << "Copy number is negative or too big!" << G4endl
392  << " Copy number: " << copyNo << G4endl
393  << " Total number of voxels: " << fNoVoxel;
394  G4Exception("G4PhantomParameterisation::CheckCopyNo()",
395  "GeomNav0002", FatalErrorInArgument, message);
396  }
397 }
G4String GetName() const
void CheckCopyNo(const G4int copyNo) const
void ComputeVoxelIndices(const G4int copyNo, size_t &nx, size_t &ny, size_t &nz) const
CLHEP::Hep3Vector G4ThreeVector
G4double GetSurfaceTolerance() 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)
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
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
std::vector< G4Material * > fMaterials
static G4GeometryTolerance * GetInstance()
G4VSolid * GetSolid() const
G4ThreeVector GetTranslation(const G4int copyNo) const