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