ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PhantomParameterisation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PhantomParameterisation.cc
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 // class G4PhantomParameterisation implementation
27 //
28 // May 2007 Pedro Arce, first version
29 //
30 // --------------------------------------------------------------------
31 
33 
34 #include "globals.hh"
35 #include "G4VSolid.hh"
36 #include "G4VPhysicalVolume.hh"
37 #include "G4LogicalVolume.hh"
39 #include "G4GeometryTolerance.hh"
40 
41 //------------------------------------------------------------------
43 {
45 }
46 
47 
48 //------------------------------------------------------------------
50 {
51 }
52 
53 
54 //------------------------------------------------------------------
57 {
58  fContainerSolid = pMotherPhysical->GetLogicalVolume()->GetSolid();
62 
63  // CheckVoxelsFillContainer();
64 }
65 
66 //------------------------------------------------------------------
68 BuildContainerSolid( G4VSolid* pMotherSolid )
69 {
70  fContainerSolid = pMotherSolid;
74 
75  // CheckVoxelsFillContainer();
76 }
77 
78 
79 //------------------------------------------------------------------
81 ComputeTransformation(const G4int copyNo, G4VPhysicalVolume* physVol ) const
82 {
83  // Voxels cannot be rotated, return translation
84  //
85  G4ThreeVector trans = GetTranslation( copyNo );
86 
87  physVol->SetTranslation( trans );
88 }
89 
90 
91 //------------------------------------------------------------------
93 GetTranslation(const G4int copyNo ) const
94 {
95  CheckCopyNo( copyNo );
96 
97  size_t nx;
98  size_t ny;
99  size_t nz;
100 
101  ComputeVoxelIndices( copyNo, nx, ny, nz );
102 
103  G4ThreeVector trans( (2*nx+1)*fVoxelHalfX - fContainerWallX,
104  (2*ny+1)*fVoxelHalfY - fContainerWallY,
105  (2*nz+1)*fVoxelHalfZ - fContainerWallZ);
106  return trans;
107 }
108 
109 
110 //------------------------------------------------------------------
112 ComputeSolid(const G4int, G4VPhysicalVolume* pPhysicalVol)
113 {
114  return pPhysicalVol->GetLogicalVolume()->GetSolid();
115 }
116 
117 
118 //------------------------------------------------------------------
121 {
122  CheckCopyNo( copyNo );
123  size_t matIndex = GetMaterialIndex(copyNo);
124 
125  return fMaterials[ matIndex ];
126 }
127 
128 
129 //------------------------------------------------------------------
131 GetMaterialIndex( size_t copyNo ) const
132 {
133  CheckCopyNo( copyNo );
134 
135  if( fMaterialIndices == nullptr ) { return 0; }
136  return *(fMaterialIndices+copyNo);
137 }
138 
139 
140 //------------------------------------------------------------------
142 GetMaterialIndex( size_t nx, size_t ny, size_t nz ) const
143 {
144  size_t copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
145  return GetMaterialIndex( copyNo );
146 }
147 
148 
149 //------------------------------------------------------------------
150 G4Material*
151 G4PhantomParameterisation::GetMaterial( size_t nx, size_t ny, size_t nz) const
152 {
153  return fMaterials[GetMaterialIndex(nx,ny,nz)];
154 }
155 
156 //------------------------------------------------------------------
158 {
159  return fMaterials[GetMaterialIndex(copyNo)];
160 }
161 
162 //------------------------------------------------------------------
164 ComputeVoxelIndices(const G4int copyNo, size_t& nx,
165  size_t& ny, size_t& nz ) const
166 {
167  CheckCopyNo( copyNo );
168  nx = size_t(copyNo%fNoVoxelX);
169  ny = size_t( (copyNo/fNoVoxelX)%fNoVoxelY );
170  nz = size_t(copyNo/fNoVoxelXY);
171 }
172 
173 
174 //------------------------------------------------------------------
177 {
178  G4double toleranceForWarning = 0.25*kCarTolerance;
179 
180  // Any bigger value than 0.25*kCarTolerance will give a warning in
181  // G4NormalNavigation::ComputeStep(), because the Inverse of a container
182  // translation that is Z+epsilon gives -Z+epsilon (and the maximum tolerance
183  // in G4Box::Inside is 0.5*kCarTolerance
184  //
185  G4double toleranceForError = 1.*kCarTolerance;
186 
187  // Any bigger value than kCarTolerance will give an error in GetReplicaNo()
188  //
189  if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForError
190  || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForError
191  || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForError )
192  {
193  std::ostringstream message;
194  message << "Voxels do not fully fill the container: "
196  << " DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
197  << " DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
198  << " DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
199  << " Maximum difference is: " << toleranceForError;
200  G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
201  "GeomNav0002", FatalException, message);
202 
203  }
204  else if( std::fabs(contX-fNoVoxelX*fVoxelHalfX) >= toleranceForWarning
205  || std::fabs(contY-fNoVoxelY*fVoxelHalfY) >= toleranceForWarning
206  || std::fabs(contZ-fNoVoxelZ*fVoxelHalfZ) >= toleranceForWarning )
207  {
208  std::ostringstream message;
209  message << "Voxels do not fully fill the container: "
211  << " DiffX= " << contX-fNoVoxelX*fVoxelHalfX << G4endl
212  << " DiffY= " << contY-fNoVoxelY*fVoxelHalfY << G4endl
213  << " DiffZ= " << contZ-fNoVoxelZ*fVoxelHalfZ << G4endl
214  << " Maximum difference is: " << toleranceForWarning;
215  G4Exception("G4PhantomParameterisation::CheckVoxelsFillContainer()",
216  "GeomNav1002", JustWarning, message);
217  }
218 }
219 
220 
221 //------------------------------------------------------------------
223 GetReplicaNo( const G4ThreeVector& localPoint, const G4ThreeVector& localDir )
224 {
225 
226  // Check first that point is really inside voxels
227  //
228  if( fContainerSolid->Inside( localPoint ) == kOutside )
229  {
230  std::ostringstream message;
231  message << "Point outside voxels!" << G4endl
232  << " localPoint - " << localPoint
233  << " - is outside container solid: "
235  << "DIFFERENCE WITH PHANTOM WALLS X: "
236  << std::fabs(localPoint.x()) - fContainerWallX
237  << " Y: " << std::fabs(localPoint.y()) - fContainerWallY
238  << " Z: " << std::fabs(localPoint.z()) - fContainerWallZ;
239  G4Exception("G4PhantomParameterisation::GetReplicaNo()", "GeomNav0003",
240  FatalErrorInArgument, message);
241  }
242 
243  // Check the voxel numbers corresponding to localPoint
244  // When a particle is on a surface, it may be between -kCarTolerance and
245  // +kCartolerance. By a simple distance as:
246  // G4int nx = G4int( (localPoint.x()+)/fVoxelHalfX/2.);
247  // those between -kCartolerance and 0 will be placed on voxel N-1 and those
248  // between 0 and kCarTolerance on voxel N.
249  // To avoid precision problems place the tracks that are on the surface on
250  // voxel N-1 if they have negative direction and on voxel N if they have
251  // positive direction.
252  // Add +kCarTolerance so that they are first placed on voxel N, and then
253  // if the direction is negative substract 1
254 
255  G4double fx = (localPoint.x()+fContainerWallX+kCarTolerance)/(fVoxelHalfX*2.);
256  G4int nx = G4int(fx);
257 
258  G4double fy = (localPoint.y()+fContainerWallY+kCarTolerance)/(fVoxelHalfY*2.);
259  G4int ny = G4int(fy);
260 
261  G4double fz = (localPoint.z()+fContainerWallZ+kCarTolerance)/(fVoxelHalfZ*2.);
262  G4int nz = G4int(fz);
263 
264  // If it is on the surface side, check the direction: if direction is
265  // negative place it in the previous voxel (if direction is positive it is
266  // already in the next voxel).
267  // Correct also cases where n = -1 or n = fNoVoxel. It is always traced to be
268  // due to multiple scattering: track is entering a voxel but multiple
269  // scattering changes the angle towards outside
270  //
271  if( fx - nx < kCarTolerance*fVoxelHalfX )
272  {
273  if( localDir.x() < 0 )
274  {
275  if( nx != 0 )
276  {
277  nx -= 1;
278  }
279  }
280  else
281  {
282  if( nx == G4int(fNoVoxelX) )
283  {
284  nx -= 1;
285  }
286  }
287  }
288  if( fy - ny < kCarTolerance*fVoxelHalfY )
289  {
290  if( localDir.y() < 0 )
291  {
292  if( ny != 0 )
293  {
294  ny -= 1;
295  }
296  }
297  else
298  {
299  if( ny == G4int(fNoVoxelY) )
300  {
301  ny -= 1;
302  }
303  }
304  }
305  if( fz - nz < kCarTolerance*fVoxelHalfZ )
306  {
307  if( localDir.z() < 0 )
308  {
309  if( nz != 0 )
310  {
311  nz -= 1;
312  }
313  }
314  else
315  {
316  if( nz == G4int(fNoVoxelZ) )
317  {
318  nz -= 1;
319  }
320  }
321  }
322 
323  G4int copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
324 
325  // Check if there are still errors
326  //
327  G4bool isOK = true;
328  if( nx < 0 )
329  {
330  nx = 0;
331  isOK = false;
332  }
333  else if( nx >= G4int(fNoVoxelX) )
334  {
335  nx = fNoVoxelX-1;
336  isOK = false;
337  }
338  if( ny < 0 )
339  {
340  ny = 0;
341  isOK = false;
342  }
343  else if( ny >= G4int(fNoVoxelY) )
344  {
345  ny = fNoVoxelY-1;
346  isOK = false;
347  }
348  if( nz < 0 )
349  {
350  nz = 0;
351  isOK = false;
352  }
353  else if( nz >= G4int(fNoVoxelZ) )
354  {
355  nz = fNoVoxelZ-1;
356  isOK = false;
357  }
358  if( !isOK )
359  {
360  std::ostringstream message;
361  message << "Corrected the copy number! It was negative or too big" << G4endl
362  << " LocalPoint: " << localPoint << G4endl
363  << " LocalDir: " << localDir << G4endl
364  << " Voxel container size: " << fContainerWallX
365  << " " << fContainerWallY << " " << fContainerWallZ << G4endl
366  << " LocalPoint - wall: "
367  << localPoint.x()-fContainerWallX << " "
368  << localPoint.y()-fContainerWallY << " "
369  << localPoint.z()-fContainerWallZ;
370  G4Exception("G4PhantomParameterisation::GetReplicaNo()",
371  "GeomNav1002", JustWarning, message);
372  copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
373  }
374 
375  return copyNo;
376 }
377 
378 
379 //------------------------------------------------------------------
381 {
382  if( copyNo < 0 || copyNo >= G4int(fNoVoxel) )
383  {
384  std::ostringstream message;
385  message << "Copy number is negative or too big!" << G4endl
386  << " Copy number: " << copyNo << G4endl
387  << " Total number of voxels: " << fNoVoxel;
388  G4Exception("G4PhantomParameterisation::CheckCopyNo()",
389  "GeomNav0002", FatalErrorInArgument, message);
390  }
391 }