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