VecGeom/volumes/CoaxialConesStruct.h

0001 /// @file CoaxialConesStruct.h
0002 /// @author Raman Sehgal (raman.sehgal@cern.ch)
0003
0004 #ifndef VECGEOM_VOLUMES_COAXIALCONESSTRUCT_H_
0005 #define VECGEOM_VOLUMES_COAXIALCONESSTRUCT_H_
0006 #include "VecGeom/base/Global.h"
0007 #include "VecGeom/base/Vector3D.h"
0008 #include "VecGeom/base/Vector.h"
0009 #include "VecGeom/volumes/ConeStruct.h"
0010
0011 namespace vecgeom {
0012
0013 inline namespace VECGEOM_IMPL_NAMESPACE {
0014
0015 template <typename T = double>
0016 struct CoaxialConesStruct {
0017
0018   // CoaxialCone parameters
0019   unsigned int fNumOfCones;
0020   T fDz;
0021   T fSPhi;
0022   T fDPhi;
0023   Vector<T> fRmin1Vect;
0024   Vector<T> fRmax1Vect;
0025   Vector<T> fRmin2Vect;
0026   Vector<T> fRmax2Vect;
0027
0028   /*
0029   ** A function to Print the sequence of Coaxial cones,
0030   ** Useful for debugging the Coaxial cones
0031   */
0032   void Print() const
0033   {
0034
0035     for (unsigned int i = 0; i < fConeStructVector.size(); i++) {
0036       std::cout << std::endl << "====== Index of cone struct : " << i << " ======" << std::endl;
0037       fConeStructVector[i]->Print();
0038       std::cout << std::endl;
0039     }
0040   }
0041
0042   VECCORE_ATT_HOST_DEVICE
0043   CoaxialConesStruct() {}
0044
0045   VECCORE_ATT_HOST_DEVICE
0046   CoaxialConesStruct(unsigned int numOfCones, T *rmin1Vect, T *rmax1Vect, T *rmin2Vect, T *rmax2Vect, T dz, T sphi,
0047                      T dphi)
0048       : fNumOfCones(numOfCones), fDz(dz), fSPhi(sphi), fDPhi(dphi)
0049   {
0050
0051     for (unsigned int i = 0; i < fNumOfCones; i++) {
0052       fRmin1Vect.push_back(rmin1Vect[i]);
0053       fRmax1Vect.push_back(rmax1Vect[i]);
0054       fRmin2Vect.push_back(rmin2Vect[i]);
0055       fRmax2Vect.push_back(rmax2Vect[i]);
0056     }
0057
0058     for (unsigned int i = 0; i < fNumOfCones; i++) {
0059       fConeStructVector.push_back(
0060           new ConeStruct<T>(fRmin1Vect[i], fRmax1Vect[i], fRmin2Vect[i], fRmax2Vect[i], dz, sphi, dphi));
0061       if (i == 0) {
0062         fMinR = fRmin1Vect[i] < fRmin2Vect[i] ? fRmin1Vect[i] : fRmin2Vect[i];
0063       }
0064
0065       if (i == fNumOfCones - 1) {
0066         fMaxR = fRmax1Vect[i] > fRmax2Vect[i] ? fRmax1Vect[i] : fRmax2Vect[i];
0067       }
0068     }
0069   }
0070
0071   VECCORE_ATT_HOST_DEVICE
0072   CoaxialConesStruct(Vector<Precision> rmin1Vect, Vector<Precision> rmax1Vect, Vector<Precision> rmin2Vect,
0073                      Vector<Precision> rmax2Vect, T dz, T sphi, T dphi)
0074       : fNumOfCones(rmin1Vect.size()), fDz(dz), fSPhi(sphi), fDPhi(dphi), fRmin1Vect(rmin1Vect), fRmax1Vect(rmax1Vect),
0075         fRmin2Vect(rmin2Vect), fRmax2Vect(rmax2Vect)
0076   {
0077
0078     for (unsigned int i = 0; i < fNumOfCones; i++) {
0079       fConeStructVector.push_back(
0080           new ConeStruct<T>(fRmin1Vect[i], fRmax1Vect[i], fRmin2Vect[i], fRmax2Vect[i], dz, sphi, dphi));
0081       if (i == 0) {
0082         fMinR = fRmin1Vect[i] < fRmin2Vect[i] ? fRmin1Vect[i] : fRmin2Vect[i];
0083       }
0084
0085       if (i == fNumOfCones - 1) {
0086         fMaxR = fRmax1Vect[i] > fRmax2Vect[i] ? fRmax1Vect[i] : fRmax2Vect[i];
0087       }
0088     }
0089   }
0090
0091   // Vector of Cones
0092   Vector<ConeStruct<T> *> fConeStructVector;
0093
0094   T fSurfaceArea; // area of the surface
0095   T fCubicVolume; // volume
0096   T fMaxR;
0097   T fMinR;
0098
0099   // Precalculated cached values
0100   VECCORE_ATT_HOST_DEVICE
0101   Precision Capacity()
0102   {
0103     Precision volume = 0.;
0104     for (unsigned int i = 0; i < fConeStructVector.size(); i++) {
0105       volume += fConeStructVector[i]->Capacity();
0106     }
0107     return volume;
0108   }
0109
0110   VECCORE_ATT_HOST_DEVICE
0111   Precision ConicalSurfaceArea()
0112   {
0113     Precision conicalSurfaceArea = 0.;
0114     Precision mmin, mmax, dmin, dmax;
0115
0116     for (unsigned int i = 0; i < fConeStructVector.size(); i++) {
0117
0118       mmin = (fConeStructVector[i]->fRmin1 + fConeStructVector[i]->fRmin2) * 0.5;
0119       mmax = (fConeStructVector[i]->fRmax1 + fConeStructVector[i]->fRmax2) * 0.5;
0120       dmin = (fConeStructVector[i]->fRmin2 - fConeStructVector[i]->fRmin1);
0121       dmax = (fConeStructVector[i]->fRmax2 - fConeStructVector[i]->fRmax1);
0122
0123       conicalSurfaceArea +=
0124           fConeStructVector[i]->fDPhi *
0125           (mmin * vecCore::math::Sqrt(dmin * dmin + 4 * fConeStructVector[i]->fDz * fConeStructVector[i]->fDz) +
0126            mmax * vecCore::math::Sqrt(dmax * dmax + 4 * fConeStructVector[i]->fDz * fConeStructVector[i]->fDz));
0127     }
0128     return conicalSurfaceArea;
0129   }
0130
0131   VECCORE_ATT_HOST_DEVICE
0132   Precision SurfaceAreaLowerZPlanes(int index)
0133   {
0134     return fConeStructVector[index]->fDPhi * 0.5 *
0135            (fConeStructVector[index]->fRmax1 * fConeStructVector[index]->fRmax1 -
0136             fConeStructVector[index]->fRmin1 * fConeStructVector[index]->fRmin1);
0137   }
0138
0139   VECCORE_ATT_HOST_DEVICE
0140   Precision SurfaceAreaUpperZPlanes(int index)
0141   {
0142     return fConeStructVector[index]->fDPhi * 0.5 *
0143            (fConeStructVector[index]->fRmax2 * fConeStructVector[index]->fRmax2 -
0144             fConeStructVector[index]->fRmin2 * fConeStructVector[index]->fRmin2);
0145   }
0146
0147   VECCORE_ATT_HOST_DEVICE
0148   Precision SurfaceAreaOfZPlanes()
0149   {
0150     Precision surfaceAreaOfZPlanes = 0.;
0151     for (unsigned int i = 0; i < fConeStructVector.size(); i++) {
0152       surfaceAreaOfZPlanes = SurfaceAreaLowerZPlanes(i) + SurfaceAreaUpperZPlanes(i);
0153     }
0154     return surfaceAreaOfZPlanes;
0155   }
0156
0157   VECCORE_ATT_HOST_DEVICE
0158   Precision SurfaceArea() { return (ConicalSurfaceArea() + SurfaceAreaOfZPlanes()); }
0159
0160   VECCORE_ATT_HOST_DEVICE
0161   bool Normal(Vector3D<Precision> const &p, Vector3D<Precision> &norm) const
0162   {
0163
0164     norm.Set(0.);
0165     bool valid = false;
0166     for (unsigned int i = 0; i < fConeStructVector.size(); i++) {
0167       bool validNormal = false;
0168       Vector3D<Precision> normal(0., 0., 0.);
0169       validNormal = fConeStructVector[i]->Normal(p, normal);
0170       if (validNormal) {
0171         norm += normal;
0172       }
0173       valid |= validNormal;
0174     }
0175     if (valid) {
0176       norm.Normalize();
0177     } else {
0178       norm.Set(0., 0., 1.);
0179     }
0180     return valid;
0181   }
0182 };
0183
0184 } // namespace VECGEOM_IMPL_NAMESPACE
0185 } // namespace vecgeom
0186
0187 #endif