Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/CSG/src/G4CutTubs.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // G4CutTubs implementation
 27 //
 28 // 01.06.11 T.Nikitina - Derived from G4Tubs
 29 // 30.10.16 E.Tcherniaev - reimplemented CalculateExtent(),
 30 //                         removed CreateRotatedVetices()
 31 // --------------------------------------------------------------------
 32 
 33 #include "G4CutTubs.hh"
 34 
 35 #if !defined(G4GEOM_USE_UCTUBS)
 36 
 37 #include "G4GeomTools.hh"
 38 #include "G4VoxelLimits.hh"
 39 #include "G4AffineTransform.hh"
 40 #include "G4GeometryTolerance.hh"
 41 #include "G4BoundingEnvelope.hh"
 42 
 43 #include "G4VPVParameterisation.hh"
 44 #include "G4QuickRand.hh"
 45 
 46 #include "G4VGraphicsScene.hh"
 47 #include "G4Polyhedron.hh"
 48 
 49 #include "G4AutoLock.hh"
 50 
 51 namespace
 52 {
 53   G4Mutex zminmaxMutex = G4MUTEX_INITIALIZER;
 54 }
 55 
 56 using namespace CLHEP;
 57 
 58 /////////////////////////////////////////////////////////////////////////
 59 //
 60 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
 61 //             - note if pdphi>2PI then reset to 2PI
 62 
 63 G4CutTubs::G4CutTubs( const G4String &pName,
 64                       G4double pRMin, G4double pRMax,
 65                       G4double pDz,
 66                       G4double pSPhi, G4double pDPhi,
 67                       G4ThreeVector pLowNorm,G4ThreeVector pHighNorm )
 68   : G4CSGSolid(pName), fRMin(pRMin), fRMax(pRMax), fDz(pDz),
 69     fSPhi(0.), fDPhi(0.), fZMin(0.), fZMax(0.)
 70 {
 71   kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
 72   kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
 73 
 74   halfCarTolerance = kCarTolerance*0.5;
 75   halfRadTolerance = kRadTolerance*0.5;
 76   halfAngTolerance = kAngTolerance*0.5;
 77 
 78   if (pDz<=0) // Check z-len
 79   {
 80     std::ostringstream message;
 81     message << "Negative Z half-length (" << pDz << ") in solid: " << GetName();
 82     G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002", FatalException, message);
 83   }
 84   if ( (pRMin >= pRMax) || (pRMin < 0) ) // Check radii
 85   {
 86     std::ostringstream message;
 87     message << "Invalid values for radii in solid: " << GetName()
 88             << G4endl
 89             << "        pRMin = " << pRMin << ", pRMax = " << pRMax;
 90     G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002", FatalException, message);
 91   }
 92 
 93   // Check angles
 94   //
 95   CheckPhiAngles(pSPhi, pDPhi);
 96 
 97   // Check on Cutted Planes Normals
 98   // If there is NO CUT, propose to use G4Tubs instead
 99   //
100   if ( ( pLowNorm.x() == 0.0) && ( pLowNorm.y() == 0.0)
101     && ( pHighNorm.x() == 0.0) && (pHighNorm.y() == 0.0) )
102   {
103     std::ostringstream message;
104     message << "Inexisting Low/High Normal to Z plane or Parallel to Z."
105             << G4endl
106             << "Normals to Z plane are " << pLowNorm << " and "
107             << pHighNorm << " in solid: " << GetName() << " \n";
108     G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids1001",
109                 JustWarning, message, "Should use G4Tubs!");
110   }
111 
112   // If Normal is (0,0,0),means parallel to R, give it value of (0,0,+/-1)
113   //
114   if (pLowNorm.mag2() == 0.)  { pLowNorm.setZ(-1.); }
115   if (pHighNorm.mag2()== 0.)  { pHighNorm.setZ(1.); }
116 
117   // Given Normals to Cut Planes have to be an unit vectors.
118   // Normalize if it is needed.
119   //
120   if (pLowNorm.mag2() != 1.)  { pLowNorm  = pLowNorm.unit();  }
121   if (pHighNorm.mag2()!= 1.)  { pHighNorm = pHighNorm.unit(); }
122 
123   // Normals to cutted planes have to point outside Solid
124   //
125   if( (pLowNorm.mag2() != 0.) && (pHighNorm.mag2()!= 0. ) )
126   {
127     if( ( pLowNorm.z()>= 0. ) || ( pHighNorm.z() <= 0.))
128     {
129       std::ostringstream message;
130       message << "Invalid Low or High Normal to Z plane; "
131                  "has to point outside Solid." << G4endl
132               << "Invalid Norm to Z plane (" << pLowNorm << " or  "
133               << pHighNorm << ") in solid: " << GetName();
134       G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002",
135                   FatalException, message);
136     }
137   }
138   fLowNorm  = pLowNorm;
139   fHighNorm = pHighNorm;
140 
141   // Check intersection of cut planes, they MUST NOT intersect
142   // each other inside the lateral surface
143   //
144   if(IsCrossingCutPlanes())
145   {
146     std::ostringstream message;
147     message << "Invalid normals to Z plane in solid : " << GetName() << G4endl
148             << "Cut planes are crossing inside lateral surface !!!\n"
149             << " Solid type: G4CutTubs\n"
150             << " Parameters: \n"
151             << "    inner radius : " << fRMin/mm << " mm \n"
152             << "    outer radius : " << fRMax/mm << " mm \n"
153             << "    half length Z: " << fDz/mm << " mm \n"
154             << "    starting phi : " << fSPhi/degree << " degrees \n"
155             << "    delta phi    : " << fDPhi/degree << " degrees \n"
156             << "    low Norm     : " << fLowNorm << "  \n"
157             << "    high Norm    : " << fHighNorm;
158     G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002",
159                 FatalException, message);
160   }
161 }
162 
163 ///////////////////////////////////////////////////////////////////////
164 //
165 // Fake default constructor - sets only member data and allocates memory
166 //                            for usage restricted to object persistency.
167 //
168 G4CutTubs::G4CutTubs( __void__& a )
169   : G4CSGSolid(a)
170 {
171 }
172 
173 //////////////////////////////////////////////////////////////////////////
174 //
175 // Destructor
176 
177 G4CutTubs::~G4CutTubs() = default;
178 
179 //////////////////////////////////////////////////////////////////////////
180 //
181 // Copy constructor
182 
183 G4CutTubs::G4CutTubs(const G4CutTubs&) = default;
184 
185 //////////////////////////////////////////////////////////////////////////
186 //
187 // Assignment operator
188 
189 G4CutTubs& G4CutTubs::operator = (const G4CutTubs& rhs)
190 {
191    // Check assignment to self
192    //
193    if (this == &rhs)  { return *this; }
194 
195    // Copy base class data
196    //
197    G4CSGSolid::operator=(rhs);
198 
199    // Copy data
200    //
201    kRadTolerance = rhs.kRadTolerance; kAngTolerance = rhs.kAngTolerance;
202    fRMin = rhs.fRMin; fRMax = rhs.fRMax; fDz = rhs.fDz;
203    fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
204    fZMin = rhs.fZMin; fZMax = rhs.fZMax;
205    sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
206    cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
207    sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
208    sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
209    fPhiFullCutTube = rhs.fPhiFullCutTube;
210    halfCarTolerance = rhs.halfCarTolerance;
211    halfRadTolerance = rhs.halfRadTolerance;
212    halfAngTolerance = rhs.halfAngTolerance;
213    fLowNorm = rhs.fLowNorm; fHighNorm = rhs.fHighNorm;
214 
215    return *this;
216 }
217 
218 //////////////////////////////////////////////////////////////////////////
219 //
220 // Get volume
221 
222 G4double G4CutTubs::GetCubicVolume()
223 {
224   constexpr G4int nphi = 200, nrho = 100;
225 
226   if (fCubicVolume == 0.)
227   {
228     // get parameters
229     G4double rmin = GetInnerRadius();
230     G4double rmax = GetOuterRadius();
231     G4double dz   = GetZHalfLength();
232     G4double sphi = GetStartPhiAngle();
233     G4double dphi = GetDeltaPhiAngle();
234 
235     // calculate volume
236     G4double volume = dz*dphi*(rmax*rmax - rmin*rmin);
237     if (dphi < twopi) // make recalculation
238     {
239       // set values for calculation of h - distance between
240       // opposite points on bases
241       G4ThreeVector nbot = GetLowNorm();
242       G4ThreeVector ntop = GetHighNorm();
243       G4double nx = nbot.x()/nbot.z() - ntop.x()/ntop.z();
244       G4double ny = nbot.y()/nbot.z() - ntop.y()/ntop.z();
245 
246       // compute volume by integration
247       G4double delrho = (rmax - rmin)/nrho;
248       G4double delphi = dphi/nphi;
249       volume = 0.;
250       for (G4int irho=0; irho<nrho; ++irho)
251       {
252         G4double r1  = rmin + delrho*irho;
253         G4double r2  = rmin + delrho*(irho + 1);
254         G4double rho = 0.5*(r1 + r2);
255         G4double sector = 0.5*delphi*(r2*r2 - r1*r1);
256         for (G4int iphi=0; iphi<nphi; ++iphi)
257         {
258           G4double phi = sphi + delphi*(iphi + 0.5);
259           G4double h = nx*rho*std::cos(phi) + ny*rho*std::sin(phi) + 2.*dz;
260           volume += sector*h;
261         }
262       }
263     }
264     fCubicVolume = volume;
265   }
266   return fCubicVolume;
267 }
268 
269 //////////////////////////////////////////////////////////////////////////
270 //
271 // Get surface area
272 
273 G4double G4CutTubs::GetSurfaceArea()
274 {
275   constexpr G4int nphi = 400;
276 
277   if (fSurfaceArea == 0.)
278   {
279     // get parameters
280     G4double rmin = GetInnerRadius();
281     G4double rmax = GetOuterRadius();
282     G4double dz   = GetZHalfLength();
283     G4double sphi = GetStartPhiAngle();
284     G4double dphi = GetDeltaPhiAngle();
285     G4ThreeVector nbot = GetLowNorm();
286     G4ThreeVector ntop = GetHighNorm();
287 
288     // calculate lateral surface area
289     G4double sinner = 2.*dz*dphi*rmin;
290     G4double souter = 2.*dz*dphi*rmax;
291     if (dphi < twopi) // make recalculation
292     {
293       // set values for calculation of h - distance between
294       // opposite points on bases
295       G4double nx = nbot.x()/nbot.z() - ntop.x()/ntop.z();
296       G4double ny = nbot.y()/nbot.z() - ntop.y()/ntop.z();
297 
298       // compute lateral surface area by integration
299       G4double delphi = dphi/nphi;
300       sinner = 0.;
301       souter = 0.;
302       for (G4int iphi=0; iphi<nphi; ++iphi)
303       {
304         G4double phi = sphi + delphi*(iphi + 0.5);
305         G4double cosphi = std::cos(phi);
306         G4double sinphi = std::sin(phi);
307         sinner += rmin*(nx*cosphi + ny*sinphi) + 2.*dz;
308         souter += rmax*(nx*cosphi + ny*sinphi) + 2.*dz;
309       }
310       sinner *= delphi*rmin;
311       souter *= delphi*rmax;
312     }
313     // set surface area
314     G4double scut  = (dphi == twopi) ? 0. : 2.*dz*(rmax - rmin);
315     G4double szero = 0.5*dphi*(rmax*rmax - rmin*rmin);
316     G4double slow  = szero/std::abs(nbot.z());
317     G4double shigh = szero/std::abs(ntop.z());
318     fSurfaceArea = sinner + souter + 2.*scut + slow + shigh;
319   }
320   return fSurfaceArea;
321 }
322 
323 //////////////////////////////////////////////////////////////////////////
324 //
325 // Get bounding box
326 
327 void G4CutTubs::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
328 {
329   G4double rmin = GetInnerRadius();
330   G4double rmax = GetOuterRadius();
331   G4double dz   = GetZHalfLength();
332   G4double dphi = GetDeltaPhiAngle();
333 
334   G4double sinSphi = GetSinStartPhi();
335   G4double cosSphi = GetCosStartPhi();
336   G4double sinEphi = GetSinEndPhi();
337   G4double cosEphi = GetCosEndPhi();
338 
339   G4ThreeVector norm;
340   G4double mag, topx, topy, dists, diste;
341   G4bool iftop;
342 
343   // Find Zmin
344   //
345   G4double zmin;
346   norm = GetLowNorm();
347   mag  = std::sqrt(norm.x()*norm.x() + norm.y()*norm.y());
348   topx = (mag == 0) ? 0 : -rmax*norm.x()/mag;
349   topy = (mag == 0) ? 0 : -rmax*norm.y()/mag;
350   dists =  sinSphi*topx - cosSphi*topy;
351   diste = -sinEphi*topx + cosEphi*topy;
352   if (dphi > pi)
353   {
354     iftop = true;
355     if (dists > 0 && diste > 0)iftop = false;
356   }
357   else
358   {
359     iftop = false;
360     if (dists <= 0 && diste <= 0) iftop = true;
361   }
362   if (iftop)
363   {
364     zmin = -(norm.x()*topx + norm.y()*topy)/norm.z() - dz;
365   }
366   else
367   {
368     G4double z1 = -rmin*(norm.x()*cosSphi + norm.y()*sinSphi)/norm.z() - dz;
369     G4double z2 = -rmin*(norm.x()*cosEphi + norm.y()*sinEphi)/norm.z() - dz;
370     G4double z3 = -rmax*(norm.x()*cosSphi + norm.y()*sinSphi)/norm.z() - dz;
371     G4double z4 = -rmax*(norm.x()*cosEphi + norm.y()*sinEphi)/norm.z() - dz;
372     zmin = std::min(std::min(std::min(z1,z2),z3),z4);
373   }
374 
375   // Find Zmax
376   //
377   G4double zmax;
378   norm = GetHighNorm();
379   mag  = std::sqrt(norm.x()*norm.x() + norm.y()*norm.y());
380   topx = (mag == 0) ? 0 : -rmax*norm.x()/mag;
381   topy = (mag == 0) ? 0 : -rmax*norm.y()/mag;
382   dists =  sinSphi*topx - cosSphi*topy;
383   diste = -sinEphi*topx + cosEphi*topy;
384   if (dphi > pi)
385   {
386     iftop = true;
387     if (dists > 0 && diste > 0) iftop = false;
388   }
389   else
390   {
391     iftop = false;
392     if (dists <= 0 && diste <= 0) iftop = true;
393   }
394   if (iftop)
395   {
396     zmax = -(norm.x()*topx + norm.y()*topy)/norm.z() + dz;
397   }
398   else
399   {
400     G4double z1 = -rmin*(norm.x()*cosSphi + norm.y()*sinSphi)/norm.z() + dz;
401     G4double z2 = -rmin*(norm.x()*cosEphi + norm.y()*sinEphi)/norm.z() + dz;
402     G4double z3 = -rmax*(norm.x()*cosSphi + norm.y()*sinSphi)/norm.z() + dz;
403     G4double z4 = -rmax*(norm.x()*cosEphi + norm.y()*sinEphi)/norm.z() + dz;
404     zmax = std::max(std::max(std::max(z1,z2),z3),z4);
405   }
406 
407   // Find bounding box
408   //
409   if (dphi < twopi)
410   {
411     G4TwoVector vmin,vmax;
412     G4GeomTools::DiskExtent(rmin,rmax,
413                             GetSinStartPhi(),GetCosStartPhi(),
414                             GetSinEndPhi(),GetCosEndPhi(),
415                             vmin,vmax);
416     pMin.set(vmin.x(),vmin.y(), zmin);
417     pMax.set(vmax.x(),vmax.y(), zmax);
418   }
419   else
420   {
421     pMin.set(-rmax,-rmax, zmin);
422     pMax.set( rmax, rmax, zmax);
423   }
424 
425   // Check correctness of the bounding box
426   //
427   if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
428   {
429     std::ostringstream message;
430     message << "Bad bounding box (min >= max) for solid: "
431             << GetName() << " !"
432             << "\npMin = " << pMin
433             << "\npMax = " << pMax;
434     G4Exception("G4CutTubs::BoundingLimits()", "GeomMgt0001",
435                 JustWarning, message);
436     DumpInfo();
437   }
438 }
439 
440 //////////////////////////////////////////////////////////////////////////
441 //
442 // Calculate extent under transform and specified limit
443 
444 G4bool G4CutTubs::CalculateExtent( const EAxis              pAxis,
445                                    const G4VoxelLimits&     pVoxelLimit,
446                                    const G4AffineTransform& pTransform,
447                                          G4double&          pMin,
448                                          G4double&          pMax    ) const
449 {
450   G4ThreeVector bmin, bmax;
451   G4bool exist;
452 
453   // Get bounding box
454   BoundingLimits(bmin,bmax);
455 
456   // Check bounding box
457   G4BoundingEnvelope bbox(bmin,bmax);
458 #ifdef G4BBOX_EXTENT
459   return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
460 #endif
461   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
462   {
463     return exist = pMin < pMax;
464   }
465 
466   // Get parameters of the solid
467   G4double rmin = GetInnerRadius();
468   G4double rmax = GetOuterRadius();
469   G4double dphi = GetDeltaPhiAngle();
470   G4double zmin = bmin.z();
471   G4double zmax = bmax.z();
472 
473   // Find bounding envelope and calculate extent
474   //
475   const G4int NSTEPS = 24;            // number of steps for whole circle
476   G4double astep  = twopi/NSTEPS;     // max angle for one step
477   G4int    ksteps = (dphi <= astep) ? 1 : (G4int)((dphi-deg)/astep) + 1;
478   G4double ang    = dphi/ksteps;
479 
480   G4double sinHalf = std::sin(0.5*ang);
481   G4double cosHalf = std::cos(0.5*ang);
482   G4double sinStep = 2.*sinHalf*cosHalf;
483   G4double cosStep = 1. - 2.*sinHalf*sinHalf;
484   G4double rext    = rmax/cosHalf;
485 
486   // bounding envelope for full cylinder consists of two polygons,
487   // in other cases it is a sequence of quadrilaterals
488   if (rmin == 0 && dphi == twopi)
489   {
490     G4double sinCur = sinHalf;
491     G4double cosCur = cosHalf;
492 
493     G4ThreeVectorList baseA(NSTEPS),baseB(NSTEPS);
494     for (G4int k=0; k<NSTEPS; ++k)
495     {
496       baseA[k].set(rext*cosCur,rext*sinCur,zmin);
497       baseB[k].set(rext*cosCur,rext*sinCur,zmax);
498 
499       G4double sinTmp = sinCur;
500       sinCur = sinCur*cosStep + cosCur*sinStep;
501       cosCur = cosCur*cosStep - sinTmp*sinStep;
502     }
503     std::vector<const G4ThreeVectorList *> polygons(2);
504     polygons[0] = &baseA;
505     polygons[1] = &baseB;
506     G4BoundingEnvelope benv(bmin,bmax,polygons);
507     exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
508   }
509   else
510   {
511     G4double sinStart = GetSinStartPhi();
512     G4double cosStart = GetCosStartPhi();
513     G4double sinEnd   = GetSinEndPhi();
514     G4double cosEnd   = GetCosEndPhi();
515     G4double sinCur   = sinStart*cosHalf + cosStart*sinHalf;
516     G4double cosCur   = cosStart*cosHalf - sinStart*sinHalf;
517 
518     // set quadrilaterals
519     G4ThreeVectorList pols[NSTEPS+2];
520     for (G4int k=0; k<ksteps+2; ++k) pols[k].resize(4);
521     pols[0][0].set(rmin*cosStart,rmin*sinStart,zmax);
522     pols[0][1].set(rmin*cosStart,rmin*sinStart,zmin);
523     pols[0][2].set(rmax*cosStart,rmax*sinStart,zmin);
524     pols[0][3].set(rmax*cosStart,rmax*sinStart,zmax);
525     for (G4int k=1; k<ksteps+1; ++k)
526     {
527       pols[k][0].set(rmin*cosCur,rmin*sinCur,zmax);
528       pols[k][1].set(rmin*cosCur,rmin*sinCur,zmin);
529       pols[k][2].set(rext*cosCur,rext*sinCur,zmin);
530       pols[k][3].set(rext*cosCur,rext*sinCur,zmax);
531 
532       G4double sinTmp = sinCur;
533       sinCur = sinCur*cosStep + cosCur*sinStep;
534       cosCur = cosCur*cosStep - sinTmp*sinStep;
535     }
536     pols[ksteps+1][0].set(rmin*cosEnd,rmin*sinEnd,zmax);
537     pols[ksteps+1][1].set(rmin*cosEnd,rmin*sinEnd,zmin);
538     pols[ksteps+1][2].set(rmax*cosEnd,rmax*sinEnd,zmin);
539     pols[ksteps+1][3].set(rmax*cosEnd,rmax*sinEnd,zmax);
540 
541     // set envelope and calculate extent
542     std::vector<const G4ThreeVectorList *> polygons;
543     polygons.resize(ksteps+2);
544     for (G4int k=0; k<ksteps+2; ++k) { polygons[k] = &pols[k]; }
545     G4BoundingEnvelope benv(bmin,bmax,polygons);
546     exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
547   }
548   return exist;
549 }
550 
551 //////////////////////////////////////////////////////////////////////////
552 //
553 // Return whether point inside/outside/on surface
554 
555 EInside G4CutTubs::Inside( const G4ThreeVector& p ) const
556 {
557   G4ThreeVector vZ = G4ThreeVector(0,0,fDz);
558   EInside in = kInside;
559 
560   // Check the lower cut plane
561   //
562   G4double zinLow =(p+vZ).dot(fLowNorm);
563   if (zinLow > halfCarTolerance)  { return kOutside; }
564 
565   // Check the higher cut plane
566   //
567   G4double zinHigh = (p-vZ).dot(fHighNorm);
568   if (zinHigh > halfCarTolerance)  { return kOutside; }
569 
570   // Check radius
571   //
572   G4double r2 = p.x()*p.x() + p.y()*p.y() ;
573 
574   G4double tolRMin = fRMin - halfRadTolerance;
575   G4double tolRMax = fRMax + halfRadTolerance;
576   if ( tolRMin < 0 )  { tolRMin = 0; }
577 
578   if (r2 < tolRMin*tolRMin || r2 > tolRMax*tolRMax) { return kOutside; }
579 
580   // Check Phi cut
581   //
582   if(!fPhiFullCutTube)
583   {
584     if ((tolRMin == 0) && (std::fabs(p.x()) <= halfCarTolerance)
585                        && (std::fabs(p.y()) <= halfCarTolerance))
586     {
587       return kSurface;
588     }
589 
590     G4double phi0 = std::atan2(p.y(),p.x());
591     G4double phi1 = phi0 - twopi;
592     G4double phi2 = phi0 + twopi;
593 
594     in = kOutside;
595     G4double sphi = fSPhi - halfAngTolerance;
596     G4double ephi = sphi + fDPhi + kAngTolerance;
597     if ((phi0  >= sphi && phi0  <= ephi) ||
598         (phi1  >= sphi && phi1  <= ephi) ||
599         (phi2  >= sphi && phi2  <= ephi)) in = kSurface;
600     if (in == kOutside)  { return kOutside; }
601 
602     sphi += kAngTolerance;
603     ephi -= kAngTolerance;
604     if ((phi0  >= sphi && phi0  <= ephi) ||
605         (phi1  >= sphi && phi1  <= ephi) ||
606         (phi2  >= sphi && phi2  <= ephi)) in = kInside;
607     if (in == kSurface)  { return kSurface; }
608   }
609 
610   // Check on the Surface for Z
611   //
612   if ((zinLow >= -halfCarTolerance) || (zinHigh >= -halfCarTolerance))
613   {
614     return kSurface;
615   }
616 
617   // Check on the Surface for R
618   //
619   if (fRMin != 0.0) { tolRMin = fRMin + halfRadTolerance; }
620   else       { tolRMin = 0; }
621   tolRMax = fRMax - halfRadTolerance;
622   if (((r2 <= tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax)) &&
623        (r2 >= halfRadTolerance*halfRadTolerance))
624   {
625     return kSurface;
626   }
627 
628   return in;
629 }
630 
631 ///////////////////////////////////////////////////////////////////////////
632 //
633 // Return unit normal of surface closest to p
634 // - note if point on z axis, ignore phi divided sides
635 // - unsafe if point close to z axis a rmin=0 - no explicit checks
636 
637 G4ThreeVector G4CutTubs::SurfaceNormal( const G4ThreeVector& p ) const
638 {
639   G4int noSurfaces = 0;
640   G4double rho, pPhi;
641   G4double distZLow,distZHigh, distRMin, distRMax;
642   G4double distSPhi = kInfinity, distEPhi = kInfinity;
643   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
644 
645   G4ThreeVector norm, sumnorm(0.,0.,0.);
646   G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
647   G4ThreeVector nR, nPs, nPe;
648 
649   rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
650 
651   distRMin = std::fabs(rho - fRMin);
652   distRMax = std::fabs(rho - fRMax);
653 
654   // dist to Low Cut
655   //
656   distZLow =std::fabs((p+vZ).dot(fLowNorm));
657 
658   // dist to High Cut
659   //
660   distZHigh = std::fabs((p-vZ).dot(fHighNorm));
661 
662   if (!fPhiFullCutTube)    // Protected against (0,0,z)
663   {
664     if ( rho > halfCarTolerance )
665     {
666       pPhi = std::atan2(p.y(),p.x());
667 
668       if(pPhi  < fSPhi- halfCarTolerance)           { pPhi += twopi; }
669       else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; }
670 
671       distSPhi = std::fabs(pPhi - fSPhi);
672       distEPhi = std::fabs(pPhi - fSPhi - fDPhi);
673     }
674     else if( fRMin == 0.0 )
675     {
676       distSPhi = 0.;
677       distEPhi = 0.;
678     }
679     nPs = G4ThreeVector( sinSPhi, -cosSPhi, 0 );
680     nPe = G4ThreeVector( -sinEPhi, cosEPhi, 0 );
681   }
682   if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
683 
684   if( distRMax <= halfCarTolerance )
685   {
686     ++noSurfaces;
687     sumnorm += nR;
688   }
689   if( (fRMin != 0.0) && (distRMin <= halfCarTolerance) )
690   {
691     ++noSurfaces;
692     sumnorm -= nR;
693   }
694   if( fDPhi < twopi )
695   {
696     if (distSPhi <= halfAngTolerance)
697     {
698       ++noSurfaces;
699       sumnorm += nPs;
700     }
701     if (distEPhi <= halfAngTolerance)
702     {
703       ++noSurfaces;
704       sumnorm += nPe;
705     }
706   }
707   if (distZLow <= halfCarTolerance)
708   {
709     ++noSurfaces;
710     sumnorm += fLowNorm;
711   }
712   if (distZHigh <= halfCarTolerance)
713   {
714     ++noSurfaces;
715     sumnorm += fHighNorm;
716   }
717   if ( noSurfaces == 0 )
718   {
719 #ifdef G4CSGDEBUG
720     G4Exception("G4CutTubs::SurfaceNormal(p)", "GeomSolids1002",
721                 JustWarning, "Point p is not on surface !?" );
722     G4int oldprc = G4cout.precision(20);
723     G4cout<< "G4CutTubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); "
724           << G4endl << G4endl;
725     G4cout.precision(oldprc) ;
726 #endif
727      norm = ApproxSurfaceNormal(p);
728   }
729   else if ( noSurfaces == 1 )  { norm = sumnorm; }
730   else                         { norm = sumnorm.unit(); }
731 
732   return norm;
733 }
734 
735 /////////////////////////////////////////////////////////////////////////////
736 //
737 // Algorithm for SurfaceNormal() following the original specification
738 // for points not on the surface
739 
740 G4ThreeVector G4CutTubs::ApproxSurfaceNormal( const G4ThreeVector& p ) const
741 {
742   enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};
743 
744   ENorm side ;
745   G4ThreeVector norm ;
746   G4double rho, phi ;
747   G4double distZLow,distZHigh,distZ;
748   G4double distRMin, distRMax, distSPhi, distEPhi, distMin ;
749   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
750 
751   rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
752 
753   distRMin = std::fabs(rho - fRMin) ;
754   distRMax = std::fabs(rho - fRMax) ;
755 
756   //dist to Low Cut
757   //
758   distZLow =std::fabs((p+vZ).dot(fLowNorm));
759 
760   //dist to High Cut
761   //
762   distZHigh = std::fabs((p-vZ).dot(fHighNorm));
763   distZ=std::min(distZLow,distZHigh);
764 
765   if (distRMin < distRMax) // First minimum
766   {
767     if ( distZ < distRMin )
768     {
769        distMin = distZ ;
770        side    = kNZ ;
771     }
772     else
773     {
774       distMin = distRMin ;
775       side    = kNRMin   ;
776     }
777   }
778   else
779   {
780     if ( distZ < distRMax )
781     {
782       distMin = distZ ;
783       side    = kNZ   ;
784     }
785     else
786     {
787       distMin = distRMax ;
788       side    = kNRMax   ;
789     }
790   }
791   if (!fPhiFullCutTube  &&  (rho != 0.0) ) // Protected against (0,0,z)
792   {
793     phi = std::atan2(p.y(),p.x()) ;
794 
795     if ( phi < 0 )  { phi += twopi; }
796 
797     if ( fSPhi < 0 )
798     {
799       distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;
800     }
801     else
802     {
803       distSPhi = std::fabs(phi - fSPhi)*rho ;
804     }
805     distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
806 
807     if (distSPhi < distEPhi) // Find new minimum
808     {
809       if ( distSPhi < distMin )
810       {
811         side = kNSPhi ;
812       }
813     }
814     else
815     {
816       if ( distEPhi < distMin )
817       {
818         side = kNEPhi ;
819       }
820     }
821   }
822   switch ( side )
823   {
824     case kNRMin : // Inner radius
825     {
826       norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ;
827       break ;
828     }
829     case kNRMax : // Outer radius
830     {
831       norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ;
832       break ;
833     }
834     case kNZ :    // + or - dz
835     {
836       if ( distZHigh > distZLow )  { norm = fHighNorm ; }
837       else                         { norm = fLowNorm; }
838       break ;
839     }
840     case kNSPhi:
841     {
842       norm = G4ThreeVector(sinSPhi, -cosSPhi, 0) ;
843       break ;
844     }
845     case kNEPhi:
846     {
847       norm = G4ThreeVector(-sinEPhi, cosEPhi, 0) ;
848       break;
849     }
850     default:      // Should never reach this case ...
851     {
852       DumpInfo();
853       G4Exception("G4CutTubs::ApproxSurfaceNormal()",
854                   "GeomSolids1002", JustWarning,
855                   "Undefined side for valid surface normal to solid.");
856       break ;
857     }
858   }
859   return norm;
860 }
861 
862 ////////////////////////////////////////////////////////////////////
863 //
864 //
865 // Calculate distance to shape from outside, along normalised vector
866 // - return kInfinity if no intersection, or intersection distance <= tolerance
867 //
868 // - Compute the intersection with the z planes
869 //        - if at valid r, phi, return
870 //
871 // -> If point is outer outer radius, compute intersection with rmax
872 //        - if at valid phi,z return
873 //
874 // -> Compute intersection with inner radius, taking largest +ve root
875 //        - if valid (in z,phi), save intersction
876 //
877 //    -> If phi segmented, compute intersections with phi half planes
878 //        - return smallest of valid phi intersections and
879 //          inner radius intersection
880 //
881 // NOTE:
882 // - 'if valid' implies tolerant checking of intersection points
883 
884 G4double G4CutTubs::DistanceToIn( const G4ThreeVector& p,
885                                   const G4ThreeVector& v  ) const
886 {
887   G4double snxt = kInfinity ;      // snxt = default return value
888   G4double tolORMin2, tolIRMax2 ;  // 'generous' radii squared
889   G4double tolORMax2, tolIRMin2;
890   const G4double dRmax = 100.*fRMax;
891   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
892 
893   // Intersection point variables
894   //
895   G4double Dist, sd=0, xi, yi, zi, rho2, inum, iden, cosPsi, Comp,calf ;
896   G4double t1, t2, t3, b, c, d ;     // Quadratic solver variables
897   G4double distZLow,distZHigh;
898   // Calculate tolerant rmin and rmax
899 
900   if (fRMin > kRadTolerance)
901   {
902     tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ;
903     tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ;
904   }
905   else
906   {
907     tolORMin2 = 0.0 ;
908     tolIRMin2 = 0.0 ;
909   }
910   tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ;
911   tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - halfRadTolerance) ;
912 
913   // Intersection with ZCut surfaces
914 
915   // dist to Low Cut
916   //
917   distZLow =(p+vZ).dot(fLowNorm);
918 
919   // dist to High Cut
920   //
921   distZHigh = (p-vZ).dot(fHighNorm);
922 
923   if ( distZLow >= -halfCarTolerance )
924   {
925     calf = v.dot(fLowNorm);
926     if (calf<0)
927     {
928       sd = -distZLow/calf;
929       if(sd < 0.0)  { sd = 0.0; }
930 
931       xi   = p.x() + sd*v.x() ;                // Intersection coords
932       yi   = p.y() + sd*v.y() ;
933       rho2 = xi*xi + yi*yi ;
934 
935       // Check validity of intersection
936 
937       if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
938       {
939         if (!fPhiFullCutTube && (rho2 != 0.0))
940         {
941           // Psi = angle made with central (average) phi of shape
942           //
943           inum   = xi*cosCPhi + yi*sinCPhi ;
944           iden   = std::sqrt(rho2) ;
945           cosPsi = inum/iden ;
946           if (cosPsi >= cosHDPhiIT)  { return sd ; }
947         }
948         else
949         {
950           return sd ;
951         }
952       }
953     }
954     else
955     {
956       if ( sd<halfCarTolerance )
957       {
958         if(calf>=0) { sd=kInfinity; }
959         return sd ;  // On/outside extent, and heading away
960       }              // -> cannot intersect
961     }
962   }
963 
964   if(distZHigh >= -halfCarTolerance )
965   {
966     calf = v.dot(fHighNorm);
967     if (calf<0)
968     {
969       sd = -distZHigh/calf;
970 
971       if(sd < 0.0)  { sd = 0.0; }
972 
973       xi   = p.x() + sd*v.x() ;                // Intersection coords
974       yi   = p.y() + sd*v.y() ;
975       rho2 = xi*xi + yi*yi ;
976 
977       // Check validity of intersection
978 
979       if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
980       {
981         if (!fPhiFullCutTube && (rho2 != 0.0))
982         {
983           // Psi = angle made with central (average) phi of shape
984           //
985           inum   = xi*cosCPhi + yi*sinCPhi ;
986           iden   = std::sqrt(rho2) ;
987           cosPsi = inum/iden ;
988           if (cosPsi >= cosHDPhiIT)  { return sd ; }
989         }
990         else
991         {
992           return sd ;
993         }
994       }
995     }
996     else
997     {
998       if ( sd<halfCarTolerance )
999       {
1000         if(calf>=0) { sd=kInfinity; }
1001         return sd ;  // On/outside extent, and heading away
1002       }              // -> cannot intersect
1003     }
1004   }
1005 
1006   // -> Can not intersect z surfaces
1007   //
1008   // Intersection with rmax (possible return) and rmin (must also check phi)
1009   //
1010   // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1011   //
1012   // Intersects with x^2+y^2=R^2
1013   //
1014   // 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
1015   //            t1                t2                t3
1016 
1017   t1 = 1.0 - v.z()*v.z() ;
1018   t2 = p.x()*v.x() + p.y()*v.y() ;
1019   t3 = p.x()*p.x() + p.y()*p.y() ;
1020   if ( t1 > 0 )        // Check not || to z axis
1021   {
1022     b = t2/t1 ;
1023     c = t3 - fRMax*fRMax ;
1024 
1025     if ((t3 >= tolORMax2) && (t2<0))   // This also handles the tangent case
1026     {
1027       // Try outer cylinder intersection, c=(t3-fRMax*fRMax)/t1;
1028 
1029       c /= t1 ;
1030       d = b*b - c ;
1031 
1032       if (d >= 0)  // If real root
1033       {
1034         sd = c/(-b+std::sqrt(d));
1035         if (sd >= 0)  // If 'forwards'
1036         {
1037           if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
1038           {               // 64 bits systems. Split long distances and recompute
1039             G4double fTerm = sd-std::fmod(sd,dRmax);
1040             sd = fTerm + DistanceToIn(p+fTerm*v,v);
1041           }
1042           // Check z intersection
1043           //
1044           zi = p.z() + sd*v.z() ;
1045           xi = p.x() + sd*v.x() ;
1046           yi = p.y() + sd*v.y() ;
1047           if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1048                -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1049           {
1050             if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1051                  +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1052             {
1053               // Z ok. Check phi intersection if reqd
1054               //
1055               if (fPhiFullCutTube)
1056               {
1057                 return sd ;
1058               }
1059               else
1060               {
1061                 xi     = p.x() + sd*v.x() ;
1062                 yi     = p.y() + sd*v.y() ;
1063                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ;
1064                 if (cosPsi >= cosHDPhiIT)  { return sd ; }
1065               }
1066             }  //  end if std::fabs(zi)
1067           }
1068         }    //  end if (sd>=0)
1069       }      //  end if (d>=0)
1070     }        //  end if (r>=fRMax)
1071     else
1072     {
1073       // Inside outer radius :
1074       // check not inside, and heading through tubs (-> 0 to in)
1075       if ((t3 > tolIRMin2) && (t2 < 0)
1076        && (std::fabs(p.z()) <= std::fabs(GetCutZ(p))-halfCarTolerance ))
1077       {
1078         // Inside both radii, delta r -ve, inside z extent
1079 
1080         if (!fPhiFullCutTube)
1081         {
1082           inum   = p.x()*cosCPhi + p.y()*sinCPhi ;
1083           iden   = std::sqrt(t3) ;
1084           cosPsi = inum/iden ;
1085           if (cosPsi >= cosHDPhiIT)
1086           {
1087             // In the old version, the small negative tangent for the point
1088             // on surface was not taken in account, and returning 0.0 ...
1089             // New version: check the tangent for the point on surface and
1090             // if no intersection, return kInfinity, if intersection instead
1091             // return sd.
1092             //
1093             c = t3-fRMax*fRMax;
1094             if ( c<=0.0 )
1095             {
1096               return 0.0;
1097             }
1098             else
1099             {
1100               c = c/t1 ;
1101               d = b*b-c;
1102               if ( d>=0.0 )
1103               {
1104                 snxt = c/(-b+std::sqrt(d)); // using safe solution
1105                                             // for quadratic equation
1106                 if ( snxt < halfCarTolerance ) { snxt=0; }
1107                 return snxt ;
1108               }
1109               else
1110               {
1111                 return kInfinity;
1112               }
1113             }
1114           }
1115         }
1116         else
1117         {
1118           // In the old version, the small negative tangent for the point
1119           // on surface was not taken in account, and returning 0.0 ...
1120           // New version: check the tangent for the point on surface and
1121           // if no intersection, return kInfinity, if intersection instead
1122           // return sd.
1123           //
1124           c = t3 - fRMax*fRMax;
1125           if ( c<=0.0 )
1126           {
1127             return 0.0;
1128           }
1129           else
1130           {
1131             c = c/t1 ;
1132             d = b*b-c;
1133             if ( d>=0.0 )
1134             {
1135               snxt= c/(-b+std::sqrt(d)); // using safe solution
1136                                          // for quadratic equation
1137               if ( snxt < halfCarTolerance ) { snxt=0; }
1138               return snxt ;
1139             }
1140             else
1141             {
1142               return kInfinity;
1143             }
1144           }
1145         } // end if   (!fPhiFullCutTube)
1146       }   // end if   (t3>tolIRMin2)
1147     }     // end if   (Inside Outer Radius)
1148 
1149     if ( fRMin != 0.0 )    // Try inner cylinder intersection
1150     {
1151       c = (t3 - fRMin*fRMin)/t1 ;
1152       d = b*b - c ;
1153       if ( d >= 0.0 )  // If real root
1154       {
1155         // Always want 2nd root - we are outside and know rmax Hit was bad
1156         // - If on surface of rmin also need farthest root
1157 
1158         sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
1159         if (sd >= -10*halfCarTolerance)  // check forwards
1160         {
1161           // Check z intersection
1162           //
1163           if (sd < 0.0)  { sd = 0.0; }
1164           if (sd>dRmax) // Avoid rounding errors due to precision issues seen
1165           {             // 64 bits systems. Split long distances and recompute
1166             G4double fTerm = sd-std::fmod(sd,dRmax);
1167             sd = fTerm + DistanceToIn(p+fTerm*v,v);
1168           }
1169           zi = p.z() + sd*v.z() ;
1170           xi = p.x() + sd*v.x() ;
1171           yi = p.y() + sd*v.y() ;
1172           if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1173                -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1174           {
1175             if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1176                  +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1177             {
1178               // Z ok. Check phi
1179               //
1180               if ( fPhiFullCutTube )
1181               {
1182                 return sd ;
1183               }
1184               else
1185               {
1186                 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMin ;
1187                 if (cosPsi >= cosHDPhiIT)
1188                 {
1189                   // Good inner radius isect
1190                   // - but earlier phi isect still possible
1191                   //
1192                   snxt = sd ;
1193                 }
1194               }
1195             }      //    end if std::fabs(zi)
1196           }
1197         }          //    end if (sd>=0)
1198       }            //    end if (d>=0)
1199     }              //    end if (fRMin)
1200   }
1201 
1202   // Phi segment intersection
1203   //
1204   // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1205   //
1206   // o NOTE: Large duplication of code between sphi & ephi checks
1207   //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1208   //            intersection check <=0 -> >=0
1209   //         -> use some form of loop Construct ?
1210   //
1211   if ( !fPhiFullCutTube )
1212   {
1213     // First phi surface (Starting phi)
1214     //
1215     Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1216 
1217     if ( Comp < 0 )  // Component in outwards normal dirn
1218     {
1219       Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1220 
1221       if ( Dist < halfCarTolerance )
1222       {
1223         sd = Dist/Comp ;
1224 
1225         if (sd < snxt)
1226         {
1227           if ( sd < 0 )  { sd = 0.0; }
1228           zi = p.z() + sd*v.z() ;
1229           xi = p.x() + sd*v.x() ;
1230           yi = p.y() + sd*v.y() ;
1231           if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1232                -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1233           {
1234             if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1235                  +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1236             {
1237               rho2 = xi*xi + yi*yi ;
1238               if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1239                 || ( (rho2 >  tolORMin2) && (rho2 <  tolIRMin2)
1240                   && ( v.y()*cosSPhi - v.x()*sinSPhi >  0 )
1241                   && ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 )     )
1242                 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1243                   && (v.y()*cosSPhi - v.x()*sinSPhi > 0)
1244                   && (v.x()*cosSPhi + v.y()*sinSPhi < 0) )    )
1245               {
1246                 // z and r intersections good
1247                 // - check intersecting with correct half-plane
1248                 //
1249                 if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = sd; }
1250               }
1251             }   //two Z conditions
1252           }
1253         }
1254       }
1255     }
1256 
1257     // Second phi surface (Ending phi)
1258     //
1259     Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1260 
1261     if (Comp < 0 )  // Component in outwards normal dirn
1262     {
1263       Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1264 
1265       if ( Dist < halfCarTolerance )
1266       {
1267         sd = Dist/Comp ;
1268 
1269         if (sd < snxt)
1270         {
1271           if ( sd < 0 )  { sd = 0; }
1272           zi = p.z() + sd*v.z() ;
1273           xi = p.x() + sd*v.x() ;
1274           yi = p.y() + sd*v.y() ;
1275           if ((-xi*fLowNorm.x()-yi*fLowNorm.y()
1276                -(zi+fDz)*fLowNorm.z())>-halfCarTolerance)
1277           {
1278             if ((-xi*fHighNorm.x()-yi*fHighNorm.y()
1279                  +(fDz-zi)*fHighNorm.z())>-halfCarTolerance)
1280             {
1281               xi   = p.x() + sd*v.x() ;
1282               yi   = p.y() + sd*v.y() ;
1283               rho2 = xi*xi + yi*yi ;
1284               if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1285                   || ( (rho2 > tolORMin2)  && (rho2 < tolIRMin2)
1286                     && (v.x()*sinEPhi - v.y()*cosEPhi >  0)
1287                     && (v.x()*cosEPhi + v.y()*sinEPhi >= 0) )
1288                   || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1289                     && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1290                     && (v.x()*cosEPhi + v.y()*sinEPhi < 0) ) )
1291               {
1292                 // z and r intersections good
1293                 // - check intersecting with correct half-plane
1294                 //
1295                 if ( (yi*cosCPhi-xi*sinCPhi) >= -halfCarTolerance )
1296                 {
1297                   snxt = sd;
1298                 }
1299               }    //?? >=-halfCarTolerance
1300             }
1301           }  // two Z conditions
1302         }
1303       }
1304     }         //  Comp < 0
1305   }           //  !fPhiFullTube
1306   if ( snxt<halfCarTolerance )  { snxt=0; }
1307 
1308   return snxt ;
1309 }
1310 
1311 //////////////////////////////////////////////////////////////////
1312 //
1313 // Calculate distance to shape from outside, along normalised vector
1314 // - return kInfinity if no intersection, or intersection distance <= tolerance
1315 //
1316 // - Compute the intersection with the z planes
1317 //        - if at valid r, phi, return
1318 //
1319 // -> If point is outer outer radius, compute intersection with rmax
1320 //        - if at valid phi,z return
1321 //
1322 // -> Compute intersection with inner radius, taking largest +ve root
1323 //        - if valid (in z,phi), save intersction
1324 //
1325 //    -> If phi segmented, compute intersections with phi half planes
1326 //        - return smallest of valid phi intersections and
1327 //          inner radius intersection
1328 //
1329 // NOTE:
1330 // - Precalculations for phi trigonometry are Done `just in time'
1331 // - `if valid' implies tolerant checking of intersection points
1332 //   Calculate distance (<= actual) to closest surface of shape from outside
1333 // - Calculate distance to z, radial planes
1334 // - Only to phi planes if outside phi extent
1335 // - Return 0 if point inside
1336 
1337 G4double G4CutTubs::DistanceToIn( const G4ThreeVector& p ) const
1338 {
1339   G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho,cosPsi;
1340   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1341 
1342   // Distance to R
1343   //
1344   rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1345 
1346   safRMin = fRMin- rho ;
1347   safRMax = rho - fRMax ;
1348 
1349   // Distances to ZCut(Low/High)
1350 
1351   // Dist to Low Cut
1352   //
1353   safZLow = (p+vZ).dot(fLowNorm);
1354 
1355   // Dist to High Cut
1356   //
1357   safZHigh = (p-vZ).dot(fHighNorm);
1358 
1359   safe = std::max(safZLow,safZHigh);
1360 
1361   if ( safRMin > safe ) { safe = safRMin; }
1362   if ( safRMax> safe )  { safe = safRMax; }
1363 
1364   // Distance to Phi
1365   //
1366   if ( (!fPhiFullCutTube) && ((rho) != 0.0) )
1367    {
1368      // Psi=angle from central phi to point
1369      //
1370      cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
1371 
1372      if ( cosPsi < cosHDPhi )
1373      {
1374        // Point lies outside phi range
1375 
1376        if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
1377        {
1378          safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ;
1379        }
1380        else
1381        {
1382          safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ;
1383        }
1384        if ( safePhi > safe )  { safe = safePhi; }
1385      }
1386    }
1387    if ( safe < 0 )  { safe = 0; }
1388 
1389    return safe ;
1390 }
1391 
1392 //////////////////////////////////////////////////////////////////////////////
1393 //
1394 // Calculate distance to surface of shape from `inside', allowing for tolerance
1395 // - Only Calc rmax intersection if no valid rmin intersection
1396 
1397 G4double G4CutTubs::DistanceToOut( const G4ThreeVector& p,
1398                                    const G4ThreeVector& v,
1399                                    const G4bool calcNorm,
1400                                          G4bool* validNorm,
1401                                          G4ThreeVector* n ) const
1402 {
1403   enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};
1404 
1405   ESide side=kNull , sider=kNull, sidephi=kNull ;
1406   G4double snxt=kInfinity, srd=kInfinity,sz=kInfinity, sphi=kInfinity ;
1407   G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
1408   G4double distZLow,distZHigh,calfH,calfL;
1409   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1410 
1411   // Vars for phi intersection:
1412   //
1413   G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1414 
1415   // Z plane intersection
1416   // Distances to ZCut(Low/High)
1417 
1418   // dist to Low Cut
1419   //
1420   distZLow =(p+vZ).dot(fLowNorm);
1421 
1422   // dist to High Cut
1423   //
1424   distZHigh = (p-vZ).dot(fHighNorm);
1425 
1426   calfH = v.dot(fHighNorm);
1427   calfL = v.dot(fLowNorm);
1428 
1429   if (calfH > 0 )
1430   {
1431     if ( distZHigh < halfCarTolerance )
1432     {
1433       snxt = -distZHigh/calfH ;
1434       side = kPZ ;
1435     }
1436     else
1437     {
1438       if (calcNorm)
1439       {
1440         *n         = G4ThreeVector(0,0,1) ;
1441         *validNorm = true ;
1442       }
1443       return snxt = 0 ;
1444     }
1445  }
1446   if ( calfL>0)
1447   {
1448 
1449     if ( distZLow < halfCarTolerance )
1450     {
1451       sz = -distZLow/calfL ;
1452       if(sz<snxt){
1453       snxt=sz;
1454       side = kMZ ;
1455       }
1456 
1457     }
1458     else
1459     {
1460       if (calcNorm)
1461       {
1462         *n         = G4ThreeVector(0,0,-1) ;
1463         *validNorm = true ;
1464       }
1465       return snxt = 0.0 ;
1466     }
1467   }
1468   if((calfH<=0)&&(calfL<=0))
1469   {
1470     snxt = kInfinity ;    // Travel perpendicular to z axis
1471     side = kNull;
1472   }
1473   // Radial Intersections
1474   //
1475   // Find intersection with cylinders at rmax/rmin
1476   // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1477   //
1478   // Intersects with x^2+y^2=R^2
1479   //
1480   // 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
1481   //
1482   //            t1                t2                    t3
1483 
1484   t1   = 1.0 - v.z()*v.z() ;      // since v normalised
1485   t2   = p.x()*v.x() + p.y()*v.y() ;
1486   t3   = p.x()*p.x() + p.y()*p.y() ;
1487 
1488   if ( snxt > 10*(fDz+fRMax) )  { roi2 = 2*fRMax*fRMax; }
1489   else  { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; }        // radius^2 on +-fDz
1490 
1491   if ( t1 > 0 ) // Check not parallel
1492   {
1493     // Calculate srd, r exit distance
1494 
1495     if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) )
1496     {
1497       // Delta r not negative => leaving via rmax
1498 
1499       deltaR = t3 - fRMax*fRMax ;
1500 
1501       // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
1502       // - avoid sqrt for efficiency
1503 
1504       if ( deltaR < -kRadTolerance*fRMax )
1505       {
1506         b     = t2/t1 ;
1507         c     = deltaR/t1 ;
1508         d2    = b*b-c;
1509         if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1510         else          { srd = 0.; }
1511         sider = kRMax ;
1512       }
1513       else
1514       {
1515         // On tolerant boundary & heading outwards (or perpendicular to)
1516         // outer radial surface -> leaving immediately
1517 
1518         if ( calcNorm )
1519         {
1520           *n         = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1521           *validNorm = true ;
1522         }
1523         return snxt = 0 ; // Leaving by rmax immediately
1524       }
1525     }
1526     else if ( t2 < 0. ) // i.e.  t2 < 0; Possible rmin intersection
1527     {
1528       roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement
1529 
1530       if ( (fRMin != 0.0) && (roMin2 < fRMin*(fRMin - kRadTolerance)) )
1531       {
1532         deltaR = t3 - fRMin*fRMin ;
1533         b      = t2/t1 ;
1534         c      = deltaR/t1 ;
1535         d2     = b*b - c ;
1536 
1537         if ( d2 >= 0 )   // Leaving via rmin
1538         {
1539           // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1540           // - avoid sqrt for efficiency
1541 
1542           if (deltaR > kRadTolerance*fRMin)
1543           {
1544             srd = c/(-b+std::sqrt(d2));
1545             sider = kRMin ;
1546           }
1547           else
1548           {
1549             if ( calcNorm ) { *validNorm = false; }  // Concave side
1550             return snxt = 0.0;
1551           }
1552         }
1553         else    // No rmin intersect -> must be rmax intersect
1554         {
1555           deltaR = t3 - fRMax*fRMax ;
1556           c     = deltaR/t1 ;
1557           d2    = b*b-c;
1558           if( d2 >=0. )
1559           {
1560             srd    = -b + std::sqrt(d2) ;
1561             sider  = kRMax ;
1562           }
1563           else // Case: On the border+t2<kRadTolerance
1564                //       (v is perpendicular to the surface)
1565           {
1566             if (calcNorm)
1567             {
1568               *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1569               *validNorm = true ;
1570             }
1571             return snxt = 0.0;
1572           }
1573         }
1574       }
1575       else if ( roi2 > fRMax*(fRMax + kRadTolerance) )
1576            // No rmin intersect -> must be rmax intersect
1577       {
1578         deltaR = t3 - fRMax*fRMax ;
1579         b      = t2/t1 ;
1580         c      = deltaR/t1;
1581         d2     = b*b-c;
1582         if( d2 >= 0 )
1583         {
1584           srd    = -b + std::sqrt(d2) ;
1585           sider  = kRMax ;
1586         }
1587         else // Case: On the border+t2<kRadTolerance
1588              //       (v is perpendicular to the surface)
1589         {
1590           if (calcNorm)
1591           {
1592             *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1593             *validNorm = true ;
1594           }
1595           return snxt = 0.0;
1596         }
1597       }
1598     }
1599     // Phi Intersection
1600 
1601     if ( !fPhiFullCutTube )
1602     {
1603       // add angle calculation with correction
1604       // of the difference in domain of atan2 and Sphi
1605       //
1606       vphi = std::atan2(v.y(),v.x()) ;
1607 
1608       if ( vphi < fSPhi - halfAngTolerance  )             { vphi += twopi; }
1609       else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1610 
1611 
1612       if ( (p.x() != 0.0) || (p.y() != 0.0) )  // Check if on z axis (rho not needed later)
1613       {
1614         // pDist -ve when inside
1615 
1616         pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1617         pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1618 
1619         // Comp -ve when in direction of outwards normal
1620 
1621         compS   = -sinSPhi*v.x() + cosSPhi*v.y() ;
1622         compE   =  sinEPhi*v.x() - cosEPhi*v.y() ;
1623 
1624         sidephi = kNull;
1625 
1626         if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1627                               && (pDistE <= halfCarTolerance) ) )
1628          || ( (fDPhi >  pi) && ((pDistS <=  halfCarTolerance)
1629                               || (pDistE <=  halfCarTolerance) ) )  )
1630         {
1631           // Inside both phi *full* planes
1632 
1633           if ( compS < 0 )
1634           {
1635             sphi = pDistS/compS ;
1636 
1637             if (sphi >= -halfCarTolerance)
1638             {
1639               xi = p.x() + sphi*v.x() ;
1640               yi = p.y() + sphi*v.y() ;
1641 
1642               // Check intersecting with correct half-plane
1643               // (if not -> no intersect)
1644               //
1645               if( (std::fabs(xi)<=kCarTolerance)
1646                && (std::fabs(yi)<=kCarTolerance) )
1647               {
1648                 sidephi = kSPhi;
1649                 if (((fSPhi-halfAngTolerance)<=vphi)
1650                    &&((fSPhi+fDPhi+halfAngTolerance)>=vphi))
1651                 {
1652                   sphi = kInfinity;
1653                 }
1654               }
1655               else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1656               {
1657                 sphi = kInfinity ;
1658               }
1659               else
1660               {
1661                 sidephi = kSPhi ;
1662                 if ( pDistS > -halfCarTolerance )
1663                 {
1664                   sphi = 0.0 ; // Leave by sphi immediately
1665                 }
1666               }
1667             }
1668             else
1669             {
1670               sphi = kInfinity ;
1671             }
1672           }
1673           else
1674           {
1675             sphi = kInfinity ;
1676           }
1677 
1678           if ( compE < 0 )
1679           {
1680             sphi2 = pDistE/compE ;
1681 
1682             // Only check further if < starting phi intersection
1683             //
1684             if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1685             {
1686               xi = p.x() + sphi2*v.x() ;
1687               yi = p.y() + sphi2*v.y() ;
1688 
1689               if ((std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance))
1690               {
1691                 // Leaving via ending phi
1692                 //
1693                 if( (fSPhi-halfAngTolerance > vphi)
1694                      ||(fSPhi+fDPhi+halfAngTolerance < vphi) )
1695                 {
1696                   sidephi = kEPhi ;
1697                   if ( pDistE <= -halfCarTolerance )  { sphi = sphi2 ; }
1698                   else                                { sphi = 0.0 ;   }
1699                 }
1700               }
1701               else    // Check intersecting with correct half-plane
1702 
1703               if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1704               {
1705                 // Leaving via ending phi
1706                 //
1707                 sidephi = kEPhi ;
1708                 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; }
1709                 else                               { sphi = 0.0 ;   }
1710               }
1711             }
1712           }
1713         }
1714         else
1715         {
1716           sphi = kInfinity ;
1717         }
1718       }
1719       else
1720       {
1721         // On z axis + travel not || to z axis -> if phi of vector direction
1722         // within phi of shape, Step limited by rmax, else Step =0
1723 
1724         if ( (fSPhi - halfAngTolerance <= vphi)
1725            && (vphi <= fSPhi + fDPhi + halfAngTolerance ) )
1726         {
1727           sphi = kInfinity ;
1728         }
1729         else
1730         {
1731           sidephi = kSPhi ; // arbitrary
1732           sphi    = 0.0 ;
1733         }
1734       }
1735       if (sphi < snxt)  // Order intersecttions
1736       {
1737         snxt = sphi ;
1738         side = sidephi ;
1739       }
1740     }
1741     if (srd < snxt)  // Order intersections
1742     {
1743       snxt = srd ;
1744       side = sider ;
1745     }
1746   }
1747   if (calcNorm)
1748   {
1749     switch(side)
1750     {
1751       case kRMax:
1752         // Note: returned vector not normalised
1753         // (divide by fRMax for unit vector)
1754         //
1755         xi = p.x() + snxt*v.x() ;
1756         yi = p.y() + snxt*v.y() ;
1757         *n = G4ThreeVector(xi/fRMax,yi/fRMax,0) ;
1758         *validNorm = true ;
1759         break ;
1760 
1761       case kRMin:
1762         *validNorm = false ;  // Rmin is inconvex
1763         break ;
1764 
1765       case kSPhi:
1766         if ( fDPhi <= pi )
1767         {
1768           *n         = G4ThreeVector(sinSPhi,-cosSPhi,0) ;
1769           *validNorm = true ;
1770         }
1771         else
1772         {
1773           *validNorm = false ;
1774         }
1775         break ;
1776 
1777       case kEPhi:
1778         if (fDPhi <= pi)
1779         {
1780           *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ;
1781           *validNorm = true ;
1782         }
1783         else
1784         {
1785           *validNorm = false ;
1786         }
1787         break ;
1788 
1789       case kPZ:
1790         *n         = fHighNorm ;
1791         *validNorm = true ;
1792         break ;
1793 
1794       case kMZ:
1795         *n         = fLowNorm ;
1796         *validNorm = true ;
1797         break ;
1798 
1799       default:
1800         G4cout << G4endl ;
1801         DumpInfo();
1802         std::ostringstream message;
1803         G4long oldprc = message.precision(16);
1804         message << "Undefined side for valid surface normal to solid."
1805                 << G4endl
1806                 << "Position:"  << G4endl << G4endl
1807                 << "p.x() = "   << p.x()/mm << " mm" << G4endl
1808                 << "p.y() = "   << p.y()/mm << " mm" << G4endl
1809                 << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl
1810                 << "Direction:" << G4endl << G4endl
1811                 << "v.x() = "   << v.x() << G4endl
1812                 << "v.y() = "   << v.y() << G4endl
1813                 << "v.z() = "   << v.z() << G4endl << G4endl
1814                 << "Proposed distance :" << G4endl << G4endl
1815                 << "snxt = "    << snxt/mm << " mm" << G4endl ;
1816         message.precision(oldprc) ;
1817         G4Exception("G4CutTubs::DistanceToOut(p,v,..)", "GeomSolids1002",
1818                     JustWarning, message);
1819         break ;
1820     }
1821   }
1822   if ( snxt<halfCarTolerance )  { snxt=0 ; }
1823   return snxt ;
1824 }
1825 
1826 //////////////////////////////////////////////////////////////////////////
1827 //
1828 // Calculate distance (<=actual) to closest surface of shape from inside
1829 
1830 G4double G4CutTubs::DistanceToOut( const G4ThreeVector& p ) const
1831 {
1832   G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho;
1833   G4ThreeVector vZ=G4ThreeVector(0,0,fDz);
1834 
1835   rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;  // Distance to R
1836 
1837   safRMin =  rho - fRMin ;
1838   safRMax =  fRMax - rho ;
1839 
1840   // Distances to ZCut(Low/High)
1841 
1842   // Dist to Low Cut
1843   //
1844   safZLow = std::fabs((p+vZ).dot(fLowNorm));
1845 
1846   // Dist to High Cut
1847   //
1848   safZHigh = std::fabs((p-vZ).dot(fHighNorm));
1849   safe = std::min(safZLow,safZHigh);
1850 
1851   if ( safRMin < safe ) { safe = safRMin; }
1852   if ( safRMax< safe )  { safe = safRMax; }
1853 
1854   // Check if phi divided, Calc distances closest phi plane
1855   //
1856   if ( !fPhiFullCutTube )
1857   {
1858     if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 )
1859     {
1860       safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
1861     }
1862     else
1863     {
1864       safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
1865     }
1866     if (safePhi < safe)  { safe = safePhi ; }
1867   }
1868   if ( safe < 0 )  { safe = 0; }
1869 
1870   return safe ;
1871 }
1872 
1873 //////////////////////////////////////////////////////////////////////////
1874 //
1875 // Stream object contents to an output stream
1876 
1877 G4GeometryType G4CutTubs::GetEntityType() const
1878 {
1879   return {"G4CutTubs"};
1880 }
1881 
1882 //////////////////////////////////////////////////////////////////////////
1883 //
1884 // Make a clone of the object
1885 //
1886 G4VSolid* G4CutTubs::Clone() const
1887 {
1888   return new G4CutTubs(*this);
1889 }
1890 
1891 //////////////////////////////////////////////////////////////////////////
1892 //
1893 // Stream object contents to an output stream
1894 
1895 std::ostream& G4CutTubs::StreamInfo( std::ostream& os ) const
1896 {
1897   G4long oldprc = os.precision(16);
1898   os << "-----------------------------------------------------------\n"
1899      << "    *** Dump for solid - " << GetName() << " ***\n"
1900      << "    ===================================================\n"
1901      << " Solid type: G4CutTubs\n"
1902      << " Parameters: \n"
1903      << "    inner radius : " << fRMin/mm << " mm \n"
1904      << "    outer radius : " << fRMax/mm << " mm \n"
1905      << "    half length Z: " << fDz/mm << " mm \n"
1906      << "    starting phi : " << fSPhi/degree << " degrees \n"
1907      << "    delta phi    : " << fDPhi/degree << " degrees \n"
1908      << "    low Norm     : " << fLowNorm     << "  \n"
1909      << "    high Norm    : "  <<fHighNorm    << "  \n"
1910      << "-----------------------------------------------------------\n";
1911   os.precision(oldprc);
1912 
1913   return os;
1914 }
1915 
1916 /////////////////////////////////////////////////////////////////////////
1917 //
1918 // GetPointOnSurface
1919 
1920 G4ThreeVector G4CutTubs::GetPointOnSurface() const
1921 {
1922   // Set min and max z
1923   if (fZMin == 0. && fZMax == 0.)
1924   {
1925     G4AutoLock l(&zminmaxMutex);
1926     G4ThreeVector bmin, bmax;
1927     BoundingLimits(bmin,bmax);
1928     fZMin = bmin.z();
1929     fZMax = bmax.z();
1930     l.unlock();
1931   }
1932 
1933   // Set parameters
1934   G4double hmax = fZMax - fZMin;
1935   G4double sphi = fSPhi;
1936   G4double dphi = fDPhi;
1937   G4double rmin = fRMin;
1938   G4double rmax = fRMax;
1939   G4double rrmax = rmax*rmax;
1940   G4double rrmin = rmin*rmin;
1941 
1942   G4ThreeVector nbot = GetLowNorm();
1943   G4ThreeVector ntop = GetHighNorm();
1944 
1945   // Set array of surface areas
1946   G4double sbase = 0.5*dphi*(rrmax - rrmin);
1947   G4double sbot = sbase/std::abs(nbot.z());
1948   G4double stop = sbase/std::abs(ntop.z());
1949   G4double scut = (dphi == twopi) ? 0. : hmax*(rmax - rmin);
1950   G4double ssurf[6] = { scut, scut, sbot, stop, dphi*rmax*hmax, dphi*rmin*hmax };
1951   ssurf[1] += ssurf[0];
1952   ssurf[2] += ssurf[1];
1953   ssurf[3] += ssurf[2];
1954   ssurf[4] += ssurf[3];
1955   ssurf[5] += ssurf[4];
1956 
1957   constexpr G4int ntry = 100000;
1958   for (G4int i=0; i<ntry; ++i)
1959   {
1960     // Select surface
1961     G4double select = ssurf[5]*G4QuickRand();
1962     G4int k = 5;
1963     k -= (G4int)(select <= ssurf[4]);
1964     k -= (G4int)(select <= ssurf[3]);
1965     k -= (G4int)(select <= ssurf[2]);
1966     k -= (G4int)(select <= ssurf[1]);
1967     k -= (G4int)(select <= ssurf[0]);
1968 
1969     // Generate point on selected surface (rejection sampling)
1970     G4ThreeVector p(0,0,0);
1971     switch(k)
1972     {
1973       case 0: // cut at start phi
1974       {
1975         G4double r = rmin + (rmax - rmin)*G4QuickRand();
1976         p.set(r*cosSPhi, r*sinSPhi, fZMin + hmax*G4QuickRand());
1977         break;
1978       }
1979       case 1: // cut at end phi
1980       {
1981         G4double r = rmin + (rmax - rmin)*G4QuickRand();
1982         p.set(r*cosEPhi, r*sinEPhi, fZMin + hmax*G4QuickRand());
1983         break;
1984       }
1985       case 2: // base at low z
1986       {
1987         G4double r = std::sqrt(rrmin + (rrmax - rrmin)*G4QuickRand());
1988         G4double phi = sphi + dphi*G4QuickRand();
1989         G4double x = r*std::cos(phi);
1990         G4double y = r*std::sin(phi);
1991         G4double z = -fDz - (x*nbot.x() + y*nbot.y())/nbot.z();
1992         return {x, y, z};
1993       }
1994       case 3: // base at high z
1995       {
1996         G4double r = std::sqrt(rrmin + (rrmax - rrmin)*G4QuickRand());
1997         G4double phi = sphi + dphi*G4QuickRand();
1998         G4double x = r*std::cos(phi);
1999         G4double y = r*std::sin(phi);
2000         G4double z = fDz - (x*ntop.x() + y*ntop.y())/ntop.z();
2001         return {x, y, z};
2002       }
2003       case 4: // external lateral surface
2004       {
2005         G4double phi = sphi + dphi*G4QuickRand();
2006         G4double z = fZMin + hmax*G4QuickRand();
2007         G4double x = rmax*std::cos(phi);
2008         G4double y = rmax*std::sin(phi);
2009         p.set(x, y, z);
2010         break;
2011       }
2012       case 5: // internal lateral surface
2013       {
2014         G4double phi = sphi + dphi*G4QuickRand();
2015         G4double z = fZMin + hmax*G4QuickRand();
2016         G4double x = rmin*std::cos(phi);
2017         G4double y = rmin*std::sin(phi);
2018         p.set(x, y, z);
2019         break;
2020       }
2021     }
2022     if ((ntop.dot(p) - fDz*ntop.z()) > 0.) continue;
2023     if ((nbot.dot(p) + fDz*nbot.z()) > 0.) continue;
2024     return p;
2025   }
2026   // Just in case, if all attempts to generate a point have failed
2027   // Normally should never happen
2028   G4double x = rmax*std::cos(sphi + 0.5*dphi);
2029   G4double y = rmax*std::sin(sphi + 0.5*dphi);
2030   G4double z = fDz - (x*ntop.x() + y*ntop.y())/ntop.z();
2031   return {x, y, z};
2032 }
2033 
2034 ///////////////////////////////////////////////////////////////////////////
2035 //
2036 // Methods for visualisation
2037 
2038 void G4CutTubs::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
2039 {
2040   scene.AddSolid (*this) ;
2041 }
2042 
2043 G4Polyhedron* G4CutTubs::CreatePolyhedron () const
2044 {
2045   typedef G4double G4double3[3];
2046   typedef G4int G4int4[4];
2047 
2048   auto ph = new G4Polyhedron;
2049   G4Polyhedron *ph1 = new G4PolyhedronTubs (fRMin, fRMax, fDz, fSPhi, fDPhi);
2050   G4int nn=ph1->GetNoVertices();
2051   G4int nf=ph1->GetNoFacets();
2052   auto xyz = new G4double3[nn];  // number of nodes
2053   auto faces = new G4int4[nf] ;  // number of faces
2054 
2055   for(G4int i=0; i<nn; ++i)
2056   {
2057     xyz[i][0]=ph1->GetVertex(i+1).x();
2058     xyz[i][1]=ph1->GetVertex(i+1).y();
2059     G4double tmpZ=ph1->GetVertex(i+1).z();
2060     if(tmpZ>=fDz-kCarTolerance)
2061     {
2062       xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],fDz));
2063     }
2064     else if(tmpZ<=-fDz+kCarTolerance)
2065     {
2066       xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],-fDz));
2067     }
2068     else
2069     {
2070       xyz[i][2]=tmpZ;
2071     }
2072   }
2073   G4int iNodes[4];
2074   G4int* iEdge = nullptr;
2075   G4int n;
2076   for(G4int i=0; i<nf ; ++i)
2077   {
2078     ph1->GetFacet(i+1,n,iNodes,iEdge);
2079     for(G4int k=0; k<n; ++k)
2080     {
2081       faces[i][k]=iNodes[k];
2082     }
2083     for(G4int k=n; k<4; ++k)
2084     {
2085       faces[i][k]=0;
2086     }
2087   }
2088   ph->createPolyhedron(nn,nf,xyz,faces);
2089 
2090   delete [] xyz;
2091   delete [] faces;
2092   delete ph1;
2093 
2094   return ph;
2095 }
2096 
2097 // Auxilary Methods for Solid
2098 
2099 //////////////////////////////////////////////////////////////////////////
2100 //
2101 // Check set of points on the outer lateral surface and return true
2102 // if the cut planes are crossing inside the surface
2103 //
2104 
2105 G4bool G4CutTubs::IsCrossingCutPlanes() const
2106 {
2107   constexpr G4int npoints = 30;
2108 
2109   // set values for calculation of h - distance between
2110   // opposite points on bases
2111   G4ThreeVector nbot = GetLowNorm();
2112   G4ThreeVector ntop = GetHighNorm();
2113   if (std::abs(nbot.z()) < kCarTolerance) return true;
2114   if (std::abs(ntop.z()) < kCarTolerance) return true;
2115   G4double nx = nbot.x()/nbot.z() - ntop.x()/ntop.z();
2116   G4double ny = nbot.y()/nbot.z() - ntop.y()/ntop.z();
2117 
2118   // check points
2119   G4double cosphi = GetCosStartPhi();
2120   G4double sinphi = GetSinStartPhi();
2121   G4double delphi = GetDeltaPhiAngle()/npoints;
2122   G4double cosdel = std::cos(delphi);
2123   G4double sindel = std::sin(delphi);
2124   G4double hzero = 2.*GetZHalfLength()/GetOuterRadius();
2125   for (G4int i=0; i<npoints+1; ++i)
2126   {
2127     G4double h = nx*cosphi + ny*sinphi + hzero;
2128     if (h < 0.) return true;
2129     G4double sintmp = sinphi;
2130     sinphi = sintmp*cosdel + cosphi*sindel;
2131     cosphi = cosphi*cosdel - sintmp*sindel;
2132   }
2133   return false;
2134 }
2135 
2136 ///////////////////////////////////////////////////////////////////////////
2137 //
2138 // Return real Z coordinate of point on Cutted +/- fDZ plane
2139 
2140 G4double G4CutTubs::GetCutZ(const G4ThreeVector& p) const
2141 {
2142   G4double newz = p.z();  // p.z() should be either +fDz or -fDz
2143   if (p.z()<0)
2144   {
2145     if(fLowNorm.z()!=0.)
2146     {
2147        newz = -fDz-(p.x()*fLowNorm.x()+p.y()*fLowNorm.y())/fLowNorm.z();
2148     }
2149   }
2150   else
2151   {
2152     if(fHighNorm.z()!=0.)
2153     {
2154        newz = fDz-(p.x()*fHighNorm.x()+p.y()*fHighNorm.y())/fHighNorm.z();
2155     }
2156   }
2157   return newz;
2158 }
2159 #endif
2160