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