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