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