VecGeom/volumes/ConeStruct.h

0001 /*
0002  * ConeStruct.h
0003  *
0004  *  Created on: May 09, 2017
0005  *      Author: Raman Sehgal
0006  */
0007 #ifndef VECGEOM_CONESTRUCT_H_
0008 #define VECGEOM_CONESTRUCT_H_
0009
0010 #include "VecGeom/base/Global.h"
0011 #include "VecGeom/volumes/Wedge_Evolution.h"
0012 #include <VecGeom/management/Logger.h>
0013
0014 namespace vecgeom {
0015
0016 VECGEOM_DEVICE_DECLARE_CONV_TEMPLATE(struct, ConeStruct, typename);
0017
0018 inline namespace VECGEOM_IMPL_NAMESPACE {
0019
0020 // a plain and lightweight struct to encapsulate data members of a Cone
0021 template <typename T = double>
0022 struct ConeStruct {
0023   // Cone defining parameters
0024   T fRmin1{0.};
0025   T fRmax1{0.};
0026   T fRmin2{0.};
0027   T fRmax2{0.};
0028   T fDz{0.};
0029   T fSPhi{0.};
0030   T fDPhi{0.};
0031
0032   /* These new data members are introduced to store the original paramters of
0033    * Cone, which may change in the case where rmin is equal to rmax.
0034    * These are basically required by the Extent functions to do more accurate
0035    * bounding box calculations.
0036    */
0037   T _frmin1{0.};
0038   T _frmin2{0.};
0039   T _frmax1{0.};
0040   T _frmax2{0.};
0041
0042   evolution::Wedge fPhiWedge;
0043
0044   // vectors characterizing the normals of phi planes
0045   // makes task to detect phi sektors very efficient
0046   Vector3D<Precision> fNormalPhi1;
0047   Vector3D<Precision> fNormalPhi2;
0048   Precision fAlongPhi1x{0.};
0049   Precision fAlongPhi1y{0.};
0050   Precision fAlongPhi2x{0.};
0051   Precision fAlongPhi2y{0.};
0052
0053   // Some Cached value, try to reduce them
0054   // Some precomputed values to avoid divisions etc
0055   Precision fInnerSlope{0.}; // "gradient" of inner surface in z direction
0056   Precision fOuterSlope{0.}; // "gradient" of outer surface in z direction
0057   Precision fInnerOffset{0.};
0058   Precision fOuterOffset{0.};
0059   Precision fInnerTolerance{0.}; // tolerance on radial direction for inner surface
0060   Precision fOuterTolerance{0.}; // tolerance on radial direction for outer surface
0061   // Values to be cached
0062   Precision fSqRmin1{0.}, fSqRmin2{0.};
0063   Precision fSqRmax1{0.}, fSqRmax2{0.};
0064   Precision fTolIz{0.}, fTolOz{0.};
0065   Precision fInnerConeApex{0.};
0066   Precision fTanInnerApexAngle{0.};
0067   Precision fOuterConeApex{0.};
0068   Precision fTanOuterApexAngle{0.};
0069
0070   Precision fSecRMin{0.};
0071   Precision fSecRMax{0.};
0072   Precision fInvSecRMin{0.};
0073   Precision fInvSecRMax{0.};
0074   Precision fTanRMin{0.};
0075   Precision fTanRMax{0.};
0076   Precision fZNormInner{0.};
0077   Precision fZNormOuter{0.};
0078
0079   /* Some additional variable to store original Rmax
0080    * for the cases when Rmax is modified because of Rmin==Rmax
0081    */
0082   Precision fOriginalRmax1{0.};
0083   Precision fOriginalRmax2{0.};
0084
0085   VECCORE_ATT_HOST_DEVICE
0086   Precision Capacity() const
0087   {
0088     return (fDz * fDPhi / 3.) *
0089            (fRmax1 * fRmax1 + fRmax2 * fRmax2 + fRmax1 * fRmax2 - fRmin1 * fRmin1 - fRmin2 * fRmin2 - fRmin1 * fRmin2);
0090   }
0091
0092   VECCORE_ATT_HOST_DEVICE
0093   void CalculateCached()
0094   {
0095     fOriginalRmax1 = fRmax1;
0096     fOriginalRmax2 = fRmax2;
0097
0098     if (fRmin1 == fRmax1) {
0099       fRmax1 += kConeTolerance;
0100     }
0101     if (fRmin2 == fRmax2) {
0102       fRmax2 += kConeTolerance;
0103     }
0104
0105     fSqRmin1 = fRmin1 * fRmin1;
0106     fSqRmax1 = fRmax1 * fRmax1;
0107     fSqRmin2 = fRmin2 * fRmin2;
0108     fSqRmax2 = fRmax2 * fRmax2;
0109
0110     fTanRMin    = (fRmin2 - fRmin1) * 0.5 / fDz;
0111     fSecRMin    = std::sqrt(1.0 + fTanRMin * fTanRMin);
0112     fInvSecRMin = 1. / NonZero(fSecRMin);
0113     fTanRMax    = (fRmax2 - fRmax1) * 0.5 / fDz;
0114
0115     fSecRMax    = std::sqrt(1.0 + fTanRMax * fTanRMax);
0116     fInvSecRMax = 1. / NonZero(fSecRMax);
0117
0118     // check this very carefully
0119     fInnerSlope     = -(fRmin1 - fRmin2) / (2. * fDz);
0120     fOuterSlope     = -(fRmax1 - fRmax2) / (2. * fDz);
0121     fInnerOffset    = fRmin2 - fInnerSlope * fDz;
0122     fOuterOffset    = fRmax2 - fOuterSlope * fDz;
0123     fInnerTolerance = kConeTolerance * fSecRMin;
0124     fOuterTolerance = kConeTolerance * fSecRMax;
0125
0126     if (fRmin2 > fRmin1) {
0127       fInnerConeApex     = 2 * fDz * fRmin1 / (fRmin2 - fRmin1);
0128       fTanInnerApexAngle = fRmin2 / (2 * fDz + fInnerConeApex);
0129     } else { // Should we add a check if(fRmin1 > fRmin2)
0130       fInnerConeApex     = 2 * fDz * fRmin2 / NonZero(fRmin1 - fRmin2);
0131       fTanInnerApexAngle = fRmin1 / (2 * fDz + fInnerConeApex);
0132     }
0133
0134     if (fRmin1 == 0. || fRmin2 == 0.) fInnerConeApex = 0.;
0135
0136     if (fRmin1 == 0.) fTanInnerApexAngle = fRmin2 / (2 * fDz);
0137     if (fRmin2 == 0.) fTanInnerApexAngle = fRmin1 / (2 * fDz);
0138
0139     if (fRmax2 > fRmax1) {
0140       fOuterConeApex     = 2 * fDz * fRmax1 / (fRmax2 - fRmax1);
0141       fTanOuterApexAngle = fRmax2 / (2 * fDz + fOuterConeApex);
0142     } else { // Should we add a check if(fRmax1 > fRmax2)
0143       fOuterConeApex     = 2 * fDz * fRmax2 / NonZero(fRmax1 - fRmax2);
0144       fTanOuterApexAngle = fRmax1 / (2 * fDz + fOuterConeApex);
0145     }
0146
0147     if (fRmax1 == 0. || fRmax2 == 0.) fOuterConeApex = 0.;
0148
0149     if (fRmax1 == 0.) fTanOuterApexAngle = fRmax2 / (2 * fDz);
0150     if (fRmax2 == 0.) fTanOuterApexAngle = fRmax1 / (2 * fDz);
0151
0152     fZNormInner = fTanRMin / NonZero(fSecRMin);
0153     fZNormOuter = -fTanRMax / NonZero(fSecRMax);
0154
0155     fTolIz = fDz - kHalfTolerance;
0156     fTolOz = fDz + kHalfTolerance;
0157
0158     // DetectConvexity();
0159   }
0160
0161   VECCORE_ATT_HOST_DEVICE
0162   void Print() const
0163   {
0164     printf("ConeStruct :  {rmin1 %.2f, rmax1 %.2f, rmin2 %.2f, "
0165            "rmax2 %.2f, dz %.2f, phistart %.2f, deltaphi %.2f}",
0166            fRmin1, fRmax1, fRmin2, fRmax2, fDz, fSPhi, fDPhi);
0167   }
0168
0169   void Print(std::ostream &os) const { os << "UnplacedCone; please implement Print to outstream\n"; }
0170
0171   VECCORE_ATT_HOST_DEVICE
0172   bool IsFullPhi() const { return fDPhi == kTwoPi; }
0173
0174   VECCORE_ATT_HOST_DEVICE
0175   bool Normal(Vector3D<Precision> const &p, Vector3D<Precision> &norm) const
0176   {
0177     int noSurfaces = 0;
0178     Precision rho, pPhi;
0179     Precision distZ, distRMin, distRMax;
0180     Precision distSPhi = kInfLength, distEPhi = kInfLength;
0181     Precision pRMin, widRMin;
0182     Precision pRMax, widRMax;
0183
0184     // const double kHalfTolerance = 0.5 * kTolerance;
0185
0186     Vector3D<Precision> sumnorm(0., 0., 0.), nZ = Vector3D<Precision>(0., 0., 1.);
0187     Vector3D<Precision> nR, nr(0., 0., 0.), nPs, nPe;
0188     norm = sumnorm;
0189
0190     // do not use an extra fabs here -- negative/positive distZ tells us when point is outside or inside
0191     distZ = vecCore::math::Abs(p.z()) - fDz;
0192     rho   = vecCore::math::Sqrt(p.x() * p.x() + p.y() * p.y());
0193
0194     pRMin   = rho - p.z() * fTanRMin;
0195     widRMin = fRmin2 - fDz * fTanRMin;
0196     if (vecCore::math::Abs(_frmin1 - _frmin2) < fInnerTolerance)
0197       distRMin = (rho - _frmin2);
0198     else
0199       distRMin = (pRMin - widRMin) / fSecRMin;
0200
0201     pRMax   = rho - p.z() * fTanRMax;
0202     widRMax = fRmax2 - fDz * fTanRMax;
0203     if (vecCore::math::Abs(_frmax1 - _frmax2) < fOuterTolerance)
0204       distRMax = (rho - _frmax2);
0205     else
0206       distRMax = (pRMax - widRMax) / fSecRMax;
0207
0208     bool inside = distZ < kTolerance && distRMax < fOuterTolerance;
0209     if (fRmin1 || fRmin2) inside &= distRMin > -fInnerTolerance;
0210
0211     distZ    = std::fabs(distZ);
0212     distRMax = std::fabs(distRMax);
0213     distRMin = std::fabs(distRMin);
0214
0215     // keep track of nearest normal, needed in case point is not on a surface
0216     Precision distNearest           = distZ;
0217     Vector3D<Precision> normNearest = nZ;
0218     if (p.z() < 0.) normNearest.Set(0, 0, -1.);
0219
0220     if (!IsFullPhi()) {
0221       if (rho) { // Protected against (0,0,z)
0222         pPhi = vecCore::math::ATan2(p.y(), p.x());
0223
0224         if (pPhi < fSPhi - kHalfTolerance)
0225           pPhi += 2 * kPi;
0226         else if (pPhi > fSPhi + fDPhi + kHalfTolerance)
0227           pPhi -= 2 * kPi;
0228
0229         distSPhi = rho * (pPhi - fSPhi);
0230         distEPhi = rho * (pPhi - fSPhi - fDPhi);
0231         inside   = inside && (distSPhi > -kTolerance) && (distEPhi < kTolerance);
0232         distSPhi = vecCore::math::Abs(distSPhi);
0233         distEPhi = vecCore::math::Abs(distEPhi);
0234       }
0235
0236       else if (!(fRmin1) || !(fRmin2)) {
0237         distSPhi = 0.;
0238         distEPhi = 0.;
0239       }
0240       nPs = Vector3D<Precision>(vecCore::math::Sin(fSPhi), -vecCore::math::Cos(fSPhi), 0);
0241       nPe = Vector3D<Precision>(-vecCore::math::Sin(fSPhi + fDPhi), vecCore::math::Cos(fSPhi + fDPhi), 0);
0242     }
0243
0244     if (rho > kHalfTolerance) {
0245       nR = Vector3D<Precision>(p.x() / rho / fSecRMax, p.y() / rho / fSecRMax, -fTanRMax / fSecRMax);
0246       if (fRmin1 || fRmin2) {
0247         nr = Vector3D<Precision>(-p.x() / rho / fSecRMin, -p.y() / rho / fSecRMin, fTanRMin / fSecRMin);
0248       }
0249     }
0250
0251     if (inside && distZ <= kHalfTolerance) {
0252       noSurfaces++;
0253       if (p.z() >= 0.)
0254         sumnorm += nZ;
0255       else
0256         sumnorm.Set(0, 0, -1.);
0257     }
0258
0259     if (inside && distRMax <= fOuterTolerance) {
0260       noSurfaces++;
0261       sumnorm += nR;
0262     } else if (noSurfaces == 0 && distRMax < distNearest) {
0263       distNearest = distRMax;
0264       normNearest = nR;
0265     }
0266
0267     if (fRmin1 || fRmin2) {
0268       if (inside && distRMin <= fInnerTolerance) {
0269         noSurfaces++;
0270         sumnorm += nr;
0271       } else if (noSurfaces == 0 && distRMin < distNearest) {
0272         distNearest = distRMin;
0273         normNearest = nr;
0274       }
0275     }
0276
0277     if (!IsFullPhi()) {
0278       if (inside && distSPhi <= kHalfTolerance) {
0279         noSurfaces++;
0280         sumnorm += nPs;
0281       } else if (noSurfaces == 0 && distSPhi < distNearest) {
0282         distNearest = distSPhi;
0283         normNearest = nPs;
0284       }
0285       if (inside && distEPhi <= kHalfTolerance) {
0286         noSurfaces++;
0287         sumnorm += nPe;
0288       } else if (noSurfaces == 0 && distEPhi < distNearest) {
0289         // No more check on distNearest, no need to assign to it.
0290         // distNearest = distEPhi;
0291         normNearest = nPe;
0292       }
0293     }
0294     // Final checks
0295     if (noSurfaces == 0)
0296       norm = normNearest;
0297     else if (noSurfaces == 1)
0298       norm = sumnorm;
0299     else
0300       norm = sumnorm.Unit();
0301
0302     bool valid = noSurfaces != 0;
0303     if (noSurfaces > 2) {
0304       // return valid=false for noSurfaces > 2
0305       valid = false;
0306     }
0307
0308     return valid;
0309   }
0310
0311   VECCORE_ATT_HOST_DEVICE
0312   void SetAndCheckSPhiAngle(Precision sPhi)
0313   {
0314     // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
0315     if (sPhi < 0) {
0316       fSPhi = kTwoPi - std::fmod(std::fabs(sPhi), kTwoPi);
0317     } else {
0318       fSPhi = std::fmod(sPhi, kTwoPi);
0319     }
0320     if (fSPhi + fDPhi > kTwoPi) {
0321       fSPhi -= kTwoPi;
0322     }
0323
0324     // Update Wedge
0325     fPhiWedge.SetStartPhi(fSPhi);
0326     // Update cached values.
0327     GetAlongVectorToPhiSector(fSPhi, fAlongPhi1x, fAlongPhi1y);
0328     GetAlongVectorToPhiSector(fSPhi + fDPhi, fAlongPhi2x, fAlongPhi2y);
0329   }
0330
0331   VECCORE_ATT_HOST_DEVICE
0332   void SetAndCheckDPhiAngle(Precision dPhi)
0333   {
0334     if (dPhi >= kTwoPi - 0.5 * kAngTolerance) {
0335       fDPhi = kTwoPi;
0336       fSPhi = 0;
0337     } else {
0338       if (dPhi > 0) {
0339         fDPhi = dPhi;
0340       } else {
0341         //        std::ostringstream message;
0342         //        message << "Invalid dphi.\n"
0343         //                << "Negative or zero delta-Phi (" << dPhi << ")\n";
0344         //        std::cerr<<"UnplacedTube::CheckDPhiAngle(): Fatal error: "<< message.str().c_str() <<"\n";
0345       }
0346     }
0347     // Update Wedge
0348     fPhiWedge.SetDeltaPhi(fDPhi);
0349     // Update cached values.
0350     GetAlongVectorToPhiSector(fSPhi, fAlongPhi1x, fAlongPhi1y);
0351     GetAlongVectorToPhiSector(fSPhi + fDPhi, fAlongPhi2x, fAlongPhi2y);
0352   }
0353
0354   VECCORE_ATT_HOST_DEVICE
0355   static void GetAlongVectorToPhiSector(Precision phi, Precision &x, Precision &y)
0356   {
0357     x = std::cos(phi);
0358     y = std::sin(phi);
0359   }
0360
0361   void SetRmin1(Precision const &arg)
0362   {
0363     fRmin1 = arg;
0364     CalculateCached();
0365   }
0366   void SetRmax1(Precision const &arg)
0367   {
0368     fRmax1 = arg;
0369     CalculateCached();
0370   }
0371   void SetRmin2(Precision const &arg)
0372   {
0373     fRmin2 = arg;
0374     CalculateCached();
0375   }
0376   void SetRmax2(Precision const &arg)
0377   {
0378     fRmax2 = arg;
0379     CalculateCached();
0380   }
0381   void SetDz(Precision const &arg)
0382   {
0383     fDz = arg;
0384     CalculateCached();
0385   }
0386   void SetSPhi(Precision const &arg)
0387   {
0388     fSPhi = arg;
0389     SetAndCheckSPhiAngle(fSPhi);
0390     // DetectConvexity();
0391   }
0392   void SetDPhi(Precision const &arg)
0393   {
0394     fDPhi = arg;
0395     SetAndCheckDPhiAngle(fDPhi);
0396     // DetectConvexity();
0397   }
0398
0399   VECCORE_ATT_HOST_DEVICE
0400   Precision GetTolIz() const { return fTolIz; }
0401   VECCORE_ATT_HOST_DEVICE
0402   Precision GetTolOz() const { return fTolOz; }
0403
0404   VECCORE_ATT_HOST_DEVICE
0405   evolution::Wedge const &GetWedge() const { return fPhiWedge; }
0406
0407   // constructors
0408   ConeStruct() = default;
0409
0410   VECCORE_ATT_HOST_DEVICE
0411   ConeStruct(T const &rmin1, T const &rmax1, T const &rmin2, T const &rmax2, T const &z, T const &sphi, T const &dphi)
0412   {
0413     Init(rmin1, rmax1, rmin2, rmax2, z, sphi, dphi);
0414   }
0415
0416   VECCORE_ATT_HOST_DEVICE
0417   void Init(T const &rmin1, T const &rmax1, T const &rmin2, T const &rmax2, T const &z, T const &sphi, T const &dphi)
0418   {
0419     fRmin1  = rmin1 < 0.0 ? 0.0 : rmin1;
0420     fRmax1  = rmax1;
0421     fRmin2  = rmin2 < 0.0 ? 0.0 : rmin2;
0422     fRmax2  = rmax2;
0423     fDz     = z;
0424     fSPhi   = sphi;
0425     fDPhi   = dphi;
0426     _frmin1 = rmin1;
0427     _frmin2 = rmin2;
0428     _frmax1 = rmax1;
0429     _frmax2 = rmax2;
0430     fPhiWedge.Init(dphi, sphi);
0431     SetAndCheckDPhiAngle(dphi);
0432     SetAndCheckSPhiAngle(sphi);
0433     CalculateCached();
0434     // DetectConvexity();
0435   }
0436 };
0437 } // namespace VECGEOM_IMPL_NAMESPACE
0438 } // namespace vecgeom
0439
0440 #endif