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
0025
0026
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
0039
0040
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
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
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
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
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
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
0241
0242
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
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;
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;
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
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);
0431 float rtHi = rtIn * (zOut - dLum) / (zIn - dLum) +
0432 rtGeom;
0433
0434
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
0650 ALPAKA_ASSERT_ACC((alpaka::getWorkDiv<alpaka::Grid, alpaka::Blocks>(acc)[0] == 1));
0651
0652
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
0660 constexpr int p08_occupancy_matrix[4][4] = {
0661 {572, 300, 183, 62},
0662 {191, 128, 0, 0},
0663 {0, 107, 102, 0},
0664 {0, 64, 79, 85}
0665 };
0666
0667
0668 constexpr int p06_occupancy_matrix[4][4] = {
0669 {936, 351, 256, 61},
0670 {1358, 763, 0, 0},
0671 {0, 210, 268, 0},
0672 {0, 60, 97, 96}
0673 };
0674
0675
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
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
0707 int matrix_cap =
0708 (category_number != -1 && eta_number != -1) ? occupancy_matrix[category_number][eta_number] : 0;
0709
0710
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
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
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
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 }
0830
0831 #endif