Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // class G4PartialPhantomParameterisation impl 26 // class G4PartialPhantomParameterisation implementation 27 // 27 // 28 // May 2007 Pedro Arce (CIEMAT), first version 28 // May 2007 Pedro Arce (CIEMAT), first version 29 // ------------------------------------------- 29 // -------------------------------------------------------------------- 30 30 31 #include "G4PartialPhantomParameterisation.hh" 31 #include "G4PartialPhantomParameterisation.hh" 32 32 33 #include "globals.hh" 33 #include "globals.hh" 34 #include "G4Material.hh" 34 #include "G4Material.hh" 35 #include "G4VSolid.hh" 35 #include "G4VSolid.hh" 36 #include "G4VPhysicalVolume.hh" 36 #include "G4VPhysicalVolume.hh" 37 #include "G4LogicalVolume.hh" 37 #include "G4LogicalVolume.hh" 38 #include "G4VVolumeMaterialScanner.hh" 38 #include "G4VVolumeMaterialScanner.hh" 39 #include "G4GeometryTolerance.hh" 39 #include "G4GeometryTolerance.hh" 40 40 41 #include <list> 41 #include <list> 42 42 43 //-------------------------------------------- 43 //------------------------------------------------------------------ >> 44 G4PartialPhantomParameterisation::G4PartialPhantomParameterisation() >> 45 { >> 46 } >> 47 >> 48 >> 49 //------------------------------------------------------------------ >> 50 G4PartialPhantomParameterisation::~G4PartialPhantomParameterisation() = default; >> 51 >> 52 //------------------------------------------------------------------ 44 void G4PartialPhantomParameterisation:: 53 void G4PartialPhantomParameterisation:: 45 ComputeTransformation( const G4int copyNo, G4V 54 ComputeTransformation( const G4int copyNo, G4VPhysicalVolume *physVol ) const 46 { 55 { 47 // Voxels cannot be rotated, return translat 56 // Voxels cannot be rotated, return translation 48 // 57 // 49 G4ThreeVector trans = GetTranslation( copyNo 58 G4ThreeVector trans = GetTranslation( copyNo ); 50 physVol->SetTranslation( trans ); 59 physVol->SetTranslation( trans ); 51 } 60 } 52 61 53 62 54 //-------------------------------------------- 63 //------------------------------------------------------------------ 55 G4ThreeVector G4PartialPhantomParameterisation 64 G4ThreeVector G4PartialPhantomParameterisation:: 56 GetTranslation(const G4int copyNo ) const 65 GetTranslation(const G4int copyNo ) const 57 { 66 { 58 CheckCopyNo( copyNo ); 67 CheckCopyNo( copyNo ); 59 68 60 std::size_t nx, ny, nz; 69 std::size_t nx, ny, nz; 61 ComputeVoxelIndices( copyNo, nx, ny, nz ); 70 ComputeVoxelIndices( copyNo, nx, ny, nz ); 62 71 63 G4ThreeVector trans( (2*nx+1)*fVoxelHalfX - 72 G4ThreeVector trans( (2*nx+1)*fVoxelHalfX - fContainerWallX, 64 (2*ny+1)*fVoxelHalfY - 73 (2*ny+1)*fVoxelHalfY - fContainerWallY, 65 (2*nz+1)*fVoxelHalfZ - 74 (2*nz+1)*fVoxelHalfZ - fContainerWallZ); 66 return trans; 75 return trans; 67 } 76 } 68 77 69 78 70 //-------------------------------------------- 79 //------------------------------------------------------------------ 71 G4Material* G4PartialPhantomParameterisation:: 80 G4Material* G4PartialPhantomParameterisation:: 72 ComputeMaterial( const G4int copyNo, G4VPhysic 81 ComputeMaterial( const G4int copyNo, G4VPhysicalVolume*, const G4VTouchable* ) 73 { 82 { 74 CheckCopyNo( copyNo ); 83 CheckCopyNo( copyNo ); 75 auto matIndex = GetMaterialIndex(copyNo); 84 auto matIndex = GetMaterialIndex(copyNo); 76 85 77 return fMaterials[ matIndex ]; 86 return fMaterials[ matIndex ]; 78 } 87 } 79 88 80 89 81 //-------------------------------------------- 90 //------------------------------------------------------------------ 82 size_t G4PartialPhantomParameterisation:: 91 size_t G4PartialPhantomParameterisation:: 83 GetMaterialIndex( std::size_t copyNo ) const 92 GetMaterialIndex( std::size_t copyNo ) const 84 { 93 { 85 CheckCopyNo( copyNo ); 94 CheckCopyNo( copyNo ); 86 95 87 if( fMaterialIndices == nullptr ) { return 0 96 if( fMaterialIndices == nullptr ) { return 0; } 88 97 89 return *(fMaterialIndices+copyNo); 98 return *(fMaterialIndices+copyNo); 90 } 99 } 91 100 92 101 93 //-------------------------------------------- 102 //------------------------------------------------------------------ 94 size_t G4PartialPhantomParameterisation:: 103 size_t G4PartialPhantomParameterisation:: 95 GetMaterialIndex( std::size_t nx, std::size_t 104 GetMaterialIndex( std::size_t nx, std::size_t ny, std::size_t nz ) const 96 { 105 { 97 std::size_t copyNo = nx + fNoVoxelsX*ny + fN 106 std::size_t copyNo = nx + fNoVoxelsX*ny + fNoVoxelsXY*nz; 98 return GetMaterialIndex( copyNo ); 107 return GetMaterialIndex( copyNo ); 99 } 108 } 100 109 101 110 102 //-------------------------------------------- 111 //------------------------------------------------------------------ 103 G4Material* G4PartialPhantomParameterisation:: 112 G4Material* G4PartialPhantomParameterisation:: 104 GetMaterial( std::size_t nx, std::size_t ny, s 113 GetMaterial( std::size_t nx, std::size_t ny, std::size_t nz) const 105 { 114 { 106 return fMaterials[GetMaterialIndex(nx,ny,nz) 115 return fMaterials[GetMaterialIndex(nx,ny,nz)]; 107 } 116 } 108 117 109 118 110 //-------------------------------------------- 119 //------------------------------------------------------------------ 111 G4Material* G4PartialPhantomParameterisation:: 120 G4Material* G4PartialPhantomParameterisation:: 112 GetMaterial( std::size_t copyNo ) const 121 GetMaterial( std::size_t copyNo ) const 113 { 122 { 114 return fMaterials[GetMaterialIndex(copyNo)]; 123 return fMaterials[GetMaterialIndex(copyNo)]; 115 } 124 } 116 125 117 126 118 //-------------------------------------------- 127 //------------------------------------------------------------------ 119 void G4PartialPhantomParameterisation:: 128 void G4PartialPhantomParameterisation:: 120 ComputeVoxelIndices(const G4int copyNo, std::s 129 ComputeVoxelIndices(const G4int copyNo, std::size_t& nx, 121 std::size_t& ny, std 130 std::size_t& ny, std::size_t& nz ) const 122 { 131 { 123 CheckCopyNo( copyNo ); 132 CheckCopyNo( copyNo ); 124 133 125 auto ite = fFilledIDs.lower_bound(copyNo); 134 auto ite = fFilledIDs.lower_bound(copyNo); 126 G4long dist = std::distance( fFilledIDs.cbeg 135 G4long dist = std::distance( fFilledIDs.cbegin(), ite ); 127 nz = std::size_t( dist/fNoVoxelsY ); 136 nz = std::size_t( dist/fNoVoxelsY ); 128 ny = std::size_t( dist%fNoVoxelsY ); 137 ny = std::size_t( dist%fNoVoxelsY ); 129 138 130 G4int ifmin = (*ite).second; 139 G4int ifmin = (*ite).second; 131 G4int nvoxXprev; 140 G4int nvoxXprev; 132 if( dist != 0 ) 141 if( dist != 0 ) 133 { 142 { 134 ite--; 143 ite--; 135 nvoxXprev = (*ite).first; 144 nvoxXprev = (*ite).first; 136 } 145 } 137 else 146 else 138 { 147 { 139 nvoxXprev = -1; 148 nvoxXprev = -1; 140 } 149 } 141 150 142 nx = ifmin+copyNo-nvoxXprev-1; 151 nx = ifmin+copyNo-nvoxXprev-1; 143 } 152 } 144 153 145 154 146 //-------------------------------------------- 155 //------------------------------------------------------------------ 147 G4int G4PartialPhantomParameterisation:: 156 G4int G4PartialPhantomParameterisation:: 148 GetReplicaNo( const G4ThreeVector& localPoint, 157 GetReplicaNo( const G4ThreeVector& localPoint, const G4ThreeVector& localDir ) 149 { 158 { 150 // Check the voxel numbers corresponding to 159 // Check the voxel numbers corresponding to localPoint 151 // When a particle is on a surface, it may b 160 // When a particle is on a surface, it may be between -kCarTolerance and 152 // +kCartolerance. By a simple distance as: 161 // +kCartolerance. By a simple distance as: 153 // G4int nx = G4int( (localPoint.x()+)/fVo 162 // G4int nx = G4int( (localPoint.x()+)/fVoxelHalfX/2.); 154 // those between -kCartolerance and 0 will b 163 // those between -kCartolerance and 0 will be placed on voxel N-1 and those 155 // between 0 and kCarTolerance on voxel N. 164 // between 0 and kCarTolerance on voxel N. 156 // To avoid precision problems place the tra 165 // To avoid precision problems place the tracks that are on the surface on 157 // voxel N-1 if they have negative direction 166 // voxel N-1 if they have negative direction and on voxel N if they have 158 // positive direction. 167 // positive direction. 159 // Add +kCarTolerance so that they are first 168 // Add +kCarTolerance so that they are first placed on voxel N, and then 160 // if the direction is negative substract 1 169 // if the direction is negative substract 1 161 170 162 G4double fx = (localPoint.x()+fContainerWall 171 G4double fx = (localPoint.x()+fContainerWallX+kCarTolerance)/(fVoxelHalfX*2.); 163 auto nx = G4int(fx); 172 auto nx = G4int(fx); 164 173 165 G4double fy = (localPoint.y()+fContainerWall 174 G4double fy = (localPoint.y()+fContainerWallY+kCarTolerance)/(fVoxelHalfY*2.); 166 auto ny = G4int(fy); 175 auto ny = G4int(fy); 167 176 168 G4double fz = (localPoint.z()+fContainerWall 177 G4double fz = (localPoint.z()+fContainerWallZ+kCarTolerance)/(fVoxelHalfZ*2.); 169 auto nz = G4int(fz); 178 auto nz = G4int(fz); 170 179 171 // If it is on the surface side, check the d 180 // If it is on the surface side, check the direction: if direction is 172 // negative place it on the previous voxel ( 181 // negative place it on the previous voxel (if direction is positive it is 173 // already in the next voxel...). 182 // already in the next voxel...). 174 // Correct also cases where n = -1 or n = fN 183 // Correct also cases where n = -1 or n = fNoVoxels. It is always traced to be 175 // due to multiple scattering: track is ente 184 // due to multiple scattering: track is entering a voxel but multiple 176 // scattering changes the angle towards outs 185 // scattering changes the angle towards outside 177 // 186 // 178 if( fx - nx < kCarTolerance/fVoxelHalfX ) 187 if( fx - nx < kCarTolerance/fVoxelHalfX ) 179 { 188 { 180 if( localDir.x() < 0 ) 189 if( localDir.x() < 0 ) 181 { 190 { 182 if( nx != 0 ) 191 if( nx != 0 ) 183 { 192 { 184 nx -= 1; 193 nx -= 1; 185 } 194 } 186 } 195 } 187 else 196 else 188 { 197 { 189 if( nx == G4int(fNoVoxelsX) ) 198 if( nx == G4int(fNoVoxelsX) ) 190 { 199 { 191 nx -= 1; 200 nx -= 1; 192 } 201 } 193 } 202 } 194 } 203 } 195 if( fy - ny < kCarTolerance/fVoxelHalfY ) 204 if( fy - ny < kCarTolerance/fVoxelHalfY ) 196 { 205 { 197 if( localDir.y() < 0 ) 206 if( localDir.y() < 0 ) 198 { 207 { 199 if( ny != 0 ) 208 if( ny != 0 ) 200 { 209 { 201 ny -= 1; 210 ny -= 1; 202 } 211 } 203 } 212 } 204 else 213 else 205 { 214 { 206 if( ny == G4int(fNoVoxelsY) ) 215 if( ny == G4int(fNoVoxelsY) ) 207 { 216 { 208 ny -= 1; 217 ny -= 1; 209 } 218 } 210 } 219 } 211 } 220 } 212 if( fz - nz < kCarTolerance/fVoxelHalfZ ) 221 if( fz - nz < kCarTolerance/fVoxelHalfZ ) 213 { 222 { 214 if( localDir.z() < 0 ) 223 if( localDir.z() < 0 ) 215 { 224 { 216 if( nz != 0 ) 225 if( nz != 0 ) 217 { 226 { 218 nz -= 1; 227 nz -= 1; 219 } 228 } 220 } 229 } 221 else 230 else 222 { 231 { 223 if( nz == G4int(fNoVoxelsZ) ) 232 if( nz == G4int(fNoVoxelsZ) ) 224 { 233 { 225 nz -= 1; 234 nz -= 1; 226 } 235 } 227 } 236 } 228 } 237 } 229 238 230 // Check if there are still errors 239 // Check if there are still errors 231 // 240 // 232 G4bool isOK = true; 241 G4bool isOK = true; 233 if( nx < 0 ) 242 if( nx < 0 ) 234 { 243 { 235 nx = 0; 244 nx = 0; 236 isOK = false; 245 isOK = false; 237 } 246 } 238 else if( nx >= G4int(fNoVoxelsX) ) 247 else if( nx >= G4int(fNoVoxelsX) ) 239 { 248 { 240 nx = G4int(fNoVoxelsX)-1; 249 nx = G4int(fNoVoxelsX)-1; 241 isOK = false; 250 isOK = false; 242 } 251 } 243 if( ny < 0 ) 252 if( ny < 0 ) 244 { 253 { 245 ny = 0; 254 ny = 0; 246 isOK = false; 255 isOK = false; 247 } 256 } 248 else if( ny >= G4int(fNoVoxelsY) ) 257 else if( ny >= G4int(fNoVoxelsY) ) 249 { 258 { 250 ny = G4int(fNoVoxelsY)-1; 259 ny = G4int(fNoVoxelsY)-1; 251 isOK = false; 260 isOK = false; 252 } 261 } 253 if( nz < 0 ) 262 if( nz < 0 ) 254 { 263 { 255 nz = 0; 264 nz = 0; 256 isOK = false; 265 isOK = false; 257 } 266 } 258 else if( nz >= G4int(fNoVoxelsZ) ) 267 else if( nz >= G4int(fNoVoxelsZ) ) 259 { 268 { 260 nz = G4int(fNoVoxelsZ)-1; 269 nz = G4int(fNoVoxelsZ)-1; 261 isOK = false; 270 isOK = false; 262 } 271 } 263 if( !isOK ) 272 if( !isOK ) 264 { 273 { 265 std::ostringstream message; 274 std::ostringstream message; 266 message << "Corrected the copy number! It 275 message << "Corrected the copy number! It was negative or too big." 267 << G4endl 276 << G4endl 268 << " LocalPoint: " << loc 277 << " LocalPoint: " << localPoint << G4endl 269 << " LocalDir: " << local 278 << " LocalDir: " << localDir << G4endl 270 << " Voxel container size 279 << " Voxel container size: " << fContainerWallX 271 << " " << fContainerWallY << " " < 280 << " " << fContainerWallY << " " << fContainerWallZ << G4endl 272 << " LocalPoint - wall: " 281 << " LocalPoint - wall: " 273 << localPoint.x()-fContainerWallX 282 << localPoint.x()-fContainerWallX << " " 274 << localPoint.y()-fContainerWallY 283 << localPoint.y()-fContainerWallY << " " 275 << localPoint.z()-fContainerWallZ; 284 << localPoint.z()-fContainerWallZ; 276 G4Exception("G4PartialPhantomParameterisat 285 G4Exception("G4PartialPhantomParameterisation::GetReplicaNo()", 277 "GeomNav1002", JustWarning, me 286 "GeomNav1002", JustWarning, message); 278 } 287 } 279 288 280 auto nyz = G4int(nz*fNoVoxelsY+ny); 289 auto nyz = G4int(nz*fNoVoxelsY+ny); 281 auto ite = fFilledIDs.cbegin(); 290 auto ite = fFilledIDs.cbegin(); 282 /* 291 /* 283 for( ite = fFilledIDs.cbegin(); ite != fFill 292 for( ite = fFilledIDs.cbegin(); ite != fFilledIDs.cend(); ++ite ) 284 { 293 { 285 G4cout << " G4PartialPhantomParameterisati 294 G4cout << " G4PartialPhantomParameterisation::GetReplicaNo filled " 286 << (*ite).first << " , " << (*ite). 295 << (*ite).first << " , " << (*ite).second << std::endl; 287 } 296 } 288 */ 297 */ 289 298 290 advance(ite,nyz); 299 advance(ite,nyz); 291 auto iteant = ite; iteant--; 300 auto iteant = ite; iteant--; 292 G4int copyNo = (*iteant).first + 1 + ( nx - 301 G4int copyNo = (*iteant).first + 1 + ( nx - (*ite).second ); 293 /* 302 /* 294 G4cout << " G4PartialPhantomParameterisation 303 G4cout << " G4PartialPhantomParameterisation::GetReplicaNo getting copyNo " 295 << copyNo << " nyz " << nyz << " (* 304 << copyNo << " nyz " << nyz << " (*iteant).first " 296 << (*iteant).first << " (*ite).secon 305 << (*iteant).first << " (*ite).second " << (*ite).second << G4endl; 297 306 298 G4cout << " G4PartialPhantomParameterisation 307 G4cout << " G4PartialPhantomParameterisation::GetReplicaNo " << copyNo 299 << " nx " << nx << " ny " << ny << " 308 << " nx " << nx << " ny " << ny << " nz " << nz 300 << " localPoint " << localPoint << " 309 << " localPoint " << localPoint << " localDir " << localDir << G4endl; 301 */ 310 */ 302 return copyNo; 311 return copyNo; 303 } 312 } 304 313 305 314 306 //-------------------------------------------- 315 //------------------------------------------------------------------ 307 void G4PartialPhantomParameterisation::CheckCo 316 void G4PartialPhantomParameterisation::CheckCopyNo( const G4long copyNo ) const 308 { 317 { 309 if( copyNo < 0 || copyNo >= G4int(fNoVoxels) 318 if( copyNo < 0 || copyNo >= G4int(fNoVoxels) ) 310 { 319 { 311 std::ostringstream message; 320 std::ostringstream message; 312 message << "Copy number is negative or too 321 message << "Copy number is negative or too big!" << G4endl 313 << " Copy number: " << copy 322 << " Copy number: " << copyNo << G4endl 314 << " Total number of voxels 323 << " Total number of voxels: " << fNoVoxels; 315 G4Exception("G4PartialPhantomParameterisat 324 G4Exception("G4PartialPhantomParameterisation::CheckCopyNo()", 316 "GeomNav0002", FatalErrorInArg 325 "GeomNav0002", FatalErrorInArgument, message); 317 } 326 } 318 } 327 } 319 328 320 329 321 //-------------------------------------------- 330 //------------------------------------------------------------------ 322 void G4PartialPhantomParameterisation::BuildCo 331 void G4PartialPhantomParameterisation::BuildContainerWalls() 323 { 332 { 324 fContainerWallX = fNoVoxelsX * fVoxelHalfX; 333 fContainerWallX = fNoVoxelsX * fVoxelHalfX; 325 fContainerWallY = fNoVoxelsY * fVoxelHalfY; 334 fContainerWallY = fNoVoxelsY * fVoxelHalfY; 326 fContainerWallZ = fNoVoxelsZ * fVoxelHalfZ; 335 fContainerWallZ = fNoVoxelsZ * fVoxelHalfZ; 327 } 336 } 328 337