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