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