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