| 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 //
26 // G4EllipticalTube implementation <<
27 // 23 //
28 // Author: David C. Williams (davidw@scipp.ucs << 24 // $Id: G4EllipticalTube.cc,v 1.11 2002/05/08 13:51:26 johna Exp $
29 // Revision: Evgueni Tcherniaev (evgueni.tcher << 25 // GEANT4 tag $Name: geant4-04-01 $
>> 26 //
>> 27 //
>> 28 // --------------------------------------------------------------------
>> 29 // GEANT 4 class source file
>> 30 //
>> 31 //
>> 32 // G4EllipticalTube.cc
>> 33 //
>> 34 // Implementation of a CSG volume representing a tube with elliptical cross
>> 35 // section (geant3 solid 'ELTU')
>> 36 //
30 // ------------------------------------------- 37 // --------------------------------------------------------------------
31 38
32 #include "G4EllipticalTube.hh" 39 #include "G4EllipticalTube.hh"
33 <<
34 #if !(defined(G4GEOM_USE_UELLIPTICALTUBE) && d <<
35 <<
36 #include "G4GeomTools.hh" <<
37 #include "G4RandomTools.hh" <<
38 #include "G4ClippablePolygon.hh" 40 #include "G4ClippablePolygon.hh"
39 #include "G4AffineTransform.hh" 41 #include "G4AffineTransform.hh"
>> 42 #include "G4SolidExtentList.hh"
40 #include "G4VoxelLimits.hh" 43 #include "G4VoxelLimits.hh"
41 #include "G4BoundingEnvelope.hh" << 44 #include "meshdefs.hh"
42 <<
43 #include "Randomize.hh" <<
44 45
45 #include "G4VGraphicsScene.hh" 46 #include "G4VGraphicsScene.hh"
>> 47 #include "G4Polyhedron.hh"
46 #include "G4VisExtent.hh" 48 #include "G4VisExtent.hh"
47 49
48 #include "G4AutoLock.hh" <<
49 <<
50 namespace <<
51 { <<
52 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE <<
53 } <<
54 <<
55 using namespace CLHEP; <<
56 <<
57 ////////////////////////////////////////////// <<
58 // 50 //
59 // Constructor 51 // Constructor
60 <<
61 G4EllipticalTube::G4EllipticalTube( const G4St <<
62 G4do <<
63 G4do <<
64 G4do <<
65 : G4VSolid(name), fDx(Dx), fDy(Dy), fDz(Dz) <<
66 { <<
67 CheckParameters(); <<
68 } <<
69 <<
70 ////////////////////////////////////////////// <<
71 // 52 //
72 // Fake default constructor - sets only member << 53 G4EllipticalTube::G4EllipticalTube( const G4String &name,
73 // for usage restri << 54 G4double theDx, G4double theDy, G4double theDz )
74 << 55 : G4VSolid( name )
75 G4EllipticalTube::G4EllipticalTube( __void__& << 56 {
76 : G4VSolid(a), halfTolerance(0.), fDx(0.), f << 57 dx = theDx;
77 fRsph(0.), fDDx(0.), fDDy(0.), fSx(0.), fS << 58 dy = theDy;
78 fQ1(0.), fQ2(0.), fScratch(0.) << 59 dz = theDz;
79 { <<
80 } 60 }
81 61
82 ////////////////////////////////////////////// << 62
83 // 63 //
84 // Destructor 64 // Destructor
85 <<
86 G4EllipticalTube::~G4EllipticalTube() <<
87 { <<
88 delete fpPolyhedron; fpPolyhedron = nullptr; <<
89 } <<
90 <<
91 ////////////////////////////////////////////// <<
92 // 65 //
93 // Copy constructor << 66 G4EllipticalTube::~G4EllipticalTube() {;}
94 67
95 G4EllipticalTube::G4EllipticalTube(const G4Ell <<
96 : G4VSolid(rhs), halfTolerance(rhs.halfToler <<
97 fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), <<
98 fCubicVolume(rhs.fCubicVolume), fSurfaceAr <<
99 fRsph(rhs.fRsph), fDDx(rhs.fDDx), fDDy(rhs <<
100 fSx(rhs.fSx), fSy(rhs.fSy), fR(rhs.fR), <<
101 fQ1(rhs.fQ1), fQ2(rhs.fQ2), fScratch(rhs.f <<
102 { <<
103 } <<
104 <<
105 ////////////////////////////////////////////// <<
106 // <<
107 // Assignment operator <<
108 <<
109 G4EllipticalTube& G4EllipticalTube::operator = <<
110 { <<
111 // Check assignment to self <<
112 // <<
113 if (this == &rhs) { return *this; } <<
114 <<
115 // Copy base class data <<
116 // <<
117 G4VSolid::operator=(rhs); <<
118 <<
119 // Copy data <<
120 // <<
121 halfTolerance = rhs.halfTolerance; <<
122 fDx = rhs.fDx; <<
123 fDy = rhs.fDy; <<
124 fDz = rhs.fDz; <<
125 fCubicVolume = rhs.fCubicVolume; <<
126 fSurfaceArea = rhs.fSurfaceArea; <<
127 <<
128 fRsph = rhs.fRsph; <<
129 fDDx = rhs.fDDx; <<
130 fDDy = rhs.fDDy; <<
131 fSx = rhs.fSx; <<
132 fSy = rhs.fSy; <<
133 fR = rhs.fR; <<
134 fQ1 = rhs.fQ1; <<
135 fQ2 = rhs.fQ2; <<
136 fScratch = rhs.fScratch; <<
137 <<
138 fRebuildPolyhedron = false; <<
139 delete fpPolyhedron; fpPolyhedron = nullptr <<
140 <<
141 return *this; <<
142 } <<
143 <<
144 ////////////////////////////////////////////// <<
145 // <<
146 // Check dimensions <<
147 <<
148 void G4EllipticalTube::CheckParameters() <<
149 { <<
150 // Check dimensions <<
151 // <<
152 halfTolerance = 0.5*kCarTolerance; // half t <<
153 G4double dmin = 2*kCarTolerance; <<
154 if (fDx < dmin || fDy < dmin || fDz < dmin) <<
155 { <<
156 std::ostringstream message; <<
157 message << "Invalid (too small or negative <<
158 << GetName() <<
159 << "\n Dx = " << fDx <<
160 << "\n Dy = " << fDy <<
161 << "\n Dz = " << fDz; <<
162 G4Exception("G4EllipticalTube::CheckParame <<
163 FatalException, message); <<
164 } <<
165 <<
166 // Set pre-calculatated values <<
167 // <<
168 halfTolerance = 0.5*kCarTolerance; // half t <<
169 fRsph = std::sqrt(fDx * fDx + fDy * fDy + fD <<
170 fDDx = fDx * fDx; // X semi-axis squared <<
171 fDDy = fDy * fDy; // Y semi-axis squared <<
172 <<
173 fR = std::min(fDx, fDy); // resulting radius <<
174 fSx = fR / fDx; // X scale factor <<
175 fSy = fR / fDy; // Y scale factor <<
176 <<
177 fQ1 = 0.5 / fR; // distance approxiamtion di <<
178 fQ2 = 0.5 * (fR + halfTolerance * halfTolera <<
179 fScratch = 2. * fR * fR * DBL_EPSILON; // sc <<
180 // fScratch = (B * B / A) * (2. + halfTolera <<
181 } <<
182 <<
183 ////////////////////////////////////////////// <<
184 // <<
185 // Get bounding box <<
186 <<
187 void G4EllipticalTube::BoundingLimits( G4Three <<
188 G4Three <<
189 { <<
190 pMin.set(-fDx,-fDy,-fDz); <<
191 pMax.set( fDx, fDy, fDz); <<
192 } <<
193 <<
194 ////////////////////////////////////////////// <<
195 // <<
196 // Calculate extent under transform and specif <<
197 <<
198 G4bool <<
199 G4EllipticalTube::CalculateExtent( const EAxis <<
200 const G4Vox <<
201 const G4Aff <<
202 G4dou <<
203 { <<
204 G4ThreeVector bmin, bmax; <<
205 G4bool exist; <<
206 <<
207 // Check bounding box (bbox) <<
208 // <<
209 BoundingLimits(bmin,bmax); <<
210 G4BoundingEnvelope bbox(bmin,bmax); <<
211 #ifdef G4BBOX_EXTENT <<
212 return bbox.CalculateExtent(pAxis,pVoxelLimi <<
213 #endif <<
214 if (bbox.BoundingBoxVsVoxelLimits(pAxis, pVo <<
215 { <<
216 return exist = pMin < pMax; <<
217 } <<
218 <<
219 G4double dx = fDx; <<
220 G4double dy = fDy; <<
221 G4double dz = fDz; <<
222 <<
223 // Set bounding envelope (benv) and calculat <<
224 // <<
225 const G4int NSTEPS = 24; // number of steps <<
226 G4double ang = twopi/NSTEPS; <<
227 <<
228 G4double sinHalf = std::sin(0.5*ang); <<
229 G4double cosHalf = std::cos(0.5*ang); <<
230 G4double sinStep = 2.*sinHalf*cosHalf; <<
231 G4double cosStep = 1. - 2.*sinHalf*sinHalf; <<
232 G4double sx = dx/cosHalf; <<
233 G4double sy = dy/cosHalf; <<
234 <<
235 G4double sinCur = sinHalf; <<
236 G4double cosCur = cosHalf; <<
237 G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS <<
238 for (G4int k=0; k<NSTEPS; ++k) <<
239 { <<
240 baseA[k].set(sx*cosCur,sy*sinCur,-dz); <<
241 baseB[k].set(sx*cosCur,sy*sinCur, dz); <<
242 <<
243 G4double sinTmp = sinCur; <<
244 sinCur = sinCur*cosStep + cosCur*sinStep; <<
245 cosCur = cosCur*cosStep - sinTmp*sinStep; <<
246 } <<
247 <<
248 std::vector<const G4ThreeVectorList *> polyg <<
249 polygons[0] = &baseA; <<
250 polygons[1] = &baseB; <<
251 G4BoundingEnvelope benv(bmin, bmax, polygons <<
252 exist = benv.CalculateExtent(pAxis, pVoxelLi <<
253 return exist; <<
254 } <<
255 <<
256 ////////////////////////////////////////////// <<
257 // <<
258 // Determine where is point: inside, outside o <<
259 // <<
260 <<
261 EInside G4EllipticalTube::Inside( const G4Thre <<
262 { <<
263 G4double x = p.x() * fSx; <<
264 G4double y = p.y() * fSy; <<
265 G4double distR = fQ1 * (x * x + y * y) - fQ2 <<
266 G4double distZ = std::abs(p.z()) - fDz; <<
267 G4double dist = std::max(distR, distZ); <<
268 <<
269 if (dist > halfTolerance) return kOutside; <<
270 return (dist > -halfTolerance) ? kSurface : <<
271 } <<
272 <<
273 ////////////////////////////////////////////// <<
274 // <<
275 // Return unit normal at surface closest to p <<
276 <<
277 G4ThreeVector G4EllipticalTube::SurfaceNormal( <<
278 { <<
279 G4ThreeVector norm(0, 0, 0); <<
280 G4int nsurf = 0; <<
281 <<
282 // check lateral surface <<
283 G4double x = p.x() * fSx; <<
284 G4double y = p.y() * fSy; <<
285 G4double distR = fQ1 * (x * x + y * y) - fQ2 <<
286 if (std::abs(distR) <= halfTolerance) <<
287 { <<
288 norm = G4ThreeVector(p.x() * fDDy, p.y() * <<
289 ++nsurf; <<
290 } <<
291 <<
292 // check lateral bases <<
293 G4double distZ = std::abs(p.z()) - fDz; <<
294 if (std::abs(distZ) <= halfTolerance) <<
295 { <<
296 norm.setZ(p.z() < 0 ? -1. : 1.); <<
297 ++nsurf; <<
298 } <<
299 <<
300 // return normal <<
301 if (nsurf == 1) return norm; <<
302 else if (nsurf > 1) return norm.unit(); // e <<
303 else <<
304 { <<
305 // Point is not on the surface <<
306 // <<
307 #ifdef G4SPECDEBUG <<
308 std::ostringstream message; <<
309 G4long oldprc = message.precision(16); <<
310 message << "Point p is not on surface (!?) <<
311 << GetName() << G4endl; <<
312 message << "Position:\n"; <<
313 message << " p.x() = " << p.x()/mm << " <<
314 message << " p.y() = " << p.y()/mm << " <<
315 message << " p.z() = " << p.z()/mm << " <<
316 G4cout.precision(oldprc); <<
317 G4Exception("G4EllipticalTube::SurfaceNorm <<
318 JustWarning, message ); <<
319 DumpInfo(); <<
320 #endif <<
321 return ApproxSurfaceNormal(p); <<
322 } <<
323 } <<
324 <<
325 ////////////////////////////////////////////// <<
326 // <<
327 // Find surface nearest to point and return co <<
328 // The algorithm is similar to the algorithm u <<
329 // This method normally should not be called. <<
330 <<
331 G4ThreeVector <<
332 G4EllipticalTube::ApproxSurfaceNormal( const G <<
333 { <<
334 G4double x = p.x() * fSx; <<
335 G4double y = p.y() * fSy; <<
336 G4double distR = fQ1 * (x * x + y * y) - fQ2 <<
337 G4double distZ = std::abs(p.z()) - fDz; <<
338 if (distR > distZ && (x * x + y * y) > 0) <<
339 return G4ThreeVector(p.x() * fDDy, p.y() * <<
340 else <<
341 return {0, 0, (p.z() < 0 ? -1. : 1.)}; <<
342 } <<
343 <<
344 ////////////////////////////////////////////// <<
345 // <<
346 // Calculate distance to shape from outside, a <<
347 // return kInfinity if no intersection, or dis <<
348 <<
349 G4double G4EllipticalTube::DistanceToIn( const <<
350 const <<
351 { <<
352 G4double offset = 0.; <<
353 G4ThreeVector pcur = p; <<
354 <<
355 // Check if point is flying away <<
356 // <<
357 G4double safex = std::abs(pcur.x()) - fDx; <<
358 G4double safey = std::abs(pcur.y()) - fDy; <<
359 G4double safez = std::abs(pcur.z()) - fDz; <<
360 <<
361 if (safez >= -halfTolerance && pcur.z() * v. <<
362 if (safey >= -halfTolerance && pcur.y() * v. <<
363 if (safex >= -halfTolerance && pcur.x() * v. <<
364 <<
365 // Relocate point, if required <<
366 // <<
367 G4double Dmax = 32. * fRsph; <<
368 if (std::max(std::max(safex, safey), safez) <<
369 { <<
370 offset = (1. - 1.e-08) * pcur.mag() - 2. * <<
371 pcur += offset * v; <<
372 G4double dist = DistanceToIn(pcur, v); <<
373 return (dist == kInfinity) ? kInfinity : d <<
374 } <<
375 <<
376 // Scale elliptical tube to cylinder <<
377 // <<
378 G4double px = pcur.x() * fSx; <<
379 G4double py = pcur.y() * fSy; <<
380 G4double pz = pcur.z(); <<
381 G4double vx = v.x() * fSx; <<
382 G4double vy = v.y() * fSy; <<
383 G4double vz = v.z(); <<
384 <<
385 // Set coefficients of quadratic equation: A <<
386 // <<
387 G4double rr = px * px + py * py; <<
388 G4double A = vx * vx + vy * vy; <<
389 G4double B = px * vx + py * vy; <<
390 G4double C = rr - fR * fR; <<
391 G4double D = B * B - A * C; <<
392 <<
393 // Check if point is flying away relative to <<
394 // <<
395 G4double distR = fQ1 * rr - fQ2; <<
396 G4bool parallelToZ = (A < DBL_EPSILON || std <<
397 if (distR >= -halfTolerance && (B >= 0. || p <<
398 <<
399 // Find intersection with Z planes <<
400 // <<
401 G4double invz = (vz == 0) ? DBL_MAX : -1./v <<
402 G4double dz = std::copysign(fDz, invz); <<
403 G4double tzmin = (pz - dz) * invz; <<
404 G4double tzmax = (pz + dz) * invz; <<
405 <<
406 // Solve qudratic equation. There are two ca <<
407 // 1) trajectory parallel to Z axis (A = 0 <<
408 // 2) touch (D = 0) or no intersection (D <<
409 // <<
410 if (parallelToZ) return (tzmin<halfTolerance <<
411 if (D <= A * A * fScratch) return kInfinity; <<
412 <<
413 // Find roots of quadratic equation <<
414 G4double tmp = -B - std::copysign(std::sqrt( <<
415 G4double t1 = tmp / A; <<
416 G4double t2 = C / tmp; <<
417 G4double trmin = std::min(t1, t2); <<
418 G4double trmax = std::max(t1, t2); <<
419 <<
420 // Return distance <<
421 G4double tin = std::max(tzmin, trmin); <<
422 G4double tout = std::min(tzmax, trmax); <<
423 68
424 if (tout <= tin + halfTolerance) return kInf <<
425 return (tin<halfTolerance) ? offset : tin + <<
426 } <<
427 <<
428 ////////////////////////////////////////////// <<
429 // 69 //
430 // Estimate distance to the surface from outsi << 70 // CalculateExtent
431 // returns 0 if point is inside << 71 //
432 << 72 G4bool G4EllipticalTube::CalculateExtent( const EAxis axis,
>> 73 const G4VoxelLimits &voxelLimit,
>> 74 const G4AffineTransform &transform,
>> 75 G4double &min, G4double &max ) const
>> 76 {
>> 77 G4SolidExtentList extentList( axis, voxelLimit );
>> 78
>> 79 //
>> 80 // We are going to divide up our elliptical face into small
>> 81 // pieces
>> 82 //
>> 83
>> 84 //
>> 85 // Choose phi size of our segment(s) based on constants as
>> 86 // defined in meshdefs.hh
>> 87 //
>> 88 G4int numPhi = kMaxMeshSections;
>> 89 G4double sigPhi = 2*M_PI/numPhi;
>> 90
>> 91 //
>> 92 // We have to be careful to keep our segments completely outside
>> 93 // of the elliptical surface. To do so we imagine we have
>> 94 // a simple (unit radius) circular cross section (as in G4Tubs)
>> 95 // and then "stretch" the dimensions as necessary to fit the ellipse.
>> 96 //
>> 97 G4double rFudge = 1.0/cos(0.5*sigPhi);
>> 98 G4double dxFudge = dx*rFudge,
>> 99 dyFudge = dy*rFudge;
>> 100
>> 101 //
>> 102 // As we work around the elliptical surface, we build
>> 103 // a "phi" segment on the way, and keep track of two
>> 104 // additional polygons for the two ends.
>> 105 //
>> 106 G4ClippablePolygon endPoly1, endPoly2, phiPoly;
>> 107
>> 108 G4double phi = 0,
>> 109 cosPhi = cos(phi),
>> 110 sinPhi = sin(phi);
>> 111 G4ThreeVector v0( dxFudge*cosPhi, dyFudge*sinPhi, +dz ),
>> 112 v1( dxFudge*cosPhi, dyFudge*sinPhi, -dz ),
>> 113 w0, w1;
>> 114 transform.ApplyPointTransform( v0 );
>> 115 transform.ApplyPointTransform( v1 );
>> 116 do {
>> 117 phi += sigPhi;
>> 118 if (numPhi == 1) phi = 0; // Try to avoid roundoff
>> 119 cosPhi = cos(phi),
>> 120 sinPhi = sin(phi);
>> 121
>> 122 w0 = G4ThreeVector( dxFudge*cosPhi, dyFudge*sinPhi, +dz );
>> 123 w1 = G4ThreeVector( dxFudge*cosPhi, dyFudge*sinPhi, -dz );
>> 124 transform.ApplyPointTransform( w0 );
>> 125 transform.ApplyPointTransform( w1 );
>> 126
>> 127 //
>> 128 // Add a point to our z ends
>> 129 //
>> 130 endPoly1.AddVertexInOrder( v0 );
>> 131 endPoly2.AddVertexInOrder( v1 );
>> 132
>> 133 //
>> 134 // Build phi polygon
>> 135 //
>> 136 phiPoly.ClearAllVertices();
>> 137
>> 138 phiPoly.AddVertexInOrder( v0 );
>> 139 phiPoly.AddVertexInOrder( v1 );
>> 140 phiPoly.AddVertexInOrder( w1 );
>> 141 phiPoly.AddVertexInOrder( w0 );
>> 142
>> 143 if (phiPoly.PartialClip( voxelLimit, axis )) {
>> 144 //
>> 145 // Get unit normal
>> 146 //
>> 147 phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
>> 148
>> 149 extentList.AddSurface( phiPoly );
>> 150 }
>> 151
>> 152 //
>> 153 // Next vertex
>> 154 //
>> 155 v0 = w0;
>> 156 v1 = w1;
>> 157 } while( --numPhi > 0 );
>> 158
>> 159 //
>> 160 // Process the end pieces
>> 161 //
>> 162 if (endPoly1.PartialClip( voxelLimit, axis )) {
>> 163 static const G4ThreeVector normal(0,0,+1);
>> 164 endPoly1.SetNormal( transform.TransformAxis(normal) );
>> 165 extentList.AddSurface( endPoly1 );
>> 166 }
>> 167
>> 168 if (endPoly2.PartialClip( voxelLimit, axis )) {
>> 169 static const G4ThreeVector normal(0,0,-1);
>> 170 endPoly2.SetNormal( transform.TransformAxis(normal) );
>> 171 extentList.AddSurface( endPoly2 );
>> 172 }
>> 173
>> 174 //
>> 175 // Return min/max value
>> 176 //
>> 177 return extentList.GetExtent( min, max );
>> 178 }
>> 179
>> 180
>> 181 //
>> 182 // Inside
>> 183 //
>> 184 // Note that for this solid, we've decided to define the tolerant
>> 185 // surface as that which is bounded by ellipses with axes
>> 186 // at +/- 0.5*kCarTolerance.
>> 187 //
>> 188 EInside G4EllipticalTube::Inside( const G4ThreeVector& p) const
>> 189 {
>> 190 static const G4double halfTol = 0.5*kCarTolerance;
>> 191
>> 192 //
>> 193 // Check z extents: are we outside?
>> 194 //
>> 195 G4double absZ = fabs(p.z());
>> 196 if (absZ > dz+halfTol) return kOutside;
>> 197
>> 198 //
>> 199 // Check x,y: are we outside?
>> 200 //
>> 201 // G4double x = p.x(), y = p.y();
>> 202
>> 203 if (CheckXY(p.x(), p.y(), +halfTol) > 1.0) return kOutside;
>> 204
>> 205 //
>> 206 // We are either inside or on the surface: recheck z extents
>> 207 //
>> 208 if (absZ > dz-halfTol) return kSurface;
>> 209
>> 210 //
>> 211 // Recheck x,y
>> 212 //
>> 213 if (CheckXY(p.x(), p.y(), -halfTol) > 1.0) return kSurface;
>> 214
>> 215 return kInside;
>> 216 }
>> 217
>> 218
>> 219 //
>> 220 // SurfaceNormal
>> 221 //
>> 222 G4ThreeVector G4EllipticalTube::SurfaceNormal( const G4ThreeVector& p) const
>> 223 {
>> 224 //
>> 225 // Which of the three surfaces are we closest to (approximately)?
>> 226 //
>> 227 G4double distZ = fabs(p.z()) - dz;
>> 228
>> 229 G4double rxy = CheckXY( p.x(), p.y() );
>> 230 G4double distR2 = (rxy < DBL_MIN) ? DBL_MAX : 1.0/rxy;
>> 231
>> 232 //
>> 233 // Closer to z?
>> 234 //
>> 235 if (distZ*distZ < distR2)
>> 236 return G4ThreeVector( 0.0, 0.0, p.z() < 0 ? -1.0 : 1.0 );
>> 237
>> 238 //
>> 239 // Closer to x/y
>> 240 //
>> 241 return G4ThreeVector( p.x()*dy*dy, p.y()*dx*dx, 0.0 ).unit();
>> 242 }
>> 243
>> 244
>> 245 //
>> 246 // DistanceToIn(p,v)
>> 247 //
>> 248 // Unlike DistanceToOut(p,v), it is possible for the trajectory
>> 249 // to miss. The geometric calculations here are quite simple.
>> 250 // More difficult is the logic required to prevent particles
>> 251 // from sneaking (or leaking) between the elliptical and end
>> 252 // surfaces.
>> 253 //
>> 254 // Keep in mind that the true distance is allowed to be
>> 255 // negative if the point is currently on the surface. For oblique
>> 256 // angles, it can be very negative.
>> 257 //
>> 258 G4double G4EllipticalTube::DistanceToIn( const G4ThreeVector& p,const G4ThreeVector& v ) const
>> 259 {
>> 260 static const G4double halfTol = 0.5*kCarTolerance;
>> 261
>> 262 //
>> 263 // Check z = -dz planer surface
>> 264 //
>> 265 G4double sigz = p.z()+dz;
>> 266
>> 267 if (sigz < halfTol) {
>> 268 //
>> 269 // We are "behind" the shape in z, and so can
>> 270 // potentially hit the rear face. Correct direction?
>> 271 //
>> 272 if (v.z() <= 0) {
>> 273 //
>> 274 // As long as we are far enough away, we know we
>> 275 // can't intersect
>> 276 //
>> 277 if (sigz < 0) return kInfinity;
>> 278
>> 279 //
>> 280 // Otherwise, we don't intersect unless we are
>> 281 // on the surface of the ellipse
>> 282 //
>> 283 if (CheckXY(p.x(),p.y(),-halfTol) <= 1.0) return kInfinity;
>> 284 }
>> 285 else {
>> 286
>> 287 //
>> 288 // How far?
>> 289 //
>> 290 G4double s = -sigz/v.z();
>> 291
>> 292 //
>> 293 // Where does that place us?
>> 294 //
>> 295 G4double xi = p.x() + s*v.x(),
>> 296 yi = p.y() + s*v.y();
>> 297
>> 298 //
>> 299 // Is this on the surface (within ellipse)?
>> 300 //
>> 301 if (CheckXY(xi,yi) <= 1.0) {
>> 302 //
>> 303 // Yup. Return s, unless we are on the surface
>> 304 //
>> 305 return (sigz < -halfTol) ? s : 0;
>> 306 }
>> 307 else if (xi*dy*dy*v.x() + yi*dx*dx*v.y() >= 0) {
>> 308 //
>> 309 // Else, if we are traveling outwards, we know
>> 310 // we must miss
>> 311 //
>> 312 return kInfinity;
>> 313 }
>> 314 }
>> 315 }
>> 316
>> 317 //
>> 318 // Check z = +dz planer surface
>> 319 //
>> 320 sigz = p.z() - dz;
>> 321
>> 322 if (sigz > -halfTol) {
>> 323 if (v.z() >= 0) {
>> 324 if (sigz > 0) return kInfinity;
>> 325 if (CheckXY(p.x(),p.y(),-halfTol) <= 1.0) return kInfinity;
>> 326 }
>> 327 else {
>> 328 G4double s = -sigz/v.z();
>> 329
>> 330 G4double xi = p.x() + s*v.x(),
>> 331 yi = p.y() + s*v.y();
>> 332
>> 333 if (CheckXY(xi,yi) <= 1.0) {
>> 334 return (sigz > -halfTol) ? s : 0;
>> 335 }
>> 336 else if (xi*dy*dy*v.x() + yi*dx*dx*v.y() >= 0) {
>> 337 return kInfinity;
>> 338 }
>> 339 }
>> 340 }
>> 341
>> 342 //
>> 343 // Check intersection with the elliptical tube
>> 344 //
>> 345 G4double s[2];
>> 346 G4int n = IntersectXY( p, v, s );
>> 347
>> 348 if (n==0) return kInfinity;
>> 349
>> 350 //
>> 351 // Is the original point on the surface?
>> 352 //
>> 353 if (fabs(p.z()) < dz+halfTol) {
>> 354 if (CheckXY( p.x(), p.y(), halfTol ) < 1.0) {
>> 355 //
>> 356 // Well, yes, but are we traveling inwards at this point?
>> 357 //
>> 358 if (p.x()*dy*dy*v.x() + p.y()*dx*dx*v.y() < 0) return 0;
>> 359 }
>> 360 }
>> 361
>> 362 //
>> 363 // We are now certain that point p is not on the surface of
>> 364 // the solid (and thus fabs(s[0]) > halfTol).
>> 365 // Return kInfinity if the intersection is "behind" the point.
>> 366 //
>> 367 if (s[0] < 0) return kInfinity;
>> 368
>> 369 //
>> 370 // Check to see if we intersect the tube within
>> 371 // dz, but only when we know it might miss
>> 372 //
>> 373 G4double zi = p.z() + s[0]*v.z();
>> 374
>> 375 if (v.z() < 0) {
>> 376 if (zi < -dz) return kInfinity;
>> 377 }
>> 378 else if (v.z() > 0) {
>> 379 if (zi > +dz) return kInfinity;
>> 380 }
>> 381
>> 382 return s[0];
>> 383 }
>> 384
>> 385
>> 386 //
>> 387 // DistanceToIn(p)
>> 388 //
>> 389 // The distance from a point to an ellipse (in 2 dimensions) is a
>> 390 // surprisingly complicated quadric expression (this is easy to
>> 391 // appreciate once one understands that there may be up to
>> 392 // four lines normal to the ellipse intersecting any point). To
>> 393 // solve it exactly would be rather time consuming. This method,
>> 394 // however, is supposed to be a quick check, and is allowed to be an
>> 395 // underestimate.
>> 396 //
>> 397 // So, I will use the following underestimate of the distance
>> 398 // from an outside point to an ellipse. First: find the intersection "A"
>> 399 // of the line from the origin to the point with the ellipse.
>> 400 // Find the line passing through "A" and tangent to the ellipse
>> 401 // at A. The distance of the point p from the ellipse will be approximated
>> 402 // as the distance to this line.
>> 403 //
433 G4double G4EllipticalTube::DistanceToIn( const 404 G4double G4EllipticalTube::DistanceToIn( const G4ThreeVector& p ) const
434 { 405 {
435 // safety distance to bounding box << 406 static const G4double halfTol = 0.5*kCarTolerance;
436 G4double distX = std::abs(p.x()) - fDx; << 407
437 G4double distY = std::abs(p.y()) - fDy; << 408 if (CheckXY( p.x(), p.y(), +halfTol ) < 1.0) {
438 G4double distZ = std::abs(p.z()) - fDz; << 409 //
439 G4double distB = std::max(std::max(distX, di << 410 // We are inside or on the surface of the
440 // return (distB < 0) ? 0 : distB; << 411 // elliptical cross section in x/y. Check z
441 << 412 //
442 // safety distance to lateral surface << 413 if (p.z() < -dz-halfTol)
443 G4double x = p.x() * fSx; << 414 return -p.z()-dz;
444 G4double y = p.y() * fSy; << 415 else if (p.z() > dz+halfTol)
445 G4double distR = std::sqrt(x * x + y * y) - << 416 return p.z()-dz;
446 << 417 else
447 // return SafetyToIn << 418 return 0; // On any surface here (or inside)
448 G4double dist = std::max(distB, distR); << 419 }
449 return (dist < 0) ? 0 : dist; << 420
450 } << 421 //
451 << 422 // Find point on ellipse
452 ////////////////////////////////////////////// << 423 //
453 // << 424 G4double qnorm = CheckXY( p.x(), p.y() );
454 // Calculate distance to shape from inside and << 425 if (qnorm < DBL_MIN) return 0; // This should never happen
455 // at exit point, if required << 426
456 // - when leaving the surface, return 0 << 427 G4double q = 1.0/sqrt(qnorm);
457 << 428
458 G4double G4EllipticalTube::DistanceToOut( cons << 429 G4double xe = q*p.x(), ye = q*p.y();
459 cons << 430
460 cons << 431 //
461 << 432 // Get tangent to ellipse
462 << 433 //
463 { << 434 G4double tx = -ye*dx*dx, ty = +xe*dy*dy;
464 // Check if point flying away relative to Z << 435 G4double tnorm = sqrt( tx*tx + ty*ty );
465 // << 436
466 G4double pz = p.z(); << 437 //
467 G4double vz = v.z(); << 438 // Calculate distance
468 G4double distZ = std::abs(pz) - fDz; << 439 //
469 if (distZ >= -halfTolerance && pz * vz > 0) << 440 G4double distR = ( (p.x()-xe)*ty - (p.y()-ye)*tx )/tnorm;
470 { << 441
471 if (calcNorm) << 442 //
472 { << 443 // Add the result in quadrature if we are, in addition,
473 *validNorm = true; << 444 // outside the z bounds of the shape
474 n->set(0, 0, (pz < 0) ? -1. : 1.); << 445 //
475 } << 446 // We could save some time by returning the maximum rather
476 return 0.; << 447 // than the quadrature sum
477 } << 448 //
478 G4double tzmax = (vz == 0) ? DBL_MAX : (std: << 449 if (p.z() < -dz)
479 << 450 return sqrt( (p.z()+dz)*(p.z()+dz) + distR*distR );
480 // Scale elliptical tube to cylinder << 451 else if (p.z() > dz)
481 // << 452 return sqrt( (p.z()-dz)*(p.z()-dz) + distR*distR );
482 G4double px = p.x() * fSx; << 453
483 G4double py = p.y() * fSy; << 454 return distR;
484 G4double vx = v.x() * fSx; << 455 }
485 G4double vy = v.y() * fSy; << 456
486 << 457
487 // Check if point is flying away relative to << 458 //
488 // << 459 // DistanceToOut(p,v)
489 G4double rr = px * px + py * py; << 460 //
490 G4double B = px * vx + py * vy; << 461 // This method can be somewhat complicated for a general shape.
491 G4double distR = fQ1 * rr - fQ2; << 462 // For a convex one, like this, there are several simplifications,
492 if (distR >= -halfTolerance && B > 0.) << 463 // the most important of which is that one can treat the surfaces
493 { << 464 // as infinite in extent when deciding if the p is on the surface.
494 if (calcNorm) << 465 //
495 { << 466 G4double G4EllipticalTube::DistanceToOut( const G4ThreeVector& p,const G4ThreeVector& v,
496 *validNorm = true; << 467 const G4bool calcNorm,
497 *n = G4ThreeVector(px * fDDy, py * fDDx, << 468 G4bool *validNorm,G4ThreeVector *norm ) const
498 } << 469 {
499 return 0.; << 470 static const G4double halfTol = 0.5*kCarTolerance;
500 } << 471
501 << 472 //
502 // Just in case check if point is outside, n << 473 // Our normal is always valid
503 // << 474 //
504 if (std::max(distZ, distR) > halfTolerance) << 475 if (calcNorm) *validNorm = true;
505 { << 476
506 #ifdef G4SPECDEBUG << 477 G4double sBest = kInfinity;
507 std::ostringstream message; << 478 const G4ThreeVector *nBest=0;
508 G4long oldprc = message.precision(16); << 479
509 message << "Point p is outside (!?) of sol << 480 //
510 << GetName() << G4endl; << 481 // Might we intersect the -dz surface?
511 message << "Position: " << p << G4endl;; << 482 //
512 message << "Direction: " << v; << 483 if (v.z() < 0) {
513 G4cout.precision(oldprc); << 484 static const G4ThreeVector normHere(0.0,0.0,-1.0);
514 G4Exception("G4EllipticalTube::DistanceToO << 485 //
515 JustWarning, message ); << 486 // Yup. What distance?
516 DumpInfo(); << 487 //
517 #endif << 488 sBest = -(p.z()+dz)/v.z();
518 if (calcNorm) << 489
519 { << 490 //
520 *validNorm = true; << 491 // Are we on the surface? If so, return zero
521 *n = ApproxSurfaceNormal(p); << 492 //
522 } << 493 if (p.z() < -dz+halfTol) {
523 return 0.; << 494 if (calcNorm) *norm = normHere;
524 } << 495 return 0;
525 << 496 }
526 // Set coefficients of quadratic equation: A << 497 else
527 // << 498 nBest = &normHere;
528 G4double A = vx * vx + vy * vy; << 499 }
529 G4double C = rr - fR * fR; << 500
530 G4double D = B * B - A * C; << 501 //
531 << 502 // How about the +dz surface?
532 // Solve qudratic equation. There are two sp << 503 //
533 // 1) trajectory parallel to Z axis (A = 0 << 504 if (v.z() > 0) {
534 // 2) touch (D = 0) or no intersection (D << 505 static const G4ThreeVector normHere(0.0,0.0,+1.0);
535 // << 506 //
536 G4bool parallelToZ = (A < DBL_EPSILON || std << 507 // Yup. What distance?
537 if (parallelToZ) // 1) << 508 //
538 { << 509 G4double s = (dz-p.z())/v.z();
539 if (calcNorm) << 510
540 { << 511 //
541 *validNorm = true; << 512 // Are we on the surface? If so, return zero
542 n->set(0, 0, (vz < 0) ? -1. : 1.); << 513 //
543 } << 514 if (p.z() > +dz-halfTol) {
544 return tzmax; << 515 if (calcNorm) *norm = normHere;
545 } << 516 return 0;
546 if (D <= A * A * fScratch) // 2) << 517 }
547 { << 518
548 if (calcNorm) << 519 //
549 { << 520 // Best so far?
550 *validNorm = true; << 521 //
551 *n = G4ThreeVector(px * fDDy, py * fDDx, << 522 if (s < sBest) { sBest = s; nBest = &normHere; }
552 } << 523 }
553 return 0.; << 524
554 } << 525 //
555 << 526 // Check furthest intersection with ellipse
556 // Find roots of quadratic equation << 527 //
557 G4double tmp = -B - std::copysign(std::sqrt( << 528 G4double s[2];
558 G4double t1 = tmp / A; << 529 G4int n = IntersectXY( p, v, s );
559 G4double t2 = C / tmp; << 530
560 G4double trmax = std::max(t1, t2); << 531 if (n == 0) {
561 << 532 if (sBest == kInfinity)
562 // Return distance << 533 G4Exception( "G4EllipticalTube::DistanceToOut - Point is outside" );
563 G4double tmax = std::min(tzmax, trmax); << 534
564 << 535 if (calcNorm) *norm = *nBest;
565 // Set normal, if required, and return dista << 536 return sBest;
566 // << 537 }
567 if (calcNorm) << 538 else if (s[n-1] > sBest) {
568 { << 539 if (calcNorm) *norm = *nBest;
569 *validNorm = true; << 540 return sBest;
570 G4ThreeVector pnew = p + tmax * v; << 541 }
571 if (tmax == tzmax) << 542 sBest = s[n-1];
572 n->set(0, 0, (pnew.z() < 0) ? -1. : 1.); << 543
573 else << 544 //
574 *n = G4ThreeVector(pnew.x() * fDDy, pnew << 545 // Intersection with ellipse. Get normal at intersection point.
575 } << 546 //
576 return tmax; << 547 if (calcNorm) {
577 } << 548 G4ThreeVector ip = p + sBest*v;
578 << 549 *norm = G4ThreeVector( ip.x()*dy*dy, ip.y()*dx*dx, 0.0 ).unit();
579 ////////////////////////////////////////////// << 550 }
580 // << 551
581 // Estimate distance to the surface from insid << 552 //
582 // returns 0 if point is outside << 553 // Do we start on the surface?
>> 554 //
>> 555 if (CheckXY( p.x(), p.y(), -halfTol ) > 1.0) {
>> 556 //
>> 557 // Well, yes, but are we traveling outwards at this point?
>> 558 //
>> 559 if (p.x()*dy*dy*v.x() + p.y()*dx*dx*v.y() > 0) return 0;
>> 560 }
>> 561
>> 562 return sBest;
>> 563 }
>> 564
>> 565
>> 566 //
>> 567 // DistanceToOut(p)
>> 568 //
>> 569 // See DistanceToIn(p) for notes on the distance from a point
>> 570 // to an ellipse in two dimensions.
>> 571 //
>> 572 // The approximation used here for a point inside the ellipse
>> 573 // is to find the intersection with the ellipse of the lines
>> 574 // through the point and parallel to the x and y axes. The
>> 575 // distance of the point from the line connecting the two
>> 576 // intersecting points is then used.
583 // 577 //
584 <<
585 G4double G4EllipticalTube::DistanceToOut( cons 578 G4double G4EllipticalTube::DistanceToOut( const G4ThreeVector& p ) const
586 { 579 {
587 #ifdef G4SPECDEBUG << 580 static const G4double halfTol = 0.5*kCarTolerance;
588 if( Inside(p) == kOutside ) << 581
589 { << 582 //
590 std::ostringstream message; << 583 // We need to calculate the distances to all surfaces,
591 G4long oldprc = message.precision(16); << 584 // and then return the smallest
592 message << "Point p is outside (!?) of sol << 585 //
593 << "Position:\n" << 586 // Check -dz and +dz surface
594 << " p.x() = " << p.x()/mm << " << 587 //
595 << " p.y() = " << p.y()/mm << " << 588 G4double sBest = dz - fabs(p.z());
596 << " p.z() = " << p.z()/mm << " << 589 if (sBest < halfTol) return 0;
597 message.precision(oldprc) ; << 590
598 G4Exception("G4ElliptocalTube::DistanceToO << 591 //
599 JustWarning, message); << 592 // Check elliptical surface: find intersection of
600 DumpInfo(); << 593 // line through p and parallel to x axis
601 } << 594 //
602 #endif << 595 G4double radical = 1.0 - p.y()*p.y()/dy/dy;
603 // safety distance to Z-bases << 596 if (radical < +DBL_MIN) return 0;
604 G4double distZ = fDz - std::abs(p.z()); << 597
605 << 598 G4double xi = dx*sqrt( radical );
606 // safety distance lateral surface << 599 if (p.x() < 0) xi = -xi;
607 G4double x = p.x() * fSx; << 600
608 G4double y = p.y() * fSy; << 601 //
609 G4double distR = fR - std::sqrt(x * x + y * << 602 // Do the same with y axis
610 << 603 //
611 // return SafetyToOut << 604 radical = 1.0 - p.x()*p.x()/dx/dx;
612 G4double dist = std::min(distZ, distR); << 605 if (radical < +DBL_MIN) return 0;
613 return (dist < 0) ? 0 : dist; << 606
>> 607 G4double yi = dy*sqrt( radical );
>> 608 if (p.y() < 0) yi = -yi;
>> 609
>> 610 //
>> 611 // Get distance from p to the line connecting
>> 612 // these two points
>> 613 //
>> 614 G4double xdi = p.x() - xi,
>> 615 ydi = yi - p.y();
>> 616
>> 617 G4double normi = sqrt( xdi*xdi + ydi*ydi );
>> 618 if (normi < halfTol) return 0;
>> 619 xdi /= normi;
>> 620 ydi /= normi;
>> 621
>> 622 G4double s = 0.5*(xdi*(p.y()-yi) - ydi*(p.x()-xi));
>> 623 if (xi*yi < 0) s = -s;
>> 624
>> 625 if (s < sBest) sBest = s;
>> 626
>> 627 //
>> 628 // Return best answer
>> 629 //
>> 630 return sBest < halfTol ? 0 : sBest;
614 } 631 }
615 632
616 ////////////////////////////////////////////// <<
617 // <<
618 // GetEntityType <<
619 633
620 G4GeometryType G4EllipticalTube::GetEntityType <<
621 { <<
622 return {"G4EllipticalTube"}; <<
623 } <<
624 <<
625 ////////////////////////////////////////////// <<
626 // 634 //
627 // Make a clone of the object << 635 // CreatePolyhedron
628 <<
629 G4VSolid* G4EllipticalTube::Clone() const <<
630 { <<
631 return new G4EllipticalTube(*this); <<
632 } <<
633 <<
634 ////////////////////////////////////////////// <<
635 // 636 //
636 // Return volume << 637 G4Polyhedron* G4EllipticalTube::CreatePolyhedron() const
637 <<
638 G4double G4EllipticalTube::GetCubicVolume() <<
639 { 638 {
640 if (fCubicVolume == 0.) << 639 // create cylinder with radius=1...
641 { << 640 G4Polyhedron* eTube = new G4PolyhedronTube(0.,1.,dz);
642 fCubicVolume = twopi * fDx * fDy * fDz; << 641 // apply non-uniform scaling...
643 } << 642 eTube->Transform(G4Scale3D(dx,dy,1.));
644 return fCubicVolume; << 643 return eTube;
645 } 644 }
646 645
647 ////////////////////////////////////////////// <<
648 // <<
649 // Return cached surface area <<
650 646
651 G4double G4EllipticalTube::GetCachedSurfaceAre <<
652 { <<
653 G4ThreadLocalStatic G4double cached_Dx = 0; <<
654 G4ThreadLocalStatic G4double cached_Dy = 0; <<
655 G4ThreadLocalStatic G4double cached_Dz = 0; <<
656 G4ThreadLocalStatic G4double cached_area = 0 <<
657 if (cached_Dx != fDx || cached_Dy != fDy || <<
658 { <<
659 cached_Dx = fDx; <<
660 cached_Dy = fDy; <<
661 cached_Dz = fDz; <<
662 cached_area = 2.*(pi*fDx*fDy + G4GeomTools <<
663 } <<
664 return cached_area; <<
665 } <<
666 <<
667 ////////////////////////////////////////////// <<
668 // 647 //
669 // Return surface area << 648 // DescribeYourselfTo
670 <<
671 G4double G4EllipticalTube::GetSurfaceArea() <<
672 { <<
673 if(fSurfaceArea == 0.) <<
674 { <<
675 fSurfaceArea = GetCachedSurfaceArea(); <<
676 } <<
677 return fSurfaceArea; <<
678 } <<
679 <<
680 ////////////////////////////////////////////// <<
681 // 649 //
682 // Stream object contents to output stream << 650 void G4EllipticalTube::DescribeYourselfTo( G4VGraphicsScene& scene ) const
683 <<
684 std::ostream& G4EllipticalTube::StreamInfo(std <<
685 { 651 {
686 G4long oldprc = os.precision(16); << 652 scene.AddThis (*this);
687 os << "------------------------------------- <<
688 << " *** Dump for solid - " << GetName <<
689 << " ================================= <<
690 << " Solid type: G4EllipticalTube\n" <<
691 << " Parameters: \n" <<
692 << " length Z: " << fDz/mm << " mm \n" <<
693 << " lateral surface equation: \n" <<
694 << " (X / " << fDx << ")^2 + (Y / " <<
695 << "------------------------------------- <<
696 os.precision(oldprc); <<
697 <<
698 return os; <<
699 } 653 }
700 654
701 655
702 ////////////////////////////////////////////// <<
703 // 656 //
704 // Pick up a random point on the surface << 657 // GetExtent
705 <<
706 G4ThreeVector G4EllipticalTube::GetPointOnSurf <<
707 { <<
708 // Select surface (0 - base at -Z, 1 - base <<
709 // <<
710 G4double sbase = pi * fDx * fDy; <<
711 G4double ssurf = GetCachedSurfaceArea(); <<
712 G4double select = ssurf * G4UniformRand(); <<
713 <<
714 G4int k = 0; <<
715 if (select > sbase) k = 1; <<
716 if (select > 2. * sbase) k = 2; <<
717 <<
718 // Pick random point on selected surface (re <<
719 // <<
720 G4ThreeVector p; <<
721 switch (k) { <<
722 case 0: // base at -Z <<
723 { <<
724 G4TwoVector rho = G4RandomPointInEllipse <<
725 p.set(rho.x(), rho.y(), -fDz); <<
726 break; <<
727 } <<
728 case 1: // base at +Z <<
729 { <<
730 G4TwoVector rho = G4RandomPointInEllipse <<
731 p.set(rho.x(), rho.y(), fDz); <<
732 break; <<
733 } <<
734 case 2: // lateral surface <<
735 { <<
736 G4TwoVector rho = G4RandomPointOnEllipse <<
737 p.set(rho.x(), rho.y(), (2. * G4UniformR <<
738 break; <<
739 } <<
740 } <<
741 return p; <<
742 } <<
743 <<
744 <<
745 ////////////////////////////////////////////// <<
746 // 658 //
747 // CreatePolyhedron << 659 G4VisExtent G4EllipticalTube::GetExtent() const
748 <<
749 G4Polyhedron* G4EllipticalTube::CreatePolyhedr <<
750 { 660 {
751 // create cylinder with radius=1... << 661 return G4VisExtent( -dx, dx, -dy, dy, -dz, dz );
752 // <<
753 G4Polyhedron* eTube = new G4PolyhedronTube(0 <<
754 <<
755 // apply non-uniform scaling... <<
756 // <<
757 eTube->Transform(G4Scale3D(fDx, fDy, 1.)); <<
758 return eTube; <<
759 } 662 }
760 663
761 ////////////////////////////////////////////// <<
762 // <<
763 // GetPolyhedron <<
764 664
765 G4Polyhedron* G4EllipticalTube::GetPolyhedron <<
766 { <<
767 if (fpPolyhedron == nullptr || <<
768 fRebuildPolyhedron || <<
769 fpPolyhedron->GetNumberOfRotationStepsAt <<
770 fpPolyhedron->GetNumberOfRotationSteps() <<
771 { <<
772 G4AutoLock l(&polyhedronMutex); <<
773 delete fpPolyhedron; <<
774 fpPolyhedron = CreatePolyhedron(); <<
775 fRebuildPolyhedron = false; <<
776 l.unlock(); <<
777 } <<
778 return fpPolyhedron; <<
779 } <<
780 <<
781 ////////////////////////////////////////////// <<
782 // 665 //
783 // DescribeYourselfTo << 666 // IntersectXY
784 <<
785 void G4EllipticalTube::DescribeYourselfTo( G4V <<
786 { <<
787 scene.AddSolid (*this); <<
788 } <<
789 <<
790 ////////////////////////////////////////////// <<
791 // 667 //
792 // GetExtent << 668 // Decide if and where the x/y trajectory hits the elliptical cross
793 << 669 // section.
794 G4VisExtent G4EllipticalTube::GetExtent() cons << 670 //
795 { << 671 // Arguments:
796 return { -fDx, fDx, -fDy, fDy, -fDz, fDz }; << 672 // p - (in) Point on trajectory
>> 673 // v - (in) Vector along trajectory
>> 674 // s - (out) Up to two points of intersection, where the
>> 675 // intersection point is p + s*v, and if there are
>> 676 // two intersections, s[0] < s[1]. May be negative.
>> 677 // Returns:
>> 678 // The number of intersections. If 0, the trajectory misses. If 1, the
>> 679 // trajectory just grazes the surface.
>> 680 //
>> 681 // Solution:
>> 682 // One needs to solve: ( (p.x + s*v.x)/dx )**2 + ( (p.y + s*v.y)/dy )**2 = 1
>> 683 //
>> 684 // The solution is quadratic: a*s**2 + b*s + c = 0
>> 685 //
>> 686 // a = (v.x/dx)**2 + (v.y/dy)**2
>> 687 // b = 2*p.x*v.x/dx**2 + 2*p.y*v.y/dy**2
>> 688 // c = (p.x/dx)**2 + (p.y/dy)**2 - 1
>> 689 //
>> 690 G4int G4EllipticalTube::IntersectXY( const G4ThreeVector &p,
>> 691 const G4ThreeVector &v, G4double s[2] ) const
>> 692 {
>> 693 G4double px = p.x(), py = p.y();
>> 694 G4double vx = v.x(), vy = v.y();
>> 695
>> 696 G4double a = (vx/dx)*(vx/dx) + (vy/dy)*(vy/dy);
>> 697 G4double b = 2.0*( px*vx/dx/dx + py*vy/dy/dy );
>> 698 G4double c = (px/dx)*(px/dx) + (py/dy)*(py/dy) - 1.0;
>> 699
>> 700 if (a < DBL_MIN) return 0; // Trajectory parallel to z axis
>> 701
>> 702 G4double radical = b*b - 4*a*c;
>> 703
>> 704 if (radical < -DBL_MIN) return 0; // No solution
>> 705
>> 706 if (radical < DBL_MIN) {
>> 707 //
>> 708 // Grazes surface
>> 709 //
>> 710 s[0] = -b/a/2.0;
>> 711 return 1;
>> 712 }
>> 713
>> 714 radical = sqrt(radical);
>> 715
>> 716 G4double q = -0.5*( b + (b < 0 ? -radical : +radical) );
>> 717 G4double sa = q/a;
>> 718 G4double sb = c/q;
>> 719 if (sa < sb) { s[0] = sa; s[1] = sb; } else { s[0] = sb; s[1] = sa; }
>> 720 return 2;
797 } 721 }
798 <<
799 #endif // !defined(G4GEOM_USE_UELLIPTICALTUBE) <<
800 722