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