| 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 // class G4Ellipsoid 26 // class G4Ellipsoid
27 // 27 //
28 // Implementation of G4Ellipsoid class 28 // Implementation of G4Ellipsoid class
29 // 29 //
30 // 10.11.99 G.Horton-Smith: first writing, bas 30 // 10.11.99 G.Horton-Smith: first writing, based on G4Sphere class
31 // 25.02.05 G.Guerrieri: Revised 31 // 25.02.05 G.Guerrieri: Revised
32 // 15.12.19 E.Tcherniaev: Complete revision 32 // 15.12.19 E.Tcherniaev: Complete revision
33 // ------------------------------------------- 33 // --------------------------------------------------------------------
34 34
35 #include "G4Ellipsoid.hh" 35 #include "G4Ellipsoid.hh"
36 36
37 #if !(defined(G4GEOM_USE_UELLIPSOID) && define 37 #if !(defined(G4GEOM_USE_UELLIPSOID) && defined(G4GEOM_USE_SYS_USOLIDS))
38 38
39 #include "globals.hh" 39 #include "globals.hh"
40 40
41 #include "G4VoxelLimits.hh" 41 #include "G4VoxelLimits.hh"
42 #include "G4AffineTransform.hh" 42 #include "G4AffineTransform.hh"
43 #include "G4GeometryTolerance.hh" 43 #include "G4GeometryTolerance.hh"
44 #include "G4BoundingEnvelope.hh" 44 #include "G4BoundingEnvelope.hh"
45 #include "G4RandomTools.hh" 45 #include "G4RandomTools.hh"
46 #include "G4QuickRand.hh" 46 #include "G4QuickRand.hh"
47 47
48 #include "G4VPVParameterisation.hh" 48 #include "G4VPVParameterisation.hh"
49 49
50 #include "G4VGraphicsScene.hh" 50 #include "G4VGraphicsScene.hh"
51 #include "G4VisExtent.hh" 51 #include "G4VisExtent.hh"
52 52
53 #include "G4AutoLock.hh" 53 #include "G4AutoLock.hh"
54 54
55 namespace 55 namespace
56 { 56 {
57 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZ 57 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
58 G4Mutex lateralareaMutex = G4MUTEX_INITIALIZ 58 G4Mutex lateralareaMutex = G4MUTEX_INITIALIZER;
59 } 59 }
60 60
61 using namespace CLHEP; 61 using namespace CLHEP;
62 62
63 ////////////////////////////////////////////// 63 //////////////////////////////////////////////////////////////////////////
64 // 64 //
65 // Constructor 65 // Constructor
66 66
67 G4Ellipsoid::G4Ellipsoid(const G4String& name, 67 G4Ellipsoid::G4Ellipsoid(const G4String& name,
68 G4double xSemiA 68 G4double xSemiAxis,
69 G4double ySemiA 69 G4double ySemiAxis,
70 G4double zSemiA 70 G4double zSemiAxis,
71 G4double zBotto 71 G4double zBottomCut,
72 G4double zTopCu 72 G4double zTopCut)
73 : G4VSolid(name), fDx(xSemiAxis), fDy(ySemiA 73 : G4VSolid(name), fDx(xSemiAxis), fDy(ySemiAxis), fDz(zSemiAxis),
74 fZBottomCut(zBottomCut), fZTopCut(zTopCut) 74 fZBottomCut(zBottomCut), fZTopCut(zTopCut)
75 { 75 {
76 CheckParameters(); 76 CheckParameters();
77 } 77 }
78 78
79 ////////////////////////////////////////////// 79 //////////////////////////////////////////////////////////////////////////
80 // 80 //
81 // Fake default constructor - sets only member 81 // Fake default constructor - sets only member data and allocates memory
82 // for usage restri 82 // for usage restricted to object persistency
83 83
84 G4Ellipsoid::G4Ellipsoid( __void__& a ) 84 G4Ellipsoid::G4Ellipsoid( __void__& a )
85 : G4VSolid(a), fDx(0.), fDy(0.), fDz(0.), fZ 85 : G4VSolid(a), fDx(0.), fDy(0.), fDz(0.), fZBottomCut(0.), fZTopCut(0.)
86 { 86 {
87 } 87 }
88 88
89 ////////////////////////////////////////////// 89 //////////////////////////////////////////////////////////////////////////
90 // 90 //
91 // Destructor 91 // Destructor
92 92
93 G4Ellipsoid::~G4Ellipsoid() 93 G4Ellipsoid::~G4Ellipsoid()
94 { 94 {
95 delete fpPolyhedron; fpPolyhedron = nullptr; 95 delete fpPolyhedron; fpPolyhedron = nullptr;
96 } 96 }
97 97
98 ////////////////////////////////////////////// 98 //////////////////////////////////////////////////////////////////////////
99 // 99 //
100 // Copy constructor 100 // Copy constructor
101 101
102 G4Ellipsoid::G4Ellipsoid(const G4Ellipsoid& rh 102 G4Ellipsoid::G4Ellipsoid(const G4Ellipsoid& rhs)
103 : G4VSolid(rhs), 103 : G4VSolid(rhs),
104 fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), 104 fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),
105 fZBottomCut(rhs.fZBottomCut), fZTopCut(rhs. 105 fZBottomCut(rhs.fZBottomCut), fZTopCut(rhs.fZTopCut),
106 halfTolerance(rhs.halfTolerance), 106 halfTolerance(rhs.halfTolerance),
107 fXmax(rhs.fXmax), fYmax(rhs.fYmax), 107 fXmax(rhs.fXmax), fYmax(rhs.fYmax),
108 fRsph(rhs.fRsph), fR(rhs.fR), 108 fRsph(rhs.fRsph), fR(rhs.fR),
109 fSx(rhs.fSx), fSy(rhs.fSy), fSz(rhs.fSz), 109 fSx(rhs.fSx), fSy(rhs.fSy), fSz(rhs.fSz),
110 fZMidCut(rhs.fZMidCut), fZDimCut(rhs.fZDimC 110 fZMidCut(rhs.fZMidCut), fZDimCut(rhs.fZDimCut),
111 fQ1(rhs.fQ1), fQ2(rhs.fQ2), 111 fQ1(rhs.fQ1), fQ2(rhs.fQ2),
112 fCubicVolume(rhs.fCubicVolume), 112 fCubicVolume(rhs.fCubicVolume),
113 fSurfaceArea(rhs.fSurfaceArea), 113 fSurfaceArea(rhs.fSurfaceArea),
114 fLateralArea(rhs.fLateralArea) 114 fLateralArea(rhs.fLateralArea)
115 { 115 {
116 } 116 }
117 117
118 ////////////////////////////////////////////// 118 //////////////////////////////////////////////////////////////////////////
119 // 119 //
120 // Assignment operator 120 // Assignment operator
121 121
122 G4Ellipsoid& G4Ellipsoid::operator = (const G4 122 G4Ellipsoid& G4Ellipsoid::operator = (const G4Ellipsoid& rhs)
123 { 123 {
124 // Check assignment to self 124 // Check assignment to self
125 // 125 //
126 if (this == &rhs) { return *this; } 126 if (this == &rhs) { return *this; }
127 127
128 // Copy base class data 128 // Copy base class data
129 // 129 //
130 G4VSolid::operator=(rhs); 130 G4VSolid::operator=(rhs);
131 131
132 // Copy data 132 // Copy data
133 // 133 //
134 fDx = rhs.fDx; 134 fDx = rhs.fDx;
135 fDy = rhs.fDy; 135 fDy = rhs.fDy;
136 fDz = rhs.fDz; 136 fDz = rhs.fDz;
137 fZBottomCut = rhs.fZBottomCut; 137 fZBottomCut = rhs.fZBottomCut;
138 fZTopCut = rhs.fZTopCut; 138 fZTopCut = rhs.fZTopCut;
139 139
140 halfTolerance = rhs.halfTolerance; 140 halfTolerance = rhs.halfTolerance;
141 fXmax = rhs.fXmax; 141 fXmax = rhs.fXmax;
142 fYmax = rhs.fYmax; 142 fYmax = rhs.fYmax;
143 fRsph = rhs.fRsph; 143 fRsph = rhs.fRsph;
144 fR = rhs.fR; 144 fR = rhs.fR;
145 fSx = rhs.fSx; 145 fSx = rhs.fSx;
146 fSy = rhs.fSy; 146 fSy = rhs.fSy;
147 fSz = rhs.fSz; 147 fSz = rhs.fSz;
148 fZMidCut = rhs.fZMidCut; 148 fZMidCut = rhs.fZMidCut;
149 fZDimCut = rhs.fZDimCut; 149 fZDimCut = rhs.fZDimCut;
150 fQ1 = rhs.fQ1; 150 fQ1 = rhs.fQ1;
151 fQ2 = rhs.fQ2; 151 fQ2 = rhs.fQ2;
152 152
153 fCubicVolume = rhs.fCubicVolume; 153 fCubicVolume = rhs.fCubicVolume;
154 fSurfaceArea = rhs.fSurfaceArea; 154 fSurfaceArea = rhs.fSurfaceArea;
155 fLateralArea = rhs.fLateralArea; 155 fLateralArea = rhs.fLateralArea;
156 156
157 fRebuildPolyhedron = false; 157 fRebuildPolyhedron = false;
158 delete fpPolyhedron; fpPolyhedron = nullptr 158 delete fpPolyhedron; fpPolyhedron = nullptr;
159 159
160 return *this; 160 return *this;
161 } 161 }
162 162
163 ////////////////////////////////////////////// 163 //////////////////////////////////////////////////////////////////////////
164 // 164 //
165 // Check parameters and make precalculation 165 // Check parameters and make precalculation
166 166
167 void G4Ellipsoid::CheckParameters() 167 void G4Ellipsoid::CheckParameters()
168 { 168 {
169 halfTolerance = 0.5 * kCarTolerance; // half 169 halfTolerance = 0.5 * kCarTolerance; // half tolerance
170 G4double dmin = 2. * kCarTolerance; 170 G4double dmin = 2. * kCarTolerance;
171 171
172 // Check dimensions 172 // Check dimensions
173 // 173 //
174 if (fDx < dmin || fDy < dmin || fDz < dmin) 174 if (fDx < dmin || fDy < dmin || fDz < dmin)
175 { 175 {
176 std::ostringstream message; 176 std::ostringstream message;
177 message << "Invalid (too small or negative 177 message << "Invalid (too small or negative) dimensions for Solid: "
178 << GetName() << "\n" 178 << GetName() << "\n"
179 << " semi-axis x: " << fDx << "\n 179 << " semi-axis x: " << fDx << "\n"
180 << " semi-axis y: " << fDy << "\n 180 << " semi-axis y: " << fDy << "\n"
181 << " semi-axis z: " << fDz; 181 << " semi-axis z: " << fDz;
182 G4Exception("G4Ellipsoid::CheckParameters( 182 G4Exception("G4Ellipsoid::CheckParameters()", "GeomSolids0002",
183 FatalException, message); 183 FatalException, message);
184 } 184 }
185 G4double A = fDx; 185 G4double A = fDx;
186 G4double B = fDy; 186 G4double B = fDy;
187 G4double C = fDz; 187 G4double C = fDz;
188 188
189 // Check cuts 189 // Check cuts
190 // 190 //
191 if (fZBottomCut == 0. && fZTopCut == 0.) 191 if (fZBottomCut == 0. && fZTopCut == 0.)
192 { 192 {
193 fZBottomCut = -C; 193 fZBottomCut = -C;
194 fZTopCut = C; 194 fZTopCut = C;
195 } 195 }
196 if (fZBottomCut >= C || fZTopCut <= -C || fZ 196 if (fZBottomCut >= C || fZTopCut <= -C || fZBottomCut >= fZTopCut)
197 { 197 {
198 std::ostringstream message; 198 std::ostringstream message;
199 message << "Invalid Z cuts for Solid: " 199 message << "Invalid Z cuts for Solid: "
200 << GetName() << "\n" 200 << GetName() << "\n"
201 << " bottom cut: " << fZBottomCut 201 << " bottom cut: " << fZBottomCut << "\n"
202 << " top cut: " << fZTopCut; 202 << " top cut: " << fZTopCut;
203 G4Exception("G4Ellipsoid::CheckParameters( 203 G4Exception("G4Ellipsoid::CheckParameters()", "GeomSolids0002",
204 FatalException, message); 204 FatalException, message);
205 205
206 } 206 }
207 fZBottomCut = std::max(fZBottomCut, -C); 207 fZBottomCut = std::max(fZBottomCut, -C);
208 fZTopCut = std::min(fZTopCut, C); 208 fZTopCut = std::min(fZTopCut, C);
209 209
210 // Set extent in x and y 210 // Set extent in x and y
211 fXmax = A; 211 fXmax = A;
212 fYmax = B; 212 fYmax = B;
213 if (fZBottomCut > 0.) 213 if (fZBottomCut > 0.)
214 { 214 {
215 G4double ratio = fZBottomCut / C; 215 G4double ratio = fZBottomCut / C;
216 G4double scale = std::sqrt((1. - ratio) * 216 G4double scale = std::sqrt((1. - ratio) * (1 + ratio));
217 fXmax *= scale; 217 fXmax *= scale;
218 fYmax *= scale; 218 fYmax *= scale;
219 } 219 }
220 if (fZTopCut < 0.) 220 if (fZTopCut < 0.)
221 { 221 {
222 G4double ratio = fZTopCut / C; 222 G4double ratio = fZTopCut / C;
223 G4double scale = std::sqrt((1. - ratio) * 223 G4double scale = std::sqrt((1. - ratio) * (1 + ratio));
224 fXmax *= scale; 224 fXmax *= scale;
225 fYmax *= scale; 225 fYmax *= scale;
226 } 226 }
227 227
228 // Set scale factors 228 // Set scale factors
229 fRsph = std::max(std::max(A, B), C); // boun 229 fRsph = std::max(std::max(A, B), C); // bounding sphere
230 fR = std::min(std::min(A, B), C); // radi 230 fR = std::min(std::min(A, B), C); // radius of sphere after scaling
231 fSx = fR / A; // X scale factor 231 fSx = fR / A; // X scale factor
232 fSy = fR / B; // Y scale factor 232 fSy = fR / B; // Y scale factor
233 fSz = fR / C; // Z scale factor 233 fSz = fR / C; // Z scale factor
234 234
235 // Scaled cuts 235 // Scaled cuts
236 fZMidCut = 0.5 * (fZTopCut + fZBottomCut) * 236 fZMidCut = 0.5 * (fZTopCut + fZBottomCut) * fSz; // middle position
237 fZDimCut = 0.5 * (fZTopCut - fZBottomCut) * 237 fZDimCut = 0.5 * (fZTopCut - fZBottomCut) * fSz; // half distance
238 238
239 // Coefficients for approximation of distanc 239 // Coefficients for approximation of distance: Q1 * (x^2 + y^2 + z^2) - Q2
240 fQ1 = 0.5 / fR; 240 fQ1 = 0.5 / fR;
241 fQ2 = 0.5 * fR + halfTolerance * halfToleran 241 fQ2 = 0.5 * fR + halfTolerance * halfTolerance * fQ1;
242 242
243 fCubicVolume = 0.; // volume 243 fCubicVolume = 0.; // volume
244 fSurfaceArea = 0.; // surface area 244 fSurfaceArea = 0.; // surface area
245 fLateralArea = 0.; // lateral surface area 245 fLateralArea = 0.; // lateral surface area
246 } 246 }
247 247
248 ////////////////////////////////////////////// 248 //////////////////////////////////////////////////////////////////////////
249 // 249 //
250 // Dispatch to parameterisation for replicatio 250 // Dispatch to parameterisation for replication mechanism dimension
251 // computation & modification 251 // computation & modification
252 252
253 void G4Ellipsoid::ComputeDimensions(G4VPVParam 253 void G4Ellipsoid::ComputeDimensions(G4VPVParameterisation* p,
254 const G4in 254 const G4int n,
255 const G4VP 255 const G4VPhysicalVolume* pRep)
256 { 256 {
257 p->ComputeDimensions(*this,n,pRep); 257 p->ComputeDimensions(*this,n,pRep);
258 } 258 }
259 259
260 ////////////////////////////////////////////// 260 //////////////////////////////////////////////////////////////////////////
261 // 261 //
262 // Get bounding box 262 // Get bounding box
263 263
264 void G4Ellipsoid::BoundingLimits(G4ThreeVector 264 void G4Ellipsoid::BoundingLimits(G4ThreeVector& pMin,
265 G4ThreeVector 265 G4ThreeVector& pMax) const
266 { 266 {
267 pMin.set(-fXmax,-fYmax, fZBottomCut); 267 pMin.set(-fXmax,-fYmax, fZBottomCut);
268 pMax.set( fXmax, fYmax, fZTopCut); 268 pMax.set( fXmax, fYmax, fZTopCut);
269 } 269 }
270 270
271 ////////////////////////////////////////////// 271 //////////////////////////////////////////////////////////////////////////
272 // 272 //
273 // Calculate extent under transform and specif 273 // Calculate extent under transform and specified limits
274 274
275 G4bool 275 G4bool
276 G4Ellipsoid::CalculateExtent(const EAxis pAxis 276 G4Ellipsoid::CalculateExtent(const EAxis pAxis,
277 const G4VoxelLimi 277 const G4VoxelLimits& pVoxelLimit,
278 const G4AffineTra 278 const G4AffineTransform& pTransform,
279 G4double& p 279 G4double& pMin, G4double& pMax) const
280 { 280 {
281 G4ThreeVector bmin, bmax; 281 G4ThreeVector bmin, bmax;
282 282
283 // Get bounding box 283 // Get bounding box
284 BoundingLimits(bmin,bmax); 284 BoundingLimits(bmin,bmax);
285 285
286 // Find extent 286 // Find extent
287 G4BoundingEnvelope bbox(bmin,bmax); 287 G4BoundingEnvelope bbox(bmin,bmax);
288 return bbox.CalculateExtent(pAxis,pVoxelLimi 288 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
289 } 289 }
290 290
291 ////////////////////////////////////////////// 291 //////////////////////////////////////////////////////////////////////////
292 // 292 //
293 // Return position of point: inside/outside/on 293 // Return position of point: inside/outside/on surface
294 294
295 EInside G4Ellipsoid::Inside(const G4ThreeVecto 295 EInside G4Ellipsoid::Inside(const G4ThreeVector& p) const
296 { 296 {
297 G4double x = p.x() * fSx; 297 G4double x = p.x() * fSx;
298 G4double y = p.y() * fSy; 298 G4double y = p.y() * fSy;
299 G4double z = p.z() * fSz; 299 G4double z = p.z() * fSz;
300 G4double rr = x * x + y * y + z * z; 300 G4double rr = x * x + y * y + z * z;
301 G4double distZ = std::abs(z - fZMidCut) - fZ 301 G4double distZ = std::abs(z - fZMidCut) - fZDimCut;
302 G4double distR = fQ1 * rr - fQ2; 302 G4double distR = fQ1 * rr - fQ2;
303 G4double dist = std::max(distZ, distR); 303 G4double dist = std::max(distZ, distR);
304 304
305 if (dist > halfTolerance) return kOutside; 305 if (dist > halfTolerance) return kOutside;
306 return (dist > -halfTolerance) ? kSurface : 306 return (dist > -halfTolerance) ? kSurface : kInside;
307 } 307 }
308 308
309 ////////////////////////////////////////////// 309 //////////////////////////////////////////////////////////////////////////
310 // 310 //
311 // Return unit normal to surface at p 311 // Return unit normal to surface at p
312 312
313 G4ThreeVector G4Ellipsoid::SurfaceNormal( cons 313 G4ThreeVector G4Ellipsoid::SurfaceNormal( const G4ThreeVector& p) const
314 { 314 {
315 G4ThreeVector norm(0., 0., 0.); 315 G4ThreeVector norm(0., 0., 0.);
316 G4int nsurf = 0; 316 G4int nsurf = 0;
317 317
318 // Check cuts 318 // Check cuts
319 G4double x = p.x() * fSx; 319 G4double x = p.x() * fSx;
320 G4double y = p.y() * fSy; 320 G4double y = p.y() * fSy;
321 G4double z = p.z() * fSz; 321 G4double z = p.z() * fSz;
322 G4double distZ = std::abs(z - fZMidCut) - fZ 322 G4double distZ = std::abs(z - fZMidCut) - fZDimCut;
323 if (std::abs(distZ) <= halfTolerance) 323 if (std::abs(distZ) <= halfTolerance)
324 { 324 {
325 norm.setZ(std::copysign(1., z - fZMidCut)) 325 norm.setZ(std::copysign(1., z - fZMidCut));
326 ++nsurf; 326 ++nsurf;
327 } 327 }
328 328
329 // Check lateral surface 329 // Check lateral surface
330 G4double distR = fQ1*(x*x + y*y + z*z) - fQ2 330 G4double distR = fQ1*(x*x + y*y + z*z) - fQ2;
331 if (std::abs(distR) <= halfTolerance) 331 if (std::abs(distR) <= halfTolerance)
332 { 332 {
333 // normal = (p.x/a^2, p.y/b^2, p.z/c^2) 333 // normal = (p.x/a^2, p.y/b^2, p.z/c^2)
334 norm += G4ThreeVector(x*fSx, y*fSy, z*fSz) 334 norm += G4ThreeVector(x*fSx, y*fSy, z*fSz).unit();
335 ++nsurf; 335 ++nsurf;
336 } 336 }
337 337
338 // Return normal 338 // Return normal
339 if (nsurf == 1) return norm; 339 if (nsurf == 1) return norm;
340 else if (nsurf > 1) return norm.unit(); // e 340 else if (nsurf > 1) return norm.unit(); // edge
341 else 341 else
342 { 342 {
343 #ifdef G4SPECSDEBUG 343 #ifdef G4SPECSDEBUG
344 std::ostringstream message; 344 std::ostringstream message;
345 G4long oldprc = message.precision(16); << 345 G4int oldprc = message.precision(16);
346 message << "Point p is not on surface (!?) 346 message << "Point p is not on surface (!?) of solid: "
347 << GetName() << "\n"; 347 << GetName() << "\n";
348 message << "Position:\n"; 348 message << "Position:\n";
349 message << " p.x() = " << p.x()/mm << " 349 message << " p.x() = " << p.x()/mm << " mm\n";
350 message << " p.y() = " << p.y()/mm << " 350 message << " p.y() = " << p.y()/mm << " mm\n";
351 message << " p.z() = " << p.z()/mm << " 351 message << " p.z() = " << p.z()/mm << " mm";
352 G4cout.precision(oldprc); 352 G4cout.precision(oldprc);
353 G4Exception("G4Ellipsoid::SurfaceNormal(p) 353 G4Exception("G4Ellipsoid::SurfaceNormal(p)", "GeomSolids1002",
354 JustWarning, message ); 354 JustWarning, message );
355 DumpInfo(); 355 DumpInfo();
356 #endif 356 #endif
357 return ApproxSurfaceNormal(p); 357 return ApproxSurfaceNormal(p);
358 } 358 }
359 } 359 }
360 360
361 ////////////////////////////////////////////// 361 //////////////////////////////////////////////////////////////////////////
362 // 362 //
363 // Find surface nearest to point and return co 363 // Find surface nearest to point and return corresponding normal.
364 // This method normally should not be called. 364 // This method normally should not be called.
365 365
366 G4ThreeVector G4Ellipsoid::ApproxSurfaceNormal 366 G4ThreeVector G4Ellipsoid::ApproxSurfaceNormal(const G4ThreeVector& p) const
367 { 367 {
368 G4double x = p.x() * fSx; 368 G4double x = p.x() * fSx;
369 G4double y = p.y() * fSy; 369 G4double y = p.y() * fSy;
370 G4double z = p.z() * fSz; 370 G4double z = p.z() * fSz;
371 G4double rr = x * x + y * y + z * z; 371 G4double rr = x * x + y * y + z * z;
372 G4double distZ = std::abs(z - fZMidCut) - fZ 372 G4double distZ = std::abs(z - fZMidCut) - fZDimCut;
373 G4double distR = std::sqrt(rr) - fR; 373 G4double distR = std::sqrt(rr) - fR;
374 if (distR > distZ && rr > 0.) // distR > di 374 if (distR > distZ && rr > 0.) // distR > distZ is correct!
375 return G4ThreeVector(x*fSx, y*fSy, z*fSz). 375 return G4ThreeVector(x*fSx, y*fSy, z*fSz).unit();
376 else 376 else
377 return { 0., 0., std::copysign(1., z - fZM << 377 return G4ThreeVector(0., 0., std::copysign(1., z - fZMidCut));
378 } 378 }
379 379
380 ////////////////////////////////////////////// 380 //////////////////////////////////////////////////////////////////////////
381 // 381 //
382 // Calculate distance to shape from outside al 382 // Calculate distance to shape from outside along normalised vector
383 383
384 G4double G4Ellipsoid::DistanceToIn(const G4Thr 384 G4double G4Ellipsoid::DistanceToIn(const G4ThreeVector& p,
385 const G4Thr 385 const G4ThreeVector& v) const
386 { 386 {
387 G4double offset = 0.; 387 G4double offset = 0.;
388 G4ThreeVector pcur = p; 388 G4ThreeVector pcur = p;
389 389
390 // Check if point is flying away, relative t 390 // Check if point is flying away, relative to bounding box
391 // 391 //
392 G4double safex = std::abs(p.x()) - fXmax; 392 G4double safex = std::abs(p.x()) - fXmax;
393 G4double safey = std::abs(p.y()) - fYmax; 393 G4double safey = std::abs(p.y()) - fYmax;
394 G4double safet = p.z() - fZTopCut; 394 G4double safet = p.z() - fZTopCut;
395 G4double safeb = fZBottomCut - p.z(); 395 G4double safeb = fZBottomCut - p.z();
396 396
397 if (safex >= -halfTolerance && p.x() * v.x() 397 if (safex >= -halfTolerance && p.x() * v.x() >= 0.) return kInfinity;
398 if (safey >= -halfTolerance && p.y() * v.y() 398 if (safey >= -halfTolerance && p.y() * v.y() >= 0.) return kInfinity;
399 if (safet >= -halfTolerance && v.z() >= 0.) 399 if (safet >= -halfTolerance && v.z() >= 0.) return kInfinity;
400 if (safeb >= -halfTolerance && v.z() <= 0.) 400 if (safeb >= -halfTolerance && v.z() <= 0.) return kInfinity;
401 401
402 // Relocate point, if required 402 // Relocate point, if required
403 // 403 //
404 G4double safe = std::max(std::max(std::max(s 404 G4double safe = std::max(std::max(std::max(safex, safey), safet), safeb);
405 if (safe > 32. * fRsph) 405 if (safe > 32. * fRsph)
406 { 406 {
407 offset = (1. - 1.e-08) * safe - 2. * fRsph 407 offset = (1. - 1.e-08) * safe - 2. * fRsph;
408 pcur += offset * v; 408 pcur += offset * v;
409 G4double dist = DistanceToIn(pcur, v); 409 G4double dist = DistanceToIn(pcur, v);
410 return (dist == kInfinity) ? kInfinity : d 410 return (dist == kInfinity) ? kInfinity : dist + offset;
411 } 411 }
412 412
413 // Scale ellipsoid to sphere 413 // Scale ellipsoid to sphere
414 // 414 //
415 G4double px = pcur.x() * fSx; 415 G4double px = pcur.x() * fSx;
416 G4double py = pcur.y() * fSy; 416 G4double py = pcur.y() * fSy;
417 G4double pz = pcur.z() * fSz; 417 G4double pz = pcur.z() * fSz;
418 G4double vx = v.x() * fSx; 418 G4double vx = v.x() * fSx;
419 G4double vy = v.y() * fSy; 419 G4double vy = v.y() * fSy;
420 G4double vz = v.z() * fSz; 420 G4double vz = v.z() * fSz;
421 421
422 // Check if point is leaving the solid 422 // Check if point is leaving the solid
423 // 423 //
424 G4double dzcut = fZDimCut; 424 G4double dzcut = fZDimCut;
425 G4double pzcut = pz - fZMidCut; 425 G4double pzcut = pz - fZMidCut;
426 G4double distZ = std::abs(pzcut) - dzcut; 426 G4double distZ = std::abs(pzcut) - dzcut;
427 if (distZ >= -halfTolerance && pzcut * vz >= 427 if (distZ >= -halfTolerance && pzcut * vz >= 0.) return kInfinity;
428 428
429 G4double rr = px * px + py * py + pz * pz; 429 G4double rr = px * px + py * py + pz * pz;
430 G4double pv = px * vx + py * vy + pz * vz; 430 G4double pv = px * vx + py * vy + pz * vz;
431 G4double distR = fQ1 * rr - fQ2; 431 G4double distR = fQ1 * rr - fQ2;
432 if (distR >= -halfTolerance && pv >= 0.) ret 432 if (distR >= -halfTolerance && pv >= 0.) return kInfinity;
433 433
434 G4double A = vx * vx + vy * vy + vz * vz; 434 G4double A = vx * vx + vy * vy + vz * vz;
435 G4double B = pv; 435 G4double B = pv;
436 G4double C = rr - fR * fR; 436 G4double C = rr - fR * fR;
437 G4double D = B * B - A * C; 437 G4double D = B * B - A * C;
438 // scratch^2 = R^2 - (R - halfTolerance)^2 = 438 // scratch^2 = R^2 - (R - halfTolerance)^2 = 2 * R * halfTolerance
439 G4double EPS = A * A * fR * kCarTolerance; / 439 G4double EPS = A * A * fR * kCarTolerance; // discriminant at scratching
440 if (D <= EPS) return kInfinity; // no inters 440 if (D <= EPS) return kInfinity; // no intersection or scratching
441 441
442 // Find intersection with Z planes 442 // Find intersection with Z planes
443 // 443 //
444 G4double invz = (vz == 0) ? DBL_MAX : -1./v 444 G4double invz = (vz == 0) ? DBL_MAX : -1./vz;
445 G4double dz = std::copysign(dzcut, invz); 445 G4double dz = std::copysign(dzcut, invz);
446 G4double tzmin = (pzcut - dz) * invz; 446 G4double tzmin = (pzcut - dz) * invz;
447 G4double tzmax = (pzcut + dz) * invz; 447 G4double tzmax = (pzcut + dz) * invz;
448 448
449 // Find intersection with lateral surface 449 // Find intersection with lateral surface
450 // 450 //
451 G4double tmp = -B - std::copysign(std::sqrt( 451 G4double tmp = -B - std::copysign(std::sqrt(D), B);
452 G4double t1 = tmp / A; 452 G4double t1 = tmp / A;
453 G4double t2 = C / tmp; 453 G4double t2 = C / tmp;
454 G4double trmin = std::min(t1, t2); 454 G4double trmin = std::min(t1, t2);
455 G4double trmax = std::max(t1, t2); 455 G4double trmax = std::max(t1, t2);
456 456
457 // Return distance 457 // Return distance
458 // 458 //
459 G4double tmin = std::max(tzmin, trmin); 459 G4double tmin = std::max(tzmin, trmin);
460 G4double tmax = std::min(tzmax, trmax); 460 G4double tmax = std::min(tzmax, trmax);
461 461
462 if (tmax - tmin <= halfTolerance) return kIn 462 if (tmax - tmin <= halfTolerance) return kInfinity; // touch or no hit
463 return (tmin < halfTolerance) ? offset : tmi 463 return (tmin < halfTolerance) ? offset : tmin + offset;
464 } 464 }
465 465
466 ////////////////////////////////////////////// 466 //////////////////////////////////////////////////////////////////////////
467 // 467 //
468 // Estimate distance to surface from outside 468 // Estimate distance to surface from outside
469 469
470 G4double G4Ellipsoid::DistanceToIn(const G4Thr 470 G4double G4Ellipsoid::DistanceToIn(const G4ThreeVector& p) const
471 { 471 {
472 G4double px = p.x(); 472 G4double px = p.x();
473 G4double py = p.y(); 473 G4double py = p.y();
474 G4double pz = p.z(); 474 G4double pz = p.z();
475 475
476 // Safety distance to bounding box 476 // Safety distance to bounding box
477 G4double distX = std::abs(px) - fXmax; 477 G4double distX = std::abs(px) - fXmax;
478 G4double distY = std::abs(py) - fYmax; 478 G4double distY = std::abs(py) - fYmax;
479 G4double distZ = std::max(pz - fZTopCut, fZB 479 G4double distZ = std::max(pz - fZTopCut, fZBottomCut - pz);
480 G4double distB = std::max(std::max(distX, di 480 G4double distB = std::max(std::max(distX, distY), distZ);
481 481
482 // Safety distance to lateral surface 482 // Safety distance to lateral surface
483 G4double x = px * fSx; 483 G4double x = px * fSx;
484 G4double y = py * fSy; 484 G4double y = py * fSy;
485 G4double z = pz * fSz; 485 G4double z = pz * fSz;
486 G4double distR = std::sqrt(x*x + y*y + z*z) 486 G4double distR = std::sqrt(x*x + y*y + z*z) - fR;
487 487
488 // Return safety to in 488 // Return safety to in
489 G4double dist = std::max(distB, distR); 489 G4double dist = std::max(distB, distR);
490 return (dist < 0.) ? 0. : dist; 490 return (dist < 0.) ? 0. : dist;
491 } 491 }
492 492
493 ////////////////////////////////////////////// 493 //////////////////////////////////////////////////////////////////////////
494 // 494 //
495 // Calculate distance to surface from inside a 495 // Calculate distance to surface from inside along normalised vector
496 496
497 G4double G4Ellipsoid::DistanceToOut(const G4Th 497 G4double G4Ellipsoid::DistanceToOut(const G4ThreeVector& p,
498 const G4Th 498 const G4ThreeVector& v,
499 const G4bo 499 const G4bool calcNorm,
500 G4bo 500 G4bool* validNorm,
501 G4Th 501 G4ThreeVector* n ) const
502 { 502 {
503 // Check if point flying away relative to Z 503 // Check if point flying away relative to Z planes
504 // 504 //
505 G4double pz = p.z() * fSz; 505 G4double pz = p.z() * fSz;
506 G4double vz = v.z() * fSz; 506 G4double vz = v.z() * fSz;
507 G4double dzcut = fZDimCut; 507 G4double dzcut = fZDimCut;
508 G4double pzcut = pz - fZMidCut; 508 G4double pzcut = pz - fZMidCut;
509 G4double distZ = std::abs(pzcut) - dzcut; 509 G4double distZ = std::abs(pzcut) - dzcut;
510 if (distZ >= -halfTolerance && pzcut * vz > 510 if (distZ >= -halfTolerance && pzcut * vz > 0.)
511 { 511 {
512 if (calcNorm) 512 if (calcNorm)
513 { 513 {
514 *validNorm = true; 514 *validNorm = true;
515 n->set(0., 0., std::copysign(1., pzcut)) 515 n->set(0., 0., std::copysign(1., pzcut));
516 } 516 }
517 return 0.; 517 return 0.;
518 } 518 }
519 519
520 // Check if point is flying away relative to 520 // Check if point is flying away relative to lateral surface
521 // 521 //
522 G4double px = p.x() * fSx; 522 G4double px = p.x() * fSx;
523 G4double py = p.y() * fSy; 523 G4double py = p.y() * fSy;
524 G4double vx = v.x() * fSx; 524 G4double vx = v.x() * fSx;
525 G4double vy = v.y() * fSy; 525 G4double vy = v.y() * fSy;
526 G4double rr = px * px + py * py + pz * pz; 526 G4double rr = px * px + py * py + pz * pz;
527 G4double pv = px * vx + py * vy + pz * vz; 527 G4double pv = px * vx + py * vy + pz * vz;
528 G4double distR = fQ1 * rr - fQ2; 528 G4double distR = fQ1 * rr - fQ2;
529 if (distR >= -halfTolerance && pv > 0.) 529 if (distR >= -halfTolerance && pv > 0.)
530 { 530 {
531 if (calcNorm) 531 if (calcNorm)
532 { 532 {
533 *validNorm = true; 533 *validNorm = true;
534 *n = G4ThreeVector(px*fSx, py*fSy, pz*fS 534 *n = G4ThreeVector(px*fSx, py*fSy, pz*fSz).unit();
535 } 535 }
536 return 0.; 536 return 0.;
537 } 537 }
538 538
539 // Just in case check if point is outside (n 539 // Just in case check if point is outside (normally it should never be)
540 // 540 //
541 if (std::max(distZ, distR) > halfTolerance) 541 if (std::max(distZ, distR) > halfTolerance)
542 { 542 {
543 #ifdef G4SPECSDEBUG 543 #ifdef G4SPECSDEBUG
544 std::ostringstream message; 544 std::ostringstream message;
545 G4long oldprc = message.precision(16); << 545 G4int oldprc = message.precision(16);
546 message << "Point p is outside (!?) of sol 546 message << "Point p is outside (!?) of solid: "
547 << GetName() << G4endl; 547 << GetName() << G4endl;
548 message << "Position: " << p << G4endl;; 548 message << "Position: " << p << G4endl;;
549 message << "Direction: " << v; 549 message << "Direction: " << v;
550 G4cout.precision(oldprc); 550 G4cout.precision(oldprc);
551 G4Exception("G4Ellipsoid::DistanceToOut(p, 551 G4Exception("G4Ellipsoid::DistanceToOut(p,v)", "GeomSolids1002",
552 JustWarning, message ); 552 JustWarning, message );
553 DumpInfo(); 553 DumpInfo();
554 #endif 554 #endif
555 if (calcNorm) 555 if (calcNorm)
556 { 556 {
557 *validNorm = true; 557 *validNorm = true;
558 *n = ApproxSurfaceNormal(p); 558 *n = ApproxSurfaceNormal(p);
559 } 559 }
560 return 0.; 560 return 0.;
561 } 561 }
562 562
563 // Set coefficients of quadratic equation: A 563 // Set coefficients of quadratic equation: A t^2 + 2B t + C = 0
564 // 564 //
565 G4double A = vx * vx + vy * vy + vz * vz; 565 G4double A = vx * vx + vy * vy + vz * vz;
566 G4double B = pv; 566 G4double B = pv;
567 G4double C = rr - fR * fR; 567 G4double C = rr - fR * fR;
568 G4double D = B * B - A * C; 568 G4double D = B * B - A * C;
569 // It is expected that the point is located 569 // It is expected that the point is located inside the sphere, so
570 // max term in the expression for discrimina 570 // max term in the expression for discriminant is A * R^2 and
571 // max calculation error can be derived as f 571 // max calculation error can be derived as follows:
572 // A * (1 + 2e) * R^2 * (1 + 2e) = A * R^2 + 572 // A * (1 + 2e) * R^2 * (1 + 2e) = A * R^2 + (4 * A * R^2 * e)
573 G4double EPS = 4. * A * fR * fR * DBL_EPSILO 573 G4double EPS = 4. * A * fR * fR * DBL_EPSILON; // calculation error
574 574
575 if (D <= EPS) // no intersection 575 if (D <= EPS) // no intersection
576 { 576 {
577 if (calcNorm) 577 if (calcNorm)
578 { 578 {
579 *validNorm = true; 579 *validNorm = true;
580 *n = G4ThreeVector(px*fSx, py*fSy, pz*fS 580 *n = G4ThreeVector(px*fSx, py*fSy, pz*fSz).unit();
581 } 581 }
582 return 0.; 582 return 0.;
583 } 583 }
584 584
585 // Find intersection with Z cuts 585 // Find intersection with Z cuts
586 // 586 //
587 G4double tzmax = (vz == 0.) ? DBL_MAX : (std 587 G4double tzmax = (vz == 0.) ? DBL_MAX : (std::copysign(dzcut, vz) - pzcut) / vz;
588 588
589 // Find intersection with lateral surface 589 // Find intersection with lateral surface
590 // 590 //
591 G4double tmp = -B - std::copysign(std::sqrt( 591 G4double tmp = -B - std::copysign(std::sqrt(D), B);
592 G4double trmax = (tmp < 0.) ? C/tmp : tmp/A; 592 G4double trmax = (tmp < 0.) ? C/tmp : tmp/A;
593 593
594 // Find distance and set normal, if required 594 // Find distance and set normal, if required
595 // 595 //
596 G4double tmax = std::min(tzmax, trmax); 596 G4double tmax = std::min(tzmax, trmax);
597 //if (tmax < halfTolerance) tmax = 0.; 597 //if (tmax < halfTolerance) tmax = 0.;
598 598
599 if (calcNorm) 599 if (calcNorm)
600 { 600 {
601 *validNorm = true; 601 *validNorm = true;
602 if (tmax == tzmax) 602 if (tmax == tzmax)
603 { 603 {
604 G4double pznew = pz + tmax * vz; 604 G4double pznew = pz + tmax * vz;
605 n->set(0., 0., (pznew > fZMidCut) ? 1. : 605 n->set(0., 0., (pznew > fZMidCut) ? 1. : -1.);
606 } 606 }
607 else 607 else
608 { 608 {
609 G4double nx = (px + tmax * vx) * fSx; 609 G4double nx = (px + tmax * vx) * fSx;
610 G4double ny = (py + tmax * vy) * fSy; 610 G4double ny = (py + tmax * vy) * fSy;
611 G4double nz = (pz + tmax * vz) * fSz; 611 G4double nz = (pz + tmax * vz) * fSz;
612 *n = G4ThreeVector(nx, ny, nz).unit(); 612 *n = G4ThreeVector(nx, ny, nz).unit();
613 } 613 }
614 } 614 }
615 return tmax; 615 return tmax;
616 } 616 }
617 617
618 ////////////////////////////////////////////// 618 //////////////////////////////////////////////////////////////////////////
619 // 619 //
620 // Estimate distance to surface from inside 620 // Estimate distance to surface from inside
621 621
622 G4double G4Ellipsoid::DistanceToOut(const G4Th 622 G4double G4Ellipsoid::DistanceToOut(const G4ThreeVector& p) const
623 { 623 {
624 // Safety distance in z direction 624 // Safety distance in z direction
625 G4double distZ = std::min(fZTopCut - p.z(), 625 G4double distZ = std::min(fZTopCut - p.z(), p.z() - fZBottomCut);
626 626
627 // Safety distance to lateral surface 627 // Safety distance to lateral surface
628 G4double x = p.x() * fSx; 628 G4double x = p.x() * fSx;
629 G4double y = p.y() * fSy; 629 G4double y = p.y() * fSy;
630 G4double z = p.z() * fSz; 630 G4double z = p.z() * fSz;
631 G4double distR = fR - std::sqrt(x*x + y*y + 631 G4double distR = fR - std::sqrt(x*x + y*y + z*z);
632 632
633 // Return safety to out 633 // Return safety to out
634 G4double dist = std::min(distZ, distR); 634 G4double dist = std::min(distZ, distR);
635 return (dist < 0.) ? 0. : dist; 635 return (dist < 0.) ? 0. : dist;
636 } 636 }
637 637
638 ////////////////////////////////////////////// 638 //////////////////////////////////////////////////////////////////////////
639 // 639 //
640 // Return entity type 640 // Return entity type
641 641
642 G4GeometryType G4Ellipsoid::GetEntityType() co 642 G4GeometryType G4Ellipsoid::GetEntityType() const
643 { 643 {
644 return {"G4Ellipsoid"}; << 644 return G4String("G4Ellipsoid");
645 } 645 }
646 646
647 ////////////////////////////////////////////// 647 //////////////////////////////////////////////////////////////////////////
648 // 648 //
649 // Make a clone of the object 649 // Make a clone of the object
650 650
651 G4VSolid* G4Ellipsoid::Clone() const 651 G4VSolid* G4Ellipsoid::Clone() const
652 { 652 {
653 return new G4Ellipsoid(*this); 653 return new G4Ellipsoid(*this);
654 } 654 }
655 655
656 ////////////////////////////////////////////// 656 //////////////////////////////////////////////////////////////////////////
657 // 657 //
658 // Stream object contents to output stream 658 // Stream object contents to output stream
659 659
660 std::ostream& G4Ellipsoid::StreamInfo( std::os 660 std::ostream& G4Ellipsoid::StreamInfo( std::ostream& os ) const
661 { 661 {
662 G4long oldprc = os.precision(16); << 662 G4int oldprc = os.precision(16);
663 os << "------------------------------------- 663 os << "-----------------------------------------------------------\n"
664 << " *** Dump for solid - " << GetName 664 << " *** Dump for solid - " << GetName() << " ***\n"
665 << " ================================= 665 << " ===================================================\n"
666 << " Solid type: " << GetEntityType() << 666 << " Solid type: " << GetEntityType() << "\n"
667 << " Parameters: \n" 667 << " Parameters: \n"
668 << " semi-axis x: " << GetDx()/mm << " 668 << " semi-axis x: " << GetDx()/mm << " mm \n"
669 << " semi-axis y: " << GetDy()/mm << " 669 << " semi-axis y: " << GetDy()/mm << " mm \n"
670 << " semi-axis z: " << GetDz()/mm << " 670 << " semi-axis z: " << GetDz()/mm << " mm \n"
671 << " lower cut in z: " << GetZBottomCu 671 << " lower cut in z: " << GetZBottomCut()/mm << " mm \n"
672 << " upper cut in z: " << GetZTopCut() 672 << " upper cut in z: " << GetZTopCut()/mm << " mm \n"
673 << "------------------------------------- 673 << "-----------------------------------------------------------\n";
674 os.precision(oldprc); 674 os.precision(oldprc);
675 return os; 675 return os;
676 } 676 }
677 677
678 ////////////////////////////////////////////// 678 //////////////////////////////////////////////////////////////////////////
679 // 679 //
680 // Return volume 680 // Return volume
681 681
682 G4double G4Ellipsoid::GetCubicVolume() 682 G4double G4Ellipsoid::GetCubicVolume()
683 { 683 {
684 if (fCubicVolume == 0.) 684 if (fCubicVolume == 0.)
685 { 685 {
686 G4double piAB_3 = CLHEP::pi * fDx * fDy / 686 G4double piAB_3 = CLHEP::pi * fDx * fDy / 3.;
687 fCubicVolume = 4. * piAB_3 * fDz; 687 fCubicVolume = 4. * piAB_3 * fDz;
688 if (fZBottomCut > -fDz) 688 if (fZBottomCut > -fDz)
689 { 689 {
690 G4double hbot = 1. + fZBottomCut / fDz; 690 G4double hbot = 1. + fZBottomCut / fDz;
691 fCubicVolume -= piAB_3 * hbot * hbot * ( 691 fCubicVolume -= piAB_3 * hbot * hbot * (2. * fDz - fZBottomCut);
692 } 692 }
693 if (fZTopCut < fDz) 693 if (fZTopCut < fDz)
694 { 694 {
695 G4double htop = 1. - fZTopCut / fDz; 695 G4double htop = 1. - fZTopCut / fDz;
696 fCubicVolume -= piAB_3 * htop * htop * ( 696 fCubicVolume -= piAB_3 * htop * htop * (2. * fDz + fZTopCut);
697 } 697 }
698 } 698 }
699 return fCubicVolume; 699 return fCubicVolume;
700 } 700 }
701 701
702 ////////////////////////////////////////////// 702 //////////////////////////////////////////////////////////////////////////
703 // 703 //
704 // Calculate area of lateral surface 704 // Calculate area of lateral surface
705 705
706 G4double G4Ellipsoid::LateralSurfaceArea() con 706 G4double G4Ellipsoid::LateralSurfaceArea() const
707 { 707 {
708 constexpr G4int NPHI = 1000.; << 708 const G4int Nphi = 100;
709 constexpr G4double dPhi = CLHEP::halfpi/NPHI << 709 const G4int Nz = 200;
710 constexpr G4double eps = 4.*DBL_EPSILON; << 710 G4double rho[Nz + 1];
711 << 711
712 G4double aa = fDx*fDx; << 712 // Set array of rho
713 G4double bb = fDy*fDy; << 713 G4double zbot = fZBottomCut / fDz;
714 G4double cc = fDz*fDz; << 714 G4double ztop = fZTopCut / fDz;
715 G4double ab = fDx*fDy; << 715 G4double dz = (ztop - zbot) / Nz;
716 G4double cc_aa = cc/aa; << 716 for (G4int iz = 0; iz < Nz; ++iz)
717 G4double cc_bb = cc/bb; << 717 {
718 G4double zmax = std::min(fZTopCut, fDz); << 718 G4double z = zbot + iz * dz;
719 G4double zmin = std::max(fZBottomCut,-fDz); << 719 rho[iz] = std::sqrt((1. + z) * (1. - z));
720 G4double zmax_c = zmax/fDz; << 720 }
721 G4double zmin_c = zmin/fDz; << 721 rho[Nz] = std::sqrt((1. + ztop) * (1. - ztop));
>> 722
>> 723 // Compute area
>> 724 zbot = fZBottomCut;
>> 725 ztop = fZTopCut;
>> 726 dz = (ztop - zbot) / Nz;
722 G4double area = 0.; 727 G4double area = 0.;
723 << 728 G4double dphi = CLHEP::halfpi / Nphi;
724 if (aa == bb) // spheroid, use analytical ex << 729 for (G4int iphi = 0; iphi < Nphi; ++iphi)
725 { 730 {
726 G4double k = fDz/fDx; << 731 G4double phi1 = iphi * dphi;
727 G4double kk = k*k; << 732 G4double phi2 = (iphi == Nphi - 1) ? CLHEP::halfpi : phi1 + dphi;
728 if (kk < 1. - eps) << 733 G4double cos1 = std::cos(phi1) * fDx;
729 { << 734 G4double cos2 = std::cos(phi2) * fDx;
730 G4double invk = fDx/fDz; << 735 G4double sin1 = std::sin(phi1) * fDy;
731 G4double root = std::sqrt(1. - kk); << 736 G4double sin2 = std::sin(phi2) * fDy;
732 G4double tmax = zmax_c*root; << 737 for (G4int iz = 0; iz < Nz; ++iz)
733 G4double tmin = zmin_c*root; << 738 {
734 area = CLHEP::pi*ab* << 739 G4double z1 = zbot + iz * dz;
735 ((zmax_c*std::sqrt(kk + tmax*tmax) - z << 740 G4double z2 = (iz == Nz - 1) ? ztop : z1 + dz;
736 (std::asinh(tmax*invk) - std::asinh(t << 741 G4double rho1 = rho[iz];
737 } << 742 G4double rho2 = rho[iz + 1];
738 else if (kk > 1. + eps) << 743 G4ThreeVector p1(rho1 * cos1, rho1 * sin1, z1);
739 { << 744 G4ThreeVector p2(rho1 * cos2, rho1 * sin2, z1);
740 G4double invk = fDx/fDz; << 745 G4ThreeVector p3(rho2 * cos1, rho2 * sin1, z2);
741 G4double root = std::sqrt(kk - 1.); << 746 G4ThreeVector p4(rho2 * cos2, rho2 * sin2, z2);
742 G4double tmax = zmax_c*root; << 747 area += ((p4 - p1).cross(p3 - p2)).mag();
743 G4double tmin = zmin_c*root; <<
744 area = CLHEP::pi*ab* <<
745 ((zmax_c*std::sqrt(kk - tmax*tmax) - z <<
746 (std::asin(tmax*invk) - std::asin(tmi <<
747 } <<
748 else <<
749 { <<
750 area = CLHEP::twopi*fDx*(zmax - zmin); <<
751 } <<
752 return area; <<
753 } <<
754 <<
755 // ellipsoid, integration along phi <<
756 for (G4int i = 0; i < NPHI; ++i) <<
757 { <<
758 G4double sinPhi = std::sin(dPhi*(i + 0.5)) <<
759 G4double kk = cc_aa + (cc_bb - cc_aa)*sinP <<
760 if (kk < 1. - eps) <<
761 { <<
762 G4double root = std::sqrt(1. - kk); <<
763 G4double tmax = zmax_c*root; <<
764 G4double tmin = zmin_c*root; <<
765 G4double invk = 1./std::sqrt(kk); <<
766 area += 2.*ab*dPhi* <<
767 ((zmax_c*std::sqrt(kk + tmax*tmax) - z <<
768 (std::asinh(tmax*invk) - std::asinh(t <<
769 } <<
770 else if (kk > 1. + eps) <<
771 { <<
772 G4double root = std::sqrt(kk - 1.); <<
773 G4double tmax = zmax_c*root; <<
774 G4double tmin = zmin_c*root; <<
775 G4double invk = 1./std::sqrt(kk); <<
776 area += 2.*ab*dPhi* <<
777 ((zmax_c*std::sqrt(kk - tmax*tmax) - z <<
778 (std::asin(tmax*invk) - std::asin(tmi <<
779 } <<
780 else <<
781 { <<
782 area += 4.*ab*dPhi*(zmax_c - zmin_c); <<
783 } 748 }
784 } 749 }
785 return area; << 750 return 2. * area;
786 } 751 }
787 752
788 ////////////////////////////////////////////// 753 //////////////////////////////////////////////////////////////////////////
789 // 754 //
790 // Return surface area 755 // Return surface area
791 756
792 G4double G4Ellipsoid::GetSurfaceArea() 757 G4double G4Ellipsoid::GetSurfaceArea()
793 { 758 {
794 if (fSurfaceArea == 0.) 759 if (fSurfaceArea == 0.)
795 { 760 {
796 G4double piAB = CLHEP::pi * fDx * fDy; 761 G4double piAB = CLHEP::pi * fDx * fDy;
797 fSurfaceArea = LateralSurfaceArea(); 762 fSurfaceArea = LateralSurfaceArea();
798 if (fZBottomCut > -fDz) 763 if (fZBottomCut > -fDz)
799 { 764 {
800 G4double hbot = 1. + fZBottomCut / fDz; 765 G4double hbot = 1. + fZBottomCut / fDz;
801 fSurfaceArea += piAB * hbot * (2. - hbot 766 fSurfaceArea += piAB * hbot * (2. - hbot);
802 } 767 }
803 if (fZTopCut < fDz) 768 if (fZTopCut < fDz)
804 { 769 {
805 G4double htop = 1. - fZTopCut / fDz; 770 G4double htop = 1. - fZTopCut / fDz;
806 fSurfaceArea += piAB * htop * (2. - htop 771 fSurfaceArea += piAB * htop * (2. - htop);
807 } 772 }
808 } 773 }
809 return fSurfaceArea; 774 return fSurfaceArea;
810 } 775 }
811 776
812 ////////////////////////////////////////////// 777 //////////////////////////////////////////////////////////////////////////
813 // 778 //
814 // Return random point on surface 779 // Return random point on surface
815 780
816 G4ThreeVector G4Ellipsoid::GetPointOnSurface() 781 G4ThreeVector G4Ellipsoid::GetPointOnSurface() const
817 { 782 {
818 G4double A = GetDx(); 783 G4double A = GetDx();
819 G4double B = GetDy(); 784 G4double B = GetDy();
820 G4double C = GetDz(); 785 G4double C = GetDz();
821 G4double Zbot = GetZBottomCut(); 786 G4double Zbot = GetZBottomCut();
822 G4double Ztop = GetZTopCut(); 787 G4double Ztop = GetZTopCut();
823 788
824 // Calculate cut areas 789 // Calculate cut areas
825 G4double Hbot = 1. + Zbot / C; 790 G4double Hbot = 1. + Zbot / C;
826 G4double Htop = 1. - Ztop / C; 791 G4double Htop = 1. - Ztop / C;
827 G4double piAB = CLHEP::pi * A * B; 792 G4double piAB = CLHEP::pi * A * B;
828 G4double Sbot = piAB * Hbot * (2. - Hbot); 793 G4double Sbot = piAB * Hbot * (2. - Hbot);
829 G4double Stop = piAB * Htop * (2. - Htop); 794 G4double Stop = piAB * Htop * (2. - Htop);
830 795
831 // Get area of lateral surface 796 // Get area of lateral surface
832 if (fLateralArea == 0.) 797 if (fLateralArea == 0.)
833 { 798 {
834 G4AutoLock l(&lateralareaMutex); 799 G4AutoLock l(&lateralareaMutex);
835 fLateralArea = LateralSurfaceArea(); 800 fLateralArea = LateralSurfaceArea();
836 l.unlock(); 801 l.unlock();
837 } 802 }
838 G4double Slat = fLateralArea; 803 G4double Slat = fLateralArea;
839 804
840 // Select surface (0 - bottom cut, 1 - later 805 // Select surface (0 - bottom cut, 1 - lateral surface, 2 - top cut)
841 G4double select = (Sbot + Slat + Stop) * G4Q 806 G4double select = (Sbot + Slat + Stop) * G4QuickRand();
842 G4int k = 0; 807 G4int k = 0;
843 if (select > Sbot) k = 1; 808 if (select > Sbot) k = 1;
844 if (select > Sbot + Slat) k = 2; 809 if (select > Sbot + Slat) k = 2;
845 810
846 // Pick random point on selected surface (re 811 // Pick random point on selected surface (rejection sampling)
847 G4ThreeVector p; 812 G4ThreeVector p;
848 switch (k) 813 switch (k)
849 { 814 {
850 case 0: // bootom z-cut 815 case 0: // bootom z-cut
851 { 816 {
852 G4double scale = std::sqrt(Hbot * (2. - 817 G4double scale = std::sqrt(Hbot * (2. - Hbot));
853 G4TwoVector rho = G4RandomPointInEllipse 818 G4TwoVector rho = G4RandomPointInEllipse(A * scale, B * scale);
854 p.set(rho.x(), rho.y(), Zbot); 819 p.set(rho.x(), rho.y(), Zbot);
855 break; 820 break;
856 } 821 }
857 case 1: // lateral surface 822 case 1: // lateral surface
858 { 823 {
859 G4double x, y, z; 824 G4double x, y, z;
860 G4double mu_max = std::max(std::max(A * 825 G4double mu_max = std::max(std::max(A * B, A * C), B * C);
861 for (G4int i = 0; i < 1000; ++i) 826 for (G4int i = 0; i < 1000; ++i)
862 { 827 {
863 // generate random point on unit spher 828 // generate random point on unit sphere
864 z = (Zbot + (Ztop - Zbot) * G4QuickRan 829 z = (Zbot + (Ztop - Zbot) * G4QuickRand()) / C;
865 G4double rho = std::sqrt((1. + z) * (1 830 G4double rho = std::sqrt((1. + z) * (1. - z));
866 G4double phi = CLHEP::twopi * G4QuickR 831 G4double phi = CLHEP::twopi * G4QuickRand();
867 x = rho * std::cos(phi); 832 x = rho * std::cos(phi);
868 y = rho * std::sin(phi); 833 y = rho * std::sin(phi);
869 // check acceptance 834 // check acceptance
870 G4double xbc = x * B * C; 835 G4double xbc = x * B * C;
871 G4double yac = y * A * C; 836 G4double yac = y * A * C;
872 G4double zab = z * A * B; 837 G4double zab = z * A * B;
873 G4double mu = std::sqrt(xbc * xbc + y 838 G4double mu = std::sqrt(xbc * xbc + yac * yac + zab * zab);
874 if (mu_max * G4QuickRand() <= mu) brea 839 if (mu_max * G4QuickRand() <= mu) break;
875 } 840 }
876 p.set(A * x, B * y, C * z); 841 p.set(A * x, B * y, C * z);
877 break; 842 break;
878 } 843 }
879 case 2: // top z-cut 844 case 2: // top z-cut
880 { 845 {
881 G4double scale = std::sqrt(Htop * (2. - 846 G4double scale = std::sqrt(Htop * (2. - Htop));
882 G4TwoVector rho = G4RandomPointInEllipse 847 G4TwoVector rho = G4RandomPointInEllipse(A * scale, B * scale);
883 p.set(rho.x(), rho.y(), Ztop); 848 p.set(rho.x(), rho.y(), Ztop);
884 break; 849 break;
885 } 850 }
886 } 851 }
887 return p; 852 return p;
888 } 853 }
889 854
890 ////////////////////////////////////////////// 855 //////////////////////////////////////////////////////////////////////////
891 // 856 //
892 // Methods for visualisation 857 // Methods for visualisation
893 858
894 void G4Ellipsoid::DescribeYourselfTo (G4VGraph 859 void G4Ellipsoid::DescribeYourselfTo (G4VGraphicsScene& scene) const
895 { 860 {
896 scene.AddSolid(*this); 861 scene.AddSolid(*this);
897 } 862 }
898 863
899 ////////////////////////////////////////////// 864 //////////////////////////////////////////////////////////////////////////
900 // 865 //
901 // Return vis extent 866 // Return vis extent
902 867
903 G4VisExtent G4Ellipsoid::GetExtent() const 868 G4VisExtent G4Ellipsoid::GetExtent() const
904 { 869 {
905 return { -fXmax, fXmax, -fYmax, fYmax, fZBot << 870 return G4VisExtent(-fXmax, fXmax, -fYmax, fYmax, fZBottomCut, fZTopCut);
906 } 871 }
907 872
908 ////////////////////////////////////////////// 873 //////////////////////////////////////////////////////////////////////////
909 // 874 //
910 // Create polyhedron 875 // Create polyhedron
911 876
912 G4Polyhedron* G4Ellipsoid::CreatePolyhedron () 877 G4Polyhedron* G4Ellipsoid::CreatePolyhedron () const
913 { 878 {
914 return new G4PolyhedronEllipsoid(fDx, fDy, f 879 return new G4PolyhedronEllipsoid(fDx, fDy, fDz, fZBottomCut, fZTopCut);
915 } 880 }
916 881
917 ////////////////////////////////////////////// 882 //////////////////////////////////////////////////////////////////////////
918 // 883 //
919 // Return pointer to polyhedron 884 // Return pointer to polyhedron
920 885
921 G4Polyhedron* G4Ellipsoid::GetPolyhedron () co 886 G4Polyhedron* G4Ellipsoid::GetPolyhedron () const
922 { 887 {
923 if (fpPolyhedron == nullptr || 888 if (fpPolyhedron == nullptr ||
924 fRebuildPolyhedron || 889 fRebuildPolyhedron ||
925 fpPolyhedron->GetNumberOfRotationStepsAt 890 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
926 fpPolyhedron->GetNumberOfRotationSteps() 891 fpPolyhedron->GetNumberOfRotationSteps())
927 { 892 {
928 G4AutoLock l(&polyhedronMutex); 893 G4AutoLock l(&polyhedronMutex);
929 delete fpPolyhedron; 894 delete fpPolyhedron;
930 fpPolyhedron = CreatePolyhedron(); 895 fpPolyhedron = CreatePolyhedron();
931 fRebuildPolyhedron = false; 896 fRebuildPolyhedron = false;
932 l.unlock(); 897 l.unlock();
933 } 898 }
934 return fpPolyhedron; 899 return fpPolyhedron;
935 } 900 }
936 901
937 #endif // !defined(G4GEOM_USE_UELLIPSOID) || ! 902 #endif // !defined(G4GEOM_USE_UELLIPSOID) || !defined(G4GEOM_USE_SYS_USOLIDS)
938 903