Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // class G4PartialPhantomParameterisation impl 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, G4V 46 { 47 // Voxels cannot be rotated, return translat 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 - 64 (2*ny+1)*fVoxelHalfY - 65 (2*nz+1)*fVoxelHalfZ - 66 return trans; 67 } 68 69 70 //-------------------------------------------- 71 G4Material* G4PartialPhantomParameterisation:: 72 ComputeMaterial( const G4int copyNo, G4VPhysic 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 96 { 97 std::size_t copyNo = nx + fNoVoxelsX*ny + fN 98 return GetMaterialIndex( copyNo ); 99 } 100 101 102 //-------------------------------------------- 103 G4Material* G4PartialPhantomParameterisation:: 104 GetMaterial( std::size_t nx, std::size_t ny, s 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::s 121 std::size_t& ny, std 122 { 123 CheckCopyNo( copyNo ); 124 125 auto ite = fFilledIDs.lower_bound(copyNo); 126 G4long dist = std::distance( fFilledIDs.cbeg 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, 149 { 150 // Check the voxel numbers corresponding to 151 // When a particle is on a surface, it may b 152 // +kCartolerance. By a simple distance as: 153 // G4int nx = G4int( (localPoint.x()+)/fVo 154 // those between -kCartolerance and 0 will b 155 // between 0 and kCarTolerance on voxel N. 156 // To avoid precision problems place the tra 157 // voxel N-1 if they have negative direction 158 // positive direction. 159 // Add +kCarTolerance so that they are first 160 // if the direction is negative substract 1 161 162 G4double fx = (localPoint.x()+fContainerWall 163 auto nx = G4int(fx); 164 165 G4double fy = (localPoint.y()+fContainerWall 166 auto ny = G4int(fy); 167 168 G4double fz = (localPoint.z()+fContainerWall 169 auto nz = G4int(fz); 170 171 // If it is on the surface side, check the d 172 // negative place it on the previous voxel ( 173 // already in the next voxel...). 174 // Correct also cases where n = -1 or n = fN 175 // due to multiple scattering: track is ente 176 // scattering changes the angle towards outs 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 267 << G4endl 268 << " LocalPoint: " << loc 269 << " LocalDir: " << local 270 << " Voxel container size 271 << " " << fContainerWallY << " " < 272 << " LocalPoint - wall: " 273 << localPoint.x()-fContainerWallX 274 << localPoint.y()-fContainerWallY 275 << localPoint.z()-fContainerWallZ; 276 G4Exception("G4PartialPhantomParameterisat 277 "GeomNav1002", JustWarning, me 278 } 279 280 auto nyz = G4int(nz*fNoVoxelsY+ny); 281 auto ite = fFilledIDs.cbegin(); 282 /* 283 for( ite = fFilledIDs.cbegin(); ite != fFill 284 { 285 G4cout << " G4PartialPhantomParameterisati 286 << (*ite).first << " , " << (*ite). 287 } 288 */ 289 290 advance(ite,nyz); 291 auto iteant = ite; iteant--; 292 G4int copyNo = (*iteant).first + 1 + ( nx - 293 /* 294 G4cout << " G4PartialPhantomParameterisation 295 << copyNo << " nyz " << nyz << " (* 296 << (*iteant).first << " (*ite).secon 297 298 G4cout << " G4PartialPhantomParameterisation 299 << " nx " << nx << " ny " << ny << " 300 << " localPoint " << localPoint << " 301 */ 302 return copyNo; 303 } 304 305 306 //-------------------------------------------- 307 void G4PartialPhantomParameterisation::CheckCo 308 { 309 if( copyNo < 0 || copyNo >= G4int(fNoVoxels) 310 { 311 std::ostringstream message; 312 message << "Copy number is negative or too 313 << " Copy number: " << copy 314 << " Total number of voxels 315 G4Exception("G4PartialPhantomParameterisat 316 "GeomNav0002", FatalErrorInArg 317 } 318 } 319 320 321 //-------------------------------------------- 322 void G4PartialPhantomParameterisation::BuildCo 323 { 324 fContainerWallX = fNoVoxelsX * fVoxelHalfX; 325 fContainerWallY = fNoVoxelsY * fVoxelHalfY; 326 fContainerWallZ = fNoVoxelsZ * fVoxelHalfZ; 327 } 328