Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-12-05 02:48:09 

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