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