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