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