Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PartialPhantomParameterisation.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$
28 //
29 //
30 // class G4PartialPhantomParameterisation implementation
31 //
32 // May 2007 Pedro Arce (CIEMAT), first version
33 //
34 // --------------------------------------------------------------------
35 
37 
38 #include "globals.hh"
39 #include "G4Material.hh"
40 #include "G4VSolid.hh"
41 #include "G4VPhysicalVolume.hh"
42 #include "G4LogicalVolume.hh"
44 #include "G4GeometryTolerance.hh"
45 
46 #include <list>
47 
48 //------------------------------------------------------------------
51 {
52 }
53 
54 
55 //------------------------------------------------------------------
57 {
58 }
59 
60 //------------------------------------------------------------------
62 ComputeTransformation(const G4int copyNo, G4VPhysicalVolume *physVol ) const
63 {
64  // Voxels cannot be rotated, return translation
65  //
66  G4ThreeVector trans = GetTranslation( copyNo );
67  physVol->SetTranslation( trans );
68 }
69 
70 
71 //------------------------------------------------------------------
73 GetTranslation(const G4int copyNo ) const
74 {
75  CheckCopyNo( copyNo );
76 
77  size_t nx;
78  size_t ny;
79  size_t nz;
80  ComputeVoxelIndices( copyNo, nx, ny, nz );
81 
82  G4ThreeVector trans( (2*nx+1)*fVoxelHalfX - fContainerWallX,
83  (2*ny+1)*fVoxelHalfY - fContainerWallY,
84  (2*nz+1)*fVoxelHalfZ - fContainerWallZ);
85  return trans;
86 }
87 
88 
89 //------------------------------------------------------------------
92 {
93  CheckCopyNo( copyNo );
94  size_t matIndex = GetMaterialIndex(copyNo);
95 
96  return fMaterials[ matIndex ];
97 }
98 
99 
100 //------------------------------------------------------------------
102 GetMaterialIndex( size_t copyNo ) const
103 {
104  CheckCopyNo( copyNo );
105 
106  if( !fMaterialIndices ) { return 0; }
107 
108  return *(fMaterialIndices+copyNo);
109 }
110 
111 
112 //------------------------------------------------------------------
114 GetMaterialIndex( size_t nx, size_t ny, size_t nz ) const
115 {
116  size_t copyNo = nx + fNoVoxelX*ny + fNoVoxelXY*nz;
117  return GetMaterialIndex( copyNo );
118 }
119 
120 
121 //------------------------------------------------------------------
123 GetMaterial( size_t nx, size_t ny, size_t nz) const
124 {
125  return fMaterials[GetMaterialIndex(nx,ny,nz)];
126 }
127 
128 
129 //------------------------------------------------------------------
131 GetMaterial( size_t copyNo ) const
132 {
133  return fMaterials[GetMaterialIndex(copyNo)];
134 }
135 
136 
137 //------------------------------------------------------------------
138 void G4PartialPhantomParameterisation::
139 ComputeVoxelIndices(const G4int copyNo, size_t& nx,
140  size_t& ny, size_t& nz ) const
141 {
142  CheckCopyNo( copyNo );
143 
144  std::multimap<G4int,G4int>::const_iterator ite =
145  fFilledIDs.lower_bound(size_t(copyNo));
146  G4int dist = std::distance( fFilledIDs.begin(), ite );
147  nz = size_t(dist/fNoVoxelY);
148  ny = size_t( dist%fNoVoxelY );
149 
150  G4int ifmin = (*ite).second;
151  G4int nvoxXprev;
152  if( dist != 0 ) {
153  ite--;
154  nvoxXprev = (*ite).first;
155  } else {
156  nvoxXprev = -1;
157  }
158 
159  nx = ifmin+copyNo-nvoxXprev-1;
160 }
161 
162 
163 //------------------------------------------------------------------
165 GetReplicaNo( const G4ThreeVector& localPoint, const G4ThreeVector& localDir )
166 {
167  // Check the voxel numbers corresponding to localPoint
168  // When a particle is on a surface, it may be between -kCarTolerance and
169  // +kCartolerance. By a simple distance as:
170  // G4int nx = G4int( (localPoint.x()+)/fVoxelHalfX/2.);
171  // those between -kCartolerance and 0 will be placed on voxel N-1 and those
172  // between 0 and kCarTolerance on voxel N.
173  // To avoid precision problems place the tracks that are on the surface on
174  // voxel N-1 if they have negative direction and on voxel N if they have
175  // positive direction.
176  // Add +kCarTolerance so that they are first placed on voxel N, and then
177  // if the direction is negative substract 1
178 
179  G4double fx = (localPoint.x()+fContainerWallX+kCarTolerance)/(fVoxelHalfX*2.);
180  G4int nx = G4int(fx);
181 
182  G4double fy = (localPoint.y()+fContainerWallY+kCarTolerance)/(fVoxelHalfY*2.);
183  G4int ny = G4int(fy);
184 
185  G4double fz = (localPoint.z()+fContainerWallZ+kCarTolerance)/(fVoxelHalfZ*2.);
186  G4int nz = G4int(fz);
187 
188  // If it is on the surface side, check the direction: if direction is
189  // negative place it on the previous voxel (if direction is positive it is
190  // already in the next voxel...).
191  // Correct also cases where n = -1 or n = fNoVoxel. It is always traced to be
192  // due to multiple scattering: track is entering a voxel but multiple
193  // scattering changes the angle towards outside
194  //
195  if( fx - nx < kCarTolerance/fVoxelHalfX )
196  {
197  if( localDir.x() < 0 )
198  {
199  if( nx != 0 )
200  {
201  nx -= 1;
202  }
203  }
204  else
205  {
206  if( nx == G4int(fNoVoxelX) )
207  {
208  nx -= 1;
209  }
210  }
211  }
212  if( fy - ny < kCarTolerance/fVoxelHalfY )
213  {
214  if( localDir.y() < 0 )
215  {
216  if( ny != 0 )
217  {
218  ny -= 1;
219  }
220  }
221  else
222  {
223  if( ny == G4int(fNoVoxelY) )
224  {
225  ny -= 1;
226  }
227  }
228  }
229  if( fz - nz < kCarTolerance/fVoxelHalfZ )
230  {
231  if( localDir.z() < 0 )
232  {
233  if( nz != 0 )
234  {
235  nz -= 1;
236  }
237  }
238  else
239  {
240  if( nz == G4int(fNoVoxelZ) )
241  {
242  nz -= 1;
243  }
244  }
245  }
246 
247  // Check if there are still errors
248  //
249  G4bool isOK = true;
250  if( nx < 0 )
251  {
252  nx = 0;
253  isOK = false;
254  }
255  else if( nx >= G4int(fNoVoxelX) )
256  {
257  nx = fNoVoxelX-1;
258  isOK = false;
259  }
260  if( ny < 0 )
261  {
262  ny = 0;
263  isOK = false;
264  }
265  else if( ny >= G4int(fNoVoxelY) )
266  {
267  ny = fNoVoxelY-1;
268  isOK = false;
269  }
270  if( nz < 0 )
271  {
272  nz = 0;
273  isOK = false;
274  }
275  else if( nz >= G4int(fNoVoxelZ) )
276  {
277  nz = fNoVoxelZ-1;
278  isOK = false;
279  }
280  if( !isOK )
281  {
282  std::ostringstream message;
283  message << "Corrected the copy number! It was negative or too big."
284  << G4endl
285  << " LocalPoint: " << localPoint << G4endl
286  << " LocalDir: " << localDir << G4endl
287  << " Voxel container size: " << fContainerWallX
288  << " " << fContainerWallY << " " << fContainerWallZ << G4endl
289  << " LocalPoint - wall: "
290  << localPoint.x()-fContainerWallX << " "
291  << localPoint.y()-fContainerWallY << " "
292  << localPoint.z()-fContainerWallZ;
293  G4Exception("G4PartialPhantomParameterisation::GetReplicaNo()",
294  "GeomNav1002", JustWarning, message);
295  }
296 
297  G4int nyz = nz*fNoVoxelY+ny;
298  std::multimap<G4int,G4int>::iterator ite = fFilledIDs.begin();
299 /*
300  for( ite = fFilledIDs.begin(); ite != fFilledIDs.end(); ite++ )
301  {
302  G4cout << " G4PartialPhantomParameterisation::GetReplicaNo filled "
303  << (*ite).first << " , " << (*ite).second << std::endl;
304  }
305 */
306  ite = fFilledIDs.begin();
307 
308  advance(ite,nyz);
309  std::multimap<G4int,G4int>::iterator iteant = ite; iteant--;
310  G4int copyNo = (*iteant).first + 1 + ( nx - (*ite).second );
311 /*
312  G4cout << " G4PartialPhantomParameterisation::GetReplicaNo getting copyNo "
313  << copyNo << " nyz " << nyz << " (*iteant).first "
314  << (*iteant).first << " (*ite).second " << (*ite).second << G4endl;
315 
316  G4cout << " G4PartialPhantomParameterisation::GetReplicaNo " << copyNo
317  << " nx " << nx << " ny " << ny << " nz " << nz
318  << " localPoint " << localPoint << " localDir " << localDir << G4endl;
319 */
320  return copyNo;
321 }
322 
323 
324 //------------------------------------------------------------------
325 void G4PartialPhantomParameterisation::CheckCopyNo( const G4int copyNo ) const
326 {
327  if( copyNo < 0 || copyNo >= G4int(fNoVoxel) )
328  {
329  std::ostringstream message;
330  message << "Copy number is negative or too big!" << G4endl
331  << " Copy number: " << copyNo << G4endl
332  << " Total number of voxels: " << fNoVoxel;
333  G4Exception("G4PartialPhantomParameterisation::CheckCopyNo()",
334  "GeomNav0002", FatalErrorInArgument, message);
335  }
336 }
337 
338 
339 //------------------------------------------------------------------
341 {
345 }