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