VecGeom/volumes/RadialConeStruct.h

0001 /// @file RadialConeStruct.h
0002 /// @brief Radial and z data for complete cone sections without phi state.
0003
0004 #ifndef VECGEOM_VOLUMES_RADIALCONESTRUCT_H_
0005 #define VECGEOM_VOLUMES_RADIALCONESTRUCT_H_
0006
0007 #include "VecGeom/base/Global.h"
0008
0009 namespace vecgeom {
0010
0011 VECGEOM_DEVICE_DECLARE_CONV_TEMPLATE(struct, RadialConeStruct, typename);
0012
0013 inline namespace VECGEOM_IMPL_NAMESPACE {
0014
0015 template <typename T = double>
0016 struct RadialConeStruct {
0017   T fRmin1{0.};
0018   T fRmax1{0.};
0019   T fRmin2{0.};
0020   T fRmax2{0.};
0021   T fDz{0.};
0022
0023   Precision fInnerSlope{0.};
0024   Precision fOuterSlope{0.};
0025   Precision fInnerOffset{0.};
0026   Precision fOuterOffset{0.};
0027   Precision fInnerTolerance{0.};
0028   Precision fOuterTolerance{0.};
0029   Precision fSqRmin1{0.}, fSqRmin2{0.};
0030   Precision fSqRmax1{0.}, fSqRmax2{0.};
0031   Precision fTolIz{0.}, fTolOz{0.};
0032   Precision fInnerConeApex{0.};
0033   Precision fTanInnerApexAngle{0.};
0034   Precision fOuterConeApex{0.};
0035   Precision fTanOuterApexAngle{0.};
0036   Precision fSecRMin{0.};
0037   Precision fSecRMax{0.};
0038   Precision fInvSecRMin{0.};
0039   Precision fInvSecRMax{0.};
0040   Precision fTanRMin{0.};
0041   Precision fTanRMax{0.};
0042   Precision fZNormInner{0.};
0043   Precision fZNormOuter{0.};
0044
0045   VECCORE_ATT_HOST_DEVICE
0046   Precision Capacity(Precision dphi) const
0047   {
0048     return (fDz * dphi / 3.) *
0049            (fRmax1 * fRmax1 + fRmax2 * fRmax2 + fRmax1 * fRmax2 - fRmin1 * fRmin1 - fRmin2 * fRmin2 - fRmin1 * fRmin2);
0050   }
0051
0052   VECCORE_ATT_HOST_DEVICE
0053   void CalculateCached()
0054   {
0055
0056     fSqRmin1 = fRmin1 * fRmin1;
0057     fSqRmax1 = fRmax1 * fRmax1;
0058     fSqRmin2 = fRmin2 * fRmin2;
0059     fSqRmax2 = fRmax2 * fRmax2;
0060
0061     fTanRMin    = (fRmin2 - fRmin1) * 0.5 / fDz;
0062     fSecRMin    = Sqrt(1.0 + fTanRMin * fTanRMin);
0063     fInvSecRMin = 1. / NonZero(fSecRMin);
0064     fTanRMax    = (fRmax2 - fRmax1) * 0.5 / fDz;
0065     fSecRMax    = Sqrt(1.0 + fTanRMax * fTanRMax);
0066     fInvSecRMax = 1. / NonZero(fSecRMax);
0067
0068     fInnerSlope     = -(fRmin1 - fRmin2) / (2. * fDz);
0069     fOuterSlope     = -(fRmax1 - fRmax2) / (2. * fDz);
0070     fInnerOffset    = fRmin2 - fInnerSlope * fDz;
0071     fOuterOffset    = fRmax2 - fOuterSlope * fDz;
0072     fInnerTolerance = kConeTolerance * fSecRMin;
0073     fOuterTolerance = kConeTolerance * fSecRMax;
0074
0075     if (fRmin2 > fRmin1) {
0076       fInnerConeApex     = 2 * fDz * fRmin1 / (fRmin2 - fRmin1);
0077       fTanInnerApexAngle = fRmin2 / (2 * fDz + fInnerConeApex);
0078     } else {
0079       fInnerConeApex     = 2 * fDz * fRmin2 / NonZero(fRmin1 - fRmin2);
0080       fTanInnerApexAngle = fRmin1 / (2 * fDz + fInnerConeApex);
0081     }
0082     if (fRmin1 == 0. || fRmin2 == 0.) fInnerConeApex = 0.;
0083     if (fRmin1 == 0.) fTanInnerApexAngle = fRmin2 / (2 * fDz);
0084     if (fRmin2 == 0.) fTanInnerApexAngle = fRmin1 / (2 * fDz);
0085
0086     if (fRmax2 > fRmax1) {
0087       fOuterConeApex     = 2 * fDz * fRmax1 / (fRmax2 - fRmax1);
0088       fTanOuterApexAngle = fRmax2 / (2 * fDz + fOuterConeApex);
0089     } else {
0090       fOuterConeApex     = 2 * fDz * fRmax2 / NonZero(fRmax1 - fRmax2);
0091       fTanOuterApexAngle = fRmax1 / (2 * fDz + fOuterConeApex);
0092     }
0093     if (fRmax1 == 0. || fRmax2 == 0.) fOuterConeApex = 0.;
0094     if (fRmax1 == 0.) fTanOuterApexAngle = fRmax2 / (2 * fDz);
0095     if (fRmax2 == 0.) fTanOuterApexAngle = fRmax1 / (2 * fDz);
0096
0097     fZNormInner = fTanRMin / NonZero(fSecRMin);
0098     fZNormOuter = -fTanRMax / NonZero(fSecRMax);
0099     fTolIz      = fDz - kHalfTolerance;
0100     fTolOz      = fDz + kHalfTolerance;
0101   }
0102
0103   VECCORE_ATT_HOST_DEVICE
0104   void Init(T rmin1, T rmax1, T rmin2, T rmax2, T dz)
0105   {
0106     fRmin1 = rmin1 < 0. ? 0. : rmin1;
0107     fRmax1 = rmax1;
0108     fRmin2 = rmin2 < 0. ? 0. : rmin2;
0109     fRmax2 = rmax2;
0110     fDz    = dz;
0111     CalculateCached();
0112   }
0113
0114   /// @brief Accumulate the outward normals of all incident radial and z surfaces.
0115   /// @details The returned vector is an unnormalized sum so callers can append
0116   /// normals from surfaces represented outside this radial section.
0117   /// @param p Point in the section's local coordinates.
0118   /// @param[out] norm Raw incident-surface sum, or the nearest-surface direction
0119   /// when the point is not on a surface.
0120   /// @param[out] count Number of incident surfaces included in `norm`; zero when
0121   /// the nearest-surface fallback is returned.
0122   /// @return True for one or two incident surfaces; false for no surface or more
0123   /// than two incident surfaces.
0124   VECCORE_ATT_HOST_DEVICE
0125   bool NormalWithCount(Vector3D<Precision> const &p, Vector3D<Precision> &norm, int &count) const
0126   {
0127     const Precision rho = Sqrt(p.Perp2());
0128     Precision distZ     = Abs(Abs(p.z()) - fDz);
0129     Precision distRMin;
0130     Precision distRMax;
0131
0132     const Precision pRMin   = rho - p.z() * fTanRMin;
0133     const Precision widRMin = fRmin2 - fDz * fTanRMin;
0134     if (Abs(fRmin1 - fRmin2) < fInnerTolerance)
0135       distRMin = Abs(rho - fRmin2);
0136     else
0137       distRMin = Abs((pRMin - widRMin) / fSecRMin);
0138
0139     const Precision pRMax   = rho - p.z() * fTanRMax;
0140     const Precision widRMax = fRmax2 - fDz * fTanRMax;
0141     if (Abs(fRmax1 - fRmax2) < fOuterTolerance)
0142       distRMax = Abs(rho - fRmax2);
0143     else
0144       distRMax = Abs((pRMax - widRMax) / fSecRMax);
0145
0146     Vector3D<Precision> radialOuter(0., 0., fZNormOuter);
0147     Vector3D<Precision> radialInner(0., 0., fZNormInner);
0148     if (rho > kHalfTolerance) {
0149       radialOuter.Set(p.x() / rho / fSecRMax, p.y() / rho / fSecRMax, -fTanRMax / fSecRMax);
0150       radialInner.Set(-p.x() / rho / fSecRMin, -p.y() / rho / fSecRMin, fTanRMin / fSecRMin);
0151     }
0152
0153     const bool radiallyInsideOuter = rho <= GetRadiusAtZ<false>(p.z()) + fOuterTolerance;
0154     const bool radiallyInsideInner =
0155         (fRmin1 == 0. && fRmin2 == 0.) || rho >= GetRadiusAtZ<true>(p.z()) - fInnerTolerance;
0156     const bool insideRadial = radiallyInsideOuter && radiallyInsideInner;
0157
0158     count = 0;
0159     Vector3D<Precision> sum(0., 0., 0.);
0160     if (insideRadial && distZ <= kHalfTolerance) {
0161       const bool upper      = p.z() >= 0.;
0162       const bool degenerate = upper ? Abs(fRmax2 - fRmin2) < kTolerance : Abs(fRmax1 - fRmin1) < kTolerance;
0163       if (!degenerate) {
0164         sum += Vector3D<Precision>(0., 0., upper ? 1. : -1.);
0165         ++count;
0166       }
0167     }
0168     if (Abs(p.z()) <= fDz + kTolerance && distRMax <= fOuterTolerance) {
0169       sum += radialOuter;
0170       ++count;
0171     }
0172     if ((fRmin1 || fRmin2) && Abs(p.z()) <= fDz + kTolerance && distRMin <= fInnerTolerance) {
0173       sum += radialInner;
0174       ++count;
0175     }
0176
0177     if (count) {
0178       norm = sum;
0179       return count <= 2;
0180     }
0181     norm = distZ < distRMax ? Vector3D<Precision>(0., 0., p.z() < 0. ? -1. : 1.) : radialOuter;
0182     if ((fRmin1 || fRmin2) && distRMin < Min(distZ, distRMax)) norm = radialInner;
0183     return false;
0184   }
0185
0186   VECCORE_ATT_HOST_DEVICE
0187   bool Normal(Vector3D<Precision> const &p, Vector3D<Precision> &norm) const
0188   {
0189     int count        = 0;
0190     const bool valid = NormalWithCount(p, norm, count);
0191     if (count > 1) norm = norm.Unit();
0192     return valid;
0193   }
0194
0195   template <bool Inner>
0196   VECCORE_ATT_HOST_DEVICE T GetRadiusAtZ(T z) const
0197   {
0198     return Inner ? T(fInnerSlope * z + fInnerOffset) : T(fOuterSlope * z + fOuterOffset);
0199   }
0200
0201   VECCORE_ATT_HOST_DEVICE
0202   void Print() const
0203   {
0204     printf("RadialConeStruct : {rmin1 %.2f, rmax1 %.2f, rmin2 %.2f, rmax2 %.2f, dz %.2f}", fRmin1, fRmax1, fRmin2,
0205            fRmax2, fDz);
0206   }
0207 };
0208
0209 } // namespace VECGEOM_IMPL_NAMESPACE
0210 } // namespace vecgeom
0211
0212 #endif