VecGeom/volumes/PolyconeStruct.h

0001 /*
0002  * ConeStruct.h
0003  *
0004  *  Created on: May 11, 2017
0005  *      Author: Raman Sehgal
0006  */
0007 #ifndef VECGEOM_POLYCONESTRUCT_H_
0008 #define VECGEOM_POLYCONESTRUCT_H_
0009
0010 #include "VecGeom/base/Global.h"
0011 #include "VecGeom/volumes/Wedge_Evolution.h"
0012 #include "VecGeom/volumes/ConeStruct.h"
0013 #include "VecGeom/base/Vector.h"
0014 #include "VecGeom/volumes/PolyconeHistorical.h"
0015 #include "VecGeom/volumes/PolyconeSection.h"
0016 #include "VecGeom/base/Array.h"
0017
0018 namespace vecgeom {
0019
0020 inline namespace VECGEOM_IMPL_NAMESPACE {
0021
0022 // a plain and lightweight struct to encapsulate data members of a Polycone
0023 template <typename T = double>
0024 struct PolyconeStruct {
0025
0026   bool fEqualRmax;
0027   bool fContinuityOverAll;
0028   bool fConvexityPossible;
0029
0030   evolution::Wedge fPhiWedge;
0031   Precision fStartPhi;
0032   Precision fDeltaPhi;
0033   unsigned int fNz;
0034
0035   Vector<PolyconeSection> fSections;
0036   Vector<Precision> fZs;
0037   PolyconeHistorical *fOriginal_parameters;
0038
0039   // Data member to hold the line segment that form Polycone boundary
0040   Vector<Vector3D<Precision>> fRMinTwoDVec;
0041   Vector<Vector3D<Precision>> fRMaxTwoDVec;
0042   Vector<Vector3D<Precision>> fTwoDVec;
0043
0044   VECCORE_ATT_HOST_DEVICE
0045   bool CheckContinuity(const Precision rOuter[], const Precision rInner[], const Precision zPlane[],
0046                        Vector<Precision> &newROuter, Vector<Precision> &newRInner, Vector<Precision> &newZPlane)
0047   {
0048     Vector<Precision> rOut, rIn;
0049     Vector<Precision> zPl;
0050     rOut.push_back(rOuter[0]);
0051     rIn.push_back(rInner[0]);
0052     zPl.push_back(zPlane[0]);
0053     for (unsigned int j = 1; j < fNz; j++) {
0054
0055       if (j == fNz - 1) {
0056         rOut.push_back(rOuter[j]);
0057         rIn.push_back(rInner[j]);
0058         zPl.push_back(zPlane[j]);
0059       } else {
0060         if ((zPlane[j] != zPlane[j + 1]) || (rOuter[j] != rOuter[j + 1])) {
0061           rOut.push_back(rOuter[j]);
0062           rOut.push_back(rOuter[j]);
0063
0064           zPl.push_back(zPlane[j]);
0065           zPl.push_back(zPlane[j]);
0066
0067           rIn.push_back(rInner[j]);
0068           rIn.push_back(rInner[j]);
0069
0070         } else {
0071           rOut.push_back(rOuter[j]);
0072           zPl.push_back(zPlane[j]);
0073           rIn.push_back(rInner[j]);
0074         }
0075       }
0076     }
0077
0078     if (rOut.size() % 2 != 0) {
0079       // fNz is odd, the adding of the last item did not happen in the loop.
0080       rOut.push_back(rOut[rOut.size() - 1]);
0081       rIn.push_back(rIn[rIn.size() - 1]);
0082       zPl.push_back(zPl[zPl.size() - 1]);
0083     }
0084
0085     /* Creating a new temporary Reduced polycone with desired data elements,
0086      *  which makes sure that denominator will never be zero (hence avoiding FPE(division by zero)),
0087      *  while calculating slope.
0088      *
0089      *  This will be the minimum polycone,i.e. no extra section which
0090      *  affect its shape
0091      */
0092
0093     for (size_t j = 0; j < rOut.size();) {
0094
0095       if (zPl[j] != zPl[j + 1]) {
0096
0097         newZPlane.push_back(zPl[j]);
0098         newZPlane.push_back(zPl[j + 1]);
0099         newROuter.push_back(rOut[j]);
0100         newROuter.push_back(rOut[j + 1]);
0101         newRInner.push_back(rIn[j]);
0102         newRInner.push_back(rIn[j + 1]);
0103       }
0104
0105       j = j + 2;
0106     }
0107     // Minimum polycone construction over
0108
0109     // Checking Slope continuity and Rmax Continuity
0110     bool contRmax  = CheckContinuityInRmax(newROuter);
0111     bool contSlope = CheckContinuityInSlope(newROuter, newZPlane);
0112
0113     // If both are true then the polycone can be convex
0114     // but still final convexity depends on Inner Radius also.
0115     return (contRmax && contSlope);
0116   }
0117
0118   VECCORE_ATT_HOST_DEVICE
0119   bool CheckContinuityInRmax(const Vector<Precision> &rOuter)
0120   {
0121     bool continuous  = true;
0122     unsigned int len = rOuter.size();
0123     if (len > 2) {
0124       for (unsigned int j = 1; j < len;) {
0125         if (j != (len - 1)) continuous &= (rOuter[j] == rOuter[j + 1]);
0126         j = j + 2;
0127       }
0128     }
0129     return continuous;
0130   }
0131
0132   VECCORE_ATT_HOST_DEVICE
0133   bool CheckContinuityInSlope(const Vector<Precision> &rOuter, const Vector<Precision> &zPlane)
0134   {
0135
0136     bool continuous      = true;
0137     Precision startSlope = kInfLength;
0138
0139     // Doing the actual slope calculation here, and checking continuity,
0140     for (size_t j = 0; j < rOuter.size(); j = j + 2) {
0141       Precision currentSlope = (rOuter[j + 1] - rOuter[j]) / (zPlane[j + 1] - zPlane[j]);
0142       continuous &= (currentSlope <= startSlope);
0143       startSlope = currentSlope;
0144     }
0145     return continuous;
0146   }
0147
0148   VECCORE_ATT_HOST_DEVICE
0149   void Init(Precision phiStart,       // initial phi starting angle
0150             Precision phiTotal,       // total phi angle
0151             unsigned int numZPlanes,  // number of z planes
0152             const Precision zPlane[], // position of z planes
0153             const Precision rInner[], // tangent distance to inner surface
0154             const Precision rOuter[])
0155   {
0156
0157     SetAndCheckDPhiAngle(phiTotal);
0158     SetAndCheckSPhiAngle(phiStart);
0159     fNz                = numZPlanes;
0160     Precision *zPlaneR = new Precision[numZPlanes];
0161     Precision *rInnerR = new Precision[numZPlanes];
0162     Precision *rOuterR = new Precision[numZPlanes];
0163     for (unsigned int i = 0; i < numZPlanes; i++) {
0164       zPlaneR[i] = zPlane[i];
0165       rInnerR[i] = rInner[i];
0166       rOuterR[i] = rOuter[i];
0167     }
0168     if (zPlane[0] > zPlane[numZPlanes - 1]) {
0169       // Reverse the arrays
0170       for (unsigned int i = 0; i < numZPlanes; i++) {
0171         zPlaneR[i] = zPlane[numZPlanes - 1 - i];
0172         rInnerR[i] = rInner[numZPlanes - 1 - i];
0173         rOuterR[i] = rOuter[numZPlanes - 1 - i];
0174       }
0175     }
0176
0177     // Conversion for angles
0178     if (phiTotal <= 0. || phiTotal > kTwoPi - kTolerance) {
0179       // phiIsOpen=false;
0180       fStartPhi = 0;
0181       fDeltaPhi = kTwoPi;
0182     } else {
0183       //
0184       // Convert phi into our convention
0185       //
0186       fStartPhi = phiStart;
0187       while (fStartPhi < 0)
0188         fStartPhi += kTwoPi;
0189     }
0190
0191     // Calculate RMax of Polycone in order to determine convexity of sections
0192     //
0193     Precision RMaxextent = rOuterR[0];
0194
0195     Vector<Precision> newROuter, newZPlane, newRInner;
0196     fContinuityOverAll &= CheckContinuity(rOuterR, rInnerR, zPlaneR, newROuter, newRInner, newZPlane);
0197     fConvexityPossible &= (newRInner[0] == 0.);
0198
0199     Precision startRmax = newROuter[0];
0200     for (unsigned int j = 1; j < newROuter.size(); j++) {
0201       fEqualRmax &= (startRmax == newROuter[j]);
0202       startRmax = newROuter[j];
0203       fConvexityPossible &= (newRInner[j] == 0.);
0204     }
0205
0206     for (unsigned int j = 1; j < numZPlanes; j++) {
0207
0208       if (rOuterR[j] > RMaxextent) RMaxextent = rOuterR[j];
0209
0210       if (rInnerR[j] > rOuterR[j]) {
0211 #ifndef VECCORE_CUDA
0212         std::cerr << "Cannot create Polycone with rInner > rOuter for the same Z"
0213                   << "\n"
0214                   << "        rInner > rOuter for the same Z !\n"
0215                   << "        rMin[" << j << "] = " << rInner[j] << " -- rMax[" << j << "] = " << rOuter[j];
0216 #endif
0217       }
0218     }
0219
0220     Precision prevZ = zPlaneR[0], prevRmax = 0, prevRmin = 0;
0221     int dirZ = 1;
0222     if (zPlaneR[1] < zPlaneR[0]) dirZ = -1;
0223
0224     for (unsigned int i = 0; i < numZPlanes; ++i) {
0225       if ((i < numZPlanes - 1) && (zPlaneR[i] == zPlaneR[i + 1])) {
0226         if ((rInnerR[i] > rOuterR[i + 1]) || (rInnerR[i + 1] > rOuterR[i])) {
0227 #ifndef VECCORE_CUDA
0228           std::cerr << "Cannot create a Polycone with no contiguous segments." << std::endl
0229                     << "                Segments are not contiguous !" << std::endl
0230                     << "                rMin[" << i << "] = " << rInnerR[i] << " -- rMax[" << i + 1
0231                     << "] = " << rOuterR[i + 1] << std::endl
0232                     << "                rMin[" << i + 1 << "] = " << rInnerR[i + 1] << " -- rMax[" << i
0233                     << "] = " << rOuterR[i];
0234 #endif
0235         }
0236       }
0237
0238       Precision rMin = rInnerR[i];
0239
0240       Precision rMax = rOuterR[i];
0241       Precision z    = zPlaneR[i];
0242
0243       // i has to be at least one to complete a section
0244       if (i > 0) {
0245         // GL: I had to add kTolerance here, otherwise this polycone would get an
0246         //     extra 0-length ZSection on the GPU -- see VECGEOM-578
0247         if (((z > prevZ + kTolerance) && (dirZ > 0)) || ((z < prevZ - kTolerance) && (dirZ < 0))) {
0248           if (dirZ * (z - prevZ) < 0) {
0249 #ifndef VECCORE_CUDA
0250             std::cerr << "Cannot create a Polycone with different Z directions.Use GenericPolycone." << std::endl
0251                       << "              ZPlane is changing direction  !" << std::endl
0252                       << "  zPlane[0] = " << zPlaneR[0] << " -- zPlane[1] = " << zPlaneR[1] << std::endl
0253                       << "  zPlane[" << i - 1 << "] = " << zPlaneR[i - 1] << " -- rPlane[" << i << "] = " << zPlaneR[i];
0254 #endif
0255           }
0256
0257           ConeStruct<Precision> *solid;
0258
0259           Precision dz = (z - prevZ) / 2;
0260
0261           solid = new ConeStruct<Precision>(prevRmin, prevRmax, rMin, rMax, dz, phiStart, phiTotal);
0262
0263           fZs.push_back(z);
0264           int zi          = fZs.size() - 1;
0265           Precision shift = fZs[zi - 1] + 0.5 * (fZs[zi] - fZs[zi - 1]);
0266
0267           PolyconeSection section;
0268           section.fShift = shift;
0269           section.fSolid = solid;
0270
0271           section.fConvex = !((rMax < prevRmax) || (rMax < RMaxextent) || (prevRmax < RMaxextent));
0272
0273           fSections.push_back(section);
0274           fRMinTwoDVec.push_back(Vector3D<Precision>(prevRmin, prevZ, 0));
0275           fRMinTwoDVec.push_back(Vector3D<Precision>(rMin, z, 0));
0276           fRMaxTwoDVec.push_back(Vector3D<Precision>(prevRmax, prevZ, 0));
0277           fRMaxTwoDVec.push_back(Vector3D<Precision>(rMax, z, 0));
0278         }
0279       } else { // for i == 0 just push back first z plane
0280         fZs.push_back(z);
0281       }
0282
0283       prevZ    = z;
0284       prevRmin = rMin;
0285       prevRmax = rMax;
0286     }
0287     for (auto val : fRMaxTwoDVec) {
0288       fTwoDVec.push_back(val);
0289     }
0290
0291     for (int k = fRMinTwoDVec.size() - 1; k >= 0; k--) {
0292       fTwoDVec.push_back(fRMinTwoDVec[k]);
0293     }
0294     fOriginal_parameters                  = new PolyconeHistorical(numZPlanes);
0295     fOriginal_parameters->fHStart_angle   = phiStart;
0296     fOriginal_parameters->fHOpening_angle = phiTotal;
0297     for (unsigned int i = 0; i < numZPlanes; i++) {
0298       fOriginal_parameters->fHZ_values[i] = zPlaneR[i];
0299       fOriginal_parameters->fHRmin[i]     = rInnerR[i];
0300       fOriginal_parameters->fHRmax[i]     = rOuterR[i];
0301     }
0302
0303     delete[] zPlaneR;
0304     delete[] rInnerR;
0305     delete[] rOuterR;
0306   }
0307
0308   VECCORE_ATT_HOST_DEVICE
0309   PolyconeHistorical *GetOriginalParameters() const { return fOriginal_parameters; }
0310
0311   VECCORE_ATT_HOST_DEVICE unsigned int GetNz() const { return fNz; }
0312
0313   VECCORE_ATT_HOST_DEVICE
0314   int GetNSections() const { return fSections.size(); }
0315
0316   VECCORE_ATT_HOST_DEVICE
0317   int GetSectionIndex(Precision zposition) const
0318   {
0319     // TODO: consider binary search
0320     // TODO: consider making these comparisons tolerant in case we need it
0321     if (zposition < fZs[0]) return -1;
0322     for (unsigned int i = 0; i < fSections.size(); ++i) {
0323       if (zposition >= fZs[i] && zposition <= fZs[i + 1]) return i;
0324     }
0325     return -2;
0326   }
0327
0328   VECCORE_ATT_HOST_DEVICE
0329   PolyconeSection const &GetSection(Precision zposition) const
0330   {
0331     // TODO: consider binary search
0332     int i = GetSectionIndex(zposition);
0333     return fSections[i];
0334   }
0335
0336   VECCORE_ATT_HOST_DEVICE
0337   // GetSection if index is known
0338   PolyconeSection const &GetSection(int index) const { return fSections[index]; }
0339
0340   VECCORE_ATT_HOST_DEVICE
0341   Precision GetRminAtPlane(int index) const
0342   {
0343     int nsect = fSections.size();
0344     assert(index >= 0 && index <= nsect);
0345     if (index == nsect)
0346       return fSections[index - 1].fSolid->fRmin2; // GetRmin2();
0347     else
0348       return fSections[index].fSolid->fRmin1; // GetRmin1();
0349   }
0350
0351   VECCORE_ATT_HOST_DEVICE
0352   Precision GetRmaxAtPlane(int index) const
0353   {
0354     int nsect = fSections.size();
0355     assert(index >= 0 || index <= nsect);
0356     if (index == nsect)
0357       return fSections[index - 1].fSolid->fRmax2; // GetRmax2();
0358     else
0359       return fSections[index].fSolid->fRmax1; // GetRmax1();
0360   }
0361
0362   VECCORE_ATT_HOST_DEVICE
0363   Precision GetZAtPlane(unsigned int index) const
0364   {
0365     assert(index <= fSections.size());
0366     return fZs[index];
0367   }
0368
0369   VECCORE_ATT_HOST_DEVICE
0370   void SetAndCheckSPhiAngle(Precision sPhi)
0371   {
0372     // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
0373     if (sPhi < 0) {
0374       fStartPhi = kTwoPi - std::fmod(std::fabs(sPhi), kTwoPi);
0375     } else {
0376       fStartPhi = std::fmod(sPhi, kTwoPi);
0377     }
0378     if (fStartPhi + fDeltaPhi > kTwoPi) {
0379       fStartPhi -= kTwoPi;
0380     }
0381
0382     // Update Wedge
0383     fPhiWedge.SetStartPhi(fStartPhi);
0384     fPhiWedge.UpdateNormals();
0385   }
0386
0387   VECCORE_ATT_HOST_DEVICE
0388   void SetAndCheckDPhiAngle(Precision dPhi)
0389   {
0390     if (dPhi >= kTwoPi - 0.5 * kAngTolerance) {
0391       fDeltaPhi = kTwoPi;
0392       fStartPhi = 0;
0393     } else {
0394       if (dPhi > 0) {
0395         fDeltaPhi = dPhi;
0396       } else {
0397         //        std::ostringstream message;
0398         //        message << "Invalid dphi.\n"
0399         //                << "Negative or zero delta-Phi (" << dPhi << ")\n";
0400         //        std::cout<<"UnplacedTube::CheckDPhiAngle(): Fatal error: "<< message.str().c_str() <<"\n";
0401       }
0402     }
0403     // Update Wedge
0404     fPhiWedge.SetDeltaPhi(fDeltaPhi);
0405     fPhiWedge.UpdateNormals();
0406   }
0407
0408   VECCORE_ATT_HOST_DEVICE
0409   PolyconeStruct() {}
0410 };
0411 } // namespace VECGEOM_IMPL_NAMESPACE
0412 } // namespace vecgeom
0413
0414 #endif