Geant4  10.00.p01
G4Tubs.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 //
27 // $Id: G4Tubs.cc 76263 2013-11-08 11:41:52Z gcosmo $
28 //
29 //
30 // class G4Tubs
31 //
32 // History:
33 //
34 // 05.04.12 M.Kelsey: Use sqrt(r) in GetPointOnSurface() for uniform points
35 // 02.08.07 T.Nikitina: bug fixed in DistanceToOut(p,v,..) for negative value under sqrt
36 // for the case: p on the surface and v is tangent to the surface
37 // 11.05.07 T.Nikitina: bug fixed in DistanceToOut(p,v,..) for phi < 2pi
38 // 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
39 // 16.03.05 V.Grichine: SurfaceNormal(p) with edges/corners for boolean
40 // 20.07.01 V.Grichine: bug fixed in Inside(p)
41 // 20.02.01 V.Grichine: bug fixed in Inside(p) and CalculateExtent was
42 // simplified base on G4Box::CalculateExtent
43 // 07.12.00 V.Grichine: phi-section algorithm was changed in Inside(p)
44 // 28.11.00 V.Grichine: bug fixed in Inside(p)
45 // 31.10.00 V.Grichine: assign srd, sphi in Distance ToOut(p,v,...)
46 // 08.08.00 V.Grichine: more stable roots of 2-equation in DistanceToOut(p,v,..)
47 // 02.08.00 V.Grichine: point is outside check in Distance ToOut(p)
48 // 17.05.00 V.Grichine: bugs (#76,#91) fixed in Distance ToOut(p,v,...)
49 // 31.03.00 V.Grichine: bug fixed in Inside(p)
50 // 19.11.99 V.Grichine: side = kNull in DistanceToOut(p,v,...)
51 // 13.10.99 V.Grichine: bugs fixed in DistanceToIn(p,v)
52 // 28.05.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...)
53 // 25.05.99 V.Grichine: bugs fixed in DistanceToIn(p,v)
54 // 23.03.99 V.Grichine: bug fixed in DistanceToIn(p,v)
55 // 09.10.98 V.Grichine: modifications in DistanceToOut(p,v,...)
56 // 18.06.98 V.Grichine: n-normalisation in DistanceToOut(p,v)
57 //
58 // 1994-95 P.Kent: implementation
59 //
61 
62 #include "G4Tubs.hh"
63 
64 #if !defined(G4GEOM_USE_UTUBS)
65 
66 #include "G4VoxelLimits.hh"
67 #include "G4AffineTransform.hh"
68 #include "G4GeometryTolerance.hh"
69 
70 #include "G4VPVParameterisation.hh"
71 
72 #include "Randomize.hh"
73 
74 #include "meshdefs.hh"
75 
76 #include "G4VGraphicsScene.hh"
77 #include "G4Polyhedron.hh"
78 
79 using namespace CLHEP;
80 
82 //
83 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
84 // - note if pdphi>2PI then reset to 2PI
85 
86 G4Tubs::G4Tubs( const G4String &pName,
87  G4double pRMin, G4double pRMax,
88  G4double pDz,
89  G4double pSPhi, G4double pDPhi )
90  : G4CSGSolid(pName), fRMin(pRMin), fRMax(pRMax), fDz(pDz), fSPhi(0), fDPhi(0)
91 {
92 
95 
99 
100  if (pDz<=0) // Check z-len
101  {
102  std::ostringstream message;
103  message << "Negative Z half-length (" << pDz << ") in solid: " << GetName();
104  G4Exception("G4Tubs::G4Tubs()", "GeomSolids0002", FatalException, message);
105  }
106  if ( (pRMin >= pRMax) || (pRMin < 0) ) // Check radii
107  {
108  std::ostringstream message;
109  message << "Invalid values for radii in solid: " << GetName()
110  << G4endl
111  << " pRMin = " << pRMin << ", pRMax = " << pRMax;
112  G4Exception("G4Tubs::G4Tubs()", "GeomSolids0002", FatalException, message);
113  }
114 
115  // Check angles
116  //
117  CheckPhiAngles(pSPhi, pDPhi);
118 }
119 
121 //
122 // Fake default constructor - sets only member data and allocates memory
123 // for usage restricted to object persistency.
124 //
125 G4Tubs::G4Tubs( __void__& a )
126  : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
127  fRMin(0.), fRMax(0.), fDz(0.), fSPhi(0.), fDPhi(0.),
128  sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
129  sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
130  fPhiFullTube(false), halfCarTolerance(0.), halfRadTolerance(0.),
131  halfAngTolerance(0.)
132 
133 {
134 }
135 
137 //
138 // Destructor
139 
141 {
142 }
143 
145 //
146 // Copy constructor
147 
149  : G4CSGSolid(rhs),
150  kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
151  fRMin(rhs.fRMin), fRMax(rhs.fRMax), fDz(rhs.fDz),
152  fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
153  sinCPhi(rhs.sinCPhi), cosCPhi(rhs.sinCPhi),
154  cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiOT),
155  sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
156  sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi), fPhiFullTube(rhs.fPhiFullTube),
157  halfCarTolerance(rhs.halfCarTolerance),
158  halfRadTolerance(rhs.halfRadTolerance),
159  halfAngTolerance(rhs.halfAngTolerance)
160 {
161 }
162 
164 //
165 // Assignment operator
166 
168 {
169  // Check assignment to self
170  //
171  if (this == &rhs) { return *this; }
172 
173  // Copy base class data
174  //
176 
177  // Copy data
178  //
180  fRMin = rhs.fRMin; fRMax = rhs.fRMax; fDz = rhs.fDz;
181  fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
182  sinCPhi = rhs.sinCPhi; cosCPhi = rhs.sinCPhi;
184  sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
185  sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
190 
191  return *this;
192 }
193 
195 //
196 // Dispatch to parameterisation for replication mechanism dimension
197 // computation & modification.
198 
200  const G4int n,
201  const G4VPhysicalVolume* pRep )
202 {
203  p->ComputeDimensions(*this,n,pRep) ;
204 }
205 
207 //
208 // Calculate extent under transform and specified limit
209 
211  const G4VoxelLimits& pVoxelLimit,
212  const G4AffineTransform& pTransform,
213  G4double& pMin,
214  G4double& pMax ) const
215 {
216 
217  if ( (!pTransform.IsRotated()) && (fDPhi == twopi) && (fRMin == 0) )
218  {
219  // Special case handling for unrotated solid tubes
220  // Compute x/y/z mins and maxs fro bounding box respecting limits,
221  // with early returns if outside limits. Then switch() on pAxis,
222  // and compute exact x and y limit for x/y case
223 
224  G4double xoffset, xMin, xMax;
225  G4double yoffset, yMin, yMax;
226  G4double zoffset, zMin, zMax;
227 
228  G4double diff1, diff2, maxDiff, newMin, newMax;
229  G4double xoff1, xoff2, yoff1, yoff2, delta;
230 
231  xoffset = pTransform.NetTranslation().x();
232  xMin = xoffset - fRMax;
233  xMax = xoffset + fRMax;
234 
235  if (pVoxelLimit.IsXLimited())
236  {
237  if ( (xMin > pVoxelLimit.GetMaxXExtent())
238  || (xMax < pVoxelLimit.GetMinXExtent()) )
239  {
240  return false;
241  }
242  else
243  {
244  if (xMin < pVoxelLimit.GetMinXExtent())
245  {
246  xMin = pVoxelLimit.GetMinXExtent();
247  }
248  if (xMax > pVoxelLimit.GetMaxXExtent())
249  {
250  xMax = pVoxelLimit.GetMaxXExtent();
251  }
252  }
253  }
254  yoffset = pTransform.NetTranslation().y();
255  yMin = yoffset - fRMax;
256  yMax = yoffset + fRMax;
257 
258  if ( pVoxelLimit.IsYLimited() )
259  {
260  if ( (yMin > pVoxelLimit.GetMaxYExtent())
261  || (yMax < pVoxelLimit.GetMinYExtent()) )
262  {
263  return false;
264  }
265  else
266  {
267  if (yMin < pVoxelLimit.GetMinYExtent())
268  {
269  yMin = pVoxelLimit.GetMinYExtent();
270  }
271  if (yMax > pVoxelLimit.GetMaxYExtent())
272  {
273  yMax=pVoxelLimit.GetMaxYExtent();
274  }
275  }
276  }
277  zoffset = pTransform.NetTranslation().z();
278  zMin = zoffset - fDz;
279  zMax = zoffset + fDz;
280 
281  if ( pVoxelLimit.IsZLimited() )
282  {
283  if ( (zMin > pVoxelLimit.GetMaxZExtent())
284  || (zMax < pVoxelLimit.GetMinZExtent()) )
285  {
286  return false;
287  }
288  else
289  {
290  if (zMin < pVoxelLimit.GetMinZExtent())
291  {
292  zMin = pVoxelLimit.GetMinZExtent();
293  }
294  if (zMax > pVoxelLimit.GetMaxZExtent())
295  {
296  zMax = pVoxelLimit.GetMaxZExtent();
297  }
298  }
299  }
300  switch ( pAxis ) // Known to cut cylinder
301  {
302  case kXAxis :
303  {
304  yoff1 = yoffset - yMin;
305  yoff2 = yMax - yoffset;
306 
307  if ( (yoff1 >= 0) && (yoff2 >= 0) ) // Y limits cross max/min x
308  { // => no change
309  pMin = xMin;
310  pMax = xMax;
311  }
312  else
313  {
314  // Y limits don't cross max/min x => compute max delta x,
315  // hence new mins/maxs
316 
317  delta = fRMax*fRMax - yoff1*yoff1;
318  diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
319  delta = fRMax*fRMax - yoff2*yoff2;
320  diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
321  maxDiff = (diff1 > diff2) ? diff1:diff2;
322  newMin = xoffset - maxDiff;
323  newMax = xoffset + maxDiff;
324  pMin = (newMin < xMin) ? xMin : newMin;
325  pMax = (newMax > xMax) ? xMax : newMax;
326  }
327  break;
328  }
329  case kYAxis :
330  {
331  xoff1 = xoffset - xMin;
332  xoff2 = xMax - xoffset;
333 
334  if ( (xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
335  { // => no change
336  pMin = yMin;
337  pMax = yMax;
338  }
339  else
340  {
341  // X limits don't cross max/min y => compute max delta y,
342  // hence new mins/maxs
343 
344  delta = fRMax*fRMax - xoff1*xoff1;
345  diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
346  delta = fRMax*fRMax - xoff2*xoff2;
347  diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
348  maxDiff = (diff1 > diff2) ? diff1 : diff2;
349  newMin = yoffset - maxDiff;
350  newMax = yoffset + maxDiff;
351  pMin = (newMin < yMin) ? yMin : newMin;
352  pMax = (newMax > yMax) ? yMax : newMax;
353  }
354  break;
355  }
356  case kZAxis:
357  {
358  pMin = zMin;
359  pMax = zMax;
360  break;
361  }
362  default:
363  break;
364  }
365  pMin -= kCarTolerance;
366  pMax += kCarTolerance;
367  return true;
368  }
369  else // Calculate rotated vertex coordinates
370  {
371  G4int i, noEntries, noBetweenSections4;
372  G4bool existsAfterClip = false;
373  G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform);
374 
375  pMin = kInfinity;
376  pMax = -kInfinity;
377 
378  noEntries = vertices->size();
379  noBetweenSections4 = noEntries - 4;
380 
381  for ( i = 0 ; i < noEntries ; i += 4 )
382  {
383  ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
384  }
385  for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
386  {
387  ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax);
388  }
389  if ( (pMin != kInfinity) || (pMax != -kInfinity) )
390  {
391  existsAfterClip = true;
392  pMin -= kCarTolerance; // Add 2*tolerance to avoid precision troubles
393  pMax += kCarTolerance;
394  }
395  else
396  {
397  // Check for case where completely enveloping clipping volume
398  // If point inside then we are confident that the solid completely
399  // envelopes the clipping volume. Hence set min/max extents according
400  // to clipping volume extents along the specified axis.
401 
402  G4ThreeVector clipCentre(
403  (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
404  (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
405  (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 );
406 
407  if ( Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
408  {
409  existsAfterClip = true;
410  pMin = pVoxelLimit.GetMinExtent(pAxis);
411  pMax = pVoxelLimit.GetMaxExtent(pAxis);
412  }
413  }
414  delete vertices;
415  return existsAfterClip;
416  }
417 }
418 
419 
421 //
422 // Return whether point inside/outside/on surface
423 
425 {
426  G4double r2,pPhi,tolRMin,tolRMax;
427  EInside in = kOutside ;
428 
429  if (std::fabs(p.z()) <= fDz - halfCarTolerance)
430  {
431  r2 = p.x()*p.x() + p.y()*p.y() ;
432 
433  if (fRMin) { tolRMin = fRMin + halfRadTolerance ; }
434  else { tolRMin = 0 ; }
435 
436  tolRMax = fRMax - halfRadTolerance ;
437 
438  if ((r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax))
439  {
440  if ( fPhiFullTube )
441  {
442  in = kInside ;
443  }
444  else
445  {
446  // Try inner tolerant phi boundaries (=>inside)
447  // if not inside, try outer tolerant phi boundaries
448 
449  if ( (tolRMin==0) && (std::fabs(p.x())<=halfCarTolerance)
450  && (std::fabs(p.y())<=halfCarTolerance) )
451  {
452  in=kSurface;
453  }
454  else
455  {
456  pPhi = std::atan2(p.y(),p.x()) ;
457  if ( pPhi < -halfAngTolerance ) { pPhi += twopi; } // 0<=pPhi<2pi
458 
459  if ( fSPhi >= 0 )
460  {
461  if ( (std::fabs(pPhi) < halfAngTolerance)
462  && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
463  {
464  pPhi += twopi ; // 0 <= pPhi < 2pi
465  }
466  if ( (pPhi >= fSPhi + halfAngTolerance)
467  && (pPhi <= fSPhi + fDPhi - halfAngTolerance) )
468  {
469  in = kInside ;
470  }
471  else if ( (pPhi >= fSPhi - halfAngTolerance)
472  && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
473  {
474  in = kSurface ;
475  }
476  }
477  else // fSPhi < 0
478  {
479  if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
480  && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;} //kOutside
481  else if ( (pPhi <= fSPhi + twopi + halfAngTolerance)
482  && (pPhi >= fSPhi + fDPhi - halfAngTolerance) )
483  {
484  in = kSurface ;
485  }
486  else
487  {
488  in = kInside ;
489  }
490  }
491  }
492  }
493  }
494  else // Try generous boundaries
495  {
496  tolRMin = fRMin - halfRadTolerance ;
497  tolRMax = fRMax + halfRadTolerance ;
498 
499  if ( tolRMin < 0 ) { tolRMin = 0; }
500 
501  if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
502  {
503  if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance) )
504  { // Continuous in phi or on z-axis
505  in = kSurface ;
506  }
507  else // Try outer tolerant phi boundaries only
508  {
509  pPhi = std::atan2(p.y(),p.x()) ;
510 
511  if ( pPhi < -halfAngTolerance) { pPhi += twopi; } // 0<=pPhi<2pi
512  if ( fSPhi >= 0 )
513  {
514  if ( (std::fabs(pPhi) < halfAngTolerance)
515  && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
516  {
517  pPhi += twopi ; // 0 <= pPhi < 2pi
518  }
519  if ( (pPhi >= fSPhi - halfAngTolerance)
520  && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
521  {
522  in = kSurface ;
523  }
524  }
525  else // fSPhi < 0
526  {
527  if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
528  && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;} // kOutside
529  else
530  {
531  in = kSurface ;
532  }
533  }
534  }
535  }
536  }
537  }
538  else if (std::fabs(p.z()) <= fDz + halfCarTolerance)
539  { // Check within tolerant r limits
540  r2 = p.x()*p.x() + p.y()*p.y() ;
541  tolRMin = fRMin - halfRadTolerance ;
542  tolRMax = fRMax + halfRadTolerance ;
543 
544  if ( tolRMin < 0 ) { tolRMin = 0; }
545 
546  if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
547  {
548  if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance))
549  { // Continuous in phi or on z-axis
550  in = kSurface ;
551  }
552  else // Try outer tolerant phi boundaries
553  {
554  pPhi = std::atan2(p.y(),p.x()) ;
555 
556  if ( pPhi < -halfAngTolerance ) { pPhi += twopi; } // 0<=pPhi<2pi
557  if ( fSPhi >= 0 )
558  {
559  if ( (std::fabs(pPhi) < halfAngTolerance)
560  && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
561  {
562  pPhi += twopi ; // 0 <= pPhi < 2pi
563  }
564  if ( (pPhi >= fSPhi - halfAngTolerance)
565  && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
566  {
567  in = kSurface;
568  }
569  }
570  else // fSPhi < 0
571  {
572  if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
573  && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;}
574  else
575  {
576  in = kSurface ;
577  }
578  }
579  }
580  }
581  }
582  return in;
583 }
584 
586 //
587 // Return unit normal of surface closest to p
588 // - note if point on z axis, ignore phi divided sides
589 // - unsafe if point close to z axis a rmin=0 - no explicit checks
590 
592 {
593  G4int noSurfaces = 0;
594  G4double rho, pPhi;
595  G4double distZ, distRMin, distRMax;
596  G4double distSPhi = kInfinity, distEPhi = kInfinity;
597 
598  G4ThreeVector norm, sumnorm(0.,0.,0.);
599  G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
600  G4ThreeVector nR, nPs, nPe;
601 
602  rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
603 
604  distRMin = std::fabs(rho - fRMin);
605  distRMax = std::fabs(rho - fRMax);
606  distZ = std::fabs(std::fabs(p.z()) - fDz);
607 
608  if (!fPhiFullTube) // Protected against (0,0,z)
609  {
610  if ( rho > halfCarTolerance )
611  {
612  pPhi = std::atan2(p.y(),p.x());
613 
614  if(pPhi < fSPhi- halfCarTolerance) { pPhi += twopi; }
615  else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; }
616 
617  distSPhi = std::fabs(pPhi - fSPhi);
618  distEPhi = std::fabs(pPhi - fSPhi - fDPhi);
619  }
620  else if( !fRMin )
621  {
622  distSPhi = 0.;
623  distEPhi = 0.;
624  }
625  nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
626  nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
627  }
628  if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
629 
630  if( distRMax <= halfCarTolerance )
631  {
632  noSurfaces ++;
633  sumnorm += nR;
634  }
635  if( fRMin && (distRMin <= halfCarTolerance) )
636  {
637  noSurfaces ++;
638  sumnorm -= nR;
639  }
640  if( fDPhi < twopi )
641  {
642  if (distSPhi <= halfAngTolerance)
643  {
644  noSurfaces ++;
645  sumnorm += nPs;
646  }
647  if (distEPhi <= halfAngTolerance)
648  {
649  noSurfaces ++;
650  sumnorm += nPe;
651  }
652  }
653  if (distZ <= halfCarTolerance)
654  {
655  noSurfaces ++;
656  if ( p.z() >= 0.) { sumnorm += nZ; }
657  else { sumnorm -= nZ; }
658  }
659  if ( noSurfaces == 0 )
660  {
661 #ifdef G4CSGDEBUG
662  G4Exception("G4Tubs::SurfaceNormal(p)", "GeomSolids1002",
663  JustWarning, "Point p is not on surface !?" );
664  G4int oldprc = G4cout.precision(20);
665  G4cout<< "G4Tubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); "
666  << G4endl << G4endl;
667  G4cout.precision(oldprc) ;
668 #endif
669  norm = ApproxSurfaceNormal(p);
670  }
671  else if ( noSurfaces == 1 ) { norm = sumnorm; }
672  else { norm = sumnorm.unit(); }
673 
674  return norm;
675 }
676 
678 //
679 // Algorithm for SurfaceNormal() following the original specification
680 // for points not on the surface
681 
683 {
684  ENorm side ;
685  G4ThreeVector norm ;
686  G4double rho, phi ;
687  G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
688 
689  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
690 
691  distRMin = std::fabs(rho - fRMin) ;
692  distRMax = std::fabs(rho - fRMax) ;
693  distZ = std::fabs(std::fabs(p.z()) - fDz) ;
694 
695  if (distRMin < distRMax) // First minimum
696  {
697  if ( distZ < distRMin )
698  {
699  distMin = distZ ;
700  side = kNZ ;
701  }
702  else
703  {
704  distMin = distRMin ;
705  side = kNRMin ;
706  }
707  }
708  else
709  {
710  if ( distZ < distRMax )
711  {
712  distMin = distZ ;
713  side = kNZ ;
714  }
715  else
716  {
717  distMin = distRMax ;
718  side = kNRMax ;
719  }
720  }
721  if (!fPhiFullTube && rho ) // Protected against (0,0,z)
722  {
723  phi = std::atan2(p.y(),p.x()) ;
724 
725  if ( phi < 0 ) { phi += twopi; }
726 
727  if ( fSPhi < 0 )
728  {
729  distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;
730  }
731  else
732  {
733  distSPhi = std::fabs(phi - fSPhi)*rho ;
734  }
735  distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
736 
737  if (distSPhi < distEPhi) // Find new minimum
738  {
739  if ( distSPhi < distMin )
740  {
741  side = kNSPhi ;
742  }
743  }
744  else
745  {
746  if ( distEPhi < distMin )
747  {
748  side = kNEPhi ;
749  }
750  }
751  }
752  switch ( side )
753  {
754  case kNRMin : // Inner radius
755  {
756  norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ;
757  break ;
758  }
759  case kNRMax : // Outer radius
760  {
761  norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ;
762  break ;
763  }
764  case kNZ : // + or - dz
765  {
766  if ( p.z() > 0 ) { norm = G4ThreeVector(0,0,1) ; }
767  else { norm = G4ThreeVector(0,0,-1); }
768  break ;
769  }
770  case kNSPhi:
771  {
772  norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
773  break ;
774  }
775  case kNEPhi:
776  {
777  norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
778  break;
779  }
780  default: // Should never reach this case ...
781  {
782  DumpInfo();
783  G4Exception("G4Tubs::ApproxSurfaceNormal()",
784  "GeomSolids1002", JustWarning,
785  "Undefined side for valid surface normal to solid.");
786  break ;
787  }
788  }
789  return norm;
790 }
791 
793 //
794 //
795 // Calculate distance to shape from outside, along normalised vector
796 // - return kInfinity if no intersection, or intersection distance <= tolerance
797 //
798 // - Compute the intersection with the z planes
799 // - if at valid r, phi, return
800 //
801 // -> If point is outer outer radius, compute intersection with rmax
802 // - if at valid phi,z return
803 //
804 // -> Compute intersection with inner radius, taking largest +ve root
805 // - if valid (in z,phi), save intersction
806 //
807 // -> If phi segmented, compute intersections with phi half planes
808 // - return smallest of valid phi intersections and
809 // inner radius intersection
810 //
811 // NOTE:
812 // - 'if valid' implies tolerant checking of intersection points
813 
815  const G4ThreeVector& v ) const
816 {
817  G4double snxt = kInfinity ; // snxt = default return value
818  G4double tolORMin2, tolIRMax2 ; // 'generous' radii squared
819  G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
820  const G4double dRmax = 100.*fRMax;
821 
822  // Intersection point variables
823  //
824  G4double Dist, sd, xi, yi, zi, rho2, inum, iden, cosPsi, Comp ;
825  G4double t1, t2, t3, b, c, d ; // Quadratic solver variables
826 
827  // Calculate tolerant rmin and rmax
828 
829  if (fRMin > kRadTolerance)
830  {
831  tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ;
832  tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ;
833  }
834  else
835  {
836  tolORMin2 = 0.0 ;
837  tolIRMin2 = 0.0 ;
838  }
839  tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ;
840  tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - halfRadTolerance) ;
841 
842  // Intersection with Z surfaces
843 
844  tolIDz = fDz - halfCarTolerance ;
845  tolODz = fDz + halfCarTolerance ;
846 
847  if (std::fabs(p.z()) >= tolIDz)
848  {
849  if ( p.z()*v.z() < 0 ) // at +Z going in -Z or visa versa
850  {
851  sd = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ; // Z intersect distance
852 
853  if(sd < 0.0) { sd = 0.0; }
854 
855  xi = p.x() + sd*v.x() ; // Intersection coords
856  yi = p.y() + sd*v.y() ;
857  rho2 = xi*xi + yi*yi ;
858 
859  // Check validity of intersection
860 
861  if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
862  {
863  if (!fPhiFullTube && rho2)
864  {
865  // Psi = angle made with central (average) phi of shape
866  //
867  inum = xi*cosCPhi + yi*sinCPhi ;
868  iden = std::sqrt(rho2) ;
869  cosPsi = inum/iden ;
870  if (cosPsi >= cosHDPhiIT) { return sd ; }
871  }
872  else
873  {
874  return sd ;
875  }
876  }
877  }
878  else
879  {
880  if ( snxt<halfCarTolerance ) { snxt=0; }
881  return snxt ; // On/outside extent, and heading away
882  // -> cannot intersect
883  }
884  }
885 
886  // -> Can not intersect z surfaces
887  //
888  // Intersection with rmax (possible return) and rmin (must also check phi)
889  //
890  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
891  //
892  // Intersects with x^2+y^2=R^2
893  //
894  // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
895  // t1 t2 t3
896 
897  t1 = 1.0 - v.z()*v.z() ;
898  t2 = p.x()*v.x() + p.y()*v.y() ;
899  t3 = p.x()*p.x() + p.y()*p.y() ;
900 
901  if ( t1 > 0 ) // Check not || to z axis
902  {
903  b = t2/t1 ;
904  c = t3 - fRMax*fRMax ;
905  if ((t3 >= tolORMax2) && (t2<0)) // This also handles the tangent case
906  {
907  // Try outer cylinder intersection
908  // c=(t3-fRMax*fRMax)/t1;
909 
910  c /= t1 ;
911  d = b*b - c ;
912 
913  if (d >= 0) // If real root
914  {
915  sd = c/(-b+std::sqrt(d));
916  if (sd >= 0) // If 'forwards'
917  {
918  if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
919  { // 64 bits systems. Split long distances and recompute
920  G4double fTerm = sd-std::fmod(sd,dRmax);
921  sd = fTerm + DistanceToIn(p+fTerm*v,v);
922  }
923  // Check z intersection
924  //
925  zi = p.z() + sd*v.z() ;
926  if (std::fabs(zi)<=tolODz)
927  {
928  // Z ok. Check phi intersection if reqd
929  //
930  if (fPhiFullTube)
931  {
932  return sd ;
933  }
934  else
935  {
936  xi = p.x() + sd*v.x() ;
937  yi = p.y() + sd*v.y() ;
938  cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ;
939  if (cosPsi >= cosHDPhiIT) { return sd ; }
940  }
941  } // end if std::fabs(zi)
942  } // end if (sd>=0)
943  } // end if (d>=0)
944  } // end if (r>=fRMax)
945  else
946  {
947  // Inside outer radius :
948  // check not inside, and heading through tubs (-> 0 to in)
949 
950  if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.z()) <= tolIDz))
951  {
952  // Inside both radii, delta r -ve, inside z extent
953 
954  if (!fPhiFullTube)
955  {
956  inum = p.x()*cosCPhi + p.y()*sinCPhi ;
957  iden = std::sqrt(t3) ;
958  cosPsi = inum/iden ;
959  if (cosPsi >= cosHDPhiIT)
960  {
961  // In the old version, the small negative tangent for the point
962  // on surface was not taken in account, and returning 0.0 ...
963  // New version: check the tangent for the point on surface and
964  // if no intersection, return kInfinity, if intersection instead
965  // return sd.
966  //
967  c = t3-fRMax*fRMax;
968  if ( c<=0.0 )
969  {
970  return 0.0;
971  }
972  else
973  {
974  c = c/t1 ;
975  d = b*b-c;
976  if ( d>=0.0 )
977  {
978  snxt = c/(-b+std::sqrt(d)); // using safe solution
979  // for quadratic equation
980  if ( snxt < halfCarTolerance ) { snxt=0; }
981  return snxt ;
982  }
983  else
984  {
985  return kInfinity;
986  }
987  }
988  }
989  }
990  else
991  {
992  // In the old version, the small negative tangent for the point
993  // on surface was not taken in account, and returning 0.0 ...
994  // New version: check the tangent for the point on surface and
995  // if no intersection, return kInfinity, if intersection instead
996  // return sd.
997  //
998  c = t3 - fRMax*fRMax;
999  if ( c<=0.0 )
1000  {
1001  return 0.0;
1002  }
1003  else
1004  {
1005  c = c/t1 ;
1006  d = b*b-c;
1007  if ( d>=0.0 )
1008  {
1009  snxt= c/(-b+std::sqrt(d)); // using safe solution
1010  // for quadratic equation
1011  if ( snxt < halfCarTolerance ) { snxt=0; }
1012  return snxt ;
1013  }
1014  else
1015  {
1016  return kInfinity;
1017  }
1018  }
1019  } // end if (!fPhiFullTube)
1020  } // end if (t3>tolIRMin2)
1021  } // end if (Inside Outer Radius)
1022  if ( fRMin ) // Try inner cylinder intersection
1023  {
1024  c = (t3 - fRMin*fRMin)/t1 ;
1025  d = b*b - c ;
1026  if ( d >= 0.0 ) // If real root
1027  {
1028  // Always want 2nd root - we are outside and know rmax Hit was bad
1029  // - If on surface of rmin also need farthest root
1030 
1031  sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
1032  if (sd >= -halfCarTolerance) // check forwards
1033  {
1034  // Check z intersection
1035  //
1036  if(sd < 0.0) { sd = 0.0; }
1037  if ( sd>dRmax ) // Avoid rounding errors due to precision issues seen
1038  { // 64 bits systems. Split long distances and recompute
1039  G4double fTerm = sd-std::fmod(sd,dRmax);
1040  sd = fTerm + DistanceToIn(p+fTerm*v,v);
1041  }
1042  zi = p.z() + sd*v.z() ;
1043  if (std::fabs(zi) <= tolODz)
1044  {
1045  // Z ok. Check phi
1046  //
1047  if ( fPhiFullTube )
1048  {
1049  return sd ;
1050  }
1051  else
1052  {
1053  xi = p.x() + sd*v.x() ;
1054  yi = p.y() + sd*v.y() ;
1055  cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMin ;
1056  if (cosPsi >= cosHDPhiIT)
1057  {
1058  // Good inner radius isect
1059  // - but earlier phi isect still possible
1060 
1061  snxt = sd ;
1062  }
1063  }
1064  } // end if std::fabs(zi)
1065  } // end if (sd>=0)
1066  } // end if (d>=0)
1067  } // end if (fRMin)
1068  }
1069 
1070  // Phi segment intersection
1071  //
1072  // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1073  //
1074  // o NOTE: Large duplication of code between sphi & ephi checks
1075  // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1076  // intersection check <=0 -> >=0
1077  // -> use some form of loop Construct ?
1078  //
1079  if ( !fPhiFullTube )
1080  {
1081  // First phi surface (Starting phi)
1082  //
1083  Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1084 
1085  if ( Comp < 0 ) // Component in outwards normal dirn
1086  {
1087  Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1088 
1089  if ( Dist < halfCarTolerance )
1090  {
1091  sd = Dist/Comp ;
1092 
1093  if (sd < snxt)
1094  {
1095  if ( sd < 0 ) { sd = 0.0; }
1096  zi = p.z() + sd*v.z() ;
1097  if ( std::fabs(zi) <= tolODz )
1098  {
1099  xi = p.x() + sd*v.x() ;
1100  yi = p.y() + sd*v.y() ;
1101  rho2 = xi*xi + yi*yi ;
1102 
1103  if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1104  || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1105  && ( v.y()*cosSPhi - v.x()*sinSPhi > 0 )
1106  && ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 ) )
1107  || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1108  && (v.y()*cosSPhi - v.x()*sinSPhi > 0)
1109  && (v.x()*cosSPhi + v.y()*sinSPhi < 0) ) )
1110  {
1111  // z and r intersections good
1112  // - check intersecting with correct half-plane
1113  //
1114  if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = sd; }
1115  }
1116  }
1117  }
1118  }
1119  }
1120 
1121  // Second phi surface (Ending phi)
1122 
1123  Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1124 
1125  if (Comp < 0 ) // Component in outwards normal dirn
1126  {
1127  Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1128 
1129  if ( Dist < halfCarTolerance )
1130  {
1131  sd = Dist/Comp ;
1132 
1133  if (sd < snxt)
1134  {
1135  if ( sd < 0 ) { sd = 0; }
1136  zi = p.z() + sd*v.z() ;
1137  if ( std::fabs(zi) <= tolODz )
1138  {
1139  xi = p.x() + sd*v.x() ;
1140  yi = p.y() + sd*v.y() ;
1141  rho2 = xi*xi + yi*yi ;
1142  if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1143  || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1144  && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1145  && (v.x()*cosEPhi + v.y()*sinEPhi >= 0) )
1146  || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1147  && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1148  && (v.x()*cosEPhi + v.y()*sinEPhi < 0) ) )
1149  {
1150  // z and r intersections good
1151  // - check intersecting with correct half-plane
1152  //
1153  if ( (yi*cosCPhi-xi*sinCPhi) >= 0 ) { snxt = sd; }
1154  } //?? >=-halfCarTolerance
1155  }
1156  }
1157  }
1158  } // Comp < 0
1159  } // !fPhiFullTube
1160  if ( snxt<halfCarTolerance ) { snxt=0; }
1161  return snxt ;
1162 }
1163 
1165 //
1166 // Calculate distance to shape from outside, along normalised vector
1167 // - return kInfinity if no intersection, or intersection distance <= tolerance
1168 //
1169 // - Compute the intersection with the z planes
1170 // - if at valid r, phi, return
1171 //
1172 // -> If point is outer outer radius, compute intersection with rmax
1173 // - if at valid phi,z return
1174 //
1175 // -> Compute intersection with inner radius, taking largest +ve root
1176 // - if valid (in z,phi), save intersction
1177 //
1178 // -> If phi segmented, compute intersections with phi half planes
1179 // - return smallest of valid phi intersections and
1180 // inner radius intersection
1181 //
1182 // NOTE:
1183 // - Precalculations for phi trigonometry are Done `just in time'
1184 // - `if valid' implies tolerant checking of intersection points
1185 // Calculate distance (<= actual) to closest surface of shape from outside
1186 // - Calculate distance to z, radial planes
1187 // - Only to phi planes if outside phi extent
1188 // - Return 0 if point inside
1189 
1191 {
1192  G4double safe=0.0, rho, safe1, safe2, safe3 ;
1193  G4double safePhi, cosPsi ;
1194 
1195  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1196  safe1 = fRMin - rho ;
1197  safe2 = rho - fRMax ;
1198  safe3 = std::fabs(p.z()) - fDz ;
1199 
1200  if ( safe1 > safe2 ) { safe = safe1; }
1201  else { safe = safe2; }
1202  if ( safe3 > safe ) { safe = safe3; }
1203 
1204  if ( (!fPhiFullTube) && (rho) )
1205  {
1206  // Psi=angle from central phi to point
1207  //
1208  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1209 
1210  if ( cosPsi < std::cos(fDPhi*0.5) )
1211  {
1212  // Point lies outside phi range
1213 
1214  if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
1215  {
1216  safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ;
1217  }
1218  else
1219  {
1220  safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
1221  }
1222  if ( safePhi > safe ) { safe = safePhi; }
1223  }
1224  }
1225  if ( safe < 0 ) { safe = 0; }
1226  return safe ;
1227 }
1228 
1230 //
1231 // Calculate distance to surface of shape from `inside', allowing for tolerance
1232 // - Only Calc rmax intersection if no valid rmin intersection
1233 
1235  const G4ThreeVector& v,
1236  const G4bool calcNorm,
1237  G4bool *validNorm,
1238  G4ThreeVector *n ) const
1239 {
1240  ESide side=kNull , sider=kNull, sidephi=kNull ;
1241  G4double snxt, srd=kInfinity, sphi=kInfinity, pdist ;
1242  G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
1243 
1244  // Vars for phi intersection:
1245 
1246  G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1247 
1248  // Z plane intersection
1249 
1250  if (v.z() > 0 )
1251  {
1252  pdist = fDz - p.z() ;
1253  if ( pdist > halfCarTolerance )
1254  {
1255  snxt = pdist/v.z() ;
1256  side = kPZ ;
1257  }
1258  else
1259  {
1260  if (calcNorm)
1261  {
1262  *n = G4ThreeVector(0,0,1) ;
1263  *validNorm = true ;
1264  }
1265  return snxt = 0 ;
1266  }
1267  }
1268  else if ( v.z() < 0 )
1269  {
1270  pdist = fDz + p.z() ;
1271 
1272  if ( pdist > halfCarTolerance )
1273  {
1274  snxt = -pdist/v.z() ;
1275  side = kMZ ;
1276  }
1277  else
1278  {
1279  if (calcNorm)
1280  {
1281  *n = G4ThreeVector(0,0,-1) ;
1282  *validNorm = true ;
1283  }
1284  return snxt = 0.0 ;
1285  }
1286  }
1287  else
1288  {
1289  snxt = kInfinity ; // Travel perpendicular to z axis
1290  side = kNull;
1291  }
1292 
1293  // Radial Intersections
1294  //
1295  // Find intersection with cylinders at rmax/rmin
1296  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1297  //
1298  // Intersects with x^2+y^2=R^2
1299  //
1300  // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
1301  //
1302  // t1 t2 t3
1303 
1304  t1 = 1.0 - v.z()*v.z() ; // since v normalised
1305  t2 = p.x()*v.x() + p.y()*v.y() ;
1306  t3 = p.x()*p.x() + p.y()*p.y() ;
1307 
1308  if ( snxt > 10*(fDz+fRMax) ) { roi2 = 2*fRMax*fRMax; }
1309  else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; } // radius^2 on +-fDz
1310 
1311  if ( t1 > 0 ) // Check not parallel
1312  {
1313  // Calculate srd, r exit distance
1314 
1315  if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) )
1316  {
1317  // Delta r not negative => leaving via rmax
1318 
1319  deltaR = t3 - fRMax*fRMax ;
1320 
1321  // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
1322  // - avoid sqrt for efficiency
1323 
1324  if ( deltaR < -kRadTolerance*fRMax )
1325  {
1326  b = t2/t1 ;
1327  c = deltaR/t1 ;
1328  d2 = b*b-c;
1329  if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1330  else { srd = 0.; }
1331  sider = kRMax ;
1332  }
1333  else
1334  {
1335  // On tolerant boundary & heading outwards (or perpendicular to)
1336  // outer radial surface -> leaving immediately
1337 
1338  if ( calcNorm )
1339  {
1340  *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1341  *validNorm = true ;
1342  }
1343  return snxt = 0 ; // Leaving by rmax immediately
1344  }
1345  }
1346  else if ( t2 < 0. ) // i.e. t2 < 0; Possible rmin intersection
1347  {
1348  roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement
1349 
1350  if ( fRMin && (roMin2 < fRMin*(fRMin - kRadTolerance)) )
1351  {
1352  deltaR = t3 - fRMin*fRMin ;
1353  b = t2/t1 ;
1354  c = deltaR/t1 ;
1355  d2 = b*b - c ;
1356 
1357  if ( d2 >= 0 ) // Leaving via rmin
1358  {
1359  // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1360  // - avoid sqrt for efficiency
1361 
1362  if (deltaR > kRadTolerance*fRMin)
1363  {
1364  srd = c/(-b+std::sqrt(d2));
1365  sider = kRMin ;
1366  }
1367  else
1368  {
1369  if ( calcNorm ) { *validNorm = false; } // Concave side
1370  return snxt = 0.0;
1371  }
1372  }
1373  else // No rmin intersect -> must be rmax intersect
1374  {
1375  deltaR = t3 - fRMax*fRMax ;
1376  c = deltaR/t1 ;
1377  d2 = b*b-c;
1378  if( d2 >=0. )
1379  {
1380  srd = -b + std::sqrt(d2) ;
1381  sider = kRMax ;
1382  }
1383  else // Case: On the border+t2<kRadTolerance
1384  // (v is perpendicular to the surface)
1385  {
1386  if (calcNorm)
1387  {
1388  *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1389  *validNorm = true ;
1390  }
1391  return snxt = 0.0;
1392  }
1393  }
1394  }
1395  else if ( roi2 > fRMax*(fRMax + kRadTolerance) )
1396  // No rmin intersect -> must be rmax intersect
1397  {
1398  deltaR = t3 - fRMax*fRMax ;
1399  b = t2/t1 ;
1400  c = deltaR/t1;
1401  d2 = b*b-c;
1402  if( d2 >= 0 )
1403  {
1404  srd = -b + std::sqrt(d2) ;
1405  sider = kRMax ;
1406  }
1407  else // Case: On the border+t2<kRadTolerance
1408  // (v is perpendicular to the surface)
1409  {
1410  if (calcNorm)
1411  {
1412  *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1413  *validNorm = true ;
1414  }
1415  return snxt = 0.0;
1416  }
1417  }
1418  }
1419 
1420  // Phi Intersection
1421 
1422  if ( !fPhiFullTube )
1423  {
1424  // add angle calculation with correction
1425  // of the difference in domain of atan2 and Sphi
1426  //
1427  vphi = std::atan2(v.y(),v.x()) ;
1428 
1429  if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1430  else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1431 
1432 
1433  if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1434  {
1435  // pDist -ve when inside
1436 
1437  pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1438  pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1439 
1440  // Comp -ve when in direction of outwards normal
1441 
1442  compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1443  compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1444 
1445  sidephi = kNull;
1446 
1447  if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1448  && (pDistE <= halfCarTolerance) ) )
1449  || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
1450  && (pDistE > halfCarTolerance) ) ) )
1451  {
1452  // Inside both phi *full* planes
1453 
1454  if ( compS < 0 )
1455  {
1456  sphi = pDistS/compS ;
1457 
1458  if (sphi >= -halfCarTolerance)
1459  {
1460  xi = p.x() + sphi*v.x() ;
1461  yi = p.y() + sphi*v.y() ;
1462 
1463  // Check intersecting with correct half-plane
1464  // (if not -> no intersect)
1465  //
1466  if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
1467  {
1468  sidephi = kSPhi;
1469  if (((fSPhi-halfAngTolerance)<=vphi)
1470  &&((fSPhi+fDPhi+halfAngTolerance)>=vphi))
1471  {
1472  sphi = kInfinity;
1473  }
1474  }
1475  else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1476  {
1477  sphi = kInfinity ;
1478  }
1479  else
1480  {
1481  sidephi = kSPhi ;
1482  if ( pDistS > -halfCarTolerance )
1483  {
1484  sphi = 0.0 ; // Leave by sphi immediately
1485  }
1486  }
1487  }
1488  else
1489  {
1490  sphi = kInfinity ;
1491  }
1492  }
1493  else
1494  {
1495  sphi = kInfinity ;
1496  }
1497 
1498  if ( compE < 0 )
1499  {
1500  sphi2 = pDistE/compE ;
1501 
1502  // Only check further if < starting phi intersection
1503  //
1504  if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1505  {
1506  xi = p.x() + sphi2*v.x() ;
1507  yi = p.y() + sphi2*v.y() ;
1508 
1509  if ((std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance))
1510  {
1511  // Leaving via ending phi
1512  //
1513  if( !((fSPhi-halfAngTolerance <= vphi)
1514  &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
1515  {
1516  sidephi = kEPhi ;
1517  if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1518  else { sphi = 0.0 ; }
1519  }
1520  }
1521  else // Check intersecting with correct half-plane
1522 
1523  if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1524  {
1525  // Leaving via ending phi
1526  //
1527  sidephi = kEPhi ;
1528  if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1529  else { sphi = 0.0 ; }
1530  }
1531  }
1532  }
1533  }
1534  else
1535  {
1536  sphi = kInfinity ;
1537  }
1538  }
1539  else
1540  {
1541  // On z axis + travel not || to z axis -> if phi of vector direction
1542  // within phi of shape, Step limited by rmax, else Step =0
1543 
1544  if ( (fSPhi - halfAngTolerance <= vphi)
1545  && (vphi <= fSPhi + fDPhi + halfAngTolerance ) )
1546  {
1547  sphi = kInfinity ;
1548  }
1549  else
1550  {
1551  sidephi = kSPhi ; // arbitrary
1552  sphi = 0.0 ;
1553  }
1554  }
1555  if (sphi < snxt) // Order intersecttions
1556  {
1557  snxt = sphi ;
1558  side = sidephi ;
1559  }
1560  }
1561  if (srd < snxt) // Order intersections
1562  {
1563  snxt = srd ;
1564  side = sider ;
1565  }
1566  }
1567  if (calcNorm)
1568  {
1569  switch(side)
1570  {
1571  case kRMax:
1572  // Note: returned vector not normalised
1573  // (divide by fRMax for unit vector)
1574  //
1575  xi = p.x() + snxt*v.x() ;
1576  yi = p.y() + snxt*v.y() ;
1577  *n = G4ThreeVector(xi/fRMax,yi/fRMax,0) ;
1578  *validNorm = true ;
1579  break ;
1580 
1581  case kRMin:
1582  *validNorm = false ; // Rmin is inconvex
1583  break ;
1584 
1585  case kSPhi:
1586  if ( fDPhi <= pi )
1587  {
1588  *n = G4ThreeVector(sinSPhi,-cosSPhi,0) ;
1589  *validNorm = true ;
1590  }
1591  else
1592  {
1593  *validNorm = false ;
1594  }
1595  break ;
1596 
1597  case kEPhi:
1598  if (fDPhi <= pi)
1599  {
1600  *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ;
1601  *validNorm = true ;
1602  }
1603  else
1604  {
1605  *validNorm = false ;
1606  }
1607  break ;
1608 
1609  case kPZ:
1610  *n = G4ThreeVector(0,0,1) ;
1611  *validNorm = true ;
1612  break ;
1613 
1614  case kMZ:
1615  *n = G4ThreeVector(0,0,-1) ;
1616  *validNorm = true ;
1617  break ;
1618 
1619  default:
1620  G4cout << G4endl ;
1621  DumpInfo();
1622  std::ostringstream message;
1623  G4int oldprc = message.precision(16);
1624  message << "Undefined side for valid surface normal to solid."
1625  << G4endl
1626  << "Position:" << G4endl << G4endl
1627  << "p.x() = " << p.x()/mm << " mm" << G4endl
1628  << "p.y() = " << p.y()/mm << " mm" << G4endl
1629  << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
1630  << "Direction:" << G4endl << G4endl
1631  << "v.x() = " << v.x() << G4endl
1632  << "v.y() = " << v.y() << G4endl
1633  << "v.z() = " << v.z() << G4endl << G4endl
1634  << "Proposed distance :" << G4endl << G4endl
1635  << "snxt = " << snxt/mm << " mm" << G4endl ;
1636  message.precision(oldprc) ;
1637  G4Exception("G4Tubs::DistanceToOut(p,v,..)", "GeomSolids1002",
1638  JustWarning, message);
1639  break ;
1640  }
1641  }
1642  if ( snxt<halfCarTolerance ) { snxt=0 ; }
1643 
1644  return snxt ;
1645 }
1646 
1648 //
1649 // Calculate distance (<=actual) to closest surface of shape from inside
1650 
1652 {
1653  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi ;
1654  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1655 
1656 #ifdef G4CSGDEBUG
1657  if( Inside(p) == kOutside )
1658  {
1659  G4int oldprc = G4cout.precision(16) ;
1660  G4cout << G4endl ;
1661  DumpInfo();
1662  G4cout << "Position:" << G4endl << G4endl ;
1663  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1664  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1665  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1666  G4cout.precision(oldprc) ;
1667  G4Exception("G4Tubs::DistanceToOut(p)", "GeomSolids1002",
1668  JustWarning, "Point p is outside !?");
1669  }
1670 #endif
1671 
1672  if ( fRMin )
1673  {
1674  safeR1 = rho - fRMin ;
1675  safeR2 = fRMax - rho ;
1676 
1677  if ( safeR1 < safeR2 ) { safe = safeR1 ; }
1678  else { safe = safeR2 ; }
1679  }
1680  else
1681  {
1682  safe = fRMax - rho ;
1683  }
1684  safeZ = fDz - std::fabs(p.z()) ;
1685 
1686  if ( safeZ < safe ) { safe = safeZ ; }
1687 
1688  // Check if phi divided, Calc distances closest phi plane
1689  //
1690  if ( !fPhiFullTube )
1691  {
1692  if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
1693  {
1694  safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
1695  }
1696  else
1697  {
1698  safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
1699  }
1700  if (safePhi < safe) { safe = safePhi ; }
1701  }
1702  if ( safe < 0 ) { safe = 0 ; }
1703 
1704  return safe ;
1705 }
1706 
1708 //
1709 // Create a List containing the transformed vertices
1710 // Ordering [0-3] -fDz cross section
1711 // [4-7] +fDz cross section such that [0] is below [4],
1712 // [1] below [5] etc.
1713 // Note:
1714 // Caller has deletion resposibility
1715 // Potential improvement: For last slice, use actual ending angle
1716 // to avoid rounding error problems.
1717 
1720 {
1721  G4ThreeVectorList* vertices ;
1722  G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
1723  G4double meshAngle, meshRMax, crossAngle,
1724  cosCrossAngle, sinCrossAngle, sAngle;
1725  G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ;
1726  G4int crossSection, noCrossSections;
1727 
1728  // Compute no of cross-sections necessary to mesh tube
1729  //
1730  noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;
1731 
1732  if ( noCrossSections < kMinMeshSections )
1733  {
1734  noCrossSections = kMinMeshSections ;
1735  }
1736  else if (noCrossSections>kMaxMeshSections)
1737  {
1738  noCrossSections = kMaxMeshSections ;
1739  }
1740  // noCrossSections = 4 ;
1741 
1742  meshAngle = fDPhi/(noCrossSections - 1) ;
1743  // meshAngle = fDPhi/(noCrossSections) ;
1744 
1745  meshRMax = (fRMax+100*kCarTolerance)/std::cos(meshAngle*0.5) ;
1746  meshRMin = fRMin - 100*kCarTolerance ;
1747 
1748  // If complete in phi, set start angle such that mesh will be at fRMax
1749  // on the x axis. Will give better extent calculations when not rotated.
1750 
1751  if (fPhiFullTube && (fSPhi == 0) ) { sAngle = -meshAngle*0.5 ; }
1752  else { sAngle = fSPhi ; }
1753 
1754  vertices = new G4ThreeVectorList();
1755 
1756  if ( vertices )
1757  {
1758  vertices->reserve(noCrossSections*4);
1759  for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ )
1760  {
1761  // Compute coordinates of cross section at section crossSection
1762 
1763  crossAngle = sAngle + crossSection*meshAngle ;
1764  cosCrossAngle = std::cos(crossAngle) ;
1765  sinCrossAngle = std::sin(crossAngle) ;
1766 
1767  rMaxX = meshRMax*cosCrossAngle ;
1768  rMaxY = meshRMax*sinCrossAngle ;
1769 
1770  if(meshRMin <= 0.0)
1771  {
1772  rMinX = 0.0 ;
1773  rMinY = 0.0 ;
1774  }
1775  else
1776  {
1777  rMinX = meshRMin*cosCrossAngle ;
1778  rMinY = meshRMin*sinCrossAngle ;
1779  }
1780  vertex0 = G4ThreeVector(rMinX,rMinY,-fDz) ;
1781  vertex1 = G4ThreeVector(rMaxX,rMaxY,-fDz) ;
1782  vertex2 = G4ThreeVector(rMaxX,rMaxY,+fDz) ;
1783  vertex3 = G4ThreeVector(rMinX,rMinY,+fDz) ;
1784 
1785  vertices->push_back(pTransform.TransformPoint(vertex0)) ;
1786  vertices->push_back(pTransform.TransformPoint(vertex1)) ;
1787  vertices->push_back(pTransform.TransformPoint(vertex2)) ;
1788  vertices->push_back(pTransform.TransformPoint(vertex3)) ;
1789  }
1790  }
1791  else
1792  {
1793  DumpInfo();
1794  G4Exception("G4Tubs::CreateRotatedVertices()",
1795  "GeomSolids0003", FatalException,
1796  "Error in allocation of vertices. Out of memory !");
1797  }
1798  return vertices ;
1799 }
1800 
1802 //
1803 // Stream object contents to an output stream
1804 
1806 {
1807  return G4String("G4Tubs");
1808 }
1809 
1811 //
1812 // Make a clone of the object
1813 //
1815 {
1816  return new G4Tubs(*this);
1817 }
1818 
1820 //
1821 // Stream object contents to an output stream
1822 
1823 std::ostream& G4Tubs::StreamInfo( std::ostream& os ) const
1824 {
1825  G4int oldprc = os.precision(16);
1826  os << "-----------------------------------------------------------\n"
1827  << " *** Dump for solid - " << GetName() << " ***\n"
1828  << " ===================================================\n"
1829  << " Solid type: G4Tubs\n"
1830  << " Parameters: \n"
1831  << " inner radius : " << fRMin/mm << " mm \n"
1832  << " outer radius : " << fRMax/mm << " mm \n"
1833  << " half length Z: " << fDz/mm << " mm \n"
1834  << " starting phi : " << fSPhi/degree << " degrees \n"
1835  << " delta phi : " << fDPhi/degree << " degrees \n"
1836  << "-----------------------------------------------------------\n";
1837  os.precision(oldprc);
1838 
1839  return os;
1840 }
1841 
1843 //
1844 // GetPointOnSurface
1845 
1847 {
1848  G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1849  aOne, aTwo, aThr, aFou;
1850  G4double rRand;
1851 
1852  aOne = 2.*fDz*fDPhi*fRMax;
1853  aTwo = 2.*fDz*fDPhi*fRMin;
1854  aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin);
1855  aFou = 2.*fDz*(fRMax-fRMin);
1856 
1857  phi = RandFlat::shoot(fSPhi, fSPhi+fDPhi);
1858  cosphi = std::cos(phi);
1859  sinphi = std::sin(phi);
1860 
1861  rRand = GetRadiusInRing(fRMin,fRMax);
1862 
1863  if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; }
1864 
1865  chose = RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou);
1866 
1867  if( (chose >=0) && (chose < aOne) )
1868  {
1869  xRand = fRMax*cosphi;
1870  yRand = fRMax*sinphi;
1871  zRand = RandFlat::shoot(-1.*fDz,fDz);
1872  return G4ThreeVector (xRand, yRand, zRand);
1873  }
1874  else if( (chose >= aOne) && (chose < aOne + aTwo) )
1875  {
1876  xRand = fRMin*cosphi;
1877  yRand = fRMin*sinphi;
1878  zRand = RandFlat::shoot(-1.*fDz,fDz);
1879  return G4ThreeVector (xRand, yRand, zRand);
1880  }
1881  else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1882  {
1883  xRand = rRand*cosphi;
1884  yRand = rRand*sinphi;
1885  zRand = fDz;
1886  return G4ThreeVector (xRand, yRand, zRand);
1887  }
1888  else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1889  {
1890  xRand = rRand*cosphi;
1891  yRand = rRand*sinphi;
1892  zRand = -1.*fDz;
1893  return G4ThreeVector (xRand, yRand, zRand);
1894  }
1895  else if( (chose >= aOne + aTwo + 2.*aThr)
1896  && (chose < aOne + aTwo + 2.*aThr + aFou) )
1897  {
1898  xRand = rRand*std::cos(fSPhi);
1899  yRand = rRand*std::sin(fSPhi);
1900  zRand = RandFlat::shoot(-1.*fDz,fDz);
1901  return G4ThreeVector (xRand, yRand, zRand);
1902  }
1903  else
1904  {
1905  xRand = rRand*std::cos(fSPhi+fDPhi);
1906  yRand = rRand*std::sin(fSPhi+fDPhi);
1907  zRand = RandFlat::shoot(-1.*fDz,fDz);
1908  return G4ThreeVector (xRand, yRand, zRand);
1909  }
1910 }
1911 
1913 //
1914 // Methods for visualisation
1915 
1917 {
1918  scene.AddSolid (*this) ;
1919 }
1920 
1922 {
1923  return new G4PolyhedronTubs (fRMin, fRMax, fDz, fSPhi, fDPhi) ;
1924 }
1925 #endif
G4String GetName() const
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:347
ThreeVector shoot(const G4int Ap, const G4int Af)
G4bool fPhiFullTube
Definition: G4Tubs.hh:225
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
G4double fDPhi
Definition: G4Tubs.hh:216
G4double halfCarTolerance
Definition: G4Tubs.hh:229
G4AffineTransform Inverse() const
G4bool IsYLimited() const
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:101
G4Tubs(const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi)
Definition: G4Tubs.cc:86
const G4double pi
Definition: G4Tubs.hh:84
G4double cosEPhi
Definition: G4Tubs.hh:220
G4bool IsRotated() const
G4ThreeVector NetTranslation() const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Tubs.cc:814
G4bool IsXLimited() const
G4double a
Definition: TRTMaterials.hh:39
const G4int kMinMeshSections
Definition: meshdefs.hh:45
G4ThreeVector GetPointOnSurface() const
Definition: G4Tubs.cc:1846
virtual void AddSolid(const G4Box &)=0
virtual G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4Tubs.cc:682
int G4int
Definition: G4Types.hh:78
G4double cosHDPhiOT
Definition: G4Tubs.hh:220
void DumpInfo() const
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Tubs.cc:1823
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Tubs.cc:1916
G4double GetMaxXExtent() const
G4double sinSPhi
Definition: G4Tubs.hh:220
G4double GetMinZExtent() const
G4GeometryType GetEntityType() const
Definition: G4Tubs.cc:1805
G4GLOB_DLL std::ostream G4cout
G4double fRMin
Definition: G4Tubs.hh:216
G4Polyhedron * CreatePolyhedron() const
Definition: G4Tubs.cc:1921
G4double GetRadialTolerance() const
bool G4bool
Definition: G4Types.hh:79
G4VSolid * Clone() const
Definition: G4Tubs.cc:1814
G4double cosCPhi
Definition: G4Tubs.hh:220
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
G4double fSPhi
Definition: G4Tubs.hh:216
const G4int n
virtual ~G4Tubs()
Definition: G4Tubs.cc:140
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Tubs.cc:210
G4double sinEPhi
Definition: G4Tubs.hh:220
G4double halfAngTolerance
Definition: G4Tubs.hh:229
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Tubs.cc:1234
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4double GetMinXExtent() const
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4double GetMaxZExtent() const
EInside Inside(const G4ThreeVector &p) const
Definition: G4Tubs.cc:424
G4double halfRadTolerance
Definition: G4Tubs.hh:229
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
static const double degree
Definition: G4SIunits.hh:125
G4double fRMax
Definition: G4Tubs.hh:216
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Tubs.cc:591
#define G4endl
Definition: G4ios.hh:61
G4double GetMaxYExtent() const
G4double kCarTolerance
Definition: G4VSolid.hh:305
G4double kRadTolerance
Definition: G4Tubs.hh:212
G4double fDz
Definition: G4Tubs.hh:216
G4double cosSPhi
Definition: G4Tubs.hh:220
double G4double
Definition: G4Types.hh:76
G4double GetMaxExtent(const EAxis pAxis) const
static const G4double d2
void CheckPhiAngles(G4double sPhi, G4double dPhi)
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:82
G4bool IsZLimited() const
G4double cosHDPhiIT
Definition: G4Tubs.hh:220
static const double mm
Definition: G4SIunits.hh:102
G4double GetAngularTolerance() const
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:378
G4double sinCPhi
Definition: G4Tubs.hh:220
G4double GetMinExtent(const EAxis pAxis) const
const G4double kMeshAngleDefault
Definition: meshdefs.hh:42
static G4GeometryTolerance * GetInstance()
G4double kAngTolerance
Definition: G4Tubs.hh:212
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Tubs.cc:199
G4Tubs & operator=(const G4Tubs &rhs)
Definition: G4Tubs.cc:167
const G4int kMaxMeshSections
Definition: meshdefs.hh:46
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
Definition: G4Tubs.cc:1719