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