Geant4  10.02.p01
G4Torus.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: G4Torus.cc 92393 2015-08-31 14:07:30Z gcosmo $
28 //
29 //
30 // class G4Torus
31 //
32 // Implementation
33 //
34 // 05.04.12 M.Kelsey: Use sqrt(r) in GetPointOnSurface() for uniform points
35 // 02.10.07 T.Nikitina: Bug fixed in SolveNumericJT(), b.969:segmentation fault.
36 // rootsrefined is used only if the number of refined roots
37 // is the same as for primary roots.
38 // 02.10.07 T.Nikitina: Bug fixed in CalculateExtent() for case of non-rotated
39 // full-phi torus:protect against negative value for sqrt,
40 // correct formula for delta.
41 // 20.11.05 V.Grichine: Bug fixed in Inside(p) for phi sections, b.810
42 // 25.08.05 O.Link: new methods for DistanceToIn/Out using JTPolynomialSolver
43 // 07.06.05 V.Grichine: SurfaceNormal(p) for rho=0, Constructor as G4Cons
44 // 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
45 // 18.03.04 V.Grichine: bug fixed in DistanceToIn(p)
46 // 11.01.01 E.Medernach: Use G4PolynomialSolver to find roots
47 // 03.10.00 E.Medernach: SafeNewton added
48 // 31.08.00 E.Medernach: numerical computation of roots wuth bounding
49 // volume technique
50 // 26.05.00 V.Grichine: new fuctions developed by O.Cremonesi were added
51 // 06.03.00 V.Grichine: modifications in Distance ToOut(p,v,...)
52 // 19.11.99 V.Grichine: side = kNull in Distance ToOut(p,v,...)
53 // 09.10.98 V.Grichine: modifications in Distance ToOut(p,v,...)
54 // 30.10.96 V.Grichine: first implementation with G4Tubs elements in Fs
55 //
56 
57 #include "G4Torus.hh"
58 
59 #if !(defined(G4GEOM_USE_UTORUS) && defined(G4GEOM_USE_SYS_USOLIDS))
60 
61 #include "G4VoxelLimits.hh"
62 #include "G4AffineTransform.hh"
63 #include "G4GeometryTolerance.hh"
64 #include "G4JTPolynomialSolver.hh"
65 
66 #include "G4VPVParameterisation.hh"
67 
68 #include "meshdefs.hh"
69 
70 #include "Randomize.hh"
71 
72 #include "G4VGraphicsScene.hh"
73 #include "G4Polyhedron.hh"
74 
75 using namespace CLHEP;
76 
78 //
79 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
80 // - note if pdphi>2PI then reset to 2PI
81 
82 G4Torus::G4Torus( const G4String &pName,
83  G4double pRmin,
84  G4double pRmax,
85  G4double pRtor,
86  G4double pSPhi,
87  G4double pDPhi)
88  : G4CSGSolid(pName)
89 {
90  SetAllParameters(pRmin, pRmax, pRtor, pSPhi, pDPhi);
91 }
92 
94 //
95 //
96 
97 void
99  G4double pRmax,
100  G4double pRtor,
101  G4double pSPhi,
102  G4double pDPhi )
103 {
104  const G4double fEpsilon = 4.e-11; // relative tolerance of radii
105 
106  fCubicVolume = 0.;
107  fSurfaceArea = 0.;
108  fRebuildPolyhedron = true;
109 
112 
115 
116  if ( pRtor >= pRmax+1.e3*kCarTolerance ) // Check swept radius, as in G4Cons
117  {
118  fRtor = pRtor ;
119  }
120  else
121  {
122  std::ostringstream message;
123  message << "Invalid swept radius for Solid: " << GetName() << G4endl
124  << " pRtor = " << pRtor << ", pRmax = " << pRmax;
125  G4Exception("G4Torus::SetAllParameters()",
126  "GeomSolids0002", FatalException, message);
127  }
128 
129  // Check radii, as in G4Cons
130  //
131  if ( pRmin < pRmax - 1.e2*kCarTolerance && pRmin >= 0 )
132  {
133  if (pRmin >= 1.e2*kCarTolerance) { fRmin = pRmin ; }
134  else { fRmin = 0.0 ; }
135  fRmax = pRmax ;
136  }
137  else
138  {
139  std::ostringstream message;
140  message << "Invalid values of radii for Solid: " << GetName() << G4endl
141  << " pRmin = " << pRmin << ", pRmax = " << pRmax;
142  G4Exception("G4Torus::SetAllParameters()",
143  "GeomSolids0002", FatalException, message);
144  }
145 
146  // Relative tolerances
147  //
149  ? 0.5*std::max( kRadTolerance, fEpsilon*(fRtor-fRmin )) : 0;
150  fRmaxTolerance = 0.5*std::max( kRadTolerance, fEpsilon*(fRtor+fRmax) );
151 
152  // Check angles
153  //
154  if ( pDPhi >= twopi ) { fDPhi = twopi ; }
155  else
156  {
157  if (pDPhi > 0) { fDPhi = pDPhi ; }
158  else
159  {
160  std::ostringstream message;
161  message << "Invalid Z delta-Phi for Solid: " << GetName() << G4endl
162  << " pDPhi = " << pDPhi;
163  G4Exception("G4Torus::SetAllParameters()",
164  "GeomSolids0002", FatalException, message);
165  }
166  }
167 
168  // Ensure psphi in 0-2PI or -2PI-0 range if shape crosses 0
169  //
170  fSPhi = pSPhi;
171 
172  if (fSPhi < 0) { fSPhi = twopi-std::fmod(std::fabs(fSPhi),twopi) ; }
173  else { fSPhi = std::fmod(fSPhi,twopi) ; }
174 
175  if (fSPhi+fDPhi > twopi) { fSPhi-=twopi ; }
176 }
177 
179 //
180 // Fake default constructor - sets only member data and allocates memory
181 // for usage restricted to object persistency.
182 //
183 G4Torus::G4Torus( __void__& a )
184  : G4CSGSolid(a), fRmin(0.), fRmax(0.), fRtor(0.), fSPhi(0.),
185  fDPhi(0.), fRminTolerance(0.), fRmaxTolerance(0. ),
186  kRadTolerance(0.), kAngTolerance(0.),
187  halfCarTolerance(0.), halfAngTolerance(0.)
188 {
189 }
190 
192 //
193 // Destructor
194 
196 {}
197 
199 //
200 // Copy constructor
201 
203  : G4CSGSolid(rhs), fRmin(rhs.fRmin),fRmax(rhs.fRmax),
204  fRtor(rhs.fRtor),fSPhi(rhs.fSPhi),fDPhi(rhs.fDPhi),
205  fRminTolerance(rhs.fRminTolerance), fRmaxTolerance(rhs.fRmaxTolerance),
206  kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
207  halfCarTolerance(rhs.halfCarTolerance),
208  halfAngTolerance(rhs.halfAngTolerance)
209 {
210 }
211 
213 //
214 // Assignment operator
215 
217 {
218  // Check assignment to self
219  //
220  if (this == &rhs) { return *this; }
221 
222  // Copy base class data
223  //
225 
226  // Copy data
227  //
228  fRmin = rhs.fRmin; fRmax = rhs.fRmax;
229  fRtor = rhs.fRtor; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
234 
235  return *this;
236 }
237 
239 //
240 // Dispatch to parameterisation for replication mechanism dimension
241 // computation & modification.
242 
244  const G4int n,
245  const G4VPhysicalVolume* pRep )
246 {
247  p->ComputeDimensions(*this,n,pRep);
248 }
249 
250 
251 
253 //
254 // Calculate the real roots to torus surface.
255 // Returns negative solutions as well.
256 
258  const G4ThreeVector& v,
259  G4double r,
260  std::vector<G4double>& roots ) const
261 {
262 
263  G4int i, num ;
264  G4double c[5], srd[4], si[4] ;
265 
266  G4double Rtor2 = fRtor*fRtor, r2 = r*r ;
267 
268  G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
269  G4double pRad2 = p.x()*p.x() + p.y()*p.y() + p.z()*p.z() ;
270 
271  c[0] = 1.0 ;
272  c[1] = 4*pDotV ;
273  c[2] = 2*(pRad2 + 2*pDotV*pDotV - Rtor2 - r2 + 2*Rtor2*v.z()*v.z()) ;
274  c[3] = 4*(pDotV*(pRad2 - Rtor2 - r2) + 2*Rtor2*p.z()*v.z()) ;
275  c[4] = pRad2*pRad2 - 2*pRad2*(Rtor2+r2)
276  + 4*Rtor2*p.z()*p.z() + (Rtor2-r2)*(Rtor2-r2) ;
277 
278  G4JTPolynomialSolver torusEq;
279 
280  num = torusEq.FindRoots( c, 4, srd, si );
281 
282  for ( i = 0; i < num; i++ )
283  {
284  if( si[i] == 0. ) { roots.push_back(srd[i]) ; } // store real roots
285  }
286 
287  std::sort(roots.begin() , roots.end() ) ; // sorting with <
288 }
289 
291 //
292 // Interface for DistanceToIn and DistanceToOut.
293 // Calls TorusRootsJT and returns the smalles possible distance to
294 // the surface.
295 // Attention: Difference in DistanceToIn/Out for points p on the surface.
296 
298  const G4ThreeVector& v,
299  G4double r,
300  G4bool IsDistanceToIn ) const
301 {
302  G4double bigdist = 10*mm ;
303  G4double tmin = kInfinity ;
304  G4double t, scal ;
305 
306  // calculate the distances to the intersections with the Torus
307  // from a given point p and direction v.
308  //
309  std::vector<G4double> roots ;
310  std::vector<G4double> rootsrefined ;
311  TorusRootsJT(p,v,r,roots) ;
312 
313  G4ThreeVector ptmp ;
314 
315  // determine the smallest non-negative solution
316  //
317  for ( size_t k = 0 ; k<roots.size() ; k++ )
318  {
319  t = roots[k] ;
320 
321  if ( t < -halfCarTolerance ) { continue ; } // skip negative roots
322 
323  if ( t > bigdist && t<kInfinity ) // problem with big distances
324  {
325  ptmp = p + t*v ;
326  TorusRootsJT(ptmp,v,r,rootsrefined) ;
327  if ( rootsrefined.size()==roots.size() )
328  {
329  t = t + rootsrefined[k] ;
330  }
331  }
332 
333  ptmp = p + t*v ; // calculate the position of the proposed intersection
334 
335  G4double theta = std::atan2(ptmp.y(),ptmp.x());
336 
337  if ( fSPhi >= 0 )
338  {
339  if ( theta < - halfAngTolerance ) { theta += twopi; }
340  if ( (std::fabs(theta) < halfAngTolerance)
341  && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
342  {
343  theta += twopi ; // 0 <= theta < 2pi
344  }
345  }
346  if ((fSPhi <= -pi )&&(theta>halfAngTolerance)) { theta = theta-twopi; }
347 
348  // We have to verify if this root is inside the region between
349  // fSPhi and fSPhi + fDPhi
350  //
351  if ( (theta - fSPhi >= - halfAngTolerance)
352  && (theta - (fSPhi + fDPhi) <= halfAngTolerance) )
353  {
354  // check if P is on the surface, and called from DistanceToIn
355  // DistanceToIn has to return 0.0 if particle is going inside the solid
356 
357  if ( IsDistanceToIn == true )
358  {
359  if (std::fabs(t) < halfCarTolerance )
360  {
361  // compute scalar product at position p : v.n
362  // ( n taken from SurfaceNormal, not normalized )
363 
364  scal = v* G4ThreeVector( p.x()*(1-fRtor/std::sqrt(p.x()*p.x()
365  + p.y()*p.y())),
366  p.y()*(1-fRtor/std::sqrt(p.x()*p.x()
367  + p.y()*p.y())),
368  p.z() );
369 
370  // change sign in case of inner radius
371  //
372  if ( r == GetRmin() ) { scal = -scal ; }
373  if ( scal < 0 ) { return 0.0 ; }
374  }
375  }
376 
377  // check if P is on the surface, and called from DistanceToOut
378  // DistanceToIn has to return 0.0 if particle is leaving the solid
379 
380  if ( IsDistanceToIn == false )
381  {
382  if (std::fabs(t) < halfCarTolerance )
383  {
384  // compute scalar product at position p : v.n
385  //
386  scal = v* G4ThreeVector( p.x()*(1-fRtor/std::sqrt(p.x()*p.x()
387  + p.y()*p.y())),
388  p.y()*(1-fRtor/std::sqrt(p.x()*p.x()
389  + p.y()*p.y())),
390  p.z() );
391 
392  // change sign in case of inner radius
393  //
394  if ( r == GetRmin() ) { scal = -scal ; }
395  if ( scal > 0 ) { return 0.0 ; }
396  }
397  }
398 
399  // check if distance is larger than 1/2 kCarTolerance
400  //
401  if( t > halfCarTolerance )
402  {
403  tmin = t ;
404  return tmin ;
405  }
406  }
407  }
408 
409  return tmin;
410 }
411 
413 //
414 // Calculate extent under transform and specified limit
415 
417  const G4VoxelLimits& pVoxelLimit,
418  const G4AffineTransform& pTransform,
419  G4double& pMin, G4double& pMax) const
420 {
421  if ((!pTransform.IsRotated()) && (fDPhi==twopi) && (fRmin==0))
422  {
423  // Special case handling for unrotated solid torus
424  // Compute x/y/z mins and maxs for bounding box respecting limits,
425  // with early returns if outside limits. Then switch() on pAxis,
426  // and compute exact x and y limit for x/y case
427 
428  G4double xoffset,xMin,xMax;
429  G4double yoffset,yMin,yMax;
430  G4double zoffset,zMin,zMax;
431 
432  G4double RTorus,delta,diff1,diff2,maxDiff,newMin,newMax;
433  G4double xoff1,xoff2,yoff1,yoff2;
434 
435  xoffset = pTransform.NetTranslation().x();
436  xMin = xoffset - fRmax - fRtor ;
437  xMax = xoffset + fRmax + fRtor ;
438 
439  if (pVoxelLimit.IsXLimited())
440  {
441  if ( (xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance)
442  || (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance) )
443  return false ;
444  else
445  {
446  if (xMin < pVoxelLimit.GetMinXExtent())
447  {
448  xMin = pVoxelLimit.GetMinXExtent() ;
449  }
450  if (xMax > pVoxelLimit.GetMaxXExtent())
451  {
452  xMax = pVoxelLimit.GetMaxXExtent() ;
453  }
454  }
455  }
456  yoffset = pTransform.NetTranslation().y();
457  yMin = yoffset - fRmax - fRtor ;
458  yMax = yoffset + fRmax + fRtor ;
459 
460  if (pVoxelLimit.IsYLimited())
461  {
462  if ( (yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance)
463  || (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance) )
464  {
465  return false ;
466  }
467  else
468  {
469  if (yMin < pVoxelLimit.GetMinYExtent() )
470  {
471  yMin = pVoxelLimit.GetMinYExtent() ;
472  }
473  if (yMax > pVoxelLimit.GetMaxYExtent() )
474  {
475  yMax = pVoxelLimit.GetMaxYExtent() ;
476  }
477  }
478  }
479  zoffset = pTransform.NetTranslation().z() ;
480  zMin = zoffset - fRmax ;
481  zMax = zoffset + fRmax ;
482 
483  if (pVoxelLimit.IsZLimited())
484  {
485  if ( (zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance)
486  || (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance) )
487  {
488  return false ;
489  }
490  else
491  {
492  if (zMin < pVoxelLimit.GetMinZExtent() )
493  {
494  zMin = pVoxelLimit.GetMinZExtent() ;
495  }
496  if (zMax > pVoxelLimit.GetMaxZExtent() )
497  {
498  zMax = pVoxelLimit.GetMaxZExtent() ;
499  }
500  }
501  }
502 
503  // Known to cut cylinder
504 
505  switch (pAxis)
506  {
507  case kXAxis:
508  yoff1=yoffset-yMin;
509  yoff2=yMax-yoffset;
510  if ( yoff1 >= 0 && yoff2 >= 0 )
511  {
512  // Y limits cross max/min x => no change
513  //
514  pMin = xMin ;
515  pMax = xMax ;
516  }
517  else
518  {
519  // Y limits don't cross max/min x => compute max delta x,
520  // hence new mins/maxs
521  //
522 
523  RTorus=fRmax+fRtor;
524  delta = RTorus*RTorus - yoff1*yoff1;
525  diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
526  delta = RTorus*RTorus - yoff2*yoff2;
527  diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
528  maxDiff = (diff1 > diff2) ? diff1:diff2 ;
529  newMin = xoffset - maxDiff ;
530  newMax = xoffset + maxDiff ;
531  pMin = (newMin < xMin) ? xMin : newMin ;
532  pMax = (newMax > xMax) ? xMax : newMax ;
533  }
534  break;
535 
536  case kYAxis:
537  xoff1 = xoffset - xMin ;
538  xoff2 = xMax - xoffset ;
539  if (xoff1 >= 0 && xoff2 >= 0 )
540  {
541  // X limits cross max/min y => no change
542  //
543  pMin = yMin ;
544  pMax = yMax ;
545  }
546  else
547  {
548  // X limits don't cross max/min y => compute max delta y,
549  // hence new mins/maxs
550  //
551  RTorus=fRmax+fRtor;
552  delta = RTorus*RTorus - xoff1*xoff1;
553  diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
554  delta = RTorus*RTorus - xoff2*xoff2;
555  diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
556  maxDiff = (diff1 > diff2) ? diff1 : diff2 ;
557  newMin = yoffset - maxDiff ;
558  newMax = yoffset + maxDiff ;
559  pMin = (newMin < yMin) ? yMin : newMin ;
560  pMax = (newMax > yMax) ? yMax : newMax ;
561  }
562  break;
563 
564  case kZAxis:
565  pMin=zMin;
566  pMax=zMax;
567  break;
568 
569  default:
570  break;
571  }
572  pMin -= kCarTolerance ;
573  pMax += kCarTolerance ;
574 
575  return true;
576  }
577  else
578  {
579  G4int i, noEntries, noBetweenSections4 ;
580  G4bool existsAfterClip = false ;
581 
582  // Calculate rotated vertex coordinates
583 
584  G4ThreeVectorList *vertices ;
585  G4int noPolygonVertices ; // will be 4
586  vertices = CreateRotatedVertices(pTransform,noPolygonVertices) ;
587 
588  pMin = +kInfinity ;
589  pMax = -kInfinity ;
590 
591  noEntries = vertices->size() ;
592  noBetweenSections4 = noEntries - noPolygonVertices ;
593 
594  for (i=0;i<noEntries;i+=noPolygonVertices)
595  {
596  ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
597  }
598  for (i=0;i<noBetweenSections4;i+=noPolygonVertices)
599  {
600  ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
601  }
602  if (pMin!=kInfinity||pMax!=-kInfinity)
603  {
604  existsAfterClip = true ; // Add 2*tolerance to avoid precision troubles
605  pMin -= kCarTolerance ;
606  pMax += kCarTolerance ;
607  }
608  else
609  {
610  // Check for case where completely enveloping clipping volume
611  // If point inside then we are confident that the solid completely
612  // envelopes the clipping volume. Hence set min/max extents according
613  // to clipping volume extents along the specified axis.
614 
615  G4ThreeVector clipCentre(
616  (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
617  (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
618  (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 ) ;
619 
620  if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
621  {
622  existsAfterClip = true ;
623  pMin = pVoxelLimit.GetMinExtent(pAxis) ;
624  pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
625  }
626  }
627  delete vertices;
628  return existsAfterClip;
629  }
630 }
631 
633 //
634 // Return whether point inside/outside/on surface
635 
637 {
638  G4double r2, pt2, pPhi, tolRMin, tolRMax ;
639 
640  EInside in = kOutside ;
641 
642  // General precals
643  //
644  r2 = p.x()*p.x() + p.y()*p.y() ;
645  pt2 = r2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*std::sqrt(r2) ;
646 
647  if (fRmin) tolRMin = fRmin + fRminTolerance ;
648  else tolRMin = 0 ;
649 
650  tolRMax = fRmax - fRmaxTolerance;
651 
652  if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax*tolRMax )
653  {
654  if ( fDPhi == twopi || pt2 == 0 ) // on torus swept axis
655  {
656  in = kInside ;
657  }
658  else
659  {
660  // Try inner tolerant phi boundaries (=>inside)
661  // if not inside, try outer tolerant phi boundaries
662 
663  pPhi = std::atan2(p.y(),p.x()) ;
664 
665  if ( pPhi < -halfAngTolerance ) { pPhi += twopi ; } // 0<=pPhi<2pi
666  if ( fSPhi >= 0 )
667  {
668  if ( (std::fabs(pPhi) < halfAngTolerance)
669  && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
670  {
671  pPhi += twopi ; // 0 <= pPhi < 2pi
672  }
673  if ( (pPhi >= fSPhi + halfAngTolerance)
674  && (pPhi <= fSPhi + fDPhi - halfAngTolerance) )
675  {
676  in = kInside ;
677  }
678  else if ( (pPhi >= fSPhi - halfAngTolerance)
679  && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
680  {
681  in = kSurface ;
682  }
683  }
684  else // fSPhi < 0
685  {
686  if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
687  && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;}
688  else
689  {
690  in = kSurface ;
691  }
692  }
693  }
694  }
695  else // Try generous boundaries
696  {
697  tolRMin = fRmin - fRminTolerance ;
698  tolRMax = fRmax + fRmaxTolerance ;
699 
700  if (tolRMin < 0 ) { tolRMin = 0 ; }
701 
702  if ( (pt2 >= tolRMin*tolRMin) && (pt2 <= tolRMax*tolRMax) )
703  {
704  if ( (fDPhi == twopi) || (pt2 == 0) ) // Continuous in phi or on z-axis
705  {
706  in = kSurface ;
707  }
708  else // Try outer tolerant phi boundaries only
709  {
710  pPhi = std::atan2(p.y(),p.x()) ;
711 
712  if ( pPhi < -halfAngTolerance ) { pPhi += twopi ; } // 0<=pPhi<2pi
713  if ( fSPhi >= 0 )
714  {
715  if ( (std::fabs(pPhi) < halfAngTolerance)
716  && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
717  {
718  pPhi += twopi ; // 0 <= pPhi < 2pi
719  }
720  if ( (pPhi >= fSPhi - halfAngTolerance)
721  && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
722  {
723  in = kSurface;
724  }
725  }
726  else // fSPhi < 0
727  {
728  if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
729  && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;}
730  else
731  {
732  in = kSurface ;
733  }
734  }
735  }
736  }
737  }
738  return in ;
739 }
740 
742 //
743 // Return unit normal of surface closest to p
744 // - note if point on z axis, ignore phi divided sides
745 // - unsafe if point close to z axis a rmin=0 - no explicit checks
746 
748 {
749  G4int noSurfaces = 0;
750  G4double rho2, rho, pt2, pt, pPhi;
751  G4double distRMin = kInfinity;
752  G4double distSPhi = kInfinity, distEPhi = kInfinity;
753 
754  // To cope with precision loss
755  //
756  const G4double delta = std::max(10.0*kCarTolerance,
757  1.0e-8*(fRtor+fRmax));
758  const G4double dAngle = 10.0*kAngTolerance;
759 
760  G4ThreeVector nR, nPs, nPe;
761  G4ThreeVector norm, sumnorm(0.,0.,0.);
762 
763  rho2 = p.x()*p.x() + p.y()*p.y();
764  rho = std::sqrt(rho2);
765  pt2 = rho2+p.z()*p.z() +fRtor * (fRtor-2*rho);
766  pt2 = std::max(pt2, 0.0); // std::fabs(pt2);
767  pt = std::sqrt(pt2) ;
768 
769  G4double distRMax = std::fabs(pt - fRmax);
770  if(fRmin) distRMin = std::fabs(pt - fRmin);
771 
772  if( rho > delta && pt != 0.0 )
773  {
774  G4double redFactor= (rho-fRtor)/rho;
775  nR = G4ThreeVector( p.x()*redFactor, // p.x()*(1.-fRtor/rho),
776  p.y()*redFactor, // p.y()*(1.-fRtor/rho),
777  p.z() );
778  nR *= 1.0/pt;
779  }
780 
781  if ( fDPhi < twopi ) // && rho ) // old limitation against (0,0,z)
782  {
783  if ( rho )
784  {
785  pPhi = std::atan2(p.y(),p.x());
786 
787  if(pPhi < fSPhi-delta) { pPhi += twopi; }
788  else if(pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; }
789 
790  distSPhi = std::fabs( pPhi - fSPhi );
791  distEPhi = std::fabs(pPhi-fSPhi-fDPhi);
792  }
793  nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
794  nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
795  }
796  if( distRMax <= delta )
797  {
798  noSurfaces ++;
799  sumnorm += nR;
800  }
801  else if( fRmin && (distRMin <= delta) ) // Must not be on both Outer and Inner
802  {
803  noSurfaces ++;
804  sumnorm -= nR;
805  }
806 
807  // To be on one of the 'phi' surfaces,
808  // it must be within the 'tube' - with tolerance
809 
810  if( (fDPhi < twopi) && (fRmin-delta <= pt) && (pt <= (fRmax+delta)) )
811  {
812  if (distSPhi <= dAngle)
813  {
814  noSurfaces ++;
815  sumnorm += nPs;
816  }
817  if (distEPhi <= dAngle)
818  {
819  noSurfaces ++;
820  sumnorm += nPe;
821  }
822  }
823  if ( noSurfaces == 0 )
824  {
825 #ifdef G4CSGDEBUG
827  ed.precision(16);
828 
829  EInside inIt= Inside( p );
830 
831  if( inIt != kSurface )
832  {
833  ed << " ERROR> Surface Normal was called for Torus,"
834  << " with point not on surface." << G4endl;
835  }
836  else
837  {
838  ed << " ERROR> Surface Normal has not found a surface, "
839  << " despite the point being on the surface. " <<G4endl;
840  }
841 
842  if( inIt != kInside)
843  {
844  ed << " Safety (Dist To In) = " << DistanceToIn(p) << G4endl;
845  }
846  if( inIt != kOutside)
847  {
848  ed << " Safety (Dist to Out) = " << DistanceToOut(p) << G4endl;
849  }
850  ed << " Coordinates of point : " << p << G4endl;
851  ed << " Parameters of solid : " << G4endl << *this << G4endl;
852 
853  if( inIt == kSurface )
854  {
855  G4Exception("G4Torus::SurfaceNormal(p)", "GeomSolids1002",
856  JustWarning, ed,
857  "Failing to find normal, even though point is on surface!" );
858  }
859  else
860  {
861  static const char* NameInside[3]= { "Inside", "Surface", "Outside" };
862  ed << " The point is " << NameInside[inIt] << " the solid. "<< G4endl;
863  G4Exception("G4Torus::SurfaceNormal(p)", "GeomSolids1002",
864  JustWarning, ed, "Point p is not on surface !?" );
865  }
866 #endif
867  norm = ApproxSurfaceNormal(p);
868  }
869  else if ( noSurfaces == 1 ) { norm = sumnorm; }
870  else { norm = sumnorm.unit(); }
871 
872  // G4cout << "G4Torus::SurfaceNormal p= " << p << " returns norm= " << norm << G4endl;
873 
874  return norm ;
875 }
876 
878 //
879 // Algorithm for SurfaceNormal() following the original specification
880 // for points not on the surface
881 
883 {
884  ENorm side ;
885  G4ThreeVector norm;
886  G4double rho2,rho,pt2,pt,phi;
887  G4double distRMin,distRMax,distSPhi,distEPhi,distMin;
888 
889  rho2 = p.x()*p.x() + p.y()*p.y();
890  rho = std::sqrt(rho2) ;
891  pt2 = std::fabs(rho2+p.z()*p.z() +fRtor*fRtor - 2*fRtor*rho) ;
892  pt = std::sqrt(pt2) ;
893 
894 #ifdef G4CSGDEBUG
895  G4cout << " G4Torus::ApproximateSurfaceNormal called for point " << p
896  << G4endl;
897 #endif
898 
899  distRMax = std::fabs(pt - fRmax) ;
900 
901  if(fRmin) // First minimum radius
902  {
903  distRMin = std::fabs(pt - fRmin) ;
904 
905  if (distRMin < distRMax)
906  {
907  distMin = distRMin ;
908  side = kNRMin ;
909  }
910  else
911  {
912  distMin = distRMax ;
913  side = kNRMax ;
914  }
915  }
916  else
917  {
918  distMin = distRMax ;
919  side = kNRMax ;
920  }
921  if ( (fDPhi < twopi) && rho )
922  {
923  phi = std::atan2(p.y(),p.x()) ; // Protected against (0,0,z) (above rho!=0)
924 
925  if (phi < 0) { phi += twopi ; }
926 
927  if (fSPhi < 0 ) { distSPhi = std::fabs(phi-(fSPhi+twopi))*rho ; }
928  else { distSPhi = std::fabs(phi-fSPhi)*rho ; }
929 
930  distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
931 
932  if (distSPhi < distEPhi) // Find new minimum
933  {
934  if (distSPhi<distMin) side = kNSPhi ;
935  }
936  else
937  {
938  if (distEPhi < distMin) { side = kNEPhi ; }
939  }
940  }
941  switch (side)
942  {
943  case kNRMin: // Inner radius
944  norm = G4ThreeVector( -p.x()*(1-fRtor/rho)/pt,
945  -p.y()*(1-fRtor/rho)/pt,
946  -p.z()/pt ) ;
947  break ;
948  case kNRMax: // Outer radius
949  norm = G4ThreeVector( p.x()*(1-fRtor/rho)/pt,
950  p.y()*(1-fRtor/rho)/pt,
951  p.z()/pt ) ;
952  break;
953  case kNSPhi:
954  norm = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
955  break;
956  case kNEPhi:
957  norm = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
958  break;
959  default: // Should never reach this case ...
960  DumpInfo();
961  G4Exception("G4Torus::ApproxSurfaceNormal()",
962  "GeomSolids1002", JustWarning,
963  "Undefined side for valid surface normal to solid.");
964  break ;
965  }
966  return norm ;
967 }
968 
970 //
971 // Calculate distance to shape from outside, along normalised vector
972 // - return kInfinity if no intersection, or intersection distance <= tolerance
973 //
974 // - Compute the intersection with the z planes
975 // - if at valid r, phi, return
976 //
977 // -> If point is outer outer radius, compute intersection with rmax
978 // - if at valid phi,z return
979 //
980 // -> Compute intersection with inner radius, taking largest +ve root
981 // - if valid (phi), save intersction
982 //
983 // -> If phi segmented, compute intersections with phi half planes
984 // - return smallest of valid phi intersections and
985 // inner radius intersection
986 //
987 // NOTE:
988 // - Precalculations for phi trigonometry are Done `just in time'
989 // - `if valid' implies tolerant checking of intersection points
990 
992  const G4ThreeVector& v ) const
993 {
994 
995  G4double snxt=kInfinity, sphi=kInfinity; // snxt = default return value
996 
997  G4double sd[4] ;
998 
999  // Precalculated trig for phi intersections - used by r,z intersections to
1000  // check validity
1001 
1002  G4bool seg; // true if segmented
1003  G4double hDPhi; // half dphi
1004  G4double cPhi,sinCPhi=0.,cosCPhi=0.; // central phi
1005 
1006  G4double tolORMin2; // `generous' radii squared
1007  G4double tolORMax2;
1008 
1009  G4double Dist,xi,yi,zi,rhoi2,it2; // Intersection point variables
1010 
1011  G4double Comp;
1012  G4double cosSPhi,sinSPhi; // Trig for phi start intersect
1013  G4double ePhi,cosEPhi,sinEPhi; // for phi end intersect
1014 
1015  // Set phi divided flag and precalcs
1016  //
1017  if ( fDPhi < twopi )
1018  {
1019  seg = true ;
1020  hDPhi = 0.5*fDPhi ; // half delta phi
1021  cPhi = fSPhi + hDPhi ;
1022  sinCPhi = std::sin(cPhi) ;
1023  cosCPhi = std::cos(cPhi) ;
1024  }
1025  else
1026  {
1027  seg = false ;
1028  }
1029 
1030  if (fRmin > fRminTolerance) // Calculate tolerant rmin and rmax
1031  {
1032  tolORMin2 = (fRmin - fRminTolerance)*(fRmin - fRminTolerance) ;
1033  }
1034  else
1035  {
1036  tolORMin2 = 0 ;
1037  }
1038  tolORMax2 = (fRmax + fRmaxTolerance)*(fRmax + fRmaxTolerance) ;
1039 
1040  // Intersection with Rmax (possible return) and Rmin (must also check phi)
1041 
1042  G4double Rtor2 = fRtor*fRtor ;
1043 
1044  snxt = SolveNumericJT(p,v,fRmax,true);
1045 
1046  if (fRmin) // Possible Rmin intersection
1047  {
1048  sd[0] = SolveNumericJT(p,v,fRmin,true);
1049  if ( sd[0] < snxt ) { snxt = sd[0] ; }
1050  }
1051 
1052  //
1053  // Phi segment intersection
1054  //
1055  // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1056  //
1057  // o NOTE: Large duplication of code between sphi & ephi checks
1058  // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1059  // intersection check <=0 -> >=0
1060  // -> use some form of loop Construct ?
1061 
1062  if (seg)
1063  {
1064  sinSPhi = std::sin(fSPhi) ; // First phi surface ('S'tarting phi)
1065  cosSPhi = std::cos(fSPhi) ;
1066  Comp = v.x()*sinSPhi - v.y()*cosSPhi ; // Component in outwards
1067  // normal direction
1068  if (Comp < 0 )
1069  {
1070  Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1071 
1072  if (Dist < halfCarTolerance)
1073  {
1074  sphi = Dist/Comp ;
1075  if (sphi < snxt)
1076  {
1077  if ( sphi < 0 ) { sphi = 0 ; }
1078 
1079  xi = p.x() + sphi*v.x() ;
1080  yi = p.y() + sphi*v.y() ;
1081  zi = p.z() + sphi*v.z() ;
1082  rhoi2 = xi*xi + yi*yi ;
1083  it2 = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ;
1084 
1085  if ( it2 >= tolORMin2 && it2 <= tolORMax2 )
1086  {
1087  // r intersection is good - check intersecting
1088  // with correct half-plane
1089  //
1090  if ((yi*cosCPhi-xi*sinCPhi)<=0) { snxt=sphi; }
1091  }
1092  }
1093  }
1094  }
1095  ePhi=fSPhi+fDPhi; // Second phi surface ('E'nding phi)
1096  sinEPhi=std::sin(ePhi);
1097  cosEPhi=std::cos(ePhi);
1098  Comp=-(v.x()*sinEPhi-v.y()*cosEPhi);
1099 
1100  if ( Comp < 0 ) // Component in outwards normal dirn
1101  {
1102  Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1103 
1104  if (Dist < halfCarTolerance )
1105  {
1106  sphi = Dist/Comp ;
1107 
1108  if (sphi < snxt )
1109  {
1110  if (sphi < 0 ) { sphi = 0 ; }
1111 
1112  xi = p.x() + sphi*v.x() ;
1113  yi = p.y() + sphi*v.y() ;
1114  zi = p.z() + sphi*v.z() ;
1115  rhoi2 = xi*xi + yi*yi ;
1116  it2 = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ;
1117 
1118  if (it2 >= tolORMin2 && it2 <= tolORMax2)
1119  {
1120  // z and r intersections good - check intersecting
1121  // with correct half-plane
1122  //
1123  if ((yi*cosCPhi-xi*sinCPhi)>=0) { snxt=sphi; }
1124  }
1125  }
1126  }
1127  }
1128  }
1129  if(snxt < halfCarTolerance) { snxt = 0.0 ; }
1130 
1131  return snxt ;
1132 }
1133 
1135 //
1136 // Calculate distance (<= actual) to closest surface of shape from outside
1137 // - Calculate distance to z, radial planes
1138 // - Only to phi planes if outside phi extent
1139 // - Return 0 if point inside
1140 
1142 {
1143  G4double safe=0.0, safe1, safe2 ;
1144  G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
1145  G4double rho2, rho, pt2, pt ;
1146 
1147  rho2 = p.x()*p.x() + p.y()*p.y() ;
1148  rho = std::sqrt(rho2) ;
1149  pt2 = std::fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ;
1150  pt = std::sqrt(pt2) ;
1151 
1152  safe1 = fRmin - pt ;
1153  safe2 = pt - fRmax ;
1154 
1155  if (safe1 > safe2) { safe = safe1; }
1156  else { safe = safe2; }
1157 
1158  if ( fDPhi < twopi && rho )
1159  {
1160  phiC = fSPhi + fDPhi*0.5 ;
1161  cosPhiC = std::cos(phiC) ;
1162  sinPhiC = std::sin(phiC) ;
1163  cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
1164 
1165  if (cosPsi < std::cos(fDPhi*0.5) ) // Psi=angle from central phi to point
1166  { // Point lies outside phi range
1167  if ((p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
1168  {
1169  safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
1170  }
1171  else
1172  {
1173  ePhi = fSPhi + fDPhi ;
1174  safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1175  }
1176  if (safePhi > safe) { safe = safePhi ; }
1177  }
1178  }
1179  if (safe < 0 ) { safe = 0 ; }
1180  return safe;
1181 }
1182 
1184 //
1185 // Calculate distance to surface of shape from `inside', allowing for tolerance
1186 // - Only Calc rmax intersection if no valid rmin intersection
1187 //
1188 
1190  const G4ThreeVector& v,
1191  const G4bool calcNorm,
1192  G4bool *validNorm,
1193  G4ThreeVector *n ) const
1194 {
1195  ESide side = kNull, sidephi = kNull ;
1196  G4double snxt = kInfinity, sphi, sd[4] ;
1197 
1198  // Vars for phi intersection
1199  //
1200  G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi;
1201  G4double cPhi, sinCPhi, cosCPhi ;
1202  G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, zi, vphi ;
1203 
1204  // Radial Intersections Defenitions & General Precals
1205 
1207 
1208 #if 1
1209 
1210  // This is the version with the calculation of CalcNorm = true
1211  // To be done: Check the precision of this calculation.
1212  // If you want return always validNorm = false, then take the version below
1213 
1214  G4double rho2 = p.x()*p.x()+p.y()*p.y();
1215  G4double rho = std::sqrt(rho2) ;
1216 
1217  G4double pt2 = rho2 + p.z()*p.z() + fRtor * (fRtor - 2.0*rho);
1218  // Regroup for slightly better FP accuracy
1219 
1220  if( pt2 < 0.0)
1221  {
1222  pt2= std::fabs( pt2 );
1223  }
1224 
1225  G4double pt = std::sqrt(pt2) ;
1226 
1227  G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
1228 
1229  G4double tolRMax = fRmax - fRmaxTolerance ;
1230 
1231  G4double vDotNmax = pDotV - fRtor*(v.x()*p.x() + v.y()*p.y())/rho ;
1232  G4double pDotxyNmax = (1 - fRtor/rho) ;
1233 
1234  if( (pt2 > tolRMax*tolRMax) && (vDotNmax >= 0) )
1235  {
1236  // On tolerant boundary & heading outwards (or perpendicular to) outer
1237  // radial surface -> leaving immediately with *n for really convex part
1238  // only
1239 
1240  if ( calcNorm && (pDotxyNmax >= -2.*fRmaxTolerance) )
1241  {
1242  *n = G4ThreeVector( p.x()*(1 - fRtor/rho)/pt,
1243  p.y()*(1 - fRtor/rho)/pt,
1244  p.z()/pt ) ;
1245  *validNorm = true ;
1246  }
1247 
1248  return snxt = 0 ; // Leaving by Rmax immediately
1249  }
1250 
1251  snxt = SolveNumericJT(p,v,fRmax,false);
1252  side = kRMax ;
1253 
1254  // rmin
1255 
1256  if ( fRmin )
1257  {
1258  G4double tolRMin = fRmin + fRminTolerance ;
1259 
1260  if ( (pt2 < tolRMin*tolRMin) && (vDotNmax < 0) )
1261  {
1262  if (calcNorm) { *validNorm = false ; } // Concave surface of the torus
1263  return snxt = 0 ; // Leaving by Rmin immediately
1264  }
1265 
1266  sd[0] = SolveNumericJT(p,v,fRmin,false);
1267  if ( sd[0] < snxt )
1268  {
1269  snxt = sd[0] ;
1270  side = kRMin ;
1271  }
1272  }
1273 
1274 #else
1275 
1276  // this is the "conservative" version which return always validnorm = false
1277  // NOTE: using this version the unit test testG4Torus will break
1278 
1279  snxt = SolveNumericJT(p,v,fRmax,false);
1280  side = kRMax ;
1281 
1282  if ( fRmin )
1283  {
1284  sd[0] = SolveNumericJT(p,v,fRmin,false);
1285  if ( sd[0] < snxt )
1286  {
1287  snxt = sd[0] ;
1288  side = kRMin ;
1289  }
1290  }
1291 
1292  if ( calcNorm && (snxt == 0.0) )
1293  {
1294  *validNorm = false ; // Leaving solid, but possible re-intersection
1295  return snxt ;
1296  }
1297 
1298 #endif
1299 
1300  if (fDPhi < twopi) // Phi Intersections
1301  {
1302  sinSPhi = std::sin(fSPhi) ;
1303  cosSPhi = std::cos(fSPhi) ;
1304  ePhi = fSPhi + fDPhi ;
1305  sinEPhi = std::sin(ePhi) ;
1306  cosEPhi = std::cos(ePhi) ;
1307  cPhi = fSPhi + fDPhi*0.5 ;
1308  sinCPhi = std::sin(cPhi) ;
1309  cosCPhi = std::cos(cPhi) ;
1310 
1311  // angle calculation with correction
1312  // of difference in domain of atan2 and Sphi
1313  //
1314  vphi = std::atan2(v.y(),v.x()) ;
1315 
1316  if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1317  else if ( vphi > ePhi + halfAngTolerance ) { vphi -= twopi; }
1318 
1319  if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1320  {
1321  pDistS = p.x()*sinSPhi - p.y()*cosSPhi ; // pDist -ve when inside
1322  pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1323 
1324  // Comp -ve when in direction of outwards normal
1325  //
1326  compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1327  compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1328  sidephi = kNull ;
1329 
1330  if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1331  && (pDistE <= halfCarTolerance) ) )
1332  || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
1333  && (pDistE > halfCarTolerance) ) ) )
1334  {
1335  // Inside both phi *full* planes
1336 
1337  if ( compS < 0 )
1338  {
1339  sphi = pDistS/compS ;
1340 
1341  if (sphi >= -halfCarTolerance)
1342  {
1343  xi = p.x() + sphi*v.x() ;
1344  yi = p.y() + sphi*v.y() ;
1345 
1346  // Check intersecting with correct half-plane
1347  // (if not -> no intersect)
1348  //
1349  if ( (std::fabs(xi)<=kCarTolerance)
1350  && (std::fabs(yi)<=kCarTolerance) )
1351  {
1352  sidephi = kSPhi;
1353  if ( ((fSPhi-halfAngTolerance)<=vphi)
1354  && ((ePhi+halfAngTolerance)>=vphi) )
1355  {
1356  sphi = kInfinity;
1357  }
1358  }
1359  else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1360  {
1361  sphi = kInfinity ;
1362  }
1363  else
1364  {
1365  sidephi = kSPhi ;
1366  }
1367  }
1368  else
1369  {
1370  sphi = kInfinity ;
1371  }
1372  }
1373  else
1374  {
1375  sphi = kInfinity ;
1376  }
1377 
1378  if ( compE < 0 )
1379  {
1380  sphi2 = pDistE/compE ;
1381 
1382  // Only check further if < starting phi intersection
1383  //
1384  if ( (sphi2 > -kCarTolerance) && (sphi2 < sphi) )
1385  {
1386  xi = p.x() + sphi2*v.x() ;
1387  yi = p.y() + sphi2*v.y() ;
1388 
1389  if ( (std::fabs(xi)<=kCarTolerance)
1390  && (std::fabs(yi)<=kCarTolerance) )
1391  {
1392  // Leaving via ending phi
1393  //
1394  if( !( (fSPhi-halfAngTolerance <= vphi)
1395  && (ePhi+halfAngTolerance >= vphi) ) )
1396  {
1397  sidephi = kEPhi ;
1398  sphi = sphi2;
1399  }
1400  }
1401  else // Check intersecting with correct half-plane
1402  {
1403  if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1404  {
1405  // Leaving via ending phi
1406  //
1407  sidephi = kEPhi ;
1408  sphi = sphi2;
1409 
1410  }
1411  }
1412  }
1413  }
1414  }
1415  else
1416  {
1417  sphi = kInfinity ;
1418  }
1419  }
1420  else
1421  {
1422  // On z axis + travel not || to z axis -> if phi of vector direction
1423  // within phi of shape, Step limited by rmax, else Step =0
1424 
1425  vphi = std::atan2(v.y(),v.x());
1426 
1427  if ( ( fSPhi-halfAngTolerance <= vphi ) &&
1428  ( vphi <= ( ePhi+halfAngTolerance ) ) )
1429  {
1430  sphi = kInfinity;
1431  }
1432  else
1433  {
1434  sidephi = kSPhi ; // arbitrary
1435  sphi=0;
1436  }
1437  }
1438 
1439  // Order intersections
1440 
1441  if (sphi<snxt)
1442  {
1443  snxt=sphi;
1444  side=sidephi;
1445  }
1446  }
1447 
1448  G4double rhoi2,rhoi,it2,it,iDotxyNmax ;
1449  // Note: by numerical computation we know where the ray hits the torus
1450  // So I propose to return the side where the ray hits
1451 
1452  if (calcNorm)
1453  {
1454  switch(side)
1455  {
1456  case kRMax: // n is unit vector
1457  xi = p.x() + snxt*v.x() ;
1458  yi =p.y() + snxt*v.y() ;
1459  zi = p.z() + snxt*v.z() ;
1460  rhoi2 = xi*xi + yi*yi ;
1461  rhoi = std::sqrt(rhoi2) ;
1462  it2 = std::fabs(rhoi2 + zi*zi + fRtor*fRtor - 2*fRtor*rhoi) ;
1463  it = std::sqrt(it2) ;
1464  iDotxyNmax = (1-fRtor/rhoi) ;
1465  if(iDotxyNmax >= -2.*fRmaxTolerance) // really convex part of Rmax
1466  {
1467  *n = G4ThreeVector( xi*(1-fRtor/rhoi)/it,
1468  yi*(1-fRtor/rhoi)/it,
1469  zi/it ) ;
1470  *validNorm = true ;
1471  }
1472  else
1473  {
1474  *validNorm = false ; // concave-convex part of Rmax
1475  }
1476  break ;
1477 
1478  case kRMin:
1479  *validNorm = false ; // Rmin is concave or concave-convex
1480  break;
1481 
1482  case kSPhi:
1483  if (fDPhi <= pi )
1484  {
1485  *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
1486  *validNorm=true;
1487  }
1488  else
1489  {
1490  *validNorm = false ;
1491  }
1492  break ;
1493 
1494  case kEPhi:
1495  if (fDPhi <= pi)
1496  {
1497  *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
1498  *validNorm=true;
1499  }
1500  else
1501  {
1502  *validNorm = false ;
1503  }
1504  break;
1505 
1506  default:
1507 
1508  // It seems we go here from time to time ...
1509 
1510  G4cout << G4endl;
1511  DumpInfo();
1512  std::ostringstream message;
1513  G4int oldprc = message.precision(16);
1514  message << "Undefined side for valid surface normal to solid."
1515  << G4endl
1516  << "Position:" << G4endl << G4endl
1517  << "p.x() = " << p.x()/mm << " mm" << G4endl
1518  << "p.y() = " << p.y()/mm << " mm" << G4endl
1519  << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
1520  << "Direction:" << G4endl << G4endl
1521  << "v.x() = " << v.x() << G4endl
1522  << "v.y() = " << v.y() << G4endl
1523  << "v.z() = " << v.z() << G4endl << G4endl
1524  << "Proposed distance :" << G4endl << G4endl
1525  << "snxt = " << snxt/mm << " mm" << G4endl;
1526  message.precision(oldprc);
1527  G4Exception("G4Torus::DistanceToOut(p,v,..)",
1528  "GeomSolids1002",JustWarning, message);
1529  break;
1530  }
1531  }
1532  if ( snxt<halfCarTolerance ) { snxt=0 ; }
1533 
1534  return snxt;
1535 }
1536 
1538 //
1539 // Calculate distance (<=actual) to closest surface of shape from inside
1540 
1542 {
1543  G4double safe=0.0,safeR1,safeR2;
1544  G4double rho2,rho,pt2,pt ;
1545  G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi;
1546  rho2 = p.x()*p.x() + p.y()*p.y() ;
1547  rho = std::sqrt(rho2) ;
1548  pt2 = std::fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ;
1549  pt = std::sqrt(pt2) ;
1550 
1551 #ifdef G4CSGDEBUG
1552  if( Inside(p) == kOutside )
1553  {
1554  G4int oldprc = G4cout.precision(16) ;
1555  G4cout << G4endl ;
1556  DumpInfo();
1557  G4cout << "Position:" << G4endl << G4endl ;
1558  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1559  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1560  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1561  G4cout.precision(oldprc);
1562  G4Exception("G4Torus::DistanceToOut(p)", "GeomSolids1002",
1563  JustWarning, "Point p is outside !?" );
1564  }
1565 #endif
1566 
1567  if (fRmin)
1568  {
1569  safeR1 = pt - fRmin ;
1570  safeR2 = fRmax - pt ;
1571 
1572  if (safeR1 < safeR2) { safe = safeR1 ; }
1573  else { safe = safeR2 ; }
1574  }
1575  else
1576  {
1577  safe = fRmax - pt ;
1578  }
1579 
1580  // Check if phi divided, Calc distances closest phi plane
1581  //
1582  if (fDPhi<twopi) // Above/below central phi of Torus?
1583  {
1584  phiC = fSPhi + fDPhi*0.5 ;
1585  cosPhiC = std::cos(phiC) ;
1586  sinPhiC = std::sin(phiC) ;
1587 
1588  if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
1589  {
1590  safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
1591  }
1592  else
1593  {
1594  ePhi = fSPhi + fDPhi ;
1595  safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1596  }
1597  if (safePhi < safe) { safe = safePhi ; }
1598  }
1599  if (safe < 0) { safe = 0 ; }
1600  return safe ;
1601 }
1602 
1604 //
1605 // Create a List containing the transformed vertices
1606 // Ordering [0-3] -fRtor cross section
1607 // [4-7] +fRtor cross section such that [0] is below [4],
1608 // [1] below [5] etc.
1609 // Note:
1610 // Caller has deletion resposibility
1611 // Potential improvement: For last slice, use actual ending angle
1612 // to avoid rounding error problems.
1613 
1616  G4int& noPolygonVertices ) const
1617 {
1618  G4ThreeVectorList *vertices;
1619  G4ThreeVector vertex0,vertex1,vertex2,vertex3;
1620  G4double meshAngle,meshRMax,crossAngle,cosCrossAngle,sinCrossAngle,sAngle;
1621  G4double rMaxX,rMaxY,rMinX,rMinY;
1622  G4int crossSection,noCrossSections;
1623 
1624  // Compute no of cross-sections necessary to mesh tube
1625  //
1626  noCrossSections = G4int (fDPhi/kMeshAngleDefault) + 1 ;
1627 
1628  if (noCrossSections < kMinMeshSections)
1629  {
1630  noCrossSections = kMinMeshSections ;
1631  }
1632  else if (noCrossSections>kMaxMeshSections)
1633  {
1634  noCrossSections=kMaxMeshSections;
1635  }
1636  meshAngle = fDPhi/(noCrossSections - 1) ;
1637  meshRMax = (fRtor + fRmax)/std::cos(meshAngle*0.5) ;
1638 
1639  // If complete in phi, set start angle such that mesh will be at fRmax
1640  // on the x axis. Will give better extent calculations when not rotated
1641 
1642  if ( (fDPhi == twopi) && (fSPhi == 0) )
1643  {
1644  sAngle = -meshAngle*0.5 ;
1645  }
1646  else
1647  {
1648  sAngle = fSPhi ;
1649  }
1650  vertices = new G4ThreeVectorList();
1651 
1652  if (vertices)
1653  {
1654  vertices->reserve(noCrossSections*4) ;
1655  for (crossSection=0;crossSection<noCrossSections;crossSection++)
1656  {
1657  // Compute coordinates of cross section at section crossSection
1658 
1659  crossAngle=sAngle+crossSection*meshAngle;
1660  cosCrossAngle=std::cos(crossAngle);
1661  sinCrossAngle=std::sin(crossAngle);
1662 
1663  rMaxX=meshRMax*cosCrossAngle;
1664  rMaxY=meshRMax*sinCrossAngle;
1665  rMinX=(fRtor-fRmax)*cosCrossAngle;
1666  rMinY=(fRtor-fRmax)*sinCrossAngle;
1667  vertex0=G4ThreeVector(rMinX,rMinY,-fRmax);
1668  vertex1=G4ThreeVector(rMaxX,rMaxY,-fRmax);
1669  vertex2=G4ThreeVector(rMaxX,rMaxY,+fRmax);
1670  vertex3=G4ThreeVector(rMinX,rMinY,+fRmax);
1671 
1672  vertices->push_back(pTransform.TransformPoint(vertex0));
1673  vertices->push_back(pTransform.TransformPoint(vertex1));
1674  vertices->push_back(pTransform.TransformPoint(vertex2));
1675  vertices->push_back(pTransform.TransformPoint(vertex3));
1676  }
1677  noPolygonVertices = 4 ;
1678  }
1679  else
1680  {
1681  DumpInfo();
1682  G4Exception("G4Torus::CreateRotatedVertices()",
1683  "GeomSolids0003", FatalException,
1684  "Error in allocation of vertices. Out of memory !");
1685  }
1686  return vertices;
1687 }
1688 
1690 //
1691 // Stream object contents to an output stream
1692 
1694 {
1695  return G4String("G4Torus");
1696 }
1697 
1699 //
1700 // Make a clone of the object
1701 //
1703 {
1704  return new G4Torus(*this);
1705 }
1706 
1708 //
1709 // Stream object contents to an output stream
1710 
1711 std::ostream& G4Torus::StreamInfo( std::ostream& os ) const
1712 {
1713  G4int oldprc = os.precision(16);
1714  os << "-----------------------------------------------------------\n"
1715  << " *** Dump for solid - " << GetName() << " ***\n"
1716  << " ===================================================\n"
1717  << " Solid type: G4Torus\n"
1718  << " Parameters: \n"
1719  << " inner radius: " << fRmin/mm << " mm \n"
1720  << " outer radius: " << fRmax/mm << " mm \n"
1721  << " swept radius: " << fRtor/mm << " mm \n"
1722  << " starting phi: " << fSPhi/degree << " degrees \n"
1723  << " delta phi : " << fDPhi/degree << " degrees \n"
1724  << "-----------------------------------------------------------\n";
1725  os.precision(oldprc);
1726 
1727  return os;
1728 }
1729 
1731 //
1732 // GetPointOnSurface
1733 
1735 {
1736  G4double cosu, sinu,cosv, sinv, aOut, aIn, aSide, chose, phi, theta, rRand;
1737 
1739  theta = RandFlat::shoot(0.,twopi);
1740 
1741  cosu = std::cos(phi); sinu = std::sin(phi);
1742  cosv = std::cos(theta); sinv = std::sin(theta);
1743 
1744  // compute the areas
1745 
1746  aOut = (fDPhi)*twopi*fRtor*fRmax;
1747  aIn = (fDPhi)*twopi*fRtor*fRmin;
1748  aSide = pi*(fRmax*fRmax-fRmin*fRmin);
1749 
1750  if ((fSPhi == 0) && (fDPhi == twopi)){ aSide = 0; }
1751  chose = RandFlat::shoot(0.,aOut + aIn + 2.*aSide);
1752 
1753  if(chose < aOut)
1754  {
1755  return G4ThreeVector ((fRtor+fRmax*cosv)*cosu,
1756  (fRtor+fRmax*cosv)*sinu, fRmax*sinv);
1757  }
1758  else if( (chose >= aOut) && (chose < aOut + aIn) )
1759  {
1760  return G4ThreeVector ((fRtor+fRmin*cosv)*cosu,
1761  (fRtor+fRmin*cosv)*sinu, fRmin*sinv);
1762  }
1763  else if( (chose >= aOut + aIn) && (chose < aOut + aIn + aSide) )
1764  {
1765  rRand = GetRadiusInRing(fRmin,fRmax);
1766  return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi),
1767  (fRtor+rRand*cosv)*std::sin(fSPhi), rRand*sinv);
1768  }
1769  else
1770  {
1771  rRand = GetRadiusInRing(fRmin,fRmax);
1772  return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi+fDPhi),
1773  (fRtor+rRand*cosv)*std::sin(fSPhi+fDPhi),
1774  rRand*sinv);
1775  }
1776 }
1777 
1779 //
1780 // Visualisation Functions
1781 
1783 {
1784  scene.AddSolid (*this);
1785 }
1786 
1788 {
1789  return new G4PolyhedronTorus (fRmin, fRmax, fRtor, fSPhi, fDPhi);
1790 }
1791 
1792 #endif // !defined(G4GEOM_USE_TORUS) || !defined(G4GEOM_USE_SYS_USOLIDS)
G4String GetName() const
EInside Inside(const G4ThreeVector &p) const
Definition: G4Torus.cc:636
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)
G4double fRminTolerance
Definition: G4Torus.hh:211
G4double fRmin
Definition: G4Torus.hh:203
G4Polyhedron * CreatePolyhedron() const
Definition: G4Torus.cc:1787
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double GetMinYExtent() const
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
CLHEP::Hep3Vector G4ThreeVector
G4bool fRebuildPolyhedron
Definition: G4CSGSolid.hh:80
G4Torus & operator=(const G4Torus &rhs)
Definition: G4Torus.cc:216
G4AffineTransform Inverse() const
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Torus.cc:1711
G4bool IsYLimited() const
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:111
G4double fRmax
Definition: G4Torus.hh:203
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Torus.cc:991
G4bool IsRotated() const
G4double SolveNumericJT(const G4ThreeVector &p, const G4ThreeVector &v, G4double r, G4bool IsDistanceToIn) const
Definition: G4Torus.cc:297
static const G4double e2
G4double fCubicVolume
Definition: G4CSGSolid.hh:78
G4ThreeVector NetTranslation() const
~G4Torus()
Definition: G4Torus.cc:195
G4bool IsXLimited() const
G4double a
Definition: TRTMaterials.hh:39
const G4int kMinMeshSections
Definition: meshdefs.hh:45
virtual void AddSolid(const G4Box &)=0
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Torus.cc:416
int G4int
Definition: G4Types.hh:78
void DumpInfo() const
G4double GetMaxXExtent() const
G4double GetMinZExtent() const
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform, G4int &noPolygonVertices) const
Definition: G4Torus.cc:1615
G4double fSurfaceArea
Definition: G4CSGSolid.hh:79
G4GLOB_DLL std::ostream G4cout
G4double kRadTolerance
Definition: G4Torus.hh:211
G4double GetRmin() const
G4double GetRadialTolerance() const
bool G4bool
Definition: G4Types.hh:79
static const double twopi
Definition: G4SIunits.hh:75
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
G4Torus(const G4String &pName, G4double pRmin, G4double pRmax, G4double pRtor, G4double pSPhi, G4double pDPhi)
Definition: G4Torus.cc:82
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Torus.cc:747
const G4int n
G4VSolid * Clone() const
Definition: G4Torus.cc:1702
G4ThreeVector GetPointOnSurface() const
Definition: G4Torus.cc:1734
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4double fSPhi
Definition: G4Torus.hh:203
G4double GetMinXExtent() const
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Torus.cc:243
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Torus.cc:1189
G4double halfCarTolerance
Definition: G4Torus.hh:214
G4double kAngTolerance
Definition: G4Torus.hh:211
G4double GetMaxZExtent() const
static const double pi
Definition: G4SIunits.hh:74
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EInside
Definition: geomdefs.hh:58
G4double fDPhi
Definition: G4Torus.hh:203
EAxis
Definition: geomdefs.hh:54
G4double halfAngTolerance
Definition: G4Torus.hh:214
static const double degree
Definition: G4SIunits.hh:143
G4GeometryType GetEntityType() const
Definition: G4Torus.cc:1693
#define G4endl
Definition: G4ios.hh:61
void TorusRootsJT(const G4ThreeVector &p, const G4ThreeVector &v, G4double r, std::vector< G4double > &roots) const
Definition: G4Torus.cc:257
G4double GetMaxYExtent() const
G4double fRtor
Definition: G4Torus.hh:203
G4double kCarTolerance
Definition: G4VSolid.hh:304
void SetAllParameters(G4double pRmin, G4double pRmax, G4double pRtor, G4double pSPhi, G4double pDPhi)
Definition: G4Torus.cc:98
double G4double
Definition: G4Types.hh:76
G4double fRmaxTolerance
Definition: G4Torus.hh:211
G4double GetMaxExtent(const EAxis pAxis) const
static const G4double e3
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Torus.cc:1782
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:91
G4bool IsZLimited() const
static const double mm
Definition: G4SIunits.hh:114
G4int FindRoots(G4double *op, G4int degree, G4double *zeror, G4double *zeroi)
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 GetMinExtent(const EAxis pAxis) const
const G4double kMeshAngleDefault
Definition: meshdefs.hh:42
static G4GeometryTolerance * GetInstance()
const G4int kMaxMeshSections
Definition: meshdefs.hh:46
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4Torus.cc:882