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