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