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