Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-06-12 23:30:07 

0001 #ifndef RecoTracker_LSTCore_src_alpaka_Segment_h
0002 #define RecoTracker_LSTCore_src_alpaka_Segment_h
0003 
0004 #include <limits>
0005 
0006 #include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
0007 #include "FWCore/Utilities/interface/CMSUnrollLoop.h"
0008 
0009 #include "RecoTracker/LSTCore/interface/alpaka/Common.h"
0010 #include "RecoTracker/LSTCore/interface/SegmentsSoA.h"
0011 #include "RecoTracker/LSTCore/interface/alpaka/SegmentsDeviceCollection.h"
0012 #include "RecoTracker/LSTCore/interface/PixelSegmentsSoA.h"
0013 #include "RecoTracker/LSTCore/interface/ModulesSoA.h"
0014 #include "RecoTracker/LSTCore/interface/HitsSoA.h"
0015 #include "RecoTracker/LSTCore/interface/MiniDoubletsSoA.h"
0016 #include "RecoTracker/LSTCore/interface/EndcapGeometry.h"
0017 #include "RecoTracker/LSTCore/interface/ObjectRangesSoA.h"
0018 
0019 #include "NeuralNetwork.h"
0020 
0021 namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
0022 
0023   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool isTighterTiltedModules_seg(ModulesConst modules, unsigned int moduleIndex) {
0024     // The "tighter" tilted modules are the subset of tilted modules that have smaller spacing
0025     // This is the same as what was previously considered as"isNormalTiltedModules"
0026     // See Figure 9.1 of https://cds.cern.ch/record/2272264/files/CMS-TDR-014.pdf
0027     short subdet = modules.subdets()[moduleIndex];
0028     short layer = modules.layers()[moduleIndex];
0029     short side = modules.sides()[moduleIndex];
0030     short rod = modules.rods()[moduleIndex];
0031 
0032     return (subdet == Barrel) && (((side != Center) && (layer == 3)) ||
0033                                   ((side == NegZ) && (((layer == 2) && (rod > 5)) || ((layer == 1) && (rod > 9)))) ||
0034                                   ((side == PosZ) && (((layer == 2) && (rod < 8)) || ((layer == 1) && (rod < 4)))));
0035   }
0036 
0037   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool isTighterTiltedModules_seg(short subdet, short layer, short side, short rod) {
0038     // The "tighter" tilted modules are the subset of tilted modules that have smaller spacing
0039     // This is the same as what was previously considered as"isNormalTiltedModules"
0040     // See Figure 9.1 of https://cds.cern.ch/record/2272264/files/CMS-TDR-014.pdf
0041     return (subdet == Barrel) && (((side != Center) && (layer == 3)) ||
0042                                   ((side == NegZ) && (((layer == 2) && (rod > 5)) || ((layer == 1) && (rod > 9)))) ||
0043                                   ((side == PosZ) && (((layer == 2) && (rod < 8)) || ((layer == 1) && (rod < 4)))));
0044   }
0045 
0046   ALPAKA_FN_ACC ALPAKA_FN_INLINE float moduleGapSize_seg(short layer, short ring, short subdet, short side, short rod) {
0047     unsigned int iL = layer - 1;
0048     unsigned int iR = ring - 1;
0049 
0050     float moduleSeparation = 0;
0051 
0052     if (subdet == Barrel and side == Center) {
0053       moduleSeparation = kMiniDeltaFlat[iL];
0054     } else if (isTighterTiltedModules_seg(subdet, layer, side, rod)) {
0055       moduleSeparation = kMiniDeltaTilted[iL];
0056     } else if (subdet == Endcap) {
0057       moduleSeparation = kMiniDeltaEndcap[iL][iR];
0058     } else  //Loose tilted modules
0059     {
0060       moduleSeparation = kMiniDeltaLooseTilted[iL];
0061     }
0062 
0063     return moduleSeparation;
0064   }
0065 
0066   ALPAKA_FN_ACC ALPAKA_FN_INLINE float moduleGapSize_seg(ModulesConst modules, unsigned int moduleIndex) {
0067     unsigned int iL = modules.layers()[moduleIndex] - 1;
0068     unsigned int iR = modules.rings()[moduleIndex] - 1;
0069     short subdet = modules.subdets()[moduleIndex];
0070     short side = modules.sides()[moduleIndex];
0071 
0072     float moduleSeparation = 0;
0073 
0074     if (subdet == Barrel and side == Center) {
0075       moduleSeparation = kMiniDeltaFlat[iL];
0076     } else if (isTighterTiltedModules_seg(modules, moduleIndex)) {
0077       moduleSeparation = kMiniDeltaTilted[iL];
0078     } else if (subdet == Endcap) {
0079       moduleSeparation = kMiniDeltaEndcap[iL][iR];
0080     } else  //Loose tilted modules
0081     {
0082       moduleSeparation = kMiniDeltaLooseTilted[iL];
0083     }
0084 
0085     return moduleSeparation;
0086   }
0087 
0088   template <typename TAcc>
0089   ALPAKA_FN_ACC ALPAKA_FN_INLINE void dAlphaThreshold(TAcc const& acc,
0090                                                       float* dAlphaThresholdValues,
0091                                                       ModulesConst modules,
0092                                                       MiniDoubletsConst mds,
0093                                                       float xIn,
0094                                                       float yIn,
0095                                                       float zIn,
0096                                                       float rtIn,
0097                                                       float xOut,
0098                                                       float yOut,
0099                                                       float zOut,
0100                                                       float rtOut,
0101                                                       uint16_t innerLowerModuleIndex,
0102                                                       uint16_t outerLowerModuleIndex,
0103                                                       unsigned int innerMDIndex,
0104                                                       unsigned int outerMDIndex,
0105                                                       const float ptCut) {
0106     float sdMuls = (modules.subdets()[innerLowerModuleIndex] == Barrel)
0107                        ? kMiniMulsPtScaleBarrel[modules.layers()[innerLowerModuleIndex] - 1] * 3.f / ptCut
0108                        : kMiniMulsPtScaleEndcap[modules.layers()[innerLowerModuleIndex] - 1] * 3.f / ptCut;
0109 
0110     //more accurate then outer rt - inner rt
0111     float segmentDr = alpaka::math::sqrt(acc, (yOut - yIn) * (yOut - yIn) + (xOut - xIn) * (xOut - xIn));
0112 
0113     const float dAlpha_Bfield =
0114         alpaka::math::asin(acc, alpaka::math::min(acc, segmentDr * k2Rinv1GeVf / ptCut, kSinAlphaMax));
0115 
0116     bool isInnerTilted =
0117         modules.subdets()[innerLowerModuleIndex] == Barrel and modules.sides()[innerLowerModuleIndex] != Center;
0118     bool isOuterTilted =
0119         modules.subdets()[outerLowerModuleIndex] == Barrel and modules.sides()[outerLowerModuleIndex] != Center;
0120 
0121     float drdzInner = modules.drdzs()[innerLowerModuleIndex];
0122     float drdzOuter = modules.drdzs()[outerLowerModuleIndex];
0123     float innerModuleGapSize = moduleGapSize_seg(modules, innerLowerModuleIndex);
0124     float outerModuleGapSize = moduleGapSize_seg(modules, outerLowerModuleIndex);
0125     const float innerminiTilt2 = isInnerTilted
0126                                      ? ((0.5f * 0.5f) * (kPixelPSZpitch * kPixelPSZpitch) * (drdzInner * drdzInner) /
0127                                         (1.f + drdzInner * drdzInner) / (innerModuleGapSize * innerModuleGapSize))
0128                                      : 0;
0129 
0130     const float outerminiTilt2 = isOuterTilted
0131                                      ? ((0.5f * 0.5f) * (kPixelPSZpitch * kPixelPSZpitch) * (drdzOuter * drdzOuter) /
0132                                         (1.f + drdzOuter * drdzOuter) / (outerModuleGapSize * outerModuleGapSize))
0133                                      : 0;
0134 
0135     float miniDelta = innerModuleGapSize;
0136 
0137     float sdLumForInnerMini2;
0138     float sdLumForOuterMini2;
0139 
0140     if (modules.subdets()[innerLowerModuleIndex] == Barrel) {
0141       sdLumForInnerMini2 = innerminiTilt2 * (dAlpha_Bfield * dAlpha_Bfield);
0142     } else {
0143       sdLumForInnerMini2 = (mds.dphis()[innerMDIndex] * mds.dphis()[innerMDIndex]) * (kDeltaZLum * kDeltaZLum) /
0144                            (mds.dzs()[innerMDIndex] * mds.dzs()[innerMDIndex]);
0145     }
0146 
0147     if (modules.subdets()[outerLowerModuleIndex] == Barrel) {
0148       sdLumForOuterMini2 = outerminiTilt2 * (dAlpha_Bfield * dAlpha_Bfield);
0149     } else {
0150       sdLumForOuterMini2 = (mds.dphis()[outerMDIndex] * mds.dphis()[outerMDIndex]) * (kDeltaZLum * kDeltaZLum) /
0151                            (mds.dzs()[outerMDIndex] * mds.dzs()[outerMDIndex]);
0152     }
0153 
0154     // Unique stuff for the segment dudes alone
0155     float dAlpha_res_inner =
0156         0.02f / miniDelta *
0157         (modules.subdets()[innerLowerModuleIndex] == Barrel ? 1.0f : alpaka::math::abs(acc, zIn) / rtIn);
0158     float dAlpha_res_outer =
0159         0.02f / miniDelta *
0160         (modules.subdets()[outerLowerModuleIndex] == Barrel ? 1.0f : alpaka::math::abs(acc, zOut) / rtOut);
0161 
0162     float dAlpha_res = dAlpha_res_inner + dAlpha_res_outer;
0163 
0164     if (modules.subdets()[innerLowerModuleIndex] == Barrel and modules.sides()[innerLowerModuleIndex] == Center) {
0165       dAlphaThresholdValues[0] = dAlpha_Bfield + alpaka::math::sqrt(acc, dAlpha_res * dAlpha_res + sdMuls * sdMuls);
0166     } else {
0167       dAlphaThresholdValues[0] =
0168           dAlpha_Bfield + alpaka::math::sqrt(acc, dAlpha_res * dAlpha_res + sdMuls * sdMuls + sdLumForInnerMini2);
0169     }
0170 
0171     if (modules.subdets()[outerLowerModuleIndex] == Barrel and modules.sides()[outerLowerModuleIndex] == Center) {
0172       dAlphaThresholdValues[1] = dAlpha_Bfield + alpaka::math::sqrt(acc, dAlpha_res * dAlpha_res + sdMuls * sdMuls);
0173     } else {
0174       dAlphaThresholdValues[1] =
0175           dAlpha_Bfield + alpaka::math::sqrt(acc, dAlpha_res * dAlpha_res + sdMuls * sdMuls + sdLumForOuterMini2);
0176     }
0177 
0178     //Inner to outer
0179     dAlphaThresholdValues[2] = dAlpha_Bfield + alpaka::math::sqrt(acc, dAlpha_res * dAlpha_res + sdMuls * sdMuls);
0180   }
0181 
0182   ALPAKA_FN_ACC ALPAKA_FN_INLINE void addSegmentToMemory(Segments segments,
0183                                                          unsigned int lowerMDIndex,
0184                                                          unsigned int upperMDIndex,
0185                                                          uint16_t innerLowerModuleIndex,
0186                                                          uint16_t outerLowerModuleIndex,
0187                                                          unsigned int innerMDAnchorHitIndex,
0188                                                          unsigned int outerMDAnchorHitIndex,
0189                                                          float dPhi,
0190                                                          float dPhiMin,
0191                                                          float dPhiMax,
0192                                                          float dPhiChange,
0193                                                          float dPhiChangeMin,
0194                                                          float dPhiChangeMax,
0195                                                          unsigned int idx) {
0196     segments.mdIndices()[idx][0] = lowerMDIndex;
0197     segments.mdIndices()[idx][1] = upperMDIndex;
0198     segments.innerLowerModuleIndices()[idx] = innerLowerModuleIndex;
0199     segments.outerLowerModuleIndices()[idx] = outerLowerModuleIndex;
0200     segments.innerMiniDoubletAnchorHitIndices()[idx] = innerMDAnchorHitIndex;
0201     segments.outerMiniDoubletAnchorHitIndices()[idx] = outerMDAnchorHitIndex;
0202 
0203     segments.dPhis()[idx] = __F2H(dPhi);
0204     segments.dPhiMins()[idx] = __F2H(dPhiMin);
0205     segments.dPhiMaxs()[idx] = __F2H(dPhiMax);
0206     segments.dPhiChanges()[idx] = __F2H(dPhiChange);
0207     segments.dPhiChangeMins()[idx] = __F2H(dPhiChangeMin);
0208     segments.dPhiChangeMaxs()[idx] = __F2H(dPhiChangeMax);
0209   }
0210 
0211   template <typename TAcc>
0212   ALPAKA_FN_ACC ALPAKA_FN_INLINE void addPixelSegmentToMemory(TAcc const& acc,
0213                                                               Segments segments,
0214                                                               PixelSegments pixelSegments,
0215                                                               PixelSeedsConst pixelSeeds,
0216                                                               MiniDoubletsConst mds,
0217                                                               unsigned int innerMDIndex,
0218                                                               unsigned int outerMDIndex,
0219                                                               uint16_t pixelModuleIndex,
0220                                                               const Params_pLS::ArrayUxHits& hitIdxs,
0221                                                               unsigned int innerAnchorHitIndex,
0222                                                               unsigned int outerAnchorHitIndex,
0223                                                               float dPhiChange,
0224                                                               unsigned int idx,
0225                                                               unsigned int pixelSegmentArrayIndex,
0226                                                               float score) {
0227     segments.mdIndices()[idx][0] = innerMDIndex;
0228     segments.mdIndices()[idx][1] = outerMDIndex;
0229     segments.innerLowerModuleIndices()[idx] = pixelModuleIndex;
0230     segments.outerLowerModuleIndices()[idx] = pixelModuleIndex;
0231     segments.innerMiniDoubletAnchorHitIndices()[idx] = innerAnchorHitIndex;
0232     segments.outerMiniDoubletAnchorHitIndices()[idx] = outerAnchorHitIndex;
0233     segments.dPhiChanges()[idx] = __F2H(dPhiChange);
0234 
0235     pixelSegments.isDup()[pixelSegmentArrayIndex] = false;
0236     pixelSegments.partOfPT5()[pixelSegmentArrayIndex] = false;
0237     pixelSegments.score()[pixelSegmentArrayIndex] = score;
0238     pixelSegments.pLSHitsIdxs()[pixelSegmentArrayIndex] = hitIdxs;
0239 
0240     //computing circle parameters
0241     /*
0242     The two anchor hits are r3PCA and r3LH. p3PCA pt, eta, phi is hitIndex1 x, y, z
0243     */
0244     float circleRadius = mds.outerX()[innerMDIndex] / (2 * k2Rinv1GeVf);
0245     float circlePhi = mds.outerZ()[innerMDIndex];
0246     float candidateCenterXs[] = {mds.anchorX()[innerMDIndex] + circleRadius * alpaka::math::sin(acc, circlePhi),
0247                                  mds.anchorX()[innerMDIndex] - circleRadius * alpaka::math::sin(acc, circlePhi)};
0248     float candidateCenterYs[] = {mds.anchorY()[innerMDIndex] - circleRadius * alpaka::math::cos(acc, circlePhi),
0249                                  mds.anchorY()[innerMDIndex] + circleRadius * alpaka::math::cos(acc, circlePhi)};
0250 
0251     //check which of the circles can accommodate r3LH better (we won't get perfect agreement)
0252     float bestChiSquared = std::numeric_limits<float>::infinity();
0253     float chiSquared;
0254     size_t bestIndex;
0255     for (size_t i = 0; i < 2; i++) {
0256       chiSquared = alpaka::math::abs(acc,
0257                                      alpaka::math::sqrt(acc,
0258                                                         (mds.anchorX()[outerMDIndex] - candidateCenterXs[i]) *
0259                                                                 (mds.anchorX()[outerMDIndex] - candidateCenterXs[i]) +
0260                                                             (mds.anchorY()[outerMDIndex] - candidateCenterYs[i]) *
0261                                                                 (mds.anchorY()[outerMDIndex] - candidateCenterYs[i])) -
0262                                          circleRadius);
0263       if (chiSquared < bestChiSquared) {
0264         bestChiSquared = chiSquared;
0265         bestIndex = i;
0266       }
0267     }
0268     pixelSegments.circleCenterX()[pixelSegmentArrayIndex] = candidateCenterXs[bestIndex];
0269     pixelSegments.circleCenterY()[pixelSegmentArrayIndex] = candidateCenterYs[bestIndex];
0270     pixelSegments.circleRadius()[pixelSegmentArrayIndex] = circleRadius;
0271 
0272     float plsEmbed[Params_pLS::kEmbed];
0273     plsembdnn::runEmbed(acc,
0274                         pixelSeeds.eta()[pixelSegmentArrayIndex],
0275                         pixelSeeds.etaErr()[pixelSegmentArrayIndex],
0276                         pixelSeeds.phi()[pixelSegmentArrayIndex],
0277                         pixelSegments.circleCenterX()[pixelSegmentArrayIndex],
0278                         pixelSegments.circleCenterY()[pixelSegmentArrayIndex],
0279                         pixelSegments.circleRadius()[pixelSegmentArrayIndex],
0280                         pixelSeeds.ptIn()[pixelSegmentArrayIndex],
0281                         pixelSeeds.ptErr()[pixelSegmentArrayIndex],
0282                         static_cast<bool>(pixelSeeds.isQuad()[pixelSegmentArrayIndex]),
0283                         plsEmbed);
0284 
0285     CMS_UNROLL_LOOP for (unsigned k = 0; k < Params_pLS::kEmbed; ++k) {
0286       pixelSegments.plsEmbed()[pixelSegmentArrayIndex][k] = plsEmbed[k];
0287     }
0288   }
0289 
0290   template <typename TAcc>
0291   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runSegmentDefaultAlgoBarrel(TAcc const& acc,
0292                                                                   ModulesConst modules,
0293                                                                   MiniDoubletsConst mds,
0294                                                                   uint16_t innerLowerModuleIndex,
0295                                                                   uint16_t outerLowerModuleIndex,
0296                                                                   unsigned int innerMDIndex,
0297                                                                   unsigned int outerMDIndex,
0298                                                                   float& dPhi,
0299                                                                   float& dPhiMin,
0300                                                                   float& dPhiMax,
0301                                                                   float& dPhiChange,
0302                                                                   float& dPhiChangeMin,
0303                                                                   float& dPhiChangeMax,
0304                                                                   const float ptCut) {
0305     float sdMuls = (modules.subdets()[innerLowerModuleIndex] == Barrel)
0306                        ? kMiniMulsPtScaleBarrel[modules.layers()[innerLowerModuleIndex] - 1] * 3.f / ptCut
0307                        : kMiniMulsPtScaleEndcap[modules.layers()[innerLowerModuleIndex] - 1] * 3.f / ptCut;
0308 
0309     float xIn, yIn, zIn, rtIn, xOut, yOut, zOut, rtOut;
0310 
0311     xIn = mds.anchorX()[innerMDIndex];
0312     yIn = mds.anchorY()[innerMDIndex];
0313     zIn = mds.anchorZ()[innerMDIndex];
0314     rtIn = mds.anchorRt()[innerMDIndex];
0315 
0316     xOut = mds.anchorX()[outerMDIndex];
0317     yOut = mds.anchorY()[outerMDIndex];
0318     zOut = mds.anchorZ()[outerMDIndex];
0319     rtOut = mds.anchorRt()[outerMDIndex];
0320 
0321     float sdSlope = alpaka::math::asin(acc, alpaka::math::min(acc, rtOut * k2Rinv1GeVf / ptCut, kSinAlphaMax));
0322     float sdPVoff = 0.1f / rtOut;
0323     float dzDrtScale = alpaka::math::tan(acc, sdSlope) / sdSlope;  //FIXME: need appropriate value
0324 
0325     const float zGeom = modules.layers()[innerLowerModuleIndex] <= 2 ? 2.f * kPixelPSZpitch : 2.f * kStrip2SZpitch;
0326 
0327     float zLo = zIn + (zIn - kDeltaZLum) * (rtOut / rtIn - 1.f) * (zIn > 0.f ? 1.f : dzDrtScale) -
0328                 zGeom;  //slope-correction only on outer end
0329     float zHi = zIn + (zIn + kDeltaZLum) * (rtOut / rtIn - 1.f) * (zIn < 0.f ? 1.f : dzDrtScale) + zGeom;
0330 
0331     if ((zOut < zLo) || (zOut > zHi))
0332       return false;
0333 
0334     float sdCut = sdSlope + alpaka::math::sqrt(acc, sdMuls * sdMuls + sdPVoff * sdPVoff);
0335 
0336     dPhi = cms::alpakatools::reducePhiRange(acc, mds.anchorPhi()[outerMDIndex] - mds.anchorPhi()[innerMDIndex]);
0337 
0338     if (alpaka::math::abs(acc, dPhi) > sdCut)
0339       return false;
0340 
0341     dPhiChange = cms::alpakatools::reducePhiRange(
0342         acc, cms::alpakatools::phi(acc, xOut - xIn, yOut - yIn) - mds.anchorPhi()[innerMDIndex]);
0343 
0344     if (alpaka::math::abs(acc, dPhiChange) > sdCut)
0345       return false;
0346 
0347     float dAlphaThresholdValues[3];
0348     dAlphaThreshold(acc,
0349                     dAlphaThresholdValues,
0350                     modules,
0351                     mds,
0352                     xIn,
0353                     yIn,
0354                     zIn,
0355                     rtIn,
0356                     xOut,
0357                     yOut,
0358                     zOut,
0359                     rtOut,
0360                     innerLowerModuleIndex,
0361                     outerLowerModuleIndex,
0362                     innerMDIndex,
0363                     outerMDIndex,
0364                     ptCut);
0365 
0366     float innerMDAlpha = mds.dphichanges()[innerMDIndex];
0367     float outerMDAlpha = mds.dphichanges()[outerMDIndex];
0368     float dAlphaInnerMDSegment = innerMDAlpha - dPhiChange;
0369     float dAlphaOuterMDSegment = outerMDAlpha - dPhiChange;
0370     float dAlphaInnerMDOuterMD = innerMDAlpha - outerMDAlpha;
0371 
0372     float dAlphaInnerMDSegmentThreshold = dAlphaThresholdValues[0];
0373     float dAlphaOuterMDSegmentThreshold = dAlphaThresholdValues[1];
0374     float dAlphaInnerMDOuterMDThreshold = dAlphaThresholdValues[2];
0375 
0376     if (alpaka::math::abs(acc, dAlphaInnerMDSegment) >= dAlphaInnerMDSegmentThreshold)
0377       return false;
0378     if (alpaka::math::abs(acc, dAlphaOuterMDSegment) >= dAlphaOuterMDSegmentThreshold)
0379       return false;
0380     return alpaka::math::abs(acc, dAlphaInnerMDOuterMD) < dAlphaInnerMDOuterMDThreshold;
0381   }
0382 
0383   template <typename TAcc>
0384   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runSegmentDefaultAlgoEndcap(TAcc const& acc,
0385                                                                   ModulesConst modules,
0386                                                                   MiniDoubletsConst mds,
0387                                                                   uint16_t innerLowerModuleIndex,
0388                                                                   uint16_t outerLowerModuleIndex,
0389                                                                   unsigned int innerMDIndex,
0390                                                                   unsigned int outerMDIndex,
0391                                                                   float& dPhi,
0392                                                                   float& dPhiMin,
0393                                                                   float& dPhiMax,
0394                                                                   float& dPhiChange,
0395                                                                   float& dPhiChangeMin,
0396                                                                   float& dPhiChangeMax,
0397                                                                   const float ptCut) {
0398     float xIn, yIn, zIn, rtIn, xOut, yOut, zOut, rtOut;
0399 
0400     xIn = mds.anchorX()[innerMDIndex];
0401     yIn = mds.anchorY()[innerMDIndex];
0402     zIn = mds.anchorZ()[innerMDIndex];
0403     rtIn = mds.anchorRt()[innerMDIndex];
0404 
0405     xOut = mds.anchorX()[outerMDIndex];
0406     yOut = mds.anchorY()[outerMDIndex];
0407     zOut = mds.anchorZ()[outerMDIndex];
0408     rtOut = mds.anchorRt()[outerMDIndex];
0409 
0410     bool outerLayerEndcapTwoS =
0411         (modules.subdets()[outerLowerModuleIndex] == Endcap) && (modules.moduleType()[outerLowerModuleIndex] == TwoS);
0412 
0413     float sdSlope = alpaka::math::asin(acc, alpaka::math::min(acc, rtOut * k2Rinv1GeVf / ptCut, kSinAlphaMax));
0414     float disks2SMinRadius = 60.f;
0415 
0416     float rtGeom = ((rtIn < disks2SMinRadius && rtOut < disks2SMinRadius)
0417                         ? (2.f * kPixelPSZpitch)
0418                         : ((rtIn < disks2SMinRadius || rtOut < disks2SMinRadius) ? (kPixelPSZpitch + kStrip2SZpitch)
0419                                                                                  : (2.f * kStrip2SZpitch)));
0420 
0421     //cut 0 - z compatibility
0422     if (zIn * zOut < 0)
0423       return false;
0424 
0425     float dz = zOut - zIn;
0426     float dLum = alpaka::math::copysign(acc, kDeltaZLum, zIn);
0427     float drtDzScale = sdSlope / alpaka::math::tan(acc, sdSlope);
0428 
0429     float rtLo = alpaka::math::max(
0430         acc, rtIn * (1.f + dz / (zIn + dLum) * drtDzScale) - rtGeom, rtIn - 0.5f * rtGeom);  //rt should increase
0431     float rtHi = rtIn * (zOut - dLum) / (zIn - dLum) +
0432                  rtGeom;  //dLum for luminous; rGeom for measurement size; no tanTheta_loc(pt) correction
0433 
0434     // Completeness
0435     if ((rtOut < rtLo) || (rtOut > rtHi))
0436       return false;
0437 
0438     dPhi = cms::alpakatools::reducePhiRange(acc, mds.anchorPhi()[outerMDIndex] - mds.anchorPhi()[innerMDIndex]);
0439 
0440     float sdCut = sdSlope;
0441     if (outerLayerEndcapTwoS) {
0442       float dPhiPos_high =
0443           cms::alpakatools::reducePhiRange(acc, mds.anchorHighEdgePhi()[outerMDIndex] - mds.anchorPhi()[innerMDIndex]);
0444       float dPhiPos_low =
0445           cms::alpakatools::reducePhiRange(acc, mds.anchorLowEdgePhi()[outerMDIndex] - mds.anchorPhi()[innerMDIndex]);
0446 
0447       dPhiMax = alpaka::math::abs(acc, dPhiPos_high) > alpaka::math::abs(acc, dPhiPos_low) ? dPhiPos_high : dPhiPos_low;
0448       dPhiMin = alpaka::math::abs(acc, dPhiPos_high) > alpaka::math::abs(acc, dPhiPos_low) ? dPhiPos_low : dPhiPos_high;
0449     } else {
0450       dPhiMax = dPhi;
0451       dPhiMin = dPhi;
0452     }
0453     if (alpaka::math::abs(acc, dPhi) > sdCut)
0454       return false;
0455 
0456     float dzFrac = dz / zIn;
0457     dPhiChange = dPhi / dzFrac * (1.f + dzFrac);
0458     dPhiChangeMin = dPhiMin / dzFrac * (1.f + dzFrac);
0459     dPhiChangeMax = dPhiMax / dzFrac * (1.f + dzFrac);
0460 
0461     if (alpaka::math::abs(acc, dPhiChange) > sdCut)
0462       return false;
0463 
0464     float dAlphaThresholdValues[3];
0465     dAlphaThreshold(acc,
0466                     dAlphaThresholdValues,
0467                     modules,
0468                     mds,
0469                     xIn,
0470                     yIn,
0471                     zIn,
0472                     rtIn,
0473                     xOut,
0474                     yOut,
0475                     zOut,
0476                     rtOut,
0477                     innerLowerModuleIndex,
0478                     outerLowerModuleIndex,
0479                     innerMDIndex,
0480                     outerMDIndex,
0481                     ptCut);
0482 
0483     float innerMDAlpha = mds.dphichanges()[innerMDIndex];
0484     float outerMDAlpha = mds.dphichanges()[outerMDIndex];
0485     float dAlphaInnerMDSegment = innerMDAlpha - dPhiChange;
0486     float dAlphaOuterMDSegment = outerMDAlpha - dPhiChange;
0487     float dAlphaInnerMDOuterMD = innerMDAlpha - outerMDAlpha;
0488 
0489     float dAlphaInnerMDSegmentThreshold = dAlphaThresholdValues[0];
0490     float dAlphaOuterMDSegmentThreshold = dAlphaThresholdValues[1];
0491     float dAlphaInnerMDOuterMDThreshold = dAlphaThresholdValues[2];
0492 
0493     if (alpaka::math::abs(acc, dAlphaInnerMDSegment) >= dAlphaInnerMDSegmentThreshold)
0494       return false;
0495     if (alpaka::math::abs(acc, dAlphaOuterMDSegment) >= dAlphaOuterMDSegmentThreshold)
0496       return false;
0497     return alpaka::math::abs(acc, dAlphaInnerMDOuterMD) < dAlphaInnerMDOuterMDThreshold;
0498   }
0499 
0500   template <typename TAcc>
0501   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runSegmentDefaultAlgo(TAcc const& acc,
0502                                                             ModulesConst modules,
0503                                                             MiniDoubletsConst mds,
0504                                                             uint16_t innerLowerModuleIndex,
0505                                                             uint16_t outerLowerModuleIndex,
0506                                                             unsigned int innerMDIndex,
0507                                                             unsigned int outerMDIndex,
0508                                                             float& dPhi,
0509                                                             float& dPhiMin,
0510                                                             float& dPhiMax,
0511                                                             float& dPhiChange,
0512                                                             float& dPhiChangeMin,
0513                                                             float& dPhiChangeMax,
0514                                                             const float ptCut) {
0515     if (modules.subdets()[innerLowerModuleIndex] == Barrel and modules.subdets()[outerLowerModuleIndex] == Barrel) {
0516       return runSegmentDefaultAlgoBarrel(acc,
0517                                          modules,
0518                                          mds,
0519                                          innerLowerModuleIndex,
0520                                          outerLowerModuleIndex,
0521                                          innerMDIndex,
0522                                          outerMDIndex,
0523                                          dPhi,
0524                                          dPhiMin,
0525                                          dPhiMax,
0526                                          dPhiChange,
0527                                          dPhiChangeMin,
0528                                          dPhiChangeMax,
0529                                          ptCut);
0530     } else {
0531       return runSegmentDefaultAlgoEndcap(acc,
0532                                          modules,
0533                                          mds,
0534                                          innerLowerModuleIndex,
0535                                          outerLowerModuleIndex,
0536                                          innerMDIndex,
0537                                          outerMDIndex,
0538                                          dPhi,
0539                                          dPhiMin,
0540                                          dPhiMax,
0541                                          dPhiChange,
0542                                          dPhiChangeMin,
0543                                          dPhiChangeMax,
0544                                          ptCut);
0545     }
0546   }
0547 
0548   struct CreateSegments {
0549     ALPAKA_FN_ACC void operator()(Acc3D const& acc,
0550                                   ModulesConst modules,
0551                                   MiniDoubletsConst mds,
0552                                   MiniDoubletsOccupancyConst mdsOccupancy,
0553                                   Segments segments,
0554                                   SegmentsOccupancy segmentsOccupancy,
0555                                   ObjectRangesConst ranges,
0556                                   const float ptCut) const {
0557       for (uint16_t innerLowerModuleIndex : cms::alpakatools::uniform_elements_z(acc, modules.nLowerModules())) {
0558         unsigned int nInnerMDs = mdsOccupancy.nMDs()[innerLowerModuleIndex];
0559         if (nInnerMDs == 0)
0560           continue;
0561 
0562         unsigned int nConnectedModules = modules.nConnectedModules()[innerLowerModuleIndex];
0563 
0564         for (uint16_t outerLowerModuleArrayIdx : cms::alpakatools::uniform_elements_y(acc, nConnectedModules)) {
0565           uint16_t outerLowerModuleIndex = modules.moduleMap()[innerLowerModuleIndex][outerLowerModuleArrayIdx];
0566 
0567           unsigned int nOuterMDs = mdsOccupancy.nMDs()[outerLowerModuleIndex];
0568 
0569           unsigned int limit = nInnerMDs * nOuterMDs;
0570 
0571           if (limit == 0)
0572             continue;
0573           for (unsigned int hitIndex : cms::alpakatools::uniform_elements_x(acc, limit)) {
0574             unsigned int innerMDArrayIdx = hitIndex / nOuterMDs;
0575             unsigned int outerMDArrayIdx = hitIndex % nOuterMDs;
0576             if (outerMDArrayIdx >= nOuterMDs)
0577               continue;
0578 
0579             unsigned int innerMDIndex = ranges.mdRanges()[innerLowerModuleIndex][0] + innerMDArrayIdx;
0580             unsigned int outerMDIndex = ranges.mdRanges()[outerLowerModuleIndex][0] + outerMDArrayIdx;
0581 
0582             float dPhi, dPhiMin, dPhiMax, dPhiChange, dPhiChangeMin, dPhiChangeMax;
0583 
0584             unsigned int innerMiniDoubletAnchorHitIndex = mds.anchorHitIndices()[innerMDIndex];
0585             unsigned int outerMiniDoubletAnchorHitIndex = mds.anchorHitIndices()[outerMDIndex];
0586             dPhiMin = 0;
0587             dPhiMax = 0;
0588             dPhiChangeMin = 0;
0589             dPhiChangeMax = 0;
0590             if (runSegmentDefaultAlgo(acc,
0591                                       modules,
0592                                       mds,
0593                                       innerLowerModuleIndex,
0594                                       outerLowerModuleIndex,
0595                                       innerMDIndex,
0596                                       outerMDIndex,
0597                                       dPhi,
0598                                       dPhiMin,
0599                                       dPhiMax,
0600                                       dPhiChange,
0601                                       dPhiChangeMin,
0602                                       dPhiChangeMax,
0603                                       ptCut)) {
0604               unsigned int totOccupancySegments =
0605                   alpaka::atomicAdd(acc,
0606                                     &segmentsOccupancy.totOccupancySegments()[innerLowerModuleIndex],
0607                                     1u,
0608                                     alpaka::hierarchy::Threads{});
0609               if (static_cast<int>(totOccupancySegments) >= ranges.segmentModuleOccupancy()[innerLowerModuleIndex]) {
0610 #ifdef WARNINGS
0611                 printf("Segment excess alert! Module index = %d, Occupancy = %d\n",
0612                        innerLowerModuleIndex,
0613                        totOccupancySegments);
0614 #endif
0615               } else {
0616                 unsigned int segmentModuleIdx = alpaka::atomicAdd(
0617                     acc, &segmentsOccupancy.nSegments()[innerLowerModuleIndex], 1u, alpaka::hierarchy::Threads{});
0618                 unsigned int segmentIdx = ranges.segmentModuleIndices()[innerLowerModuleIndex] + segmentModuleIdx;
0619 
0620                 addSegmentToMemory(segments,
0621                                    innerMDIndex,
0622                                    outerMDIndex,
0623                                    innerLowerModuleIndex,
0624                                    outerLowerModuleIndex,
0625                                    innerMiniDoubletAnchorHitIndex,
0626                                    outerMiniDoubletAnchorHitIndex,
0627                                    dPhi,
0628                                    dPhiMin,
0629                                    dPhiMax,
0630                                    dPhiChange,
0631                                    dPhiChangeMin,
0632                                    dPhiChangeMax,
0633                                    segmentIdx);
0634               }
0635             }
0636           }
0637         }
0638       }
0639     }
0640   };
0641 
0642   struct CreateSegmentArrayRanges {
0643     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
0644                                   ModulesConst modules,
0645                                   ObjectRanges ranges,
0646                                   MiniDoubletsConst mds,
0647                                   MiniDoubletsOccupancyConst mdsOccupancy,
0648                                   const float ptCut) const {
0649       // implementation is 1D with a single block
0650       ALPAKA_ASSERT_ACC((alpaka::getWorkDiv<alpaka::Grid, alpaka::Blocks>(acc)[0] == 1));
0651 
0652       // Initialize variables in shared memory and set to 0
0653       int& nTotalSegments = alpaka::declareSharedVar<int, __COUNTER__>(acc);
0654       if (cms::alpakatools::once_per_block(acc)) {
0655         nTotalSegments = 0;
0656       }
0657       alpaka::syncBlockThreads(acc);
0658 
0659       // Occupancy matrix for 0.8 GeV pT Cut
0660       constexpr int p08_occupancy_matrix[4][4] = {
0661           {572, 300, 183, 62},  // category 0
0662           {191, 128, 0, 0},     // category 1
0663           {0, 107, 102, 0},     // category 2
0664           {0, 64, 79, 85}       // category 3
0665       };
0666 
0667       // Occupancy matrix for 0.6 GeV pT Cut, 99.9%
0668       constexpr int p06_occupancy_matrix[4][4] = {
0669           {936, 351, 256, 61},  // category 0
0670           {1358, 763, 0, 0},    // category 1
0671           {0, 210, 268, 0},     // category 2
0672           {0, 60, 97, 96}       // category 3
0673       };
0674 
0675       // Select the appropriate occupancy matrix based on ptCut
0676       const auto& occupancy_matrix = (ptCut < 0.8f) ? p06_occupancy_matrix : p08_occupancy_matrix;
0677 
0678       for (uint16_t i : cms::alpakatools::uniform_elements(acc, modules.nLowerModules())) {
0679         if (modules.nConnectedModules()[i] == 0) {
0680           ranges.segmentModuleIndices()[i] = nTotalSegments;
0681           ranges.segmentModuleOccupancy()[i] = 0;
0682           continue;
0683         }
0684 
0685         short module_rings = modules.rings()[i];
0686         short module_layers = modules.layers()[i];
0687         short module_subdets = modules.subdets()[i];
0688         float module_eta = alpaka::math::abs(acc, modules.eta()[i]);
0689 
0690         int category_number = getCategoryNumber(module_layers, module_subdets, module_rings);
0691         int eta_number = getEtaBin(module_eta);
0692 
0693         int dynamic_count = 0;
0694         // Calculate dynamic limit based on connected modules
0695         unsigned int nInnerMDs = mdsOccupancy.nMDs()[i];
0696         for (unsigned int c = 0; c < modules.nConnectedModules()[i]; ++c) {
0697           uint16_t connectedModule = modules.moduleMap()[i][c];
0698           dynamic_count += nInnerMDs * mdsOccupancy.nMDs()[connectedModule];
0699         }
0700 
0701 #ifdef WARNINGS
0702         if (category_number == -1 || eta_number == -1) {
0703           printf("Unhandled case in createSegmentArrayRanges! Module index = %i\n", i);
0704         }
0705 #endif
0706         // Get matrix-based cap
0707         int matrix_cap =
0708             (category_number != -1 && eta_number != -1) ? occupancy_matrix[category_number][eta_number] : 0;
0709 
0710         // Cap occupancy at minimum of dynamic count and matrix value
0711         int occupancy = alpaka::math::min(acc, dynamic_count, matrix_cap);
0712 
0713         int nTotSegs = alpaka::atomicAdd(acc, &nTotalSegments, occupancy, alpaka::hierarchy::Threads{});
0714         ranges.segmentModuleIndices()[i] = nTotSegs;
0715         ranges.segmentModuleOccupancy()[i] = occupancy;
0716       }
0717 
0718       // Wait for all threads to finish before reporting final values
0719       alpaka::syncBlockThreads(acc);
0720       if (cms::alpakatools::once_per_block(acc)) {
0721         ranges.segmentModuleIndices()[modules.nLowerModules()] = nTotalSegments;
0722         ranges.nTotalSegs() = nTotalSegments;
0723       }
0724     }
0725   };
0726 
0727   struct AddSegmentRangesToEventExplicit {
0728     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
0729                                   ModulesConst modules,
0730                                   SegmentsOccupancyConst segmentsOccupancy,
0731                                   ObjectRanges ranges) const {
0732       // implementation is 1D with a single block
0733       ALPAKA_ASSERT_ACC((alpaka::getWorkDiv<alpaka::Grid, alpaka::Blocks>(acc)[0] == 1));
0734 
0735       for (uint16_t i : cms::alpakatools::uniform_elements(acc, modules.nLowerModules())) {
0736         if (segmentsOccupancy.nSegments()[i] == 0) {
0737           ranges.segmentRanges()[i][0] = -1;
0738           ranges.segmentRanges()[i][1] = -1;
0739         } else {
0740           ranges.segmentRanges()[i][0] = ranges.segmentModuleIndices()[i];
0741           ranges.segmentRanges()[i][1] = ranges.segmentModuleIndices()[i] + segmentsOccupancy.nSegments()[i] - 1;
0742         }
0743       }
0744     }
0745   };
0746 
0747   struct AddPixelSegmentToEventKernel {
0748     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
0749                                   ModulesConst modules,
0750                                   ObjectRangesConst ranges,
0751                                   HitsBaseConst hitsBase,
0752                                   HitsExtendedConst hitsExtended,
0753                                   PixelSeedsConst pixelSeeds,
0754                                   MiniDoublets mds,
0755                                   Segments segments,
0756                                   PixelSegments pixelSegments,
0757                                   uint16_t pixelModuleIndex,
0758                                   int size) const {
0759       for (int tid : cms::alpakatools::uniform_elements(acc, size)) {
0760         unsigned int innerMDIndex = ranges.miniDoubletModuleIndices()[pixelModuleIndex] + 2 * (tid);
0761         unsigned int outerMDIndex = ranges.miniDoubletModuleIndices()[pixelModuleIndex] + 2 * (tid) + 1;
0762         unsigned int pixelSegmentIndex = ranges.segmentModuleIndices()[pixelModuleIndex] + tid;
0763 
0764         addMDToMemory(acc,
0765                       mds,
0766                       hitsBase,
0767                       hitsExtended,
0768                       modules,
0769                       pixelSeeds.hitIndices()[tid][0],
0770                       pixelSeeds.hitIndices()[tid][1],
0771                       pixelModuleIndex,
0772                       0,
0773                       0,
0774                       0,
0775                       0,
0776                       0,
0777                       0,
0778                       0,
0779                       0,
0780                       innerMDIndex);
0781         addMDToMemory(acc,
0782                       mds,
0783                       hitsBase,
0784                       hitsExtended,
0785                       modules,
0786                       pixelSeeds.hitIndices()[tid][2],
0787                       pixelSeeds.hitIndices()[tid][3],
0788                       pixelModuleIndex,
0789                       0,
0790                       0,
0791                       0,
0792                       0,
0793                       0,
0794                       0,
0795                       0,
0796                       0,
0797                       outerMDIndex);
0798 
0799         //in outer hits - pt, eta, phi
0800         float slope = alpaka::math::sinh(acc, hitsBase.ys()[mds.outerHitIndices()[innerMDIndex]]);
0801         float intercept = hitsBase.zs()[mds.anchorHitIndices()[innerMDIndex]] -
0802                           slope * hitsExtended.rts()[mds.anchorHitIndices()[innerMDIndex]];
0803         float score_lsq = (hitsExtended.rts()[mds.anchorHitIndices()[outerMDIndex]] * slope + intercept) -
0804                           (hitsBase.zs()[mds.anchorHitIndices()[outerMDIndex]]);
0805         score_lsq = score_lsq * score_lsq;
0806 
0807         const Params_pLS::ArrayUxHits hits1{{hitsBase.idxs()[mds.anchorHitIndices()[innerMDIndex]],
0808                                              hitsBase.idxs()[mds.anchorHitIndices()[outerMDIndex]],
0809                                              hitsBase.idxs()[mds.outerHitIndices()[innerMDIndex]],
0810                                              hitsBase.idxs()[mds.outerHitIndices()[outerMDIndex]]}};
0811         addPixelSegmentToMemory(acc,
0812                                 segments,
0813                                 pixelSegments,
0814                                 pixelSeeds,
0815                                 mds,
0816                                 innerMDIndex,
0817                                 outerMDIndex,
0818                                 pixelModuleIndex,
0819                                 hits1,
0820                                 pixelSeeds.hitIndices()[tid][0],
0821                                 pixelSeeds.hitIndices()[tid][2],
0822                                 pixelSeeds.deltaPhi()[tid],
0823                                 pixelSegmentIndex,
0824                                 tid,
0825                                 score_lsq);
0826       }
0827     }
0828   };
0829 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE::lst
0830 
0831 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!