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