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