1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // $Id: G4Ellipsoid.cc,v 1.14 2007/05/18 07:39:56 gcosmo Exp $ |
---|
27 | // GEANT4 tag $Name: geant4-09-02-ref-02 $ |
---|
28 | // |
---|
29 | // class G4Ellipsoid |
---|
30 | // |
---|
31 | // Implementation for G4Ellipsoid class |
---|
32 | // |
---|
33 | // History: |
---|
34 | // |
---|
35 | // 10.11.99 G.Horton-Smith -- first writing, based on G4Sphere class |
---|
36 | // 25.02.05 G.Guerrieri -- Modified for future Geant4 release |
---|
37 | // |
---|
38 | // -------------------------------------------------------------------- |
---|
39 | |
---|
40 | #include "globals.hh" |
---|
41 | |
---|
42 | #include "G4Ellipsoid.hh" |
---|
43 | |
---|
44 | #include "G4VoxelLimits.hh" |
---|
45 | #include "G4AffineTransform.hh" |
---|
46 | #include "G4GeometryTolerance.hh" |
---|
47 | |
---|
48 | #include "meshdefs.hh" |
---|
49 | |
---|
50 | #include "Randomize.hh" |
---|
51 | |
---|
52 | #include "G4VGraphicsScene.hh" |
---|
53 | #include "G4Polyhedron.hh" |
---|
54 | #include "G4NURBS.hh" |
---|
55 | #include "G4NURBSbox.hh" |
---|
56 | #include "G4VisExtent.hh" |
---|
57 | |
---|
58 | using namespace CLHEP; |
---|
59 | |
---|
60 | /////////////////////////////////////////////////////////////////////////////// |
---|
61 | // |
---|
62 | // constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI |
---|
63 | // - note if pDPhi>2PI then reset to 2PI |
---|
64 | |
---|
65 | G4Ellipsoid::G4Ellipsoid(const G4String& pName, |
---|
66 | G4double pxSemiAxis, |
---|
67 | G4double pySemiAxis, |
---|
68 | G4double pzSemiAxis, |
---|
69 | G4double pzBottomCut, |
---|
70 | G4double pzTopCut) |
---|
71 | : G4VSolid(pName), fpPolyhedron(0), fCubicVolume(0.), fSurfaceArea(0.), |
---|
72 | zBottomCut(0.), zTopCut(0.) |
---|
73 | { |
---|
74 | // note: for users that want to use the full ellipsoid it is useful |
---|
75 | // to include a default for the cuts |
---|
76 | |
---|
77 | kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance(); |
---|
78 | |
---|
79 | // Check Semi-Axis |
---|
80 | if ( (pxSemiAxis>0.) && (pySemiAxis>0.) && (pzSemiAxis>0.) ) |
---|
81 | { |
---|
82 | SetSemiAxis(pxSemiAxis, pySemiAxis, pzSemiAxis); |
---|
83 | } |
---|
84 | else |
---|
85 | { |
---|
86 | G4cerr << "ERROR - G4Ellipsoid::G4Ellipsoid(): " << GetName() << G4endl |
---|
87 | << " Invalid semi-axis !" |
---|
88 | << G4endl; |
---|
89 | G4Exception("G4Ellipsoid::G4Ellipsoid()", "InvalidSetup", |
---|
90 | FatalException, "Invalid semi-axis."); |
---|
91 | } |
---|
92 | |
---|
93 | if ( pzBottomCut == 0 && pzTopCut == 0 ) |
---|
94 | { |
---|
95 | SetZCuts(-pzSemiAxis, pzSemiAxis); |
---|
96 | } |
---|
97 | else if ( (pzBottomCut < pzSemiAxis) && (pzTopCut > -pzSemiAxis) |
---|
98 | && (pzBottomCut < pzTopCut) ) |
---|
99 | { |
---|
100 | SetZCuts(pzBottomCut, pzTopCut); |
---|
101 | } |
---|
102 | else |
---|
103 | { |
---|
104 | G4cerr << "ERROR - G4Ellipsoid::G4Ellipsoid(): " << GetName() << G4endl |
---|
105 | << " Invalid z-coordinate for cutting plane !" |
---|
106 | << G4endl; |
---|
107 | G4Exception("G4Ellipsoid::G4Ellipsoid()", "InvalidSetup", |
---|
108 | FatalException, "Invalid z-coordinate for cutting plane."); |
---|
109 | } |
---|
110 | } |
---|
111 | |
---|
112 | /////////////////////////////////////////////////////////////////////////////// |
---|
113 | // |
---|
114 | // Fake default constructor - sets only member data and allocates memory |
---|
115 | // for usage restricted to object persistency. |
---|
116 | // |
---|
117 | G4Ellipsoid::G4Ellipsoid( __void__& a ) |
---|
118 | : G4VSolid(a), fpPolyhedron(0), fCubicVolume(0.), fSurfaceArea(0.) |
---|
119 | { |
---|
120 | } |
---|
121 | |
---|
122 | /////////////////////////////////////////////////////////////////////////////// |
---|
123 | // |
---|
124 | // Destructor |
---|
125 | |
---|
126 | G4Ellipsoid::~G4Ellipsoid() |
---|
127 | { |
---|
128 | } |
---|
129 | |
---|
130 | /////////////////////////////////////////////////////////////////////////////// |
---|
131 | // |
---|
132 | // Calculate extent under transform and specified limit |
---|
133 | |
---|
134 | G4bool |
---|
135 | G4Ellipsoid::CalculateExtent(const EAxis pAxis, |
---|
136 | const G4VoxelLimits& pVoxelLimit, |
---|
137 | const G4AffineTransform& pTransform, |
---|
138 | G4double& pMin, G4double& pMax) const |
---|
139 | { |
---|
140 | if (!pTransform.IsRotated()) |
---|
141 | { |
---|
142 | // Special case handling for unrotated solid ellipsoid |
---|
143 | // Compute x/y/z mins and maxs for bounding box respecting limits, |
---|
144 | // with early returns if outside limits. Then switch() on pAxis, |
---|
145 | // and compute exact x and y limit for x/y case |
---|
146 | |
---|
147 | G4double xoffset,xMin,xMax; |
---|
148 | G4double yoffset,yMin,yMax; |
---|
149 | G4double zoffset,zMin,zMax; |
---|
150 | |
---|
151 | G4double maxDiff,newMin,newMax; |
---|
152 | G4double xoff,yoff; |
---|
153 | |
---|
154 | xoffset=pTransform.NetTranslation().x(); |
---|
155 | xMin=xoffset - xSemiAxis; |
---|
156 | xMax=xoffset + xSemiAxis; |
---|
157 | if (pVoxelLimit.IsXLimited()) |
---|
158 | { |
---|
159 | if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance) |
---|
160 | || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) ) |
---|
161 | { |
---|
162 | return false; |
---|
163 | } |
---|
164 | else |
---|
165 | { |
---|
166 | if (xMin<pVoxelLimit.GetMinXExtent()) |
---|
167 | { |
---|
168 | xMin=pVoxelLimit.GetMinXExtent(); |
---|
169 | } |
---|
170 | if (xMax>pVoxelLimit.GetMaxXExtent()) |
---|
171 | { |
---|
172 | xMax=pVoxelLimit.GetMaxXExtent(); |
---|
173 | } |
---|
174 | } |
---|
175 | } |
---|
176 | |
---|
177 | yoffset=pTransform.NetTranslation().y(); |
---|
178 | yMin=yoffset - ySemiAxis; |
---|
179 | yMax=yoffset + ySemiAxis; |
---|
180 | if (pVoxelLimit.IsYLimited()) |
---|
181 | { |
---|
182 | if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance) |
---|
183 | || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) ) |
---|
184 | { |
---|
185 | return false; |
---|
186 | } |
---|
187 | else |
---|
188 | { |
---|
189 | if (yMin<pVoxelLimit.GetMinYExtent()) |
---|
190 | { |
---|
191 | yMin=pVoxelLimit.GetMinYExtent(); |
---|
192 | } |
---|
193 | if (yMax>pVoxelLimit.GetMaxYExtent()) |
---|
194 | { |
---|
195 | yMax=pVoxelLimit.GetMaxYExtent(); |
---|
196 | } |
---|
197 | } |
---|
198 | } |
---|
199 | |
---|
200 | zoffset=pTransform.NetTranslation().z(); |
---|
201 | zMin=zoffset + (-zSemiAxis > zBottomCut ? -zSemiAxis : zBottomCut); |
---|
202 | zMax=zoffset + ( zSemiAxis < zTopCut ? zSemiAxis : zTopCut); |
---|
203 | if (pVoxelLimit.IsZLimited()) |
---|
204 | { |
---|
205 | if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance) |
---|
206 | || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) ) |
---|
207 | { |
---|
208 | return false; |
---|
209 | } |
---|
210 | else |
---|
211 | { |
---|
212 | if (zMin<pVoxelLimit.GetMinZExtent()) |
---|
213 | { |
---|
214 | zMin=pVoxelLimit.GetMinZExtent(); |
---|
215 | } |
---|
216 | if (zMax>pVoxelLimit.GetMaxZExtent()) |
---|
217 | { |
---|
218 | zMax=pVoxelLimit.GetMaxZExtent(); |
---|
219 | } |
---|
220 | } |
---|
221 | } |
---|
222 | |
---|
223 | // if here, then known to cut bounding box around ellipsoid |
---|
224 | // |
---|
225 | xoff = (xoffset < xMin) ? (xMin-xoffset) |
---|
226 | : (xoffset > xMax) ? (xoffset-xMax) : 0.0; |
---|
227 | yoff = (yoffset < yMin) ? (yMin-yoffset) |
---|
228 | : (yoffset > yMax) ? (yoffset-yMax) : 0.0; |
---|
229 | |
---|
230 | // detailed calculations |
---|
231 | // NOTE: does not use X or Y offsets to adjust Z range, |
---|
232 | // and does not use Z offset to adjust X or Y range, |
---|
233 | // which is consistent with G4Sphere::CalculateExtent behavior |
---|
234 | // |
---|
235 | switch (pAxis) |
---|
236 | { |
---|
237 | case kXAxis: |
---|
238 | if (yoff==0.) |
---|
239 | { |
---|
240 | // YZ limits cross max/min x => no change |
---|
241 | // |
---|
242 | pMin=xMin; |
---|
243 | pMax=xMax; |
---|
244 | } |
---|
245 | else |
---|
246 | { |
---|
247 | // YZ limits don't cross max/min x => compute max delta x, |
---|
248 | // hence new mins/maxs |
---|
249 | // |
---|
250 | maxDiff= 1.0-sqr(yoff/ySemiAxis); |
---|
251 | if (maxDiff < 0.0) { return false; } |
---|
252 | maxDiff= xSemiAxis * std::sqrt(maxDiff); |
---|
253 | newMin=xoffset-maxDiff; |
---|
254 | newMax=xoffset+maxDiff; |
---|
255 | pMin=(newMin<xMin) ? xMin : newMin; |
---|
256 | pMax=(newMax>xMax) ? xMax : newMax; |
---|
257 | } |
---|
258 | break; |
---|
259 | case kYAxis: |
---|
260 | if (xoff==0.) |
---|
261 | { |
---|
262 | // XZ limits cross max/min y => no change |
---|
263 | // |
---|
264 | pMin=yMin; |
---|
265 | pMax=yMax; |
---|
266 | } |
---|
267 | else |
---|
268 | { |
---|
269 | // XZ limits don't cross max/min y => compute max delta y, |
---|
270 | // hence new mins/maxs |
---|
271 | // |
---|
272 | maxDiff= 1.0-sqr(xoff/xSemiAxis); |
---|
273 | if (maxDiff < 0.0) { return false; } |
---|
274 | maxDiff= ySemiAxis * std::sqrt(maxDiff); |
---|
275 | newMin=yoffset-maxDiff; |
---|
276 | newMax=yoffset+maxDiff; |
---|
277 | pMin=(newMin<yMin) ? yMin : newMin; |
---|
278 | pMax=(newMax>yMax) ? yMax : newMax; |
---|
279 | } |
---|
280 | break; |
---|
281 | case kZAxis: |
---|
282 | pMin=zMin; |
---|
283 | pMax=zMax; |
---|
284 | break; |
---|
285 | default: |
---|
286 | break; |
---|
287 | } |
---|
288 | |
---|
289 | pMin-=kCarTolerance; |
---|
290 | pMax+=kCarTolerance; |
---|
291 | return true; |
---|
292 | } |
---|
293 | else // not rotated |
---|
294 | { |
---|
295 | G4int i,j,noEntries,noBetweenSections; |
---|
296 | G4bool existsAfterClip=false; |
---|
297 | |
---|
298 | // Calculate rotated vertex coordinates |
---|
299 | |
---|
300 | G4int noPolygonVertices=0; |
---|
301 | G4ThreeVectorList* vertices = |
---|
302 | CreateRotatedVertices(pTransform,noPolygonVertices); |
---|
303 | |
---|
304 | pMin=+kInfinity; |
---|
305 | pMax=-kInfinity; |
---|
306 | |
---|
307 | noEntries=vertices->size(); // noPolygonVertices*noPhiCrossSections |
---|
308 | noBetweenSections=noEntries-noPolygonVertices; |
---|
309 | |
---|
310 | G4ThreeVectorList ThetaPolygon; |
---|
311 | for (i=0;i<noEntries;i+=noPolygonVertices) |
---|
312 | { |
---|
313 | for(j=0;j<(noPolygonVertices/2)-1;j++) |
---|
314 | { |
---|
315 | ThetaPolygon.push_back((*vertices)[i+j]); |
---|
316 | ThetaPolygon.push_back((*vertices)[i+j+1]); |
---|
317 | ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]); |
---|
318 | ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]); |
---|
319 | CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax); |
---|
320 | ThetaPolygon.clear(); |
---|
321 | } |
---|
322 | } |
---|
323 | for (i=0;i<noBetweenSections;i+=noPolygonVertices) |
---|
324 | { |
---|
325 | for(j=0;j<noPolygonVertices-1;j++) |
---|
326 | { |
---|
327 | ThetaPolygon.push_back((*vertices)[i+j]); |
---|
328 | ThetaPolygon.push_back((*vertices)[i+j+1]); |
---|
329 | ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]); |
---|
330 | ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]); |
---|
331 | CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax); |
---|
332 | ThetaPolygon.clear(); |
---|
333 | } |
---|
334 | ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]); |
---|
335 | ThetaPolygon.push_back((*vertices)[i]); |
---|
336 | ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]); |
---|
337 | ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]); |
---|
338 | CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax); |
---|
339 | ThetaPolygon.clear(); |
---|
340 | } |
---|
341 | if ( (pMin!=kInfinity) || (pMax!=-kInfinity) ) |
---|
342 | { |
---|
343 | existsAfterClip=true; |
---|
344 | |
---|
345 | // Add 2*tolerance to avoid precision troubles |
---|
346 | // |
---|
347 | pMin-=kCarTolerance; |
---|
348 | pMax+=kCarTolerance; |
---|
349 | |
---|
350 | } |
---|
351 | else |
---|
352 | { |
---|
353 | // Check for case where completely enveloping clipping volume |
---|
354 | // If point inside then we are confident that the solid completely |
---|
355 | // envelopes the clipping volume. Hence set min/max extents according |
---|
356 | // to clipping volume extents along the specified axis. |
---|
357 | // |
---|
358 | G4ThreeVector |
---|
359 | clipCentre((pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5, |
---|
360 | (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5, |
---|
361 | (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5); |
---|
362 | |
---|
363 | if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside) |
---|
364 | { |
---|
365 | existsAfterClip=true; |
---|
366 | pMin=pVoxelLimit.GetMinExtent(pAxis); |
---|
367 | pMax=pVoxelLimit.GetMaxExtent(pAxis); |
---|
368 | } |
---|
369 | } |
---|
370 | delete vertices; |
---|
371 | return existsAfterClip; |
---|
372 | } |
---|
373 | } |
---|
374 | |
---|
375 | /////////////////////////////////////////////////////////////////////////////// |
---|
376 | // |
---|
377 | // Return whether point inside/outside/on surface |
---|
378 | // Split into radius, phi, theta checks |
---|
379 | // Each check modifies `in', or returns as approprate |
---|
380 | |
---|
381 | EInside G4Ellipsoid::Inside(const G4ThreeVector& p) const |
---|
382 | { |
---|
383 | G4double rad2oo, // outside surface outer tolerance |
---|
384 | rad2oi; // outside surface inner tolerance |
---|
385 | EInside in; |
---|
386 | |
---|
387 | // check this side of z cut first, because that's fast |
---|
388 | // |
---|
389 | if (p.z() < zBottomCut-kRadTolerance/2.0) |
---|
390 | { return in=kOutside; } |
---|
391 | if (p.z() > zTopCut+kRadTolerance/2.0) |
---|
392 | { return in=kOutside; } |
---|
393 | |
---|
394 | rad2oo= sqr(p.x()/(xSemiAxis+kRadTolerance/2.)) |
---|
395 | + sqr(p.y()/(ySemiAxis+kRadTolerance/2.)) |
---|
396 | + sqr(p.z()/(zSemiAxis+kRadTolerance/2.)); |
---|
397 | |
---|
398 | if (rad2oo > 1.0) |
---|
399 | { return in=kOutside; } |
---|
400 | |
---|
401 | rad2oi= sqr(p.x()*(1.0+kRadTolerance/2./xSemiAxis)/xSemiAxis) |
---|
402 | + sqr(p.y()*(1.0+kRadTolerance/2./ySemiAxis)/ySemiAxis) |
---|
403 | + sqr(p.z()*(1.0+kRadTolerance/2./zSemiAxis)/zSemiAxis); |
---|
404 | |
---|
405 | // Check radial surfaces |
---|
406 | // sets `in' (already checked for rad2oo > 1.0) |
---|
407 | // |
---|
408 | if (rad2oi < 1.0) |
---|
409 | { |
---|
410 | in = ( (p.z() < zBottomCut+kRadTolerance/2.0) |
---|
411 | || (p.z() > zTopCut-kRadTolerance/2.0) ) ? kSurface : kInside; |
---|
412 | } |
---|
413 | else |
---|
414 | { |
---|
415 | in = kSurface; |
---|
416 | } |
---|
417 | |
---|
418 | return in; |
---|
419 | } |
---|
420 | |
---|
421 | /////////////////////////////////////////////////////////////////////////////// |
---|
422 | // |
---|
423 | // Return unit normal of surface closest to p not protected against p=0 |
---|
424 | |
---|
425 | G4ThreeVector G4Ellipsoid::SurfaceNormal( const G4ThreeVector& p) const |
---|
426 | { |
---|
427 | G4double distR, distZBottom, distZTop; |
---|
428 | |
---|
429 | // normal vector with special magnitude: parallel to normal, units 1/length |
---|
430 | // norm*p == 1.0 if on surface, >1.0 if outside, <1.0 if inside |
---|
431 | // |
---|
432 | G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis), |
---|
433 | p.y()/(ySemiAxis*ySemiAxis), |
---|
434 | p.z()/(zSemiAxis*zSemiAxis)); |
---|
435 | G4double radius = 1.0/norm.mag(); |
---|
436 | |
---|
437 | // approximate distance to curved surface |
---|
438 | // |
---|
439 | distR = std::fabs( (p*norm - 1.0) * radius ) / 2.0; |
---|
440 | |
---|
441 | // Distance to z-cut plane |
---|
442 | // |
---|
443 | distZBottom = std::fabs( p.z() - zBottomCut ); |
---|
444 | distZTop = std::fabs( p.z() - zTopCut ); |
---|
445 | |
---|
446 | if ( (distZBottom < distR) || (distZTop < distR) ) |
---|
447 | { |
---|
448 | return G4ThreeVector(0.,0.,(distZBottom < distZTop) ? -1.0 : 1.0); |
---|
449 | } |
---|
450 | return ( norm *= radius ); |
---|
451 | } |
---|
452 | |
---|
453 | /////////////////////////////////////////////////////////////////////////////// |
---|
454 | // |
---|
455 | // Calculate distance to shape from outside, along normalised vector |
---|
456 | // - return kInfinity if no intersection, or intersection distance <= tolerance |
---|
457 | // |
---|
458 | |
---|
459 | G4double G4Ellipsoid::DistanceToIn( const G4ThreeVector& p, |
---|
460 | const G4ThreeVector& v ) const |
---|
461 | { |
---|
462 | G4double distMin; |
---|
463 | |
---|
464 | distMin= kInfinity; |
---|
465 | |
---|
466 | // check to see if Z plane is relevant |
---|
467 | if (p.z() < zBottomCut) { |
---|
468 | if (v.z() <= 0.0) |
---|
469 | return distMin; |
---|
470 | G4double distZ = (zBottomCut - p.z()) / v.z(); |
---|
471 | if (distZ > kRadTolerance/2.0 && Inside(p+distZ*v) != kOutside ) |
---|
472 | { |
---|
473 | // early exit since can't intercept curved surface if we reach here |
---|
474 | return distMin= distZ; |
---|
475 | } |
---|
476 | } |
---|
477 | if (p.z() > zTopCut) { |
---|
478 | if (v.z() >= 0.0) |
---|
479 | return distMin; |
---|
480 | G4double distZ = (zTopCut - p.z()) / v.z(); |
---|
481 | if (distZ > kRadTolerance/2.0 && Inside(p+distZ*v) != kOutside ) |
---|
482 | { |
---|
483 | // early exit since can't intercept curved surface if we reach here |
---|
484 | return distMin= distZ; |
---|
485 | } |
---|
486 | } |
---|
487 | // if fZCut1 <= p.z() <= fZCut2, then must hit curved surface |
---|
488 | |
---|
489 | // now check curved surface intercept |
---|
490 | G4double A,B,C; |
---|
491 | |
---|
492 | A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis); |
---|
493 | C= sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) + sqr(p.z()/zSemiAxis) - 1.0; |
---|
494 | B= 2.0 * ( p.x()*v.x()/(xSemiAxis*xSemiAxis) + p.y()*v.y()/(ySemiAxis*ySemiAxis) |
---|
495 | + p.z()*v.z()/(zSemiAxis*zSemiAxis) ); |
---|
496 | |
---|
497 | C= B*B - 4.0*A*C; |
---|
498 | if (C > 0.0) |
---|
499 | { |
---|
500 | G4double distR= (-B - std::sqrt(C) ) / (2.0*A); |
---|
501 | G4double intZ= p.z()+distR*v.z(); |
---|
502 | if (distR > kRadTolerance/2.0 |
---|
503 | && intZ >= zBottomCut-kRadTolerance/2.0 |
---|
504 | && intZ <= zTopCut+kRadTolerance/2.0) |
---|
505 | { |
---|
506 | distMin = distR; |
---|
507 | } |
---|
508 | else |
---|
509 | { |
---|
510 | distR= (-B + std::sqrt(C) ) / (2.0*A); |
---|
511 | intZ= p.z()+distR*v.z(); |
---|
512 | if (distR > kRadTolerance/2.0 |
---|
513 | && intZ >= zBottomCut-kRadTolerance/2.0 |
---|
514 | && intZ <= zTopCut+kRadTolerance/2.0) |
---|
515 | { |
---|
516 | distMin = distR; |
---|
517 | } |
---|
518 | } |
---|
519 | } |
---|
520 | |
---|
521 | return distMin; |
---|
522 | } |
---|
523 | |
---|
524 | /////////////////////////////////////////////////////////////////////////////// |
---|
525 | // |
---|
526 | // Calculate distance (<= actual) to closest surface of shape from outside |
---|
527 | // - Return 0 if point inside |
---|
528 | |
---|
529 | G4double G4Ellipsoid::DistanceToIn(const G4ThreeVector& p) const |
---|
530 | { |
---|
531 | G4double distR, distZ; |
---|
532 | |
---|
533 | // normal vector: parallel to normal, magnitude 1/(characteristic radius) |
---|
534 | // |
---|
535 | G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis), |
---|
536 | p.y()/(ySemiAxis*ySemiAxis), |
---|
537 | p.z()/(zSemiAxis*zSemiAxis)); |
---|
538 | G4double radius= 1.0/norm.mag(); |
---|
539 | |
---|
540 | // approximate distance to curved surface ( <= actual distance ) |
---|
541 | // |
---|
542 | distR= (p*norm - 1.0) * radius / 2.0; |
---|
543 | |
---|
544 | // Distance to z-cut plane |
---|
545 | // |
---|
546 | distZ= zBottomCut - p.z(); |
---|
547 | if (distZ < 0.0) |
---|
548 | { |
---|
549 | distZ = p.z() - zTopCut; |
---|
550 | } |
---|
551 | |
---|
552 | // Distance to closest surface from outside |
---|
553 | // |
---|
554 | if (distZ < 0.0) |
---|
555 | { |
---|
556 | return (distR < 0.0) ? 0.0 : distR; |
---|
557 | } |
---|
558 | else if (distR < 0.0) |
---|
559 | { |
---|
560 | return distZ; |
---|
561 | } |
---|
562 | else |
---|
563 | { |
---|
564 | return (distZ < distR) ? distZ : distR; |
---|
565 | } |
---|
566 | } |
---|
567 | |
---|
568 | /////////////////////////////////////////////////////////////////////////////// |
---|
569 | // |
---|
570 | // Calculate distance to surface of shape from `inside', allowing for tolerance |
---|
571 | |
---|
572 | G4double G4Ellipsoid::DistanceToOut(const G4ThreeVector& p, |
---|
573 | const G4ThreeVector& v, |
---|
574 | const G4bool calcNorm, |
---|
575 | G4bool *validNorm, |
---|
576 | G4ThreeVector *n ) const |
---|
577 | { |
---|
578 | G4double distMin; |
---|
579 | enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface; |
---|
580 | |
---|
581 | distMin= kInfinity; |
---|
582 | surface= kNoSurf; |
---|
583 | |
---|
584 | // check to see if Z plane is relevant |
---|
585 | // |
---|
586 | if (v.z() < 0.0) |
---|
587 | { |
---|
588 | G4double distZ = (zBottomCut - p.z()) / v.z(); |
---|
589 | if (distZ < 0.0) |
---|
590 | { |
---|
591 | distZ= 0.0; |
---|
592 | if (!calcNorm) {return 0.0;} |
---|
593 | } |
---|
594 | distMin= distZ; |
---|
595 | surface= kPlaneSurf; |
---|
596 | } |
---|
597 | if (v.z() > 0.0) |
---|
598 | { |
---|
599 | G4double distZ = (zTopCut - p.z()) / v.z(); |
---|
600 | if (distZ < 0.0) |
---|
601 | { |
---|
602 | distZ= 0.0; |
---|
603 | if (!calcNorm) {return 0.0;} |
---|
604 | } |
---|
605 | distMin= distZ; |
---|
606 | surface= kPlaneSurf; |
---|
607 | } |
---|
608 | |
---|
609 | // normal vector: parallel to normal, magnitude 1/(characteristic radius) |
---|
610 | // |
---|
611 | G4ThreeVector nearnorm(p.x()/(xSemiAxis*xSemiAxis), |
---|
612 | p.y()/(ySemiAxis*ySemiAxis), |
---|
613 | p.z()/(zSemiAxis*zSemiAxis)); |
---|
614 | |
---|
615 | // now check curved surface intercept |
---|
616 | // |
---|
617 | G4double A,B,C; |
---|
618 | |
---|
619 | A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis); |
---|
620 | C= (p * nearnorm) - 1.0; |
---|
621 | B= 2.0 * (v * nearnorm); |
---|
622 | |
---|
623 | C= B*B - 4.0*A*C; |
---|
624 | if (C > 0.0) |
---|
625 | { |
---|
626 | G4double distR= (-B + std::sqrt(C) ) / (2.0*A); |
---|
627 | if (distR < 0.0) |
---|
628 | { |
---|
629 | distR= 0.0; |
---|
630 | if (!calcNorm) {return 0.0;} |
---|
631 | } |
---|
632 | if (distR < distMin) |
---|
633 | { |
---|
634 | distMin= distR; |
---|
635 | surface= kCurvedSurf; |
---|
636 | } |
---|
637 | } |
---|
638 | |
---|
639 | // set normal if requested |
---|
640 | // |
---|
641 | if (calcNorm) |
---|
642 | { |
---|
643 | if (surface == kNoSurf) |
---|
644 | { |
---|
645 | *validNorm = false; |
---|
646 | } |
---|
647 | else |
---|
648 | { |
---|
649 | *validNorm = true; |
---|
650 | switch (surface) |
---|
651 | { |
---|
652 | case kPlaneSurf: |
---|
653 | *n= G4ThreeVector(0.,0.,(v.z() > 1.0 ? 1. : -1.)); |
---|
654 | break; |
---|
655 | case kCurvedSurf: |
---|
656 | { |
---|
657 | G4ThreeVector pexit= p + distMin*v; |
---|
658 | G4ThreeVector truenorm(pexit.x()/(xSemiAxis*xSemiAxis), |
---|
659 | pexit.y()/(ySemiAxis*ySemiAxis), |
---|
660 | pexit.z()/(zSemiAxis*zSemiAxis)); |
---|
661 | truenorm *= 1.0/truenorm.mag(); |
---|
662 | *n= truenorm; |
---|
663 | } break; |
---|
664 | default: |
---|
665 | G4cout.precision(16); |
---|
666 | G4cout << G4endl; |
---|
667 | DumpInfo(); |
---|
668 | G4cout << "Position:" << G4endl << G4endl; |
---|
669 | G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl; |
---|
670 | G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl; |
---|
671 | G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl; |
---|
672 | G4cout << "Direction:" << G4endl << G4endl; |
---|
673 | G4cout << "v.x() = " << v.x() << G4endl; |
---|
674 | G4cout << "v.y() = " << v.y() << G4endl; |
---|
675 | G4cout << "v.z() = " << v.z() << G4endl << G4endl; |
---|
676 | G4cout << "Proposed distance :" << G4endl << G4endl; |
---|
677 | G4cout << "distMin = " << distMin/mm << " mm" << G4endl << G4endl; |
---|
678 | G4Exception("G4Ellipsoid::DistanceToOut(p,v,..)", |
---|
679 | "Notification", JustWarning, |
---|
680 | "Undefined side for valid surface normal to solid."); |
---|
681 | break; |
---|
682 | } |
---|
683 | } |
---|
684 | } |
---|
685 | return distMin; |
---|
686 | } |
---|
687 | |
---|
688 | /////////////////////////////////////////////////////////////////////////////// |
---|
689 | // |
---|
690 | // Calculate distance (<=actual) to closest surface of shape from inside |
---|
691 | |
---|
692 | G4double G4Ellipsoid::DistanceToOut(const G4ThreeVector& p) const |
---|
693 | { |
---|
694 | G4double distR, distZ; |
---|
695 | |
---|
696 | #ifdef G4SPECSDEBUG |
---|
697 | if( Inside(p) == kOutside ) |
---|
698 | { |
---|
699 | G4cout.precision(16) ; |
---|
700 | G4cout << G4endl ; |
---|
701 | DumpInfo(); |
---|
702 | G4cout << "Position:" << G4endl << G4endl ; |
---|
703 | G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; |
---|
704 | G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; |
---|
705 | G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; |
---|
706 | G4Exception("G4Ellipsoid::DistanceToOut(p)", "Notification", JustWarning, |
---|
707 | "Point p is outside !?" ); |
---|
708 | } |
---|
709 | #endif |
---|
710 | |
---|
711 | // Normal vector: parallel to normal, magnitude 1/(characteristic radius) |
---|
712 | // |
---|
713 | G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis), |
---|
714 | p.y()/(ySemiAxis*ySemiAxis), |
---|
715 | p.z()/(zSemiAxis*zSemiAxis)); |
---|
716 | |
---|
717 | // the following is a safe inlined "radius= min(1.0/norm.mag(),p.mag()) |
---|
718 | // |
---|
719 | G4double radius= p.mag(); |
---|
720 | G4double tmp= norm.mag(); |
---|
721 | if ( (tmp > 0.0) && (1.0 < radius*tmp) ) {radius = 1.0/tmp;} |
---|
722 | |
---|
723 | // Approximate distance to curved surface ( <= actual distance ) |
---|
724 | // |
---|
725 | distR = (1.0 - p*norm) * radius / 2.0; |
---|
726 | |
---|
727 | // Distance to z-cut plane |
---|
728 | // |
---|
729 | distZ = p.z() - zBottomCut; |
---|
730 | if (distZ < 0.0) {distZ= zTopCut - p.z();} |
---|
731 | |
---|
732 | // Distance to closest surface from inside |
---|
733 | // |
---|
734 | if ( (distZ < 0.0) || (distR < 0.0) ) |
---|
735 | { |
---|
736 | return 0.0; |
---|
737 | } |
---|
738 | else |
---|
739 | { |
---|
740 | return (distZ < distR) ? distZ : distR; |
---|
741 | } |
---|
742 | } |
---|
743 | |
---|
744 | /////////////////////////////////////////////////////////////////////////////// |
---|
745 | // |
---|
746 | // Create a List containing the transformed vertices |
---|
747 | // Ordering [0-3] -fDz cross section |
---|
748 | // [4-7] +fDz cross section such that [0] is below [4], |
---|
749 | // [1] below [5] etc. |
---|
750 | // Note: |
---|
751 | // Caller has deletion resposibility |
---|
752 | // Potential improvement: For last slice, use actual ending angle |
---|
753 | // to avoid rounding error problems. |
---|
754 | |
---|
755 | G4ThreeVectorList* |
---|
756 | G4Ellipsoid::CreateRotatedVertices(const G4AffineTransform& pTransform, |
---|
757 | G4int& noPolygonVertices) const |
---|
758 | { |
---|
759 | G4ThreeVectorList *vertices; |
---|
760 | G4ThreeVector vertex; |
---|
761 | G4double meshAnglePhi, meshRMaxFactor, |
---|
762 | crossAnglePhi, coscrossAnglePhi, sincrossAnglePhi, sAnglePhi; |
---|
763 | G4double meshTheta, crossTheta, startTheta; |
---|
764 | G4double rMaxX, rMaxY, rMaxZ, rMaxMax, rx, ry, rz; |
---|
765 | G4int crossSectionPhi, noPhiCrossSections, crossSectionTheta, noThetaSections; |
---|
766 | |
---|
767 | // Phi cross sections |
---|
768 | // |
---|
769 | noPhiCrossSections=G4int (twopi/kMeshAngleDefault)+1; |
---|
770 | |
---|
771 | if (noPhiCrossSections<kMinMeshSections) |
---|
772 | { |
---|
773 | noPhiCrossSections=kMinMeshSections; |
---|
774 | } |
---|
775 | else if (noPhiCrossSections>kMaxMeshSections) |
---|
776 | { |
---|
777 | noPhiCrossSections=kMaxMeshSections; |
---|
778 | } |
---|
779 | meshAnglePhi=twopi/(noPhiCrossSections-1); |
---|
780 | |
---|
781 | // Set start angle such that mesh will be at fRMax |
---|
782 | // on the x axis. Will give better extent calculations when not rotated. |
---|
783 | |
---|
784 | sAnglePhi = -meshAnglePhi*0.5; |
---|
785 | |
---|
786 | // Theta cross sections |
---|
787 | |
---|
788 | noThetaSections = G4int(pi/kMeshAngleDefault)+3; |
---|
789 | |
---|
790 | if (noThetaSections<kMinMeshSections) |
---|
791 | { |
---|
792 | noThetaSections=kMinMeshSections; |
---|
793 | } |
---|
794 | else if (noThetaSections>kMaxMeshSections) |
---|
795 | { |
---|
796 | noThetaSections=kMaxMeshSections; |
---|
797 | } |
---|
798 | meshTheta= pi/(noThetaSections-2); |
---|
799 | |
---|
800 | // Set start angle such that mesh will be at fRMax |
---|
801 | // on the z axis. Will give better extent calculations when not rotated. |
---|
802 | |
---|
803 | startTheta = -meshTheta*0.5; |
---|
804 | |
---|
805 | meshRMaxFactor = 1.0/std::cos(0.5* |
---|
806 | std::sqrt(meshAnglePhi*meshAnglePhi+meshTheta*meshTheta)); |
---|
807 | rMaxMax= (xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis); |
---|
808 | if (zSemiAxis > rMaxMax) rMaxMax= zSemiAxis; |
---|
809 | rMaxX= xSemiAxis + rMaxMax*(meshRMaxFactor-1.0); |
---|
810 | rMaxY= ySemiAxis + rMaxMax*(meshRMaxFactor-1.0); |
---|
811 | rMaxZ= zSemiAxis + rMaxMax*(meshRMaxFactor-1.0); |
---|
812 | G4double* cosCrossTheta = new G4double[noThetaSections]; |
---|
813 | G4double* sinCrossTheta = new G4double[noThetaSections]; |
---|
814 | vertices=new G4ThreeVectorList(noPhiCrossSections*noThetaSections); |
---|
815 | if (vertices && cosCrossTheta && sinCrossTheta) |
---|
816 | { |
---|
817 | for (crossSectionTheta=0; crossSectionTheta<noThetaSections; |
---|
818 | crossSectionTheta++) |
---|
819 | { |
---|
820 | // Compute sine and cosine table (for historical reasons) |
---|
821 | // |
---|
822 | crossTheta=startTheta+crossSectionTheta*meshTheta; |
---|
823 | cosCrossTheta[crossSectionTheta]=std::cos(crossTheta); |
---|
824 | sinCrossTheta[crossSectionTheta]=std::sin(crossTheta); |
---|
825 | } |
---|
826 | for (crossSectionPhi=0; crossSectionPhi<noPhiCrossSections; |
---|
827 | crossSectionPhi++) |
---|
828 | { |
---|
829 | crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi; |
---|
830 | coscrossAnglePhi=std::cos(crossAnglePhi); |
---|
831 | sincrossAnglePhi=std::sin(crossAnglePhi); |
---|
832 | for (crossSectionTheta=0; crossSectionTheta<noThetaSections; |
---|
833 | crossSectionTheta++) |
---|
834 | { |
---|
835 | // Compute coordinates of cross section at section crossSectionPhi |
---|
836 | // |
---|
837 | rx= sinCrossTheta[crossSectionTheta]*coscrossAnglePhi*rMaxX; |
---|
838 | ry= sinCrossTheta[crossSectionTheta]*sincrossAnglePhi*rMaxY; |
---|
839 | rz= cosCrossTheta[crossSectionTheta]*rMaxZ; |
---|
840 | if (rz < zBottomCut) |
---|
841 | { rz= zBottomCut; } |
---|
842 | if (rz > zTopCut) |
---|
843 | { rz= zTopCut; } |
---|
844 | vertex= G4ThreeVector(rx,ry,rz); |
---|
845 | vertices->push_back(pTransform.TransformPoint(vertex)); |
---|
846 | } // Theta forward |
---|
847 | } // Phi |
---|
848 | noPolygonVertices = noThetaSections ; |
---|
849 | } |
---|
850 | else |
---|
851 | { |
---|
852 | DumpInfo(); |
---|
853 | G4Exception("G4Ellipsoid::CreateRotatedVertices()", |
---|
854 | "FatalError", FatalException, |
---|
855 | "Error in allocation of vertices. Out of memory !"); |
---|
856 | } |
---|
857 | |
---|
858 | delete[] cosCrossTheta; |
---|
859 | delete[] sinCrossTheta; |
---|
860 | |
---|
861 | return vertices; |
---|
862 | } |
---|
863 | |
---|
864 | ////////////////////////////////////////////////////////////////////////// |
---|
865 | // |
---|
866 | // G4EntityType |
---|
867 | |
---|
868 | G4GeometryType G4Ellipsoid::GetEntityType() const |
---|
869 | { |
---|
870 | return G4String("G4Ellipsoid"); |
---|
871 | } |
---|
872 | |
---|
873 | ////////////////////////////////////////////////////////////////////////// |
---|
874 | // |
---|
875 | // Stream object contents to an output stream |
---|
876 | |
---|
877 | std::ostream& G4Ellipsoid::StreamInfo( std::ostream& os ) const |
---|
878 | { |
---|
879 | os << "-----------------------------------------------------------\n" |
---|
880 | << " *** Dump for solid - " << GetName() << " ***\n" |
---|
881 | << " ===================================================\n" |
---|
882 | << " Solid type: G4Ellipsoid\n" |
---|
883 | << " Parameters: \n" |
---|
884 | |
---|
885 | << " semi-axis x: " << xSemiAxis/mm << " mm \n" |
---|
886 | << " semi-axis y: " << ySemiAxis/mm << " mm \n" |
---|
887 | << " semi-axis z: " << zSemiAxis/mm << " mm \n" |
---|
888 | << " max semi-axis: " << semiAxisMax/mm << " mm \n" |
---|
889 | << " lower cut plane level z: " << zBottomCut/mm << " mm \n" |
---|
890 | << " upper cut plane level z: " << zTopCut/mm << " mm \n" |
---|
891 | << "-----------------------------------------------------------\n"; |
---|
892 | |
---|
893 | return os; |
---|
894 | } |
---|
895 | |
---|
896 | //////////////////////////////////////////////////////////////////// |
---|
897 | // |
---|
898 | // GetPointOnSurface |
---|
899 | |
---|
900 | G4ThreeVector G4Ellipsoid::GetPointOnSurface() const |
---|
901 | { |
---|
902 | G4double aTop, aBottom, aCurved, chose, xRand, yRand, zRand, phi, theta; |
---|
903 | G4double cosphi, sinphi, costheta, sintheta, alpha, beta, max1, max2, max3; |
---|
904 | |
---|
905 | max1 = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis; |
---|
906 | max1 = max1 > zSemiAxis ? max1 : zSemiAxis; |
---|
907 | if(max1 == xSemiAxis){max2 = ySemiAxis; max3 = zSemiAxis;} |
---|
908 | else if(max1 == ySemiAxis){max2 = xSemiAxis; max3 = zSemiAxis;} |
---|
909 | else {max2 = xSemiAxis; max3 = ySemiAxis; } |
---|
910 | |
---|
911 | phi = RandFlat::shoot(0.,2.*pi); |
---|
912 | theta = RandFlat::shoot(0.,pi); |
---|
913 | |
---|
914 | cosphi = std::cos(phi); sinphi = std::sin(phi); |
---|
915 | costheta = RandFlat::shoot(zBottomCut,zTopCut)/zSemiAxis; |
---|
916 | sintheta = std::sqrt(1.-sqr(costheta)); |
---|
917 | |
---|
918 | alpha = 1.-sqr(max2/max1); beta = 1.-sqr(max3/max1); |
---|
919 | |
---|
920 | aTop = pi*xSemiAxis*ySemiAxis*(1 - sqr(zTopCut/zSemiAxis)); |
---|
921 | aBottom = pi*xSemiAxis*ySemiAxis*(1 - sqr(zBottomCut/zSemiAxis)); |
---|
922 | |
---|
923 | // approximation |
---|
924 | // from:" http://www.citr.auckland.ac.nz/techreports/2004/CITR-TR-139.pdf" |
---|
925 | aCurved = 4.*pi*max1*max2*(1.-1./6.*(alpha+beta)- |
---|
926 | 1./120.*(3.*sqr(alpha)+2.*alpha*beta+3.*sqr(beta))); |
---|
927 | |
---|
928 | aCurved *= 0.5*(1.2*zTopCut/zSemiAxis - 1.2*zBottomCut/zSemiAxis); |
---|
929 | |
---|
930 | if( ( zTopCut >= zSemiAxis && zBottomCut <= -1.*zSemiAxis ) |
---|
931 | || ( zTopCut == 0 && zBottomCut ==0 ) ) |
---|
932 | { |
---|
933 | aTop = 0; aBottom = 0; |
---|
934 | } |
---|
935 | |
---|
936 | chose = RandFlat::shoot(0.,aTop + aBottom + aCurved); |
---|
937 | |
---|
938 | if(chose < aCurved) |
---|
939 | { |
---|
940 | xRand = xSemiAxis*sintheta*cosphi; |
---|
941 | yRand = ySemiAxis*sintheta*sinphi; |
---|
942 | zRand = zSemiAxis*costheta; |
---|
943 | return G4ThreeVector (xRand,yRand,zRand); |
---|
944 | } |
---|
945 | else if(chose >= aCurved && chose < aCurved + aTop) |
---|
946 | { |
---|
947 | xRand = RandFlat::shoot(-1.,1.)*xSemiAxis |
---|
948 | * std::sqrt(1-sqr(zTopCut/zSemiAxis)); |
---|
949 | yRand = RandFlat::shoot(-1.,1.)*ySemiAxis |
---|
950 | * std::sqrt(1.-sqr(zTopCut/zSemiAxis)-sqr(xRand/xSemiAxis)); |
---|
951 | zRand = zTopCut; |
---|
952 | return G4ThreeVector (xRand,yRand,zRand); |
---|
953 | } |
---|
954 | else |
---|
955 | { |
---|
956 | xRand = RandFlat::shoot(-1.,1.)*xSemiAxis |
---|
957 | * std::sqrt(1-sqr(zBottomCut/zSemiAxis)); |
---|
958 | yRand = RandFlat::shoot(-1.,1.)*ySemiAxis |
---|
959 | * std::sqrt(1.-sqr(zBottomCut/zSemiAxis)-sqr(xRand/xSemiAxis)); |
---|
960 | zRand = zBottomCut; |
---|
961 | return G4ThreeVector (xRand,yRand,zRand); |
---|
962 | } |
---|
963 | } |
---|
964 | |
---|
965 | ///////////////////////////////////////////////////////////////////////////// |
---|
966 | // |
---|
967 | // Methods for visualisation |
---|
968 | |
---|
969 | void G4Ellipsoid::DescribeYourselfTo (G4VGraphicsScene& scene) const |
---|
970 | { |
---|
971 | scene.AddSolid(*this); |
---|
972 | } |
---|
973 | |
---|
974 | G4VisExtent G4Ellipsoid::GetExtent() const |
---|
975 | { |
---|
976 | // Define the sides of the box into which the G4Ellipsoid instance would fit. |
---|
977 | // |
---|
978 | return G4VisExtent (-semiAxisMax, semiAxisMax, |
---|
979 | -semiAxisMax, semiAxisMax, |
---|
980 | -semiAxisMax, semiAxisMax); |
---|
981 | } |
---|
982 | |
---|
983 | G4NURBS* G4Ellipsoid::CreateNURBS () const |
---|
984 | { |
---|
985 | // Box for now!!! |
---|
986 | // |
---|
987 | return new G4NURBSbox(semiAxisMax, semiAxisMax, semiAxisMax); |
---|
988 | } |
---|
989 | |
---|
990 | G4Polyhedron* G4Ellipsoid::CreatePolyhedron () const |
---|
991 | { |
---|
992 | return new G4PolyhedronEllipsoid(xSemiAxis, ySemiAxis, zSemiAxis, |
---|
993 | zBottomCut, zTopCut); |
---|
994 | } |
---|
995 | |
---|
996 | G4Polyhedron* G4Ellipsoid::GetPolyhedron () const |
---|
997 | { |
---|
998 | if (!fpPolyhedron || |
---|
999 | fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() != |
---|
1000 | fpPolyhedron->GetNumberOfRotationSteps()) |
---|
1001 | { |
---|
1002 | delete fpPolyhedron; |
---|
1003 | fpPolyhedron = CreatePolyhedron(); |
---|
1004 | } |
---|
1005 | return fpPolyhedron; |
---|
1006 | } |
---|