VecGeom/volumes/SphereStruct.h

0001 /*
0002  * SphereStruct.h
0003  *
0004  *  Created on: Jul 10, 2017
0005  *      Author: rsehgal
0006  */
0007
0008 #ifndef VECGEOM_SPHERESTRUCT_H_
0009 #define VECGEOM_SPHERESTRUCT_H_
0010
0011 #include "VecGeom/volumes/Wedge_Evolution.h"
0012 //#include "VecGeom/volumes/Wedge.h"
0013 //#include "VecGeom/volumes/ThetaCone.h"
0014 #include "VecGeom/volumes/ThetaCone_Evolution.h"
0015 #include "VecGeom/base/Global.h"
0016
0017 namespace vecgeom {
0018
0019 inline namespace VECGEOM_IMPL_NAMESPACE {
0020
0021 template <typename T = double>
0022 struct SphereStruct {
0023   T fRmin;
0024   T fRmax;
0025   T fSPhi;
0026   T fDPhi;
0027   T fSTheta;
0028   T fDTheta;
0029
0030   // Radial and angular tolerances
0031   Precision fRminTolerance, mkTolerance, //, kAngTolerance, kRadTolerance,
0032       fEpsilon;
0033
0034   // Cached trigonometric values for Phi angle
0035   Precision sinCPhi, cosCPhi, cosHDPhiOT, cosHDPhiIT, sinSPhi, cosSPhi, sinEPhi, cosEPhi, hDPhi, cPhi, ePhi;
0036
0037   // Cached trigonometric values for Theta angle
0038   Precision sinSTheta, cosSTheta, sinETheta, cosETheta, tanSTheta, tanSTheta2, tanETheta, tanETheta2, eTheta;
0039
0040   Precision fabsTanSTheta, fabsTanETheta;
0041
0042   // Flags for identification of section, shell or full sphere
0043   bool fFullPhiSphere, fFullThetaSphere, fFullSphere;
0044
0045   Precision fCubicVolume, fSurfaceArea;
0046
0047   // data members for Theta and Phi
0048   evolution::Wedge fPhiWedge;
0049   ThetaCone fThetaCone;
0050
0051   Precision kAngTolerance;
0052
0053   VECCORE_ATT_HOST_DEVICE
0054   VECGEOM_FORCE_INLINE
0055   void InitializePhiTrigonometry()
0056   {
0057     hDPhi = 0.5 * fDPhi; // half delta phi
0058     cPhi  = fSPhi + hDPhi;
0059     ePhi  = fSPhi + fDPhi;
0060
0061     sinCPhi    = std::sin(cPhi);
0062     cosCPhi    = vecCore::math::Cos(cPhi);
0063     cosHDPhiIT = vecCore::math::Cos(hDPhi - 0.5 * kAngTolerance); // inner/outer tol half dphi
0064     cosHDPhiOT = vecCore::math::Cos(hDPhi + 0.5 * kAngTolerance);
0065     sinSPhi    = std::sin(fSPhi);
0066     cosSPhi    = vecCore::math::Cos(fSPhi);
0067     sinEPhi    = std::sin(ePhi);
0068     cosEPhi    = vecCore::math::Cos(ePhi);
0069   }
0070
0071   VECCORE_ATT_HOST_DEVICE
0072   VECGEOM_FORCE_INLINE
0073   void InitializeThetaTrigonometry()
0074   {
0075     eTheta = fSTheta + fDTheta;
0076
0077     sinSTheta = std::sin(fSTheta);
0078     cosSTheta = vecCore::math::Cos(fSTheta);
0079     sinETheta = std::sin(eTheta);
0080     cosETheta = vecCore::math::Cos(eTheta);
0081
0082     tanSTheta     = sinSTheta / cosSTheta;
0083     fabsTanSTheta = std::fabs(tanSTheta);
0084     tanSTheta2    = tanSTheta * tanSTheta;
0085     tanETheta     = sinETheta / cosETheta;
0086     fabsTanETheta = std::fabs(tanETheta);
0087     tanETheta2    = tanETheta * tanETheta;
0088   }
0089
0090   VECCORE_ATT_HOST_DEVICE
0091   VECGEOM_FORCE_INLINE
0092   void CheckThetaAngles(Precision sTheta, Precision dTheta)
0093   {
0094     if ((sTheta < 0) || (sTheta > kPi)) {
0095       // std::ostringstream message;
0096       // message << "sTheta outside 0-PI range." << std::endl
0097       //       << "Invalid starting Theta angle for solid: " ;//<< GetName();
0098       // return;
0099       // UUtils::Exception("USphere::CheckThetaAngles()", "GeomSolids0002",
0100       //                FatalError, 1, message.str().c_str());
0101     } else {
0102       fSTheta = sTheta;
0103     }
0104     if (dTheta + sTheta >= kPi) {
0105       fDTheta = kPi - sTheta;
0106     } else if (dTheta > 0) {
0107       fDTheta = dTheta;
0108     } else {
0109       /*
0110         std::ostringstream message;
0111       message << "Invalid dTheta." << std::endl
0112               << "Negative delta-Theta (" << dTheta << "), for solid: ";
0113       return;
0114        */
0115       //<< GetName();
0116       // UUtils::Exception("USphere::CheckThetaAngles()", "GeomSolids0002",
0117       //                FatalError, 1, message.str().c_str());
0118     }
0119     if (fDTheta - fSTheta < kPi) {
0120       fFullThetaSphere = false;
0121     } else {
0122       fFullThetaSphere = true;
0123     }
0124     fFullSphere = fFullPhiSphere && fFullThetaSphere;
0125
0126     InitializeThetaTrigonometry();
0127   }
0128
0129   VECCORE_ATT_HOST_DEVICE
0130   VECGEOM_FORCE_INLINE
0131   void CheckSPhiAngle(Precision sPhi)
0132   {
0133     // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
0134
0135     if (sPhi < 0) {
0136       fSPhi = 2 * kPi - std::fmod(std::fabs(sPhi), 2 * kPi);
0137     } else {
0138       fSPhi = std::fmod(sPhi, 2 * kPi);
0139     }
0140     if (fSPhi + fDPhi > 2 * kPi) {
0141       fSPhi -= 2 * kPi;
0142     }
0143   }
0144
0145   VECCORE_ATT_HOST_DEVICE
0146   VECGEOM_FORCE_INLINE
0147   void CheckDPhiAngle(Precision dPhi)
0148   {
0149     if (dPhi >= 2 * kPi - kAngTolerance * 0.5) {
0150
0151       fDPhi          = 2 * kPi;
0152       fSPhi          = 0;
0153       fFullPhiSphere = true;
0154
0155     } else {
0156       fFullPhiSphere = false;
0157       if (dPhi > 0) {
0158         fDPhi = dPhi;
0159       } else {
0160         /*
0161       std::ostringstream message;
0162       message << "Invalid dphi." << std::endl
0163               << "Negative delta-Phi (" << dPhi << "), for solid: ";
0164       return;
0165          */
0166
0167         // << GetName();
0168         // UUtils::Exception("USphere::CheckDPhiAngle()", "GeomSolids0002",
0169         //                FatalError, 1, message.str().c_str());
0170       }
0171     }
0172   }
0173
0174   VECCORE_ATT_HOST_DEVICE
0175   VECGEOM_FORCE_INLINE
0176   void CheckPhiAngles(Precision sPhi, Precision dPhi)
0177   {
0178     CheckDPhiAngle(dPhi);
0179     // if (!fFullPhiSphere && sPhi) { CheckSPhiAngle(sPhi); }
0180     if (!fFullPhiSphere) {
0181       CheckSPhiAngle(sPhi);
0182     }
0183     fFullSphere = fFullPhiSphere && fFullThetaSphere;
0184
0185     InitializePhiTrigonometry();
0186   }
0187
0188   // VECCORE_ATT_HOST_DEVICE
0189   // VECGEOM_FORCE_INLINE
0190   void SetInsideRadius(Precision newRmin)
0191   {
0192     fRmin          = newRmin;
0193     fRminTolerance = (fRmin) ? std::max(kRadTolerance, fEpsilon * fRmin) : 0;
0194     Initialize();
0195 #ifndef VECCORE_CUDA
0196     CalcCapacity();
0197     CalcSurfaceArea();
0198 #endif
0199   }
0200
0201   // VECCORE_ATT_HOST_DEVICE
0202   // VECGEOM_FORCE_INLINE
0203   void SetInnerRadius(Precision newRmin) { SetInsideRadius(newRmin); }
0204
0205   // VECCORE_ATT_HOST_DEVICE
0206   // VECGEOM_FORCE_INLINE
0207   void SetOuterRadius(Precision newRmax)
0208   {
0209     fRmax       = newRmax;
0210     mkTolerance = std::max(kRadTolerance,
0211                            fEpsilon * fRmax); // RELOOK at kTolerance, may be we will take directly from base/global.h
0212     Initialize();
0213 #ifndef VECCORE_CUDA
0214     CalcCapacity();
0215     CalcSurfaceArea();
0216 #endif
0217   }
0218
0219   // VECCORE_ATT_HOST_DEVICE
0220   // VECGEOM_FORCE_INLINE
0221   void SetStartPhiAngle(Precision newSPhi, bool compute = true)
0222   {
0223     // Flag 'compute' can be used to explicitely avoid recomputation of
0224     // trigonometry in case SetDeltaPhiAngle() is invoked afterwards
0225
0226     CheckSPhiAngle(newSPhi);
0227     fFullPhiSphere = false;
0228     if (compute) {
0229       InitializePhiTrigonometry();
0230     }
0231     Initialize();
0232 #ifndef VECCORE_CUDA
0233     CalcCapacity();
0234     CalcSurfaceArea();
0235 #endif
0236   }
0237
0238   // VECCORE_ATT_HOST_DEVICE
0239   // VECGEOM_FORCE_INLINE
0240   void SetDeltaPhiAngle(Precision newDPhi)
0241   {
0242     CheckPhiAngles(fSPhi, newDPhi);
0243     Initialize();
0244 #ifndef VECCORE_CUDA
0245     CalcCapacity();
0246     CalcSurfaceArea();
0247 #endif
0248   }
0249
0250   // VECCORE_ATT_HOST_DEVICE
0251   // VECGEOM_FORCE_INLINE
0252   void SetStartThetaAngle(Precision newSTheta)
0253   {
0254     CheckThetaAngles(newSTheta, fDTheta);
0255     Initialize();
0256 #ifndef VECCORE_CUDA
0257     CalcCapacity();
0258     CalcSurfaceArea();
0259 #endif
0260   }
0261
0262   // VECCORE_ATT_HOST_DEVICE
0263   // VECGEOM_FORCE_INLINE
0264   void SetDeltaThetaAngle(Precision newDTheta)
0265   {
0266     CheckThetaAngles(fSTheta, newDTheta);
0267     Initialize();
0268 #ifndef VECCORE_CUDA
0269     CalcCapacity();
0270     CalcSurfaceArea();
0271 #endif
0272   }
0273
0274   VECCORE_ATT_HOST_DEVICE
0275   VECGEOM_FORCE_INLINE
0276   void Initialize()
0277   {
0278     fCubicVolume = 0.;
0279     fSurfaceArea = 0.;
0280   }
0281
0282   // Constructor
0283   VECCORE_ATT_HOST_DEVICE
0284   VECGEOM_FORCE_INLINE
0285   SphereStruct(T pRmin, T pRmax, T pSPhi, T pDPhi, T pSTheta, T pDTheta)
0286       : fRmin(pRmin), fRmax(pRmax), fSPhi(pSPhi), fDPhi(pDPhi), fSTheta(pSTheta), fDTheta(pDTheta), fRminTolerance(0),
0287         mkTolerance(0), fEpsilon(kEpsilon), sinCPhi(0), cosCPhi(0), cosHDPhiOT(0), cosHDPhiIT(0), sinSPhi(0),
0288         cosSPhi(0), sinEPhi(0), cosEPhi(0), hDPhi(0), cPhi(0), ePhi(0), sinSTheta(0), cosSTheta(0), sinETheta(0),
0289         cosETheta(0), tanSTheta(0), tanSTheta2(0), tanETheta(0), tanETheta2(0), eTheta(0), fFullPhiSphere(true),
0290         fFullThetaSphere(true), fFullSphere(true), fCubicVolume(0.), fSurfaceArea(0.), fPhiWedge(pDPhi, pSPhi),
0291         fThetaCone(pSTheta, pDTheta)
0292   {
0293
0294     kAngTolerance  = 1e-9;
0295     fRminTolerance = (fRmin) ? Max(kRadTolerance, fEpsilon * fRmin) : 0;
0296     mkTolerance    = Max(kRadTolerance, fEpsilon * fRmax);
0297
0298     CheckPhiAngles(pSPhi, pDPhi);
0299     CheckThetaAngles(pSTheta, pDTheta);
0300 #ifndef VECCORE_CUDA
0301     CalcCapacity();
0302     CalcSurfaceArea();
0303 #endif
0304   }
0305
0306   VECCORE_ATT_HOST_DEVICE
0307   VECGEOM_FORCE_INLINE
0308   SphereStruct() {}
0309
0310   //#ifndef VECCORE_CUDA
0311   VECCORE_ATT_HOST_DEVICE
0312   VECGEOM_FORCE_INLINE
0313   void CalcCapacity()
0314   {
0315     if (fCubicVolume != 0.) {
0316       ;
0317     } else {
0318       fCubicVolume = fDPhi * (vecCore::math::Cos(fSTheta) - vecCore::math::Cos(fSTheta + fDTheta)) *
0319                      (fRmax * fRmax * fRmax - fRmin * fRmin * fRmin) / 3.;
0320     }
0321   }
0322   VECCORE_ATT_HOST_DEVICE
0323   VECGEOM_FORCE_INLINE
0324   void CalcSurfaceArea()
0325   {
0326
0327     if (fSurfaceArea != 0.) {
0328       ;
0329     } else {
0330       Precision Rsq = fRmax * fRmax;
0331       Precision rsq = fRmin * fRmin;
0332
0333       fSurfaceArea = fDPhi * (rsq + Rsq) * (cosSTheta - cosETheta);
0334       if (!fFullPhiSphere) {
0335         fSurfaceArea = fSurfaceArea + fDTheta * (Rsq - rsq);
0336       }
0337       if (fSTheta > 0) {
0338         Precision acos1 = 0.;
0339         if (fDPhi != kTwoPi)
0340           acos1 = vecCore::math::ACos((sinSTheta * sinSTheta) * vecCore::math::Cos(fDPhi) + (cosSTheta * cosSTheta));
0341         if (fDPhi > kPi) {
0342           fSurfaceArea = fSurfaceArea + 0.5 * (Rsq - rsq) * (2 * kPi - acos1);
0343         } else {
0344           fSurfaceArea = fSurfaceArea + 0.5 * (Rsq - rsq) * acos1;
0345         }
0346       }
0347       if (eTheta < kPi) {
0348         Precision acos2 = 0.;
0349         if (fDPhi != kTwoPi)
0350           acos2 = vecCore::math::ACos((sinETheta * sinETheta) * vecCore::math::Cos(fDPhi) + (cosETheta * cosETheta));
0351         if (fDPhi > kPi) {
0352           fSurfaceArea = fSurfaceArea + 0.5 * (Rsq - rsq) * (2 * kPi - acos2);
0353         } else {
0354           fSurfaceArea = fSurfaceArea + 0.5 * (Rsq - rsq) * acos2;
0355         }
0356       }
0357     }
0358   }
0359   //#endif
0360 };
0361 } // namespace VECGEOM_IMPL_NAMESPACE
0362 } // namespace vecgeom
0363
0364 #endif /* VOLUMES_SPHERESTRUCT_H_ */