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