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