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 // $Id: G4EllipticalCone.cc,v 1.15 2007/08/21 12:58:36 gcosmo Exp $ >> 27 // GEANT4 tag $Name: geant4-09-01 $ >> 28 // 26 // Implementation of G4EllipticalCone class 29 // Implementation of G4EllipticalCone class 27 // 30 // 28 // This code implements an Elliptical Cone giv 31 // This code implements an Elliptical Cone given explicitly by the 29 // equation: 32 // equation: 30 // x^2/a^2 + y^2/b^2 = (z-h)^2 33 // x^2/a^2 + y^2/b^2 = (z-h)^2 31 // and specified by the parameters (a,b,h) and 34 // and specified by the parameters (a,b,h) and a cut parallel to the 32 // xy plane above z = 0. 35 // xy plane above z = 0. 33 // 36 // 34 // Author: Dionysios Anninos 37 // Author: Dionysios Anninos 35 // Revised: Evgueni Tcherniaev << 38 // 36 // ------------------------------------------- 39 // -------------------------------------------------------------------- 37 40 38 #if !(defined(G4GEOM_USE_UELLIPTICALCONE) && d << 39 << 40 #include "globals.hh" 41 #include "globals.hh" 41 42 42 #include "G4EllipticalCone.hh" 43 #include "G4EllipticalCone.hh" 43 44 44 #include "G4RandomTools.hh" << 45 #include "G4GeomTools.hh" << 46 #include "G4ClippablePolygon.hh" 45 #include "G4ClippablePolygon.hh" >> 46 #include "G4SolidExtentList.hh" 47 #include "G4VoxelLimits.hh" 47 #include "G4VoxelLimits.hh" 48 #include "G4AffineTransform.hh" 48 #include "G4AffineTransform.hh" 49 #include "G4BoundingEnvelope.hh" << 50 #include "G4GeometryTolerance.hh" 49 #include "G4GeometryTolerance.hh" 51 50 52 #include "meshdefs.hh" 51 #include "meshdefs.hh" 53 52 54 #include "Randomize.hh" 53 #include "Randomize.hh" 55 54 56 #include "G4VGraphicsScene.hh" 55 #include "G4VGraphicsScene.hh" >> 56 #include "G4Polyhedron.hh" >> 57 #include "G4NURBS.hh" >> 58 #include "G4NURBSbox.hh" 57 #include "G4VisExtent.hh" 59 #include "G4VisExtent.hh" 58 60 59 #include "G4AutoLock.hh" << 61 //#define G4SPECSDEBUG 1 60 << 61 namespace << 62 { << 63 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE << 64 } << 65 62 66 using namespace CLHEP; 63 using namespace CLHEP; 67 64 68 ////////////////////////////////////////////// << 65 ////////////////////////////////////////////////////////////////////// 69 // 66 // 70 // Constructor - check parameters 67 // Constructor - check parameters 71 << 68 // 72 G4EllipticalCone::G4EllipticalCone(const G4Str 69 G4EllipticalCone::G4EllipticalCone(const G4String& pName, 73 G4dou 70 G4double pxSemiAxis, 74 G4dou 71 G4double pySemiAxis, 75 G4dou 72 G4double pzMax, 76 G4dou 73 G4double pzTopCut) 77 : G4VSolid(pName), zTopCut(0.) << 74 : G4VSolid(pName), fpPolyhedron(0), fCubicVolume(0.), fSurfaceArea(0.), >> 75 zTopCut(0.) 78 { 76 { 79 halfCarTol = 0.5*kCarTolerance; << 80 77 81 // Check Semi-Axis & Z-cut << 78 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance(); >> 79 >> 80 // Check Semi-Axis 82 // 81 // 83 if ( (pxSemiAxis <= 0.) || (pySemiAxis <= 0. << 82 if ( (pxSemiAxis > 0.) && (pySemiAxis > 0.) && (pzMax > 0.) ) 84 { 83 { 85 std::ostringstream message; << 84 SetSemiAxis( pxSemiAxis, pySemiAxis, pzMax ); 86 message << "Invalid semi-axis or height f << 85 } 87 << "\n X semi-axis, Y semi-axis << 86 else 88 << pxSemiAxis << ", " << pySemiAx << 89 G4Exception("G4EllipticalCone::G4Elliptic << 90 FatalErrorInArgument, message << 91 } << 92 << 93 if ( pzTopCut <= 0 ) << 94 { 87 { 95 std::ostringstream message; << 88 G4cerr << "ERROR - G4EllipticalCone::G4EllipticalCone(): " 96 message << "Invalid z-coordinate for cutt << 89 << GetName() << G4endl 97 << "\n Z top cut = " << pzTopCu << 90 << " Invalid semi-axis or height!" << G4endl; 98 G4Exception("G4EllipticalCone::G4Elliptic << 91 G4Exception("G4EllipticalCone::G4EllipticalCone()", "InvalidSetup", 99 FatalErrorInArgument, message << 92 FatalException, "Invalid semi-axis or height."); 100 } 93 } 101 94 102 SetSemiAxis( pxSemiAxis, pySemiAxis, pzMax ) << 95 if ( pzTopCut > 0 ) 103 SetZCut(pzTopCut); << 96 { >> 97 SetZCut(pzTopCut); >> 98 } >> 99 else >> 100 { >> 101 G4cerr << "ERROR - G4EllipticalCone::G4EllipticalCone(): " >> 102 << GetName() << G4endl >> 103 << " Invalid z-coordinate for cutting plane !" << G4endl; >> 104 G4Exception("G4EllipticalCone::G4EllipticalCone()", "InvalidSetup", >> 105 FatalException, "Invalid z-coordinate for cutting plane."); >> 106 } 104 } 107 } 105 108 106 ////////////////////////////////////////////// << 109 /////////////////////////////////////////////////////////////////////////////// 107 // 110 // 108 // Fake default constructor - sets only member 111 // Fake default constructor - sets only member data and allocates memory 109 // for usage restri 112 // for usage restricted to object persistency. 110 << 113 // 111 G4EllipticalCone::G4EllipticalCone( __void__& 114 G4EllipticalCone::G4EllipticalCone( __void__& a ) 112 : G4VSolid(a), halfCarTol(0.), << 115 : G4VSolid(a), fpPolyhedron(0), fCubicVolume(0.), fSurfaceArea(0.), 113 xSemiAxis(0.), ySemiAxis(0.), zheight(0.), << 116 zTopCut(0.) 114 cosAxisMin(0.), invXX(0.), invYY(0.) << 115 { 117 { 116 } 118 } 117 119 118 ////////////////////////////////////////////// << 120 /////////////////////////////////////////////////////////////////////////////// 119 // 121 // 120 // Destructor 122 // Destructor 121 << 123 // 122 G4EllipticalCone::~G4EllipticalCone() 124 G4EllipticalCone::~G4EllipticalCone() 123 { 125 { 124 delete fpPolyhedron; fpPolyhedron = nullptr; << 125 } 126 } 126 127 127 ////////////////////////////////////////////// << 128 /////////////////////////////////////////////////////////////////////////////// 128 // 129 // 129 // Copy constructor << 130 // Calculate extent under transform and specified limit 130 << 131 G4EllipticalCone::G4EllipticalCone(const G4Ell << 132 : G4VSolid(rhs), halfCarTol(rhs.halfCarTol), << 133 fCubicVolume(rhs.fCubicVolume), fSurfaceAr << 134 xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.yS << 135 zheight(rhs.zheight), zTopCut(rhs.zTopCut) << 136 cosAxisMin(rhs.cosAxisMin), invXX(rhs.invX << 137 { << 138 } << 139 << 140 ////////////////////////////////////////////// << 141 // 131 // 142 // Assignment operator << 132 G4bool 143 << 133 G4EllipticalCone::CalculateExtent( const EAxis axis, 144 G4EllipticalCone& G4EllipticalCone::operator = << 134 const G4VoxelLimits &voxelLimit, >> 135 const G4AffineTransform &transform, >> 136 G4double &min, G4double &max ) const 145 { 137 { 146 // Check assignment to self << 138 G4SolidExtentList extentList( axis, voxelLimit ); 147 // << 139 148 if (this == &rhs) { return *this; } << 140 // 149 << 141 // We are going to divide up our elliptical face into small pieces 150 // Copy base class data << 142 // 151 // << 143 152 G4VSolid::operator=(rhs); << 144 // 153 << 145 // Choose phi size of our segment(s) based on constants as 154 // Copy data << 146 // defined in meshdefs.hh 155 // << 147 // 156 halfCarTol = rhs.halfCarTol; << 148 G4int numPhi = kMaxMeshSections; 157 fCubicVolume = rhs.fCubicVolume; fSurfaceAr << 149 G4double sigPhi = twopi/numPhi; 158 xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs. << 150 159 zheight = rhs.zheight; zTopCut = rhs.zTopCu << 151 // 160 cosAxisMin = rhs.cosAxisMin; invXX = rhs.in << 152 // We have to be careful to keep our segments completely outside >> 153 // of the elliptical surface. To do so we imagine we have >> 154 // a simple (unit radius) circular cross section (as in G4Tubs) >> 155 // and then "stretch" the dimensions as necessary to fit the ellipse. >> 156 // >> 157 G4double rFudge = 1.0/std::cos(0.5*sigPhi); >> 158 G4double dxFudgeBot = xSemiAxis*2.*zheight*rFudge, >> 159 dyFudgeBot = ySemiAxis*2.*zheight*rFudge; >> 160 G4double dxFudgeTop = xSemiAxis*(zheight-zTopCut)*rFudge, >> 161 dyFudgeTop = ySemiAxis*(zheight-zTopCut)*rFudge; >> 162 >> 163 // >> 164 // As we work around the elliptical surface, we build >> 165 // a "phi" segment on the way, and keep track of two >> 166 // additional polygons for the two ends. >> 167 // >> 168 G4ClippablePolygon endPoly1, endPoly2, phiPoly; >> 169 >> 170 G4double phi = 0, >> 171 cosPhi = std::cos(phi), >> 172 sinPhi = std::sin(phi); >> 173 G4ThreeVector v0( dxFudgeTop*cosPhi, dyFudgeTop*sinPhi, +zTopCut ), >> 174 v1( dxFudgeBot*cosPhi, dyFudgeBot*sinPhi, -zTopCut ), >> 175 w0, w1; >> 176 transform.ApplyPointTransform( v0 ); >> 177 transform.ApplyPointTransform( v1 ); >> 178 do >> 179 { >> 180 phi += sigPhi; >> 181 if (numPhi == 1) phi = 0; // Try to avoid roundoff >> 182 cosPhi = std::cos(phi), >> 183 sinPhi = std::sin(phi); >> 184 >> 185 w0 = G4ThreeVector( dxFudgeTop*cosPhi, dyFudgeTop*sinPhi, +zTopCut ); >> 186 w1 = G4ThreeVector( dxFudgeBot*cosPhi, dyFudgeBot*sinPhi, -zTopCut ); >> 187 transform.ApplyPointTransform( w0 ); >> 188 transform.ApplyPointTransform( w1 ); >> 189 >> 190 // >> 191 // Add a point to our z ends >> 192 // >> 193 endPoly1.AddVertexInOrder( v0 ); >> 194 endPoly2.AddVertexInOrder( v1 ); >> 195 >> 196 // >> 197 // Build phi polygon >> 198 // >> 199 phiPoly.ClearAllVertices(); >> 200 >> 201 phiPoly.AddVertexInOrder( v0 ); >> 202 phiPoly.AddVertexInOrder( v1 ); >> 203 phiPoly.AddVertexInOrder( w1 ); >> 204 phiPoly.AddVertexInOrder( w0 ); >> 205 >> 206 if (phiPoly.PartialClip( voxelLimit, axis )) >> 207 { >> 208 // >> 209 // Get unit normal >> 210 // >> 211 phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() ); >> 212 >> 213 extentList.AddSurface( phiPoly ); >> 214 } 161 215 162 fRebuildPolyhedron = false; << 216 // 163 delete fpPolyhedron; fpPolyhedron = nullptr << 217 // Next vertex >> 218 // >> 219 v0 = w0; >> 220 v1 = w1; >> 221 } while( --numPhi > 0 ); 164 222 165 return *this; << 223 // >> 224 // Process the end pieces >> 225 // >> 226 if (endPoly1.PartialClip( voxelLimit, axis )) >> 227 { >> 228 static const G4ThreeVector normal(0,0,+1); >> 229 endPoly1.SetNormal( transform.TransformAxis(normal) ); >> 230 extentList.AddSurface( endPoly1 ); >> 231 } >> 232 >> 233 if (endPoly2.PartialClip( voxelLimit, axis )) >> 234 { >> 235 static const G4ThreeVector normal(0,0,-1); >> 236 endPoly2.SetNormal( transform.TransformAxis(normal) ); >> 237 extentList.AddSurface( endPoly2 ); >> 238 } >> 239 >> 240 // >> 241 // Return min/max value >> 242 // >> 243 return extentList.GetExtent( min, max ); 166 } 244 } 167 245 168 ////////////////////////////////////////////// << 246 //////////////////////////////////////////////////////////////////////// 169 // 247 // 170 // Get bounding box << 248 // Return whether point inside/outside/on surface 171 << 249 // Split into radius, phi, theta checks 172 void G4EllipticalCone::BoundingLimits(G4ThreeV << 250 // Each check modifies `in', or returns as approprate 173 G4ThreeV << 251 // >> 252 EInside G4EllipticalCone::Inside(const G4ThreeVector& p) const 174 { 253 { 175 G4double zcut = GetZTopCut(); << 254 G4double rad2oo, // outside surface outer tolerance 176 G4double height = GetZMax(); << 255 rad2oi; // outside surface inner tolerance 177 G4double xmax = GetSemiAxisX()*(height+zcu << 256 178 G4double ymax = GetSemiAxisY()*(height+zcu << 257 EInside in; 179 pMin.set(-xmax,-ymax,-zcut); << 180 pMax.set( xmax, ymax, zcut); << 181 << 182 // Check correctness of the bounding box << 183 // << 184 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 185 { << 186 std::ostringstream message; << 187 message << "Bad bounding box (min >= max) << 188 << GetName() << " !" << 189 << "\npMin = " << pMin << 190 << "\npMax = " << pMax; << 191 G4Exception("G4EllipticalCone::BoundingLim << 192 JustWarning, message); << 193 DumpInfo(); << 194 } << 195 } << 196 258 197 ////////////////////////////////////////////// << 259 // check this side of z cut first, because that's fast 198 // << 260 // 199 // Calculate extent under transform and specif << 200 261 201 G4bool << 262 if ( (p.z() < -zTopCut - 0.5*kCarTolerance) 202 G4EllipticalCone::CalculateExtent(const EAxis << 263 || (p.z() > zTopCut + 0.5*kCarTolerance ) ) 203 const G4Voxe << 204 const G4Affi << 205 G4doub << 206 { << 207 G4ThreeVector bmin,bmax; << 208 G4bool exist; << 209 << 210 // Check bounding box (bbox) << 211 // << 212 BoundingLimits(bmin,bmax); << 213 G4BoundingEnvelope bbox(bmin,bmax); << 214 #ifdef G4BBOX_EXTENT << 215 return bbox.CalculateExtent(pAxis,pVoxelLimi << 216 #endif << 217 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 218 { 264 { 219 return exist = pMin < pMax; << 265 return in = kOutside; 220 } 266 } 221 267 222 // Set bounding envelope (benv) and calculat << 268 rad2oo= sqr(p.x()/( xSemiAxis + 0.5*kRadTolerance )) 223 // << 269 + sqr(p.y()/( ySemiAxis + 0.5*kRadTolerance )); 224 static const G4int NSTEPS = 48; // number of << 225 static const G4double ang = twopi/NSTEPS; << 226 static const G4double sinHalf = std::sin(0.5 << 227 static const G4double cosHalf = std::cos(0.5 << 228 static const G4double sinStep = 2.*sinHalf*c << 229 static const G4double cosStep = 1. - 2.*sinH << 230 G4double zcut = bmax.z(); << 231 G4double height = GetZMax(); << 232 G4double sxmin = GetSemiAxisX()*(height-zcu << 233 G4double symin = GetSemiAxisY()*(height-zcu << 234 G4double sxmax = bmax.x()/cosHalf; << 235 G4double symax = bmax.y()/cosHalf; << 236 << 237 G4double sinCur = sinHalf; << 238 G4double cosCur = cosHalf; << 239 G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS << 240 for (G4int k=0; k<NSTEPS; ++k) << 241 { << 242 baseA[k].set(sxmax*cosCur,symax*sinCur,-zc << 243 baseB[k].set(sxmin*cosCur,symin*sinCur, zc << 244 << 245 G4double sinTmp = sinCur; << 246 sinCur = sinCur*cosStep + cosCur*sinStep; << 247 cosCur = cosCur*cosStep - sinTmp*sinStep; << 248 } << 249 << 250 std::vector<const G4ThreeVectorList *> polyg << 251 polygons[0] = &baseA; << 252 polygons[1] = &baseB; << 253 G4BoundingEnvelope benv(bmin,bmax,polygons); << 254 exist = benv.CalculateExtent(pAxis,pVoxelLim << 255 return exist; << 256 } << 257 270 258 ////////////////////////////////////////////// << 271 if ( rad2oo > sqr( zheight-p.z() ) ) 259 // << 272 { 260 // Determine where is point: inside, outside o << 273 return in = kOutside; >> 274 } 261 275 262 EInside G4EllipticalCone::Inside(const G4Three << 276 // rad2oi= sqr( p.x()*(1.0 + 0.5*kRadTolerance/(xSemiAxis*xSemiAxis)) ) 263 { << 277 // + sqr( p.y()*(1.0 + 0.5*kRadTolerance/(ySemiAxis*ySemiAxis)) ); 264 G4double hp = std::sqrt(p.x()*p.x()*invXX + << 278 rad2oi = sqr(p.x()/( xSemiAxis - 0.5*kRadTolerance )) 265 G4double ds = (hp - zheight)*cosAxisMin; << 279 + sqr(p.y()/( ySemiAxis - 0.5*kRadTolerance )); 266 G4double dz = std::abs(p.z()) - zTopCut; << 280 267 G4double dist = std::max(ds,dz); << 281 if (rad2oi < sqr( zheight-p.z() ) ) >> 282 { >> 283 in = ( ( p.z() < -zTopCut + 0.5*kRadTolerance ) >> 284 || ( p.z() > zTopCut - 0.5*kRadTolerance ) ) ? kSurface : kInside; >> 285 } >> 286 else >> 287 { >> 288 in = kSurface; >> 289 } 268 290 269 if (dist > halfCarTol) return kOutside; << 291 return in; 270 return (dist > -halfCarTol) ? kSurface : kIn << 271 } 292 } 272 293 273 ////////////////////////////////////////////// 294 ///////////////////////////////////////////////////////////////////////// 274 // 295 // 275 // Return unit normal at surface closest to p << 296 // Return unit normal of surface closest to p not protected against p=0 276 << 297 // 277 G4ThreeVector G4EllipticalCone::SurfaceNormal( 298 G4ThreeVector G4EllipticalCone::SurfaceNormal( const G4ThreeVector& p) const 278 { 299 { 279 G4ThreeVector norm(0,0,0); << 280 G4int nsurf = 0; // number of surfaces wher << 281 300 282 G4double hp = std::sqrt(p.x()*p.x()*invXX + << 301 G4double rx = sqr(p.x()/xSemiAxis), 283 G4double ds = (hp - zheight)*cosAxisMin; << 302 ry = sqr(p.y()/ySemiAxis); 284 if (std::abs(ds) <= halfCarTol) << 303 >> 304 G4double rad = std::sqrt(rx + ry); >> 305 >> 306 G4ThreeVector norm; >> 307 >> 308 if( (p.z() < -zTopCut) && ((rx+ry) < sqr(zTopCut + zheight)) ) 285 { 309 { 286 norm = G4ThreeVector(p.x()*invXX, p.y()*in << 310 return G4ThreeVector( 0., 0., -1. ); 287 G4double mag = norm.mag(); << 288 if (mag == 0) return {0,0,1}; // apex << 289 norm *= (1/mag); << 290 ++nsurf; << 291 } 311 } 292 G4double dz = std::abs(p.z()) - zTopCut; << 312 293 if (std::abs(dz) <= halfCarTol) << 313 if( (p.z() > (zheight > zTopCut ? zheight : zTopCut)) && >> 314 ((rx+ry) < sqr(zheight-zTopCut)) ) 294 { 315 { 295 norm += G4ThreeVector(0., 0.,(p.z() < 0) ? << 316 return G4ThreeVector( 0., 0., 1. ); 296 ++nsurf; << 297 } 317 } 298 318 299 if (nsurf == 1) return norm; << 319 if( p.z() > rad + 2.*zTopCut - zheight ) 300 else if (nsurf > 1) return norm.unit(); // << 301 else << 302 { 320 { 303 // Point is not on the surface << 321 if ( p.z() > zTopCut ) 304 // << 322 { 305 #ifdef G4CSGDEBUG << 323 if( p.x() == 0. ) 306 std::ostringstream message; << 324 { 307 G4long oldprc = message.precision(16); << 325 norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., 1. ); 308 message << "Point p is not on surface (!?) << 326 return norm /= norm.mag(); 309 << GetName() << G4endl; << 327 } 310 message << "Position:\n"; << 328 if( p.y() == 0. ) 311 message << " p.x() = " << p.x()/mm << " << 329 { 312 message << " p.y() = " << p.y()/mm << " << 330 norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., 1. ); 313 message << " p.z() = " << p.z()/mm << " << 331 return norm /= norm.mag(); 314 G4cout.precision(oldprc); << 332 } 315 G4Exception("G4EllipticalCone::SurfaceNorm << 333 316 JustWarning, message ); << 334 G4double m = std::fabs(p.x()/p.y()); 317 DumpInfo(); << 335 G4double c2 = sqr(zheight-zTopCut)/(1./sqr(xSemiAxis)+sqr(m/ySemiAxis)); 318 #endif << 336 G4double x = std::sqrt(c2); 319 return ApproxSurfaceNormal(p); << 337 G4double y = m*x; >> 338 >> 339 x /= sqr(xSemiAxis); >> 340 y /= sqr(ySemiAxis); >> 341 >> 342 norm = G4ThreeVector( p.x() < 0. ? -x : x, >> 343 p.y() < 0. ? -y : y, >> 344 zheight - zTopCut ); >> 345 norm /= norm.mag(); >> 346 norm += G4ThreeVector( 0., 0., 1. ); >> 347 return norm /= norm.mag(); >> 348 } >> 349 >> 350 return G4ThreeVector( 0., 0., 1. ); 320 } 351 } 321 } << 352 >> 353 if( p.z() < rad - 2.*zTopCut - zheight ) >> 354 { >> 355 if( p.x() == 0. ) >> 356 { >> 357 norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., -1. ); >> 358 return norm /= norm.mag(); >> 359 } >> 360 if( p.y() == 0. ) >> 361 { >> 362 norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., -1. ); >> 363 return norm /= norm.mag(); >> 364 } >> 365 >> 366 G4double m = std::fabs(p.x()/p.y()); >> 367 G4double c2 = sqr(zheight+zTopCut)/(1./sqr(xSemiAxis)+sqr(m/ySemiAxis)); >> 368 G4double x = std::sqrt(c2); >> 369 G4double y = m*x; >> 370 >> 371 x /= sqr(xSemiAxis); >> 372 y /= sqr(ySemiAxis); >> 373 >> 374 norm = G4ThreeVector( p.x() < 0. ? -x : x, >> 375 p.y() < 0. ? -y : y, >> 376 zheight - zTopCut ); >> 377 norm /= norm.mag(); >> 378 norm += G4ThreeVector( 0., 0., -1. ); >> 379 return norm /= norm.mag(); >> 380 } >> 381 >> 382 norm = G4ThreeVector(p.x()/sqr(xSemiAxis), p.y()/sqr(ySemiAxis), rad); >> 383 >> 384 G4double m = std::tan(pi/8.); >> 385 G4double c = -zTopCut - m*(zTopCut + zheight); 322 386 323 ////////////////////////////////////////////// << 387 if( p.z() < -m*rad + c ) 324 // << 388 return G4ThreeVector (0.,0.,-1.); 325 // Find surface nearest to point and return co << 389 326 // The algorithm is similar to the algorithm u << 390 return norm /= norm.mag(); 327 // This method normally should not be called. << 328 << 329 G4ThreeVector << 330 G4EllipticalCone::ApproxSurfaceNormal(const G4 << 331 { << 332 G4double hp = std::sqrt(p.x()*p.x()*invXX + << 333 G4double ds = (hp - zheight)*cosAxisMin; << 334 G4double dz = std::abs(p.z()) - zTopCut; << 335 if (ds > dz && std::abs(hp - p.z()) > halfCa << 336 return G4ThreeVector(p.x()*invXX, p.y()*in << 337 else << 338 return { 0., 0., (G4double)((p.z() < 0) ? << 339 } 391 } 340 392 341 ////////////////////////////////////////////// << 393 ////////////////////////////////////////////////////////////////////////// 342 // 394 // 343 // Calculate distance to shape from outside, a 395 // Calculate distance to shape from outside, along normalised vector 344 // return kInfinity if no intersection, or int 396 // return kInfinity if no intersection, or intersection distance <= tolerance 345 << 397 // 346 G4double G4EllipticalCone::DistanceToIn( const 398 G4double G4EllipticalCone::DistanceToIn( const G4ThreeVector& p, 347 const 399 const G4ThreeVector& v ) const 348 { 400 { >> 401 >> 402 static const G4double halfTol = 0.5*kCarTolerance; >> 403 349 G4double distMin = kInfinity; 404 G4double distMin = kInfinity; 350 405 351 // code from EllipticalTube 406 // code from EllipticalTube 352 407 353 G4double sigz = p.z()+zTopCut; 408 G4double sigz = p.z()+zTopCut; 354 409 355 // 410 // 356 // Check z = -dz planer surface 411 // Check z = -dz planer surface 357 // 412 // 358 413 359 if (sigz < halfCarTol) << 414 if (sigz < halfTol) 360 { 415 { 361 // 416 // 362 // We are "behind" the shape in z, and so 417 // We are "behind" the shape in z, and so can 363 // potentially hit the rear face. Correct 418 // potentially hit the rear face. Correct direction? 364 // 419 // 365 if (v.z() <= 0) 420 if (v.z() <= 0) 366 { 421 { 367 // 422 // 368 // As long as we are far enough away, we 423 // As long as we are far enough away, we know we 369 // can't intersect 424 // can't intersect 370 // 425 // 371 if (sigz < 0) return kInfinity; 426 if (sigz < 0) return kInfinity; 372 427 373 // 428 // 374 // Otherwise, we don't intersect unless 429 // Otherwise, we don't intersect unless we are 375 // on the surface of the ellipse 430 // on the surface of the ellipse 376 // 431 // 377 432 378 if ( sqr(p.x()/( xSemiAxis - halfCarTol << 433 if ( sqr(p.x()/( xSemiAxis - halfTol )) 379 + sqr(p.y()/( ySemiAxis - halfCarTol << 434 + sqr(p.y()/( ySemiAxis - halfTol )) <= sqr( zheight+zTopCut ) ) 380 return kInfinity; 435 return kInfinity; 381 436 382 } 437 } 383 else 438 else 384 { 439 { 385 // 440 // 386 // How far? 441 // How far? 387 // 442 // 388 G4double q = -sigz/v.z(); << 443 G4double s = -sigz/v.z(); 389 444 390 // 445 // 391 // Where does that place us? 446 // Where does that place us? 392 // 447 // 393 G4double xi = p.x() + q*v.x(), << 448 G4double xi = p.x() + s*v.x(), 394 yi = p.y() + q*v.y(); << 449 yi = p.y() + s*v.y(); 395 450 396 // 451 // 397 // Is this on the surface (within ellips 452 // Is this on the surface (within ellipse)? 398 // 453 // 399 if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxi 454 if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight + zTopCut ) ) 400 { 455 { 401 // 456 // 402 // Yup. Return q, unless we are on the << 457 // Yup. Return s, unless we are on the surface 403 // 458 // 404 return (sigz < -halfCarTol) ? q : 0; << 459 return (sigz < -halfTol) ? s : 0; 405 } 460 } 406 else if (xi/(xSemiAxis*xSemiAxis)*v.x() 461 else if (xi/(xSemiAxis*xSemiAxis)*v.x() 407 + yi/(ySemiAxis*ySemiAxis)*v.y() 462 + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0) 408 { 463 { 409 // 464 // 410 // Else, if we are traveling outwards, 465 // Else, if we are traveling outwards, we know 411 // we must miss 466 // we must miss 412 // 467 // 413 // return kInfinity; 468 // return kInfinity; 414 } 469 } 415 } 470 } 416 } 471 } 417 472 418 // 473 // 419 // Check z = +dz planer surface 474 // Check z = +dz planer surface 420 // 475 // >> 476 421 sigz = p.z() - zTopCut; 477 sigz = p.z() - zTopCut; 422 478 423 if (sigz > -halfCarTol) << 479 if (sigz > -halfTol) 424 { 480 { 425 if (v.z() >= 0) 481 if (v.z() >= 0) 426 { 482 { 427 483 428 if (sigz > 0) return kInfinity; 484 if (sigz > 0) return kInfinity; 429 485 430 if ( sqr(p.x()/( xSemiAxis - halfCarTol << 486 if ( sqr(p.x()/( xSemiAxis - halfTol )) 431 + sqr(p.y()/( ySemiAxis - halfCarTol << 487 + sqr(p.y()/( ySemiAxis - halfTol )) <= sqr( zheight-zTopCut ) ) 432 return kInfinity; 488 return kInfinity; 433 489 434 } 490 } 435 else { 491 else { 436 G4double q = -sigz/v.z(); << 492 G4double s = -sigz/v.z(); 437 493 438 G4double xi = p.x() + q*v.x(), << 494 G4double xi = p.x() + s*v.x(), 439 yi = p.y() + q*v.y(); << 495 yi = p.y() + s*v.y(); 440 496 441 if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxi 497 if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight - zTopCut ) ) 442 { 498 { 443 return (sigz > -halfCarTol) ? q : 0; << 499 return (sigz > -halfTol) ? s : 0; 444 } 500 } 445 else if (xi/(xSemiAxis*xSemiAxis)*v.x() 501 else if (xi/(xSemiAxis*xSemiAxis)*v.x() 446 + yi/(ySemiAxis*ySemiAxis)*v.y() 502 + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0) 447 { 503 { 448 // return kInfinity; 504 // return kInfinity; 449 } 505 } 450 } 506 } 451 } 507 } 452 508 453 509 454 #if 0 510 #if 0 455 511 456 // check to see if Z plane is relevant 512 // check to see if Z plane is relevant 457 // 513 // 458 if (p.z() < -zTopCut - halfCarTol) << 514 if (p.z() < -zTopCut - 0.5*kCarTolerance) 459 { 515 { 460 if (v.z() <= 0.0) 516 if (v.z() <= 0.0) 461 return distMin; 517 return distMin; 462 518 463 G4double lambda = (-zTopCut - p.z())/v.z() 519 G4double lambda = (-zTopCut - p.z())/v.z(); 464 520 465 if ( sqr((lambda*v.x()+p.x())/xSemiAxis) + 521 if ( sqr((lambda*v.x()+p.x())/xSemiAxis) + 466 sqr((lambda*v.y()+p.y())/ySemiAxis) < 522 sqr((lambda*v.y()+p.y())/ySemiAxis) <= 467 sqr(zTopCut + zheight + halfCarTol) ) << 523 sqr(zTopCut + zheight + 0.5*kRadTolerance) ) 468 { 524 { 469 return distMin = std::fabs(lambda); 525 return distMin = std::fabs(lambda); 470 } 526 } 471 } 527 } 472 528 473 if (p.z() > zTopCut + halfCarTol) << 529 if (p.z() > zTopCut+0.5*kCarTolerance) 474 { 530 { 475 if (v.z() >= 0.0) 531 if (v.z() >= 0.0) 476 { return distMin; } 532 { return distMin; } 477 533 478 G4double lambda = (zTopCut - p.z()) / v.z 534 G4double lambda = (zTopCut - p.z()) / v.z(); 479 535 480 if ( sqr((lambda*v.x() + p.x())/xSemiAxis) 536 if ( sqr((lambda*v.x() + p.x())/xSemiAxis) + 481 sqr((lambda*v.y() + p.y())/ySemiAxis) 537 sqr((lambda*v.y() + p.y())/ySemiAxis) <= 482 sqr(zheight - zTopCut + halfCarTol) ) << 538 sqr(zheight - zTopCut + 0.5*kRadTolerance) ) 483 { 539 { 484 return distMin = std::fabs(lambda); 540 return distMin = std::fabs(lambda); 485 } 541 } 486 } 542 } 487 543 488 if (p.z() > zTopCut - halfCarTol << 544 if (p.z() > zTopCut - 0.5*kCarTolerance 489 && p.z() < zTopCut + halfCarTol ) << 545 && p.z() < zTopCut + 0.5*kCarTolerance ) 490 { 546 { 491 if (v.z() > 0.) 547 if (v.z() > 0.) 492 { return kInfinity; } 548 { return kInfinity; } 493 549 494 return distMin = 0.; 550 return distMin = 0.; 495 } 551 } 496 552 497 if (p.z() < -zTopCut + halfCarTol << 553 if (p.z() < -zTopCut + 0.5*kCarTolerance 498 && p.z() > -zTopCut - halfCarTol) << 554 && p.z() > -zTopCut - 0.5*kCarTolerance) 499 { 555 { 500 if (v.z() < 0.) 556 if (v.z() < 0.) 501 { return distMin = kInfinity; } 557 { return distMin = kInfinity; } 502 558 503 return distMin = 0.; 559 return distMin = 0.; 504 } 560 } 505 561 506 #endif 562 #endif 507 563 508 // if we are here then it either intersects 564 // if we are here then it either intersects or grazes the curved surface 509 // or it does not intersect at all 565 // or it does not intersect at all 510 // 566 // >> 567 511 G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y( 568 G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z()); 512 G4double B = 2*(v.x()*p.x()/sqr(xSemiAxis) + 569 G4double B = 2*(v.x()*p.x()/sqr(xSemiAxis) + 513 v.y()*p.y()/sqr(ySemiAxis) + 570 v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z())); 514 G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y( 571 G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) - 515 sqr(zheight - p.z()); 572 sqr(zheight - p.z()); 516 573 517 G4double discr = B*B - 4.*A*C; 574 G4double discr = B*B - 4.*A*C; 518 575 519 // if the discriminant is negative it never 576 // if the discriminant is negative it never hits the curved object 520 // 577 // 521 if ( discr < -halfCarTol ) << 578 if ( discr < -0.5*kCarTolerance ) 522 { return distMin; } 579 { return distMin; } 523 580 524 // case below is when it hits or grazes the << 581 //case below is when it hits or grazes the surface 525 // 582 // 526 if ( (discr >= -halfCarTol ) && (discr < hal << 583 if ( (discr >= - 0.5*kCarTolerance ) && (discr < 0.5*kCarTolerance ) ) 527 { 584 { 528 return distMin = std::fabs(-B/(2.*A)); 585 return distMin = std::fabs(-B/(2.*A)); 529 } 586 } 530 587 531 G4double plus = (-B+std::sqrt(discr))/(2.*A 588 G4double plus = (-B+std::sqrt(discr))/(2.*A); 532 G4double minus = (-B-std::sqrt(discr))/(2.*A 589 G4double minus = (-B-std::sqrt(discr))/(2.*A); 533 << 590 // G4double lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus; 534 // Special case::Point on Surface, Check nor << 535 << 536 if ( ( std::fabs(plus) < halfCarTol )||( std << 537 { << 538 G4ThreeVector truenorm(p.x()/(xSemiAxis*xS << 539 p.y()/(ySemiAxis*yS << 540 -( p.z() - zheight << 541 if ( truenorm*v >= 0) // going outside t << 542 { << 543 return kInfinity; << 544 } << 545 else << 546 { << 547 return 0; << 548 } << 549 } << 550 591 551 // G4double lambda = std::fabs(plus) < std:: << 552 G4double lambda = 0; 592 G4double lambda = 0; 553 593 554 if ( minus > halfCarTol && minus < distMin ) << 594 if ( minus > halfTol && minus < distMin ) 555 { 595 { 556 lambda = minus ; 596 lambda = minus ; 557 // check normal vector n * v < 0 597 // check normal vector n * v < 0 558 G4ThreeVector pin = p + lambda*v; 598 G4ThreeVector pin = p + lambda*v; 559 if(std::fabs(pin.z())< zTopCut + halfCarTo << 599 560 { << 600 G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis), 561 G4ThreeVector truenorm(pin.x()/(xSemiAxi << 601 pin.y()/(ySemiAxis*ySemiAxis), 562 pin.y()/(ySemiAxi << 602 - ( pin.z() - zheight )); 563 - ( pin.z() - zhe << 603 if ( truenorm*v < 0) 564 if ( truenorm*v < 0) << 604 { // yes, going inside the solid 565 { // yes, going inside the solid << 605 distMin = lambda; 566 distMin = lambda; << 567 } << 568 } 606 } 569 } 607 } 570 if ( plus > halfCarTol && plus < distMin ) << 608 >> 609 if ( plus > halfTol && plus < distMin ) 571 { 610 { 572 lambda = plus ; 611 lambda = plus ; 573 // check normal vector n * v < 0 612 // check normal vector n * v < 0 574 G4ThreeVector pin = p + lambda*v; 613 G4ThreeVector pin = p + lambda*v; 575 if(std::fabs(pin.z()) < zTopCut + halfCarT << 614 576 { << 615 G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis), 577 G4ThreeVector truenorm(pin.x()/(xSemiAxi << 616 pin.y()/(ySemiAxis*ySemiAxis), 578 pin.y()/(ySemiAxi << 617 - ( pin.z() - zheight ) ); 579 - ( pin.z() - zhe << 618 if ( truenorm*v < 0) 580 if ( truenorm*v < 0) << 619 { // yes, going inside the solid 581 { // yes, going inside the solid << 620 distMin = lambda; 582 distMin = lambda; << 583 } << 584 } 621 } 585 } 622 } 586 if (distMin < halfCarTol) distMin=0.; << 623 >> 624 #ifdef G4SPECSDEBUG >> 625 // G4cout << "DToIn: plus,minus, lambda = " << plus >> 626 // << ", " << minus << ", " << lambda << G4endl ; >> 627 // G4cout << "DToIn: distMin = " << distMin << G4endl ; >> 628 #endif >> 629 >> 630 587 return distMin ; 631 return distMin ; 588 } 632 } 589 633 590 ////////////////////////////////////////////// << 634 ////////////////////////////////////////////////////////////////////////// 591 // 635 // 592 // Calculate distance (<= actual) to closest s 636 // Calculate distance (<= actual) to closest surface of shape from outside 593 // Return 0 if point inside 637 // Return 0 if point inside 594 << 638 // 595 G4double G4EllipticalCone::DistanceToIn(const 639 G4double G4EllipticalCone::DistanceToIn(const G4ThreeVector& p) const 596 { 640 { 597 G4double hp = std::sqrt(p.x()*p.x()*invXX + << 641 G4double distR, distR2, distZ, maxDim; 598 G4double ds = (hp - zheight)*cosAxisMin; << 642 G4double distRad; 599 G4double dz = std::abs(p.z()) - zTopCut; << 643 600 G4double dist = std::max(ds,dz); << 644 // check if the point lies either below z=-zTopCut in bottom elliptical 601 return (dist > 0) ? dist : 0.; << 645 // region or on top within cut elliptical region >> 646 // >> 647 if( (p.z() <= -zTopCut) && (sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) >> 648 <= sqr(zTopCut + zheight + 0.5*kCarTolerance )) ) >> 649 { >> 650 //return distZ = std::fabs(zTopCut - p.z()); >> 651 return distZ = std::fabs(zTopCut + p.z()); >> 652 } >> 653 >> 654 if( (p.z() >= zTopCut) && (sqr(p.x()/xSemiAxis)+sqr(p.y()/ySemiAxis) >> 655 <= sqr(zheight - zTopCut + kCarTolerance/2.0 )) ) >> 656 { >> 657 return distZ = std::fabs(p.z() - zTopCut); >> 658 } >> 659 >> 660 // below we use the following approximation: we take the largest of the >> 661 // axes and find the shortest distance to the circular (cut) cone of that >> 662 // radius. >> 663 // >> 664 maxDim = xSemiAxis >= ySemiAxis ? xSemiAxis:ySemiAxis; >> 665 distRad = std::sqrt(p.x()*p.x()+p.y()*p.y()); >> 666 >> 667 if( p.z() > maxDim*distRad + zTopCut*(1.+maxDim)-sqr(maxDim)*zheight ) >> 668 { >> 669 distR2 = sqr(p.z() - zTopCut) + sqr(distRad - maxDim*(zheight - zTopCut)); >> 670 return std::sqrt( distR2 ); >> 671 } >> 672 >> 673 if( distRad > maxDim*( zheight - p.z() ) ) >> 674 { >> 675 if( p.z() > maxDim*distRad - (zTopCut*(1.+maxDim)+sqr(maxDim)*zheight) ) >> 676 { >> 677 G4double zVal = (p.z()-maxDim*(distRad-maxDim*zheight))/(1.+sqr(maxDim)); >> 678 G4double rVal = maxDim*(zheight - zVal); >> 679 return distR = std::sqrt(sqr(p.z() - zVal) + sqr(distRad - rVal)); >> 680 } >> 681 } >> 682 >> 683 if( distRad <= maxDim*(zheight - p.z()) ) >> 684 { >> 685 distR2 = sqr(distRad - maxDim*(zheight + zTopCut)) + sqr(p.z() + zTopCut); >> 686 return std::sqrt( distR2 ); >> 687 } >> 688 >> 689 return distR = 0; 602 } 690 } 603 691 604 ////////////////////////////////////////////// << 692 ///////////////////////////////////////////////////////////////////////// 605 // 693 // 606 // Calculate distance to surface of shape from 694 // Calculate distance to surface of shape from `inside', 607 // allowing for tolerance 695 // allowing for tolerance 608 << 696 // 609 G4double G4EllipticalCone::DistanceToOut(const 697 G4double G4EllipticalCone::DistanceToOut(const G4ThreeVector& p, 610 const 698 const G4ThreeVector& v, 611 const 699 const G4bool calcNorm, 612 << 700 G4bool *validNorm, 613 << 701 G4ThreeVector *n ) const 614 { 702 { 615 G4double distMin, lambda; 703 G4double distMin, lambda; 616 enum surface_e {kPlaneSurf, kCurvedSurf, kNo 704 enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface; 617 705 618 distMin = kInfinity; 706 distMin = kInfinity; 619 surface = kNoSurf; 707 surface = kNoSurf; 620 708 >> 709 #ifdef G4SPECSDEBUG >> 710 G4cout << "DToOut: vz < 0" << G4endl ; >> 711 #endif >> 712 621 if (v.z() < 0.0) 713 if (v.z() < 0.0) 622 { 714 { 623 lambda = (-p.z() - zTopCut)/v.z(); 715 lambda = (-p.z() - zTopCut)/v.z(); 624 716 625 if ( (sqr((p.x() + lambda*v.x())/xSemiAxis 717 if ( (sqr((p.x() + lambda*v.x())/xSemiAxis) + 626 sqr((p.y() + lambda*v.y())/ySemiAxis 718 sqr((p.y() + lambda*v.y())/ySemiAxis)) < 627 sqr(zheight + zTopCut + halfCarTol) << 719 sqr(zheight + zTopCut + 0.5*kCarTolerance) ) 628 { 720 { 629 distMin = std::fabs(lambda); 721 distMin = std::fabs(lambda); 630 722 631 if (!calcNorm) { return distMin; } 723 if (!calcNorm) { return distMin; } 632 } 724 } 633 distMin = std::fabs(lambda); 725 distMin = std::fabs(lambda); 634 surface = kPlaneSurf; 726 surface = kPlaneSurf; 635 } 727 } 636 728 >> 729 #ifdef G4SPECSDEBUG >> 730 G4cout << "DToOut: vz > 0" << G4endl ; >> 731 #endif >> 732 637 if (v.z() > 0.0) 733 if (v.z() > 0.0) 638 { 734 { 639 lambda = (zTopCut - p.z()) / v.z(); 735 lambda = (zTopCut - p.z()) / v.z(); 640 736 641 if ( (sqr((p.x() + lambda*v.x())/xSemiAxis 737 if ( (sqr((p.x() + lambda*v.x())/xSemiAxis) 642 + sqr((p.y() + lambda*v.y())/ySemiAxis 738 + sqr((p.y() + lambda*v.y())/ySemiAxis) ) 643 < (sqr(zheight - zTopCut + halfCarTol)) << 739 < (sqr(zheight - zTopCut + 0.5*kCarTolerance)) ) 644 { 740 { 645 distMin = std::fabs(lambda); 741 distMin = std::fabs(lambda); 646 if (!calcNorm) { return distMin; } 742 if (!calcNorm) { return distMin; } 647 } 743 } 648 distMin = std::fabs(lambda); 744 distMin = std::fabs(lambda); 649 surface = kPlaneSurf; 745 surface = kPlaneSurf; 650 } 746 } 651 747 652 // if we are here then it either intersects 748 // if we are here then it either intersects or grazes the 653 // curved surface... 749 // curved surface... 654 // 750 // 655 G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y( 751 G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z()); 656 G4double B = 2.*(v.x()*p.x()/sqr(xSemiAxis) 752 G4double B = 2.*(v.x()*p.x()/sqr(xSemiAxis) + 657 v.y()*p.y()/sqr(ySemiAxis) 753 v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z())); 658 G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y( 754 G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) 659 - sqr(zheight - p.z()); 755 - sqr(zheight - p.z()); 660 756 661 G4double discr = B*B - 4.*A*C; 757 G4double discr = B*B - 4.*A*C; 662 758 663 if ( discr >= - halfCarTol && discr < halfCa << 759 if ( discr >= - 0.5*kCarTolerance && discr < 0.5*kCarTolerance ) 664 { 760 { 665 if(!calcNorm) { return distMin = std::fabs 761 if(!calcNorm) { return distMin = std::fabs(-B/(2.*A)); } 666 } 762 } 667 763 668 else if ( discr > halfCarTol ) << 764 else if ( discr > 0.5*kCarTolerance ) 669 { 765 { 670 G4double plus = (-B+std::sqrt(discr))/(2. 766 G4double plus = (-B+std::sqrt(discr))/(2.*A); 671 G4double minus = (-B-std::sqrt(discr))/(2. 767 G4double minus = (-B-std::sqrt(discr))/(2.*A); 672 768 673 if ( plus > halfCarTol && minus > halfCarT << 769 if ( plus > 0.5*kCarTolerance && minus > 0.5*kCarTolerance ) 674 { 770 { 675 // take the shorter distance 771 // take the shorter distance 676 // 772 // 677 lambda = std::fabs(plus) < std::fabs(m 773 lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus; 678 } 774 } 679 else 775 else 680 { 776 { 681 // at least one solution is close to zer 777 // at least one solution is close to zero or negative 682 // so, take small positive solution or z 778 // so, take small positive solution or zero 683 // 779 // 684 lambda = plus > -halfCarTol ? plus : 0 << 780 lambda = plus > -0.5*kCarTolerance ? plus : 0; 685 } 781 } 686 782 687 if ( std::fabs(lambda) < distMin ) 783 if ( std::fabs(lambda) < distMin ) 688 { 784 { 689 if( std::fabs(lambda) > halfCarTol) << 785 distMin = std::fabs(lambda); 690 { << 786 surface = kCurvedSurf; 691 distMin = std::fabs(lambda); << 692 surface = kCurvedSurf; << 693 } << 694 else // Point is On the Surface, Check << 695 { << 696 G4ThreeVector truenorm(p.x()/(xSemiAxi << 697 p.y()/(ySemiAxi << 698 -( p.z() - zhei << 699 if( truenorm.dot(v) > 0 ) << 700 { << 701 distMin = 0.0; << 702 surface = kCurvedSurf; << 703 } << 704 } << 705 } 787 } 706 } 788 } 707 789 708 // set normal if requested 790 // set normal if requested 709 // 791 // 710 if (calcNorm) 792 if (calcNorm) 711 { 793 { 712 if (surface == kNoSurf) 794 if (surface == kNoSurf) 713 { 795 { 714 *validNorm = false; 796 *validNorm = false; 715 } 797 } 716 else 798 else 717 { 799 { 718 *validNorm = true; 800 *validNorm = true; 719 switch (surface) 801 switch (surface) 720 { 802 { 721 case kPlaneSurf: 803 case kPlaneSurf: 722 { 804 { 723 *n = G4ThreeVector(0.,0.,(v.z() > 0. 805 *n = G4ThreeVector(0.,0.,(v.z() > 0.0 ? 1. : -1.)); 724 } 806 } 725 break; 807 break; 726 808 727 case kCurvedSurf: 809 case kCurvedSurf: 728 { 810 { 729 G4ThreeVector pexit = p + distMin*v; 811 G4ThreeVector pexit = p + distMin*v; 730 G4ThreeVector truenorm( pexit.x()/(x << 812 G4ThreeVector truenorm(pexit.x()/(xSemiAxis*xSemiAxis), 731 pexit.y()/(y << 813 pexit.y()/(ySemiAxis*ySemiAxis), 732 -( pexit.z() << 814 pexit.z() - zheight ); 733 truenorm /= truenorm.mag(); 815 truenorm /= truenorm.mag(); 734 *n= truenorm; 816 *n= truenorm; 735 } 817 } 736 break; 818 break; 737 819 738 default: // Should never re << 820 default: >> 821 G4cout.precision(16); >> 822 G4cout << G4endl; 739 DumpInfo(); 823 DumpInfo(); 740 std::ostringstream message; << 824 G4cout << "Position:" << G4endl << G4endl; 741 G4long oldprc = message.precision(16 << 825 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl; 742 message << "Undefined side for valid << 826 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl; 743 << G4endl << 827 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl; 744 << "Position:" << G4endl << 828 G4cout << "Direction:" << G4endl << G4endl; 745 << " p.x() = " << p.x()/ << 829 G4cout << "v.x() = " << v.x() << G4endl; 746 << " p.y() = " << p.y()/ << 830 G4cout << "v.y() = " << v.y() << G4endl; 747 << " p.z() = " << p.z()/ << 831 G4cout << "v.z() = " << v.z() << G4endl << G4endl; 748 << "Direction:" << G4endl << 832 G4cout << "Proposed distance :" << G4endl << G4endl; 749 << " v.x() = " << v.x() << 833 G4cout << "distMin = " << distMin/mm << " mm" << G4endl << G4endl; 750 << " v.y() = " << v.y() << 751 << " v.z() = " << v.z() << 752 << "Proposed distance :" << << 753 << " distMin = " << dis << 754 message.precision(oldprc); << 755 G4Exception("G4EllipticalCone::Dista 834 G4Exception("G4EllipticalCone::DistanceToOut(p,v,..)", 756 "GeomSolids1002", JustWa << 835 "Notification", JustWarning, >> 836 "Undefined side for valid surface normal to solid."); 757 break; 837 break; 758 } 838 } 759 } 839 } 760 } 840 } 761 841 762 if (distMin < halfCarTol) { distMin=0; } << 763 << 764 return distMin; 842 return distMin; 765 } 843 } 766 844 767 ////////////////////////////////////////////// 845 ///////////////////////////////////////////////////////////////////////// 768 // 846 // 769 // Calculate distance (<=actual) to closest su 847 // Calculate distance (<=actual) to closest surface of shape from inside 770 << 848 // 771 G4double G4EllipticalCone::DistanceToOut(const 849 G4double G4EllipticalCone::DistanceToOut(const G4ThreeVector& p) const 772 { 850 { >> 851 G4double rad,roo,roo1, distR, distZ, distMin=0.; >> 852 G4double minAxis = xSemiAxis < ySemiAxis ? xSemiAxis : ySemiAxis; >> 853 773 #ifdef G4SPECSDEBUG 854 #ifdef G4SPECSDEBUG 774 if( Inside(p) == kOutside ) 855 if( Inside(p) == kOutside ) 775 { 856 { 776 std::ostringstream message; << 857 G4cout.precision(16) ; 777 G4long oldprc = message.precision(16); << 858 G4cout << G4endl ; 778 message << "Point p is outside (!?) of so << 779 << "Position:\n" << 780 << " p.x() = " << p.x()/mm << << 781 << " p.y() = " << p.y()/mm << << 782 << " p.z() = " << p.z()/mm << << 783 message.precision(oldprc) ; << 784 G4Exception("G4Ellipsoid::DistanceToOut(p << 785 JustWarning, message); << 786 DumpInfo(); 859 DumpInfo(); >> 860 G4cout << "Position:" << G4endl << G4endl ; >> 861 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; >> 862 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; >> 863 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; >> 864 G4Exception("G4Ellipsoid::DistanceToOut(p)", "Notification", JustWarning, >> 865 "Point p is outside !?" ); 787 } 866 } 788 #endif 867 #endif 789 G4double hp = std::sqrt(p.x()*p.x()*invXX + << 868 790 G4double ds = (zheight - hp)*cosAxisMin; << 869 // since we have made the above warning, below we are working assuming p 791 G4double dz = zTopCut - std::abs(p.z()); << 870 // is inside check how close it is to the circular cone with radius equal 792 G4double dist = std::min(ds,dz); << 871 // to the smaller of the axes 793 return (dist > 0) ? dist : 0.; << 872 // 794 } << 873 if( sqr(p.x()/minAxis)+sqr(p.y()/minAxis) < sqr(zheight - p.z()) ) >> 874 { >> 875 rad = std::sqrt(sqr(p.x()) + sqr(p.y())); >> 876 roo = minAxis*(zheight-p.z()); // radius of cone at z= p.z() >> 877 roo1 = minAxis*(zheight-zTopCut); // radius of cone at z=+zTopCut 795 878 796 ////////////////////////////////////////////// << 879 distZ=zTopCut - std::fabs(p.z()) ; 797 // << 880 distR=(roo-rad)/(std::sqrt(1+sqr(minAxis))); 798 // GetEntityType << 799 881 800 G4GeometryType G4EllipticalCone::GetEntityType << 882 if(rad>roo1) 801 { << 883 { 802 return {"G4EllipticalCone"}; << 884 distMin=(zTopCut-p.z())*(roo-rad)/(roo-roo1); >> 885 distMin=std::min(distMin,distR); >> 886 } >> 887 distMin=std::min(distR,distZ); >> 888 } >> 889 >> 890 return distMin; 803 } 891 } 804 892 805 ////////////////////////////////////////////// << 893 ////////////////////////////////////////////////////////////////////////// 806 // 894 // 807 // Make a clone of the object << 895 // GetEntityType 808 << 896 // 809 G4VSolid* G4EllipticalCone::Clone() const << 897 G4GeometryType G4EllipticalCone::GetEntityType() const 810 { 898 { 811 return new G4EllipticalCone(*this); << 899 return G4String("G4EllipticalCone"); 812 } 900 } 813 901 814 ////////////////////////////////////////////// << 902 ////////////////////////////////////////////////////////////////////////// 815 // 903 // 816 // Stream object contents to an output stream 904 // Stream object contents to an output stream 817 << 905 // 818 std::ostream& G4EllipticalCone::StreamInfo( st 906 std::ostream& G4EllipticalCone::StreamInfo( std::ostream& os ) const 819 { 907 { 820 G4long oldprc = os.precision(16); << 821 os << "------------------------------------- 908 os << "-----------------------------------------------------------\n" 822 << " *** Dump for solid - " << GetName 909 << " *** Dump for solid - " << GetName() << " ***\n" 823 << " ================================= 910 << " ===================================================\n" 824 << " Solid type: G4EllipticalCone\n" 911 << " Solid type: G4EllipticalCone\n" 825 << " Parameters: \n" 912 << " Parameters: \n" 826 913 827 << " semi-axis x: " << xSemiAxis/mm << 914 << " semi-axis x: " << xSemiAxis/mm << " mm \n" 828 << " semi-axis y: " << ySemiAxis/mm << 915 << " semi-axis y: " << ySemiAxis/mm << " mm \n" 829 << " height z: " << zheight/mm << " 916 << " height z: " << zheight/mm << " mm \n" 830 << " half length in z: " << zTopCut/m 917 << " half length in z: " << zTopCut/mm << " mm \n" 831 << "------------------------------------- 918 << "-----------------------------------------------------------\n"; 832 os.precision(oldprc); << 833 919 834 return os; 920 return os; 835 } 921 } 836 922 837 ////////////////////////////////////////////// 923 ///////////////////////////////////////////////////////////////////////// 838 // 924 // 839 // Return random point on the surface of the s << 925 // GetPointOnSurface 840 << 926 // >> 927 // returns quasi-uniformly distributed point on surface of elliptical cone >> 928 // 841 G4ThreeVector G4EllipticalCone::GetPointOnSurf 929 G4ThreeVector G4EllipticalCone::GetPointOnSurface() const 842 { 930 { 843 G4double x0 = xSemiAxis*zheight; // x semi a << 844 G4double y0 = ySemiAxis*zheight; // y semi a << 845 G4double s0 = G4GeomTools::EllipticConeLater << 846 G4double kmin = (zTopCut >= zheight ) ? 0. : << 847 G4double kmax = (zTopCut >= zheight ) ? 2. : << 848 << 849 // Set areas (base at -Z, side surface, base << 850 // << 851 G4double szmin = pi*x0*y0*kmax*kmax; << 852 G4double szmax = pi*x0*y0*kmin*kmin; << 853 G4double sside = s0*(kmax*kmax - kmin*kmin) << 854 G4double ssurf[3] = { szmin, sside, szmax }; << 855 for (auto i=1; i<3; ++i) { ssurf[i] += ssurf << 856 931 857 // Select surface << 932 G4double phi, sinphi, cosphi, aOne, aTwo, aThree, 858 // << 933 chose, zRand, rRand1, rRand2; 859 G4double select = ssurf[2]*G4UniformRand(); << 934 860 G4int k = 2; << 935 G4double rOne = std::sqrt(sqr(xSemiAxis) 861 if (select <= ssurf[1]) k = 1; << 936 + sqr(ySemiAxis))*(zheight - zTopCut); 862 if (select <= ssurf[0]) k = 0; << 937 G4double rTwo = std::sqrt(sqr(xSemiAxis) >> 938 + sqr(ySemiAxis))*(zheight + zTopCut); >> 939 >> 940 aOne = pi*(rOne + rTwo)*std::sqrt(sqr(rOne - rTwo)+sqr(2.*zTopCut)); >> 941 aTwo = pi*xSemiAxis*ySemiAxis*sqr(zheight+zTopCut); >> 942 aThree = pi*xSemiAxis*ySemiAxis*sqr(zheight-zTopCut); >> 943 >> 944 phi = RandFlat::shoot(0.,twopi); >> 945 cosphi = std::cos(phi); >> 946 sinphi = std::sin(phi); >> 947 >> 948 if(zTopCut >= zheight) aThree = 0.; 863 949 864 // Pick random point on selected surface << 950 chose = RandFlat::shoot(0.,aOne+aTwo+aThree); 865 // << 951 if((chose>=0.) && (chose<aOne)) 866 G4ThreeVector p; << 867 switch(k) << 868 { 952 { 869 case 0: // base at -Z, uniform distributio << 953 zRand = RandFlat::shoot(-zTopCut,zTopCut); 870 { << 954 return G4ThreeVector(xSemiAxis*(zheight-zRand)*cosphi, 871 G4double zh = zheight + zTopCut; << 955 ySemiAxis*(zheight-zRand)*sinphi,zRand); 872 G4TwoVector rho = G4RandomPointInEllipse << 956 } 873 p.set(rho.x(),rho.y(),-zTopCut); << 957 else if((chose>=aOne) && (chose<aOne+aTwo)) 874 break; << 958 { 875 } << 959 do 876 case 1: // side surface, uniform distribut << 877 { 960 { 878 G4double zh = G4RandomRadiusInRing(zheig << 961 rRand1 = RandFlat::shoot(0.,1.) ; 879 G4double a = x0; << 962 rRand2 = RandFlat::shoot(0.,1.) ; 880 G4double b = y0; << 963 } while ( rRand2 >= rRand1 ) ; 881 964 882 G4double hh = zheight*zheight; << 965 // rRand2 = RandFlat::shoot(0.,std::sqrt(1.-sqr(rRand1))); 883 G4double aa = a*a; << 966 return G4ThreeVector(rRand1*xSemiAxis*(zheight+zTopCut)*cosphi, 884 G4double bb = b*b; << 967 rRand1*ySemiAxis*(zheight+zTopCut)*sinphi, -zTopCut); 885 G4double R = std::max(a,b); << 886 G4double mu_max = R*std::sqrt(hh + R*R); << 887 968 888 G4double x,y; << 889 for (auto i=0; i<1000; ++i) << 890 { << 891 G4double phi = CLHEP::twopi*G4UniformRand(); << 892 x = std::cos(phi); << 893 y = std::sin(phi); << 894 G4double xx = x*x; << 895 G4double yy = y*y; << 896 G4double E = hh + aa*xx + bb*yy; << 897 G4double F = (aa-bb)*x*y; << 898 G4double G = aa*yy + bb*xx; << 899 G4double mu = std::sqrt(E*G - F*F); << 900 if (mu_max*G4UniformRand() <= mu) brea << 901 } << 902 p.set(zh*xSemiAxis*x,zh*ySemiAxis*y,zhei << 903 break; << 904 } << 905 case 2: // base at +Z, uniform distributio << 906 { << 907 G4double zh = zheight - zTopCut; << 908 G4TwoVector rho = G4RandomPointInEllipse << 909 p.set(rho.x(),rho.y(),zTopCut); << 910 break; << 911 } << 912 } 969 } 913 return p; << 970 // else 914 } << 971 // 915 << 916 ////////////////////////////////////////////// << 917 // << 918 // Get cubic volume << 919 972 920 G4double G4EllipticalCone::GetCubicVolume() << 973 do 921 { << 922 if (fCubicVolume == 0.0) << 923 { 974 { 924 G4double x0 = xSemiAxis*zheight; // x semi << 975 rRand1 = RandFlat::shoot(0.,1.) ; 925 G4double y0 = ySemiAxis*zheight; // y semi << 976 rRand2 = RandFlat::shoot(0.,1.) ; 926 G4double v0 = CLHEP::pi*x0*y0*zheight/3.; << 977 } while ( rRand2 >= rRand1 ) ; 927 G4double kmin = (zTopCut >= zheight ) ? 0. << 928 G4double kmax = (zTopCut >= zheight ) ? 2. << 929 fCubicVolume = (kmax - kmin)*(kmax*kmax + << 930 } << 931 return fCubicVolume; << 932 } << 933 << 934 ////////////////////////////////////////////// << 935 // << 936 // Get surface area << 937 978 938 G4double G4EllipticalCone::GetSurfaceArea() << 979 return G4ThreeVector(rRand1*xSemiAxis*(zheight-zTopCut)*cosphi, 939 { << 980 rRand1*ySemiAxis*(zheight-zTopCut)*sinphi, zTopCut); 940 if (fSurfaceArea == 0.0) << 941 { << 942 G4double x0 = xSemiAxis*zheight; // x semi << 943 G4double y0 = ySemiAxis*zheight; // y semi << 944 G4double s0 = G4GeomTools::EllipticConeLat << 945 G4double kmin = (zTopCut >= zheight ) ? 0. << 946 G4double kmax = (zTopCut >= zheight ) ? 2. << 947 fSurfaceArea = (kmax - kmin)*(kmax + kmin) << 948 + CLHEP::pi*x0*y0*(kmin*kmin << 949 } << 950 return fSurfaceArea; << 951 } 981 } 952 982 953 ////////////////////////////////////////////// << 954 // 983 // 955 // Methods for visualisation 984 // Methods for visualisation >> 985 // 956 986 957 void G4EllipticalCone::DescribeYourselfTo (G4V 987 void G4EllipticalCone::DescribeYourselfTo (G4VGraphicsScene& scene) const 958 { 988 { 959 scene.AddSolid(*this); 989 scene.AddSolid(*this); 960 } 990 } 961 991 962 G4VisExtent G4EllipticalCone::GetExtent() cons 992 G4VisExtent G4EllipticalCone::GetExtent() const 963 { 993 { 964 // Define the sides of the box into which th 994 // Define the sides of the box into which the solid instance would fit. 965 // 995 // 966 G4ThreeVector pmin,pmax; << 996 G4double maxDim; 967 BoundingLimits(pmin,pmax); << 997 maxDim = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis; 968 return { pmin.x(), pmax.x(), pmin.y(), pmax. << 998 maxDim = maxDim > zTopCut ? maxDim : zTopCut; >> 999 >> 1000 return G4VisExtent (-maxDim, maxDim, >> 1001 -maxDim, maxDim, >> 1002 -maxDim, maxDim); >> 1003 } >> 1004 >> 1005 G4NURBS* G4EllipticalCone::CreateNURBS () const >> 1006 { >> 1007 // Box for now!!! >> 1008 // >> 1009 return new G4NURBSbox(xSemiAxis, ySemiAxis,zheight); 969 } 1010 } 970 1011 971 G4Polyhedron* G4EllipticalCone::CreatePolyhedr 1012 G4Polyhedron* G4EllipticalCone::CreatePolyhedron () const 972 { 1013 { 973 return new G4PolyhedronEllipticalCone(xSemiA 1014 return new G4PolyhedronEllipticalCone(xSemiAxis, ySemiAxis, zheight, zTopCut); 974 } 1015 } 975 1016 976 G4Polyhedron* G4EllipticalCone::GetPolyhedron 1017 G4Polyhedron* G4EllipticalCone::GetPolyhedron () const 977 { 1018 { 978 if ( (fpPolyhedron == nullptr) << 1019 if ( (!fpPolyhedron) 979 || fRebuildPolyhedron << 980 || (fpPolyhedron->GetNumberOfRotationSteps 1020 || (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() != 981 fpPolyhedron->GetNumberOfRotationSteps 1021 fpPolyhedron->GetNumberOfRotationSteps()) ) 982 { 1022 { 983 G4AutoLock l(&polyhedronMutex); << 984 delete fpPolyhedron; 1023 delete fpPolyhedron; 985 fpPolyhedron = CreatePolyhedron(); 1024 fpPolyhedron = CreatePolyhedron(); 986 fRebuildPolyhedron = false; << 987 l.unlock(); << 988 } 1025 } 989 return fpPolyhedron; 1026 return fpPolyhedron; 990 } 1027 } 991 << 992 #endif // !defined(G4GEOM_USE_UELLIPTICALCONE) << 993 1028