Geant4  10.03
G4Tubs.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: G4Tubs.cc 101121 2016-11-07 09:18:01Z gcosmo $
28 //
29 //
30 // class G4Tubs
31 //
32 // History:
33 //
34 // 24.08.16 E.Tcherniaev: reimplemented CalculateExtent() to make use
35 // of G4BoundingEnvelope
36 // 05.04.12 M.Kelsey: Use sqrt(r) in GetPointOnSurface() for uniform points
37 // 02.08.07 T.Nikitina: bug fixed in DistanceToOut(p,v,..) for negative value under sqrt
38 // for the case: p on the surface and v is tangent to the surface
39 // 11.05.07 T.Nikitina: bug fixed in DistanceToOut(p,v,..) for phi < 2pi
40 // 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
41 // 16.03.05 V.Grichine: SurfaceNormal(p) with edges/corners for boolean
42 // 20.07.01 V.Grichine: bug fixed in Inside(p)
43 // 20.02.01 V.Grichine: bug fixed in Inside(p) and CalculateExtent was
44 // simplified base on G4Box::CalculateExtent
45 // 07.12.00 V.Grichine: phi-section algorithm was changed in Inside(p)
46 // 28.11.00 V.Grichine: bug fixed in Inside(p)
47 // 31.10.00 V.Grichine: assign srd, sphi in Distance ToOut(p,v,...)
48 // 08.08.00 V.Grichine: more stable roots of 2-equation in DistanceToOut(p,v,..)
49 // 02.08.00 V.Grichine: point is outside check in Distance ToOut(p)
50 // 17.05.00 V.Grichine: bugs (#76,#91) fixed in Distance ToOut(p,v,...)
51 // 31.03.00 V.Grichine: bug fixed in Inside(p)
52 // 19.11.99 V.Grichine: side = kNull in DistanceToOut(p,v,...)
53 // 13.10.99 V.Grichine: bugs fixed in DistanceToIn(p,v)
54 // 28.05.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...)
55 // 25.05.99 V.Grichine: bugs fixed in DistanceToIn(p,v)
56 // 23.03.99 V.Grichine: bug fixed in DistanceToIn(p,v)
57 // 09.10.98 V.Grichine: modifications in DistanceToOut(p,v,...)
58 // 18.06.98 V.Grichine: n-normalisation in DistanceToOut(p,v)
59 //
60 // 1994-95 P.Kent: implementation
61 //
63 
64 #include "G4Tubs.hh"
65 
66 #if !defined(G4GEOM_USE_UTUBS)
67 
68 #include "G4GeomTools.hh"
69 #include "G4VoxelLimits.hh"
70 #include "G4AffineTransform.hh"
71 #include "G4GeometryTolerance.hh"
72 #include "G4BoundingEnvelope.hh"
73 
74 #include "G4VPVParameterisation.hh"
75 
76 #include "Randomize.hh"
77 
78 #include "meshdefs.hh"
79 
80 #include "G4VGraphicsScene.hh"
81 
82 using namespace CLHEP;
83 
85 //
86 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
87 // - note if pdphi>2PI then reset to 2PI
88 
89 G4Tubs::G4Tubs( const G4String &pName,
90  G4double pRMin, G4double pRMax,
91  G4double pDz,
92  G4double pSPhi, G4double pDPhi )
93  : G4CSGSolid(pName), fRMin(pRMin), fRMax(pRMax), fDz(pDz), fSPhi(0), fDPhi(0)
94 {
95 
98 
102 
103  if (pDz<=0) // Check z-len
104  {
105  std::ostringstream message;
106  message << "Negative Z half-length (" << pDz << ") in solid: " << GetName();
107  G4Exception("G4Tubs::G4Tubs()", "GeomSolids0002", FatalException, message);
108  }
109  if ( (pRMin >= pRMax) || (pRMin < 0) ) // Check radii
110  {
111  std::ostringstream message;
112  message << "Invalid values for radii in solid: " << GetName()
113  << G4endl
114  << " pRMin = " << pRMin << ", pRMax = " << pRMax;
115  G4Exception("G4Tubs::G4Tubs()", "GeomSolids0002", FatalException, message);
116  }
117 
118  // Check angles
119  //
120  CheckPhiAngles(pSPhi, pDPhi);
121 }
122 
124 //
125 // Fake default constructor - sets only member data and allocates memory
126 // for usage restricted to object persistency.
127 //
128 G4Tubs::G4Tubs( __void__& a )
129  : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
130  fRMin(0.), fRMax(0.), fDz(0.), fSPhi(0.), fDPhi(0.),
131  sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
132  sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
133  fPhiFullTube(false), halfCarTolerance(0.), halfRadTolerance(0.),
134  halfAngTolerance(0.)
135 {
136 }
137 
139 //
140 // Destructor
141 
143 {
144 }
145 
147 //
148 // Copy constructor
149 
151  : G4CSGSolid(rhs),
152  kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
153  fRMin(rhs.fRMin), fRMax(rhs.fRMax), fDz(rhs.fDz),
154  fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
155  sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
156  cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
157  sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
158  sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi), fPhiFullTube(rhs.fPhiFullTube),
159  halfCarTolerance(rhs.halfCarTolerance),
160  halfRadTolerance(rhs.halfRadTolerance),
161  halfAngTolerance(rhs.halfAngTolerance)
162 {
163 }
164 
166 //
167 // Assignment operator
168 
170 {
171  // Check assignment to self
172  //
173  if (this == &rhs) { return *this; }
174 
175  // Copy base class data
176  //
178 
179  // Copy data
180  //
182  fRMin = rhs.fRMin; fRMax = rhs.fRMax; fDz = rhs.fDz;
183  fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
184  sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
186  sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
187  sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
192 
193  return *this;
194 }
195 
197 //
198 // Dispatch to parameterisation for replication mechanism dimension
199 // computation & modification.
200 
202  const G4int n,
203  const G4VPhysicalVolume* pRep )
204 {
205  p->ComputeDimensions(*this,n,pRep) ;
206 }
207 
209 //
210 // Get bounding box
211 
213 {
214  G4double rmin = GetInnerRadius();
215  G4double rmax = GetOuterRadius();
216  G4double dz = GetZHalfLength();
217 
218  // Find bounding box
219  //
220  if (GetDeltaPhiAngle() < twopi)
221  {
222  G4TwoVector vmin,vmax;
223  G4GeomTools::DiskExtent(rmin,rmax,
226  vmin,vmax);
227  pMin.set(vmin.x(),vmin.y(),-dz);
228  pMax.set(vmax.x(),vmax.y(), dz);
229  }
230  else
231  {
232  pMin.set(-rmax,-rmax,-dz);
233  pMax.set( rmax, rmax, dz);
234  }
235 
236  // Check correctness of the bounding box
237  //
238  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
239  {
240  std::ostringstream message;
241  message << "Bad bounding box (min >= max) for solid: "
242  << GetName() << " !"
243  << "\npMin = " << pMin
244  << "\npMax = " << pMax;
245  G4Exception("G4Tubs::Extent()", "GeomMgt0001", JustWarning, message);
246  DumpInfo();
247  }
248 }
249 
251 //
252 // Calculate extent under transform and specified limit
253 
255  const G4VoxelLimits& pVoxelLimit,
256  const G4AffineTransform& pTransform,
257  G4double& pMin,
258  G4double& pMax ) const
259 {
260  G4ThreeVector bmin, bmax;
261  G4bool exist;
262 
263  // Get bounding box
264  Extent(bmin,bmax);
265 
266  // Check bounding box
267  G4BoundingEnvelope bbox(bmin,bmax);
268 #ifdef G4BBOX_EXTENT
269  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
270 #endif
271  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
272  {
273  return exist = (pMin < pMax) ? true : false;
274  }
275 
276  // Get parameters of the solid
277  G4double rmin = GetInnerRadius();
278  G4double rmax = GetOuterRadius();
279  G4double dz = GetZHalfLength();
280  G4double dphi = GetDeltaPhiAngle();
281 
282  // Find bounding envelope and calculate extent
283  //
284  const G4int NSTEPS = 24; // number of steps for whole circle
285  G4double astep = (360/NSTEPS)*deg; // max angle for one step
286  G4int ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
287  G4double ang = dphi/ksteps;
288 
289  G4double sinHalf = std::sin(0.5*ang);
290  G4double cosHalf = std::cos(0.5*ang);
291  G4double sinStep = 2.*sinHalf*cosHalf;
292  G4double cosStep = 1. - 2.*sinHalf*sinHalf;
293  G4double rext = rmax/cosHalf;
294 
295  // bounding envelope for full cylinder consists of two polygons,
296  // in other cases it is a sequence of quadrilaterals
297  if (rmin == 0 && dphi == twopi)
298  {
299  G4double sinCur = sinHalf;
300  G4double cosCur = cosHalf;
301 
302  G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS);
303  for (G4int k=0; k<NSTEPS; ++k)
304  {
305  baseA[k].set(rext*cosCur,rext*sinCur,-dz);
306  baseB[k].set(rext*cosCur,rext*sinCur, dz);
307 
308  G4double sinTmp = sinCur;
309  sinCur = sinCur*cosStep + cosCur*sinStep;
310  cosCur = cosCur*cosStep - sinTmp*sinStep;
311  }
312  std::vector<const G4ThreeVectorList *> polygons(2);
313  polygons[0] = &baseA;
314  polygons[1] = &baseB;
315  G4BoundingEnvelope benv(bmin,bmax,polygons);
316  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
317  }
318  else
319  {
320  G4double sinStart = GetSinStartPhi();
321  G4double cosStart = GetCosStartPhi();
322  G4double sinEnd = GetSinEndPhi();
323  G4double cosEnd = GetCosEndPhi();
324  G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
325  G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
326 
327  // set quadrilaterals
328  G4ThreeVectorList pols[NSTEPS+2];
329  for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(4);
330  pols[0][0].set(rmin*cosStart,rmin*sinStart, dz);
331  pols[0][1].set(rmin*cosStart,rmin*sinStart,-dz);
332  pols[0][2].set(rmax*cosStart,rmax*sinStart,-dz);
333  pols[0][3].set(rmax*cosStart,rmax*sinStart, dz);
334  for (G4int k=1; k<ksteps+1; ++k)
335  {
336  pols[k][0].set(rmin*cosCur,rmin*sinCur, dz);
337  pols[k][1].set(rmin*cosCur,rmin*sinCur,-dz);
338  pols[k][2].set(rext*cosCur,rext*sinCur,-dz);
339  pols[k][3].set(rext*cosCur,rext*sinCur, dz);
340 
341  G4double sinTmp = sinCur;
342  sinCur = sinCur*cosStep + cosCur*sinStep;
343  cosCur = cosCur*cosStep - sinTmp*sinStep;
344  }
345  pols[ksteps+1][0].set(rmin*cosEnd,rmin*sinEnd, dz);
346  pols[ksteps+1][1].set(rmin*cosEnd,rmin*sinEnd,-dz);
347  pols[ksteps+1][2].set(rmax*cosEnd,rmax*sinEnd,-dz);
348  pols[ksteps+1][3].set(rmax*cosEnd,rmax*sinEnd, dz);
349 
350  // set envelope and calculate extent
351  std::vector<const G4ThreeVectorList *> polygons;
352  polygons.resize(ksteps+2);
353  for (G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
354  G4BoundingEnvelope benv(bmin,bmax,polygons);
355  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
356  }
357  return exist;
358 }
359 
361 //
362 // Return whether point inside/outside/on surface
363 
365 {
366  G4double r2,pPhi,tolRMin,tolRMax;
367  EInside in = kOutside ;
368 
369  if (std::fabs(p.z()) <= fDz - halfCarTolerance)
370  {
371  r2 = p.x()*p.x() + p.y()*p.y() ;
372 
373  if (fRMin) { tolRMin = fRMin + halfRadTolerance ; }
374  else { tolRMin = 0 ; }
375 
376  tolRMax = fRMax - halfRadTolerance ;
377 
378  if ((r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax))
379  {
380  if ( fPhiFullTube )
381  {
382  in = kInside ;
383  }
384  else
385  {
386  // Try inner tolerant phi boundaries (=>inside)
387  // if not inside, try outer tolerant phi boundaries
388 
389  if ( (tolRMin==0) && (std::fabs(p.x())<=halfCarTolerance)
390  && (std::fabs(p.y())<=halfCarTolerance) )
391  {
392  in=kSurface;
393  }
394  else
395  {
396  pPhi = std::atan2(p.y(),p.x()) ;
397  if ( pPhi < -halfAngTolerance ) { pPhi += twopi; } // 0<=pPhi<2pi
398 
399  if ( fSPhi >= 0 )
400  {
401  if ( (std::fabs(pPhi) < halfAngTolerance)
402  && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
403  {
404  pPhi += twopi ; // 0 <= pPhi < 2pi
405  }
406  if ( (pPhi >= fSPhi + halfAngTolerance)
407  && (pPhi <= fSPhi + fDPhi - halfAngTolerance) )
408  {
409  in = kInside ;
410  }
411  else if ( (pPhi >= fSPhi - halfAngTolerance)
412  && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
413  {
414  in = kSurface ;
415  }
416  }
417  else // fSPhi < 0
418  {
419  if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
420  && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;} //kOutside
421  else if ( (pPhi <= fSPhi + twopi + halfAngTolerance)
422  && (pPhi >= fSPhi + fDPhi - halfAngTolerance) )
423  {
424  in = kSurface ;
425  }
426  else
427  {
428  in = kInside ;
429  }
430  }
431  }
432  }
433  }
434  else // Try generous boundaries
435  {
436  tolRMin = fRMin - halfRadTolerance ;
437  tolRMax = fRMax + halfRadTolerance ;
438 
439  if ( tolRMin < 0 ) { tolRMin = 0; }
440 
441  if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
442  {
443  if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance) )
444  { // Continuous in phi or on z-axis
445  in = kSurface ;
446  }
447  else // Try outer tolerant phi boundaries only
448  {
449  pPhi = std::atan2(p.y(),p.x()) ;
450 
451  if ( pPhi < -halfAngTolerance) { pPhi += twopi; } // 0<=pPhi<2pi
452  if ( fSPhi >= 0 )
453  {
454  if ( (std::fabs(pPhi) < halfAngTolerance)
455  && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
456  {
457  pPhi += twopi ; // 0 <= pPhi < 2pi
458  }
459  if ( (pPhi >= fSPhi - halfAngTolerance)
460  && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
461  {
462  in = kSurface ;
463  }
464  }
465  else // fSPhi < 0
466  {
467  if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
468  && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;} // kOutside
469  else
470  {
471  in = kSurface ;
472  }
473  }
474  }
475  }
476  }
477  }
478  else if (std::fabs(p.z()) <= fDz + halfCarTolerance)
479  { // Check within tolerant r limits
480  r2 = p.x()*p.x() + p.y()*p.y() ;
481  tolRMin = fRMin - halfRadTolerance ;
482  tolRMax = fRMax + halfRadTolerance ;
483 
484  if ( tolRMin < 0 ) { tolRMin = 0; }
485 
486  if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
487  {
488  if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance))
489  { // Continuous in phi or on z-axis
490  in = kSurface ;
491  }
492  else // Try outer tolerant phi boundaries
493  {
494  pPhi = std::atan2(p.y(),p.x()) ;
495 
496  if ( pPhi < -halfAngTolerance ) { pPhi += twopi; } // 0<=pPhi<2pi
497  if ( fSPhi >= 0 )
498  {
499  if ( (std::fabs(pPhi) < halfAngTolerance)
500  && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
501  {
502  pPhi += twopi ; // 0 <= pPhi < 2pi
503  }
504  if ( (pPhi >= fSPhi - halfAngTolerance)
505  && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
506  {
507  in = kSurface;
508  }
509  }
510  else // fSPhi < 0
511  {
512  if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
513  && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;}
514  else
515  {
516  in = kSurface ;
517  }
518  }
519  }
520  }
521  }
522  return in;
523 }
524 
526 //
527 // Return unit normal of surface closest to p
528 // - note if point on z axis, ignore phi divided sides
529 // - unsafe if point close to z axis a rmin=0 - no explicit checks
530 
532 {
533  G4int noSurfaces = 0;
534  G4double rho, pPhi;
535  G4double distZ, distRMin, distRMax;
536  G4double distSPhi = kInfinity, distEPhi = kInfinity;
537 
538  G4ThreeVector norm, sumnorm(0.,0.,0.);
539  G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
540  G4ThreeVector nR, nPs, nPe;
541 
542  rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
543 
544  distRMin = std::fabs(rho - fRMin);
545  distRMax = std::fabs(rho - fRMax);
546  distZ = std::fabs(std::fabs(p.z()) - fDz);
547 
548  if (!fPhiFullTube) // Protected against (0,0,z)
549  {
550  if ( rho > halfCarTolerance )
551  {
552  pPhi = std::atan2(p.y(),p.x());
553 
554  if(pPhi < fSPhi- halfCarTolerance) { pPhi += twopi; }
555  else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; }
556 
557  distSPhi = std::fabs(pPhi - fSPhi);
558  distEPhi = std::fabs(pPhi - fSPhi - fDPhi);
559  }
560  else if( !fRMin )
561  {
562  distSPhi = 0.;
563  distEPhi = 0.;
564  }
565  nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
566  nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
567  }
568  if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
569 
570  if( distRMax <= halfCarTolerance )
571  {
572  noSurfaces ++;
573  sumnorm += nR;
574  }
575  if( fRMin && (distRMin <= halfCarTolerance) )
576  {
577  noSurfaces ++;
578  sumnorm -= nR;
579  }
580  if( fDPhi < twopi )
581  {
582  if (distSPhi <= halfAngTolerance)
583  {
584  noSurfaces ++;
585  sumnorm += nPs;
586  }
587  if (distEPhi <= halfAngTolerance)
588  {
589  noSurfaces ++;
590  sumnorm += nPe;
591  }
592  }
593  if (distZ <= halfCarTolerance)
594  {
595  noSurfaces ++;
596  if ( p.z() >= 0.) { sumnorm += nZ; }
597  else { sumnorm -= nZ; }
598  }
599  if ( noSurfaces == 0 )
600  {
601 #ifdef G4CSGDEBUG
602  G4Exception("G4Tubs::SurfaceNormal(p)", "GeomSolids1002",
603  JustWarning, "Point p is not on surface !?" );
604  G4int oldprc = G4cout.precision(20);
605  G4cout<< "G4Tubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); "
606  << G4endl << G4endl;
607  G4cout.precision(oldprc) ;
608 #endif
609  norm = ApproxSurfaceNormal(p);
610  }
611  else if ( noSurfaces == 1 ) { norm = sumnorm; }
612  else { norm = sumnorm.unit(); }
613 
614  return norm;
615 }
616 
618 //
619 // Algorithm for SurfaceNormal() following the original specification
620 // for points not on the surface
621 
623 {
624  ENorm side ;
625  G4ThreeVector norm ;
626  G4double rho, phi ;
627  G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
628 
629  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
630 
631  distRMin = std::fabs(rho - fRMin) ;
632  distRMax = std::fabs(rho - fRMax) ;
633  distZ = std::fabs(std::fabs(p.z()) - fDz) ;
634 
635  if (distRMin < distRMax) // First minimum
636  {
637  if ( distZ < distRMin )
638  {
639  distMin = distZ ;
640  side = kNZ ;
641  }
642  else
643  {
644  distMin = distRMin ;
645  side = kNRMin ;
646  }
647  }
648  else
649  {
650  if ( distZ < distRMax )
651  {
652  distMin = distZ ;
653  side = kNZ ;
654  }
655  else
656  {
657  distMin = distRMax ;
658  side = kNRMax ;
659  }
660  }
661  if (!fPhiFullTube && rho ) // Protected against (0,0,z)
662  {
663  phi = std::atan2(p.y(),p.x()) ;
664 
665  if ( phi < 0 ) { phi += twopi; }
666 
667  if ( fSPhi < 0 )
668  {
669  distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;
670  }
671  else
672  {
673  distSPhi = std::fabs(phi - fSPhi)*rho ;
674  }
675  distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
676 
677  if (distSPhi < distEPhi) // Find new minimum
678  {
679  if ( distSPhi < distMin )
680  {
681  side = kNSPhi ;
682  }
683  }
684  else
685  {
686  if ( distEPhi < distMin )
687  {
688  side = kNEPhi ;
689  }
690  }
691  }
692  switch ( side )
693  {
694  case kNRMin : // Inner radius
695  {
696  norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ;
697  break ;
698  }
699  case kNRMax : // Outer radius
700  {
701  norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ;
702  break ;
703  }
704  case kNZ : // + or - dz
705  {
706  if ( p.z() > 0 ) { norm = G4ThreeVector(0,0,1) ; }
707  else { norm = G4ThreeVector(0,0,-1); }
708  break ;
709  }
710  case kNSPhi:
711  {
712  norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
713  break ;
714  }
715  case kNEPhi:
716  {
717  norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
718  break;
719  }
720  default: // Should never reach this case ...
721  {
722  DumpInfo();
723  G4Exception("G4Tubs::ApproxSurfaceNormal()",
724  "GeomSolids1002", JustWarning,
725  "Undefined side for valid surface normal to solid.");
726  break ;
727  }
728  }
729  return norm;
730 }
731 
733 //
734 //
735 // Calculate distance to shape from outside, along normalised vector
736 // - return kInfinity if no intersection, or intersection distance <= tolerance
737 //
738 // - Compute the intersection with the z planes
739 // - if at valid r, phi, return
740 //
741 // -> If point is outer outer radius, compute intersection with rmax
742 // - if at valid phi,z return
743 //
744 // -> Compute intersection with inner radius, taking largest +ve root
745 // - if valid (in z,phi), save intersction
746 //
747 // -> If phi segmented, compute intersections with phi half planes
748 // - return smallest of valid phi intersections and
749 // inner radius intersection
750 //
751 // NOTE:
752 // - 'if valid' implies tolerant checking of intersection points
753 
755  const G4ThreeVector& v ) const
756 {
757  G4double snxt = kInfinity ; // snxt = default return value
758  G4double tolORMin2, tolIRMax2 ; // 'generous' radii squared
759  G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
760  const G4double dRmax = 100.*fRMax;
761 
762  // Intersection point variables
763  //
764  G4double Dist, sd, xi, yi, zi, rho2, inum, iden, cosPsi, Comp ;
765  G4double t1, t2, t3, b, c, d ; // Quadratic solver variables
766 
767  // Calculate tolerant rmin and rmax
768 
769  if (fRMin > kRadTolerance)
770  {
771  tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ;
772  tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ;
773  }
774  else
775  {
776  tolORMin2 = 0.0 ;
777  tolIRMin2 = 0.0 ;
778  }
779  tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ;
780  tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - halfRadTolerance) ;
781 
782  // Intersection with Z surfaces
783 
784  tolIDz = fDz - halfCarTolerance ;
785  tolODz = fDz + halfCarTolerance ;
786 
787  if (std::fabs(p.z()) >= tolIDz)
788  {
789  if ( p.z()*v.z() < 0 ) // at +Z going in -Z or visa versa
790  {
791  sd = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ; // Z intersect distance
792 
793  if(sd < 0.0) { sd = 0.0; }
794 
795  xi = p.x() + sd*v.x() ; // Intersection coords
796  yi = p.y() + sd*v.y() ;
797  rho2 = xi*xi + yi*yi ;
798 
799  // Check validity of intersection
800 
801  if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
802  {
803  if (!fPhiFullTube && rho2)
804  {
805  // Psi = angle made with central (average) phi of shape
806  //
807  inum = xi*cosCPhi + yi*sinCPhi ;
808  iden = std::sqrt(rho2) ;
809  cosPsi = inum/iden ;
810  if (cosPsi >= cosHDPhiIT) { return sd ; }
811  }
812  else
813  {
814  return sd ;
815  }
816  }
817  }
818  else
819  {
820  if ( snxt<halfCarTolerance ) { snxt=0; }
821  return snxt ; // On/outside extent, and heading away
822  // -> cannot intersect
823  }
824  }
825 
826  // -> Can not intersect z surfaces
827  //
828  // Intersection with rmax (possible return) and rmin (must also check phi)
829  //
830  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
831  //
832  // Intersects with x^2+y^2=R^2
833  //
834  // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
835  // t1 t2 t3
836 
837  t1 = 1.0 - v.z()*v.z() ;
838  t2 = p.x()*v.x() + p.y()*v.y() ;
839  t3 = p.x()*p.x() + p.y()*p.y() ;
840 
841  if ( t1 > 0 ) // Check not || to z axis
842  {
843  b = t2/t1 ;
844  c = t3 - fRMax*fRMax ;
845  if ((t3 >= tolORMax2) && (t2<0)) // This also handles the tangent case
846  {
847  // Try outer cylinder intersection
848  // c=(t3-fRMax*fRMax)/t1;
849 
850  c /= t1 ;
851  d = b*b - c ;
852 
853  if (d >= 0) // If real root
854  {
855  sd = c/(-b+std::sqrt(d));
856  if (sd >= 0) // If 'forwards'
857  {
858  if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
859  { // 64 bits systems. Split long distances and recompute
860  G4double fTerm = sd-std::fmod(sd,dRmax);
861  sd = fTerm + DistanceToIn(p+fTerm*v,v);
862  }
863  // Check z intersection
864  //
865  zi = p.z() + sd*v.z() ;
866  if (std::fabs(zi)<=tolODz)
867  {
868  // Z ok. Check phi intersection if reqd
869  //
870  if (fPhiFullTube)
871  {
872  return sd ;
873  }
874  else
875  {
876  xi = p.x() + sd*v.x() ;
877  yi = p.y() + sd*v.y() ;
878  cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ;
879  if (cosPsi >= cosHDPhiIT) { return sd ; }
880  }
881  } // end if std::fabs(zi)
882  } // end if (sd>=0)
883  } // end if (d>=0)
884  } // end if (r>=fRMax)
885  else
886  {
887  // Inside outer radius :
888  // check not inside, and heading through tubs (-> 0 to in)
889 
890  if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.z()) <= tolIDz))
891  {
892  // Inside both radii, delta r -ve, inside z extent
893 
894  if (!fPhiFullTube)
895  {
896  inum = p.x()*cosCPhi + p.y()*sinCPhi ;
897  iden = std::sqrt(t3) ;
898  cosPsi = inum/iden ;
899  if (cosPsi >= cosHDPhiIT)
900  {
901  // In the old version, the small negative tangent for the point
902  // on surface was not taken in account, and returning 0.0 ...
903  // New version: check the tangent for the point on surface and
904  // if no intersection, return kInfinity, if intersection instead
905  // return sd.
906  //
907  c = t3-fRMax*fRMax;
908  if ( c<=0.0 )
909  {
910  return 0.0;
911  }
912  else
913  {
914  c = c/t1 ;
915  d = b*b-c;
916  if ( d>=0.0 )
917  {
918  snxt = c/(-b+std::sqrt(d)); // using safe solution
919  // for quadratic equation
920  if ( snxt < halfCarTolerance ) { snxt=0; }
921  return snxt ;
922  }
923  else
924  {
925  return kInfinity;
926  }
927  }
928  }
929  }
930  else
931  {
932  // In the old version, the small negative tangent for the point
933  // on surface was not taken in account, and returning 0.0 ...
934  // New version: check the tangent for the point on surface and
935  // if no intersection, return kInfinity, if intersection instead
936  // return sd.
937  //
938  c = t3 - fRMax*fRMax;
939  if ( c<=0.0 )
940  {
941  return 0.0;
942  }
943  else
944  {
945  c = c/t1 ;
946  d = b*b-c;
947  if ( d>=0.0 )
948  {
949  snxt= c/(-b+std::sqrt(d)); // using safe solution
950  // for quadratic equation
951  if ( snxt < halfCarTolerance ) { snxt=0; }
952  return snxt ;
953  }
954  else
955  {
956  return kInfinity;
957  }
958  }
959  } // end if (!fPhiFullTube)
960  } // end if (t3>tolIRMin2)
961  } // end if (Inside Outer Radius)
962  if ( fRMin ) // Try inner cylinder intersection
963  {
964  c = (t3 - fRMin*fRMin)/t1 ;
965  d = b*b - c ;
966  if ( d >= 0.0 ) // If real root
967  {
968  // Always want 2nd root - we are outside and know rmax Hit was bad
969  // - If on surface of rmin also need farthest root
970 
971  sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
972  if (sd >= -halfCarTolerance) // check forwards
973  {
974  // Check z intersection
975  //
976  if(sd < 0.0) { sd = 0.0; }
977  if ( sd>dRmax ) // Avoid rounding errors due to precision issues seen
978  { // 64 bits systems. Split long distances and recompute
979  G4double fTerm = sd-std::fmod(sd,dRmax);
980  sd = fTerm + DistanceToIn(p+fTerm*v,v);
981  }
982  zi = p.z() + sd*v.z() ;
983  if (std::fabs(zi) <= tolODz)
984  {
985  // Z ok. Check phi
986  //
987  if ( fPhiFullTube )
988  {
989  return sd ;
990  }
991  else
992  {
993  xi = p.x() + sd*v.x() ;
994  yi = p.y() + sd*v.y() ;
995  cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMin ;
996  if (cosPsi >= cosHDPhiIT)
997  {
998  // Good inner radius isect
999  // - but earlier phi isect still possible
1000 
1001  snxt = sd ;
1002  }
1003  }
1004  } // end if std::fabs(zi)
1005  } // end if (sd>=0)
1006  } // end if (d>=0)
1007  } // end if (fRMin)
1008  }
1009 
1010  // Phi segment intersection
1011  //
1012  // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1013  //
1014  // o NOTE: Large duplication of code between sphi & ephi checks
1015  // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1016  // intersection check <=0 -> >=0
1017  // -> use some form of loop Construct ?
1018  //
1019  if ( !fPhiFullTube )
1020  {
1021  // First phi surface (Starting phi)
1022  //
1023  Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1024 
1025  if ( Comp < 0 ) // Component in outwards normal dirn
1026  {
1027  Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1028 
1029  if ( Dist < halfCarTolerance )
1030  {
1031  sd = Dist/Comp ;
1032 
1033  if (sd < snxt)
1034  {
1035  if ( sd < 0 ) { sd = 0.0; }
1036  zi = p.z() + sd*v.z() ;
1037  if ( std::fabs(zi) <= tolODz )
1038  {
1039  xi = p.x() + sd*v.x() ;
1040  yi = p.y() + sd*v.y() ;
1041  rho2 = xi*xi + yi*yi ;
1042 
1043  if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1044  || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1045  && ( v.y()*cosSPhi - v.x()*sinSPhi > 0 )
1046  && ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 ) )
1047  || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1048  && (v.y()*cosSPhi - v.x()*sinSPhi > 0)
1049  && (v.x()*cosSPhi + v.y()*sinSPhi < 0) ) )
1050  {
1051  // z and r intersections good
1052  // - check intersecting with correct half-plane
1053  //
1054  if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = sd; }
1055  }
1056  }
1057  }
1058  }
1059  }
1060 
1061  // Second phi surface (Ending phi)
1062 
1063  Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1064 
1065  if (Comp < 0 ) // Component in outwards normal dirn
1066  {
1067  Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1068 
1069  if ( Dist < halfCarTolerance )
1070  {
1071  sd = Dist/Comp ;
1072 
1073  if (sd < snxt)
1074  {
1075  if ( sd < 0 ) { sd = 0; }
1076  zi = p.z() + sd*v.z() ;
1077  if ( std::fabs(zi) <= tolODz )
1078  {
1079  xi = p.x() + sd*v.x() ;
1080  yi = p.y() + sd*v.y() ;
1081  rho2 = xi*xi + yi*yi ;
1082  if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1083  || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1084  && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1085  && (v.x()*cosEPhi + v.y()*sinEPhi >= 0) )
1086  || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1087  && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1088  && (v.x()*cosEPhi + v.y()*sinEPhi < 0) ) )
1089  {
1090  // z and r intersections good
1091  // - check intersecting with correct half-plane
1092  //
1093  if ( (yi*cosCPhi-xi*sinCPhi) >= 0 ) { snxt = sd; }
1094  } //?? >=-halfCarTolerance
1095  }
1096  }
1097  }
1098  } // Comp < 0
1099  } // !fPhiFullTube
1100  if ( snxt<halfCarTolerance ) { snxt=0; }
1101  return snxt ;
1102 }
1103 
1105 //
1106 // Calculate distance to shape from outside, along normalised vector
1107 // - return kInfinity if no intersection, or intersection distance <= tolerance
1108 //
1109 // - Compute the intersection with the z planes
1110 // - if at valid r, phi, return
1111 //
1112 // -> If point is outer outer radius, compute intersection with rmax
1113 // - if at valid phi,z return
1114 //
1115 // -> Compute intersection with inner radius, taking largest +ve root
1116 // - if valid (in z,phi), save intersction
1117 //
1118 // -> If phi segmented, compute intersections with phi half planes
1119 // - return smallest of valid phi intersections and
1120 // inner radius intersection
1121 //
1122 // NOTE:
1123 // - Precalculations for phi trigonometry are Done `just in time'
1124 // - `if valid' implies tolerant checking of intersection points
1125 // Calculate distance (<= actual) to closest surface of shape from outside
1126 // - Calculate distance to z, radial planes
1127 // - Only to phi planes if outside phi extent
1128 // - Return 0 if point inside
1129 
1131 {
1132  G4double safe=0.0, rho, safe1, safe2, safe3 ;
1133  G4double safePhi, cosPsi ;
1134 
1135  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1136  safe1 = fRMin - rho ;
1137  safe2 = rho - fRMax ;
1138  safe3 = std::fabs(p.z()) - fDz ;
1139 
1140  if ( safe1 > safe2 ) { safe = safe1; }
1141  else { safe = safe2; }
1142  if ( safe3 > safe ) { safe = safe3; }
1143 
1144  if ( (!fPhiFullTube) && (rho) )
1145  {
1146  // Psi=angle from central phi to point
1147  //
1148  cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1149 
1150  if ( cosPsi < std::cos(fDPhi*0.5) )
1151  {
1152  // Point lies outside phi range
1153 
1154  if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
1155  {
1156  safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ;
1157  }
1158  else
1159  {
1160  safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
1161  }
1162  if ( safePhi > safe ) { safe = safePhi; }
1163  }
1164  }
1165  if ( safe < 0 ) { safe = 0; }
1166  return safe ;
1167 }
1168 
1170 //
1171 // Calculate distance to surface of shape from `inside', allowing for tolerance
1172 // - Only Calc rmax intersection if no valid rmin intersection
1173 
1175  const G4ThreeVector& v,
1176  const G4bool calcNorm,
1177  G4bool *validNorm,
1178  G4ThreeVector *n ) const
1179 {
1180  ESide side=kNull , sider=kNull, sidephi=kNull ;
1181  G4double snxt, srd=kInfinity, sphi=kInfinity, pdist ;
1182  G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
1183 
1184  // Vars for phi intersection:
1185 
1186  G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1187 
1188  // Z plane intersection
1189 
1190  if (v.z() > 0 )
1191  {
1192  pdist = fDz - p.z() ;
1193  if ( pdist > halfCarTolerance )
1194  {
1195  snxt = pdist/v.z() ;
1196  side = kPZ ;
1197  }
1198  else
1199  {
1200  if (calcNorm)
1201  {
1202  *n = G4ThreeVector(0,0,1) ;
1203  *validNorm = true ;
1204  }
1205  return snxt = 0 ;
1206  }
1207  }
1208  else if ( v.z() < 0 )
1209  {
1210  pdist = fDz + p.z() ;
1211 
1212  if ( pdist > halfCarTolerance )
1213  {
1214  snxt = -pdist/v.z() ;
1215  side = kMZ ;
1216  }
1217  else
1218  {
1219  if (calcNorm)
1220  {
1221  *n = G4ThreeVector(0,0,-1) ;
1222  *validNorm = true ;
1223  }
1224  return snxt = 0.0 ;
1225  }
1226  }
1227  else
1228  {
1229  snxt = kInfinity ; // Travel perpendicular to z axis
1230  side = kNull;
1231  }
1232 
1233  // Radial Intersections
1234  //
1235  // Find intersection with cylinders at rmax/rmin
1236  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1237  //
1238  // Intersects with x^2+y^2=R^2
1239  //
1240  // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
1241  //
1242  // t1 t2 t3
1243 
1244  t1 = 1.0 - v.z()*v.z() ; // since v normalised
1245  t2 = p.x()*v.x() + p.y()*v.y() ;
1246  t3 = p.x()*p.x() + p.y()*p.y() ;
1247 
1248  if ( snxt > 10*(fDz+fRMax) ) { roi2 = 2*fRMax*fRMax; }
1249  else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; } // radius^2 on +-fDz
1250 
1251  if ( t1 > 0 ) // Check not parallel
1252  {
1253  // Calculate srd, r exit distance
1254 
1255  if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) )
1256  {
1257  // Delta r not negative => leaving via rmax
1258 
1259  deltaR = t3 - fRMax*fRMax ;
1260 
1261  // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
1262  // - avoid sqrt for efficiency
1263 
1264  if ( deltaR < -kRadTolerance*fRMax )
1265  {
1266  b = t2/t1 ;
1267  c = deltaR/t1 ;
1268  d2 = b*b-c;
1269  if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1270  else { srd = 0.; }
1271  sider = kRMax ;
1272  }
1273  else
1274  {
1275  // On tolerant boundary & heading outwards (or perpendicular to)
1276  // outer radial surface -> leaving immediately
1277 
1278  if ( calcNorm )
1279  {
1280  *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1281  *validNorm = true ;
1282  }
1283  return snxt = 0 ; // Leaving by rmax immediately
1284  }
1285  }
1286  else if ( t2 < 0. ) // i.e. t2 < 0; Possible rmin intersection
1287  {
1288  roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement
1289 
1290  if ( fRMin && (roMin2 < fRMin*(fRMin - kRadTolerance)) )
1291  {
1292  deltaR = t3 - fRMin*fRMin ;
1293  b = t2/t1 ;
1294  c = deltaR/t1 ;
1295  d2 = b*b - c ;
1296 
1297  if ( d2 >= 0 ) // Leaving via rmin
1298  {
1299  // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1300  // - avoid sqrt for efficiency
1301 
1302  if (deltaR > kRadTolerance*fRMin)
1303  {
1304  srd = c/(-b+std::sqrt(d2));
1305  sider = kRMin ;
1306  }
1307  else
1308  {
1309  if ( calcNorm ) { *validNorm = false; } // Concave side
1310  return snxt = 0.0;
1311  }
1312  }
1313  else // No rmin intersect -> must be rmax intersect
1314  {
1315  deltaR = t3 - fRMax*fRMax ;
1316  c = deltaR/t1 ;
1317  d2 = b*b-c;
1318  if( d2 >=0. )
1319  {
1320  srd = -b + std::sqrt(d2) ;
1321  sider = kRMax ;
1322  }
1323  else // Case: On the border+t2<kRadTolerance
1324  // (v is perpendicular to the surface)
1325  {
1326  if (calcNorm)
1327  {
1328  *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1329  *validNorm = true ;
1330  }
1331  return snxt = 0.0;
1332  }
1333  }
1334  }
1335  else if ( roi2 > fRMax*(fRMax + kRadTolerance) )
1336  // No rmin intersect -> must be rmax intersect
1337  {
1338  deltaR = t3 - fRMax*fRMax ;
1339  b = t2/t1 ;
1340  c = deltaR/t1;
1341  d2 = b*b-c;
1342  if( d2 >= 0 )
1343  {
1344  srd = -b + std::sqrt(d2) ;
1345  sider = kRMax ;
1346  }
1347  else // Case: On the border+t2<kRadTolerance
1348  // (v is perpendicular to the surface)
1349  {
1350  if (calcNorm)
1351  {
1352  *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1353  *validNorm = true ;
1354  }
1355  return snxt = 0.0;
1356  }
1357  }
1358  }
1359 
1360  // Phi Intersection
1361 
1362  if ( !fPhiFullTube )
1363  {
1364  // add angle calculation with correction
1365  // of the difference in domain of atan2 and Sphi
1366  //
1367  vphi = std::atan2(v.y(),v.x()) ;
1368 
1369  if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1370  else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1371 
1372 
1373  if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
1374  {
1375  // pDist -ve when inside
1376 
1377  pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1378  pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1379 
1380  // Comp -ve when in direction of outwards normal
1381 
1382  compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1383  compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1384 
1385  sidephi = kNull;
1386 
1387  if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1388  && (pDistE <= halfCarTolerance) ) )
1389  || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
1390  && (pDistE > halfCarTolerance) ) ) )
1391  {
1392  // Inside both phi *full* planes
1393 
1394  if ( compS < 0 )
1395  {
1396  sphi = pDistS/compS ;
1397 
1398  if (sphi >= -halfCarTolerance)
1399  {
1400  xi = p.x() + sphi*v.x() ;
1401  yi = p.y() + sphi*v.y() ;
1402 
1403  // Check intersecting with correct half-plane
1404  // (if not -> no intersect)
1405  //
1406  if((std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance))
1407  {
1408  sidephi = kSPhi;
1409  if (((fSPhi-halfAngTolerance)<=vphi)
1410  &&((fSPhi+fDPhi+halfAngTolerance)>=vphi))
1411  {
1412  sphi = kInfinity;
1413  }
1414  }
1415  else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1416  {
1417  sphi = kInfinity ;
1418  }
1419  else
1420  {
1421  sidephi = kSPhi ;
1422  if ( pDistS > -halfCarTolerance )
1423  {
1424  sphi = 0.0 ; // Leave by sphi immediately
1425  }
1426  }
1427  }
1428  else
1429  {
1430  sphi = kInfinity ;
1431  }
1432  }
1433  else
1434  {
1435  sphi = kInfinity ;
1436  }
1437 
1438  if ( compE < 0 )
1439  {
1440  sphi2 = pDistE/compE ;
1441 
1442  // Only check further if < starting phi intersection
1443  //
1444  if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1445  {
1446  xi = p.x() + sphi2*v.x() ;
1447  yi = p.y() + sphi2*v.y() ;
1448 
1449  if((std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance))
1450  {
1451  // Leaving via ending phi
1452  //
1453  if( !((fSPhi-halfAngTolerance <= vphi)
1454  &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
1455  {
1456  sidephi = kEPhi ;
1457  if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1458  else { sphi = 0.0 ; }
1459  }
1460  }
1461  else // Check intersecting with correct half-plane
1462 
1463  if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1464  {
1465  // Leaving via ending phi
1466  //
1467  sidephi = kEPhi ;
1468  if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1469  else { sphi = 0.0 ; }
1470  }
1471  }
1472  }
1473  }
1474  else
1475  {
1476  sphi = kInfinity ;
1477  }
1478  }
1479  else
1480  {
1481  // On z axis + travel not || to z axis -> if phi of vector direction
1482  // within phi of shape, Step limited by rmax, else Step =0
1483 
1484  if ( (fSPhi - halfAngTolerance <= vphi)
1485  && (vphi <= fSPhi + fDPhi + halfAngTolerance ) )
1486  {
1487  sphi = kInfinity ;
1488  }
1489  else
1490  {
1491  sidephi = kSPhi ; // arbitrary
1492  sphi = 0.0 ;
1493  }
1494  }
1495  if (sphi < snxt) // Order intersecttions
1496  {
1497  snxt = sphi ;
1498  side = sidephi ;
1499  }
1500  }
1501  if (srd < snxt) // Order intersections
1502  {
1503  snxt = srd ;
1504  side = sider ;
1505  }
1506  }
1507  if (calcNorm)
1508  {
1509  switch(side)
1510  {
1511  case kRMax:
1512  // Note: returned vector not normalised
1513  // (divide by fRMax for unit vector)
1514  //
1515  xi = p.x() + snxt*v.x() ;
1516  yi = p.y() + snxt*v.y() ;
1517  *n = G4ThreeVector(xi/fRMax,yi/fRMax,0) ;
1518  *validNorm = true ;
1519  break ;
1520 
1521  case kRMin:
1522  *validNorm = false ; // Rmin is inconvex
1523  break ;
1524 
1525  case kSPhi:
1526  if ( fDPhi <= pi )
1527  {
1528  *n = G4ThreeVector(sinSPhi,-cosSPhi,0) ;
1529  *validNorm = true ;
1530  }
1531  else
1532  {
1533  *validNorm = false ;
1534  }
1535  break ;
1536 
1537  case kEPhi:
1538  if (fDPhi <= pi)
1539  {
1540  *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ;
1541  *validNorm = true ;
1542  }
1543  else
1544  {
1545  *validNorm = false ;
1546  }
1547  break ;
1548 
1549  case kPZ:
1550  *n = G4ThreeVector(0,0,1) ;
1551  *validNorm = true ;
1552  break ;
1553 
1554  case kMZ:
1555  *n = G4ThreeVector(0,0,-1) ;
1556  *validNorm = true ;
1557  break ;
1558 
1559  default:
1560  G4cout << G4endl ;
1561  DumpInfo();
1562  std::ostringstream message;
1563  G4int oldprc = message.precision(16);
1564  message << "Undefined side for valid surface normal to solid."
1565  << G4endl
1566  << "Position:" << G4endl << G4endl
1567  << "p.x() = " << p.x()/mm << " mm" << G4endl
1568  << "p.y() = " << p.y()/mm << " mm" << G4endl
1569  << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
1570  << "Direction:" << G4endl << G4endl
1571  << "v.x() = " << v.x() << G4endl
1572  << "v.y() = " << v.y() << G4endl
1573  << "v.z() = " << v.z() << G4endl << G4endl
1574  << "Proposed distance :" << G4endl << G4endl
1575  << "snxt = " << snxt/mm << " mm" << G4endl ;
1576  message.precision(oldprc) ;
1577  G4Exception("G4Tubs::DistanceToOut(p,v,..)", "GeomSolids1002",
1578  JustWarning, message);
1579  break ;
1580  }
1581  }
1582  if ( snxt<halfCarTolerance ) { snxt=0 ; }
1583 
1584  return snxt ;
1585 }
1586 
1588 //
1589 // Calculate distance (<=actual) to closest surface of shape from inside
1590 
1592 {
1593  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi ;
1594  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1595 
1596 #ifdef G4CSGDEBUG
1597  if( Inside(p) == kOutside )
1598  {
1599  G4int oldprc = G4cout.precision(16) ;
1600  G4cout << G4endl ;
1601  DumpInfo();
1602  G4cout << "Position:" << G4endl << G4endl ;
1603  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1604  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1605  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1606  G4cout.precision(oldprc) ;
1607  G4Exception("G4Tubs::DistanceToOut(p)", "GeomSolids1002",
1608  JustWarning, "Point p is outside !?");
1609  }
1610 #endif
1611 
1612  if ( fRMin )
1613  {
1614  safeR1 = rho - fRMin ;
1615  safeR2 = fRMax - rho ;
1616 
1617  if ( safeR1 < safeR2 ) { safe = safeR1 ; }
1618  else { safe = safeR2 ; }
1619  }
1620  else
1621  {
1622  safe = fRMax - rho ;
1623  }
1624  safeZ = fDz - std::fabs(p.z()) ;
1625 
1626  if ( safeZ < safe ) { safe = safeZ ; }
1627 
1628  // Check if phi divided, Calc distances closest phi plane
1629  //
1630  if ( !fPhiFullTube )
1631  {
1632  if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
1633  {
1634  safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
1635  }
1636  else
1637  {
1638  safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
1639  }
1640  if (safePhi < safe) { safe = safePhi ; }
1641  }
1642  if ( safe < 0 ) { safe = 0 ; }
1643 
1644  return safe ;
1645 }
1646 
1648 //
1649 // Stream object contents to an output stream
1650 
1652 {
1653  return G4String("G4Tubs");
1654 }
1655 
1657 //
1658 // Make a clone of the object
1659 //
1661 {
1662  return new G4Tubs(*this);
1663 }
1664 
1666 //
1667 // Stream object contents to an output stream
1668 
1669 std::ostream& G4Tubs::StreamInfo( std::ostream& os ) const
1670 {
1671  G4int oldprc = os.precision(16);
1672  os << "-----------------------------------------------------------\n"
1673  << " *** Dump for solid - " << GetName() << " ***\n"
1674  << " ===================================================\n"
1675  << " Solid type: G4Tubs\n"
1676  << " Parameters: \n"
1677  << " inner radius : " << fRMin/mm << " mm \n"
1678  << " outer radius : " << fRMax/mm << " mm \n"
1679  << " half length Z: " << fDz/mm << " mm \n"
1680  << " starting phi : " << fSPhi/degree << " degrees \n"
1681  << " delta phi : " << fDPhi/degree << " degrees \n"
1682  << "-----------------------------------------------------------\n";
1683  os.precision(oldprc);
1684 
1685  return os;
1686 }
1687 
1689 //
1690 // GetPointOnSurface
1691 
1693 {
1694  G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1695  aOne, aTwo, aThr, aFou;
1696  G4double rRand;
1697 
1698  aOne = 2.*fDz*fDPhi*fRMax;
1699  aTwo = 2.*fDz*fDPhi*fRMin;
1700  aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin);
1701  aFou = 2.*fDz*(fRMax-fRMin);
1702 
1704  cosphi = std::cos(phi);
1705  sinphi = std::sin(phi);
1706 
1707  rRand = GetRadiusInRing(fRMin,fRMax);
1708 
1709  if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; }
1710 
1711  chose = G4RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou);
1712 
1713  if( (chose >=0) && (chose < aOne) )
1714  {
1715  xRand = fRMax*cosphi;
1716  yRand = fRMax*sinphi;
1717  zRand = G4RandFlat::shoot(-1.*fDz,fDz);
1718  return G4ThreeVector (xRand, yRand, zRand);
1719  }
1720  else if( (chose >= aOne) && (chose < aOne + aTwo) )
1721  {
1722  xRand = fRMin*cosphi;
1723  yRand = fRMin*sinphi;
1724  zRand = G4RandFlat::shoot(-1.*fDz,fDz);
1725  return G4ThreeVector (xRand, yRand, zRand);
1726  }
1727  else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1728  {
1729  xRand = rRand*cosphi;
1730  yRand = rRand*sinphi;
1731  zRand = fDz;
1732  return G4ThreeVector (xRand, yRand, zRand);
1733  }
1734  else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1735  {
1736  xRand = rRand*cosphi;
1737  yRand = rRand*sinphi;
1738  zRand = -1.*fDz;
1739  return G4ThreeVector (xRand, yRand, zRand);
1740  }
1741  else if( (chose >= aOne + aTwo + 2.*aThr)
1742  && (chose < aOne + aTwo + 2.*aThr + aFou) )
1743  {
1744  xRand = rRand*std::cos(fSPhi);
1745  yRand = rRand*std::sin(fSPhi);
1746  zRand = G4RandFlat::shoot(-1.*fDz,fDz);
1747  return G4ThreeVector (xRand, yRand, zRand);
1748  }
1749  else
1750  {
1751  xRand = rRand*std::cos(fSPhi+fDPhi);
1752  yRand = rRand*std::sin(fSPhi+fDPhi);
1753  zRand = G4RandFlat::shoot(-1.*fDz,fDz);
1754  return G4ThreeVector (xRand, yRand, zRand);
1755  }
1756 }
1757 
1759 //
1760 // Methods for visualisation
1761 
1763 {
1764  scene.AddSolid (*this) ;
1765 }
1766 
1768 {
1769  return new G4PolyhedronTubs (fRMin, fRMax, fDz, fSPhi, fDPhi) ;
1770 }
1771 #endif
G4String GetName() const
ThreeVector shoot(const G4int Ap, const G4int Af)
G4bool fPhiFullTube
Definition: G4Tubs.hh:226
static constexpr double mm
Definition: G4SIunits.hh:115
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
G4double fDPhi
Definition: G4Tubs.hh:217
G4double halfCarTolerance
Definition: G4Tubs.hh:230
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:111
G4Tubs(const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi)
Definition: G4Tubs.cc:89
static const G4double d2
Definition: G4Tubs.hh:85
G4double cosEPhi
Definition: G4Tubs.hh:221
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4Tubs.cc:212
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Tubs.cc:754
G4double GetCosEndPhi() const
G4ThreeVector GetPointOnSurface() const
Definition: G4Tubs.cc:1692
virtual void AddSolid(const G4Box &)=0
virtual G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4Tubs.cc:622
int G4int
Definition: G4Types.hh:78
G4double cosHDPhiOT
Definition: G4Tubs.hh:221
void DumpInfo() const
static constexpr double twopi
Definition: G4SIunits.hh:76
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Tubs.cc:1669
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Tubs.cc:1762
G4double sinSPhi
Definition: G4Tubs.hh:221
G4GeometryType GetEntityType() const
Definition: G4Tubs.cc:1651
G4GLOB_DLL std::ostream G4cout
G4double GetDeltaPhiAngle() const
G4double fRMin
Definition: G4Tubs.hh:217
static constexpr double degree
Definition: G4SIunits.hh:144
G4double GetSinEndPhi() const
G4Polyhedron * CreatePolyhedron() const
Definition: G4Tubs.cc:1767
G4double GetRadialTolerance() const
bool G4bool
Definition: G4Types.hh:79
G4VSolid * Clone() const
Definition: G4Tubs.cc:1660
G4double GetSinStartPhi() const
G4double cosCPhi
Definition: G4Tubs.hh:221
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4double fSPhi
Definition: G4Tubs.hh:217
const G4int n
G4double GetInnerRadius() const
std::vector< G4ThreeVector > G4ThreeVectorList
virtual ~G4Tubs()
Definition: G4Tubs.cc:142
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Tubs.cc:254
G4double sinEPhi
Definition: G4Tubs.hh:221
G4double halfAngTolerance
Definition: G4Tubs.hh:230
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Tubs.cc:1174
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
EInside Inside(const G4ThreeVector &p) const
Definition: G4Tubs.cc:364
G4double halfRadTolerance
Definition: G4Tubs.hh:230
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
G4double GetZHalfLength() const
G4double fRMax
Definition: G4Tubs.hh:217
CLHEP::Hep2Vector G4TwoVector
Definition: G4TwoVector.hh:42
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Tubs.cc:531
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4double kRadTolerance
Definition: G4Tubs.hh:213
static constexpr double pi
Definition: G4SIunits.hh:75
G4double fDz
Definition: G4Tubs.hh:217
G4double cosSPhi
Definition: G4Tubs.hh:221
double G4double
Definition: G4Types.hh:76
static constexpr double deg
Definition: G4SIunits.hh:152
void CheckPhiAngles(G4double sPhi, G4double dPhi)
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:91
G4double cosHDPhiIT
Definition: G4Tubs.hh:221
G4double GetAngularTolerance() const
G4double sinCPhi
Definition: G4Tubs.hh:221
G4double GetCosStartPhi() const
static G4bool DiskExtent(G4double rmin, G4double rmax, G4double startPhi, G4double delPhi, G4TwoVector &pmin, G4TwoVector &pmax)
Definition: G4GeomTools.cc:378
static G4GeometryTolerance * GetInstance()
G4double GetOuterRadius() const
G4double kAngTolerance
Definition: G4Tubs.hh:213
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Tubs.cc:201
G4Tubs & operator=(const G4Tubs &rhs)
Definition: G4Tubs.cc:169