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