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