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