1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // $Id: G4Paraboloid.cc,v 1.5.2.1 2008/04/23 08:10:24 gcosmo Exp $ |
---|
27 | // GEANT4 tag $Name: geant4-09-01-patch-02 $ |
---|
28 | // |
---|
29 | // class G4Paraboloid |
---|
30 | // |
---|
31 | // Implementation for G4Paraboloid class |
---|
32 | // |
---|
33 | // History: |
---|
34 | // |
---|
35 | // Author : Lukas Lindroos (CERN), July 2007 |
---|
36 | // Revised: Tatiana Nikitina (CERN) |
---|
37 | // -------------------------------------------------------------------- |
---|
38 | |
---|
39 | #include "globals.hh" |
---|
40 | |
---|
41 | #include "G4Paraboloid.hh" |
---|
42 | |
---|
43 | #include "G4VoxelLimits.hh" |
---|
44 | #include "G4AffineTransform.hh" |
---|
45 | |
---|
46 | #include "meshdefs.hh" |
---|
47 | |
---|
48 | #include "Randomize.hh" |
---|
49 | |
---|
50 | #include "G4VGraphicsScene.hh" |
---|
51 | #include "G4Polyhedron.hh" |
---|
52 | #include "G4NURBS.hh" |
---|
53 | #include "G4NURBSbox.hh" |
---|
54 | #include "G4VisExtent.hh" |
---|
55 | |
---|
56 | using namespace CLHEP; |
---|
57 | |
---|
58 | /////////////////////////////////////////////////////////////////////////////// |
---|
59 | // |
---|
60 | // constructor - check parameters |
---|
61 | |
---|
62 | G4Paraboloid::G4Paraboloid(const G4String& pName, |
---|
63 | G4double pDz, |
---|
64 | G4double pR1, |
---|
65 | G4double pR2) |
---|
66 | : G4VSolid(pName),fpPolyhedron(0), fCubicVolume(0.) |
---|
67 | |
---|
68 | { |
---|
69 | if(pDz > 0. && pR2 > pR1 && pR1 >= 0.) |
---|
70 | { |
---|
71 | r1 = pR1; |
---|
72 | r2 = pR2; |
---|
73 | dz = pDz; |
---|
74 | } |
---|
75 | else |
---|
76 | { |
---|
77 | G4cerr << "Error - G4Paraboloid::G4Paraboloid(): " << GetName() << G4endl |
---|
78 | << "Z half-length must be larger than zero or R1>=R2 " << G4endl; |
---|
79 | G4Exception("G4Paraboloid::G4Paraboloid()", "InvalidSetup", |
---|
80 | FatalException, |
---|
81 | "Invalid dimensions. Negative Input Values or R1>=R2."); |
---|
82 | } |
---|
83 | |
---|
84 | // r1^2 = k1 * (-dz) + k2 |
---|
85 | // r2^2 = k1 * ( dz) + k2 |
---|
86 | // => r1^2 + r2^2 = k2 + k2 => k2 = (r2^2 + r1^2) / 2 |
---|
87 | // and r2^2 - r1^2 = k1 * dz - k1 * (-dz) => k1 = (r2^2 - r1^2) / 2 / dz |
---|
88 | |
---|
89 | k1 = (r2 * r2 - r1 * r1) / 2 / dz; |
---|
90 | k2 = (r2 * r2 + r1 * r1) / 2; |
---|
91 | |
---|
92 | fSurfaceArea = 0.; |
---|
93 | } |
---|
94 | |
---|
95 | /////////////////////////////////////////////////////////////////////////////// |
---|
96 | // |
---|
97 | // Fake default constructor - sets only member data and allocates memory |
---|
98 | // for usage restricted to object persistency. |
---|
99 | // |
---|
100 | G4Paraboloid::G4Paraboloid( __void__& a ) |
---|
101 | : G4VSolid(a), fpPolyhedron(0), fCubicVolume(0.) |
---|
102 | { |
---|
103 | fSurfaceArea = 0.; |
---|
104 | } |
---|
105 | |
---|
106 | /////////////////////////////////////////////////////////////////////////////// |
---|
107 | // |
---|
108 | // Destructor |
---|
109 | |
---|
110 | G4Paraboloid::~G4Paraboloid() |
---|
111 | { |
---|
112 | } |
---|
113 | |
---|
114 | ///////////////////////////////////////////////////////////////////////// |
---|
115 | // |
---|
116 | // Dispatch to parameterisation for replication mechanism dimension |
---|
117 | // computation & modification. |
---|
118 | |
---|
119 | //void ComputeDimensions( G4VPVParamerisation p, |
---|
120 | // const G4Int n, |
---|
121 | // const G4VPhysicalVolume* pRep ) |
---|
122 | //{ |
---|
123 | // p->ComputeDimensions(*this,n,pRep) ; |
---|
124 | //} |
---|
125 | |
---|
126 | |
---|
127 | /////////////////////////////////////////////////////////////////////////////// |
---|
128 | // |
---|
129 | // Calculate extent under transform and specified limit |
---|
130 | |
---|
131 | G4bool |
---|
132 | G4Paraboloid::CalculateExtent(const EAxis pAxis, |
---|
133 | const G4VoxelLimits& pVoxelLimit, |
---|
134 | const G4AffineTransform& pTransform, |
---|
135 | G4double& pMin, G4double& pMax) const |
---|
136 | { |
---|
137 | G4double xMin = -r2 + pTransform.NetTranslation().x(), |
---|
138 | xMax = r2 + pTransform.NetTranslation().x(), |
---|
139 | yMin = -r2 + pTransform.NetTranslation().y(), |
---|
140 | yMax = r2 + pTransform.NetTranslation().y(), |
---|
141 | zMin = -dz + pTransform.NetTranslation().z(), |
---|
142 | zMax = dz + pTransform.NetTranslation().z(); |
---|
143 | |
---|
144 | if(!pTransform.IsRotated() |
---|
145 | || pTransform.NetRotation()(G4ThreeVector(0, 0, 1)) == G4ThreeVector(0, 0, 1)) |
---|
146 | { |
---|
147 | if(pVoxelLimit.IsXLimited()) |
---|
148 | { |
---|
149 | if(pVoxelLimit.GetMaxXExtent() < xMin - 0.5 * kCarTolerance |
---|
150 | || pVoxelLimit.GetMinXExtent() > xMax + 0.5 * kCarTolerance) |
---|
151 | { |
---|
152 | return false; |
---|
153 | } |
---|
154 | else |
---|
155 | { |
---|
156 | if(pVoxelLimit.GetMinXExtent() > xMin) |
---|
157 | { |
---|
158 | xMin = pVoxelLimit.GetMinXExtent(); |
---|
159 | } |
---|
160 | if(pVoxelLimit.GetMaxXExtent() < xMax) |
---|
161 | { |
---|
162 | xMax = pVoxelLimit.GetMaxXExtent(); |
---|
163 | } |
---|
164 | } |
---|
165 | } |
---|
166 | if(pVoxelLimit.IsYLimited()) |
---|
167 | { |
---|
168 | if(pVoxelLimit.GetMaxYExtent() < yMin - 0.5 * kCarTolerance |
---|
169 | || pVoxelLimit.GetMinYExtent() > yMax + 0.5 * kCarTolerance) |
---|
170 | { |
---|
171 | return false; |
---|
172 | } |
---|
173 | else |
---|
174 | { |
---|
175 | if(pVoxelLimit.GetMinYExtent() > yMin) |
---|
176 | { |
---|
177 | yMin = pVoxelLimit.GetMinYExtent(); |
---|
178 | } |
---|
179 | if(pVoxelLimit.GetMaxYExtent() < yMax) |
---|
180 | { |
---|
181 | yMax = pVoxelLimit.GetMaxYExtent(); |
---|
182 | } |
---|
183 | } |
---|
184 | } |
---|
185 | if(pVoxelLimit.IsZLimited()) |
---|
186 | { |
---|
187 | if(pVoxelLimit.GetMaxZExtent() < zMin - 0.5 * kCarTolerance |
---|
188 | || pVoxelLimit.GetMinZExtent() > zMax + 0.5 * kCarTolerance) |
---|
189 | { |
---|
190 | return false; |
---|
191 | } |
---|
192 | else |
---|
193 | { |
---|
194 | if(pVoxelLimit.GetMinZExtent() > zMin) |
---|
195 | { |
---|
196 | zMin = pVoxelLimit.GetMinZExtent(); |
---|
197 | } |
---|
198 | if(pVoxelLimit.GetMaxZExtent() < zMax) |
---|
199 | { |
---|
200 | zMax = pVoxelLimit.GetMaxZExtent(); |
---|
201 | } |
---|
202 | } |
---|
203 | } |
---|
204 | switch(pAxis) |
---|
205 | { |
---|
206 | case kXAxis: |
---|
207 | pMin = xMin; |
---|
208 | pMax = xMax; |
---|
209 | break; |
---|
210 | case kYAxis: |
---|
211 | pMin = yMin; |
---|
212 | pMax = yMax; |
---|
213 | break; |
---|
214 | case kZAxis: |
---|
215 | pMin = zMin; |
---|
216 | pMax = zMax; |
---|
217 | break; |
---|
218 | default: |
---|
219 | pMin = 0; |
---|
220 | pMax = 0; |
---|
221 | return false; |
---|
222 | } |
---|
223 | } |
---|
224 | else |
---|
225 | { |
---|
226 | G4bool existsAfterClip=true; |
---|
227 | |
---|
228 | // Calculate rotated vertex coordinates |
---|
229 | |
---|
230 | G4int noPolygonVertices=0; |
---|
231 | G4ThreeVectorList* vertices |
---|
232 | = CreateRotatedVertices(pTransform,noPolygonVertices); |
---|
233 | |
---|
234 | if(pAxis == kXAxis || pAxis == kYAxis || pAxis == kZAxis) |
---|
235 | { |
---|
236 | pMin = kInfinity; |
---|
237 | pMax = -kInfinity; |
---|
238 | |
---|
239 | for(G4ThreeVectorList::iterator it = vertices->begin(); |
---|
240 | it < vertices->end(); it++) |
---|
241 | { |
---|
242 | if(pMin > (*it)[pAxis]) pMin = (*it)[pAxis]; |
---|
243 | if((*it)[pAxis] < pVoxelLimit.GetMinExtent(pAxis)) |
---|
244 | { |
---|
245 | pMin = pVoxelLimit.GetMinExtent(pAxis); |
---|
246 | } |
---|
247 | if(pMax < (*it)[pAxis]) |
---|
248 | { |
---|
249 | pMax = (*it)[pAxis]; |
---|
250 | } |
---|
251 | if((*it)[pAxis] > pVoxelLimit.GetMaxExtent(pAxis)) |
---|
252 | { |
---|
253 | pMax = pVoxelLimit.GetMaxExtent(pAxis); |
---|
254 | } |
---|
255 | } |
---|
256 | |
---|
257 | if(pMin > pVoxelLimit.GetMaxExtent(pAxis) |
---|
258 | || pMax < pVoxelLimit.GetMinExtent(pAxis)) { existsAfterClip = false; } |
---|
259 | } |
---|
260 | else |
---|
261 | { |
---|
262 | pMin = 0; |
---|
263 | pMax = 0; |
---|
264 | existsAfterClip = false; |
---|
265 | } |
---|
266 | delete vertices; |
---|
267 | return existsAfterClip; |
---|
268 | } |
---|
269 | return true; |
---|
270 | } |
---|
271 | |
---|
272 | /////////////////////////////////////////////////////////////////////////////// |
---|
273 | // |
---|
274 | // Return whether point inside/outside/on surface |
---|
275 | |
---|
276 | EInside G4Paraboloid::Inside(const G4ThreeVector& p) const |
---|
277 | { |
---|
278 | // First check is the point is above or below the solid. |
---|
279 | // |
---|
280 | if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance) { return kOutside; } |
---|
281 | |
---|
282 | G4double rho2 = p.perp2(), |
---|
283 | rhoSurfTimesTol2 = (k1 * p.z() + k2) * sqr(kCarTolerance), |
---|
284 | A = rho2 - ((k1 *p.z() + k2) + 0.25 * kCarTolerance * kCarTolerance); |
---|
285 | |
---|
286 | if(A < 0 && sqr(A) > rhoSurfTimesTol2) |
---|
287 | { |
---|
288 | // Actually checking rho < radius of paraboloid at z = p.z(). |
---|
289 | // We're either inside or in lower/upper cutoff area. |
---|
290 | |
---|
291 | if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance) |
---|
292 | { |
---|
293 | // We're in the upper/lower cutoff area, sides have a paraboloid shape |
---|
294 | // maybe further checks should be made to make these nicer |
---|
295 | |
---|
296 | return kSurface; |
---|
297 | } |
---|
298 | else |
---|
299 | { |
---|
300 | return kInside; |
---|
301 | } |
---|
302 | } |
---|
303 | else if(A <= 0 || sqr(A) < rhoSurfTimesTol2) |
---|
304 | { |
---|
305 | // We're in the parabolic surface. |
---|
306 | |
---|
307 | return kSurface; |
---|
308 | } |
---|
309 | else |
---|
310 | { |
---|
311 | return kOutside; |
---|
312 | } |
---|
313 | } |
---|
314 | |
---|
315 | /////////////////////////////////////////////////////////////////////////////// |
---|
316 | // |
---|
317 | |
---|
318 | G4ThreeVector G4Paraboloid::SurfaceNormal( const G4ThreeVector& p) const |
---|
319 | { |
---|
320 | G4ThreeVector n(0, 0, 0); |
---|
321 | if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance) |
---|
322 | { |
---|
323 | // If above or below just return normal vector for the cutoff plane. |
---|
324 | |
---|
325 | n = G4ThreeVector(0, 0, p.z()/std::fabs(p.z())); |
---|
326 | } |
---|
327 | else if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance) |
---|
328 | { |
---|
329 | // This means we're somewhere in the plane z = dz or z = -dz. |
---|
330 | // (As far as the program is concerned anyway. |
---|
331 | |
---|
332 | if(p.z() < 0) // Are we in upper or lower plane? |
---|
333 | { |
---|
334 | if(p.perp2() > sqr(r1 + 0.5 * kCarTolerance)) |
---|
335 | { |
---|
336 | n = G4ThreeVector(p.x(), p.y(), -k1 / 2).unit(); |
---|
337 | } |
---|
338 | else if(r1 < 0.5 * kCarTolerance |
---|
339 | || p.perp2() > sqr(r1 - 0.5 * kCarTolerance)) |
---|
340 | { |
---|
341 | n = G4ThreeVector(p.x(), p.y(), 0.).unit() |
---|
342 | + G4ThreeVector(0., 0., -1.).unit(); |
---|
343 | n = n.unit(); |
---|
344 | } |
---|
345 | else |
---|
346 | { |
---|
347 | n = G4ThreeVector(0., 0., -1.); |
---|
348 | } |
---|
349 | } |
---|
350 | else |
---|
351 | { |
---|
352 | if(p.perp2() > sqr(r2 + 0.5 * kCarTolerance)) |
---|
353 | { |
---|
354 | n = G4ThreeVector(p.x(), p.y(), 0.).unit(); |
---|
355 | } |
---|
356 | else if(r2 < 0.5 * kCarTolerance |
---|
357 | || p.perp2() > sqr(r2 - 0.5 * kCarTolerance)) |
---|
358 | { |
---|
359 | n = G4ThreeVector(p.x(), p.y(), 0.).unit() |
---|
360 | + G4ThreeVector(0., 0., 1.).unit(); |
---|
361 | n = n.unit(); |
---|
362 | } |
---|
363 | else |
---|
364 | { |
---|
365 | n = G4ThreeVector(0., 0., 1.); |
---|
366 | } |
---|
367 | } |
---|
368 | } |
---|
369 | else |
---|
370 | { |
---|
371 | G4double rho2 = p.perp2(); |
---|
372 | G4double rhoSurfTimesTol2 = (k1 * p.z() + k2) * sqr(kCarTolerance); |
---|
373 | G4double A = rho2 - ((k1 *p.z() + k2) |
---|
374 | + 0.25 * kCarTolerance * kCarTolerance); |
---|
375 | |
---|
376 | if(A < 0 && sqr(A) > rhoSurfTimesTol2) |
---|
377 | { |
---|
378 | // Actually checking rho < radius of paraboloid at z = p.z(). |
---|
379 | // We're inside. |
---|
380 | |
---|
381 | if(p.mag2() != 0) { n = p.unit(); } |
---|
382 | } |
---|
383 | else if(A <= 0 || sqr(A) < rhoSurfTimesTol2) |
---|
384 | { |
---|
385 | // We're in the parabolic surface. |
---|
386 | |
---|
387 | n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit(); |
---|
388 | } |
---|
389 | else |
---|
390 | { |
---|
391 | n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit(); |
---|
392 | } |
---|
393 | } |
---|
394 | |
---|
395 | if(n.mag2() == 0) |
---|
396 | { |
---|
397 | G4cerr << "WARNING - G4Paraboloid::SurfaceNormal(p)" << G4endl |
---|
398 | << " p = " << 1 / mm * p << " mm" << G4endl; |
---|
399 | G4Exception("G4Paraboloid::SurfaceNormal(p)", "Notification", |
---|
400 | JustWarning, "No normal defined for this point p."); |
---|
401 | } |
---|
402 | return n; |
---|
403 | } |
---|
404 | |
---|
405 | /////////////////////////////////////////////////////////////////////////////// |
---|
406 | // |
---|
407 | // Calculate distance to shape from outside, along normalised vector |
---|
408 | // - return kInfinity if no intersection |
---|
409 | // |
---|
410 | |
---|
411 | G4double G4Paraboloid::DistanceToIn( const G4ThreeVector& p, |
---|
412 | const G4ThreeVector& v ) const |
---|
413 | { |
---|
414 | G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2); |
---|
415 | G4double tol2 = kCarTolerance*kCarTolerance; |
---|
416 | G4double tolh = 0.5*kCarTolerance; |
---|
417 | |
---|
418 | if(r2 && p.z() > - tolh + dz) |
---|
419 | { |
---|
420 | // If the points is above check for intersection with upper edge. |
---|
421 | |
---|
422 | if(v.z() < 0) |
---|
423 | { |
---|
424 | G4double intersection = (dz - p.z()) / v.z(); // With plane z = dz. |
---|
425 | if(sqr(p.x() + v.x()*intersection) |
---|
426 | + sqr(p.y() + v.y()*intersection) < sqr(r2 + 0.5 * kCarTolerance)) |
---|
427 | { |
---|
428 | if(p.z() < tolh + dz) |
---|
429 | { return 0; } |
---|
430 | else |
---|
431 | { return intersection; } |
---|
432 | } |
---|
433 | } |
---|
434 | else // Direction away, no posibility of intersection |
---|
435 | { return kInfinity; } |
---|
436 | } |
---|
437 | else if(r1 && p.z() < tolh - dz) |
---|
438 | { |
---|
439 | // If the points is belove check for intersection with lower edge. |
---|
440 | |
---|
441 | if(v.z() > 0) |
---|
442 | { |
---|
443 | G4double intersection = (-dz - p.z()) / v.z(); // With plane z = -dz. |
---|
444 | if(sqr(p.x() + v.x()*intersection) |
---|
445 | + sqr(p.y() + v.y()*intersection) < sqr(r1 + 0.5 * kCarTolerance)) |
---|
446 | { |
---|
447 | if(p.z() > -tolh - dz) |
---|
448 | { |
---|
449 | return 0; |
---|
450 | } |
---|
451 | else |
---|
452 | { |
---|
453 | return intersection; |
---|
454 | } |
---|
455 | } |
---|
456 | } |
---|
457 | else// Direction away, no posibility of intersection |
---|
458 | { return kInfinity; } |
---|
459 | } |
---|
460 | |
---|
461 | G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y(), |
---|
462 | vRho2 = v.perp2(), intersection, |
---|
463 | B = (k1 * p.z() + k2 - rho2) * vRho2; |
---|
464 | |
---|
465 | if ( ( (rho2 > paraRho2) && (sqr(rho2-paraRho2-0.25*tol2) > tol2*paraRho2) ) |
---|
466 | || (p.z() < - dz+kCarTolerance) |
---|
467 | || (p.z() > dz-kCarTolerance) ) // Make sure it's safely outside. |
---|
468 | { |
---|
469 | // Is there a problem with squaring rho twice? |
---|
470 | |
---|
471 | if(!vRho2) // Needs to be treated seperately. |
---|
472 | { |
---|
473 | intersection = ((rho2 - k2)/k1 - p.z())/v.z(); |
---|
474 | if(intersection < 0) { return kInfinity; } |
---|
475 | else if(std::fabs(p.z() + v.z() * intersection) <= dz) |
---|
476 | { |
---|
477 | return intersection; |
---|
478 | } |
---|
479 | else |
---|
480 | { |
---|
481 | return kInfinity; |
---|
482 | } |
---|
483 | } |
---|
484 | else if(A*A + B < 0) // No real intersections. |
---|
485 | { |
---|
486 | return kInfinity; |
---|
487 | } |
---|
488 | else |
---|
489 | { |
---|
490 | intersection = (A - std::sqrt(B + sqr(A))) / vRho2; |
---|
491 | if(intersection < 0) |
---|
492 | { |
---|
493 | return kInfinity; |
---|
494 | } |
---|
495 | else if(std::fabs(p.z() + intersection * v.z()) < dz + tolh) |
---|
496 | { |
---|
497 | return intersection; |
---|
498 | } |
---|
499 | else |
---|
500 | { |
---|
501 | return kInfinity; |
---|
502 | } |
---|
503 | } |
---|
504 | } |
---|
505 | else if(sqr(rho2 - paraRho2 - .25 * tol2) <= tol2 * paraRho2) |
---|
506 | { |
---|
507 | // If this is true we're somewhere in the border. |
---|
508 | |
---|
509 | G4ThreeVector normal(p.x(), p.y(), -k1/2); |
---|
510 | if(normal.dot(v) <= 0) |
---|
511 | { return 0; } |
---|
512 | else |
---|
513 | { return kInfinity; } |
---|
514 | } |
---|
515 | else |
---|
516 | { |
---|
517 | G4cerr << "WARNING - G4Paraboloid::DistanceToIn(p,v)" << G4endl |
---|
518 | << " p = " << p * (1/mm) << " mm" << G4endl |
---|
519 | << " v = " << v * (1/mm) << " mm" << G4endl; |
---|
520 | if(Inside(p) == kInside) |
---|
521 | { |
---|
522 | G4Exception("G4Paraboloid::DistanceToIn(p,v)", "Notification", |
---|
523 | JustWarning, "Point p is inside!"); |
---|
524 | } |
---|
525 | else |
---|
526 | { |
---|
527 | G4Exception("G4Paraboloid::DistanceToIn(p,v)", "Notification", |
---|
528 | JustWarning, "There's a bug in this function (apa)!"); |
---|
529 | } |
---|
530 | return 0; |
---|
531 | } |
---|
532 | } |
---|
533 | |
---|
534 | /////////////////////////////////////////////////////////////////////////////// |
---|
535 | // |
---|
536 | // Calculate distance (<= actual) to closest surface of shape from outside |
---|
537 | // - Return 0 if point inside |
---|
538 | |
---|
539 | G4double G4Paraboloid::DistanceToIn(const G4ThreeVector& p) const |
---|
540 | { |
---|
541 | G4double safe = 0; |
---|
542 | if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance) |
---|
543 | { |
---|
544 | // If we're below or above the paraboloid treat it as a cylinder with |
---|
545 | // radius r2. |
---|
546 | |
---|
547 | if(p.perp2() > sqr(r2 + 0.5 * kCarTolerance)) |
---|
548 | { |
---|
549 | // This algorithm is exact for now but contains 2 sqrt calculations. |
---|
550 | // Might be better to replace with approximated version |
---|
551 | |
---|
552 | G4double dRho = p.perp() - r2; |
---|
553 | safe = std::sqrt(sqr(dRho) + sqr(std::fabs(p.z()) - dz)); |
---|
554 | } |
---|
555 | else |
---|
556 | { |
---|
557 | safe = std::fabs(p.z()) - dz; |
---|
558 | } |
---|
559 | } |
---|
560 | else |
---|
561 | { |
---|
562 | G4double paraRho = std::sqrt(k1 * p.z() + k2); |
---|
563 | G4double rho = p.perp(); |
---|
564 | |
---|
565 | if(rho > paraRho + 0.5 * kCarTolerance) |
---|
566 | { |
---|
567 | // Should check the value of paraRho here, |
---|
568 | // for small values the following algorithm is bad. |
---|
569 | |
---|
570 | safe = k1 / 2 * (-paraRho + rho) / rho; |
---|
571 | if(safe < 0) { safe = 0; } |
---|
572 | } |
---|
573 | } |
---|
574 | return safe; |
---|
575 | } |
---|
576 | |
---|
577 | /////////////////////////////////////////////////////////////////////////////// |
---|
578 | // |
---|
579 | // Calculate distance to surface of shape from `inside' |
---|
580 | |
---|
581 | G4double G4Paraboloid::DistanceToOut(const G4ThreeVector& p, |
---|
582 | const G4ThreeVector& v, |
---|
583 | const G4bool calcNorm, |
---|
584 | G4bool *validNorm, |
---|
585 | G4ThreeVector *n ) const |
---|
586 | { |
---|
587 | G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2); |
---|
588 | G4double vRho2 = v.perp2(), intersection; |
---|
589 | G4double tol2 = kCarTolerance*kCarTolerance; |
---|
590 | G4double tolh = 0.5*kCarTolerance; |
---|
591 | |
---|
592 | if(calcNorm) { *validNorm = false; } |
---|
593 | |
---|
594 | // We have that the particle p follows the line x = p + s * v |
---|
595 | // meaning x = p.x() + s * v.x(), y = p.y() + s * v.y() and |
---|
596 | // z = p.z() + s * v.z() |
---|
597 | // The equation for all points on the surface (surface expanded for |
---|
598 | // to include all z) x^2 + y^2 = k1 * z + k2 => .. => |
---|
599 | // => s = (A +- std::sqrt(A^2 + B)) / vRho2 |
---|
600 | // where |
---|
601 | G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y(); |
---|
602 | // and |
---|
603 | G4double B = (-rho2 + paraRho2) * vRho2; |
---|
604 | |
---|
605 | if ( rho2 < paraRho2 && sqr(rho2 - paraRho2 - 0.25 * tol2) > tol2 * paraRho2 |
---|
606 | && std::fabs(p.z()) < dz - kCarTolerance) |
---|
607 | { |
---|
608 | // Make sure it's safely inside. |
---|
609 | |
---|
610 | if(v.z() > 0) |
---|
611 | { |
---|
612 | // It's heading upwards, check where it colides with the plane z = dz. |
---|
613 | // When it does, is that in the surface of the paraboloid. |
---|
614 | // z = p.z() + variable * v.z() for all points the particle can go. |
---|
615 | // => variable = (z - p.z()) / v.z() so intersection must be: |
---|
616 | |
---|
617 | intersection = (dz - p.z()) / v.z(); |
---|
618 | G4ThreeVector ip = p + intersection * v; // Point of intersection. |
---|
619 | |
---|
620 | if(ip.perp2() < sqr(r2 + kCarTolerance)) |
---|
621 | { |
---|
622 | if(calcNorm) |
---|
623 | { |
---|
624 | *n = G4ThreeVector(0, 0, 1); |
---|
625 | if(r2 < tolh || ip.perp2() > sqr(r2 - tolh)) |
---|
626 | { |
---|
627 | *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit(); |
---|
628 | *n = n->unit(); |
---|
629 | } |
---|
630 | *validNorm = true; |
---|
631 | } |
---|
632 | return intersection; |
---|
633 | } |
---|
634 | } |
---|
635 | else if(v.z() < 0) |
---|
636 | { |
---|
637 | // It's heading downwards, check were it colides with the plane z = -dz. |
---|
638 | // When it does, is that in the surface of the paraboloid. |
---|
639 | // z = p.z() + variable * v.z() for all points the particle can go. |
---|
640 | // => variable = (z - p.z()) / v.z() so intersection must be: |
---|
641 | |
---|
642 | intersection = (-dz - p.z()) / v.z(); |
---|
643 | G4ThreeVector ip = p + intersection * v; // Point of intersection. |
---|
644 | |
---|
645 | if(ip.perp2() < sqr(r1 + tolh)) |
---|
646 | { |
---|
647 | if(calcNorm) |
---|
648 | { |
---|
649 | *n = G4ThreeVector(0, 0, -1); |
---|
650 | if(r1 < tolh || ip.perp2() > sqr(r1 - tolh)) |
---|
651 | { |
---|
652 | *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit(); |
---|
653 | *n = n->unit(); |
---|
654 | } |
---|
655 | *validNorm = true; |
---|
656 | } |
---|
657 | return intersection; |
---|
658 | } |
---|
659 | } |
---|
660 | |
---|
661 | // Now check for collisions with paraboloid surface. |
---|
662 | |
---|
663 | if(vRho2 == 0) // Needs to be treated seperately. |
---|
664 | { |
---|
665 | intersection = ((rho2 - k2)/k1 - p.z())/v.z(); |
---|
666 | if(calcNorm) |
---|
667 | { |
---|
668 | G4ThreeVector intersectionP = p + v * intersection; |
---|
669 | *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2); |
---|
670 | *n = n->unit(); |
---|
671 | |
---|
672 | *validNorm = true; |
---|
673 | } |
---|
674 | return intersection; |
---|
675 | } |
---|
676 | else if( ((A <= 0) && (B >= sqr(A) * (sqr(vRho2) - 1))) || (A >= 0)) |
---|
677 | { |
---|
678 | intersection = (A + std::sqrt(B + sqr(A))) / vRho2; |
---|
679 | if(calcNorm) |
---|
680 | { |
---|
681 | G4ThreeVector intersectionP = p + v * intersection; |
---|
682 | *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2); |
---|
683 | *n = n->unit(); |
---|
684 | *validNorm = true; |
---|
685 | } |
---|
686 | return intersection; |
---|
687 | } |
---|
688 | G4cerr << "WARNING - G4Paraboloid::DistanceToOut(p,v,...)" << G4endl |
---|
689 | << " p = " << p << G4endl |
---|
690 | << " v = " << v << G4endl; |
---|
691 | G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "Notification", |
---|
692 | JustWarning, |
---|
693 | "There is no intersection between given line and solid!"); |
---|
694 | |
---|
695 | return kInfinity; |
---|
696 | } |
---|
697 | else if ( (rho2 < paraRho2 + kCarTolerance |
---|
698 | || sqr(rho2 - paraRho2 - 0.25 * tol2) < tol2 * paraRho2 ) |
---|
699 | && std::fabs(p.z()) < dz + tolh) |
---|
700 | { |
---|
701 | // If this is true we're somewhere in the border. |
---|
702 | |
---|
703 | G4ThreeVector normal = G4ThreeVector (p.x(), p.y(), -k1/2); |
---|
704 | |
---|
705 | if(std::fabs(p.z()) > dz - tolh) |
---|
706 | { |
---|
707 | // We're in the lower or upper edge |
---|
708 | if( ((v.z() > 0) && (p.z() > 0)) || ((v.z() < 0) && (p.z() < 0)) ) |
---|
709 | // If we're headig out of the object that is treated here |
---|
710 | { |
---|
711 | if(calcNorm) |
---|
712 | { |
---|
713 | *validNorm = true; |
---|
714 | if(p.z() > 0) |
---|
715 | { *n = G4ThreeVector(0, 0, 1); } |
---|
716 | else |
---|
717 | { *n = G4ThreeVector(0, 0, -1); } |
---|
718 | } |
---|
719 | return 0; |
---|
720 | } |
---|
721 | |
---|
722 | if(v.z() == 0) |
---|
723 | { |
---|
724 | // Case where we're moving inside the surface needs to be |
---|
725 | // treated separately. |
---|
726 | // Distance until it goes out through a side is returned. |
---|
727 | |
---|
728 | G4double r = (p.z() > 0)? r2 : r1; |
---|
729 | G4double pDotV = p.dot(v); |
---|
730 | G4double A = vRho2 * ( sqr(r) - sqr(p.x()) - sqr(p.y())); |
---|
731 | intersection = (-pDotV + std::sqrt(A + sqr(pDotV))) / vRho2; |
---|
732 | |
---|
733 | if(calcNorm) |
---|
734 | { |
---|
735 | *validNorm = true; |
---|
736 | |
---|
737 | *n = (G4ThreeVector(0, 0, p.z()/std::fabs(p.z())) |
---|
738 | + G4ThreeVector(p.x() + v.x() * intersection, p.y() + v.y() |
---|
739 | * intersection, -k1/2).unit()).unit(); |
---|
740 | } |
---|
741 | return intersection; |
---|
742 | } |
---|
743 | } |
---|
744 | else if(normal.dot(v) >= 0) |
---|
745 | { |
---|
746 | if(calcNorm) |
---|
747 | { |
---|
748 | *validNorm = true; |
---|
749 | *n = normal.unit(); |
---|
750 | } |
---|
751 | return 0; |
---|
752 | } |
---|
753 | |
---|
754 | if(v.z() > 0) |
---|
755 | { |
---|
756 | // Check for collision with upper edge. |
---|
757 | |
---|
758 | intersection = (dz - p.z()) / v.z(); |
---|
759 | G4ThreeVector ip = p + intersection * v; |
---|
760 | |
---|
761 | if(ip.perp2() < sqr(r2 - tolh)) |
---|
762 | { |
---|
763 | if(calcNorm) |
---|
764 | { |
---|
765 | *validNorm = true; |
---|
766 | *n = G4ThreeVector(0, 0, 1); |
---|
767 | } |
---|
768 | return intersection; |
---|
769 | } |
---|
770 | else if(ip.perp2() < sqr(r2 + tolh)) |
---|
771 | { |
---|
772 | if(calcNorm) |
---|
773 | { |
---|
774 | *validNorm = true; |
---|
775 | *n = G4ThreeVector(0, 0, 1) |
---|
776 | + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit(); |
---|
777 | *n = n->unit(); |
---|
778 | } |
---|
779 | return intersection; |
---|
780 | } |
---|
781 | } |
---|
782 | if(r1 && v.z() < 0) |
---|
783 | { |
---|
784 | // Check for collision with lower edge. |
---|
785 | |
---|
786 | intersection = (-dz - p.z()) / v.z(); |
---|
787 | G4ThreeVector ip = p + intersection * v; |
---|
788 | |
---|
789 | if(ip.perp2() < sqr(r1 - tolh)) |
---|
790 | { |
---|
791 | if(calcNorm) |
---|
792 | { |
---|
793 | *validNorm = true; |
---|
794 | *n = G4ThreeVector(0, 0, -1); |
---|
795 | } |
---|
796 | return intersection; |
---|
797 | } |
---|
798 | else if(ip.perp2() < sqr(r1 + tolh)) |
---|
799 | { |
---|
800 | if(calcNorm) |
---|
801 | { |
---|
802 | *validNorm = true; |
---|
803 | *n = G4ThreeVector(0, 0, -1) |
---|
804 | + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit(); |
---|
805 | *n = n->unit(); |
---|
806 | } |
---|
807 | return intersection; |
---|
808 | } |
---|
809 | } |
---|
810 | |
---|
811 | if(vRho2 != 0) |
---|
812 | { intersection = (A + std::sqrt(B + sqr(A))) / vRho2; } |
---|
813 | else |
---|
814 | { intersection = ((rho2 - k2) / k1 - p.z()) / v.z(); } |
---|
815 | |
---|
816 | if(calcNorm) |
---|
817 | { |
---|
818 | *validNorm = true; |
---|
819 | *n = G4ThreeVector(p.x() + intersection * v.x(), p.y() |
---|
820 | + intersection * v.y(), - k1 / 2); |
---|
821 | } |
---|
822 | return intersection; |
---|
823 | } |
---|
824 | else |
---|
825 | { |
---|
826 | #ifdef G4SPECSDEBUG |
---|
827 | if(kOutside == Inside(p)) |
---|
828 | { |
---|
829 | G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "Notification", |
---|
830 | JustWarning, "Point p is outside!"); |
---|
831 | } |
---|
832 | else |
---|
833 | G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "Notification", |
---|
834 | JustWarning, "There's an error in this functions code."); |
---|
835 | #endif |
---|
836 | return kInfinity; |
---|
837 | } |
---|
838 | return 0; |
---|
839 | } |
---|
840 | |
---|
841 | /////////////////////////////////////////////////////////////////////////////// |
---|
842 | // |
---|
843 | // Calculate distance (<=actual) to closest surface of shape from inside |
---|
844 | |
---|
845 | G4double G4Paraboloid::DistanceToOut(const G4ThreeVector& p) const |
---|
846 | { |
---|
847 | G4double safe=0.0,rho,safeR,safeZ ; |
---|
848 | G4double tanRMax,secRMax,pRMax ; |
---|
849 | |
---|
850 | #ifdef G4SPECSDEBUG |
---|
851 | if( Inside(p) == kOutside ) |
---|
852 | { |
---|
853 | G4cout.precision(16) ; |
---|
854 | G4cout << G4endl ; |
---|
855 | DumpInfo(); |
---|
856 | G4cout << "Position:" << G4endl << G4endl ; |
---|
857 | G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; |
---|
858 | G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; |
---|
859 | G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; |
---|
860 | G4Exception("G4Paraboloid::DistanceToOut(p)", "Notification", JustWarning, |
---|
861 | "Point p is outside !?" ); |
---|
862 | } |
---|
863 | #endif |
---|
864 | |
---|
865 | rho = p.perp(); |
---|
866 | safeZ = dz - std::fabs(p.z()) ; |
---|
867 | |
---|
868 | tanRMax = (r2 - r1)*0.5/dz ; |
---|
869 | secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; |
---|
870 | pRMax = tanRMax*p.z() + (r1+r2)*0.5 ; |
---|
871 | safeR = (pRMax - rho)/secRMax ; |
---|
872 | |
---|
873 | if (safeZ < safeR) { safe = safeZ; } |
---|
874 | else { safe = safeR; } |
---|
875 | if ( safe < 0.5 * kCarTolerance ) { safe = 0; } |
---|
876 | return safe ; |
---|
877 | } |
---|
878 | |
---|
879 | ////////////////////////////////////////////////////////////////////////// |
---|
880 | // |
---|
881 | // G4EntityType |
---|
882 | |
---|
883 | G4GeometryType G4Paraboloid::GetEntityType() const |
---|
884 | { |
---|
885 | return G4String("G4Paraboloid"); |
---|
886 | } |
---|
887 | |
---|
888 | ////////////////////////////////////////////////////////////////////////// |
---|
889 | // |
---|
890 | // Stream object contents to an output stream |
---|
891 | |
---|
892 | std::ostream& G4Paraboloid::StreamInfo( std::ostream& os ) const |
---|
893 | { |
---|
894 | os << "-----------------------------------------------------------\n" |
---|
895 | << " *** Dump for solid - " << GetName() << " ***\n" |
---|
896 | << " ===================================================\n" |
---|
897 | << " Solid type: G4Paraboloid\n" |
---|
898 | << " Parameters: \n" |
---|
899 | << " z half-axis: " << dz/mm << " mm \n" |
---|
900 | << " radius at -dz: " << r1/mm << " mm \n" |
---|
901 | << " radius at dz: " << r2/mm << " mm \n" |
---|
902 | << "-----------------------------------------------------------\n"; |
---|
903 | |
---|
904 | return os; |
---|
905 | } |
---|
906 | |
---|
907 | //////////////////////////////////////////////////////////////////// |
---|
908 | // |
---|
909 | // GetPointOnSurface |
---|
910 | |
---|
911 | G4ThreeVector G4Paraboloid::GetPointOnSurface() const |
---|
912 | { |
---|
913 | G4double A = (fSurfaceArea == 0)? CalculateSurfaceArea(): fSurfaceArea; |
---|
914 | G4double z = RandFlat::shoot(0.,1.); |
---|
915 | G4double phi = RandFlat::shoot(0., twopi); |
---|
916 | if(pi*(sqr(r1) + sqr(r2))/A >= z) |
---|
917 | { |
---|
918 | G4double rho; |
---|
919 | if(pi * sqr(r1) / A > z) |
---|
920 | { |
---|
921 | rho = RandFlat::shoot(0., 1.); |
---|
922 | rho = std::sqrt(rho); |
---|
923 | rho *= r1; |
---|
924 | return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), -dz); |
---|
925 | } |
---|
926 | else |
---|
927 | { |
---|
928 | rho = RandFlat::shoot(0., 1); |
---|
929 | rho = std::sqrt(rho); |
---|
930 | rho *= r2; |
---|
931 | return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), dz); |
---|
932 | } |
---|
933 | } |
---|
934 | else |
---|
935 | { |
---|
936 | z = RandFlat::shoot(0., 1.)*2*dz - dz; |
---|
937 | return G4ThreeVector(std::sqrt(z*k1 + k2)*std::cos(phi), |
---|
938 | std::sqrt(z*k1 + k2)*std::sin(phi), z); |
---|
939 | } |
---|
940 | } |
---|
941 | |
---|
942 | G4ThreeVectorList* |
---|
943 | G4Paraboloid::CreateRotatedVertices(const G4AffineTransform& pTransform, |
---|
944 | G4int& noPolygonVertices) const |
---|
945 | { |
---|
946 | G4ThreeVectorList *vertices; |
---|
947 | G4ThreeVector vertex; |
---|
948 | G4double meshAnglePhi, cosMeshAnglePhiPer2, |
---|
949 | crossAnglePhi, coscrossAnglePhi, sincrossAnglePhi, sAnglePhi, |
---|
950 | sRho, dRho, rho, lastRho = 0., swapRho; |
---|
951 | G4double rx, ry, rz, k3, k4, zm; |
---|
952 | G4int crossSectionPhi, noPhiCrossSections, noRhoSections; |
---|
953 | |
---|
954 | // Phi cross sections |
---|
955 | // |
---|
956 | noPhiCrossSections = G4int(twopi/kMeshAngleDefault)+1; |
---|
957 | |
---|
958 | if (noPhiCrossSections<kMinMeshSections) |
---|
959 | { |
---|
960 | noPhiCrossSections=kMinMeshSections; |
---|
961 | } |
---|
962 | else if (noPhiCrossSections>kMaxMeshSections) |
---|
963 | { |
---|
964 | noPhiCrossSections=kMaxMeshSections; |
---|
965 | } |
---|
966 | meshAnglePhi=twopi/(noPhiCrossSections-1); |
---|
967 | |
---|
968 | sAnglePhi = -meshAnglePhi*0.5*0; |
---|
969 | cosMeshAnglePhiPer2 = std::cos(meshAnglePhi / 2.); |
---|
970 | |
---|
971 | noRhoSections = G4int(pi/2/kMeshAngleDefault) + 1; |
---|
972 | |
---|
973 | // There is no obvious value for noRhoSections, at the moment the parabola is |
---|
974 | // viewed as a quarter circle mean this formula for it. |
---|
975 | |
---|
976 | // An alternetive would be to calculate max deviation from parabola and |
---|
977 | // keep adding new vertices there until it was under a decided constant. |
---|
978 | |
---|
979 | // maxDeviation on a line between points (rho1, z1) and (rho2, z2) is given |
---|
980 | // by rhoMax = sqrt(k1 * z + k2) - z * (rho2 - rho1) |
---|
981 | // / (z2 - z1) - (rho1 * z2 - rho2 * z1) / (z2 - z1) |
---|
982 | // where z is k1 / 2 * (rho1 + rho2) - k2 / k1 |
---|
983 | |
---|
984 | sRho = r1; |
---|
985 | dRho = (r2 - r1) / double(noRhoSections - 1); |
---|
986 | |
---|
987 | vertices=new G4ThreeVectorList(); |
---|
988 | |
---|
989 | if (vertices) |
---|
990 | { |
---|
991 | for (crossSectionPhi=0; crossSectionPhi<noPhiCrossSections; |
---|
992 | crossSectionPhi++) |
---|
993 | { |
---|
994 | crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi; |
---|
995 | coscrossAnglePhi=std::cos(crossAnglePhi); |
---|
996 | sincrossAnglePhi=std::sin(crossAnglePhi); |
---|
997 | lastRho = 0; |
---|
998 | for (int iRho=0; iRho < noRhoSections; |
---|
999 | iRho++) |
---|
1000 | { |
---|
1001 | // Compute coordinates of cross section at section crossSectionPhi |
---|
1002 | // |
---|
1003 | if(iRho == noRhoSections - 1) |
---|
1004 | { |
---|
1005 | rho = r2; |
---|
1006 | } |
---|
1007 | else |
---|
1008 | { |
---|
1009 | rho = iRho * dRho + sRho; |
---|
1010 | |
---|
1011 | // This part is to ensure that the vertices |
---|
1012 | // will form a volume larger than the paraboloid |
---|
1013 | |
---|
1014 | k3 = k1 / (2*rho + dRho); |
---|
1015 | k4 = rho - k3 * (sqr(rho) - k2) / k1; |
---|
1016 | zm = (sqr(k1 / (2 * k3)) - k2) / k1; |
---|
1017 | rho += std::sqrt(k1 * zm + k2) - zm * k3 - k4; |
---|
1018 | } |
---|
1019 | |
---|
1020 | rho += (1 / cosMeshAnglePhiPer2 - 1) * (iRho * dRho + sRho); |
---|
1021 | |
---|
1022 | if(rho < lastRho) |
---|
1023 | { |
---|
1024 | swapRho = lastRho; |
---|
1025 | lastRho = rho + dRho; |
---|
1026 | rho = swapRho; |
---|
1027 | } |
---|
1028 | else |
---|
1029 | { |
---|
1030 | lastRho = rho + dRho; |
---|
1031 | } |
---|
1032 | |
---|
1033 | rx = coscrossAnglePhi*rho; |
---|
1034 | ry = sincrossAnglePhi*rho; |
---|
1035 | rz = (sqr(iRho * dRho + sRho) - k2) / k1; |
---|
1036 | vertex = G4ThreeVector(rx,ry,rz); |
---|
1037 | vertices->push_back(pTransform.TransformPoint(vertex)); |
---|
1038 | } |
---|
1039 | } // Phi |
---|
1040 | noPolygonVertices = noRhoSections ; |
---|
1041 | } |
---|
1042 | else |
---|
1043 | { |
---|
1044 | DumpInfo(); |
---|
1045 | G4Exception("G4Paraboloid::CreateRotatedVertices()", |
---|
1046 | "FatalError", FatalException, |
---|
1047 | "Error in allocation of vertices. Out of memory !"); |
---|
1048 | } |
---|
1049 | return vertices; |
---|
1050 | } |
---|
1051 | |
---|
1052 | ///////////////////////////////////////////////////////////////////////////// |
---|
1053 | // |
---|
1054 | // Methods for visualisation |
---|
1055 | |
---|
1056 | void G4Paraboloid::DescribeYourselfTo (G4VGraphicsScene& scene) const |
---|
1057 | { |
---|
1058 | scene.AddSolid(*this); |
---|
1059 | } |
---|
1060 | |
---|
1061 | G4NURBS* G4Paraboloid::CreateNURBS () const |
---|
1062 | { |
---|
1063 | // Box for now!!! |
---|
1064 | // |
---|
1065 | return new G4NURBSbox(r1, r1, dz); |
---|
1066 | } |
---|
1067 | |
---|
1068 | G4Polyhedron* G4Paraboloid::CreatePolyhedron () const |
---|
1069 | { |
---|
1070 | return new G4PolyhedronParaboloid(r1, r2, dz, 0., twopi); |
---|
1071 | } |
---|
1072 | |
---|
1073 | |
---|
1074 | G4Polyhedron* G4Paraboloid::GetPolyhedron () const |
---|
1075 | { |
---|
1076 | if (!fpPolyhedron || |
---|
1077 | fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() != |
---|
1078 | fpPolyhedron->GetNumberOfRotationSteps()) |
---|
1079 | { |
---|
1080 | delete fpPolyhedron; |
---|
1081 | fpPolyhedron = CreatePolyhedron(); |
---|
1082 | } |
---|
1083 | return fpPolyhedron; |
---|
1084 | } |
---|