Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-30 10:26:21 

0001 /// \file SecondOrderSurfaceShell.h
0002 /// \author: swenzel
0003 ///  Created on: Aug 1, 2014
0004 ///  Modified and completed: mihaela.gheata@cern.ch
0005 
0006 #ifndef VECGEOM_SECONDORDERSURFACESHELL_H_
0007 #define VECGEOM_SECONDORDERSURFACESHELL_H_
0008 
0009 #include "VecGeom/base/Vector3D.h"
0010 #include "VecGeom/volumes/kernel/GenericKernels.h"
0011 
0012 namespace vecgeom {
0013 /**
0014  * Templated class providing a (SOA) encapsulation of
0015  * second order curved surfaces used in generic trap
0016  */
0017 VECGEOM_DEVICE_FORWARD_DECLARE(class SecondOrderSurfaceShell;);
0018 
0019 inline namespace VECGEOM_IMPL_NAMESPACE {
0020 
0021 template <int N>
0022 class SecondOrderSurfaceShell {
0023 
0024   using Vertex_t = Vector3D<Precision>;
0025 
0026 private:
0027   // caching some important values for each of the curved planes
0028   Precision fxa[N], fya[N], fxb[N], fyb[N], fxc[N], fyc[N], fxd[N], fyd[N]; /** Coordinates of vertices */
0029   Precision ftx1[N], fty1[N], ftx2[N], fty2[N];                             /** Connecting components */
0030   Precision ft1crosst2[N]; /** Cross term ftx1[i]*fty2[i] - ftx2[i]*fty1[i] */
0031   Precision fDeltatx[N];   /** Term ftx2[i] - ftx1[i] */
0032   Precision fDeltaty[N];   /** Term fty2[i] - fty1[i] */
0033 
0034   Precision fDz;          /** height of surface (coming from the height of the GenTrap) */
0035   Precision fDz2;         /** 0.5/fDz */
0036   Precision fiscurved[N]; /** Indicate which surfaces are planar */
0037   bool fisplanar;         /** Flag that all surfaces are planar */
0038   bool fdegenerated[N];   /** Flags for each top-bottom edge marking that this is degenerated */
0039 
0040   Vertex_t fNormals[N]; /** Pre-computed normals for the planar case */
0041 
0042   // pre-computed cross products for normal computation
0043   Vertex_t fViCrossHi0[N];  /** Pre-computed vi X hi0 */
0044   Vertex_t fViCrossVj[N];   /** Pre-computed vi X vj */
0045   Vertex_t fHi1CrossHi0[N]; /** Pre-computed hi1 X hi0 */
0046 
0047 public:
0048   /** @brief SecondOrderSurfaceShell dummy constructor */
0049   VECCORE_ATT_HOST_DEVICE
0050   SecondOrderSurfaceShell() : fDz(0), fDz2(0) {}
0051 
0052   /** @brief SecondOrderSurfaceShell constructor
0053    * @param verticesx X positions of vertices in array form
0054    * @param verticesy Y positions of vertices in array form
0055    * @param dz The half-height of the GenTrap
0056    */
0057   VECCORE_ATT_HOST_DEVICE
0058   SecondOrderSurfaceShell(const Precision *verticesx, const Precision *verticesy, Precision dz) : fDz(0.), fDz2(0.)
0059   {
0060     // Constructor
0061     Initialize(verticesx, verticesy, dz);
0062   }
0063 
0064   VECCORE_ATT_HOST_DEVICE
0065   void Initialize(const Precision *verticesx, const Precision *verticesy, Precision dz)
0066   {
0067 #ifndef VECCORE_CUDA
0068     if (dz <= 0.) throw std::runtime_error("Half-length of generic trapezoid must be positive");
0069 #endif
0070     fDz  = dz;
0071     fDz2 = 0.5 / dz;
0072     // Store vertex coordinates
0073     Vertex_t va, vb, vc, vd;
0074     for (int i = 0; i < N; ++i) {
0075       int j = (i + 1) % N;
0076       va.Set(verticesx[i], verticesy[i], -dz);
0077       fxa[i] = verticesx[i];
0078       fya[i] = verticesy[i];
0079       vb.Set(verticesx[i + N], verticesy[i + N], dz);
0080       fxb[i] = verticesx[i + N];
0081       fyb[i] = verticesy[i + N];
0082       vc.Set(verticesx[j], verticesy[j], -dz);
0083       fxc[i] = verticesx[j];
0084       fyc[i] = verticesy[j];
0085       vd.Set(verticesx[j + N], verticesy[j + N], dz);
0086       fxd[i]          = verticesx[N + j];
0087       fyd[i]          = verticesy[N + j];
0088       fdegenerated[i] = (Abs(fxa[i] - fxc[i]) < kTolerance) && (Abs(fya[i] - fyc[i]) < kTolerance) &&
0089                         (Abs(fxb[i] - fxd[i]) < kTolerance) && (Abs(fyb[i] - fyd[i]) < kTolerance);
0090       ftx1[i] = fDz2 * (fxb[i] - fxa[i]);
0091       fty1[i] = fDz2 * (fyb[i] - fya[i]);
0092       ftx2[i] = fDz2 * (fxd[i] - fxc[i]);
0093       fty2[i] = fDz2 * (fyd[i] - fyc[i]);
0094 
0095       ft1crosst2[i] = ftx1[i] * fty2[i] - ftx2[i] * fty1[i];
0096       fDeltatx[i]   = ftx2[i] - ftx1[i];
0097       fDeltaty[i]   = fty2[i] - fty1[i];
0098       fNormals[i]   = Vertex_t::Cross(vb - va, vc - va);
0099       // The computation of normals is done also for the curved surfaces, even if they will not be used
0100       if (fNormals[i].Mag2() < kTolerance) {
0101         // points i and i+1/N are overlapping - use i+N and j+N instead
0102         fNormals[i] = Vertex_t::Cross(vb - va, vd - vb);
0103         if (fNormals[i].Mag2() < kTolerance) fNormals[i].Set(0., 0., 1.); // No surface, just a line
0104       }
0105       fNormals[i].Normalize();
0106       // Cross products used for normal computation
0107       fViCrossHi0[i]  = Vertex_t::Cross(vb - va, vc - va);
0108       fViCrossVj[i]   = Vertex_t::Cross(vb - va, vd - vc);
0109       fHi1CrossHi0[i] = Vertex_t::Cross(vd - vb, vc - va);
0110     }
0111 
0112     // Analyze planarity and precompute normals
0113     fisplanar = true;
0114     for (int i = 0; i < N; ++i) {
0115       fiscurved[i] = (((Abs(fxc[i] - fxa[i]) < kTolerance) && (Abs(fyc[i] - fya[i]) < kTolerance)) ||
0116                       ((Abs(fxd[i] - fxb[i]) < kTolerance) && (Abs(fyd[i] - fyb[i]) < kTolerance)) ||
0117                       (Abs((fxc[i] - fxa[i]) * (fyd[i] - fyb[i]) - (fxd[i] - fxb[i]) * (fyc[i] - fya[i])) < kTolerance))
0118                          ? 0
0119                          : 1;
0120       if (fiscurved[i]) fisplanar = false;
0121     }
0122   }
0123 
0124   //______________________________________________________________________________
0125   template <typename Real_v>
0126   VECGEOM_FORCE_INLINE
0127   VECCORE_ATT_HOST_DEVICE
0128   /** @brief Stripped down version of Inside method used for boundary and wrong side detection
0129    * @param point Starting point in the local frame
0130    * @param completelyinside Inside flag
0131    * @param completelyoutside Onside flag
0132    * @param onsurf On boundary flag
0133    */
0134   void CheckInside(Vector3D<Real_v> const &point, vecCore::Mask_v<Real_v> &completelyinside,
0135                    vecCore::Mask_v<Real_v> &completelyoutside, vecCore::Mask_v<Real_v> &onsurf) const
0136   {
0137 
0138     using Bool_v                        = vecCore::Mask_v<Real_v>;
0139     VECGEOM_CONST Precision tolerancesq = 10000. * kTolerance * kTolerance;
0140 
0141     onsurf            = Bool_v(false);
0142     completelyinside  = (Abs(point.z()) < Real_v(MakeMinusTolerant<true>(fDz)));
0143     completelyoutside = (Abs(point.z()) > Real_v(MakePlusTolerant<true>(fDz)));
0144     //  if (vecCore::EarlyReturnAllowed()) {
0145     if (vecCore::MaskFull(completelyoutside)) return;
0146     //  }
0147 
0148     Real_v cross;
0149     Real_v vertexX[N];
0150     Real_v vertexY[N];
0151     Real_v dzp = fDz + point[2];
0152     // vectorizes for scalar backend
0153     for (int i = 0; i < N; i++) {
0154       // calculate x-y positions of vertex i at this z-height
0155       vertexX[i] = fxa[i] + ftx1[i] * dzp;
0156       vertexY[i] = fya[i] + fty1[i] * dzp;
0157     }
0158     Bool_v degenerated = Bool_v(true); // Can only happen if |zpoint| = fDz
0159     for (int i = 0; i < N; i++) {
0160       //      if (fdegenerated[i])
0161       //        continue;
0162       int j          = (i + 1) % 4;
0163       Real_v DeltaX  = vertexX[j] - vertexX[i];
0164       Real_v DeltaY  = vertexY[j] - vertexY[i];
0165       Real_v deltasq = DeltaX * DeltaX + DeltaY * DeltaY;
0166       // If the current vertex is degenerated, ignore the check
0167       Bool_v samevertex = deltasq < Real_v(MakePlusTolerant<true>(0.));
0168       degenerated       = degenerated && samevertex;
0169       // Cross product to check if point is right side or left side
0170       // If vertices are same, this will be 0
0171       cross = (point.x() - vertexX[i]) * DeltaY - (point.y() - vertexY[i]) * DeltaX;
0172 
0173       onsurf            = (cross * cross < tolerancesq * deltasq) && (!samevertex);
0174       completelyoutside = completelyoutside || (((cross < Real_v(MakeMinusTolerant<true>(0.))) && (!onsurf)));
0175       completelyinside =
0176           completelyinside && (samevertex || ((cross > Real_v(MakePlusTolerant<true>(0.))) && (!onsurf)));
0177     }
0178     onsurf = (!completelyoutside) && (!completelyinside);
0179     // In fully degenerated case consider the point outside always
0180     completelyoutside = completelyoutside || degenerated;
0181   }
0182 
0183   //______________________________________________________________________________
0184   template <typename Real_v>
0185   VECGEOM_FORCE_INLINE
0186   VECCORE_ATT_HOST_DEVICE
0187   /** @brief Compute distance to a set of curved/planar surfaces
0188    * @param point Starting point in the local frame
0189    * @param dir Direction in the local frame
0190    */
0191   Real_v DistanceToOut(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir) const
0192   {
0193 
0194     using Bool_v = vecCore::Mask_v<Real_v>;
0195 
0196     // Planar case
0197     if (fisplanar) return DistanceToOutPlanar<Real_v>(point, dir);
0198 
0199     // The algorithmic tolerance in distance
0200     const Real_v tolerance(100. * kTolerance);
0201 
0202     Real_v dist(InfinityLength<Real_v>());
0203     Real_v smin[N], smax[N];
0204     Vector3D<Real_v> unorm;
0205     Real_v r(-1.0);
0206     Real_v rz;
0207     Bool_v completelyinside, completelyoutside, onsurf;
0208     CheckInside<Real_v>(point, completelyinside, completelyoutside, onsurf);
0209 
0210     // If on the wrong side, return -1.
0211     Real_v wrongsidedist(-1.0);
0212     vecCore::MaskedAssign(dist, completelyoutside, wrongsidedist);
0213     if (vecCore::MaskFull(completelyoutside)) return dist;
0214 
0215     // Solve the second order equation and return distance solutions for each surface
0216     ComputeSminSmax<Real_v>(point, dir, smin, smax);
0217     for (int i = 0; i < N; ++i) {
0218       // Check if point(s) is(are) on boundary, and in this case compute normal
0219       Bool_v crtbound = (Abs(smin[i]) < tolerance || Abs(smax[i]) < tolerance);
0220       if (!vecCore::MaskEmpty(crtbound)) {
0221         if (fiscurved[i])
0222           UNormal<Real_v>(point, i, unorm, rz, r);
0223         else
0224           unorm = fNormals[i];
0225       }
0226       // Starting point may be propagated close to boundary
0227       // === MaskedMultipleAssign needed
0228       vecCore__MaskedAssignFunc(smin[i], !completelyoutside && Abs(smin[i]) < tolerance && dir.Dot(unorm) < 0,
0229                                 InfinityLength<Real_v>());
0230       vecCore__MaskedAssignFunc(smax[i], !completelyoutside && Abs(smax[i]) < tolerance && dir.Dot(unorm) < 0,
0231                                 InfinityLength<Real_v>());
0232 
0233       vecCore__MaskedAssignFunc(dist, !completelyoutside && (smin[i] > -tolerance) && (smin[i] < dist),
0234                                 Max(smin[i], Real_v(0.)));
0235       vecCore__MaskedAssignFunc(dist, !completelyoutside && (smax[i] > -tolerance) && (smax[i] < dist),
0236                                 Max(smax[i], Real_v(0.)));
0237     }
0238     vecCore__MaskedAssignFunc(dist, dist < tolerance && onsurf, Real_v(0.));
0239     return (dist);
0240   } // end of function
0241 
0242   //______________________________________________________________________________
0243   /** @brief Compute distance to exiting the set of surfaces in the planar case.
0244    * @param point Starting point in the local frame
0245    * @param dir Direction in the local frame
0246    */
0247   template <typename Real_v>
0248   VECGEOM_FORCE_INLINE
0249   VECCORE_ATT_HOST_DEVICE
0250   Real_v DistanceToOutPlanar(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir) const
0251   {
0252 
0253     using Bool_v = vecCore::Mask_v<Real_v>;
0254     //    const Real_v tolerance = 100. * kTolerance;
0255 
0256     Vertex_t va;         // vertex i of lower base
0257     Vector3D<Real_v> pa; // same vertex converted to backend type
0258     Real_v distance = InfinityLength<Real_v>();
0259 
0260     // Check every surface
0261     Bool_v outside = (Abs(point.z()) > MakePlusTolerant<true>(fDz)); // If point is outside, we need to know
0262     for (int i = 0; i < N && (!vecCore::MaskFull(outside)); ++i) {
0263       if (fdegenerated[i]) continue;
0264       // Point A is the current vertex on lower Z. P is the point we come from.
0265       pa.Set(Real_v(fxa[i]), Real_v(fya[i]), Real_v(-fDz));
0266       Vector3D<Real_v> vecAP = point - pa;
0267       // Dot product between AP vector and normal to surface has to be negative
0268       Real_v dotAPNorm = vecAP.Dot(fNormals[i]);
0269       Bool_v otherside = (dotAPNorm > Real_v(10.) * kTolerance);
0270       // Dot product between direction and normal to surface has to be positive
0271       Real_v dotDirNorm = dir.Dot(fNormals[i]);
0272       Bool_v outgoing   = (dotDirNorm > Real_v(0.));
0273       // Update globally outside flag
0274       outside      = outside || otherside;
0275       Bool_v valid = outgoing && (!otherside);
0276       if (vecCore::MaskEmpty(valid)) continue;
0277       Real_v snext = -dotAPNorm / NonZero(dotDirNorm);
0278       vecCore__MaskedAssignFunc(distance, valid && snext < distance, Max(snext, Real_v(0.)));
0279     }
0280     // Return -1 for points actually outside
0281     vecCore__MaskedAssignFunc(distance, distance < kTolerance, Real_v(0.));
0282     vecCore__MaskedAssignFunc(distance, outside, Real_v(-1.));
0283     return distance;
0284   }
0285 
0286   //______________________________________________________________________________
0287   /** @brief Compute distance to entering the set of surfaces in the planar case.
0288    * @param point Starting point in the local frame
0289    * @param dir Direction in the local frame
0290    */
0291   template <typename Real_v>
0292   VECGEOM_FORCE_INLINE
0293   VECCORE_ATT_HOST_DEVICE
0294   Real_v DistanceToInPlanar(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir,
0295                             vecCore::Mask_v<Real_v> &done) const
0296   {
0297     //    const Real_v tolerance = 100. * kTolerance;
0298 
0299     using Bool_v = vecCore::Mask_v<Real_v>;
0300     Vertex_t va;         // vertex i of lower base
0301     Vector3D<Real_v> pa; // same vertex converted to backend type
0302     Real_v distance = InfinityLength<Real_v>();
0303 
0304     // Check every surface
0305     Bool_v inside = (Abs(point.z()) < MakeMinusTolerant<true>(fDz)); // If point is inside, we need to know
0306     for (int i = 0; i < N; ++i) {
0307       if (fdegenerated[i]) continue;
0308       // Point A is the current vertex on lower Z. P is the point we come from.
0309       pa.Set(Real_v(fxa[i]), Real_v(fya[i]), Real_v(-fDz));
0310       Vector3D<Real_v> vecAP = point - pa;
0311       // Dot product between AP vector and normal to surface has to be positive
0312       Real_v dotAPNorm = vecAP.Dot(fNormals[i]);
0313       Bool_v otherside = (dotAPNorm < Real_v(-10.) * kTolerance);
0314       // Dot product between direction and normal to surface has to be negative
0315       Real_v dotDirNorm = dir.Dot(fNormals[i]);
0316       Bool_v ingoing    = (dotDirNorm < Real_v(0.));
0317       // Update globally outside flag
0318       inside       = inside && otherside;
0319       Bool_v valid = ingoing && (!otherside);
0320       if (vecCore::MaskEmpty(valid)) continue;
0321       Real_v snext = -dotAPNorm / NonZero(dotDirNorm);
0322       // Now propagate the point to surface and check if in range
0323       Vector3D<Real_v> psurf = point + snext * dir;
0324       valid                  = valid && InSurfLimits<Real_v>(psurf, i);
0325       vecCore__MaskedAssignFunc(distance, (!done) && valid && snext < distance, Max(snext, Real_v(0.)));
0326     }
0327     // Return -1 for points actually inside
0328     vecCore__MaskedAssignFunc(distance, (!done) && (distance < kTolerance), Real_v(0.));
0329     vecCore__MaskedAssignFunc(distance, (!done) && inside, Real_v(-1.));
0330     return distance;
0331   }
0332 
0333   //______________________________________________________________________________
0334   template <typename Real_v>
0335   VECGEOM_FORCE_INLINE
0336   VECCORE_ATT_HOST_DEVICE
0337   /**
0338    * A generic function calculation the distance to a set of curved/planar surfaces
0339    *
0340    * things to improve: error handling for boundary cases
0341    *
0342    * another possibility relies on the following idea:
0343    * we always have an even number of planar/curved surfaces. We could organize them in separate substructures...
0344    */
0345   Real_v DistanceToIn(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir, vecCore::Mask_v<Real_v> &done) const
0346   {
0347     // Planar case
0348     using Bool_v = vecCore::Mask_v<Real_v>;
0349     if (fisplanar) return DistanceToInPlanar<Real_v>(point, dir, done);
0350     Real_v crtdist;
0351     Vector3D<Real_v> hit;
0352     Real_v distance = InfinityLength<Real_v>();
0353     Real_v tolerance(100. * kTolerance);
0354     Vector3D<Real_v> unorm;
0355     Real_v r(-1.0);
0356     Real_v rz;
0357     Bool_v completelyinside, completelyoutside, onsurf;
0358     CheckInside<Real_v>(point, completelyinside, completelyoutside, onsurf);
0359     // If on the wrong side, return -1.
0360     Real_v wrongsidedist(-1.0);
0361     vecCore::MaskedAssign(distance, Bool_v(completelyinside && (!done)), wrongsidedist);
0362     Bool_v checked = completelyinside || done;
0363     if (vecCore::MaskFull(checked)) return (distance);
0364     // Now solve the second degree equation to find crossings
0365     Real_v smin[N], smax[N];
0366     ComputeSminSmax<Real_v>(point, dir, smin, smax);
0367 
0368     // now we need to analyse which of those distances is good
0369     // does not vectorize
0370     for (int i = 0; i < N; ++i) {
0371       crtdist = smin[i];
0372       // Extrapolate with hit distance candidate
0373       hit             = point + crtdist * dir;
0374       Bool_v crossing = (crtdist > -tolerance) && (Abs(hit.z()) < fDz + kTolerance);
0375       // Early skip surface if not in Z range
0376       if (!vecCore::MaskEmpty(Bool_v(crossing && (!checked)))) {
0377         // Compute local un-normalized outwards normal direction and hit ratio factors
0378         UNormal<Real_v>(hit, i, unorm, rz, r);
0379         // Distance have to be positive within tolerance, and crossing must be inwards
0380         crossing = crossing && dir.Dot(unorm) < Real_v(0.);
0381         // Propagated hitpoint must be on surface (rz in [0,1] checked already)
0382         crossing = crossing && (r >= Real_v(0.)) && (r <= Real_v(1.));
0383         vecCore__MaskedAssignFunc(distance, crossing && (!checked) && crtdist < distance, Max(crtdist, Real_v(0.)));
0384       }
0385       // For the particle(s) not crossing at smin, try smax
0386       if (!vecCore::MaskFull(Bool_v(crossing || checked))) {
0387         // Treat only particles not crossing at smin
0388         crossing = !crossing;
0389         crtdist  = smax[i];
0390         hit      = point + crtdist * dir;
0391         crossing = crossing && (crtdist > -tolerance) && (Abs(hit.z()) < fDz + kTolerance);
0392         if (vecCore::MaskEmpty(crossing)) continue;
0393         if (fiscurved[i])
0394           UNormal<Real_v>(hit, i, unorm, rz, r);
0395         else
0396           unorm = fNormals[i];
0397         crossing = crossing && (dir.Dot(unorm) < Real_v(0.));
0398         crossing = crossing && (r >= Real_v(0.)) && (r <= Real_v(1.));
0399         vecCore__MaskedAssignFunc(distance, crossing && (!checked) && crtdist < distance, Max(crtdist, Real_v(0.)));
0400       }
0401     }
0402     vecCore__MaskedAssignFunc(distance, distance < tolerance && onsurf, Real_v(0.));
0403     return (distance);
0404 
0405   } // end distanceToIn function
0406 
0407   //______________________________________________________________________________
0408   /**
0409    * @brief A generic function calculation for the safety to a set of curved/planar surfaces.
0410            Should be smaller than safmax
0411    * @param point Starting point inside, in the local frame
0412    * @param safmax current safety value
0413    */
0414   template <typename Real_v>
0415   VECGEOM_FORCE_INLINE
0416   VECCORE_ATT_HOST_DEVICE
0417   Real_v SafetyToOut(Vector3D<Real_v> const &point, Real_v const &safmax) const
0418   {
0419 
0420     using Bool_v                = vecCore::Mask_v<Real_v>;
0421     VECGEOM_CONST Precision eps = 100. * kTolerance;
0422 
0423     Real_v safety = safmax;
0424     Bool_v done   = (Abs(safety) < eps);
0425     if (vecCore::MaskFull(done)) return (safety);
0426     Real_v safetyface = InfinityLength<Real_v>();
0427 
0428     // loop lateral surfaces
0429     // We can use the surface normals to get safety for non-curved surfaces
0430     Vertex_t va;         // vertex i of lower base
0431     Vector3D<Real_v> pa; // same vertex converted to backend type
0432     if (fisplanar) {
0433       for (int i = 0; i < N; ++i) {
0434         if (fdegenerated[i]) continue;
0435         va.Set(fxa[i], fya[i], -fDz);
0436         pa         = va;
0437         safetyface = (pa - point).Dot(fNormals[i]);
0438         vecCore::MaskedAssign(safety, (safetyface < safety) && (!done), safetyface);
0439       }
0440       vecCore__MaskedAssignFunc(safety, Abs(safety) < eps, Real_v(0.));
0441       return safety;
0442     }
0443 
0444     // Not fully planar - use mixed case
0445     safetyface = SafetyCurved<Real_v>(point, Bool_v(true));
0446     //  std::cout << "safetycurved = " << safetyface << std::endl;
0447     vecCore::MaskedAssign(safety, (safetyface < safety) && (!done), safetyface);
0448     //  std::cout << "safety = " << safety << std::endl;
0449     vecCore__MaskedAssignFunc(safety, safety > Real_v(0.) && safety < eps, Real_v(0.));
0450     return safety;
0451 
0452   } // end SafetyToOut
0453 
0454   //______________________________________________________________________________
0455   /**
0456    * @brief A generic function calculation for the safety to a set of curved/planar surfaces.
0457            Should be smaller than safmax
0458    * @param point Starting point outside, in the local frame
0459    * @param safmax current safety value
0460    */
0461   template <typename Real_v>
0462   VECGEOM_FORCE_INLINE
0463   VECCORE_ATT_HOST_DEVICE
0464   Real_v SafetyToIn(Vector3D<Real_v> const &point, Real_v const &safmax) const
0465   {
0466 
0467     using Bool_v                = vecCore::Mask_v<Real_v>;
0468     VECGEOM_CONST Precision eps = 100. * kTolerance;
0469 
0470     Real_v safety = safmax;
0471     Bool_v done   = (Abs(safety) < eps);
0472     if (vecCore::MaskFull(done)) return (safety);
0473     Real_v safetyface = InfinityLength<Real_v>();
0474 
0475     // loop lateral surfaces
0476     // We can use the surface normals to get safety for non-curved surfaces
0477     Vertex_t va;         // vertex i of lower base
0478     Vector3D<Real_v> pa; // same vertex converted to backend type
0479     if (fisplanar) {
0480       for (int i = 0; i < N; ++i) {
0481         if (fdegenerated[i]) continue;
0482         va.Set(fxa[i], fya[i], -fDz);
0483         pa         = va;
0484         safetyface = (point - pa).Dot(fNormals[i]);
0485         vecCore::MaskedAssign(safety, (safetyface > safety) && (!done), safetyface);
0486       }
0487       vecCore__MaskedAssignFunc(safety, Abs(safety) < eps, Real_v(0.));
0488       return safety;
0489     }
0490 
0491     // Not fully planar - use mixed case
0492     safetyface = SafetyCurved<Real_v>(point, Bool_v(false));
0493     //  std::cout << "safetycurved = " << safetyface << std::endl;
0494     vecCore::MaskedAssign(safety, (safetyface > safety) && (!done), safetyface);
0495     //  std::cout << "safety = " << safety << std::endl;
0496     vecCore__MaskedAssignFunc(safety, safety > Real_v(0.) && safety < eps, Real_v(0.));
0497     return (safety);
0498 
0499   } // end SafetyToIn
0500 
0501   //______________________________________________________________________________
0502   /**
0503    * @brief A generic function calculation for the safety to a set of curved/planar surfaces.
0504            Should be smaller than safmax
0505    * @param point Starting point
0506    * @param in Inside value for starting point
0507    */
0508   template <typename Real_v>
0509   VECGEOM_FORCE_INLINE
0510   VECCORE_ATT_HOST_DEVICE
0511   Real_v SafetyCurved(Vector3D<Real_v> const &point, vecCore::Mask_v<Real_v> in) const
0512   {
0513     using Bool_v     = vecCore::Mask_v<Real_v>;
0514     Real_v safety    = InfinityLength<Real_v>();
0515     Real_v safplanar = InfinityLength<Real_v>();
0516     Real_v tolerance = Real_v(100 * kTolerance);
0517     vecCore__MaskedAssignFunc(tolerance, !in, -tolerance);
0518 
0519     //  loop over edges connecting points i with i+4
0520     Real_v dx, dy, dpx, dpy, lsq, u;
0521     Real_v dx1 = Real_v(0.);
0522     Real_v dx2 = Real_v(0.);
0523     Real_v dy1 = Real_v(0.);
0524     Real_v dy2 = Real_v(0.);
0525 
0526     Bool_v completelyinside, completelyoutside, onsurf;
0527     CheckInside<Real_v>(point, completelyinside, completelyoutside, onsurf);
0528     if (vecCore::MaskFull(onsurf)) return (Real_v(0.));
0529 
0530     Bool_v wrong = in && (completelyoutside);
0531     wrong        = wrong || ((!in) && completelyinside);
0532     if (vecCore::MaskFull(wrong)) {
0533       return (Real_v(-1.));
0534     }
0535     Real_v vertexX[N];
0536     Real_v vertexY[N];
0537     Real_v dzp = fDz + point[2];
0538     for (int i = 0; i < N; i++) {
0539       // calculate x-y positions of vertex i at this z-height
0540       vertexX[i] = fxa[i] + ftx1[i] * dzp;
0541       vertexY[i] = fya[i] + fty1[i] * dzp;
0542     }
0543     Real_v umin = Real_v(0.);
0544     for (int i = 0; i < N; i++) {
0545       if (fiscurved[i] == 0) {
0546         // We can use the surface normals to get safety for non-curved surfaces
0547         Vector3D<Real_v> pa; // same vertex converted to backend type
0548         pa.Set(Real_v(fxa[i]), Real_v(fya[i]), Real_v(-fDz));
0549         Real_v sface = Abs((point - pa).Dot(fNormals[i]));
0550         vecCore::MaskedAssign(safplanar, sface < safplanar, sface);
0551         continue;
0552       }
0553       int j = (i + 1) % N;
0554       dx    = vertexX[j] - vertexX[i];
0555       dy    = vertexY[j] - vertexY[i];
0556       dpx   = point[0] - vertexX[i];
0557       dpy   = point[1] - vertexY[i];
0558       lsq   = dx * dx + dy * dy;
0559       u     = (dpx * dx + dpy * dy) / NonZero(lsq);
0560       vecCore__MaskedAssignFunc(dpx, u > 1, point[0] - vertexX[j]);
0561       vecCore__MaskedAssignFunc(dpy, u > 1, point[1] - vertexY[j]);
0562       vecCore__MaskedAssignFunc(dpx, u >= 0 && u <= 1, dpx - u * dx);
0563       vecCore__MaskedAssignFunc(dpy, u >= 0 && u <= 1, dpy - u * dy);
0564       Real_v ssq = dpx * dpx + dpy * dpy; // safety squared
0565       vecCore__MaskedAssignFunc(dx1, ssq < safety, Real_v(fxc[i] - fxa[i]));
0566       vecCore__MaskedAssignFunc(dx2, ssq < safety, Real_v(fxd[i] - fxb[i]));
0567       vecCore__MaskedAssignFunc(dy1, ssq < safety, Real_v(fyc[i] - fya[i]));
0568       vecCore__MaskedAssignFunc(dy2, ssq < safety, Real_v(fyd[i] - fyb[i]));
0569       vecCore::MaskedAssign(umin, ssq < safety, u);
0570       vecCore::MaskedAssign(safety, ssq < safety, ssq);
0571     }
0572     vecCore__MaskedAssignFunc(umin, umin < 0 || umin > 1, Real_v(0.));
0573     dx = dx1 + umin * (dx2 - dx1);
0574     dy = dy1 + umin * (dy2 - dy1);
0575     // Denominator below always positive as fDz>0
0576     safety *= Real_v(1.) - Real_v(4.) * fDz * fDz / (dx * dx + dy * dy + Real_v(4.) * fDz * fDz);
0577     safety = Sqrt(safety);
0578     vecCore::MaskedAssign(safety, safplanar < safety, safplanar);
0579     vecCore__MaskedAssignFunc(safety, wrong, -safety);
0580     vecCore__MaskedAssignFunc(safety, onsurf, Real_v(0.));
0581     return safety;
0582   } // end SafetyCurved
0583 
0584   VECCORE_ATT_HOST_DEVICE
0585   VECGEOM_FORCE_INLINE
0586   Vertex_t const *GetNormals() const { return fNormals; }
0587 
0588   //______________________________________________________________________________
0589   /**
0590    * @brief Computes if point on surface isurf is within surface limits
0591    * @param point Starting point
0592    * @param isurf Surface index
0593    */
0594   template <typename Real_v>
0595   VECGEOM_FORCE_INLINE
0596   VECCORE_ATT_HOST_DEVICE
0597   vecCore::Mask_v<Real_v> InSurfLimits(Vector3D<Real_v> const &point, int isurf) const
0598   {
0599 
0600     using Bool_v = vecCore::Mask_v<Real_v>;
0601     // Check first if Z is in range
0602     Real_v rz     = fDz2 * (point.z() + fDz);
0603     Bool_v insurf = (rz > Real_v(MakeMinusTolerant<true>(0.))) && (rz < Real_v(MakePlusTolerant<true>(1.)));
0604     if (vecCore::MaskEmpty(insurf)) return insurf;
0605 
0606     Real_v r     = InfinityLength<Real_v>();
0607     Real_v num   = (point.x() - fxa[isurf]) - rz * (fxb[isurf] - fxa[isurf]);
0608     Real_v denom = (fxc[isurf] - fxa[isurf]) + rz * (fxd[isurf] - fxc[isurf] - fxb[isurf] + fxa[isurf]);
0609     vecCore__MaskedAssignFunc(r, (Abs(denom) > Real_v(1.e-6)), num / NonZero(denom));
0610     num   = (point.y() - fya[isurf]) - rz * (fyb[isurf] - fya[isurf]);
0611     denom = (fyc[isurf] - fya[isurf]) + rz * (fyd[isurf] - fyc[isurf] - fyb[isurf] + fya[isurf]);
0612     vecCore__MaskedAssignFunc(r, (Abs(denom) > Real_v(1.e-6)), num / NonZero(denom));
0613     insurf = insurf && (r > Real_v(MakeMinusTolerant<true>(0.))) && (r < Real_v(MakePlusTolerant<true>(1.)));
0614     return insurf;
0615   }
0616 
0617   //______________________________________________________________________________
0618   /** @brief Computes un-normalized normal to surface isurf, on the input point */
0619   template <typename Real_v>
0620   VECGEOM_FORCE_INLINE
0621   VECCORE_ATT_HOST_DEVICE
0622   void UNormal(Vector3D<Real_v> const &point, int isurf, Vector3D<Real_v> &unorm, Real_v &rz, Real_v &r) const
0623   {
0624 
0625     // unorm = (vi X hi0) + rz*(vi X vj) + r*(hi1 X hi0)
0626     //    where: vi, vj are the vectors (AB) and (CD) (see constructor)
0627     //           hi0 = (AC) and hi1 = (BD)
0628     //           rz = 0.5*(point.z()+dz)/dz is the vertical ratio
0629     //           r = ((AP)-rz*vi) / (hi0+rz(vj-vi)) is the horizontal ratio
0630     // Any point within the surface range should reurn r and rz in the range [0,1]
0631     // These can be used as surface crossing criteria
0632     rz = fDz2 * (point.z() + fDz);
0633     /*
0634       Vector3D<Real_v> a(fxa[isurf], fya[isurf], -fDz);
0635       Vector3D<Real_v> vi(fxb[isurf]-fxa[isurf], fyb[isurf]-fya[isurf], 2*fDz);
0636       Vector3D<Real_v> vj(fxd[isurf]-fxc[isurf], fyd[isurf]-fyc[isurf], 2*fDz);
0637       Vector3D<Real_v> hi0(fxc[isurf]-fxa[isurf], fyc[isurf]-fya[isurf], 0.);
0638     */
0639     Real_v num   = (point.x() - fxa[isurf]) - rz * (fxb[isurf] - fxa[isurf]);
0640     Real_v denom = (fxc[isurf] - fxa[isurf]) + rz * (fxd[isurf] - fxc[isurf] - fxb[isurf] + fxa[isurf]);
0641     vecCore__MaskedAssignFunc(r, Abs(denom) > Real_v(1.e-6), num / NonZero(denom));
0642     num   = (point.y() - fya[isurf]) - rz * (fyb[isurf] - fya[isurf]);
0643     denom = (fyc[isurf] - fya[isurf]) + rz * (fyd[isurf] - fyc[isurf] - fyb[isurf] + fya[isurf]);
0644     vecCore__MaskedAssignFunc(r, Abs(denom) > Real_v(1.e-6), num / NonZero(denom));
0645 
0646     unorm = (Vector3D<Real_v>)fViCrossHi0[isurf] + rz * (Vector3D<Real_v>)fViCrossVj[isurf] +
0647             r * (Vector3D<Real_v>)fHi1CrossHi0[isurf];
0648   } // end UNormal
0649 
0650   //______________________________________________________________________________
0651   /** @brief Solver for the second degree equation for curved surface crossing */
0652   template <typename Real_v>
0653   VECGEOM_FORCE_INLINE
0654   VECCORE_ATT_HOST_DEVICE
0655   /**
0656    * Function to compute smin and smax crossings with the N lateral surfaces.
0657    */
0658   void ComputeSminSmax(Vector3D<Real_v> const &point, Vector3D<Real_v> const &dir, Real_v smin[N], Real_v smax[N]) const
0659   {
0660 
0661     const Real_v big = InfinityLength<Real_v>();
0662 
0663     Real_v dzp = fDz + point[2];
0664     // calculate everything needed to solve the second order equation
0665     Real_v a[N], b[N], c[N], d[N];
0666     Real_v signa[N], inva[N];
0667 
0668     // vectorizes
0669     for (int i = 0; i < N; ++i) {
0670       Real_v xs1 = fxa[i] + ftx1[i] * dzp;
0671       Real_v ys1 = fya[i] + fty1[i] * dzp;
0672       Real_v xs2 = fxc[i] + ftx2[i] * dzp;
0673       Real_v ys2 = fyc[i] + fty2[i] * dzp;
0674       Real_v dxs = xs2 - xs1;
0675       Real_v dys = ys2 - ys1;
0676       a[i]       = (fDeltatx[i] * dir[1] - fDeltaty[i] * dir[0] + ft1crosst2[i] * dir[2]) * dir[2];
0677       b[i]       = dxs * dir[1] - dys * dir[0] +
0678              (fDeltatx[i] * point[1] - fDeltaty[i] * point[0] + fty2[i] * xs1 - fty1[i] * xs2 + ftx1[i] * ys2 -
0679               ftx2[i] * ys1) *
0680                  dir[2];
0681       c[i] = dxs * point[1] - dys * point[0] + xs1 * ys2 - xs2 * ys1;
0682       d[i] = b[i] * b[i] - 4 * a[i] * c[i];
0683     }
0684 
0685     // does not vectorize
0686     for (int i = 0; i < N; ++i) {
0687       // zero or one to start with
0688       signa[i] = Real_v(0.);
0689       vecCore__MaskedAssignFunc(signa[i], a[i] < -kTolerance, Real_v(-1.));
0690       vecCore__MaskedAssignFunc(signa[i], a[i] > kTolerance, Real_v(1.));
0691       inva[i] = c[i] / NonZero(b[i] * b[i]);
0692       vecCore__MaskedAssignFunc(inva[i], Abs(a[i]) > kTolerance, Real_v(1.) / NonZero(Real_v(2.) * a[i]));
0693     }
0694 
0695     // vectorizes
0696     for (int i = 0; i < N; ++i) {
0697       // treatment for curved surfaces. Invalid solutions will be excluded.
0698       Real_v sqrtd = signa[i] * Sqrt(Abs(d[i]));
0699       vecCore::MaskedAssign(sqrtd, d[i] < Real_v(0.), big);
0700       // what is the meaning of this??
0701       smin[i] = (-b[i] - sqrtd) * inva[i];
0702       smax[i] = (-b[i] + sqrtd) * inva[i];
0703     }
0704     // For the planar surfaces, the above may be wrong, redo the work using
0705     // just the normal. This does not vectorize
0706     for (int i = 0; i < N; ++i) {
0707       if (fiscurved[i]) continue;
0708       Vertex_t va;         // vertex i of lower base
0709       Vector3D<Real_v> pa; // same vertex converted to backend type
0710       // Point A is the current vertex on lower Z. P is the point we come from.
0711       pa.Set(Real_v(fxa[i]), Real_v(fya[i]), Real_v(-fDz));
0712       Vector3D<Real_v> vecAP = point - pa;
0713       // Dot product between AP vector and normal to surface has to be negative
0714       Real_v dotAPNorm = vecAP.Dot(fNormals[i]);
0715       // Dot product between direction and normal to surface has to be positive
0716       Real_v dotDirNorm = dir.Dot(fNormals[i]);
0717       smin[i]           = -dotAPNorm / NonZero(dotDirNorm);
0718       smax[i]           = InfinityLength<Real_v>(); // not to be checked
0719     }
0720   } // end ComputeSminSmax
0721 
0722 }; // end class definition
0723 
0724 } // namespace VECGEOM_IMPL_NAMESPACE
0725 
0726 } // namespace vecgeom
0727 
0728 #endif /* SECONDORDERSURFACESHELL_H_ */
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!