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 // G4VSolid implementation for solid base clas 26 // G4VSolid implementation for solid base class 27 // 27 // 28 // 10.10.18 E.Tcherniaev, more robust Estimate 28 // 10.10.18 E.Tcherniaev, more robust EstimateSurfaceArea() based on distance 29 // 30.06.95 P.Kent, Created. 29 // 30.06.95 P.Kent, Created. 30 // ------------------------------------------- 30 // -------------------------------------------------------------------- 31 31 32 #include "G4VSolid.hh" 32 #include "G4VSolid.hh" 33 #include "G4SolidStore.hh" 33 #include "G4SolidStore.hh" 34 #include "globals.hh" 34 #include "globals.hh" 35 #include "G4QuickRand.hh" 35 #include "G4QuickRand.hh" 36 #include "G4GeometryTolerance.hh" 36 #include "G4GeometryTolerance.hh" 37 37 38 #include "G4VoxelLimits.hh" 38 #include "G4VoxelLimits.hh" 39 #include "G4AffineTransform.hh" 39 #include "G4AffineTransform.hh" 40 #include "G4VisExtent.hh" 40 #include "G4VisExtent.hh" 41 41 42 ////////////////////////////////////////////// 42 ////////////////////////////////////////////////////////////////////////// 43 // 43 // 44 // Streaming operator dumping solid contents << 45 << 46 std::ostream& operator<< ( std::ostream& os, c << 47 { << 48 return e.StreamInfo(os); << 49 } << 50 << 51 ////////////////////////////////////////////// << 52 // << 53 // Constructor 44 // Constructor 54 // - Copies name 45 // - Copies name 55 // - Add ourselves to solid Store 46 // - Add ourselves to solid Store 56 47 57 G4VSolid::G4VSolid(const G4String& name) 48 G4VSolid::G4VSolid(const G4String& name) 58 : fshapeName(name) 49 : fshapeName(name) 59 { 50 { 60 kCarTolerance = G4GeometryTolerance::GetIn 51 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 61 52 62 // Register to store 53 // Register to store 63 // 54 // 64 G4SolidStore::GetInstance()->Register(this 55 G4SolidStore::GetInstance()->Register(this); 65 } 56 } 66 57 67 ////////////////////////////////////////////// 58 ////////////////////////////////////////////////////////////////////////// 68 // 59 // 69 // Copy constructor 60 // Copy constructor 70 // 61 // 71 62 72 G4VSolid::G4VSolid(const G4VSolid& rhs) 63 G4VSolid::G4VSolid(const G4VSolid& rhs) 73 : kCarTolerance(rhs.kCarTolerance), fshapeNa 64 : kCarTolerance(rhs.kCarTolerance), fshapeName(rhs.fshapeName) 74 { 65 { 75 // Register to store 66 // Register to store 76 // 67 // 77 G4SolidStore::GetInstance()->Register(this 68 G4SolidStore::GetInstance()->Register(this); 78 } 69 } 79 70 80 ////////////////////////////////////////////// 71 ////////////////////////////////////////////////////////////////////////// 81 // 72 // 82 // Fake default constructor - sets only member 73 // Fake default constructor - sets only member data and allocates memory 83 // for usage restri 74 // for usage restricted to object persistency. 84 // 75 // 85 G4VSolid::G4VSolid( __void__& ) 76 G4VSolid::G4VSolid( __void__& ) 86 : fshapeName("") 77 : fshapeName("") 87 { 78 { 88 // Register to store 79 // Register to store 89 // 80 // 90 G4SolidStore::GetInstance()->Register(this 81 G4SolidStore::GetInstance()->Register(this); 91 } 82 } 92 83 93 ////////////////////////////////////////////// 84 ////////////////////////////////////////////////////////////////////////// 94 // 85 // 95 // Destructor (virtual) 86 // Destructor (virtual) 96 // - Remove ourselves from solid Store 87 // - Remove ourselves from solid Store 97 88 98 G4VSolid::~G4VSolid() 89 G4VSolid::~G4VSolid() 99 { 90 { 100 G4SolidStore::GetInstance()->DeRegister(th 91 G4SolidStore::GetInstance()->DeRegister(this); 101 } 92 } 102 93 103 ////////////////////////////////////////////// 94 ////////////////////////////////////////////////////////////////////////// 104 // 95 // 105 // Assignment operator 96 // Assignment operator 106 97 107 G4VSolid& G4VSolid::operator = (const G4VSolid << 98 G4VSolid& G4VSolid::operator = (const G4VSolid& rhs) 108 { 99 { 109 // Check assignment to self 100 // Check assignment to self 110 // 101 // 111 if (this == &rhs) { return *this; } 102 if (this == &rhs) { return *this; } 112 103 113 // Copy data 104 // Copy data 114 // 105 // 115 kCarTolerance = rhs.kCarTolerance; 106 kCarTolerance = rhs.kCarTolerance; 116 fshapeName = rhs.fshapeName; 107 fshapeName = rhs.fshapeName; 117 108 118 return *this; 109 return *this; 119 } << 110 } 120 << 121 << 122 111 123 ////////////////////////////////////////////// 112 ////////////////////////////////////////////////////////////////////////// 124 // 113 // 125 // Set solid name and notify store of the chan << 114 // Streaming operator dumping solid contents 126 115 127 void G4VSolid::SetName(const G4String& name) << 116 std::ostream& operator<< ( std::ostream& os, const G4VSolid& e ) 128 { 117 { 129 fshapeName = name; << 118 return e.StreamInfo(os); 130 G4SolidStore::GetInstance()->SetMapValid(fal << 131 } 119 } 132 120 133 ////////////////////////////////////////////// 121 ////////////////////////////////////////////////////////////////////////// 134 // 122 // 135 // Throw exception if ComputeDimensions called 123 // Throw exception if ComputeDimensions called for illegal derived class 136 124 137 void G4VSolid::ComputeDimensions(G4VPVParamete 125 void G4VSolid::ComputeDimensions(G4VPVParameterisation*, 138 const G4int, 126 const G4int, 139 const G4VPhys 127 const G4VPhysicalVolume*) 140 { 128 { 141 std::ostringstream message; 129 std::ostringstream message; 142 message << "Illegal call to G4VSolid::Comp 130 message << "Illegal call to G4VSolid::ComputeDimensions()" << G4endl 143 << "Method not overloaded by deriv 131 << "Method not overloaded by derived class !"; 144 G4Exception("G4VSolid::ComputeDimensions() 132 G4Exception("G4VSolid::ComputeDimensions()", "GeomMgt0003", 145 FatalException, message); 133 FatalException, message); 146 } 134 } 147 135 148 ////////////////////////////////////////////// 136 ////////////////////////////////////////////////////////////////////////// 149 // 137 // 150 // Throw exception (warning) for solids not im 138 // Throw exception (warning) for solids not implementing the method 151 139 152 G4ThreeVector G4VSolid::GetPointOnSurface() co 140 G4ThreeVector G4VSolid::GetPointOnSurface() const 153 { 141 { 154 std::ostringstream message; 142 std::ostringstream message; 155 message << "Not implemented for solid: " 143 message << "Not implemented for solid: " 156 << GetEntityType() << " !" << G4en 144 << GetEntityType() << " !" << G4endl 157 << "Returning origin."; 145 << "Returning origin."; 158 G4Exception("G4VSolid::GetPointOnSurface() 146 G4Exception("G4VSolid::GetPointOnSurface()", "GeomMgt1001", 159 JustWarning, message); 147 JustWarning, message); 160 return {0,0,0}; << 148 return G4ThreeVector(0,0,0); 161 } 149 } 162 150 163 ////////////////////////////////////////////// 151 ////////////////////////////////////////////////////////////////////////// 164 // 152 // 165 // Returns total number of constituents that w << 166 // of the solid. For non-Boolean solids the re << 167 << 168 G4int G4VSolid::GetNumOfConstituents() const << 169 { return 1; } << 170 << 171 ////////////////////////////////////////////// << 172 // << 173 // Returns true if the solid has only planar f << 174 << 175 G4bool G4VSolid::IsFaceted() const << 176 { return false; } << 177 << 178 ////////////////////////////////////////////// << 179 // << 180 // Dummy implementations ... 153 // Dummy implementations ... 181 154 182 const G4VSolid* G4VSolid::GetConstituentSolid( 155 const G4VSolid* G4VSolid::GetConstituentSolid(G4int) const 183 { return nullptr; } << 156 { return nullptr; } 184 157 185 G4VSolid* G4VSolid::GetConstituentSolid(G4int) 158 G4VSolid* G4VSolid::GetConstituentSolid(G4int) 186 { return nullptr; } << 159 { return nullptr; } 187 160 188 const G4DisplacedSolid* G4VSolid::GetDisplaced 161 const G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr() const 189 { return nullptr; } << 162 { return nullptr; } 190 163 191 G4DisplacedSolid* G4VSolid::GetDisplacedSolidP << 164 G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr() 192 { return nullptr; } << 165 { return nullptr; } 193 166 194 ////////////////////////////////////////////// 167 //////////////////////////////////////////////////////////////// 195 // 168 // 196 // Returns an estimation of the solid volume i 169 // Returns an estimation of the solid volume in internal units. 197 // The number of statistics and error accuracy 170 // The number of statistics and error accuracy is fixed. 198 // This method may be overloaded by derived cl 171 // This method may be overloaded by derived classes to compute the 199 // exact geometrical quantity for solids where 172 // exact geometrical quantity for solids where this is possible. 200 // or anyway to cache the computed value. 173 // or anyway to cache the computed value. 201 // This implementation does NOT cache the comp 174 // This implementation does NOT cache the computed value. 202 175 203 G4double G4VSolid::GetCubicVolume() 176 G4double G4VSolid::GetCubicVolume() 204 { 177 { 205 G4int cubVolStatistics = 1000000; 178 G4int cubVolStatistics = 1000000; 206 G4double cubVolEpsilon = 0.001; 179 G4double cubVolEpsilon = 0.001; 207 return EstimateCubicVolume(cubVolStatistics, 180 return EstimateCubicVolume(cubVolStatistics, cubVolEpsilon); 208 } 181 } 209 182 210 ////////////////////////////////////////////// 183 //////////////////////////////////////////////////////////////// 211 // 184 // 212 // Calculate cubic volume based on Inside() me 185 // Calculate cubic volume based on Inside() method. 213 // Accuracy is limited by the second argument 186 // Accuracy is limited by the second argument or the statistics 214 // expressed by the first argument. 187 // expressed by the first argument. 215 // Implementation is courtesy of Vasiliki Desp 188 // Implementation is courtesy of Vasiliki Despoina Mitsou, 216 // University of Athens. 189 // University of Athens. 217 190 218 G4double G4VSolid::EstimateCubicVolume(G4int n 191 G4double G4VSolid::EstimateCubicVolume(G4int nStat, G4double epsilon) const 219 { 192 { 220 G4int iInside=0; 193 G4int iInside=0; 221 G4double px,py,pz,minX,maxX,minY,maxY,minZ,m 194 G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,volume,halfepsilon; 222 G4ThreeVector p; 195 G4ThreeVector p; 223 EInside in; 196 EInside in; 224 197 225 // values needed for CalculateExtent signatu 198 // values needed for CalculateExtent signature 226 199 227 G4VoxelLimits limit; // Unlim 200 G4VoxelLimits limit; // Unlimited 228 G4AffineTransform origin; 201 G4AffineTransform origin; 229 202 230 // min max extents of pSolid along X,Y,Z 203 // min max extents of pSolid along X,Y,Z 231 204 232 CalculateExtent(kXAxis,limit,origin,minX,max 205 CalculateExtent(kXAxis,limit,origin,minX,maxX); 233 CalculateExtent(kYAxis,limit,origin,minY,max 206 CalculateExtent(kYAxis,limit,origin,minY,maxY); 234 CalculateExtent(kZAxis,limit,origin,minZ,max 207 CalculateExtent(kZAxis,limit,origin,minZ,maxZ); 235 208 236 // limits 209 // limits 237 210 238 if(nStat < 100) nStat = 100; 211 if(nStat < 100) nStat = 100; 239 if(epsilon > 0.01) epsilon = 0.01; 212 if(epsilon > 0.01) epsilon = 0.01; 240 halfepsilon = 0.5*epsilon; 213 halfepsilon = 0.5*epsilon; 241 214 242 for(auto i = 0; i < nStat; ++i ) 215 for(auto i = 0; i < nStat; ++i ) 243 { 216 { 244 px = minX-halfepsilon+(maxX-minX+epsilon)* 217 px = minX-halfepsilon+(maxX-minX+epsilon)*G4QuickRand(); 245 py = minY-halfepsilon+(maxY-minY+epsilon)* 218 py = minY-halfepsilon+(maxY-minY+epsilon)*G4QuickRand(); 246 pz = minZ-halfepsilon+(maxZ-minZ+epsilon)* 219 pz = minZ-halfepsilon+(maxZ-minZ+epsilon)*G4QuickRand(); 247 p = G4ThreeVector(px,py,pz); 220 p = G4ThreeVector(px,py,pz); 248 in = Inside(p); 221 in = Inside(p); 249 if(in != kOutside) ++iInside; << 222 if(in != kOutside) ++iInside; 250 } 223 } 251 volume = (maxX-minX+epsilon)*(maxY-minY+epsi 224 volume = (maxX-minX+epsilon)*(maxY-minY+epsilon) 252 * (maxZ-minZ+epsilon)*iInside/nStat; 225 * (maxZ-minZ+epsilon)*iInside/nStat; 253 return volume; 226 return volume; 254 } 227 } 255 228 256 ////////////////////////////////////////////// 229 //////////////////////////////////////////////////////////////// 257 // 230 // 258 // Returns an estimation of the solid surface 231 // Returns an estimation of the solid surface area in internal units. 259 // The number of statistics and error accuracy 232 // The number of statistics and error accuracy is fixed. 260 // This method may be overloaded by derived cl 233 // This method may be overloaded by derived classes to compute the 261 // exact geometrical quantity for solids where 234 // exact geometrical quantity for solids where this is possible. 262 // or anyway to cache the computed value. 235 // or anyway to cache the computed value. 263 // This implementation does NOT cache the comp 236 // This implementation does NOT cache the computed value. 264 237 265 G4double G4VSolid::GetSurfaceArea() 238 G4double G4VSolid::GetSurfaceArea() 266 { 239 { 267 G4int stat = 1000000; 240 G4int stat = 1000000; 268 G4double ell = -1.; 241 G4double ell = -1.; 269 return EstimateSurfaceArea(stat,ell); 242 return EstimateSurfaceArea(stat,ell); 270 } 243 } 271 244 272 ////////////////////////////////////////////// 245 ////////////////////////////////////////////////////////////////////////// 273 // 246 // 274 // Calculate surface area by estimating volume 247 // Calculate surface area by estimating volume of a thin shell 275 // surrounding the surface using Monte-Carlo m 248 // surrounding the surface using Monte-Carlo method. 276 // Input parameters: 249 // Input parameters: 277 // nstat - statistics (number of random poi 250 // nstat - statistics (number of random points) 278 // eps - shell thinkness 251 // eps - shell thinkness 279 252 280 G4double G4VSolid::EstimateSurfaceArea(G4int n 253 G4double G4VSolid::EstimateSurfaceArea(G4int nstat, G4double ell) const 281 { 254 { 282 static const G4double s2 = 1./std::sqrt(2.); 255 static const G4double s2 = 1./std::sqrt(2.); 283 static const G4double s3 = 1./std::sqrt(3.); 256 static const G4double s3 = 1./std::sqrt(3.); 284 static const G4ThreeVector directions[64] = 257 static const G4ThreeVector directions[64] = 285 { 258 { 286 G4ThreeVector( 0, 0, 0), G4ThreeVector( 259 G4ThreeVector( 0, 0, 0), G4ThreeVector( -1, 0, 0), // ( , , ) ( -, , ) 287 G4ThreeVector( 1, 0, 0), G4ThreeVector( 260 G4ThreeVector( 1, 0, 0), G4ThreeVector( -1, 0, 0), // ( +, , ) (-+, , ) 288 G4ThreeVector( 0, -1, 0), G4ThreeVector( 261 G4ThreeVector( 0, -1, 0), G4ThreeVector(-s2,-s2, 0), // ( , -, ) ( -, -, ) 289 G4ThreeVector( s2, -s2, 0), G4ThreeVector( 262 G4ThreeVector( s2, -s2, 0), G4ThreeVector( 0, -1, 0), // ( +, -, ) (-+, -, ) 290 263 291 G4ThreeVector( 0, 1, 0), G4ThreeVector( 264 G4ThreeVector( 0, 1, 0), G4ThreeVector( -s2, s2, 0), // ( , +, ) ( -, +, ) 292 G4ThreeVector( s2, s2, 0), G4ThreeVector( 265 G4ThreeVector( s2, s2, 0), G4ThreeVector( 0, 1, 0), // ( +, +, ) (-+, +, ) 293 G4ThreeVector( 0, -1, 0), G4ThreeVector( 266 G4ThreeVector( 0, -1, 0), G4ThreeVector( -1, 0, 0), // ( ,-+, ) ( -,-+, ) 294 G4ThreeVector( 1, 0, 0), G4ThreeVector( 267 G4ThreeVector( 1, 0, 0), G4ThreeVector( -1, 0, 0), // ( +,-+, ) (-+,-+, ) 295 268 296 G4ThreeVector( 0, 0, -1), G4ThreeVector( 269 G4ThreeVector( 0, 0, -1), G4ThreeVector(-s2, 0,-s2), // ( , , -) ( -, , -) 297 G4ThreeVector( s2, 0,-s2), G4ThreeVector( 270 G4ThreeVector( s2, 0,-s2), G4ThreeVector( 0, 0, -1), // ( +, , -) (-+, , -) 298 G4ThreeVector( 0,-s2,-s2), G4ThreeVector( 271 G4ThreeVector( 0,-s2,-s2), G4ThreeVector(-s3,-s3,-s3), // ( , -, -) ( -, -, -) 299 G4ThreeVector( s3,-s3,-s3), G4ThreeVector( 272 G4ThreeVector( s3,-s3,-s3), G4ThreeVector( 0,-s2,-s2), // ( +, -, -) (-+, -, -) 300 273 301 G4ThreeVector( 0, s2,-s2), G4ThreeVector( 274 G4ThreeVector( 0, s2,-s2), G4ThreeVector(-s3, s3,-s3), // ( , +, -) ( -, +, -) 302 G4ThreeVector( s3, s3,-s3), G4ThreeVector( 275 G4ThreeVector( s3, s3,-s3), G4ThreeVector( 0, s2,-s2), // ( +, +, -) (-+, +, -) 303 G4ThreeVector( 0, 0, -1), G4ThreeVector( 276 G4ThreeVector( 0, 0, -1), G4ThreeVector(-s2, 0,-s2), // ( ,-+, -) ( -,-+, -) 304 G4ThreeVector( s2, 0,-s2), G4ThreeVector( 277 G4ThreeVector( s2, 0,-s2), G4ThreeVector( 0, 0, -1), // ( +,-+, -) (-+,-+, -) 305 278 306 G4ThreeVector( 0, 0, 1), G4ThreeVector( 279 G4ThreeVector( 0, 0, 1), G4ThreeVector(-s2, 0, s2), // ( , , +) ( -, , +) 307 G4ThreeVector( s2, 0, s2), G4ThreeVector( 280 G4ThreeVector( s2, 0, s2), G4ThreeVector( 0, 0, 1), // ( +, , +) (-+, , +) 308 G4ThreeVector( 0,-s2, s2), G4ThreeVector( 281 G4ThreeVector( 0,-s2, s2), G4ThreeVector(-s3,-s3, s3), // ( , -, +) ( -, -, +) 309 G4ThreeVector( s3,-s3, s3), G4ThreeVector( 282 G4ThreeVector( s3,-s3, s3), G4ThreeVector( 0,-s2, s2), // ( +, -, +) (-+, -, +) 310 283 311 G4ThreeVector( 0, s2, s2), G4ThreeVector( 284 G4ThreeVector( 0, s2, s2), G4ThreeVector(-s3, s3, s3), // ( , +, +) ( -, +, +) 312 G4ThreeVector( s3, s3, s3), G4ThreeVector( 285 G4ThreeVector( s3, s3, s3), G4ThreeVector( 0, s2, s2), // ( +, +, +) (-+, +, +) 313 G4ThreeVector( 0, 0, 1), G4ThreeVector( 286 G4ThreeVector( 0, 0, 1), G4ThreeVector(-s2, 0, s2), // ( ,-+, +) ( -,-+, +) 314 G4ThreeVector( s2, 0, s2), G4ThreeVector( 287 G4ThreeVector( s2, 0, s2), G4ThreeVector( 0, 0, 1), // ( +,-+, +) (-+,-+, +) 315 288 316 G4ThreeVector( 0, 0, -1), G4ThreeVector( 289 G4ThreeVector( 0, 0, -1), G4ThreeVector( -1, 0, 0), // ( , ,-+) ( -, ,-+) 317 G4ThreeVector( 1, 0, 0), G4ThreeVector( 290 G4ThreeVector( 1, 0, 0), G4ThreeVector( -1, 0, 0), // ( +, ,-+) (-+, ,-+) 318 G4ThreeVector( 0, -1, 0), G4ThreeVector( 291 G4ThreeVector( 0, -1, 0), G4ThreeVector(-s2,-s2, 0), // ( , -,-+) ( -, -,-+) 319 G4ThreeVector( s2, -s2, 0), G4ThreeVector( 292 G4ThreeVector( s2, -s2, 0), G4ThreeVector( 0, -1, 0), // ( +, -,-+) (-+, -,-+) 320 293 321 G4ThreeVector( 0, 1, 0), G4ThreeVector( 294 G4ThreeVector( 0, 1, 0), G4ThreeVector( -s2, s2, 0), // ( , +,-+) ( -, +,-+) 322 G4ThreeVector( s2, s2, 0), G4ThreeVector( 295 G4ThreeVector( s2, s2, 0), G4ThreeVector( 0, 1, 0), // ( +, +,-+) (-+, +,-+) 323 G4ThreeVector( 0, -1, 0), G4ThreeVector( 296 G4ThreeVector( 0, -1, 0), G4ThreeVector( -1, 0, 0), // ( ,-+,-+) ( -,-+,-+) 324 G4ThreeVector( 1, 0, 0), G4ThreeVector( 297 G4ThreeVector( 1, 0, 0), G4ThreeVector( -1, 0, 0), // ( +,-+,-+) (-+,-+,-+) 325 }; 298 }; 326 299 327 G4ThreeVector bmin, bmax; 300 G4ThreeVector bmin, bmax; 328 BoundingLimits(bmin, bmax); 301 BoundingLimits(bmin, bmax); 329 302 330 G4double dX = bmax.x() - bmin.x(); 303 G4double dX = bmax.x() - bmin.x(); 331 G4double dY = bmax.y() - bmin.y(); 304 G4double dY = bmax.y() - bmin.y(); 332 G4double dZ = bmax.z() - bmin.z(); 305 G4double dZ = bmax.z() - bmin.z(); 333 << 306 334 // Define statistics and shell thickness 307 // Define statistics and shell thickness 335 // 308 // 336 G4int npoints = (nstat < 1000) ? 1000 : nsta 309 G4int npoints = (nstat < 1000) ? 1000 : nstat; 337 G4double coeff = 0.5 / std::cbrt(G4double(np 310 G4double coeff = 0.5 / std::cbrt(G4double(npoints)); 338 G4double eps = (ell > 0) ? ell : coeff * std 311 G4double eps = (ell > 0) ? ell : coeff * std::min(std::min(dX, dY), dZ); 339 G4double del = 1.8 * eps; // shold be more t 312 G4double del = 1.8 * eps; // shold be more than sqrt(3.) 340 313 341 G4double minX = bmin.x() - eps; 314 G4double minX = bmin.x() - eps; 342 G4double minY = bmin.y() - eps; 315 G4double minY = bmin.y() - eps; 343 G4double minZ = bmin.z() - eps; 316 G4double minZ = bmin.z() - eps; 344 317 345 G4double dd = 2. * eps; 318 G4double dd = 2. * eps; 346 dX += dd; 319 dX += dd; 347 dY += dd; 320 dY += dd; 348 dZ += dd; 321 dZ += dd; 349 322 350 // Calculate surface area 323 // Calculate surface area 351 // 324 // 352 G4int icount = 0; 325 G4int icount = 0; 353 for(auto i = 0; i < npoints; ++i) 326 for(auto i = 0; i < npoints; ++i) 354 { 327 { 355 G4double px = minX + dX*G4QuickRand(); 328 G4double px = minX + dX*G4QuickRand(); 356 G4double py = minY + dY*G4QuickRand(); 329 G4double py = minY + dY*G4QuickRand(); 357 G4double pz = minZ + dZ*G4QuickRand(); 330 G4double pz = minZ + dZ*G4QuickRand(); 358 G4ThreeVector p = G4ThreeVector(px, py, p 331 G4ThreeVector p = G4ThreeVector(px, py, pz); 359 EInside in = Inside(p); 332 EInside in = Inside(p); 360 G4double dist = 0; 333 G4double dist = 0; 361 if (in == kInside) 334 if (in == kInside) 362 { 335 { 363 if (DistanceToOut(p) >= eps) continue; 336 if (DistanceToOut(p) >= eps) continue; 364 G4int icase = 0; 337 G4int icase = 0; 365 if (Inside(G4ThreeVector(px-del, py, pz) 338 if (Inside(G4ThreeVector(px-del, py, pz)) != kInside) icase += 1; 366 if (Inside(G4ThreeVector(px+del, py, pz) 339 if (Inside(G4ThreeVector(px+del, py, pz)) != kInside) icase += 2; 367 if (Inside(G4ThreeVector(px, py-del, pz) 340 if (Inside(G4ThreeVector(px, py-del, pz)) != kInside) icase += 4; 368 if (Inside(G4ThreeVector(px, py+del, pz) 341 if (Inside(G4ThreeVector(px, py+del, pz)) != kInside) icase += 8; 369 if (Inside(G4ThreeVector(px, py, pz-del) 342 if (Inside(G4ThreeVector(px, py, pz-del)) != kInside) icase += 16; 370 if (Inside(G4ThreeVector(px, py, pz+del) 343 if (Inside(G4ThreeVector(px, py, pz+del)) != kInside) icase += 32; 371 if (icase == 0) continue; 344 if (icase == 0) continue; 372 G4ThreeVector v = directions[icase]; 345 G4ThreeVector v = directions[icase]; 373 dist = DistanceToOut(p, v); 346 dist = DistanceToOut(p, v); 374 G4ThreeVector n = SurfaceNormal(p + v*di 347 G4ThreeVector n = SurfaceNormal(p + v*dist); 375 dist *= v.dot(n); 348 dist *= v.dot(n); 376 } 349 } 377 else if (in == kOutside) 350 else if (in == kOutside) 378 { 351 { 379 if (DistanceToIn(p) >= eps) continue; 352 if (DistanceToIn(p) >= eps) continue; 380 G4int icase = 0; 353 G4int icase = 0; 381 if (Inside(G4ThreeVector(px-del, py, pz) 354 if (Inside(G4ThreeVector(px-del, py, pz)) != kOutside) icase += 1; 382 if (Inside(G4ThreeVector(px+del, py, pz) 355 if (Inside(G4ThreeVector(px+del, py, pz)) != kOutside) icase += 2; 383 if (Inside(G4ThreeVector(px, py-del, pz) 356 if (Inside(G4ThreeVector(px, py-del, pz)) != kOutside) icase += 4; 384 if (Inside(G4ThreeVector(px, py+del, pz) 357 if (Inside(G4ThreeVector(px, py+del, pz)) != kOutside) icase += 8; 385 if (Inside(G4ThreeVector(px, py, pz-del) 358 if (Inside(G4ThreeVector(px, py, pz-del)) != kOutside) icase += 16; 386 if (Inside(G4ThreeVector(px, py, pz+del) 359 if (Inside(G4ThreeVector(px, py, pz+del)) != kOutside) icase += 32; 387 if (icase == 0) continue; 360 if (icase == 0) continue; 388 G4ThreeVector v = directions[icase]; 361 G4ThreeVector v = directions[icase]; 389 dist = DistanceToIn(p, v); 362 dist = DistanceToIn(p, v); 390 if (dist == kInfinity) continue; 363 if (dist == kInfinity) continue; 391 G4ThreeVector n = SurfaceNormal(p + v*di 364 G4ThreeVector n = SurfaceNormal(p + v*dist); 392 dist *= -(v.dot(n)); 365 dist *= -(v.dot(n)); 393 } 366 } 394 if (dist < eps) ++icount; 367 if (dist < eps) ++icount; 395 } 368 } 396 return dX*dY*dZ*icount/npoints/dd; 369 return dX*dY*dZ*icount/npoints/dd; 397 } 370 } 398 371 399 ////////////////////////////////////////////// 372 /////////////////////////////////////////////////////////////////////////// 400 // << 373 // 401 // Returns a pointer of a dynamically allocate 374 // Returns a pointer of a dynamically allocated copy of the solid. 402 // Returns NULL pointer with warning in case t 375 // Returns NULL pointer with warning in case the concrete solid does not 403 // implement this method. The caller has respo 376 // implement this method. The caller has responsibility for ownership. 404 // 377 // 405 378 406 G4VSolid* G4VSolid::Clone() const 379 G4VSolid* G4VSolid::Clone() const 407 { 380 { 408 std::ostringstream message; 381 std::ostringstream message; 409 message << "Clone() method not implemented f 382 message << "Clone() method not implemented for type: " 410 << GetEntityType() << "!" << G4endl 383 << GetEntityType() << "!" << G4endl 411 << "Returning NULL pointer!"; 384 << "Returning NULL pointer!"; 412 G4Exception("G4VSolid::Clone()", "GeomMgt100 385 G4Exception("G4VSolid::Clone()", "GeomMgt1001", JustWarning, message); 413 return nullptr; 386 return nullptr; 414 } 387 } 415 388 416 ////////////////////////////////////////////// 389 /////////////////////////////////////////////////////////////////////////// 417 // << 390 // 418 // Calculate the maximum and minimum extents o 391 // Calculate the maximum and minimum extents of the polygon described 419 // by the vertices: pSectionIndex->pSectionInd 392 // by the vertices: pSectionIndex->pSectionIndex+1-> 420 // pSectionIndex+2->pSection 393 // pSectionIndex+2->pSectionIndex+3->pSectionIndex 421 // in the List pVertices 394 // in the List pVertices 422 // 395 // 423 // If the minimum is <pMin pMin is set to the 396 // If the minimum is <pMin pMin is set to the new minimum 424 // If the maximum is >pMax pMax is set to the 397 // If the maximum is >pMax pMax is set to the new maximum 425 // 398 // 426 // No modifications are made to pVertices 399 // No modifications are made to pVertices 427 // 400 // 428 401 429 void G4VSolid::ClipCrossSection( G4Three 402 void G4VSolid::ClipCrossSection( G4ThreeVectorList* pVertices, 430 const G4int p 403 const G4int pSectionIndex, 431 const G4Voxel 404 const G4VoxelLimits& pVoxelLimit, 432 const EAxis p << 405 const EAxis pAxis, 433 G4doubl 406 G4double& pMin, G4double& pMax) const 434 { 407 { 435 408 436 G4ThreeVectorList polygon; 409 G4ThreeVectorList polygon; 437 polygon.reserve(4); 410 polygon.reserve(4); 438 polygon.push_back((*pVertices)[pSectionIndex 411 polygon.push_back((*pVertices)[pSectionIndex]); 439 polygon.push_back((*pVertices)[pSectionIndex 412 polygon.push_back((*pVertices)[pSectionIndex+1]); 440 polygon.push_back((*pVertices)[pSectionIndex 413 polygon.push_back((*pVertices)[pSectionIndex+2]); 441 polygon.push_back((*pVertices)[pSectionIndex 414 polygon.push_back((*pVertices)[pSectionIndex+3]); 442 CalculateClippedPolygonExtent(polygon,pVoxel 415 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax); 443 return; 416 return; 444 } 417 } 445 418 446 ////////////////////////////////////////////// 419 ////////////////////////////////////////////////////////////////////////////////// 447 // 420 // 448 // Calculate the maximum and minimum extents o 421 // Calculate the maximum and minimum extents of the polygons 449 // joining the CrossSections at pSectionIndex- 422 // joining the CrossSections at pSectionIndex->pSectionIndex+3 and 450 // pSectionIndex+ 423 // pSectionIndex+4->pSectionIndex7 451 // 424 // 452 // in the List pVertices, within the boundarie 425 // in the List pVertices, within the boundaries of the voxel limits pVoxelLimit 453 // 426 // 454 // If the minimum is <pMin pMin is set to the 427 // If the minimum is <pMin pMin is set to the new minimum 455 // If the maximum is >pMax pMax is set to the 428 // If the maximum is >pMax pMax is set to the new maximum 456 // 429 // 457 // No modifications are made to pVertices 430 // No modifications are made to pVertices 458 431 459 void G4VSolid::ClipBetweenSections( G4Thr 432 void G4VSolid::ClipBetweenSections( G4ThreeVectorList* pVertices, 460 const G4int 433 const G4int pSectionIndex, 461 const G4Vox 434 const G4VoxelLimits& pVoxelLimit, 462 const EAxis << 435 const EAxis pAxis, 463 G4dou 436 G4double& pMin, G4double& pMax) const 464 { 437 { 465 G4ThreeVectorList polygon; 438 G4ThreeVectorList polygon; 466 polygon.reserve(4); 439 polygon.reserve(4); 467 polygon.push_back((*pVertices)[pSectionIndex 440 polygon.push_back((*pVertices)[pSectionIndex]); 468 polygon.push_back((*pVertices)[pSectionIndex 441 polygon.push_back((*pVertices)[pSectionIndex+4]); 469 polygon.push_back((*pVertices)[pSectionIndex 442 polygon.push_back((*pVertices)[pSectionIndex+5]); 470 polygon.push_back((*pVertices)[pSectionIndex 443 polygon.push_back((*pVertices)[pSectionIndex+1]); 471 CalculateClippedPolygonExtent(polygon,pVoxel 444 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax); 472 polygon.clear(); 445 polygon.clear(); 473 446 474 polygon.push_back((*pVertices)[pSectionIndex 447 polygon.push_back((*pVertices)[pSectionIndex+1]); 475 polygon.push_back((*pVertices)[pSectionIndex 448 polygon.push_back((*pVertices)[pSectionIndex+5]); 476 polygon.push_back((*pVertices)[pSectionIndex 449 polygon.push_back((*pVertices)[pSectionIndex+6]); 477 polygon.push_back((*pVertices)[pSectionIndex 450 polygon.push_back((*pVertices)[pSectionIndex+2]); 478 CalculateClippedPolygonExtent(polygon,pVoxel 451 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax); 479 polygon.clear(); 452 polygon.clear(); 480 453 481 polygon.push_back((*pVertices)[pSectionIndex 454 polygon.push_back((*pVertices)[pSectionIndex+2]); 482 polygon.push_back((*pVertices)[pSectionIndex 455 polygon.push_back((*pVertices)[pSectionIndex+6]); 483 polygon.push_back((*pVertices)[pSectionIndex 456 polygon.push_back((*pVertices)[pSectionIndex+7]); 484 polygon.push_back((*pVertices)[pSectionIndex 457 polygon.push_back((*pVertices)[pSectionIndex+3]); 485 CalculateClippedPolygonExtent(polygon,pVoxel 458 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax); 486 polygon.clear(); 459 polygon.clear(); 487 460 488 polygon.push_back((*pVertices)[pSectionIndex 461 polygon.push_back((*pVertices)[pSectionIndex+3]); 489 polygon.push_back((*pVertices)[pSectionIndex 462 polygon.push_back((*pVertices)[pSectionIndex+7]); 490 polygon.push_back((*pVertices)[pSectionIndex 463 polygon.push_back((*pVertices)[pSectionIndex+4]); 491 polygon.push_back((*pVertices)[pSectionIndex 464 polygon.push_back((*pVertices)[pSectionIndex]); 492 CalculateClippedPolygonExtent(polygon,pVoxel 465 CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax); 493 return; 466 return; 494 } 467 } 495 468 496 469 497 ////////////////////////////////////////////// 470 /////////////////////////////////////////////////////////////////////////////// 498 // 471 // 499 // Calculate the maximum and minimum extents o 472 // Calculate the maximum and minimum extents of the convex polygon pPolygon 500 // along the axis pAxis, within the limits pVo 473 // along the axis pAxis, within the limits pVoxelLimit 501 // 474 // 502 475 503 void 476 void 504 G4VSolid::CalculateClippedPolygonExtent(G4Thre 477 G4VSolid::CalculateClippedPolygonExtent(G4ThreeVectorList& pPolygon, 505 const G4Voxe 478 const G4VoxelLimits& pVoxelLimit, 506 const EAxis << 479 const EAxis pAxis, 507 G4doub 480 G4double& pMin, 508 G4doub 481 G4double& pMax) const 509 { 482 { 510 G4int noLeft,i; 483 G4int noLeft,i; 511 G4double component; 484 G4double component; 512 485 513 ClipPolygon(pPolygon,pVoxelLimit,pAxis); 486 ClipPolygon(pPolygon,pVoxelLimit,pAxis); 514 noLeft = (G4int)pPolygon.size(); << 487 noLeft = pPolygon.size(); 515 488 516 if ( noLeft != 0 ) << 489 if ( noLeft ) 517 { 490 { 518 for (i=0; i<noLeft; ++i) 491 for (i=0; i<noLeft; ++i) 519 { 492 { 520 component = pPolygon[i].operator()(pAxis 493 component = pPolygon[i].operator()(pAxis); 521 << 494 522 if (component < pMin) << 495 if (component < pMin) 523 { << 496 { 524 pMin = component; << 497 pMin = component; 525 } 498 } 526 if (component > pMax) 499 if (component > pMax) 527 { << 500 { 528 pMax = component; << 501 pMax = component; 529 } << 502 } 530 } 503 } 531 } 504 } 532 } 505 } 533 506 534 ////////////////////////////////////////////// 507 ///////////////////////////////////////////////////////////////////////////// 535 // 508 // 536 // Clip the convex polygon described by the ve 509 // Clip the convex polygon described by the vertices at 537 // pSectionIndex ->pSectionIndex+3 within pVer 510 // pSectionIndex ->pSectionIndex+3 within pVertices to the limits pVoxelLimit 538 // 511 // 539 // Set pMin to the smallest 512 // Set pMin to the smallest 540 // 513 // 541 // Calculate the extent of the polygon along p 514 // Calculate the extent of the polygon along pAxis, when clipped to the 542 // limits pVoxelLimit. If the polygon exists a 515 // limits pVoxelLimit. If the polygon exists after clippin, set pMin to 543 // the polygon's minimum extent along the axis 516 // the polygon's minimum extent along the axis if <pMin, and set pMax to 544 // the polygon's maximum extent along the axis 517 // the polygon's maximum extent along the axis if >pMax. 545 // 518 // 546 // The polygon is described by a set of vector 519 // The polygon is described by a set of vectors, where each vector represents 547 // a vertex, so that the polygon is described 520 // a vertex, so that the polygon is described by the vertex sequence: 548 // 0th->1st 1st->2nd 2nd->... nth->0th 521 // 0th->1st 1st->2nd 2nd->... nth->0th 549 // 522 // 550 // Modifications to the polygon are made 523 // Modifications to the polygon are made 551 // 524 // 552 // NOTE: Execessive copying during clipping 525 // NOTE: Execessive copying during clipping 553 526 554 void G4VSolid::ClipPolygon( G4ThreeVector 527 void G4VSolid::ClipPolygon( G4ThreeVectorList& pPolygon, 555 const G4VoxelLimits 528 const G4VoxelLimits& pVoxelLimit, 556 const EAxis 529 const EAxis ) const 557 { 530 { 558 G4ThreeVectorList outputPolygon; 531 G4ThreeVectorList outputPolygon; 559 532 560 if ( pVoxelLimit.IsLimited() ) 533 if ( pVoxelLimit.IsLimited() ) 561 { 534 { 562 if (pVoxelLimit.IsXLimited() ) // && pAxis 535 if (pVoxelLimit.IsXLimited() ) // && pAxis != kXAxis) 563 { 536 { 564 G4VoxelLimits simpleLimit1; 537 G4VoxelLimits simpleLimit1; 565 simpleLimit1.AddLimit(kXAxis,pVoxelLimit 538 simpleLimit1.AddLimit(kXAxis,pVoxelLimit.GetMinXExtent(),kInfinity); 566 ClipPolygonToSimpleLimits(pPolygon,outpu 539 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1); 567 << 540 568 pPolygon.clear(); 541 pPolygon.clear(); 569 542 570 if ( outputPolygon.empty() ) return; << 543 if ( !outputPolygon.size() ) return; 571 544 572 G4VoxelLimits simpleLimit2; 545 G4VoxelLimits simpleLimit2; 573 simpleLimit2.AddLimit(kXAxis,-kInfinity, 546 simpleLimit2.AddLimit(kXAxis,-kInfinity,pVoxelLimit.GetMaxXExtent()); 574 ClipPolygonToSimpleLimits(outputPolygon, 547 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2); 575 548 576 if ( pPolygon.empty() ) return; << 549 if ( !pPolygon.size() ) return; 577 else outputPoly 550 else outputPolygon.clear(); 578 } 551 } 579 if ( pVoxelLimit.IsYLimited() ) // && pAxi 552 if ( pVoxelLimit.IsYLimited() ) // && pAxis != kYAxis) 580 { 553 { 581 G4VoxelLimits simpleLimit1; 554 G4VoxelLimits simpleLimit1; 582 simpleLimit1.AddLimit(kYAxis,pVoxelLimit 555 simpleLimit1.AddLimit(kYAxis,pVoxelLimit.GetMinYExtent(),kInfinity); 583 ClipPolygonToSimpleLimits(pPolygon,outpu 556 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1); 584 557 585 // Must always clear pPolygon - for clip 558 // Must always clear pPolygon - for clip to simpleLimit2 and in case of 586 // early exit 559 // early exit 587 560 588 pPolygon.clear(); 561 pPolygon.clear(); 589 562 590 if ( outputPolygon.empty() ) return; << 563 if ( !outputPolygon.size() ) return; 591 564 592 G4VoxelLimits simpleLimit2; 565 G4VoxelLimits simpleLimit2; 593 simpleLimit2.AddLimit(kYAxis,-kInfinity, 566 simpleLimit2.AddLimit(kYAxis,-kInfinity,pVoxelLimit.GetMaxYExtent()); 594 ClipPolygonToSimpleLimits(outputPolygon, 567 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2); 595 568 596 if ( pPolygon.empty() ) return; << 569 if ( !pPolygon.size() ) return; 597 else outputPoly 570 else outputPolygon.clear(); 598 } 571 } 599 if ( pVoxelLimit.IsZLimited() ) // && pAxi 572 if ( pVoxelLimit.IsZLimited() ) // && pAxis != kZAxis) 600 { 573 { 601 G4VoxelLimits simpleLimit1; 574 G4VoxelLimits simpleLimit1; 602 simpleLimit1.AddLimit(kZAxis,pVoxelLimit 575 simpleLimit1.AddLimit(kZAxis,pVoxelLimit.GetMinZExtent(),kInfinity); 603 ClipPolygonToSimpleLimits(pPolygon,outpu 576 ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1); 604 577 605 // Must always clear pPolygon - for clip 578 // Must always clear pPolygon - for clip to simpleLimit2 and in case of 606 // early exit 579 // early exit 607 580 608 pPolygon.clear(); 581 pPolygon.clear(); 609 582 610 if ( outputPolygon.empty() ) return; << 583 if ( !outputPolygon.size() ) return; 611 584 612 G4VoxelLimits simpleLimit2; 585 G4VoxelLimits simpleLimit2; 613 simpleLimit2.AddLimit(kZAxis,-kInfinity, 586 simpleLimit2.AddLimit(kZAxis,-kInfinity,pVoxelLimit.GetMaxZExtent()); 614 ClipPolygonToSimpleLimits(outputPolygon, 587 ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2); 615 588 616 // Return after final clip - no cleanup 589 // Return after final clip - no cleanup 617 } 590 } 618 } 591 } 619 } 592 } 620 593 621 ////////////////////////////////////////////// 594 //////////////////////////////////////////////////////////////////////////// 622 // 595 // 623 // pVoxelLimits must be only limited along one 596 // pVoxelLimits must be only limited along one axis, and either the maximum 624 // along the axis must be +kInfinity, or the m 597 // along the axis must be +kInfinity, or the minimum -kInfinity 625 598 626 void 599 void 627 G4VSolid::ClipPolygonToSimpleLimits( G4ThreeVe 600 G4VSolid::ClipPolygonToSimpleLimits( G4ThreeVectorList& pPolygon, 628 G4ThreeVe 601 G4ThreeVectorList& outputPolygon, 629 const G4VoxelLi 602 const G4VoxelLimits& pVoxelLimit ) const 630 { 603 { 631 G4int i; 604 G4int i; 632 auto noVertices = (G4int)pPolygon.size(); << 605 G4int noVertices=pPolygon.size(); 633 G4ThreeVector vEnd,vStart; 606 G4ThreeVector vEnd,vStart; 634 607 635 for (i = 0 ; i < noVertices ; ++i ) 608 for (i = 0 ; i < noVertices ; ++i ) 636 { 609 { 637 vStart = pPolygon[i]; 610 vStart = pPolygon[i]; 638 if ( i == noVertices-1 ) vEnd = pPolygo 611 if ( i == noVertices-1 ) vEnd = pPolygon[0]; 639 else vEnd = pPolygo 612 else vEnd = pPolygon[i+1]; 640 613 641 if ( pVoxelLimit.Inside(vStart) ) 614 if ( pVoxelLimit.Inside(vStart) ) 642 { 615 { 643 if (pVoxelLimit.Inside(vEnd)) 616 if (pVoxelLimit.Inside(vEnd)) 644 { 617 { 645 // vStart and vEnd inside -> output en 618 // vStart and vEnd inside -> output end point 646 // 619 // 647 outputPolygon.push_back(vEnd); 620 outputPolygon.push_back(vEnd); 648 } 621 } 649 else 622 else 650 { 623 { 651 // vStart inside, vEnd outside -> outp 624 // vStart inside, vEnd outside -> output crossing point 652 // 625 // 653 pVoxelLimit.ClipToLimits(vStart,vEnd); 626 pVoxelLimit.ClipToLimits(vStart,vEnd); 654 outputPolygon.push_back(vEnd); 627 outputPolygon.push_back(vEnd); 655 } << 628 } 656 } 629 } 657 else 630 else 658 { 631 { 659 if (pVoxelLimit.Inside(vEnd)) 632 if (pVoxelLimit.Inside(vEnd)) 660 { 633 { 661 // vStart outside, vEnd inside -> outp 634 // vStart outside, vEnd inside -> output inside section 662 // 635 // 663 pVoxelLimit.ClipToLimits(vStart,vEnd); 636 pVoxelLimit.ClipToLimits(vStart,vEnd); 664 outputPolygon.push_back(vStart); 637 outputPolygon.push_back(vStart); 665 outputPolygon.push_back(vEnd); << 638 outputPolygon.push_back(vEnd); 666 } 639 } 667 else // Both point outside -> no output 640 else // Both point outside -> no output 668 { 641 { 669 // outputPolygon.push_back(vStart); 642 // outputPolygon.push_back(vStart); 670 // outputPolygon.push_back(vEnd); << 643 // outputPolygon.push_back(vEnd); 671 } 644 } 672 } 645 } 673 } 646 } 674 } 647 } 675 648 676 ////////////////////////////////////////////// 649 ////////////////////////////////////////////////////////////////////////// 677 // 650 // 678 // Throw exception (warning) for solids not im 651 // Throw exception (warning) for solids not implementing the method 679 652 680 void G4VSolid::BoundingLimits(G4ThreeVector& p 653 void G4VSolid::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const 681 { 654 { 682 std::ostringstream message; 655 std::ostringstream message; 683 message << "Not implemented for solid: " 656 message << "Not implemented for solid: " 684 << GetEntityType() << " !" 657 << GetEntityType() << " !" 685 << "\nReturning infinite boundinx bo 658 << "\nReturning infinite boundinx box."; 686 G4Exception("G4VSolid::BoundingLimits()", "G 659 G4Exception("G4VSolid::BoundingLimits()", "GeomMgt1001", 687 JustWarning, message); 660 JustWarning, message); 688 661 689 pMin.set(-kInfinity,-kInfinity,-kInfinity); 662 pMin.set(-kInfinity,-kInfinity,-kInfinity); 690 pMax.set( kInfinity, kInfinity, kInfinity); 663 pMax.set( kInfinity, kInfinity, kInfinity); 691 } 664 } 692 665 693 ////////////////////////////////////////////// 666 ////////////////////////////////////////////////////////////////////////// 694 // 667 // 695 // Get G4VisExtent - bounding box for graphics 668 // Get G4VisExtent - bounding box for graphics 696 669 697 G4VisExtent G4VSolid::GetExtent () const << 670 G4VisExtent G4VSolid::GetExtent () const 698 { 671 { 699 G4VisExtent extent; 672 G4VisExtent extent; 700 G4VoxelLimits voxelLimits; // Defaults to " 673 G4VoxelLimits voxelLimits; // Defaults to "infinite" limits. 701 G4AffineTransform affineTransform; 674 G4AffineTransform affineTransform; 702 G4double vmin, vmax; 675 G4double vmin, vmax; 703 CalculateExtent(kXAxis,voxelLimits,affineTra 676 CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax); 704 extent.SetXmin (vmin); 677 extent.SetXmin (vmin); 705 extent.SetXmax (vmax); 678 extent.SetXmax (vmax); 706 CalculateExtent(kYAxis,voxelLimits,affineTra 679 CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax); 707 extent.SetYmin (vmin); 680 extent.SetYmin (vmin); 708 extent.SetYmax (vmax); 681 extent.SetYmax (vmax); 709 CalculateExtent(kZAxis,voxelLimits,affineTra 682 CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax); 710 extent.SetZmin (vmin); 683 extent.SetZmin (vmin); 711 extent.SetZmax (vmax); 684 extent.SetZmax (vmax); 712 return extent; 685 return extent; 713 } 686 } 714 687 715 G4Polyhedron* G4VSolid::CreatePolyhedron () co 688 G4Polyhedron* G4VSolid::CreatePolyhedron () const 716 { 689 { 717 return nullptr; 690 return nullptr; 718 } 691 } 719 692 720 G4Polyhedron* G4VSolid::GetPolyhedron () const 693 G4Polyhedron* G4VSolid::GetPolyhedron () const 721 { 694 { 722 return nullptr; 695 return nullptr; 723 } 696 } 724 697