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 // G4EllipticalTube implementation 26 // G4EllipticalTube implementation 27 // 27 // 28 // Author: David C. Williams (davidw@scipp.ucs 28 // Author: David C. Williams (davidw@scipp.ucsc.edu) 29 // Revision: Evgueni Tcherniaev (evgueni.tcher 29 // Revision: Evgueni Tcherniaev (evgueni.tcherniaev@cern.ch), 23.12.2019 30 // ------------------------------------------- 30 // -------------------------------------------------------------------- 31 31 32 #include "G4EllipticalTube.hh" 32 #include "G4EllipticalTube.hh" 33 33 34 #if !(defined(G4GEOM_USE_UELLIPTICALTUBE) && d 34 #if !(defined(G4GEOM_USE_UELLIPTICALTUBE) && defined(G4GEOM_USE_SYS_USOLIDS)) 35 35 36 #include "G4GeomTools.hh" 36 #include "G4GeomTools.hh" 37 #include "G4RandomTools.hh" 37 #include "G4RandomTools.hh" 38 #include "G4ClippablePolygon.hh" 38 #include "G4ClippablePolygon.hh" 39 #include "G4AffineTransform.hh" 39 #include "G4AffineTransform.hh" 40 #include "G4VoxelLimits.hh" 40 #include "G4VoxelLimits.hh" 41 #include "G4BoundingEnvelope.hh" 41 #include "G4BoundingEnvelope.hh" 42 42 43 #include "Randomize.hh" 43 #include "Randomize.hh" 44 44 45 #include "G4VGraphicsScene.hh" 45 #include "G4VGraphicsScene.hh" 46 #include "G4VisExtent.hh" 46 #include "G4VisExtent.hh" 47 47 48 #include "G4AutoLock.hh" 48 #include "G4AutoLock.hh" 49 49 50 namespace 50 namespace 51 { 51 { 52 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE 52 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER; 53 } 53 } 54 54 55 using namespace CLHEP; 55 using namespace CLHEP; 56 56 57 ////////////////////////////////////////////// 57 ////////////////////////////////////////////////////////////////////////// 58 // 58 // 59 // Constructor 59 // Constructor 60 60 61 G4EllipticalTube::G4EllipticalTube( const G4St 61 G4EllipticalTube::G4EllipticalTube( const G4String &name, 62 G4do 62 G4double Dx, 63 G4do 63 G4double Dy, 64 G4do 64 G4double Dz ) 65 : G4VSolid(name), fDx(Dx), fDy(Dy), fDz(Dz) 65 : G4VSolid(name), fDx(Dx), fDy(Dy), fDz(Dz) 66 { 66 { 67 CheckParameters(); 67 CheckParameters(); 68 } 68 } 69 69 70 ////////////////////////////////////////////// 70 ////////////////////////////////////////////////////////////////////////// 71 // 71 // 72 // Fake default constructor - sets only member 72 // Fake default constructor - sets only member data and allocates memory 73 // for usage restri 73 // for usage restricted to object persistency. 74 74 75 G4EllipticalTube::G4EllipticalTube( __void__& 75 G4EllipticalTube::G4EllipticalTube( __void__& a ) 76 : G4VSolid(a), halfTolerance(0.), fDx(0.), f 76 : G4VSolid(a), halfTolerance(0.), fDx(0.), fDy(0.), fDz(0.), 77 fRsph(0.), fDDx(0.), fDDy(0.), fSx(0.), fS 77 fRsph(0.), fDDx(0.), fDDy(0.), fSx(0.), fSy(0.), fR(0.), 78 fQ1(0.), fQ2(0.), fScratch(0.) 78 fQ1(0.), fQ2(0.), fScratch(0.) 79 { 79 { 80 } 80 } 81 81 82 ////////////////////////////////////////////// 82 ////////////////////////////////////////////////////////////////////////// 83 // 83 // 84 // Destructor 84 // Destructor 85 85 86 G4EllipticalTube::~G4EllipticalTube() 86 G4EllipticalTube::~G4EllipticalTube() 87 { 87 { 88 delete fpPolyhedron; fpPolyhedron = nullptr; 88 delete fpPolyhedron; fpPolyhedron = nullptr; 89 } 89 } 90 90 91 ////////////////////////////////////////////// 91 ////////////////////////////////////////////////////////////////////////// 92 // 92 // 93 // Copy constructor 93 // Copy constructor 94 94 95 G4EllipticalTube::G4EllipticalTube(const G4Ell 95 G4EllipticalTube::G4EllipticalTube(const G4EllipticalTube& rhs) 96 : G4VSolid(rhs), halfTolerance(rhs.halfToler 96 : G4VSolid(rhs), halfTolerance(rhs.halfTolerance), 97 fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), 97 fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), 98 fCubicVolume(rhs.fCubicVolume), fSurfaceAr 98 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea), 99 fRsph(rhs.fRsph), fDDx(rhs.fDDx), fDDy(rhs 99 fRsph(rhs.fRsph), fDDx(rhs.fDDx), fDDy(rhs.fDDy), 100 fSx(rhs.fSx), fSy(rhs.fSy), fR(rhs.fR), 100 fSx(rhs.fSx), fSy(rhs.fSy), fR(rhs.fR), 101 fQ1(rhs.fQ1), fQ2(rhs.fQ2), fScratch(rhs.f 101 fQ1(rhs.fQ1), fQ2(rhs.fQ2), fScratch(rhs.fScratch) 102 { 102 { 103 } 103 } 104 104 105 ////////////////////////////////////////////// 105 ////////////////////////////////////////////////////////////////////////// 106 // 106 // 107 // Assignment operator 107 // Assignment operator 108 108 109 G4EllipticalTube& G4EllipticalTube::operator = 109 G4EllipticalTube& G4EllipticalTube::operator = (const G4EllipticalTube& rhs) 110 { 110 { 111 // Check assignment to self 111 // Check assignment to self 112 // 112 // 113 if (this == &rhs) { return *this; } 113 if (this == &rhs) { return *this; } 114 114 115 // Copy base class data 115 // Copy base class data 116 // 116 // 117 G4VSolid::operator=(rhs); 117 G4VSolid::operator=(rhs); 118 118 119 // Copy data 119 // Copy data 120 // 120 // 121 halfTolerance = rhs.halfTolerance; 121 halfTolerance = rhs.halfTolerance; 122 fDx = rhs.fDx; 122 fDx = rhs.fDx; 123 fDy = rhs.fDy; 123 fDy = rhs.fDy; 124 fDz = rhs.fDz; 124 fDz = rhs.fDz; 125 fCubicVolume = rhs.fCubicVolume; 125 fCubicVolume = rhs.fCubicVolume; 126 fSurfaceArea = rhs.fSurfaceArea; 126 fSurfaceArea = rhs.fSurfaceArea; 127 127 128 fRsph = rhs.fRsph; 128 fRsph = rhs.fRsph; 129 fDDx = rhs.fDDx; 129 fDDx = rhs.fDDx; 130 fDDy = rhs.fDDy; 130 fDDy = rhs.fDDy; 131 fSx = rhs.fSx; 131 fSx = rhs.fSx; 132 fSy = rhs.fSy; 132 fSy = rhs.fSy; 133 fR = rhs.fR; 133 fR = rhs.fR; 134 fQ1 = rhs.fQ1; 134 fQ1 = rhs.fQ1; 135 fQ2 = rhs.fQ2; 135 fQ2 = rhs.fQ2; 136 fScratch = rhs.fScratch; 136 fScratch = rhs.fScratch; 137 137 138 fRebuildPolyhedron = false; 138 fRebuildPolyhedron = false; 139 delete fpPolyhedron; fpPolyhedron = nullptr 139 delete fpPolyhedron; fpPolyhedron = nullptr; 140 140 141 return *this; 141 return *this; 142 } 142 } 143 143 144 ////////////////////////////////////////////// 144 ////////////////////////////////////////////////////////////////////////// 145 // 145 // 146 // Check dimensions 146 // Check dimensions 147 147 148 void G4EllipticalTube::CheckParameters() 148 void G4EllipticalTube::CheckParameters() 149 { 149 { 150 // Check dimensions 150 // Check dimensions 151 // 151 // 152 halfTolerance = 0.5*kCarTolerance; // half t 152 halfTolerance = 0.5*kCarTolerance; // half tolerance 153 G4double dmin = 2*kCarTolerance; 153 G4double dmin = 2*kCarTolerance; 154 if (fDx < dmin || fDy < dmin || fDz < dmin) 154 if (fDx < dmin || fDy < dmin || fDz < dmin) 155 { 155 { 156 std::ostringstream message; 156 std::ostringstream message; 157 message << "Invalid (too small or negative 157 message << "Invalid (too small or negative) dimensions for Solid: " 158 << GetName() 158 << GetName() 159 << "\n Dx = " << fDx 159 << "\n Dx = " << fDx 160 << "\n Dy = " << fDy 160 << "\n Dy = " << fDy 161 << "\n Dz = " << fDz; 161 << "\n Dz = " << fDz; 162 G4Exception("G4EllipticalTube::CheckParame 162 G4Exception("G4EllipticalTube::CheckParameters()", "GeomSolids0002", 163 FatalException, message); 163 FatalException, message); 164 } 164 } 165 165 166 // Set pre-calculatated values 166 // Set pre-calculatated values 167 // 167 // 168 halfTolerance = 0.5*kCarTolerance; // half t 168 halfTolerance = 0.5*kCarTolerance; // half tolerance 169 fRsph = std::sqrt(fDx * fDx + fDy * fDy + fD 169 fRsph = std::sqrt(fDx * fDx + fDy * fDy + fDz * fDz); // radius of surrounding sphere 170 fDDx = fDx * fDx; // X semi-axis squared 170 fDDx = fDx * fDx; // X semi-axis squared 171 fDDy = fDy * fDy; // Y semi-axis squared 171 fDDy = fDy * fDy; // Y semi-axis squared 172 172 173 fR = std::min(fDx, fDy); // resulting radius 173 fR = std::min(fDx, fDy); // resulting radius, after scaling elipse to circle 174 fSx = fR / fDx; // X scale factor 174 fSx = fR / fDx; // X scale factor 175 fSy = fR / fDy; // Y scale factor 175 fSy = fR / fDy; // Y scale factor 176 176 177 fQ1 = 0.5 / fR; // distance approxiamtion di 177 fQ1 = 0.5 / fR; // distance approxiamtion dist = Q1 * (x^2 + y^2) - Q2 178 fQ2 = 0.5 * (fR + halfTolerance * halfTolera 178 fQ2 = 0.5 * (fR + halfTolerance * halfTolerance / fR); 179 fScratch = 2. * fR * fR * DBL_EPSILON; // sc 179 fScratch = 2. * fR * fR * DBL_EPSILON; // scratch within calculation error thickness 180 // fScratch = (B * B / A) * (2. + halfTolera 180 // fScratch = (B * B / A) * (2. + halfTolerance / A) * halfTolerance; // alternative 181 } 181 } 182 182 183 ////////////////////////////////////////////// 183 ////////////////////////////////////////////////////////////////////////// 184 // 184 // 185 // Get bounding box 185 // Get bounding box 186 186 187 void G4EllipticalTube::BoundingLimits( G4Three 187 void G4EllipticalTube::BoundingLimits( G4ThreeVector& pMin, 188 G4Three 188 G4ThreeVector& pMax ) const 189 { 189 { 190 pMin.set(-fDx,-fDy,-fDz); 190 pMin.set(-fDx,-fDy,-fDz); 191 pMax.set( fDx, fDy, fDz); 191 pMax.set( fDx, fDy, fDz); 192 } 192 } 193 193 194 ////////////////////////////////////////////// 194 ////////////////////////////////////////////////////////////////////////// 195 // 195 // 196 // Calculate extent under transform and specif 196 // Calculate extent under transform and specified limit 197 197 198 G4bool 198 G4bool 199 G4EllipticalTube::CalculateExtent( const EAxis 199 G4EllipticalTube::CalculateExtent( const EAxis pAxis, 200 const G4Vox 200 const G4VoxelLimits& pVoxelLimit, 201 const G4Aff 201 const G4AffineTransform& pTransform, 202 G4dou 202 G4double& pMin, G4double& pMax ) const 203 { 203 { 204 G4ThreeVector bmin, bmax; 204 G4ThreeVector bmin, bmax; 205 G4bool exist; 205 G4bool exist; 206 206 207 // Check bounding box (bbox) 207 // Check bounding box (bbox) 208 // 208 // 209 BoundingLimits(bmin,bmax); 209 BoundingLimits(bmin,bmax); 210 G4BoundingEnvelope bbox(bmin,bmax); 210 G4BoundingEnvelope bbox(bmin,bmax); 211 #ifdef G4BBOX_EXTENT 211 #ifdef G4BBOX_EXTENT 212 return bbox.CalculateExtent(pAxis,pVoxelLimi 212 return bbox.CalculateExtent(pAxis,pVoxelLimit, pTransform, pMin, pMax); 213 #endif 213 #endif 214 if (bbox.BoundingBoxVsVoxelLimits(pAxis, pVo 214 if (bbox.BoundingBoxVsVoxelLimits(pAxis, pVoxelLimit, pTransform, pMin, pMax)) 215 { 215 { 216 return exist = pMin < pMax; << 216 return exist = (pMin < pMax) ? true : false; 217 } 217 } 218 218 219 G4double dx = fDx; 219 G4double dx = fDx; 220 G4double dy = fDy; 220 G4double dy = fDy; 221 G4double dz = fDz; 221 G4double dz = fDz; 222 222 223 // Set bounding envelope (benv) and calculat 223 // Set bounding envelope (benv) and calculate extent 224 // 224 // 225 const G4int NSTEPS = 24; // number of steps 225 const G4int NSTEPS = 24; // number of steps for whole circle 226 G4double ang = twopi/NSTEPS; 226 G4double ang = twopi/NSTEPS; 227 227 228 G4double sinHalf = std::sin(0.5*ang); 228 G4double sinHalf = std::sin(0.5*ang); 229 G4double cosHalf = std::cos(0.5*ang); 229 G4double cosHalf = std::cos(0.5*ang); 230 G4double sinStep = 2.*sinHalf*cosHalf; 230 G4double sinStep = 2.*sinHalf*cosHalf; 231 G4double cosStep = 1. - 2.*sinHalf*sinHalf; 231 G4double cosStep = 1. - 2.*sinHalf*sinHalf; 232 G4double sx = dx/cosHalf; 232 G4double sx = dx/cosHalf; 233 G4double sy = dy/cosHalf; 233 G4double sy = dy/cosHalf; 234 234 235 G4double sinCur = sinHalf; 235 G4double sinCur = sinHalf; 236 G4double cosCur = cosHalf; 236 G4double cosCur = cosHalf; 237 G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS 237 G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS); 238 for (G4int k=0; k<NSTEPS; ++k) 238 for (G4int k=0; k<NSTEPS; ++k) 239 { 239 { 240 baseA[k].set(sx*cosCur,sy*sinCur,-dz); 240 baseA[k].set(sx*cosCur,sy*sinCur,-dz); 241 baseB[k].set(sx*cosCur,sy*sinCur, dz); 241 baseB[k].set(sx*cosCur,sy*sinCur, dz); 242 242 243 G4double sinTmp = sinCur; 243 G4double sinTmp = sinCur; 244 sinCur = sinCur*cosStep + cosCur*sinStep; 244 sinCur = sinCur*cosStep + cosCur*sinStep; 245 cosCur = cosCur*cosStep - sinTmp*sinStep; 245 cosCur = cosCur*cosStep - sinTmp*sinStep; 246 } 246 } 247 247 248 std::vector<const G4ThreeVectorList *> polyg 248 std::vector<const G4ThreeVectorList *> polygons(2); 249 polygons[0] = &baseA; 249 polygons[0] = &baseA; 250 polygons[1] = &baseB; 250 polygons[1] = &baseB; 251 G4BoundingEnvelope benv(bmin, bmax, polygons 251 G4BoundingEnvelope benv(bmin, bmax, polygons); 252 exist = benv.CalculateExtent(pAxis, pVoxelLi 252 exist = benv.CalculateExtent(pAxis, pVoxelLimit, pTransform, pMin, pMax); 253 return exist; 253 return exist; 254 } 254 } 255 255 256 ////////////////////////////////////////////// 256 ////////////////////////////////////////////////////////////////////////// 257 // 257 // 258 // Determine where is point: inside, outside o 258 // Determine where is point: inside, outside or on surface 259 // 259 // 260 260 261 EInside G4EllipticalTube::Inside( const G4Thre 261 EInside G4EllipticalTube::Inside( const G4ThreeVector& p ) const 262 { 262 { 263 G4double x = p.x() * fSx; 263 G4double x = p.x() * fSx; 264 G4double y = p.y() * fSy; 264 G4double y = p.y() * fSy; 265 G4double distR = fQ1 * (x * x + y * y) - fQ2 265 G4double distR = fQ1 * (x * x + y * y) - fQ2; 266 G4double distZ = std::abs(p.z()) - fDz; 266 G4double distZ = std::abs(p.z()) - fDz; 267 G4double dist = std::max(distR, distZ); 267 G4double dist = std::max(distR, distZ); 268 268 269 if (dist > halfTolerance) return kOutside; 269 if (dist > halfTolerance) return kOutside; 270 return (dist > -halfTolerance) ? kSurface : 270 return (dist > -halfTolerance) ? kSurface : kInside; 271 } 271 } 272 272 273 ////////////////////////////////////////////// 273 ////////////////////////////////////////////////////////////////////////// 274 // 274 // 275 // Return unit normal at surface closest to p 275 // Return unit normal at surface closest to p 276 276 277 G4ThreeVector G4EllipticalTube::SurfaceNormal( 277 G4ThreeVector G4EllipticalTube::SurfaceNormal( const G4ThreeVector& p ) const 278 { 278 { 279 G4ThreeVector norm(0, 0, 0); 279 G4ThreeVector norm(0, 0, 0); 280 G4int nsurf = 0; 280 G4int nsurf = 0; 281 281 282 // check lateral surface 282 // check lateral surface 283 G4double x = p.x() * fSx; 283 G4double x = p.x() * fSx; 284 G4double y = p.y() * fSy; 284 G4double y = p.y() * fSy; 285 G4double distR = fQ1 * (x * x + y * y) - fQ2 285 G4double distR = fQ1 * (x * x + y * y) - fQ2; 286 if (std::abs(distR) <= halfTolerance) 286 if (std::abs(distR) <= halfTolerance) 287 { 287 { 288 norm = G4ThreeVector(p.x() * fDDy, p.y() * 288 norm = G4ThreeVector(p.x() * fDDy, p.y() * fDDx, 0.).unit(); 289 ++nsurf; 289 ++nsurf; 290 } 290 } 291 291 292 // check lateral bases 292 // check lateral bases 293 G4double distZ = std::abs(p.z()) - fDz; 293 G4double distZ = std::abs(p.z()) - fDz; 294 if (std::abs(distZ) <= halfTolerance) 294 if (std::abs(distZ) <= halfTolerance) 295 { 295 { 296 norm.setZ(p.z() < 0 ? -1. : 1.); 296 norm.setZ(p.z() < 0 ? -1. : 1.); 297 ++nsurf; 297 ++nsurf; 298 } 298 } 299 299 300 // return normal 300 // return normal 301 if (nsurf == 1) return norm; 301 if (nsurf == 1) return norm; 302 else if (nsurf > 1) return norm.unit(); // e 302 else if (nsurf > 1) return norm.unit(); // edge 303 else 303 else 304 { 304 { 305 // Point is not on the surface 305 // Point is not on the surface 306 // 306 // 307 #ifdef G4SPECDEBUG 307 #ifdef G4SPECDEBUG 308 std::ostringstream message; 308 std::ostringstream message; 309 G4long oldprc = message.precision(16); << 309 G4int oldprc = message.precision(16); 310 message << "Point p is not on surface (!?) 310 message << "Point p is not on surface (!?) of solid: " 311 << GetName() << G4endl; 311 << GetName() << G4endl; 312 message << "Position:\n"; 312 message << "Position:\n"; 313 message << " p.x() = " << p.x()/mm << " 313 message << " p.x() = " << p.x()/mm << " mm\n"; 314 message << " p.y() = " << p.y()/mm << " 314 message << " p.y() = " << p.y()/mm << " mm\n"; 315 message << " p.z() = " << p.z()/mm << " 315 message << " p.z() = " << p.z()/mm << " mm"; 316 G4cout.precision(oldprc); 316 G4cout.precision(oldprc); 317 G4Exception("G4EllipticalTube::SurfaceNorm 317 G4Exception("G4EllipticalTube::SurfaceNormal(p)", "GeomSolids1002", 318 JustWarning, message ); 318 JustWarning, message ); 319 DumpInfo(); 319 DumpInfo(); 320 #endif 320 #endif 321 return ApproxSurfaceNormal(p); 321 return ApproxSurfaceNormal(p); 322 } 322 } 323 } 323 } 324 324 325 ////////////////////////////////////////////// 325 ////////////////////////////////////////////////////////////////////////// 326 // 326 // 327 // Find surface nearest to point and return co 327 // Find surface nearest to point and return corresponding normal. 328 // The algorithm is similar to the algorithm u 328 // The algorithm is similar to the algorithm used in Inside(). 329 // This method normally should not be called. 329 // This method normally should not be called. 330 330 331 G4ThreeVector 331 G4ThreeVector 332 G4EllipticalTube::ApproxSurfaceNormal( const G 332 G4EllipticalTube::ApproxSurfaceNormal( const G4ThreeVector& p ) const 333 { 333 { 334 G4double x = p.x() * fSx; 334 G4double x = p.x() * fSx; 335 G4double y = p.y() * fSy; 335 G4double y = p.y() * fSy; 336 G4double distR = fQ1 * (x * x + y * y) - fQ2 336 G4double distR = fQ1 * (x * x + y * y) - fQ2; 337 G4double distZ = std::abs(p.z()) - fDz; 337 G4double distZ = std::abs(p.z()) - fDz; 338 if (distR > distZ && (x * x + y * y) > 0) 338 if (distR > distZ && (x * x + y * y) > 0) 339 return G4ThreeVector(p.x() * fDDy, p.y() * 339 return G4ThreeVector(p.x() * fDDy, p.y() * fDDx, 0.).unit(); 340 else 340 else 341 return {0, 0, (p.z() < 0 ? -1. : 1.)}; << 341 return G4ThreeVector(0, 0, (p.z() < 0 ? -1. : 1.)); 342 } 342 } 343 343 344 ////////////////////////////////////////////// 344 ////////////////////////////////////////////////////////////////////////// 345 // 345 // 346 // Calculate distance to shape from outside, a 346 // Calculate distance to shape from outside, along normalised vector, 347 // return kInfinity if no intersection, or dis 347 // return kInfinity if no intersection, or distance < halfTolerance 348 348 349 G4double G4EllipticalTube::DistanceToIn( const 349 G4double G4EllipticalTube::DistanceToIn( const G4ThreeVector& p, 350 const 350 const G4ThreeVector& v ) const 351 { 351 { 352 G4double offset = 0.; 352 G4double offset = 0.; 353 G4ThreeVector pcur = p; 353 G4ThreeVector pcur = p; 354 354 355 // Check if point is flying away 355 // Check if point is flying away 356 // 356 // 357 G4double safex = std::abs(pcur.x()) - fDx; 357 G4double safex = std::abs(pcur.x()) - fDx; 358 G4double safey = std::abs(pcur.y()) - fDy; 358 G4double safey = std::abs(pcur.y()) - fDy; 359 G4double safez = std::abs(pcur.z()) - fDz; 359 G4double safez = std::abs(pcur.z()) - fDz; 360 360 361 if (safez >= -halfTolerance && pcur.z() * v. 361 if (safez >= -halfTolerance && pcur.z() * v.z() >= 0.) return kInfinity; 362 if (safey >= -halfTolerance && pcur.y() * v. 362 if (safey >= -halfTolerance && pcur.y() * v.y() >= 0.) return kInfinity; 363 if (safex >= -halfTolerance && pcur.x() * v. 363 if (safex >= -halfTolerance && pcur.x() * v.x() >= 0.) return kInfinity; 364 364 365 // Relocate point, if required 365 // Relocate point, if required 366 // 366 // 367 G4double Dmax = 32. * fRsph; 367 G4double Dmax = 32. * fRsph; 368 if (std::max(std::max(safex, safey), safez) 368 if (std::max(std::max(safex, safey), safez) > Dmax) 369 { 369 { 370 offset = (1. - 1.e-08) * pcur.mag() - 2. * 370 offset = (1. - 1.e-08) * pcur.mag() - 2. * fRsph; 371 pcur += offset * v; 371 pcur += offset * v; 372 G4double dist = DistanceToIn(pcur, v); 372 G4double dist = DistanceToIn(pcur, v); 373 return (dist == kInfinity) ? kInfinity : d 373 return (dist == kInfinity) ? kInfinity : dist + offset; 374 } 374 } 375 375 376 // Scale elliptical tube to cylinder 376 // Scale elliptical tube to cylinder 377 // 377 // 378 G4double px = pcur.x() * fSx; 378 G4double px = pcur.x() * fSx; 379 G4double py = pcur.y() * fSy; 379 G4double py = pcur.y() * fSy; 380 G4double pz = pcur.z(); 380 G4double pz = pcur.z(); 381 G4double vx = v.x() * fSx; 381 G4double vx = v.x() * fSx; 382 G4double vy = v.y() * fSy; 382 G4double vy = v.y() * fSy; 383 G4double vz = v.z(); 383 G4double vz = v.z(); 384 384 385 // Set coefficients of quadratic equation: A 385 // Set coefficients of quadratic equation: A t^2 + 2B t + C = 0 386 // 386 // 387 G4double rr = px * px + py * py; 387 G4double rr = px * px + py * py; 388 G4double A = vx * vx + vy * vy; 388 G4double A = vx * vx + vy * vy; 389 G4double B = px * vx + py * vy; 389 G4double B = px * vx + py * vy; 390 G4double C = rr - fR * fR; 390 G4double C = rr - fR * fR; 391 G4double D = B * B - A * C; 391 G4double D = B * B - A * C; 392 392 393 // Check if point is flying away relative to 393 // Check if point is flying away relative to lateral surface 394 // 394 // 395 G4double distR = fQ1 * rr - fQ2; 395 G4double distR = fQ1 * rr - fQ2; 396 G4bool parallelToZ = (A < DBL_EPSILON || std 396 G4bool parallelToZ = (A < DBL_EPSILON || std::abs(vz) >= 1.); 397 if (distR >= -halfTolerance && (B >= 0. || p 397 if (distR >= -halfTolerance && (B >= 0. || parallelToZ)) return kInfinity; 398 398 399 // Find intersection with Z planes 399 // Find intersection with Z planes 400 // 400 // 401 G4double invz = (vz == 0) ? DBL_MAX : -1./v 401 G4double invz = (vz == 0) ? DBL_MAX : -1./vz; 402 G4double dz = std::copysign(fDz, invz); 402 G4double dz = std::copysign(fDz, invz); 403 G4double tzmin = (pz - dz) * invz; 403 G4double tzmin = (pz - dz) * invz; 404 G4double tzmax = (pz + dz) * invz; 404 G4double tzmax = (pz + dz) * invz; 405 405 406 // Solve qudratic equation. There are two ca 406 // Solve qudratic equation. There are two cases special where D <= 0: 407 // 1) trajectory parallel to Z axis (A = 0 407 // 1) trajectory parallel to Z axis (A = 0, B = 0, C - any, D = 0) 408 // 2) touch (D = 0) or no intersection (D 408 // 2) touch (D = 0) or no intersection (D < 0) with lateral surface 409 // 409 // 410 if (parallelToZ) return (tzmin<halfTolerance 410 if (parallelToZ) return (tzmin<halfTolerance) ? offset : tzmin + offset; // 1) 411 if (D <= A * A * fScratch) return kInfinity; 411 if (D <= A * A * fScratch) return kInfinity; // 2) 412 412 413 // Find roots of quadratic equation 413 // Find roots of quadratic equation 414 G4double tmp = -B - std::copysign(std::sqrt( 414 G4double tmp = -B - std::copysign(std::sqrt(D), B); 415 G4double t1 = tmp / A; 415 G4double t1 = tmp / A; 416 G4double t2 = C / tmp; 416 G4double t2 = C / tmp; 417 G4double trmin = std::min(t1, t2); 417 G4double trmin = std::min(t1, t2); 418 G4double trmax = std::max(t1, t2); 418 G4double trmax = std::max(t1, t2); 419 419 420 // Return distance 420 // Return distance 421 G4double tin = std::max(tzmin, trmin); 421 G4double tin = std::max(tzmin, trmin); 422 G4double tout = std::min(tzmax, trmax); 422 G4double tout = std::min(tzmax, trmax); 423 423 424 if (tout <= tin + halfTolerance) return kInf 424 if (tout <= tin + halfTolerance) return kInfinity; // touch or no hit 425 return (tin<halfTolerance) ? offset : tin + 425 return (tin<halfTolerance) ? offset : tin + offset; 426 } 426 } 427 427 428 ////////////////////////////////////////////// 428 ////////////////////////////////////////////////////////////////////////// 429 // 429 // 430 // Estimate distance to the surface from outsi 430 // Estimate distance to the surface from outside, 431 // returns 0 if point is inside 431 // returns 0 if point is inside 432 432 433 G4double G4EllipticalTube::DistanceToIn( const 433 G4double G4EllipticalTube::DistanceToIn( const G4ThreeVector& p ) const 434 { 434 { 435 // safety distance to bounding box 435 // safety distance to bounding box 436 G4double distX = std::abs(p.x()) - fDx; 436 G4double distX = std::abs(p.x()) - fDx; 437 G4double distY = std::abs(p.y()) - fDy; 437 G4double distY = std::abs(p.y()) - fDy; 438 G4double distZ = std::abs(p.z()) - fDz; 438 G4double distZ = std::abs(p.z()) - fDz; 439 G4double distB = std::max(std::max(distX, di 439 G4double distB = std::max(std::max(distX, distY), distZ); 440 // return (distB < 0) ? 0 : distB; 440 // return (distB < 0) ? 0 : distB; 441 441 442 // safety distance to lateral surface 442 // safety distance to lateral surface 443 G4double x = p.x() * fSx; 443 G4double x = p.x() * fSx; 444 G4double y = p.y() * fSy; 444 G4double y = p.y() * fSy; 445 G4double distR = std::sqrt(x * x + y * y) - 445 G4double distR = std::sqrt(x * x + y * y) - fR; 446 446 447 // return SafetyToIn 447 // return SafetyToIn 448 G4double dist = std::max(distB, distR); 448 G4double dist = std::max(distB, distR); 449 return (dist < 0) ? 0 : dist; 449 return (dist < 0) ? 0 : dist; 450 } 450 } 451 451 452 ////////////////////////////////////////////// 452 ////////////////////////////////////////////////////////////////////////// 453 // 453 // 454 // Calculate distance to shape from inside and 454 // Calculate distance to shape from inside and find normal 455 // at exit point, if required 455 // at exit point, if required 456 // - when leaving the surface, return 0 456 // - when leaving the surface, return 0 457 457 458 G4double G4EllipticalTube::DistanceToOut( cons 458 G4double G4EllipticalTube::DistanceToOut( const G4ThreeVector& p, 459 cons 459 const G4ThreeVector& v, 460 cons 460 const G4bool calcNorm, 461 461 G4bool* validNorm, 462 462 G4ThreeVector* n ) const 463 { 463 { 464 // Check if point flying away relative to Z 464 // Check if point flying away relative to Z planes 465 // 465 // 466 G4double pz = p.z(); 466 G4double pz = p.z(); 467 G4double vz = v.z(); 467 G4double vz = v.z(); 468 G4double distZ = std::abs(pz) - fDz; 468 G4double distZ = std::abs(pz) - fDz; 469 if (distZ >= -halfTolerance && pz * vz > 0) 469 if (distZ >= -halfTolerance && pz * vz > 0) 470 { 470 { 471 if (calcNorm) 471 if (calcNorm) 472 { 472 { 473 *validNorm = true; 473 *validNorm = true; 474 n->set(0, 0, (pz < 0) ? -1. : 1.); 474 n->set(0, 0, (pz < 0) ? -1. : 1.); 475 } 475 } 476 return 0.; 476 return 0.; 477 } 477 } 478 G4double tzmax = (vz == 0) ? DBL_MAX : (std: 478 G4double tzmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz, vz) - pz) / vz; 479 479 480 // Scale elliptical tube to cylinder 480 // Scale elliptical tube to cylinder 481 // 481 // 482 G4double px = p.x() * fSx; 482 G4double px = p.x() * fSx; 483 G4double py = p.y() * fSy; 483 G4double py = p.y() * fSy; 484 G4double vx = v.x() * fSx; 484 G4double vx = v.x() * fSx; 485 G4double vy = v.y() * fSy; 485 G4double vy = v.y() * fSy; 486 486 487 // Check if point is flying away relative to 487 // Check if point is flying away relative to lateral surface 488 // 488 // 489 G4double rr = px * px + py * py; 489 G4double rr = px * px + py * py; 490 G4double B = px * vx + py * vy; 490 G4double B = px * vx + py * vy; 491 G4double distR = fQ1 * rr - fQ2; 491 G4double distR = fQ1 * rr - fQ2; 492 if (distR >= -halfTolerance && B > 0.) 492 if (distR >= -halfTolerance && B > 0.) 493 { 493 { 494 if (calcNorm) 494 if (calcNorm) 495 { 495 { 496 *validNorm = true; 496 *validNorm = true; 497 *n = G4ThreeVector(px * fDDy, py * fDDx, 497 *n = G4ThreeVector(px * fDDy, py * fDDx, 0.).unit(); 498 } 498 } 499 return 0.; 499 return 0.; 500 } 500 } 501 501 502 // Just in case check if point is outside, n 502 // Just in case check if point is outside, normally it should never be 503 // 503 // 504 if (std::max(distZ, distR) > halfTolerance) 504 if (std::max(distZ, distR) > halfTolerance) 505 { 505 { 506 #ifdef G4SPECDEBUG 506 #ifdef G4SPECDEBUG 507 std::ostringstream message; 507 std::ostringstream message; 508 G4long oldprc = message.precision(16); << 508 G4int oldprc = message.precision(16); 509 message << "Point p is outside (!?) of sol 509 message << "Point p is outside (!?) of solid: " 510 << GetName() << G4endl; 510 << GetName() << G4endl; 511 message << "Position: " << p << G4endl;; 511 message << "Position: " << p << G4endl;; 512 message << "Direction: " << v; 512 message << "Direction: " << v; 513 G4cout.precision(oldprc); 513 G4cout.precision(oldprc); 514 G4Exception("G4EllipticalTube::DistanceToO 514 G4Exception("G4EllipticalTube::DistanceToOut(p,v)", "GeomSolids1002", 515 JustWarning, message ); 515 JustWarning, message ); 516 DumpInfo(); 516 DumpInfo(); 517 #endif 517 #endif 518 if (calcNorm) 518 if (calcNorm) 519 { 519 { 520 *validNorm = true; 520 *validNorm = true; 521 *n = ApproxSurfaceNormal(p); 521 *n = ApproxSurfaceNormal(p); 522 } 522 } 523 return 0.; 523 return 0.; 524 } 524 } 525 525 526 // Set coefficients of quadratic equation: A 526 // Set coefficients of quadratic equation: A t^2 + 2B t + C = 0 527 // 527 // 528 G4double A = vx * vx + vy * vy; 528 G4double A = vx * vx + vy * vy; 529 G4double C = rr - fR * fR; 529 G4double C = rr - fR * fR; 530 G4double D = B * B - A * C; 530 G4double D = B * B - A * C; 531 531 532 // Solve qudratic equation. There are two sp 532 // Solve qudratic equation. There are two special cases where D <= 0: 533 // 1) trajectory parallel to Z axis (A = 0 533 // 1) trajectory parallel to Z axis (A = 0, B = 0, C - any, D = 0) 534 // 2) touch (D = 0) or no intersection (D 534 // 2) touch (D = 0) or no intersection (D < 0) with lateral surface 535 // 535 // 536 G4bool parallelToZ = (A < DBL_EPSILON || std 536 G4bool parallelToZ = (A < DBL_EPSILON || std::abs(vz) >= 1.); 537 if (parallelToZ) // 1) 537 if (parallelToZ) // 1) 538 { 538 { 539 if (calcNorm) 539 if (calcNorm) 540 { 540 { 541 *validNorm = true; 541 *validNorm = true; 542 n->set(0, 0, (vz < 0) ? -1. : 1.); 542 n->set(0, 0, (vz < 0) ? -1. : 1.); 543 } 543 } 544 return tzmax; 544 return tzmax; 545 } 545 } 546 if (D <= A * A * fScratch) // 2) 546 if (D <= A * A * fScratch) // 2) 547 { 547 { 548 if (calcNorm) 548 if (calcNorm) 549 { 549 { 550 *validNorm = true; 550 *validNorm = true; 551 *n = G4ThreeVector(px * fDDy, py * fDDx, 551 *n = G4ThreeVector(px * fDDy, py * fDDx, 0.).unit(); 552 } 552 } 553 return 0.; 553 return 0.; 554 } 554 } 555 555 556 // Find roots of quadratic equation 556 // Find roots of quadratic equation 557 G4double tmp = -B - std::copysign(std::sqrt( 557 G4double tmp = -B - std::copysign(std::sqrt(D), B); 558 G4double t1 = tmp / A; 558 G4double t1 = tmp / A; 559 G4double t2 = C / tmp; 559 G4double t2 = C / tmp; 560 G4double trmax = std::max(t1, t2); 560 G4double trmax = std::max(t1, t2); 561 561 562 // Return distance 562 // Return distance 563 G4double tmax = std::min(tzmax, trmax); 563 G4double tmax = std::min(tzmax, trmax); 564 564 565 // Set normal, if required, and return dista 565 // Set normal, if required, and return distance 566 // 566 // 567 if (calcNorm) 567 if (calcNorm) 568 { 568 { 569 *validNorm = true; 569 *validNorm = true; 570 G4ThreeVector pnew = p + tmax * v; 570 G4ThreeVector pnew = p + tmax * v; 571 if (tmax == tzmax) 571 if (tmax == tzmax) 572 n->set(0, 0, (pnew.z() < 0) ? -1. : 1.); 572 n->set(0, 0, (pnew.z() < 0) ? -1. : 1.); 573 else 573 else 574 *n = G4ThreeVector(pnew.x() * fDDy, pnew 574 *n = G4ThreeVector(pnew.x() * fDDy, pnew.y() * fDDx, 0.).unit(); 575 } 575 } 576 return tmax; 576 return tmax; 577 } 577 } 578 578 579 ////////////////////////////////////////////// 579 ////////////////////////////////////////////////////////////////////////// 580 // 580 // 581 // Estimate distance to the surface from insid 581 // Estimate distance to the surface from inside, 582 // returns 0 if point is outside 582 // returns 0 if point is outside 583 // 583 // 584 584 585 G4double G4EllipticalTube::DistanceToOut( cons 585 G4double G4EllipticalTube::DistanceToOut( const G4ThreeVector& p ) const 586 { 586 { 587 #ifdef G4SPECDEBUG 587 #ifdef G4SPECDEBUG 588 if( Inside(p) == kOutside ) 588 if( Inside(p) == kOutside ) 589 { 589 { 590 std::ostringstream message; 590 std::ostringstream message; 591 G4long oldprc = message.precision(16); << 591 G4int oldprc = message.precision(16); 592 message << "Point p is outside (!?) of sol 592 message << "Point p is outside (!?) of solid: " << GetName() << "\n" 593 << "Position:\n" 593 << "Position:\n" 594 << " p.x() = " << p.x()/mm << " 594 << " p.x() = " << p.x()/mm << " mm\n" 595 << " p.y() = " << p.y()/mm << " 595 << " p.y() = " << p.y()/mm << " mm\n" 596 << " p.z() = " << p.z()/mm << " 596 << " p.z() = " << p.z()/mm << " mm"; 597 message.precision(oldprc) ; 597 message.precision(oldprc) ; 598 G4Exception("G4ElliptocalTube::DistanceToO 598 G4Exception("G4ElliptocalTube::DistanceToOut(p)", "GeomSolids1002", 599 JustWarning, message); 599 JustWarning, message); 600 DumpInfo(); 600 DumpInfo(); 601 } 601 } 602 #endif 602 #endif 603 // safety distance to Z-bases 603 // safety distance to Z-bases 604 G4double distZ = fDz - std::abs(p.z()); 604 G4double distZ = fDz - std::abs(p.z()); 605 605 606 // safety distance lateral surface 606 // safety distance lateral surface 607 G4double x = p.x() * fSx; 607 G4double x = p.x() * fSx; 608 G4double y = p.y() * fSy; 608 G4double y = p.y() * fSy; 609 G4double distR = fR - std::sqrt(x * x + y * 609 G4double distR = fR - std::sqrt(x * x + y * y); 610 610 611 // return SafetyToOut 611 // return SafetyToOut 612 G4double dist = std::min(distZ, distR); 612 G4double dist = std::min(distZ, distR); 613 return (dist < 0) ? 0 : dist; 613 return (dist < 0) ? 0 : dist; 614 } 614 } 615 615 616 ////////////////////////////////////////////// 616 ////////////////////////////////////////////////////////////////////////// 617 // 617 // 618 // GetEntityType 618 // GetEntityType 619 619 620 G4GeometryType G4EllipticalTube::GetEntityType 620 G4GeometryType G4EllipticalTube::GetEntityType() const 621 { 621 { 622 return {"G4EllipticalTube"}; << 622 return G4String("G4EllipticalTube"); 623 } 623 } 624 624 625 ////////////////////////////////////////////// 625 ////////////////////////////////////////////////////////////////////////// 626 // 626 // 627 // Make a clone of the object 627 // Make a clone of the object 628 628 629 G4VSolid* G4EllipticalTube::Clone() const 629 G4VSolid* G4EllipticalTube::Clone() const 630 { 630 { 631 return new G4EllipticalTube(*this); 631 return new G4EllipticalTube(*this); 632 } 632 } 633 633 634 ////////////////////////////////////////////// 634 ////////////////////////////////////////////////////////////////////////// 635 // 635 // 636 // Return volume 636 // Return volume 637 637 638 G4double G4EllipticalTube::GetCubicVolume() 638 G4double G4EllipticalTube::GetCubicVolume() 639 { 639 { 640 if (fCubicVolume == 0.) 640 if (fCubicVolume == 0.) 641 { 641 { 642 fCubicVolume = twopi * fDx * fDy * fDz; 642 fCubicVolume = twopi * fDx * fDy * fDz; 643 } 643 } 644 return fCubicVolume; 644 return fCubicVolume; 645 } 645 } 646 646 647 ////////////////////////////////////////////// 647 ////////////////////////////////////////////////////////////////////////// 648 // 648 // 649 // Return cached surface area 649 // Return cached surface area 650 650 651 G4double G4EllipticalTube::GetCachedSurfaceAre 651 G4double G4EllipticalTube::GetCachedSurfaceArea() const 652 { 652 { 653 G4ThreadLocalStatic G4double cached_Dx = 0; 653 G4ThreadLocalStatic G4double cached_Dx = 0; 654 G4ThreadLocalStatic G4double cached_Dy = 0; 654 G4ThreadLocalStatic G4double cached_Dy = 0; 655 G4ThreadLocalStatic G4double cached_Dz = 0; 655 G4ThreadLocalStatic G4double cached_Dz = 0; 656 G4ThreadLocalStatic G4double cached_area = 0 656 G4ThreadLocalStatic G4double cached_area = 0; 657 if (cached_Dx != fDx || cached_Dy != fDy || 657 if (cached_Dx != fDx || cached_Dy != fDy || cached_Dz != fDz) 658 { 658 { 659 cached_Dx = fDx; 659 cached_Dx = fDx; 660 cached_Dy = fDy; 660 cached_Dy = fDy; 661 cached_Dz = fDz; 661 cached_Dz = fDz; 662 cached_area = 2.*(pi*fDx*fDy + G4GeomTools 662 cached_area = 2.*(pi*fDx*fDy + G4GeomTools::EllipsePerimeter(fDx, fDy)*fDz); 663 } 663 } 664 return cached_area; 664 return cached_area; 665 } 665 } 666 666 667 ////////////////////////////////////////////// 667 ////////////////////////////////////////////////////////////////////////// 668 // 668 // 669 // Return surface area 669 // Return surface area 670 670 671 G4double G4EllipticalTube::GetSurfaceArea() 671 G4double G4EllipticalTube::GetSurfaceArea() 672 { 672 { 673 if(fSurfaceArea == 0.) 673 if(fSurfaceArea == 0.) 674 { 674 { 675 fSurfaceArea = GetCachedSurfaceArea(); 675 fSurfaceArea = GetCachedSurfaceArea(); 676 } 676 } 677 return fSurfaceArea; 677 return fSurfaceArea; 678 } 678 } 679 679 680 ////////////////////////////////////////////// 680 ////////////////////////////////////////////////////////////////////////// 681 // 681 // 682 // Stream object contents to output stream 682 // Stream object contents to output stream 683 683 684 std::ostream& G4EllipticalTube::StreamInfo(std 684 std::ostream& G4EllipticalTube::StreamInfo(std::ostream& os) const 685 { 685 { 686 G4long oldprc = os.precision(16); << 686 G4int oldprc = os.precision(16); 687 os << "------------------------------------- 687 os << "-----------------------------------------------------------\n" 688 << " *** Dump for solid - " << GetName 688 << " *** Dump for solid - " << GetName() << " ***\n" 689 << " ================================= 689 << " ===================================================\n" 690 << " Solid type: G4EllipticalTube\n" 690 << " Solid type: G4EllipticalTube\n" 691 << " Parameters: \n" 691 << " Parameters: \n" 692 << " length Z: " << fDz/mm << " mm \n" 692 << " length Z: " << fDz/mm << " mm \n" 693 << " lateral surface equation: \n" 693 << " lateral surface equation: \n" 694 << " (X / " << fDx << ")^2 + (Y / " 694 << " (X / " << fDx << ")^2 + (Y / " << fDy << ")^2 = 1 \n" 695 << "------------------------------------- 695 << "-----------------------------------------------------------\n"; 696 os.precision(oldprc); 696 os.precision(oldprc); 697 697 698 return os; 698 return os; 699 } 699 } 700 700 701 701 702 ////////////////////////////////////////////// 702 ////////////////////////////////////////////////////////////////////////// 703 // 703 // 704 // Pick up a random point on the surface 704 // Pick up a random point on the surface 705 705 706 G4ThreeVector G4EllipticalTube::GetPointOnSurf 706 G4ThreeVector G4EllipticalTube::GetPointOnSurface() const 707 { 707 { 708 // Select surface (0 - base at -Z, 1 - base 708 // Select surface (0 - base at -Z, 1 - base at +Z, 2 - lateral surface) 709 // 709 // 710 G4double sbase = pi * fDx * fDy; 710 G4double sbase = pi * fDx * fDy; 711 G4double ssurf = GetCachedSurfaceArea(); 711 G4double ssurf = GetCachedSurfaceArea(); 712 G4double select = ssurf * G4UniformRand(); 712 G4double select = ssurf * G4UniformRand(); 713 713 714 G4int k = 0; 714 G4int k = 0; 715 if (select > sbase) k = 1; 715 if (select > sbase) k = 1; 716 if (select > 2. * sbase) k = 2; 716 if (select > 2. * sbase) k = 2; 717 717 718 // Pick random point on selected surface (re 718 // Pick random point on selected surface (rejection sampling) 719 // 719 // 720 G4ThreeVector p; 720 G4ThreeVector p; 721 switch (k) { 721 switch (k) { 722 case 0: // base at -Z 722 case 0: // base at -Z 723 { 723 { 724 G4TwoVector rho = G4RandomPointInEllipse 724 G4TwoVector rho = G4RandomPointInEllipse(fDx, fDy); 725 p.set(rho.x(), rho.y(), -fDz); 725 p.set(rho.x(), rho.y(), -fDz); 726 break; 726 break; 727 } 727 } 728 case 1: // base at +Z 728 case 1: // base at +Z 729 { 729 { 730 G4TwoVector rho = G4RandomPointInEllipse 730 G4TwoVector rho = G4RandomPointInEllipse(fDx, fDy); 731 p.set(rho.x(), rho.y(), fDz); 731 p.set(rho.x(), rho.y(), fDz); 732 break; 732 break; 733 } 733 } 734 case 2: // lateral surface 734 case 2: // lateral surface 735 { 735 { 736 G4TwoVector rho = G4RandomPointOnEllipse 736 G4TwoVector rho = G4RandomPointOnEllipse(fDx, fDy); 737 p.set(rho.x(), rho.y(), (2. * G4UniformR 737 p.set(rho.x(), rho.y(), (2. * G4UniformRand() - 1.) * fDz); 738 break; 738 break; 739 } 739 } 740 } 740 } 741 return p; 741 return p; 742 } 742 } 743 743 744 744 745 ////////////////////////////////////////////// 745 ////////////////////////////////////////////////////////////////////////// 746 // 746 // 747 // CreatePolyhedron 747 // CreatePolyhedron 748 748 749 G4Polyhedron* G4EllipticalTube::CreatePolyhedr 749 G4Polyhedron* G4EllipticalTube::CreatePolyhedron() const 750 { 750 { 751 // create cylinder with radius=1... 751 // create cylinder with radius=1... 752 // 752 // 753 G4Polyhedron* eTube = new G4PolyhedronTube(0 753 G4Polyhedron* eTube = new G4PolyhedronTube(0., 1., fDz); 754 754 755 // apply non-uniform scaling... 755 // apply non-uniform scaling... 756 // 756 // 757 eTube->Transform(G4Scale3D(fDx, fDy, 1.)); 757 eTube->Transform(G4Scale3D(fDx, fDy, 1.)); 758 return eTube; 758 return eTube; 759 } 759 } 760 760 761 ////////////////////////////////////////////// 761 ////////////////////////////////////////////////////////////////////////// 762 // 762 // 763 // GetPolyhedron 763 // GetPolyhedron 764 764 765 G4Polyhedron* G4EllipticalTube::GetPolyhedron 765 G4Polyhedron* G4EllipticalTube::GetPolyhedron () const 766 { 766 { 767 if (fpPolyhedron == nullptr || 767 if (fpPolyhedron == nullptr || 768 fRebuildPolyhedron || 768 fRebuildPolyhedron || 769 fpPolyhedron->GetNumberOfRotationStepsAt 769 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() != 770 fpPolyhedron->GetNumberOfRotationSteps() 770 fpPolyhedron->GetNumberOfRotationSteps()) 771 { 771 { 772 G4AutoLock l(&polyhedronMutex); 772 G4AutoLock l(&polyhedronMutex); 773 delete fpPolyhedron; 773 delete fpPolyhedron; 774 fpPolyhedron = CreatePolyhedron(); 774 fpPolyhedron = CreatePolyhedron(); 775 fRebuildPolyhedron = false; 775 fRebuildPolyhedron = false; 776 l.unlock(); 776 l.unlock(); 777 } 777 } 778 return fpPolyhedron; 778 return fpPolyhedron; 779 } 779 } 780 780 781 ////////////////////////////////////////////// 781 ////////////////////////////////////////////////////////////////////////// 782 // 782 // 783 // DescribeYourselfTo 783 // DescribeYourselfTo 784 784 785 void G4EllipticalTube::DescribeYourselfTo( G4V 785 void G4EllipticalTube::DescribeYourselfTo( G4VGraphicsScene& scene ) const 786 { 786 { 787 scene.AddSolid (*this); 787 scene.AddSolid (*this); 788 } 788 } 789 789 790 ////////////////////////////////////////////// 790 ////////////////////////////////////////////////////////////////////////// 791 // 791 // 792 // GetExtent 792 // GetExtent 793 793 794 G4VisExtent G4EllipticalTube::GetExtent() cons 794 G4VisExtent G4EllipticalTube::GetExtent() const 795 { 795 { 796 return { -fDx, fDx, -fDy, fDy, -fDz, fDz }; << 796 return G4VisExtent( -fDx, fDx, -fDy, fDy, -fDz, fDz ); 797 } 797 } 798 798 799 #endif // !defined(G4GEOM_USE_UELLIPTICALTUBE) 799 #endif // !defined(G4GEOM_USE_UELLIPTICALTUBE) || !defined(G4GEOM_USE_SYS_USOLIDS) 800 800