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