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