Loading...
Searching...
No Matches
G4Cons Class Reference
#include <G4Cons.hh>
Inheritance diagram for G4Cons:
Detailed Description
Constructor & Destructor Documentation
◆ G4Cons() [1/3]
| G4Cons::G4Cons | ( | const G4String & | pName, |
| G4double | pRmin1, | ||
| G4double | pRmax1, | ||
| G4double | pRmin2, | ||
| G4double | pRmax2, | ||
| G4double | pDz, | ||
| G4double | pSPhi, | ||
| G4double | pDPhi | ||
| ) |
Definition at line 79 of file G4Cons.cc.
84 : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2),
85 fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
86{
89
90 // Check z-len
91 //
92 if ( pDz < 0 )
93 {
94 std::ostringstream message;
96 << " hZ = " << pDz;
98 FatalException, message);
99 }
100
101 // Check radii
102 //
103 if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
104 {
105 std::ostringstream message;
107 << " pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2
108 << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2;
110 FatalException, message) ;
111 }
112 if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
113 if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }
114
115 // Check angles
116 //
117 CheckPhiAngles(pSPhi, pDPhi);
118}
Definition: G4CSGSolid.hh:50
G4double GetRadialTolerance() const
Definition: G4GeometryTolerance.cc:97
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:73
G4double GetAngularTolerance() const
Definition: G4GeometryTolerance.cc:92
G4String GetName() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
◆ ~G4Cons()
◆ G4Cons() [2/3]
| G4Cons::G4Cons | ( | __void__ & | a | ) |
Definition at line 125 of file G4Cons.cc.
126 : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
127 fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.),
128 fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
129 cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
130 fPhiFullCone(false)
131{
132}
◆ G4Cons() [3/3]
Definition at line 146 of file G4Cons.cc.
147 : G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance),
148 kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
149 fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi),
150 fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
151 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
152 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
153 cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone)
154{
155}
Member Function Documentation
◆ CalculateExtent()
|
virtual |
Implements G4VSolid.
Definition at line 261 of file G4Cons.cc.
266{
268 && (fRmin1 == 0) && (fRmin2 == 0) )
269 {
270 // Special case handling for unrotated solid cones
271 // Compute z/x/y mins and maxs for bounding box respecting limits,
272 // with early returns if outside limits. Then switch() on pAxis,
273 // and compute exact x and y limit for x/y case
274
275 G4double xoffset, xMin, xMax ;
276 G4double yoffset, yMin, yMax ;
277 G4double zoffset, zMin, zMax ;
278
279 G4double diff1, diff2, delta, maxDiff, newMin, newMax, RMax ;
280 G4double xoff1, xoff2, yoff1, yoff2 ;
281
283 zMin = zoffset - fDz ;
284 zMax = zoffset + fDz ;
285
287 {
290 {
291 return false ;
292 }
293 else
294 {
296 {
297 zMin = pVoxelLimit.GetMinZExtent() ;
298 }
300 {
301 zMax = pVoxelLimit.GetMaxZExtent() ;
302 }
303 }
304 }
306 RMax = (fRmax2 >= fRmax1) ? zMax : zMin ;
307 xMax = xoffset + (fRmax1 + fRmax2)*0.5 +
308 (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
309 xMin = 2*xoffset-xMax ;
310
312 {
315 {
316 return false ;
317 }
318 else
319 {
321 {
322 xMin = pVoxelLimit.GetMinXExtent() ;
323 }
325 {
326 xMax=pVoxelLimit.GetMaxXExtent() ;
327 }
328 }
329 }
331 yMax = yoffset + (fRmax1 + fRmax2)*0.5 +
332 (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
333 yMin = 2*yoffset-yMax ;
334 RMax = yMax - yoffset ; // = max radius due to Zmax/Zmin cuttings
335
337 {
340 {
341 return false ;
342 }
343 else
344 {
346 {
347 yMin = pVoxelLimit.GetMinYExtent() ;
348 }
350 {
351 yMax = pVoxelLimit.GetMaxYExtent() ;
352 }
353 }
354 }
355 switch (pAxis) // Known to cut cones
356 {
358 yoff1 = yoffset - yMin ;
359 yoff2 = yMax - yoffset ;
360
361 if ((yoff1 >= 0) && (yoff2 >= 0)) // Y limits cross max/min x
362 { // => no change
363 pMin = xMin ;
364 pMax = xMax ;
365 }
366 else
367 {
368 // Y limits don't cross max/min x => compute max delta x,
369 // hence new mins/maxs
370 delta=RMax*RMax-yoff1*yoff1;
371 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
372 delta=RMax*RMax-yoff2*yoff2;
373 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
374 maxDiff = (diff1>diff2) ? diff1:diff2 ;
375 newMin = xoffset - maxDiff ;
376 newMax = xoffset + maxDiff ;
377 pMin = ( newMin < xMin ) ? xMin : newMin ;
378 pMax = ( newMax > xMax) ? xMax : newMax ;
379 }
380 break ;
381
383 xoff1 = xoffset - xMin ;
384 xoff2 = xMax - xoffset ;
385
386 if ((xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
387 { // => no change
388 pMin = yMin ;
389 pMax = yMax ;
390 }
391 else
392 {
393 // X limits don't cross max/min y => compute max delta y,
394 // hence new mins/maxs
395 delta=RMax*RMax-xoff1*xoff1;
396 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
397 delta=RMax*RMax-xoff2*xoff2;
398 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
399 maxDiff = (diff1 > diff2) ? diff1:diff2 ;
400 newMin = yoffset - maxDiff ;
401 newMax = yoffset + maxDiff ;
402 pMin = (newMin < yMin) ? yMin : newMin ;
403 pMax = (newMax > yMax) ? yMax : newMax ;
404 }
405 break ;
406
408 pMin = zMin ;
409 pMax = zMax ;
410 break ;
411
412 default:
413 break ;
414 }
415 pMin -= kCarTolerance ;
416 pMax += kCarTolerance ;
417
418 return true ;
419 }
420 else // Calculate rotated vertex coordinates
421 {
422 G4int i, noEntries, noBetweenSections4 ;
423 G4bool existsAfterClip = false ;
424 G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform) ;
425
426 pMin = +kInfinity ;
427 pMax = -kInfinity ;
428
429 noEntries = vertices->size() ;
430 noBetweenSections4 = noEntries-4 ;
431
432 for ( i = 0 ; i < noEntries ; i += 4 )
433 {
434 ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
435 }
436 for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
437 {
438 ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
439 }
440 if ( (pMin != kInfinity) || (pMax != -kInfinity) )
441 {
442 existsAfterClip = true ;
443
444 // Add 2*tolerance to avoid precision troubles
445
446 pMin -= kCarTolerance ;
447 pMax += kCarTolerance ;
448 }
449 else
450 {
451 // Check for case where completely enveloping clipping volume
452 // If point inside then we are confident that the solid completely
453 // envelopes the clipping volume. Hence set min/max extents according
454 // to clipping volume extents along the specified axis.
455
456 G4ThreeVector clipCentre(
460
462 {
463 existsAfterClip = true ;
464 pMin = pVoxelLimit.GetMinExtent(pAxis) ;
465 pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
466 }
467 }
468 delete vertices ;
469 return existsAfterClip ;
470 }
471}
Definition: ThreeVector.h:41
double z() const
double x() const
double y() const
G4bool IsRotated() const
G4AffineTransform Inverse() const
G4ThreeVector NetTranslation() const
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:376
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:345
G4double GetMinExtent(const EAxis pAxis) const
G4bool IsYLimited() const
G4double GetMinZExtent() const
G4bool IsXLimited() const
G4double GetMaxExtent(const EAxis pAxis) const
G4double GetMaxYExtent() const
G4double GetMaxZExtent() const
G4double GetMinYExtent() const
G4double GetMinXExtent() const
G4bool IsZLimited() const
G4double GetMaxXExtent() const
◆ Clone()
|
virtual |
◆ ComputeDimensions()
|
virtual |
Reimplemented from G4VSolid.
Definition at line 249 of file G4Cons.cc.
252{
254}
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
Definition: G4VPVParameterisation.hh:94
◆ CreateNURBS()
|
virtual |
Reimplemented from G4VSolid.
Definition at line 2386 of file G4Cons.cc.
2387{
2388 G4double RMax = (fRmax2 >= fRmax1) ? fRmax2 : fRmax1 ;
2390}
Definition: G4NURBSbox.hh:45
◆ CreatePolyhedron()
|
virtual |
Reimplemented from G4VSolid.
◆ DescribeYourselfTo()
|
virtual |
◆ DistanceToIn() [1/2]
|
virtual |
Implements G4VSolid.
Definition at line 1383 of file G4Cons.cc.
1384{
1385 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
1386 G4double tanRMin, secRMin, pRMin ;
1387 G4double tanRMax, secRMax, pRMax ;
1388
1390 safeZ = std::fabs(p.z()) - fDz ;
1391
1392 if ( fRmin1 || fRmin2 )
1393 {
1394 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1395 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1396 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
1397 safeR1 = (pRMin - rho)/secRMin ;
1398
1399 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1400 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1401 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1402 safeR2 = (rho - pRMax)/secRMax ;
1403
1404 if ( safeR1 > safeR2) { safe = safeR1; }
1405 else { safe = safeR2; }
1406 }
1407 else
1408 {
1409 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1410 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1411 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
1412 safe = (rho - pRMax)/secRMax ;
1413 }
1414 if ( safeZ > safe ) { safe = safeZ; }
1415
1416 if ( !fPhiFullCone && rho )
1417 {
1418 // Psi=angle from central phi to point
1419
1421
1422 if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range
1423 {
1425 {
1427 }
1428 else
1429 {
1431 }
1432 if ( safePhi > safe ) { safe = safePhi; }
1433 }
1434 }
1435 if ( safe < 0.0 ) { safe = 0.0; }
1436
1437 return safe ;
1438}
◆ DistanceToIn() [2/2]
|
virtual |
Implements G4VSolid.
Definition at line 712 of file G4Cons.cc.
714{
719
721 G4double tanRMin,secRMin,rMinAv,rMinOAv ;
722 G4double rout,rin ;
723
725 G4double tolORMax2,tolIRMax,tolIRMax2 ;
726 G4double tolODz,tolIDz ;
727
729
731 G4double nt1,nt2,nt3 ;
732 G4double Comp ;
733
734 G4ThreeVector Normal;
735
736 // Cone Precalcs
737
738 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
739 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
740 rMinAv = (fRmin1 + fRmin2)*0.5 ;
741
742 if (rMinAv > halfRadTolerance)
743 {
744 rMinOAv = rMinAv - halfRadTolerance ;
745 }
746 else
747 {
748 rMinOAv = 0.0 ;
749 }
750 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
751 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
752 rMaxAv = (fRmax1 + fRmax2)*0.5 ;
753 rMaxOAv = rMaxAv + halfRadTolerance ;
754
755 // Intersection with z-surfaces
756
757 tolIDz = fDz - halfCarTolerance ;
758 tolODz = fDz + halfCarTolerance ;
759
761 {
763 {
765
766 if( sd < 0.0 ) { sd = 0.0; } // negative dist -> zero
767
770 rhoi2 = xi*xi + yi*yi ;
771
772 // Check validity of intersection
773 // Calculate (outer) tolerant radi^2 at intersecion
774
776 {
777 tolORMin = fRmin1 - halfRadTolerance*secRMin ;
778 tolIRMin = fRmin1 + halfRadTolerance*secRMin ;
779 tolIRMax = fRmax1 - halfRadTolerance*secRMin ;
780 tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)*
781 (fRmax1 + halfRadTolerance*secRMax) ;
782 }
783 else
784 {
785 tolORMin = fRmin2 - halfRadTolerance*secRMin ;
786 tolIRMin = fRmin2 + halfRadTolerance*secRMin ;
787 tolIRMax = fRmax2 - halfRadTolerance*secRMin ;
788 tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)*
789 (fRmax2 + halfRadTolerance*secRMax) ;
790 }
791 if ( tolORMin > 0 )
792 {
793 tolORMin2 = tolORMin*tolORMin ;
794 tolIRMin2 = tolIRMin*tolIRMin ;
795 }
796 else
797 {
798 tolORMin2 = 0.0 ;
799 tolIRMin2 = 0.0 ;
800 }
801 if ( tolIRMax > 0 ) { tolIRMax2 = tolIRMax*tolIRMax; }
802 else { tolIRMax2 = 0.0; }
803
804 if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) )
805 {
806 if ( !fPhiFullCone && rhoi2 )
807 {
808 // Psi = angle made with central (average) phi of shape
809
810 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
811
812 if (cosPsi >= cosHDPhiIT) { return sd; }
813 }
814 else
815 {
816 return sd;
817 }
818 }
819 }
820 else // On/outside extent, and heading away -> cannot intersect
821 {
822 return snxt ;
823 }
824 }
825
826// ----> Can not intersect z surfaces
827
828
829// Intersection with outer cone (possible return) and
830// inner cone (must also check phi)
831//
832// Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
833//
834// Intersects with x^2+y^2=(a*z+b)^2
835//
836// where a=tanRMax or tanRMin
837// b=rMaxAv or rMinAv
838//
839// (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
840// t1 t2 t3
841//
842// \--------u-------/ \-----------v----------/ \---------w--------/
843//
844
848 rin = tanRMin*p.z() + rMinAv ;
849 rout = tanRMax*p.z() + rMaxAv ;
850
851 // Outer Cone Intersection
852 // Must be outside/on outer cone for valid intersection
853
855 nt2 = t2 - tanRMax*v.z()*rout ;
856 nt3 = t3 - rout*rout ;
857
858 if (std::fabs(nt1) > kRadTolerance) // Equation quadratic => 2 roots
859 {
860 b = nt2/nt1;
861 c = nt3/nt1;
862 d = b*b-c ;
863 if ( (nt3 > rout*rout*kRadTolerance*kRadTolerance*secRMax*secRMax)
864 || (rout < 0) )
865 {
866 // If outside real cone (should be rho-rout>kRadTolerance*0.5
867 // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy
868
869 if (d >= 0)
870 {
871
872 if ((rout < 0) && (nt3 <= 0))
873 {
874 // Inside `shadow cone' with -ve radius
875 // -> 2nd root could be on real cone
876
877 if (b>0) { sd = c/(-b-std::sqrt(d)); }
878 else { sd = -b + std::sqrt(d); }
879 }
880 else
881 {
882 if ((b <= 0) && (c >= 0)) // both >=0, try smaller root
883 {
884 sd=c/(-b+std::sqrt(d));
885 }
886 else
887 {
888 if ( c <= 0 ) // second >=0
889 {
890 sd = -b + std::sqrt(d) ;
891 }
892 else // both negative, travel away
893 {
894 return kInfinity ;
895 }
896 }
897 }
898 if ( sd > 0 ) // If 'forwards'. Check z intersection
899 {
900 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
901 { // 64 bits systems. Split long distances and recompute
902 G4double fTerm = sd-std::fmod(sd,dRmax);
903 sd = fTerm + DistanceToIn(p+fTerm*v,v);
904 }
906
907 if (std::fabs(zi) <= tolODz)
908 {
909 // Z ok. Check phi intersection if reqd
910
911 if ( fPhiFullCone ) { return sd; }
912 else
913 {
916 ri = rMaxAv + zi*tanRMax ;
917 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
918
919 if ( cosPsi >= cosHDPhiIT ) { return sd; }
920 }
921 }
922 } // end if (sd>0)
923 }
924 }
925 else
926 {
927 // Inside outer cone
928 // check not inside, and heading through G4Cons (-> 0 to in)
929
930 if ( ( t3 > (rin + halfRadTolerance*secRMin)*
931 (rin + halfRadTolerance*secRMin) )
932 && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) )
933 {
934 // Inside cones, delta r -ve, inside z extent
935 // Point is on the Surface => check Direction using Normal.dot(v)
936
937 xi = p.x() ;
938 yi = p.y() ;
939 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
940 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
941 if ( !fPhiFullCone )
942 {
944 if ( cosPsi >= cosHDPhiIT )
945 {
947 }
948 }
949 else
950 {
952 }
953 }
954 }
955 }
956 else // Single root case
957 {
958 if ( std::fabs(nt2) > kRadTolerance )
959 {
960 sd = -0.5*nt3/nt2 ;
961
962 if ( sd < 0 ) { return kInfinity; } // travel away
963 else // sd >= 0, If 'forwards'. Check z intersection
964 {
966
967 if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
968 {
969 // Z ok. Check phi intersection if reqd
970
971 if ( fPhiFullCone ) { return sd; }
972 else
973 {
976 ri = rMaxAv + zi*tanRMax ;
977 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
978
979 if (cosPsi >= cosHDPhiIT) { return sd; }
980 }
981 }
982 }
983 }
984 else // travel || cone surface from its origin
985 {
986 sd = kInfinity ;
987 }
988 }
989
990 // Inner Cone Intersection
991 // o Space is divided into 3 areas:
992 // 1) Radius greater than real inner cone & imaginary cone & outside
993 // tolerance
994 // 2) Radius less than inner or imaginary cone & outside tolarance
995 // 3) Within tolerance of real or imaginary cones
996 // - Extra checks needed for 3's intersections
997 // => lots of duplicated code
998
999 if (rMinAv)
1000 {
1002 nt2 = t2 - tanRMin*v.z()*rin ;
1003 nt3 = t3 - rin*rin ;
1004
1005 if ( nt1 )
1006 {
1007 if ( nt3 > rin*kRadTolerance*secRMin )
1008 {
1009 // At radius greater than real & imaginary cones
1010 // -> 2nd root, with zi check
1011
1012 b = nt2/nt1 ;
1013 c = nt3/nt1 ;
1014 d = b*b-c ;
1015 if (d >= 0) // > 0
1016 {
1017 if(b>0){sd = c/( -b-std::sqrt(d));}
1018 else {sd = -b + std::sqrt(d) ;}
1019
1020 if ( sd >= 0 ) // > 0
1021 {
1022 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
1023 { // 64 bits systems. Split long distance and recompute
1024 G4double fTerm = sd-std::fmod(sd,dRmax);
1025 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1026 }
1028
1029 if ( std::fabs(zi) <= tolODz )
1030 {
1031 if ( !fPhiFullCone )
1032 {
1035 ri = rMinAv + zi*tanRMin ;
1036 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1037
1038 if (cosPsi >= cosHDPhiIT)
1039 {
1040 if ( sd > halfRadTolerance ) { snxt=sd; }
1041 else
1042 {
1043 // Calculate a normal vector in order to check Direction
1044
1045 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1046 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1048 }
1049 }
1050 }
1051 else
1052 {
1053 if ( sd > halfRadTolerance ) { return sd; }
1054 else
1055 {
1056 // Calculate a normal vector in order to check Direction
1057
1060 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1061 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1063 }
1064 }
1065 }
1066 }
1067 }
1068 }
1069 else if ( nt3 < -rin*kRadTolerance*secRMin )
1070 {
1071 // Within radius of inner cone (real or imaginary)
1072 // -> Try 2nd root, with checking intersection is with real cone
1073 // -> If check fails, try 1st root, also checking intersection is
1074 // on real cone
1075
1076 b = nt2/nt1 ;
1077 c = nt3/nt1 ;
1078 d = b*b - c ;
1079
1080 if ( d >= 0 ) // > 0
1081 {
1082 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1083 else { sd = -b + std::sqrt(d); }
1085 ri = rMinAv + zi*tanRMin ;
1086
1087 if ( ri > 0 )
1088 {
1089 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd > 0
1090 {
1091 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1092 { // seen on 64 bits systems. Split and recompute
1093 G4double fTerm = sd-std::fmod(sd,dRmax);
1094 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1095 }
1096 if ( !fPhiFullCone )
1097 {
1100 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1101
1102 if (cosPsi >= cosHDPhiOT)
1103 {
1104 if ( sd > halfRadTolerance ) { snxt=sd; }
1105 else
1106 {
1107 // Calculate a normal vector in order to check Direction
1108
1109 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1110 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1112 }
1113 }
1114 }
1115 else
1116 {
1117 if( sd > halfRadTolerance ) { return sd; }
1118 else
1119 {
1120 // Calculate a normal vector in order to check Direction
1121
1124 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1125 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1127 }
1128 }
1129 }
1130 }
1131 else
1132 {
1133 if (b>0) { sd = -b - std::sqrt(d); }
1134 else { sd = c/(-b+std::sqrt(d)); }
1136 ri = rMinAv + zi*tanRMin ;
1137
1138 if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1139 {
1140 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1141 { // seen on 64 bits systems. Split and recompute
1142 G4double fTerm = sd-std::fmod(sd,dRmax);
1143 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1144 }
1145 if ( !fPhiFullCone )
1146 {
1149 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1150
1151 if (cosPsi >= cosHDPhiIT)
1152 {
1153 if ( sd > halfRadTolerance ) { snxt=sd; }
1154 else
1155 {
1156 // Calculate a normal vector in order to check Direction
1157
1158 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1159 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
1161 }
1162 }
1163 }
1164 else
1165 {
1166 if ( sd > halfRadTolerance ) { return sd; }
1167 else
1168 {
1169 // Calculate a normal vector in order to check Direction
1170
1173 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1174 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1176 }
1177 }
1178 }
1179 }
1180 }
1181 }
1182 else
1183 {
1184 // Within kRadTol*0.5 of inner cone (real OR imaginary)
1185 // ----> Check not travelling through (=>0 to in)
1186 // ----> if not:
1187 // -2nd root with validity check
1188
1190 {
1191 if ( nt2 > 0 )
1192 {
1193 // Inside inner real cone, heading outwards, inside z range
1194
1195 if ( !fPhiFullCone )
1196 {
1198
1199 if (cosPsi >= cosHDPhiIT) { return 0.0; }
1200 }
1201 else { return 0.0; }
1202 }
1203 else
1204 {
1205 // Within z extent, but not travelling through
1206 // -> 2nd root or kInfinity if 1st root on imaginary cone
1207
1208 b = nt2/nt1 ;
1209 c = nt3/nt1 ;
1210 d = b*b - c ;
1211
1212 if ( d >= 0 ) // > 0
1213 {
1214 if (b>0) { sd = -b - std::sqrt(d); }
1215 else { sd = c/(-b+std::sqrt(d)); }
1217 ri = rMinAv + zi*tanRMin ;
1218
1219 if ( ri > 0 ) // 2nd root
1220 {
1221 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1222 else { sd = -b + std::sqrt(d); }
1223
1225
1226 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1227 {
1228 if ( sd>dRmax ) // Avoid rounding errors due to precision issue
1229 { // seen on 64 bits systems. Split and recompute
1230 G4double fTerm = sd-std::fmod(sd,dRmax);
1231 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1232 }
1233 if ( !fPhiFullCone )
1234 {
1237 ri = rMinAv + zi*tanRMin ;
1238 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
1239
1240 if ( cosPsi >= cosHDPhiIT ) { snxt = sd; }
1241 }
1242 else { return sd; }
1243 }
1244 }
1245 else { return kInfinity; }
1246 }
1247 }
1248 }
1249 else // 2nd root
1250 {
1251 b = nt2/nt1 ;
1252 c = nt3/nt1 ;
1253 d = b*b - c ;
1254
1255 if ( d > 0 )
1256 {
1257 if (b>0) { sd = c/(-b-std::sqrt(d)); }
1258 else { sd = -b + std::sqrt(d) ; }
1260
1261 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0
1262 {
1263 if ( sd>dRmax ) // Avoid rounding errors due to precision issues
1264 { // seen on 64 bits systems. Split and recompute
1265 G4double fTerm = sd-std::fmod(sd,dRmax);
1266 sd = fTerm + DistanceToIn(p+fTerm*v,v);
1267 }
1268 if ( !fPhiFullCone )
1269 {
1272 ri = rMinAv + zi*tanRMin ;
1273 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;
1274
1275 if (cosPsi >= cosHDPhiIT) { snxt = sd; }
1276 }
1277 else { return sd; }
1278 }
1279 }
1280 }
1281 }
1282 }
1283 }
1284
1285 // Phi segment intersection
1286 //
1287 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1288 //
1289 // o NOTE: Large duplication of code between sphi & ephi checks
1290 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1291 // intersection check <=0 -> >=0
1292 // -> Should use some form of loop Construct
1293
1294 if ( !fPhiFullCone )
1295 {
1296 // First phi surface (starting phi)
1297
1299
1300 if ( Comp < 0 ) // Component in outwards normal dirn
1301 {
1303
1304 if (Dist < halfCarTolerance)
1305 {
1306 sd = Dist/Comp ;
1307
1308 if ( sd < snxt )
1309 {
1310 if ( sd < 0 ) { sd = 0.0; }
1311
1313
1314 if ( std::fabs(zi) <= tolODz )
1315 {
1318 rhoi2 = xi*xi + yi*yi ;
1319 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1320 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1321
1322 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1323 {
1324 // z and r intersections good - check intersecting with
1325 // correct half-plane
1326
1327 if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) { snxt = sd; }
1328 }
1329 }
1330 }
1331 }
1332 }
1333
1334 // Second phi surface (Ending phi)
1335
1337
1338 if ( Comp < 0 ) // Component in outwards normal dirn
1339 {
1341 if (Dist < halfCarTolerance)
1342 {
1343 sd = Dist/Comp ;
1344
1345 if ( sd < snxt )
1346 {
1347 if ( sd < 0 ) { sd = 0.0; }
1348
1350
1351 if (std::fabs(zi) <= tolODz)
1352 {
1355 rhoi2 = xi*xi + yi*yi ;
1356 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
1357 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
1358
1359 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
1360 {
1361 // z and r intersections good - check intersecting with
1362 // correct half-plane
1363
1364 if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) { snxt = sd; }
1365 }
1366 }
1367 }
1368 }
1369 }
1370 }
1371 if (snxt < halfCarTolerance) { snxt = 0.; }
1372
1373 return snxt ;
1374}
double dot(const Hep3Vector &) const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Cons.cc:712
Referenced by DistanceToIn().
◆ DistanceToOut() [1/2]
|
virtual |
Implements G4VSolid.
Definition at line 2072 of file G4Cons.cc.
2073{
2074 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
2075 G4double tanRMin, secRMin, pRMin;
2076 G4double tanRMax, secRMax, pRMax;
2077
2078#ifdef G4CSGDEBUG
2080 {
2083 DumpInfo();
2091 {
2094 }
2095 G4cout.precision(oldprc) ;
2098 }
2099#endif
2100
2102 safeZ = fDz - std::fabs(p.z()) ;
2103
2104 if (fRmin1 || fRmin2)
2105 {
2106 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
2107 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
2108 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
2109 safeR1 = (rho - pRMin)/secRMin ;
2110 }
2111 else
2112 {
2113 safeR1 = kInfinity ;
2114 }
2115
2116 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
2117 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
2118 pRMax = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ;
2119 safeR2 = (pRMax - rho)/secRMax ;
2120
2121 if (safeR1 < safeR2) { safe = safeR1; }
2122 else { safe = safeR2; }
2123 if (safeZ < safe) { safe = safeZ ; }
2124
2125 // Check if phi divided, Calc distances closest phi plane
2126
2127 if (!fPhiFullCone)
2128 {
2129 // Above/below central phi of G4Cons?
2130
2132 {
2134 }
2135 else
2136 {
2138 }
2139 if (safePhi < safe) { safe = safePhi; }
2140 }
2141 if ( safe < 0 ) { safe = 0; }
2142
2143 return safe ;
2144}
G4DLLIMPORT std::ostream G4cout
void DumpInfo() const
◆ DistanceToOut() [2/2]
|
virtual |
Implements G4VSolid.
Definition at line 1445 of file G4Cons.cc.
1450{
1451 ESide side = kNull, sider = kNull, sidephi = kNull;
1452
1456
1457 G4double snxt,srd,sphi,pdist ;
1458
1461
1462 G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ;
1463 G4double b, c, d, sr2, sr3 ;
1464
1465 // Vars for intersection within tolerance
1466
1467 ESide sidetol = kNull ;
1468 G4double slentol = kInfinity ;
1469
1470 // Vars for phi intersection:
1471
1472 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ;
1473 G4double zi, ri, deltaRoi2 ;
1474
1475 // Z plane intersection
1476
1478 {
1479 pdist = fDz - p.z() ;
1480
1481 if (pdist > halfCarTolerance)
1482 {
1483 snxt = pdist/v.z() ;
1484 side = kPZ ;
1485 }
1486 else
1487 {
1488 if (calcNorm)
1489 {
1491 *validNorm = true ;
1492 }
1493 return snxt = 0.0;
1494 }
1495 }
1497 {
1498 pdist = fDz + p.z() ;
1499
1500 if ( pdist > halfCarTolerance)
1501 {
1502 snxt = -pdist/v.z() ;
1503 side = kMZ ;
1504 }
1505 else
1506 {
1507 if ( calcNorm )
1508 {
1510 *validNorm = true ;
1511 }
1512 return snxt = 0.0 ;
1513 }
1514 }
1515 else // Travel perpendicular to z axis
1516 {
1517 snxt = kInfinity ;
1518 side = kNull ;
1519 }
1520
1521 // Radial Intersections
1522 //
1523 // Intersection with outer cone (possible return) and
1524 // inner cone (must also check phi)
1525 //
1526 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1527 //
1528 // Intersects with x^2+y^2=(a*z+b)^2
1529 //
1530 // where a=tanRMax or tanRMin
1531 // b=rMaxAv or rMinAv
1532 //
1533 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
1534 // t1 t2 t3
1535 //
1536 // \--------u-------/ \-----------v----------/ \---------w--------/
1537
1538 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
1539 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
1540 rMaxAv = (fRmax1 + fRmax2)*0.5 ;
1541
1542
1546 rout = tanRMax*p.z() + rMaxAv ;
1547
1549 nt2 = t2 - tanRMax*v.z()*rout ;
1550 nt3 = t3 - rout*rout ;
1551
1553 {
1554 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1555 - fRmax2*(fRmax2 + kRadTolerance*secRMax);
1556 }
1558 {
1559 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
1560 - fRmax1*(fRmax1 + kRadTolerance*secRMax);
1561 }
1562 else
1563 {
1564 deltaRoi2 = 1.0;
1565 }
1566
1567 if ( nt1 && (deltaRoi2 > 0.0) )
1568 {
1569 // Equation quadratic => 2 roots : second root must be leaving
1570
1571 b = nt2/nt1 ;
1572 c = nt3/nt1 ;
1573 d = b*b - c ;
1574
1575 if ( d >= 0 )
1576 {
1577 // Check if on outer cone & heading outwards
1578 // NOTE: Should use rho-rout>-kRadTolerance*0.5
1579
1580 if (nt3 > -halfRadTolerance && nt2 >= 0 )
1581 {
1582 if (calcNorm)
1583 {
1584 risec = std::sqrt(t3)*secRMax ;
1585 *validNorm = true ;
1587 }
1588 return snxt=0 ;
1589 }
1590 else
1591 {
1592 sider = kRMax ;
1593 if (b>0) { srd = -b - std::sqrt(d); }
1594 else { srd = c/(-b+std::sqrt(d)) ; }
1595
1597 ri = tanRMax*zi + rMaxAv ;
1598
1599 if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance))
1600 {
1601 // An intersection within the tolerance
1602 // we will Store it in case it is good -
1603 //
1604 slentol = srd ;
1605 sidetol = kRMax ;
1606 }
1607 if ( (ri < 0) || (srd < halfRadTolerance) )
1608 {
1609 // Safety: if both roots -ve ensure that srd cannot `win'
1610 // distance to out
1611
1612 if (b>0) { sr2 = c/(-b-std::sqrt(d)); }
1613 else { sr2 = -b + std::sqrt(d); }
1615 ri = tanRMax*zi + rMaxAv ;
1616
1617 if ((ri >= 0) && (sr2 > halfRadTolerance))
1618 {
1619 srd = sr2;
1620 }
1621 else
1622 {
1623 srd = kInfinity ;
1624
1625 if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) )
1626 {
1627 // An intersection within the tolerance.
1628 // Storing it in case it is good.
1629
1630 slentol = sr2 ;
1631 sidetol = kRMax ;
1632 }
1633 }
1634 }
1635 }
1636 }
1637 else
1638 {
1639 // No intersection with outer cone & not parallel
1640 // -> already outside, no intersection
1641
1642 if ( calcNorm )
1643 {
1644 risec = std::sqrt(t3)*secRMax;
1645 *validNorm = true;
1647 }
1648 return snxt = 0.0 ;
1649 }
1650 }
1651 else if ( nt2 && (deltaRoi2 > 0.0) )
1652 {
1653 // Linear case (only one intersection) => point outside outer cone
1654
1655 if ( calcNorm )
1656 {
1657 risec = std::sqrt(t3)*secRMax;
1658 *validNorm = true;
1660 }
1661 return snxt = 0.0 ;
1662 }
1663 else
1664 {
1665 // No intersection -> parallel to outer cone
1666 // => Z or inner cone intersection
1667
1668 srd = kInfinity ;
1669 }
1670
1671 // Check possible intersection within tolerance
1672
1673 if ( slentol <= halfCarTolerance )
1674 {
1675 // An intersection within the tolerance was found.
1676 // We must accept it only if the momentum points outwards.
1677 //
1678 // G4ThreeVector ptTol ; // The point of the intersection
1679 // ptTol= p + slentol*v ;
1680 // ri=tanRMax*zi+rMaxAv ;
1681 //
1682 // Calculate a normal vector, as below
1683
1686 risec = std::sqrt(xi*xi + yi*yi)*secRMax;
1688
1690 {
1691 if ( calcNorm )
1692 {
1694 *validNorm = true ;
1695 }
1696 return snxt = 0.0 ;
1697 }
1698 else // On the surface, but not heading out so we ignore this intersection
1699 { // (as it is within tolerance).
1700 slentol = kInfinity ;
1701 }
1702 }
1703
1704 // Inner Cone intersection
1705
1706 if ( fRmin1 || fRmin2 )
1707 {
1708 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
1710
1711 if ( nt1 )
1712 {
1713 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
1714 rMinAv = (fRmin1 + fRmin2)*0.5 ;
1715 rin = tanRMin*p.z() + rMinAv ;
1716 nt2 = t2 - tanRMin*v.z()*rin ;
1717 nt3 = t3 - rin*rin ;
1718
1719 // Equation quadratic => 2 roots : first root must be leaving
1720
1721 b = nt2/nt1 ;
1722 c = nt3/nt1 ;
1723 d = b*b - c ;
1724
1725 if ( d >= 0.0 )
1726 {
1727 // NOTE: should be rho-rin<kRadTolerance*0.5,
1728 // but using squared versions for efficiency
1729
1730 if (nt3 < kRadTolerance*(rin + kRadTolerance*0.25))
1731 {
1732 if ( nt2 < 0.0 )
1733 {
1734 if (calcNorm) { *validNorm = false; }
1735 return snxt = 0.0;
1736 }
1737 }
1738 else
1739 {
1740 if (b>0) { sr2 = -b - std::sqrt(d); }
1741 else { sr2 = c/(-b+std::sqrt(d)); }
1743 ri = tanRMin*zi + rMinAv ;
1744
1745 if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) )
1746 {
1747 // An intersection within the tolerance
1748 // storing it in case it is good.
1749
1750 slentol = sr2 ;
1751 sidetol = kRMax ;
1752 }
1753 if( (ri<0) || (sr2 < halfRadTolerance) )
1754 {
1755 if (b>0) { sr3 = c/(-b-std::sqrt(d)); }
1756 else { sr3 = -b + std::sqrt(d) ; }
1757
1758 // Safety: if both roots -ve ensure that srd cannot `win'
1759 // distancetoout
1760
1761 if ( sr3 > halfRadTolerance )
1762 {
1763 if( sr3 < srd )
1764 {
1766 ri = tanRMin*zi + rMinAv ;
1767
1768 if ( ri >= 0.0 )
1769 {
1770 srd=sr3 ;
1771 sider=kRMin ;
1772 }
1773 }
1774 }
1775 else if ( sr3 > -halfRadTolerance )
1776 {
1777 // Intersection in tolerance. Store to check if it's good
1778
1779 slentol = sr3 ;
1780 sidetol = kRMin ;
1781 }
1782 }
1783 else if ( (sr2 < srd) && (sr2 > halfCarTolerance) )
1784 {
1785 srd = sr2 ;
1786 sider = kRMin ;
1787 }
1788 else if (sr2 > -halfCarTolerance)
1789 {
1790 // Intersection in tolerance. Store to check if it's good
1791
1792 slentol = sr2 ;
1793 sidetol = kRMin ;
1794 }
1795 if( slentol <= halfCarTolerance )
1796 {
1797 // An intersection within the tolerance was found.
1798 // We must accept it only if the momentum points outwards.
1799
1800 G4ThreeVector Normal ;
1801
1802 // Calculate a normal vector, as below
1803
1806 if( sidetol==kRMax )
1807 {
1808 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1809 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
1810 }
1811 else
1812 {
1813 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
1814 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
1815 }
1817 {
1818 // We will leave the cone immediately
1819
1820 if( calcNorm )
1821 {
1823 *validNorm = true ;
1824 }
1825 return snxt = 0.0 ;
1826 }
1827 else
1828 {
1829 // On the surface, but not heading out so we ignore this
1830 // intersection (as it is within tolerance).
1831
1832 slentol = kInfinity ;
1833 }
1834 }
1835 }
1836 }
1837 }
1838 }
1839
1840 // Linear case => point outside inner cone ---> outer cone intersect
1841 //
1842 // Phi Intersection
1843
1844 if ( !fPhiFullCone )
1845 {
1846 // add angle calculation with correction
1847 // of the difference in domain of atan2 and Sphi
1848
1850
1851 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
1852 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
1853
1855 {
1856 // pDist -ve when inside
1857
1860
1861 // Comp -ve when in direction of outwards normal
1862
1865
1866 sidephi = kNull ;
1867
1868 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
1869 && (pDistE <= halfCarTolerance) ) )
1870 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
1871 && (pDistE > halfCarTolerance) ) ) )
1872 {
1873 // Inside both phi *full* planes
1874 if ( compS < 0 )
1875 {
1876 sphi = pDistS/compS ;
1877 if (sphi >= -halfCarTolerance)
1878 {
1881
1882 // Check intersecting with correct half-plane
1883 // (if not -> no intersect)
1884 //
1886 && (std::fabs(yi)<=kCarTolerance) )
1887 {
1888 sidephi= kSPhi;
1889 if ( ( fSPhi-halfAngTolerance <= vphi )
1890 && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) )
1891 {
1892 sphi = kInfinity;
1893 }
1894 }
1895 else
1896 if ( (yi*cosCPhi-xi*sinCPhi)>=0 )
1897 {
1898 sphi = kInfinity ;
1899 }
1900 else
1901 {
1902 sidephi = kSPhi ;
1903 if ( pDistS > -halfCarTolerance )
1904 {
1905 sphi = 0.0 ; // Leave by sphi immediately
1906 }
1907 }
1908 }
1909 else
1910 {
1911 sphi = kInfinity ;
1912 }
1913 }
1914 else
1915 {
1916 sphi = kInfinity ;
1917 }
1918
1919 if ( compE < 0 )
1920 {
1921 sphi2 = pDistE/compE ;
1922
1923 // Only check further if < starting phi intersection
1924 //
1925 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
1926 {
1929
1930 // Check intersecting with correct half-plane
1931
1933 && (std::fabs(yi)<=kCarTolerance) )
1934 {
1935 // Leaving via ending phi
1936
1937 if(!( (fSPhi-halfAngTolerance <= vphi)
1938 && (fSPhi+fDPhi+halfAngTolerance >= vphi) ) )
1939 {
1940 sidephi = kEPhi ;
1941 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1942 else { sphi = 0.0; }
1943 }
1944 }
1945 else // Check intersecting with correct half-plane
1946 if ( yi*cosCPhi-xi*sinCPhi >= 0 )
1947 {
1948 // Leaving via ending phi
1949
1950 sidephi = kEPhi ;
1951 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
1952 else { sphi = 0.0; }
1953 }
1954 }
1955 }
1956 }
1957 else
1958 {
1959 sphi = kInfinity ;
1960 }
1961 }
1962 else
1963 {
1964 // On z axis + travel not || to z axis -> if phi of vector direction
1965 // within phi of shape, Step limited by rmax, else Step =0
1966
1967 if ( (fSPhi-halfAngTolerance <= vphi)
1968 && (vphi <= fSPhi+fDPhi+halfAngTolerance) )
1969 {
1970 sphi = kInfinity ;
1971 }
1972 else
1973 {
1974 sidephi = kSPhi ; // arbitrary
1975 sphi = 0.0 ;
1976 }
1977 }
1978 if ( sphi < snxt ) // Order intersecttions
1979 {
1980 snxt=sphi ;
1981 side=sidephi ;
1982 }
1983 }
1984 if ( srd < snxt ) // Order intersections
1985 {
1986 snxt = srd ;
1987 side = sider ;
1988 }
1989 if (calcNorm)
1990 {
1991 switch(side)
1992 { // Note: returned vector not normalised
1993 case kRMax: // (divide by frmax for unit vector)
1996 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
1998 *validNorm = true ;
1999 break ;
2000 case kRMin:
2001 *validNorm = false ; // Rmin is inconvex
2002 break ;
2003 case kSPhi:
2004 if ( fDPhi <= pi )
2005 {
2007 *validNorm = true ;
2008 }
2009 else
2010 {
2011 *validNorm = false ;
2012 }
2013 break ;
2014 case kEPhi:
2015 if ( fDPhi <= pi )
2016 {
2018 *validNorm = true ;
2019 }
2020 else
2021 {
2022 *validNorm = false ;
2023 }
2024 break ;
2025 case kPZ:
2027 *validNorm = true ;
2028 break ;
2029 case kMZ:
2031 *validNorm = true ;
2032 break ;
2033 default:
2035 DumpInfo();
2036 std::ostringstream message;
2037 G4int oldprc = message.precision(16) ;
2038 message << "Undefined side for valid surface normal to solid."
2039 << G4endl
2047 {
2050 }
2057 message.precision(oldprc) ;
2059 JustWarning, message) ;
2060 break ;
2061 }
2062 }
2063 if (snxt < halfCarTolerance) { snxt = 0.; }
2064
2065 return snxt ;
2066}
Hep3Vector unit() const
◆ GetCubicVolume()
◆ GetDeltaPhiAngle()
|
inline |
Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsPhi::G4ParameterisationConsPhi(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsPhi::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().
◆ GetDPhi()
|
inline |
Referenced by G4tgbVolume::BuildSolidForDivision().
◆ GetDz()
|
inline |
◆ GetEntityType()
|
virtual |
◆ GetInnerRadiusMinusZ()
|
inline |
Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), G4tgbVolume::BuildSolidForDivision(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsRho::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().
◆ GetInnerRadiusPlusZ()
|
inline |
Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), G4tgbVolume::BuildSolidForDivision(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), and G4tgbGeometryDumper::GetSolidParams().
◆ GetOuterRadiusMinusZ()
|
inline |
Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), G4tgbVolume::BuildSolidForDivision(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsRho::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().
◆ GetOuterRadiusPlusZ()
|
inline |
Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), G4tgbVolume::BuildSolidForDivision(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4VParameterisationCons::G4VParameterisationCons(), and G4tgbGeometryDumper::GetSolidParams().
◆ GetPointOnSurface()
|
virtual |
Reimplemented from G4VSolid.
Definition at line 2289 of file G4Cons.cc.
2290{
2291 // declare working variables
2292 //
2293 G4double Aone, Atwo, Athree, Afour, Afive, slin, slout, phi;
2294 G4double zRand, cosu, sinu, rRand1, rRand2, chose, rone, rtwo, qone, qtwo;
2295 rone = (fRmax1-fRmax2)/(2.*fDz);
2296 rtwo = (fRmin1-fRmin2)/(2.*fDz);
2297 qone=0.; qtwo=0.;
2298 if(fRmax1!=fRmax2) { qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2); }
2299 if(fRmin1!=fRmin2) { qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2); }
2302 Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*slout;
2303 Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*slin;
2304 Athree = 0.5*fDPhi*(fRmax1*fRmax1-fRmin1*fRmin1);
2305 Afour = 0.5*fDPhi*(fRmax2*fRmax2-fRmin2*fRmin2);
2306 Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
2307
2308 phi = RandFlat::shoot(fSPhi,fSPhi+fDPhi);
2309 cosu = std::cos(phi); sinu = std::sin(phi);
2310 rRand1 = GetRadiusInRing(fRmin1, fRmin2);
2311 rRand2 = GetRadiusInRing(fRmax1, fRmax2);
2312
2313 if ( (fSPhi == 0.) && fPhiFullCone ) { Afive = 0.; }
2314 chose = RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
2315
2316 if( (chose >= 0.) && (chose < Aone) )
2317 {
2318 if(fRmin1 != fRmin2)
2319 {
2320 zRand = RandFlat::shoot(-1.*fDz,fDz);
2322 rtwo*sinu*(qtwo-zRand), zRand);
2323 }
2324 else
2325 {
2327 RandFlat::shoot(-1.*fDz,fDz));
2328 }
2329 }
2330 else if( (chose >= Aone) && (chose <= Aone + Atwo) )
2331 {
2332 if(fRmax1 != fRmax2)
2333 {
2334 zRand = RandFlat::shoot(-1.*fDz,fDz);
2336 rone*sinu*(qone-zRand), zRand);
2337 }
2338 else
2339 {
2341 RandFlat::shoot(-1.*fDz,fDz));
2342 }
2343 }
2344 else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) )
2345 {
2347 }
2348 else if( (chose >= Aone + Atwo + Athree)
2349 && (chose < Aone + Atwo + Athree + Afour) )
2350 {
2352 }
2353 else if( (chose >= Aone + Atwo + Athree + Afour)
2354 && (chose < Aone + Atwo + Athree + Afour + Afive) )
2355 {
2356 zRand = RandFlat::shoot(-1.*fDz,fDz);
2357 rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2358 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2360 rRand1*std::sin(fSPhi), zRand);
2361 }
2362 else
2363 {
2364 zRand = RandFlat::shoot(-1.*fDz,fDz);
2365 rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
2366 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
2368 rRand1*std::sin(fSPhi+fDPhi), zRand);
2369 }
2370}
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
Definition: G4CSGSolid.cc:101
◆ GetRmax1()
|
inline |
◆ GetRmax2()
|
inline |
◆ GetRmin1()
|
inline |
◆ GetRmin2()
|
inline |
◆ GetSPhi()
|
inline |
Referenced by G4tgbVolume::BuildSolidForDivision().
◆ GetStartPhiAngle()
|
inline |
Referenced by G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4VParameterisationCons::G4VParameterisationCons(), and G4tgbGeometryDumper::GetSolidParams().
◆ GetSurfaceArea()
◆ GetZHalfLength()
|
inline |
Referenced by G4HepRepFileSceneHandler::AddSolid(), G4HepRepSceneHandler::AddSolid(), G4tgbVolume::BuildSolidForDivision(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsZ::ComputeTransformation(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsZ::G4ParameterisationConsZ(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsZ::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().
◆ Inside()
|
virtual |
Implements G4VSolid.
Definition at line 191 of file G4Cons.cc.
192{
194 EInside in;
198
202
206
207 // rh2 = rh*rh;
208
209 tolRMin = rl - halfRadTolerance;
210 if ( tolRMin < 0 ) { tolRMin = 0; }
211 tolRMax = rh + halfRadTolerance;
212
214
215 if (rl) { tolRMin = rl + halfRadTolerance; }
216 else { tolRMin = 0.0; }
217 tolRMax = rh - halfRadTolerance;
218
220 {
222 }
224 {
226
227 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
228 else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
229
230 if ( (pPhi < fSPhi - halfAngTolerance) ||
232
234 {
235 if ( (pPhi < fSPhi + halfAngTolerance) ||
236 (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; }
237 }
238 }
240
241 return in ;
242}
Referenced by CalculateExtent(), and DistanceToOut().
◆ operator=()
Definition at line 161 of file G4Cons.cc.
162{
163 // Check assignment to self
164 //
165 if (this == &rhs) { return *this; }
166
167 // Copy base class data
168 //
169 G4CSGSolid::operator=(rhs);
170
171 // Copy data
172 //
173 kRadTolerance = rhs.kRadTolerance;
174 kAngTolerance = rhs.kAngTolerance;
175 fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
176 fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
177 fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
178 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
179 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
180 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
181 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
182 fPhiFullCone = rhs.fPhiFullCone;
183
184 return *this;
185}
◆ SetDeltaPhiAngle()
|
inline |
◆ SetInnerRadiusMinusZ()
|
inline |
◆ SetInnerRadiusPlusZ()
|
inline |
◆ SetOuterRadiusMinusZ()
|
inline |
◆ SetOuterRadiusPlusZ()
|
inline |
◆ SetStartPhiAngle()
◆ SetZHalfLength()
|
inline |
◆ StreamInfo()
|
virtual |
Reimplemented from G4CSGSolid.
Definition at line 2262 of file G4Cons.cc.
2263{
2264 G4int oldprc = os.precision(16);
2265 os << "-----------------------------------------------------------\n"
2267 << " ===================================================\n"
2268 << " Solid type: G4Cons\n"
2269 << " Parameters: \n"
2270 << " inside -fDz radius: " << fRmin1/mm << " mm \n"
2271 << " outside -fDz radius: " << fRmax1/mm << " mm \n"
2272 << " inside +fDz radius: " << fRmin2/mm << " mm \n"
2273 << " outside +fDz radius: " << fRmax2/mm << " mm \n"
2274 << " half length in Z : " << fDz/mm << " mm \n"
2275 << " starting angle of segment: " << fSPhi/degree << " degrees \n"
2276 << " delta angle of segment : " << fDPhi/degree << " degrees \n"
2277 << "-----------------------------------------------------------\n";
2278 os.precision(oldprc);
2279
2280 return os;
2281}
◆ SurfaceNormal()
|
virtual |
Implements G4VSolid.
Definition at line 479 of file G4Cons.cc.
480{
481 G4int noSurfaces = 0;
482 G4double rho, pPhi;
483 G4double distZ, distRMin, distRMax;
484 G4double distSPhi = kInfinity, distEPhi = kInfinity;
485 G4double tanRMin, secRMin, pRMin, widRMin;
486 G4double tanRMax, secRMax, pRMax, widRMax;
487
490
492 G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;
493
494 distZ = std::fabs(std::fabs(p.z()) - fDz);
496
497 tanRMin = (fRmin2 - fRmin1)*0.5/fDz;
498 secRMin = std::sqrt(1 + tanRMin*tanRMin);
499 pRMin = rho - p.z()*tanRMin;
500 widRMin = fRmin2 - fDz*tanRMin;
501 distRMin = std::fabs(pRMin - widRMin)/secRMin;
502
503 tanRMax = (fRmax2 - fRmax1)*0.5/fDz;
504 secRMax = std::sqrt(1+tanRMax*tanRMax);
505 pRMax = rho - p.z()*tanRMax;
506 widRMax = fRmax2 - fDz*tanRMax;
507 distRMax = std::fabs(pRMax - widRMax)/secRMax;
508
509 if (!fPhiFullCone) // Protected against (0,0,z)
510 {
511 if ( rho )
512 {
514
515 if (pPhi < fSPhi-delta) { pPhi += twopi; }
516 else if (pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; }
517
518 distSPhi = std::fabs( pPhi - fSPhi );
519 distEPhi = std::fabs( pPhi - fSPhi - fDPhi );
520 }
521 else if( !(fRmin1) || !(fRmin2) )
522 {
523 distSPhi = 0.;
524 distEPhi = 0.;
525 }
526 nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0);
527 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0);
528 }
529 if ( rho > delta )
530 {
532 if (fRmin1 || fRmin2)
533 {
535 }
536 }
537
538 if( distRMax <= delta )
539 {
540 noSurfaces ++;
541 sumnorm += nR;
542 }
543 if( (fRmin1 || fRmin2) && (distRMin <= delta) )
544 {
545 noSurfaces ++;
546 sumnorm += nr;
547 }
548 if( !fPhiFullCone )
549 {
550 if (distSPhi <= dAngle)
551 {
552 noSurfaces ++;
553 sumnorm += nPs;
554 }
555 if (distEPhi <= dAngle)
556 {
557 noSurfaces ++;
558 sumnorm += nPe;
559 }
560 }
561 if (distZ <= delta)
562 {
563 noSurfaces ++;
565 else { sumnorm -= nZ; }
566 }
567 if ( noSurfaces == 0 )
568 {
569#ifdef G4CSGDEBUG
572#endif
573 norm = ApproxSurfaceNormal(p);
574 }
575 else if ( noSurfaces == 1 ) { norm = sumnorm; }
577
578 return norm ;
579}
The documentation for this class was generated from the following files:
