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