Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 //
27 // $Id: G4Torus.cc 104316 2017-05-24 13:04:23Z 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 // 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 
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  //
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;
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 
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 
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::BoundingLimits()", "GeomMgt0001",
453  JustWarning, message);
454  DumpInfo();
455  }
456 }
457 
459 //
460 // Calculate extent under transform and specified limit
461 
463  const G4VoxelLimits& pVoxelLimit,
464  const G4AffineTransform& pTransform,
465  G4double& pMin, G4double& pMax) const
466 {
467  G4ThreeVector bmin, bmax;
468  G4bool exist;
469 
470  // Get bounding box
471  BoundingLimits(bmin,bmax);
472 
473  // Check bounding box
474  G4BoundingEnvelope bbox(bmin,bmax);
475 #ifdef G4BBOX_EXTENT
476  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
477 #endif
478  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
479  {
480  return exist = (pMin < pMax) ? true : false;
481  }
482 
483  // Get parameters of the solid
484  G4double rmin = GetRmin();
485  G4double rmax = GetRmax();
486  G4double rtor = GetRtor();
487  G4double dphi = GetDPhi();
488  G4double sinStart = GetSinStartPhi();
489  G4double cosStart = GetCosStartPhi();
490  G4double sinEnd = GetSinEndPhi();
491  G4double cosEnd = GetCosEndPhi();
492  G4double rint = rtor - rmax;
493  G4double rext = rtor + rmax;
494 
495  // Find bounding envelope and calculate extent
496  //
497  static const G4int NPHI = 24; // number of steps for whole torus
498  static const G4int NDISK = 16; // number of steps for disk
499  static const G4double sinHalfDisk = std::sin(pi/NDISK);
500  static const G4double cosHalfDisk = std::cos(pi/NDISK);
501  static const G4double sinStepDisk = 2.*sinHalfDisk*cosHalfDisk;
502  static const G4double cosStepDisk = 1. - 2.*sinHalfDisk*sinHalfDisk;
503 
504  G4double astep = (360/NPHI)*deg; // max angle for one slice in phi
505  G4int kphi = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
506  G4double ang = dphi/kphi;
507 
508  G4double sinHalf = std::sin(0.5*ang);
509  G4double cosHalf = std::cos(0.5*ang);
510  G4double sinStep = 2.*sinHalf*cosHalf;
511  G4double cosStep = 1. - 2.*sinHalf*sinHalf;
512 
513  // define vectors for bounding envelope
514  G4ThreeVectorList pols[NDISK+1];
515  for (G4int k=0; k<NDISK+1; ++k) pols[k].resize(4);
516 
517  std::vector<const G4ThreeVectorList *> polygons;
518  polygons.resize(NDISK+1);
519  for (G4int k=0; k<NDISK+1; ++k) polygons[k] = &pols[k];
520 
521  // set internal and external reference circles
522  G4TwoVector rzmin[NDISK];
523  G4TwoVector rzmax[NDISK];
524 
525  if ((rtor-rmin*sinHalfDisk)/cosHalf > (rtor+rmin*sinHalfDisk)) rmin = 0;
526  rmax /= cosHalfDisk;
527  G4double sinCurDisk = sinHalfDisk;
528  G4double cosCurDisk = cosHalfDisk;
529  for (G4int k=0; k<NDISK; ++k)
530  {
531  G4double rmincur = rtor + rmin*cosCurDisk;
532  if (cosCurDisk < 0 && rmin > 0) rmincur /= cosHalf;
533  rzmin[k].set(rmincur,rmin*sinCurDisk);
534 
535  G4double rmaxcur = rtor + rmax*cosCurDisk;
536  if (cosCurDisk > 0) rmaxcur /= cosHalf;
537  rzmax[k].set(rmaxcur,rmax*sinCurDisk);
538 
539  G4double sinTmpDisk = sinCurDisk;
540  sinCurDisk = sinCurDisk*cosStepDisk + cosCurDisk*sinStepDisk;
541  cosCurDisk = cosCurDisk*cosStepDisk - sinTmpDisk*sinStepDisk;
542  }
543 
544  // Loop along slices in Phi. The extent is calculated as cumulative
545  // extent of the slices
546  pMin = kInfinity;
547  pMax = -kInfinity;
548  G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
549  G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
550  G4double sinCur1 = 0, cosCur1 = 0, sinCur2 = 0, cosCur2 = 0;
551  for (G4int i=0; i<kphi+1; ++i)
552  {
553  if (i == 0)
554  {
555  sinCur1 = sinStart;
556  cosCur1 = cosStart;
557  sinCur2 = sinCur1*cosHalf + cosCur1*sinHalf;
558  cosCur2 = cosCur1*cosHalf - sinCur1*sinHalf;
559  }
560  else
561  {
562  sinCur1 = sinCur2;
563  cosCur1 = cosCur2;
564  sinCur2 = (i == kphi) ? sinEnd : sinCur1*cosStep + cosCur1*sinStep;
565  cosCur2 = (i == kphi) ? cosEnd : cosCur1*cosStep - sinCur1*sinStep;
566  }
567  for (G4int k=0; k<NDISK; ++k)
568  {
569  G4double r1 = rzmin[k].x(), r2 = rzmax[k].x();
570  G4double z1 = rzmin[k].y(), z2 = rzmax[k].y();
571  pols[k][0].set(r1*cosCur1,r1*sinCur1,z1);
572  pols[k][1].set(r2*cosCur1,r2*sinCur1,z2);
573  pols[k][2].set(r2*cosCur2,r2*sinCur2,z2);
574  pols[k][3].set(r1*cosCur2,r1*sinCur2,z1);
575  }
576  pols[NDISK] = pols[0];
577 
578  // get bounding box of current slice
579  G4TwoVector vmin,vmax;
581  DiskExtent(rint,rext,sinCur1,cosCur1,sinCur2,cosCur2,vmin,vmax);
582  bmin.setX(vmin.x()); bmin.setY(vmin.y());
583  bmax.setX(vmax.x()); bmax.setY(vmax.y());
584 
585  // set bounding envelope for current slice and adjust extent
586  G4double emin,emax;
587  G4BoundingEnvelope benv(bmin,bmax,polygons);
588  if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
589  if (emin < pMin) pMin = emin;
590  if (emax > pMax) pMax = emax;
591  if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
592  }
593  return (pMin < pMax);
594 }
595 
597 //
598 // Return whether point inside/outside/on surface
599 
601 {
602  G4double r, pt2, pPhi, tolRMin, tolRMax ;
603 
604  EInside in = kOutside ;
605 
606  // General precals
607  //
608  r = std::hypot(p.x(),p.y());
609  pt2 = p.z()*p.z() + (r-fRtor)*(r-fRtor);
610 
611  if (fRmin) tolRMin = fRmin + fRminTolerance ;
612  else tolRMin = 0 ;
613 
614  tolRMax = fRmax - fRmaxTolerance;
615 
616  if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax*tolRMax )
617  {
618  if ( fDPhi == twopi || pt2 == 0 ) // on torus swept axis
619  {
620  in = kInside ;
621  }
622  else
623  {
624  // Try inner tolerant phi boundaries (=>inside)
625  // if not inside, try outer tolerant phi boundaries
626 
627  pPhi = std::atan2(p.y(),p.x()) ;
628 
629  if ( pPhi < -halfAngTolerance ) { pPhi += twopi ; } // 0<=pPhi<2pi
630  if ( fSPhi >= 0 )
631  {
632  if ( (std::fabs(pPhi) < halfAngTolerance)
633  && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
634  {
635  pPhi += twopi ; // 0 <= pPhi < 2pi
636  }
637  if ( (pPhi >= fSPhi + halfAngTolerance)
638  && (pPhi <= fSPhi + fDPhi - halfAngTolerance) )
639  {
640  in = kInside ;
641  }
642  else if ( (pPhi >= fSPhi - halfAngTolerance)
643  && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
644  {
645  in = kSurface ;
646  }
647  }
648  else // fSPhi < 0
649  {
650  if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
651  && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;}
652  else
653  {
654  in = kSurface ;
655  }
656  }
657  }
658  }
659  else // Try generous boundaries
660  {
661  tolRMin = fRmin - fRminTolerance ;
662  tolRMax = fRmax + fRmaxTolerance ;
663 
664  if (tolRMin < 0 ) { tolRMin = 0 ; }
665 
666  if ( (pt2 >= tolRMin*tolRMin) && (pt2 <= tolRMax*tolRMax) )
667  {
668  if ( (fDPhi == twopi) || (pt2 == 0) ) // Continuous in phi or on z-axis
669  {
670  in = kSurface ;
671  }
672  else // Try outer tolerant phi boundaries only
673  {
674  pPhi = std::atan2(p.y(),p.x()) ;
675 
676  if ( pPhi < -halfAngTolerance ) { pPhi += twopi ; } // 0<=pPhi<2pi
677  if ( fSPhi >= 0 )
678  {
679  if ( (std::fabs(pPhi) < halfAngTolerance)
680  && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
681  {
682  pPhi += twopi ; // 0 <= pPhi < 2pi
683  }
684  if ( (pPhi >= fSPhi - halfAngTolerance)
685  && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
686  {
687  in = kSurface;
688  }
689  }
690  else // fSPhi < 0
691  {
692  if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
693  && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;}
694  else
695  {
696  in = kSurface ;
697  }
698  }
699  }
700  }
701  }
702  return in ;
703 }
704 
706 //
707 // Return unit normal of surface closest to p
708 // - note if point on z axis, ignore phi divided sides
709 // - unsafe if point close to z axis a rmin=0 - no explicit checks
710 
712 {
713  G4int noSurfaces = 0;
714  G4double rho, pt, pPhi;
715  G4double distRMin = kInfinity;
716  G4double distSPhi = kInfinity, distEPhi = kInfinity;
717 
718  // To cope with precision loss
719  //
720  const G4double delta = std::max(10.0*kCarTolerance,
721  1.0e-8*(fRtor+fRmax));
722  const G4double dAngle = 10.0*kAngTolerance;
723 
724  G4ThreeVector nR, nPs, nPe;
725  G4ThreeVector norm, sumnorm(0.,0.,0.);
726 
727  rho = std::hypot(p.x(),p.y());
728  pt = std::hypot(p.z(),rho-fRtor);
729 
730  G4double distRMax = std::fabs(pt - fRmax);
731  if(fRmin) distRMin = std::fabs(pt - fRmin);
732 
733  if( rho > delta && pt != 0.0 )
734  {
735  G4double redFactor= (rho-fRtor)/rho;
736  nR = G4ThreeVector( p.x()*redFactor, // p.x()*(1.-fRtor/rho),
737  p.y()*redFactor, // p.y()*(1.-fRtor/rho),
738  p.z() );
739  nR *= 1.0/pt;
740  }
741 
742  if ( fDPhi < twopi ) // && rho ) // old limitation against (0,0,z)
743  {
744  if ( rho )
745  {
746  pPhi = std::atan2(p.y(),p.x());
747 
748  if(pPhi < fSPhi-delta) { pPhi += twopi; }
749  else if(pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; }
750 
751  distSPhi = std::fabs( pPhi - fSPhi );
752  distEPhi = std::fabs(pPhi-fSPhi-fDPhi);
753  }
754  nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
755  nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
756  }
757  if( distRMax <= delta )
758  {
759  noSurfaces ++;
760  sumnorm += nR;
761  }
762  else if( fRmin && (distRMin <= delta) ) // Must not be on both Outer and Inner
763  {
764  noSurfaces ++;
765  sumnorm -= nR;
766  }
767 
768  // To be on one of the 'phi' surfaces,
769  // it must be within the 'tube' - with tolerance
770 
771  if( (fDPhi < twopi) && (fRmin-delta <= pt) && (pt <= (fRmax+delta)) )
772  {
773  if (distSPhi <= dAngle)
774  {
775  noSurfaces ++;
776  sumnorm += nPs;
777  }
778  if (distEPhi <= dAngle)
779  {
780  noSurfaces ++;
781  sumnorm += nPe;
782  }
783  }
784  if ( noSurfaces == 0 )
785  {
786 #ifdef G4CSGDEBUG
788  ed.precision(16);
789 
790  EInside inIt= Inside( p );
791 
792  if( inIt != kSurface )
793  {
794  ed << " ERROR> Surface Normal was called for Torus,"
795  << " with point not on surface." << G4endl;
796  }
797  else
798  {
799  ed << " ERROR> Surface Normal has not found a surface, "
800  << " despite the point being on the surface. " <<G4endl;
801  }
802 
803  if( inIt != kInside)
804  {
805  ed << " Safety (Dist To In) = " << DistanceToIn(p) << G4endl;
806  }
807  if( inIt != kOutside)
808  {
809  ed << " Safety (Dist to Out) = " << DistanceToOut(p) << G4endl;
810  }
811  ed << " Coordinates of point : " << p << G4endl;
812  ed << " Parameters of solid : " << G4endl << *this << G4endl;
813 
814  if( inIt == kSurface )
815  {
816  G4Exception("G4Torus::SurfaceNormal(p)", "GeomSolids1002",
817  JustWarning, ed,
818  "Failing to find normal, even though point is on surface!");
819  }
820  else
821  {
822  static const char* NameInside[3]= { "Inside", "Surface", "Outside" };
823  ed << " The point is " << NameInside[inIt] << " the solid. "<< G4endl;
824  G4Exception("G4Torus::SurfaceNormal(p)", "GeomSolids1002",
825  JustWarning, ed, "Point p is not on surface !?" );
826  }
827 #endif
828  norm = ApproxSurfaceNormal(p);
829  }
830  else if ( noSurfaces == 1 ) { norm = sumnorm; }
831  else { norm = sumnorm.unit(); }
832 
833  return norm ;
834 }
835 
837 //
838 // Algorithm for SurfaceNormal() following the original specification
839 // for points not on the surface
840 
842 {
843  ENorm side ;
845  G4double rho,pt,phi;
846  G4double distRMin,distRMax,distSPhi,distEPhi,distMin;
847 
848  rho = std::hypot(p.x(),p.y());
849  pt = std::hypot(p.z(),rho-fRtor);
850 
851 #ifdef G4CSGDEBUG
852  G4cout << " G4Torus::ApproximateSurfaceNormal called for point " << p
853  << G4endl;
854 #endif
855 
856  distRMax = std::fabs(pt - fRmax) ;
857 
858  if(fRmin) // First minimum radius
859  {
860  distRMin = std::fabs(pt - fRmin) ;
861 
862  if (distRMin < distRMax)
863  {
864  distMin = distRMin ;
865  side = kNRMin ;
866  }
867  else
868  {
869  distMin = distRMax ;
870  side = kNRMax ;
871  }
872  }
873  else
874  {
875  distMin = distRMax ;
876  side = kNRMax ;
877  }
878  if ( (fDPhi < twopi) && rho )
879  {
880  phi = std::atan2(p.y(),p.x()) ; // Protected against (0,0,z) (above rho!=0)
881 
882  if (phi < 0) { phi += twopi ; }
883 
884  if (fSPhi < 0 ) { distSPhi = std::fabs(phi-(fSPhi+twopi))*rho ; }
885  else { distSPhi = std::fabs(phi-fSPhi)*rho ; }
886 
887  distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
888 
889  if (distSPhi < distEPhi) // Find new minimum
890  {
891  if (distSPhi<distMin) side = kNSPhi ;
892  }
893  else
894  {
895  if (distEPhi < distMin) { side = kNEPhi ; }
896  }
897  }
898  switch (side)
899  {
900  case kNRMin: // Inner radius
901  norm = G4ThreeVector( -p.x()*(1-fRtor/rho)/pt,
902  -p.y()*(1-fRtor/rho)/pt,
903  -p.z()/pt ) ;
904  break ;
905  case kNRMax: // Outer radius
906  norm = G4ThreeVector( p.x()*(1-fRtor/rho)/pt,
907  p.y()*(1-fRtor/rho)/pt,
908  p.z()/pt ) ;
909  break;
910  case kNSPhi:
911  norm = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
912  break;
913  case kNEPhi:
914  norm = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
915  break;
916  default: // Should never reach this case ...
917  DumpInfo();
918  G4Exception("G4Torus::ApproxSurfaceNormal()",
919  "GeomSolids1002", JustWarning,
920  "Undefined side for valid surface normal to solid.");
921  break ;
922  }
923  return norm ;
924 }
925 
927 //
928 // Calculate distance to shape from outside, along normalised vector
929 // - return kInfinity if no intersection, or intersection distance <= tolerance
930 //
931 // - Compute the intersection with the z planes
932 // - if at valid r, phi, return
933 //
934 // -> If point is outer outer radius, compute intersection with rmax
935 // - if at valid phi,z return
936 //
937 // -> Compute intersection with inner radius, taking largest +ve root
938 // - if valid (phi), save intersction
939 //
940 // -> If phi segmented, compute intersections with phi half planes
941 // - return smallest of valid phi intersections and
942 // inner radius intersection
943 //
944 // NOTE:
945 // - Precalculations for phi trigonometry are Done `just in time'
946 // - `if valid' implies tolerant checking of intersection points
947 
949  const G4ThreeVector& v ) const
950 {
951 
952  G4double snxt=kInfinity, sphi=kInfinity; // snxt = default return value
953 
954  G4double sd[4] ;
955 
956  // Precalculated trig for phi intersections - used by r,z intersections to
957  // check validity
958 
959  G4bool seg; // true if segmented
960  G4double hDPhi; // half dphi
961  G4double cPhi,sinCPhi=0.,cosCPhi=0.; // central phi
962 
963  G4double tolORMin2; // `generous' radii squared
964  G4double tolORMax2;
965 
966  G4double Dist,xi,yi,zi,rhoi,it2; // Intersection point variables
967 
968  G4double Comp;
969  G4double cosSPhi,sinSPhi; // Trig for phi start intersect
970  G4double ePhi,cosEPhi,sinEPhi; // for phi end intersect
971 
972  // Set phi divided flag and precalcs
973  //
974  if ( fDPhi < twopi )
975  {
976  seg = true ;
977  hDPhi = 0.5*fDPhi ; // half delta phi
978  cPhi = fSPhi + hDPhi ;
979  sinCPhi = std::sin(cPhi) ;
980  cosCPhi = std::cos(cPhi) ;
981  }
982  else
983  {
984  seg = false ;
985  }
986 
987  if (fRmin > fRminTolerance) // Calculate tolerant rmin and rmax
988  {
989  tolORMin2 = (fRmin - fRminTolerance)*(fRmin - fRminTolerance) ;
990  }
991  else
992  {
993  tolORMin2 = 0 ;
994  }
995  tolORMax2 = (fRmax + fRmaxTolerance)*(fRmax + fRmaxTolerance) ;
996 
997  // Intersection with Rmax (possible return) and Rmin (must also check phi)
998 
999  snxt = SolveNumericJT(p,v,fRmax,true);
1000 
1001  if (fRmin) // Possible Rmin intersection
1002  {
1003  sd[0] = SolveNumericJT(p,v,fRmin,true);
1004  if ( sd[0] < snxt ) { snxt = sd[0] ; }
1005  }
1006 
1007  //
1008  // Phi segment intersection
1009  //
1010  // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1011  //
1012  // o NOTE: Large duplication of code between sphi & ephi checks
1013  // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1014  // intersection check <=0 -> >=0
1015  // -> use some form of loop Construct ?
1016 
1017  if (seg)
1018  {
1019  sinSPhi = std::sin(fSPhi) ; // First phi surface ('S'tarting phi)
1020  cosSPhi = std::cos(fSPhi) ;
1021  Comp = v.x()*sinSPhi - v.y()*cosSPhi ; // Component in outwards
1022  // normal direction
1023  if (Comp < 0 )
1024  {
1025  Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1026 
1027  if (Dist < halfCarTolerance)
1028  {
1029  sphi = Dist/Comp ;
1030  if (sphi < snxt)
1031  {
1032  if ( sphi < 0 ) { sphi = 0 ; }
1033 
1034  xi = p.x() + sphi*v.x() ;
1035  yi = p.y() + sphi*v.y() ;
1036  zi = p.z() + sphi*v.z() ;
1037  rhoi = std::hypot(xi,yi);
1038  it2 = zi*zi + (rhoi-fRtor)*(rhoi-fRtor);
1039 
1040  if ( it2 >= tolORMin2 && it2 <= tolORMax2 )
1041  {
1042  // r intersection is good - check intersecting
1043  // with correct half-plane
1044  //
1045  if ((yi*cosCPhi-xi*sinCPhi)<=0) { snxt=sphi; }
1046  }
1047  }
1048  }
1049  }
1050  ePhi=fSPhi+fDPhi; // Second phi surface ('E'nding phi)
1051  sinEPhi=std::sin(ePhi);
1052  cosEPhi=std::cos(ePhi);
1053  Comp=-(v.x()*sinEPhi-v.y()*cosEPhi);
1054 
1055  if ( Comp < 0 ) // Component in outwards normal dirn
1056  {
1057  Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1058 
1059  if (Dist < halfCarTolerance )
1060  {
1061  sphi = Dist/Comp ;
1062 
1063  if (sphi < snxt )
1064  {
1065  if (sphi < 0 ) { sphi = 0 ; }
1066 
1067  xi = p.x() + sphi*v.x() ;
1068  yi = p.y() + sphi*v.y() ;
1069  zi = p.z() + sphi*v.z() ;
1070  rhoi = std::hypot(xi,yi);
1071  it2 = zi*zi + (rhoi-fRtor)*(rhoi-fRtor);
1072 
1073  if (it2 >= tolORMin2 && it2 <= tolORMax2)
1074  {
1075  // z and r intersections good - check intersecting
1076  // with correct half-plane
1077  //
1078  if ((yi*cosCPhi-xi*sinCPhi)>=0) { snxt=sphi; }
1079  }
1080  }
1081  }
1082  }
1083  }
1084  if(snxt < halfCarTolerance) { snxt = 0.0 ; }
1085 
1086  return snxt ;
1087 }
1088 
1090 //
1091 // Calculate distance (<= actual) to closest surface of shape from outside
1092 // - Calculate distance to z, radial planes
1093 // - Only to phi planes if outside phi extent
1094 // - Return 0 if point inside
1095 
1097 {
1098  G4double safe=0.0, safe1, safe2 ;
1099  G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
1100  G4double rho, pt ;
1101 
1102  rho = std::hypot(p.x(),p.y());
1103  pt = std::hypot(p.z(),rho-fRtor);
1104  safe1 = fRmin - pt ;
1105  safe2 = pt - fRmax ;
1106 
1107  if (safe1 > safe2) { safe = safe1; }
1108  else { safe = safe2; }
1109 
1110  if ( fDPhi < twopi && rho )
1111  {
1112  phiC = fSPhi + fDPhi*0.5 ;
1113  cosPhiC = std::cos(phiC) ;
1114  sinPhiC = std::sin(phiC) ;
1115  cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
1116 
1117  if (cosPsi < std::cos(fDPhi*0.5) ) // Psi=angle from central phi to point
1118  { // Point lies outside phi range
1119  if ((p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
1120  {
1121  safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
1122  }
1123  else
1124  {
1125  ePhi = fSPhi + fDPhi ;
1126  safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1127  }
1128  if (safePhi > safe) { safe = safePhi ; }
1129  }
1130  }
1131  if (safe < 0 ) { safe = 0 ; }
1132  return safe;
1133 }
1134 
1136 //
1137 // Calculate distance to surface of shape from `inside', allowing for tolerance
1138 // - Only Calc rmax intersection if no valid rmin intersection
1139 //
1140 
1142  const G4ThreeVector& v,
1143  const G4bool calcNorm,
1144  G4bool *validNorm,
1145  G4ThreeVector *n ) const
1146 {
1147  ESide side = kNull, sidephi = kNull ;
1148  G4double snxt = kInfinity, sphi, sd[4] ;
1149 
1150  // Vars for phi intersection
1151  //
1152  G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi;
1153  G4double cPhi, sinCPhi, cosCPhi ;
1154  G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, zi, vphi ;
1155 
1156  // Radial Intersections Defenitions & General Precals
1157 
1159 
1160 #if 1
1161 
1162  // This is the version with the calculation of CalcNorm = true
1163  // To be done: Check the precision of this calculation.
1164  // If you want return always validNorm = false, then take the version below
1165 
1166 
1167  G4double rho = std::hypot(p.x(),p.y());
1168  G4double pt = hypot(p.z(),rho-fRtor);
1169 
1170  G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
1171 
1172  G4double tolRMax = fRmax - fRmaxTolerance ;
1173 
1174  G4double vDotNmax = pDotV - fRtor*(v.x()*p.x() + v.y()*p.y())/rho ;
1175  G4double pDotxyNmax = (1 - fRtor/rho) ;
1176 
1177  if( (pt*pt > tolRMax*tolRMax) && (vDotNmax >= 0) )
1178  {
1179  // On tolerant boundary & heading outwards (or perpendicular to) outer
1180  // radial surface -> leaving immediately with *n for really convex part
1181  // only
1182 
1183  if ( calcNorm && (pDotxyNmax >= -2.*fRmaxTolerance) )
1184  {
1185  *n = G4ThreeVector( p.x()*(1 - fRtor/rho)/pt,
1186  p.y()*(1 - fRtor/rho)/pt,
1187  p.z()/pt ) ;
1188  *validNorm = true ;
1189  }
1190 
1191  return snxt = 0 ; // Leaving by Rmax immediately
1192  }
1193 
1194  snxt = SolveNumericJT(p,v,fRmax,false);
1195  side = kRMax ;
1196 
1197  // rmin
1198 
1199  if ( fRmin )
1200  {
1201  G4double tolRMin = fRmin + fRminTolerance ;
1202 
1203  if ( (pt*pt < tolRMin*tolRMin) && (vDotNmax < 0) )
1204  {
1205  if (calcNorm) { *validNorm = false ; } // Concave surface of the torus
1206  return snxt = 0 ; // Leaving by Rmin immediately
1207  }
1208 
1209  sd[0] = SolveNumericJT(p,v,fRmin,false);
1210  if ( sd[0] < snxt )
1211  {
1212  snxt = sd[0] ;
1213  side = kRMin ;
1214  }
1215  }
1216 
1217 #else
1218 
1219  // this is the "conservative" version which return always validnorm = false
1220  // NOTE: using this version the unit test testG4Torus will break
1221 
1222  snxt = SolveNumericJT(p,v,fRmax,false);
1223  side = kRMax ;
1224 
1225  if ( fRmin )
1226  {
1227  sd[0] = SolveNumericJT(p,v,fRmin,false);
1228  if ( sd[0] < snxt )
1229  {
1230  snxt = sd[0] ;
1231  side = kRMin ;
1232  }
1233  }
1234 
1235  if ( calcNorm && (snxt == 0.0) )
1236  {
1237  *validNorm = false ; // Leaving solid, but possible re-intersection
1238  return snxt ;
1239  }
1240 
1241 #endif
1242 
1243  if (fDPhi < twopi) // Phi Intersections
1244  {
1245  sinSPhi = std::sin(fSPhi) ;
1246  cosSPhi = std::cos(fSPhi) ;
1247  ePhi = fSPhi + fDPhi ;
1248  sinEPhi = std::sin(ePhi) ;
1249  cosEPhi = std::cos(ePhi) ;
1250  cPhi = fSPhi + fDPhi*0.5 ;
1251  sinCPhi = std::sin(cPhi) ;
1252  cosCPhi = std::cos(cPhi) ;
1253 
1254  // angle calculation with correction
1255  // of difference in domain of atan2 and Sphi
1256  //
1257  vphi = std::atan2(v.y(),v.x()) ;
1258 
1259  if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1260  else if ( vphi > ePhi + halfAngTolerance ) { vphi -= twopi; }
1261 
1262  if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1263  {
1264  pDistS = p.x()*sinSPhi - p.y()*cosSPhi ; // pDist -ve when inside
1265  pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1266 
1267  // Comp -ve when in direction of outwards normal
1268  //
1269  compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1270  compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1271  sidephi = kNull ;
1272 
1273  if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1274  && (pDistE <= halfCarTolerance) ) )
1275  || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
1276  && (pDistE > halfCarTolerance) ) ) )
1277  {
1278  // Inside both phi *full* planes
1279 
1280  if ( compS < 0 )
1281  {
1282  sphi = pDistS/compS ;
1283 
1284  if (sphi >= -halfCarTolerance)
1285  {
1286  xi = p.x() + sphi*v.x() ;
1287  yi = p.y() + sphi*v.y() ;
1288 
1289  // Check intersecting with correct half-plane
1290  // (if not -> no intersect)
1291  //
1292  if ( (std::fabs(xi)<=kCarTolerance)
1293  && (std::fabs(yi)<=kCarTolerance) )
1294  {
1295  sidephi = kSPhi;
1296  if ( ((fSPhi-halfAngTolerance)<=vphi)
1297  && ((ePhi+halfAngTolerance)>=vphi) )
1298  {
1299  sphi = kInfinity;
1300  }
1301  }
1302  else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1303  {
1304  sphi = kInfinity ;
1305  }
1306  else
1307  {
1308  sidephi = kSPhi ;
1309  }
1310  }
1311  else
1312  {
1313  sphi = kInfinity ;
1314  }
1315  }
1316  else
1317  {
1318  sphi = kInfinity ;
1319  }
1320 
1321  if ( compE < 0 )
1322  {
1323  sphi2 = pDistE/compE ;
1324 
1325  // Only check further if < starting phi intersection
1326  //
1327  if ( (sphi2 > -kCarTolerance) && (sphi2 < sphi) )
1328  {
1329  xi = p.x() + sphi2*v.x() ;
1330  yi = p.y() + sphi2*v.y() ;
1331 
1332  if ( (std::fabs(xi)<=kCarTolerance)
1333  && (std::fabs(yi)<=kCarTolerance) )
1334  {
1335  // Leaving via ending phi
1336  //
1337  if( !( (fSPhi-halfAngTolerance <= vphi)
1338  && (ePhi+halfAngTolerance >= vphi) ) )
1339  {
1340  sidephi = kEPhi ;
1341  sphi = sphi2;
1342  }
1343  }
1344  else // Check intersecting with correct half-plane
1345  {
1346  if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1347  {
1348  // Leaving via ending phi
1349  //
1350  sidephi = kEPhi ;
1351  sphi = sphi2;
1352 
1353  }
1354  }
1355  }
1356  }
1357  }
1358  else
1359  {
1360  sphi = kInfinity ;
1361  }
1362  }
1363  else
1364  {
1365  // On z axis + travel not || to z axis -> if phi of vector direction
1366  // within phi of shape, Step limited by rmax, else Step =0
1367 
1368  vphi = std::atan2(v.y(),v.x());
1369 
1370  if ( ( fSPhi-halfAngTolerance <= vphi ) &&
1371  ( vphi <= ( ePhi+halfAngTolerance ) ) )
1372  {
1373  sphi = kInfinity;
1374  }
1375  else
1376  {
1377  sidephi = kSPhi ; // arbitrary
1378  sphi=0;
1379  }
1380  }
1381 
1382  // Order intersections
1383 
1384  if (sphi<snxt)
1385  {
1386  snxt=sphi;
1387  side=sidephi;
1388  }
1389  }
1390 
1391  G4double rhoi,it,iDotxyNmax ;
1392  // Note: by numerical computation we know where the ray hits the torus
1393  // So I propose to return the side where the ray hits
1394 
1395  if (calcNorm)
1396  {
1397  switch(side)
1398  {
1399  case kRMax: // n is unit vector
1400  xi = p.x() + snxt*v.x() ;
1401  yi = p.y() + snxt*v.y() ;
1402  zi = p.z() + snxt*v.z() ;
1403  rhoi = std::hypot(xi,yi);
1404  it = hypot(zi,rhoi-fRtor);
1405 
1406  iDotxyNmax = (1-fRtor/rhoi) ;
1407  if(iDotxyNmax >= -2.*fRmaxTolerance) // really convex part of Rmax
1408  {
1409  *n = G4ThreeVector( xi*(1-fRtor/rhoi)/it,
1410  yi*(1-fRtor/rhoi)/it,
1411  zi/it ) ;
1412  *validNorm = true ;
1413  }
1414  else
1415  {
1416  *validNorm = false ; // concave-convex part of Rmax
1417  }
1418  break ;
1419 
1420  case kRMin:
1421  *validNorm = false ; // Rmin is concave or concave-convex
1422  break;
1423 
1424  case kSPhi:
1425  if (fDPhi <= pi )
1426  {
1427  *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
1428  *validNorm=true;
1429  }
1430  else
1431  {
1432  *validNorm = false ;
1433  }
1434  break ;
1435 
1436  case kEPhi:
1437  if (fDPhi <= pi)
1438  {
1439  *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
1440  *validNorm=true;
1441  }
1442  else
1443  {
1444  *validNorm = false ;
1445  }
1446  break;
1447 
1448  default:
1449 
1450  // It seems we go here from time to time ...
1451 
1452  G4cout << G4endl;
1453  DumpInfo();
1454  std::ostringstream message;
1455  G4int oldprc = message.precision(16);
1456  message << "Undefined side for valid surface normal to solid."
1457  << G4endl
1458  << "Position:" << G4endl << G4endl
1459  << "p.x() = " << p.x()/mm << " mm" << G4endl
1460  << "p.y() = " << p.y()/mm << " mm" << G4endl
1461  << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
1462  << "Direction:" << G4endl << G4endl
1463  << "v.x() = " << v.x() << G4endl
1464  << "v.y() = " << v.y() << G4endl
1465  << "v.z() = " << v.z() << G4endl << G4endl
1466  << "Proposed distance :" << G4endl << G4endl
1467  << "snxt = " << snxt/mm << " mm" << G4endl;
1468  message.precision(oldprc);
1469  G4Exception("G4Torus::DistanceToOut(p,v,..)",
1470  "GeomSolids1002",JustWarning, message);
1471  break;
1472  }
1473  }
1474  if ( snxt<halfCarTolerance ) { snxt=0 ; }
1475 
1476  return snxt;
1477 }
1478 
1480 //
1481 // Calculate distance (<=actual) to closest surface of shape from inside
1482 
1484 {
1485  G4double safe=0.0,safeR1,safeR2;
1486  G4double rho,pt ;
1487  G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi;
1488 
1489  rho = std::hypot(p.x(),p.y());
1490  pt = std::hypot(p.z(),rho-fRtor);
1491 
1492 #ifdef G4CSGDEBUG
1493  if( Inside(p) == kOutside )
1494  {
1495  G4int oldprc = G4cout.precision(16) ;
1496  G4cout << G4endl ;
1497  DumpInfo();
1498  G4cout << "Position:" << G4endl << G4endl ;
1499  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1500  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1501  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1502  G4cout.precision(oldprc);
1503  G4Exception("G4Torus::DistanceToOut(p)", "GeomSolids1002",
1504  JustWarning, "Point p is outside !?" );
1505  }
1506 #endif
1507 
1508  if (fRmin)
1509  {
1510  safeR1 = pt - fRmin ;
1511  safeR2 = fRmax - pt ;
1512 
1513  if (safeR1 < safeR2) { safe = safeR1 ; }
1514  else { safe = safeR2 ; }
1515  }
1516  else
1517  {
1518  safe = fRmax - pt ;
1519  }
1520 
1521  // Check if phi divided, Calc distances closest phi plane
1522  //
1523  if (fDPhi<twopi) // Above/below central phi of Torus?
1524  {
1525  phiC = fSPhi + fDPhi*0.5 ;
1526  cosPhiC = std::cos(phiC) ;
1527  sinPhiC = std::sin(phiC) ;
1528 
1529  if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
1530  {
1531  safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
1532  }
1533  else
1534  {
1535  ePhi = fSPhi + fDPhi ;
1536  safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1537  }
1538  if (safePhi < safe) { safe = safePhi ; }
1539  }
1540  if (safe < 0) { safe = 0 ; }
1541  return safe ;
1542 }
1543 
1545 //
1546 // Stream object contents to an output stream
1547 
1549 {
1550  return G4String("G4Torus");
1551 }
1552 
1554 //
1555 // Make a clone of the object
1556 //
1558 {
1559  return new G4Torus(*this);
1560 }
1561 
1563 //
1564 // Stream object contents to an output stream
1565 
1566 std::ostream& G4Torus::StreamInfo( std::ostream& os ) const
1567 {
1568  G4int oldprc = os.precision(16);
1569  os << "-----------------------------------------------------------\n"
1570  << " *** Dump for solid - " << GetName() << " ***\n"
1571  << " ===================================================\n"
1572  << " Solid type: G4Torus\n"
1573  << " Parameters: \n"
1574  << " inner radius: " << fRmin/mm << " mm \n"
1575  << " outer radius: " << fRmax/mm << " mm \n"
1576  << " swept radius: " << fRtor/mm << " mm \n"
1577  << " starting phi: " << fSPhi/degree << " degrees \n"
1578  << " delta phi : " << fDPhi/degree << " degrees \n"
1579  << "-----------------------------------------------------------\n";
1580  os.precision(oldprc);
1581 
1582  return os;
1583 }
1584 
1586 //
1587 // GetPointOnSurface
1588 
1590 {
1591  G4double cosu, sinu,cosv, sinv, aOut, aIn, aSide, chose, phi, theta, rRand;
1592 
1594  theta = G4RandFlat::shoot(0.,twopi);
1595 
1596  cosu = std::cos(phi); sinu = std::sin(phi);
1597  cosv = std::cos(theta); sinv = std::sin(theta);
1598 
1599  // compute the areas
1600 
1601  aOut = (fDPhi)*twopi*fRtor*fRmax;
1602  aIn = (fDPhi)*twopi*fRtor*fRmin;
1603  aSide = pi*(fRmax*fRmax-fRmin*fRmin);
1604 
1605  if ((fSPhi == 0) && (fDPhi == twopi)){ aSide = 0; }
1606  chose = G4RandFlat::shoot(0.,aOut + aIn + 2.*aSide);
1607 
1608  if(chose < aOut)
1609  {
1610  return G4ThreeVector ((fRtor+fRmax*cosv)*cosu,
1611  (fRtor+fRmax*cosv)*sinu, fRmax*sinv);
1612  }
1613  else if( (chose >= aOut) && (chose < aOut + aIn) )
1614  {
1615  return G4ThreeVector ((fRtor+fRmin*cosv)*cosu,
1616  (fRtor+fRmin*cosv)*sinu, fRmin*sinv);
1617  }
1618  else if( (chose >= aOut + aIn) && (chose < aOut + aIn + aSide) )
1619  {
1620  rRand = GetRadiusInRing(fRmin,fRmax);
1621  return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi),
1622  (fRtor+rRand*cosv)*std::sin(fSPhi), rRand*sinv);
1623  }
1624  else
1625  {
1626  rRand = GetRadiusInRing(fRmin,fRmax);
1627  return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi+fDPhi),
1628  (fRtor+rRand*cosv)*std::sin(fSPhi+fDPhi),
1629  rRand*sinv);
1630  }
1631 }
1632 
1634 //
1635 // Visualisation Functions
1636 
1638 {
1639  scene.AddSolid (*this);
1640 }
1641 
1643 {
1644  return new G4PolyhedronTorus (fRmin, fRmax, fRtor, fSPhi, fDPhi);
1645 }
1646 
1647 #endif // !defined(G4GEOM_USE_TORUS) || !defined(G4GEOM_USE_SYS_USOLIDS)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
void SetAllParameters(G4double pRmin, G4double pRmax, G4double pRtor, G4double pSPhi, G4double pDPhi)
Definition: G4Torus.cc:103
void set(double x, double y, double z)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double GetCosEndPhi() const
G4double GetRtor() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetMaxExtent(const EAxis pAxis) const
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
G4double fSPhi
Definition: G4Torus.hh:195
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void TorusRootsJT(const G4ThreeVector &p, const G4ThreeVector &v, G4double r, std::vector< G4double > &roots) const
Definition: G4Torus.cc:262
virtual void AddSolid(const G4Box &)=0
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Torus.cc:948
static constexpr double mm
Definition: G4SIunits.hh:115
EInside Inside(const G4ThreeVector &p) const
Definition: G4Torus.cc:600
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:92
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Torus.cc:711
#define G4endl
Definition: G4ios.hh:61
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4Torus.cc:417
G4GeometryType GetEntityType() const
Definition: G4Torus.cc:1548
const char * p
Definition: xmltok.h:285
void message(RunManager *runmanager)
Definition: ts_scorers.cc:72
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:112
G4VSolid * Clone() const
Definition: G4Torus.cc:1557
G4double fDPhi
Definition: G4Torus.hh:195
G4double GetAngularTolerance() const
double z() const
G4double GetRadialTolerance() const
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4Torus.cc:841
G4double GetRmax() const
G4double fRmin
Definition: G4Torus.hh:195
static const G4double emax
G4double fRmaxTolerance
Definition: G4Torus.hh:203
void setX(double)
G4double halfCarTolerance
Definition: G4Torus.hh:206
static G4bool DiskExtent(G4double rmin, G4double rmax, G4double startPhi, G4double delPhi, G4TwoVector &pmin, G4TwoVector &pmax)
Definition: G4GeomTools.cc:396
G4Polyhedron * CreatePolyhedron() const
Definition: G4Torus.cc:1642
G4double SolveNumericJT(const G4ThreeVector &p, const G4ThreeVector &v, G4double r, G4bool IsDistanceToIn) const
Definition: G4Torus.cc:302
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4double kAngTolerance
Definition: G4Torus.hh:203
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
std::vector< G4ThreeVector > G4ThreeVectorList
G4double halfAngTolerance
Definition: G4Torus.hh:206
G4Torus(const G4String &pName, G4double pRmin, G4double pRmax, G4double pRtor, G4double pSPhi, G4double pDPhi)
Definition: G4Torus.cc:87
G4int FindRoots(G4double *op, G4int degree, G4double *zeror, G4double *zeroi)
G4double GetMinExtent(const EAxis pAxis) const
double x() const
static constexpr double deg
Definition: G4SIunits.hh:152
G4double kCarTolerance
Definition: G4VSolid.hh:307
static constexpr double twopi
Definition: G4SIunits.hh:76
TMarker * pt
Definition: egs.C:25
G4ThreeVector GetPointOnSurface() const
Definition: G4Torus.cc:1589
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Torus.cc:1566
Float_t d
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Torus.cc:248
G4double kRadTolerance
Definition: G4Torus.hh:203
ThreeVector shoot(const G4int Ap, const G4int Af)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Torus.cc:462
G4double GetSinStartPhi() const
Hep3Vector unit() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
G4double GetSinEndPhi() const
int G4int
Definition: G4Types.hh:78
EInside
Definition: geomdefs.hh:58
ifstream in
Definition: comparison.C:7
~G4Torus()
Definition: G4Torus.cc:200
EAxis
Definition: geomdefs.hh:54
void set(double x, double y)
G4String GetName() const
void DumpInfo() const
Float_t norm
double y() const
G4double GetRmin() const
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Torus.cc:1637
static constexpr double degree
Definition: G4SIunits.hh:144
G4GLOB_DLL std::ostream G4cout
double x() const
G4double GetDPhi() const
static G4GeometryTolerance * GetInstance()
G4double fRmax
Definition: G4Torus.hh:195
Char_t n[5]
static constexpr double pi
Definition: G4SIunits.hh:75
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Torus.cc:1141
double y() const
G4double fRtor
Definition: G4Torus.hh:195
G4double fSurfaceArea
Definition: G4CSGSolid.hh:79
G4double fCubicVolume
Definition: G4CSGSolid.hh:78
G4bool fRebuildPolyhedron
Definition: G4CSGSolid.hh:80
void setY(double)
G4double fRminTolerance
Definition: G4Torus.hh:203
G4double GetCosStartPhi() const
G4Torus & operator=(const G4Torus &rhs)
Definition: G4Torus.cc:221