src/alpaka/MiniDoublet.h

0001 #ifndef RecoTracker_LSTCore_src_alpaka_MiniDoublet_h
0002 #define RecoTracker_LSTCore_src_alpaka_MiniDoublet_h
0003
0004 #include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
0005 #include "FWCore/Utilities/interface/isFinite.h"
0006
0007 #include "RecoTracker/LSTCore/interface/alpaka/Common.h"
0008 #include "RecoTracker/LSTCore/interface/HitsSoA.h"
0009 #include "RecoTracker/LSTCore/interface/MiniDoubletsSoA.h"
0010 #include "RecoTracker/LSTCore/interface/alpaka/MiniDoubletsDeviceCollection.h"
0011 #include "RecoTracker/LSTCore/interface/ModulesSoA.h"
0012 #include "RecoTracker/LSTCore/interface/EndcapGeometry.h"
0013 #include "RecoTracker/LSTCore/interface/ObjectRangesSoA.h"
0014
0015 namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
0016   template <typename TAcc>
0017   ALPAKA_FN_ACC ALPAKA_FN_INLINE void addMDToMemory(TAcc const& acc,
0018                                                     MiniDoublets mds,
0019                                                     HitsBaseConst hitsBase,
0020                                                     HitsExtendedConst hitsExtended,
0021                                                     ModulesConst modules,
0022                                                     unsigned int lowerHitIdx,
0023                                                     unsigned int upperHitIdx,
0024                                                     uint16_t lowerModuleIdx,
0025                                                     float dz,
0026                                                     float dPhi,
0027                                                     float dPhiChange,
0028                                                     float shiftedX,
0029                                                     float shiftedY,
0030                                                     float shiftedZ,
0031                                                     float noShiftedDphi,
0032                                                     float noShiftedDPhiChange,
0033                                                     unsigned int idx) {
0034     //the index into which this MD needs to be written will be computed in the kernel
0035     //nMDs variable will be incremented in the kernel, no need to worry about that here
0036
0037     mds.moduleIndices()[idx] = lowerModuleIdx;
0038     unsigned int anchorHitIndex, outerHitIndex;
0039     if (modules.moduleType()[lowerModuleIdx] == PS and modules.moduleLayerType()[lowerModuleIdx] == Strip) {
0040       mds.anchorHitIndices()[idx] = upperHitIdx;
0041       mds.outerHitIndices()[idx] = lowerHitIdx;
0042
0043       anchorHitIndex = upperHitIdx;
0044       outerHitIndex = lowerHitIdx;
0045     } else {
0046       mds.anchorHitIndices()[idx] = lowerHitIdx;
0047       mds.outerHitIndices()[idx] = upperHitIdx;
0048
0049       anchorHitIndex = lowerHitIdx;
0050       outerHitIndex = upperHitIdx;
0051     }
0052
0053     mds.dphichanges()[idx] = dPhiChange;
0054
0055     mds.dphis()[idx] = dPhi;
0056     mds.dzs()[idx] = dz;
0057     mds.shiftedXs()[idx] = shiftedX;
0058     mds.shiftedYs()[idx] = shiftedY;
0059     mds.shiftedZs()[idx] = shiftedZ;
0060
0061     mds.noShiftedDphis()[idx] = noShiftedDphi;
0062     mds.noShiftedDphiChanges()[idx] = noShiftedDPhiChange;
0063
0064     mds.anchorX()[idx] = hitsBase.xs()[anchorHitIndex];
0065     mds.anchorY()[idx] = hitsBase.ys()[anchorHitIndex];
0066     mds.anchorZ()[idx] = hitsBase.zs()[anchorHitIndex];
0067     mds.anchorRt()[idx] = hitsExtended.rts()[anchorHitIndex];
0068     mds.anchorPhi()[idx] = hitsExtended.phis()[anchorHitIndex];
0069     mds.anchorEta()[idx] = hitsExtended.etas()[anchorHitIndex];
0070     mds.anchorHighEdgeX()[idx] = hitsExtended.highEdgeXs()[anchorHitIndex];
0071     mds.anchorHighEdgeY()[idx] = hitsExtended.highEdgeYs()[anchorHitIndex];
0072     mds.anchorLowEdgeX()[idx] = hitsExtended.lowEdgeXs()[anchorHitIndex];
0073     mds.anchorLowEdgeY()[idx] = hitsExtended.lowEdgeYs()[anchorHitIndex];
0074     mds.anchorHighEdgePhi()[idx] = alpaka::math::atan2(acc, mds.anchorHighEdgeY()[idx], mds.anchorHighEdgeX()[idx]);
0075     mds.anchorLowEdgePhi()[idx] = alpaka::math::atan2(acc, mds.anchorLowEdgeY()[idx], mds.anchorLowEdgeX()[idx]);
0076
0077     mds.outerX()[idx] = hitsBase.xs()[outerHitIndex];
0078     mds.outerY()[idx] = hitsBase.ys()[outerHitIndex];
0079     mds.outerZ()[idx] = hitsBase.zs()[outerHitIndex];
0080     mds.outerRt()[idx] = hitsExtended.rts()[outerHitIndex];
0081     mds.outerPhi()[idx] = hitsExtended.phis()[outerHitIndex];
0082     mds.outerEta()[idx] = hitsExtended.etas()[outerHitIndex];
0083     mds.outerHighEdgeX()[idx] = hitsExtended.highEdgeXs()[outerHitIndex];
0084     mds.outerHighEdgeY()[idx] = hitsExtended.highEdgeYs()[outerHitIndex];
0085     mds.outerLowEdgeX()[idx] = hitsExtended.lowEdgeXs()[outerHitIndex];
0086     mds.outerLowEdgeY()[idx] = hitsExtended.lowEdgeYs()[outerHitIndex];
0087   }
0088
0089   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool isTighterTiltedModules(ModulesConst modules, uint16_t moduleIndex) {
0090     // The "tighter" tilted modules are the subset of tilted modules that have smaller spacing
0091     // This is the same as what was previously considered as"isNormalTiltedModules"
0092     // See Figure 9.1 of https://cds.cern.ch/record/2272264/files/CMS-TDR-014.pdf
0093     short subdet = modules.subdets()[moduleIndex];
0094     short layer = modules.layers()[moduleIndex];
0095     short side = modules.sides()[moduleIndex];
0096     short rod = modules.rods()[moduleIndex];
0097
0098     if (subdet == Barrel) {
0099       if ((side != Center and layer == 3) or (side == NegZ and layer == 2 and rod > 5) or
0100           (side == PosZ and layer == 2 and rod < 8) or (side == NegZ and layer == 1 and rod > 9) or
0101           (side == PosZ and layer == 1 and rod < 4))
0102         return true;
0103       else
0104         return false;
0105     } else
0106       return false;
0107   }
0108
0109   ALPAKA_FN_ACC ALPAKA_FN_INLINE float moduleGapSize(ModulesConst modules, uint16_t moduleIndex) {
0110     unsigned int iL = modules.layers()[moduleIndex] - 1;
0111     unsigned int iR = modules.rings()[moduleIndex] - 1;
0112     short subdet = modules.subdets()[moduleIndex];
0113     short side = modules.sides()[moduleIndex];
0114
0115     float moduleSeparation = 0;
0116
0117     if (subdet == Barrel and side == Center) {
0118       moduleSeparation = kMiniDeltaFlat[iL];
0119     } else if (isTighterTiltedModules(modules, moduleIndex)) {
0120       moduleSeparation = kMiniDeltaTilted[iL];
0121     } else if (subdet == Endcap) {
0122       moduleSeparation = kMiniDeltaEndcap[iL][iR];
0123     } else  //Loose tilted modules
0124     {
0125       moduleSeparation = kMiniDeltaLooseTilted[iL];
0126     }
0127
0128     return moduleSeparation;
0129   }
0130
0131   template <typename TAcc>
0132   ALPAKA_FN_ACC ALPAKA_FN_INLINE float dPhiThreshold(TAcc const& acc,
0133                                                      float rt,
0134                                                      ModulesConst modules,
0135                                                      uint16_t moduleIndex,
0136                                                      const float ptCut,
0137                                                      float dPhi = 0,
0138                                                      float dz = 0) {
0139     // =================================================================
0140     // Various constants
0141     // =================================================================
0142     //mean of the horizontal layer position in y; treat this as R below
0143
0144     // =================================================================
0145     // Computing some components that make up the cut threshold
0146     // =================================================================
0147
0148     unsigned int iL = modules.layers()[moduleIndex] - 1;
0149     const float miniSlope = alpaka::math::asin(acc, alpaka::math::min(acc, rt * k2Rinv1GeVf / ptCut, kSinAlphaMax));
0150     const float rLayNominal =
0151         ((modules.subdets()[moduleIndex] == Barrel) ? kMiniRminMeanBarrel[iL] : kMiniRminMeanEndcap[iL]);
0152     const float miniPVoff = 0.1f / rLayNominal;
0153     const float miniMuls = ((modules.subdets()[moduleIndex] == Barrel) ? kMiniMulsPtScaleBarrel[iL] * 3.f / ptCut
0154                                                                        : kMiniMulsPtScaleEndcap[iL] * 3.f / ptCut);
0155     const bool isTilted = modules.subdets()[moduleIndex] == Barrel and modules.sides()[moduleIndex] != Center;
0156     //the lower module is sent in irrespective of its layer type. We need to fetch the drdz properly
0157
0158     float drdz;
0159     if (isTilted) {
0160       if (modules.moduleType()[moduleIndex] == PS and modules.moduleLayerType()[moduleIndex] == Strip) {
0161         drdz = modules.drdzs()[moduleIndex];
0162       } else {
0163         drdz = modules.drdzs()[modules.partnerModuleIndices()[moduleIndex]];
0164       }
0165     } else {
0166       drdz = 0;
0167     }
0168     const float miniTilt2 = ((isTilted) ? (0.5f * 0.5f) * (kPixelPSZpitch * kPixelPSZpitch) * (drdz * drdz) /
0169                                               (1.f + drdz * drdz) / moduleGapSize(modules, moduleIndex)
0170                                         : 0);
0171
0172     // Compute luminous region requirement for endcap
0173     const float miniLum = alpaka::math::abs(acc, dPhi * kDeltaZLum / dz);  // Balaji's new error
0174
0175     // =================================================================
0176     // Return the threshold value
0177     // =================================================================
0178     // Following condition is met if the module is central and flatly lying
0179     if (modules.subdets()[moduleIndex] == Barrel and modules.sides()[moduleIndex] == Center) {
0180       return miniSlope + alpaka::math::sqrt(acc, miniMuls * miniMuls + miniPVoff * miniPVoff);
0181     }
0182     // Following condition is met if the module is central and tilted
0183     else if (modules.subdets()[moduleIndex] == Barrel and
0184              modules.sides()[moduleIndex] != Center)  //all types of tilted modules
0185     {
0186       return miniSlope +
0187              alpaka::math::sqrt(acc, miniMuls * miniMuls + miniPVoff * miniPVoff + miniTilt2 * miniSlope * miniSlope);
0188     }
0189     // If not barrel, it is Endcap
0190     else {
0191       return miniSlope + alpaka::math::sqrt(acc, miniMuls * miniMuls + miniPVoff * miniPVoff + miniLum * miniLum);
0192     }
0193   }
0194
0195   template <typename TAcc>
0196   ALPAKA_FN_INLINE ALPAKA_FN_ACC void shiftStripHits(TAcc const& acc,
0197                                                      ModulesConst modules,
0198                                                      uint16_t lowerModuleIndex,
0199                                                      uint16_t upperModuleIndex,
0200                                                      unsigned int lowerHitIndex,
0201                                                      unsigned int upperHitIndex,
0202                                                      float* shiftedCoords,
0203                                                      float xLower,
0204                                                      float yLower,
0205                                                      float zLower,
0206                                                      float rtLower,
0207                                                      float xUpper,
0208                                                      float yUpper,
0209                                                      float zUpper,
0210                                                      float rtUpper) {
0211     // This is the strip shift scheme that is explained in http://uaf-10.t2.ucsd.edu/~phchang/talks/PhilipChang20190607_SDL_Update.pdf (see backup slides)
0212     // The main feature of this shifting is that the strip hits are shifted to be "aligned" in the line of sight from interaction point to the the pixel hit.
0213     // (since pixel hit is well defined in 3-d)
0214     // The strip hit is shifted along the strip detector to be placed in a guessed position where we think they would have actually crossed
0215     // The size of the radial direction shift due to module separation gap is computed in "radial" size, while the shift is done along the actual strip orientation
0216     // This means that there may be very very subtle edge effects coming from whether the strip hit is center of the module or the at the edge of the module
0217     // But this should be relatively minor effect
0218
0219     // dependent variables for this if statement
0220     // lowerModule
0221     // lowerHit
0222     // upperHit
0223     // endcapGeometry
0224     // tiltedGeometry
0225
0226     // Some variables relevant to the function
0227     float xp;       // pixel x (pixel hit x)
0228     float yp;       // pixel y (pixel hit y)
0229     float zp;       // pixel y (pixel hit y)
0230     float rtp;      // pixel y (pixel hit y)
0231     float xa;       // "anchor" x (the anchor position on the strip module plane from pixel hit)
0232     float ya;       // "anchor" y (the anchor position on the strip module plane from pixel hit)
0233     float xo;       // old x (before the strip hit is moved up or down)
0234     float yo;       // old y (before the strip hit is moved up or down)
0235     float xn;       // new x (after the strip hit is moved up or down)
0236     float yn;       // new y (after the strip hit is moved up or down)
0237     float abszn;    // new z in absolute value
0238     float zn;       // new z with the sign (+/-) accounted
0239     float angleA;   // in r-z plane the theta of the pixel hit in polar coordinate is the angleA
0240     float angleB;   // this is the angle of tilted module in r-z plane ("drdz"), for endcap this is 90 degrees
0241     bool isEndcap;  // If endcap, drdz = infinity
0242     float moduleSeparation;
0243     float drprime;    // The radial shift size in x-y plane projection
0244     float drprime_x;  // x-component of drprime
0245     float drprime_y;  // y-component of drprime
0246     const float& slope =
0247         modules.dxdys()[lowerModuleIndex];  // The slope of the possible strip hits for a given module in x-y plane
0248     float absArctanSlope;
0249     float angleM;  // the angle M is the angle of rotation of the module in x-y plane if the possible strip hits are along the x-axis, then angleM = 0, and if the possible strip hits are along y-axis angleM = 90 degrees
0250     float absdzprime;  // The distance between the two points after shifting
0251     const float& drdz_ = modules.drdzs()[lowerModuleIndex];
0252     // Assign hit pointers based on their hit type
0253     if (modules.moduleType()[lowerModuleIndex] == PS) {
0254       // TODO: This is somewhat of an mystery.... somewhat confused why this is the case
0255       if (modules.subdets()[lowerModuleIndex] == Barrel ? modules.moduleLayerType()[lowerModuleIndex] != Pixel
0256                                                         : modules.moduleLayerType()[lowerModuleIndex] == Pixel) {
0257         xo = xUpper;
0258         yo = yUpper;
0259         xp = xLower;
0260         yp = yLower;
0261         zp = zLower;
0262         rtp = rtLower;
0263       } else {
0264         xo = xLower;
0265         yo = yLower;
0266         xp = xUpper;
0267         yp = yUpper;
0268         zp = zUpper;
0269         rtp = rtUpper;
0270       }
0271     } else {
0272       xo = xUpper;
0273       yo = yUpper;
0274       xp = xLower;
0275       yp = yLower;
0276       zp = zLower;
0277       rtp = rtLower;
0278     }
0279
0280     // If it is endcap some of the math gets simplified (and also computers don't like infinities)
0281     isEndcap = modules.subdets()[lowerModuleIndex] == Endcap;
0282
0283     // NOTE: TODO: Keep in mind that the sin(atan) function can be simplified to something like x / sqrt(1 + x^2) and similar for cos
0284     // I am not sure how slow sin, atan, cos, functions are in c++. If x / sqrt(1 + x^2) are faster change this later to reduce arithmetic computation time
0285     angleA = alpaka::math::abs(acc, alpaka::math::atan(acc, rtp / zp));
0286     angleB =
0287         ((isEndcap)
0288              ? kPi / 2.f
0289              : alpaka::math::atan(
0290                    acc,
0291                    drdz_));  // The tilt module on the positive z-axis has negative drdz slope in r-z plane and vice versa
0292
0293     moduleSeparation = moduleGapSize(modules, lowerModuleIndex);
0294
0295     // Sign flips if the pixel is later layer
0296     if (modules.moduleType()[lowerModuleIndex] == PS and modules.moduleLayerType()[lowerModuleIndex] != Pixel) {
0297       moduleSeparation *= -1;
0298     }
0299
0300     drprime = (moduleSeparation / alpaka::math::sin(acc, angleA + angleB)) * alpaka::math::sin(acc, angleA);
0301
0302     // Compute arctan of the slope and take care of the slope = infinity case
0303     absArctanSlope =
0304         ((slope != kVerticalModuleSlope && edm::isFinite(slope)) ? fabs(alpaka::math::atan(acc, slope)) : kPi / 2.f);
0305
0306     // Depending on which quadrant the pixel hit lies, we define the angleM by shifting them slightly differently
0307     if (xp > 0 and yp > 0) {
0308       angleM = absArctanSlope;
0309     } else if (xp > 0 and yp < 0) {
0310       angleM = kPi - absArctanSlope;
0311     } else if (xp < 0 and yp < 0) {
0312       angleM = kPi + absArctanSlope;
0313     } else  // if (xp < 0 and yp > 0)
0314     {
0315       angleM = 2.f * kPi - absArctanSlope;
0316     }
0317
0318     // Then since the angleM sign is taken care of properly
0319     drprime_x = drprime * alpaka::math::sin(acc, angleM);
0320     drprime_y = drprime * alpaka::math::cos(acc, angleM);
0321
0322     // The new anchor position is
0323     xa = xp + drprime_x;
0324     ya = yp + drprime_y;
0325
0326     // Compute the new strip hit position (if the slope value is in special condition take care of the exceptions)
0327     if (slope == kVerticalModuleSlope ||
0328         edm::isNotFinite(slope))  // Designated for tilted module when the slope is infinity (module lying along y-axis)
0329     {
0330       xn = xa;  // New x point is simply where the anchor is
0331       yn = yo;  // No shift in y
0332     } else if (slope == 0) {
0333       xn = xo;  // New x point is simply where the anchor is
0334       yn = ya;  // No shift in y
0335     } else {
0336       xn = (slope * xa + (1.f / slope) * xo - ya + yo) / (slope + (1.f / slope));  // new xn
0337       yn = (xn - xa) * slope + ya;                                                 // new yn
0338     }
0339
0340     // Computing new Z position
0341     absdzprime = alpaka::math::abs(
0342         acc,
0343         moduleSeparation / alpaka::math::sin(acc, angleA + angleB) *
0344             alpaka::math::cos(
0345                 acc,
0346                 angleA));  // module separation sign is for shifting in radial direction for z-axis direction take care of the sign later
0347
0348     // Depending on which one as closer to the interactin point compute the new z wrt to the pixel properly
0349     if (modules.moduleLayerType()[lowerModuleIndex] == Pixel) {
0350       abszn = alpaka::math::abs(acc, zp) + absdzprime;
0351     } else {
0352       abszn = alpaka::math::abs(acc, zp) - absdzprime;
0353     }
0354
0355     zn = abszn * ((zp > 0) ? 1 : -1);  // Apply the sign of the zn
0356
0357     shiftedCoords[0] = xn;
0358     shiftedCoords[1] = yn;
0359     shiftedCoords[2] = zn;
0360   }
0361
0362   template <typename TAcc>
0363   ALPAKA_FN_ACC bool runMiniDoubletDefaultAlgoBarrel(TAcc const& acc,
0364                                                      ModulesConst modules,
0365                                                      uint16_t lowerModuleIndex,
0366                                                      uint16_t upperModuleIndex,
0367                                                      unsigned int lowerHitIndex,
0368                                                      unsigned int upperHitIndex,
0369                                                      float& dz,
0370                                                      float& dPhi,
0371                                                      float& dPhiChange,
0372                                                      float& shiftedX,
0373                                                      float& shiftedY,
0374                                                      float& shiftedZ,
0375                                                      float& noShiftedDphi,
0376                                                      float& noShiftedDphiChange,
0377                                                      float xLower,
0378                                                      float yLower,
0379                                                      float zLower,
0380                                                      float rtLower,
0381                                                      float xUpper,
0382                                                      float yUpper,
0383                                                      float zUpper,
0384                                                      float rtUpper,
0385                                                      const float ptCut) {
0386     dz = zLower - zUpper;
0387     const float dzCut = modules.moduleType()[lowerModuleIndex] == PS ? 2.f : 10.f;
0388     const float sign = ((dz > 0) - (dz < 0)) * ((zLower > 0) - (zLower < 0));
0389     const float invertedcrossercut = (alpaka::math::abs(acc, dz) > 2) * sign;
0390
0391     if ((alpaka::math::abs(acc, dz) >= dzCut) || (invertedcrossercut > 0))
0392       return false;
0393
0394     float miniCut = 0;
0395
0396     miniCut = modules.moduleLayerType()[lowerModuleIndex] == Pixel
0397                   ? dPhiThreshold(acc, rtLower, modules, lowerModuleIndex, ptCut)
0398                   : dPhiThreshold(acc, rtUpper, modules, lowerModuleIndex, ptCut);
0399
0400     // Cut #2: dphi difference
0401     // Ref to original code: https://github.com/slava77/cms-tkph2-ntuple/blob/184d2325147e6930030d3d1f780136bc2dd29ce6/doubletAnalysis.C#L3085
0402     float xn = 0.f, yn = 0.f;
0403     float shiftedRt2 = 0.f;
0404     if (modules.sides()[lowerModuleIndex] != Center)  // If barrel and not center it is tilted
0405     {
0406       // Shift the hits and calculate new xn, yn position
0407       float shiftedCoords[3];
0408       shiftStripHits(acc,
0409                      modules,
0410                      lowerModuleIndex,
0411                      upperModuleIndex,
0412                      lowerHitIndex,
0413                      upperHitIndex,
0414                      shiftedCoords,
0415                      xLower,
0416                      yLower,
0417                      zLower,
0418                      rtLower,
0419                      xUpper,
0420                      yUpper,
0421                      zUpper,
0422                      rtUpper);
0423       xn = shiftedCoords[0];
0424       yn = shiftedCoords[1];
0425
0426       // Lower or the upper hit needs to be modified depending on which one was actually shifted
0427       if (modules.moduleLayerType()[lowerModuleIndex] == Pixel) {
0428         shiftedX = xn;
0429         shiftedY = yn;
0430         shiftedZ = zUpper;
0431         shiftedRt2 = xn * xn + yn * yn;
0432
0433         dPhi = cms::alpakatools::deltaPhi(acc, xLower, yLower, shiftedX, shiftedY);  //function from Hit.cc
0434         noShiftedDphi = cms::alpakatools::deltaPhi(acc, xLower, yLower, xUpper, yUpper);
0435       } else {
0436         shiftedX = xn;
0437         shiftedY = yn;
0438         shiftedZ = zLower;
0439         shiftedRt2 = xn * xn + yn * yn;
0440         dPhi = cms::alpakatools::deltaPhi(acc, shiftedX, shiftedY, xUpper, yUpper);
0441         noShiftedDphi = cms::alpakatools::deltaPhi(acc, xLower, yLower, xUpper, yUpper);
0442       }
0443     } else {
0444       shiftedX = 0.f;
0445       shiftedY = 0.f;
0446       shiftedZ = 0.f;
0447       shiftedRt2 = 0.f;
0448       dPhi = cms::alpakatools::deltaPhi(acc, xLower, yLower, xUpper, yUpper);
0449       noShiftedDphi = dPhi;
0450     }
0451
0452     if (alpaka::math::abs(acc, dPhi) >= miniCut)
0453       return false;
0454
0455     // Cut #3: The dphi change going from lower Hit to upper Hit
0456     // Ref to original code: https://github.com/slava77/cms-tkph2-ntuple/blob/184d2325147e6930030d3d1f780136bc2dd29ce6/doubletAnalysis.C#L3076
0457     if (modules.sides()[lowerModuleIndex] != Center) {
0458       // When it is tilted, use the new shifted positions
0459       if (modules.moduleLayerType()[lowerModuleIndex] != Pixel) {
0460         // dPhi Change should be calculated so that the upper hit has higher rt.
0461         // In principle, this kind of check rt_lower < rt_upper should not be necessary because the hit shifting should have taken care of this.
0462         // (i.e. the strip hit is shifted to be aligned in the line of sight from interaction point to pixel hit of PS module guaranteeing rt ordering)
0463         // But I still placed this check for safety. (TODO: After checking explicitly if not needed remove later?)
0464         // setdeltaPhiChange(lowerHit.rt() < upperHitMod.rt() ? lowerHit.deltaPhiChange(upperHitMod) : upperHitMod.deltaPhiChange(lowerHit));
0465
0466         dPhiChange = (rtLower * rtLower < shiftedRt2) ? deltaPhiChange(acc, xLower, yLower, shiftedX, shiftedY)
0467                                                       : deltaPhiChange(acc, shiftedX, shiftedY, xLower, yLower);
0468         noShiftedDphiChange = rtLower < rtUpper ? deltaPhiChange(acc, xLower, yLower, xUpper, yUpper)
0469                                                 : deltaPhiChange(acc, xUpper, yUpper, xLower, yLower);
0470       } else {
0471         // dPhi Change should be calculated so that the upper hit has higher rt.
0472         // In principle, this kind of check rt_lower < rt_upper should not be necessary because the hit shifting should have taken care of this.
0473         // (i.e. the strip hit is shifted to be aligned in the line of sight from interaction point to pixel hit of PS module guaranteeing rt ordering)
0474         // But I still placed this check for safety. (TODO: After checking explicitly if not needed remove later?)
0475
0476         dPhiChange = (shiftedRt2 < rtUpper * rtUpper) ? deltaPhiChange(acc, shiftedX, shiftedY, xUpper, yUpper)
0477                                                       : deltaPhiChange(acc, xUpper, yUpper, shiftedX, shiftedY);
0478         noShiftedDphiChange = rtLower < rtUpper ? deltaPhiChange(acc, xLower, yLower, xUpper, yUpper)
0479                                                 : deltaPhiChange(acc, xUpper, yUpper, xLower, yLower);
0480       }
0481     } else {
0482       // When it is flat lying module, whichever is the lowerSide will always have rt lower
0483       dPhiChange = deltaPhiChange(acc, xLower, yLower, xUpper, yUpper);
0484       noShiftedDphiChange = dPhiChange;
0485     }
0486
0487     return alpaka::math::abs(acc, dPhiChange) < miniCut;
0488   }
0489
0490   template <typename TAcc>
0491   ALPAKA_FN_ACC bool runMiniDoubletDefaultAlgoEndcap(TAcc const& acc,
0492                                                      ModulesConst modules,
0493                                                      uint16_t lowerModuleIndex,
0494                                                      uint16_t upperModuleIndex,
0495                                                      unsigned int lowerHitIndex,
0496                                                      unsigned int upperHitIndex,
0497                                                      float& drt,
0498                                                      float& dPhi,
0499                                                      float& dPhiChange,
0500                                                      float& shiftedX,
0501                                                      float& shiftedY,
0502                                                      float& shiftedZ,
0503                                                      float& noShiftedDphi,
0504                                                      float& noShiftedDphichange,
0505                                                      float xLower,
0506                                                      float yLower,
0507                                                      float zLower,
0508                                                      float rtLower,
0509                                                      float xUpper,
0510                                                      float yUpper,
0511                                                      float zUpper,
0512                                                      float rtUpper,
0513                                                      const float ptCut) {
0514     // There are series of cuts that applies to mini-doublet in a "endcap" region
0515     // Cut #1 : dz cut. The dz difference can't be larger than 1cm. (max separation is 4mm for modules in the endcap)
0516     // Ref to original code: https://github.com/slava77/cms-tkph2-ntuple/blob/184d2325147e6930030d3d1f780136bc2dd29ce6/doubletAnalysis.C#L3093
0517     // For PS module in case when it is tilted a different dz (after the strip hit shift) is calculated later.
0518
0519     float dz = zLower - zUpper;  // Not const since later it might change depending on the type of module
0520
0521     const float dzCut = 1.f;
0522
0523     if (alpaka::math::abs(acc, dz) >= dzCut)
0524       return false;
0525     // Cut #2 : drt cut. The dz difference can't be larger than 1cm. (max separation is 4mm for modules in the endcap)
0526     // Ref to original code: https://github.com/slava77/cms-tkph2-ntuple/blob/184d2325147e6930030d3d1f780136bc2dd29ce6/doubletAnalysis.C#L3100
0527     const float drtCut = modules.moduleType()[lowerModuleIndex] == PS ? 2.f : 10.f;
0528     drt = rtLower - rtUpper;
0529     if (alpaka::math::abs(acc, drt) >= drtCut)
0530       return false;
0531     // The new scheme shifts strip hits to be "aligned" along the line of sight from interaction point to the pixel hit (if it is PS modules)
0532     float xn = 0, yn = 0, zn = 0;
0533
0534     float shiftedCoords[3];
0535     shiftStripHits(acc,
0536                    modules,
0537                    lowerModuleIndex,
0538                    upperModuleIndex,
0539                    lowerHitIndex,
0540                    upperHitIndex,
0541                    shiftedCoords,
0542                    xLower,
0543                    yLower,
0544                    zLower,
0545                    rtLower,
0546                    xUpper,
0547                    yUpper,
0548                    zUpper,
0549                    rtUpper);
0550
0551     xn = shiftedCoords[0];
0552     yn = shiftedCoords[1];
0553     zn = shiftedCoords[2];
0554
0555     if (modules.moduleType()[lowerModuleIndex] == PS) {
0556       // Appropriate lower or upper hit is modified after checking which one was actually shifted
0557       if (modules.moduleLayerType()[lowerModuleIndex] == Pixel) {
0558         shiftedX = xn;
0559         shiftedY = yn;
0560         shiftedZ = zUpper;
0561         dPhi = cms::alpakatools::deltaPhi(acc, xLower, yLower, shiftedX, shiftedY);
0562         noShiftedDphi = cms::alpakatools::deltaPhi(acc, xLower, yLower, xUpper, yUpper);
0563       } else {
0564         shiftedX = xn;
0565         shiftedY = yn;
0566         shiftedZ = zLower;
0567         dPhi = cms::alpakatools::deltaPhi(acc, shiftedX, shiftedY, xUpper, yUpper);
0568         noShiftedDphi = cms::alpakatools::deltaPhi(acc, xLower, yLower, xUpper, yUpper);
0569       }
0570     } else {
0571       shiftedX = xn;
0572       shiftedY = yn;
0573       shiftedZ = zUpper;
0574       dPhi = cms::alpakatools::deltaPhi(acc, xLower, yLower, xn, yn);
0575       noShiftedDphi = cms::alpakatools::deltaPhi(acc, xLower, yLower, xUpper, yUpper);
0576     }
0577
0578     // dz needs to change if it is a PS module where the strip hits are shifted in order to properly account for the case when a tilted module falls under "endcap logic"
0579     // if it was an endcap it will have zero effect
0580     if (modules.moduleType()[lowerModuleIndex] == PS) {
0581       dz = modules.moduleLayerType()[lowerModuleIndex] == Pixel ? zLower - zn : zUpper - zn;
0582     }
0583
0584     float miniCut = 0;
0585     miniCut = modules.moduleLayerType()[lowerModuleIndex] == Pixel
0586                   ? dPhiThreshold(acc, rtLower, modules, lowerModuleIndex, ptCut, dPhi, dz)
0587                   : dPhiThreshold(acc, rtUpper, modules, lowerModuleIndex, ptCut, dPhi, dz);
0588
0589     if (alpaka::math::abs(acc, dPhi) >= miniCut)
0590       return false;
0591
0592     // Cut #4: Another cut on the dphi after some modification
0593     // Ref to original code: https://github.com/slava77/cms-tkph2-ntuple/blob/184d2325147e6930030d3d1f780136bc2dd29ce6/doubletAnalysis.C#L3119-L3124
0594
0595     float dzFrac = alpaka::math::abs(acc, dz) / alpaka::math::abs(acc, zLower);
0596     dPhiChange = dPhi / dzFrac * (1.f + dzFrac);
0597     noShiftedDphichange = noShiftedDphi / dzFrac * (1.f + dzFrac);
0598
0599     return alpaka::math::abs(acc, dPhiChange) < miniCut;
0600   }
0601
0602   template <typename TAcc>
0603   ALPAKA_FN_ACC bool runMiniDoubletDefaultAlgo(TAcc const& acc,
0604                                                ModulesConst modules,
0605                                                uint16_t lowerModuleIndex,
0606                                                uint16_t upperModuleIndex,
0607                                                unsigned int lowerHitIndex,
0608                                                unsigned int upperHitIndex,
0609                                                float& dz,
0610                                                float& dPhi,
0611                                                float& dPhiChange,
0612                                                float& shiftedX,
0613                                                float& shiftedY,
0614                                                float& shiftedZ,
0615                                                float& noShiftedDphi,
0616                                                float& noShiftedDphiChange,
0617                                                float xLower,
0618                                                float yLower,
0619                                                float zLower,
0620                                                float rtLower,
0621                                                float xUpper,
0622                                                float yUpper,
0623                                                float zUpper,
0624                                                float rtUpper,
0625                                                const float ptCut) {
0626     if (modules.subdets()[lowerModuleIndex] == Barrel) {
0627       return runMiniDoubletDefaultAlgoBarrel(acc,
0628                                              modules,
0629                                              lowerModuleIndex,
0630                                              upperModuleIndex,
0631                                              lowerHitIndex,
0632                                              upperHitIndex,
0633                                              dz,
0634                                              dPhi,
0635                                              dPhiChange,
0636                                              shiftedX,
0637                                              shiftedY,
0638                                              shiftedZ,
0639                                              noShiftedDphi,
0640                                              noShiftedDphiChange,
0641                                              xLower,
0642                                              yLower,
0643                                              zLower,
0644                                              rtLower,
0645                                              xUpper,
0646                                              yUpper,
0647                                              zUpper,
0648                                              rtUpper,
0649                                              ptCut);
0650     } else {
0651       return runMiniDoubletDefaultAlgoEndcap(acc,
0652                                              modules,
0653                                              lowerModuleIndex,
0654                                              upperModuleIndex,
0655                                              lowerHitIndex,
0656                                              upperHitIndex,
0657                                              dz,
0658                                              dPhi,
0659                                              dPhiChange,
0660                                              shiftedX,
0661                                              shiftedY,
0662                                              shiftedZ,
0663                                              noShiftedDphi,
0664                                              noShiftedDphiChange,
0665                                              xLower,
0666                                              yLower,
0667                                              zLower,
0668                                              rtLower,
0669                                              xUpper,
0670                                              yUpper,
0671                                              zUpper,
0672                                              rtUpper,
0673                                              ptCut);
0674     }
0675   }
0676
0677   struct CreateMiniDoublets {
0678     ALPAKA_FN_ACC void operator()(Acc2D const& acc,
0679                                   ModulesConst modules,
0680                                   HitsBaseConst hitsBase,
0681                                   HitsExtendedConst hitsExtended,
0682                                   HitsRangesConst hitsRanges,
0683                                   MiniDoublets mds,
0684                                   MiniDoubletsOccupancy mdsOccupancy,
0685                                   ObjectRangesConst ranges,
0686                                   const float ptCut) const {
0687       for (uint16_t lowerModuleIndex : cms::alpakatools::uniform_elements_y(acc, modules.nLowerModules())) {
0688         uint16_t upperModuleIndex = modules.partnerModuleIndices()[lowerModuleIndex];
0689         int nLowerHits = hitsRanges.hitRangesnLower()[lowerModuleIndex];
0690         int nUpperHits = hitsRanges.hitRangesnUpper()[lowerModuleIndex];
0691         if (hitsRanges.hitRangesLower()[lowerModuleIndex] == -1)
0692           continue;
0693         unsigned int upHitArrayIndex = hitsRanges.hitRangesUpper()[lowerModuleIndex];
0694         unsigned int loHitArrayIndex = hitsRanges.hitRangesLower()[lowerModuleIndex];
0695         int limit = nUpperHits * nLowerHits;
0696
0697         for (int hitIndex : cms::alpakatools::uniform_elements_x(acc, limit)) {
0698           int lowerHitIndex = hitIndex / nUpperHits;
0699           int upperHitIndex = hitIndex % nUpperHits;
0700           if (upperHitIndex >= nUpperHits)
0701             continue;
0702           if (lowerHitIndex >= nLowerHits)
0703             continue;
0704           unsigned int lowerHitArrayIndex = loHitArrayIndex + lowerHitIndex;
0705           float xLower = hitsBase.xs()[lowerHitArrayIndex];
0706           float yLower = hitsBase.ys()[lowerHitArrayIndex];
0707           float zLower = hitsBase.zs()[lowerHitArrayIndex];
0708           float rtLower = hitsExtended.rts()[lowerHitArrayIndex];
0709           unsigned int upperHitArrayIndex = upHitArrayIndex + upperHitIndex;
0710           float xUpper = hitsBase.xs()[upperHitArrayIndex];
0711           float yUpper = hitsBase.ys()[upperHitArrayIndex];
0712           float zUpper = hitsBase.zs()[upperHitArrayIndex];
0713           float rtUpper = hitsExtended.rts()[upperHitArrayIndex];
0714
0715           float dz, dphi, dphichange, shiftedX, shiftedY, shiftedZ, noShiftedDphi, noShiftedDphiChange;
0716           bool success = runMiniDoubletDefaultAlgo(acc,
0717                                                    modules,
0718                                                    lowerModuleIndex,
0719                                                    upperModuleIndex,
0720                                                    lowerHitArrayIndex,
0721                                                    upperHitArrayIndex,
0722                                                    dz,
0723                                                    dphi,
0724                                                    dphichange,
0725                                                    shiftedX,
0726                                                    shiftedY,
0727                                                    shiftedZ,
0728                                                    noShiftedDphi,
0729                                                    noShiftedDphiChange,
0730                                                    xLower,
0731                                                    yLower,
0732                                                    zLower,
0733                                                    rtLower,
0734                                                    xUpper,
0735                                                    yUpper,
0736                                                    zUpper,
0737                                                    rtUpper,
0738                                                    ptCut);
0739           if (success) {
0740             int totOccupancyMDs = alpaka::atomicAdd(
0741                 acc, &mdsOccupancy.totOccupancyMDs()[lowerModuleIndex], 1u, alpaka::hierarchy::Threads{});
0742             if (totOccupancyMDs >= (ranges.miniDoubletModuleOccupancy()[lowerModuleIndex])) {
0743 #ifdef WARNINGS
0744               printf(
0745                   "Mini-doublet excess alert! Module index = %d, Occupancy = %d\n", lowerModuleIndex, totOccupancyMDs);
0746 #endif
0747             } else {
0748               int mdModuleIndex =
0749                   alpaka::atomicAdd(acc, &mdsOccupancy.nMDs()[lowerModuleIndex], 1u, alpaka::hierarchy::Threads{});
0750               unsigned int mdIndex = ranges.miniDoubletModuleIndices()[lowerModuleIndex] + mdModuleIndex;
0751
0752               addMDToMemory(acc,
0753                             mds,
0754                             hitsBase,
0755                             hitsExtended,
0756                             modules,
0757                             lowerHitArrayIndex,
0758                             upperHitArrayIndex,
0759                             lowerModuleIndex,
0760                             dz,
0761                             dphi,
0762                             dphichange,
0763                             shiftedX,
0764                             shiftedY,
0765                             shiftedZ,
0766                             noShiftedDphi,
0767                             noShiftedDphiChange,
0768                             mdIndex);
0769             }
0770           }
0771         }
0772       }
0773     }
0774   };
0775
0776   // Helper function to determine eta bin for occupancies
0777   ALPAKA_FN_ACC ALPAKA_FN_INLINE int getEtaBin(const float module_eta) {
0778     if (module_eta < 0.75f)
0779       return 0;
0780     else if (module_eta < 1.5f)
0781       return 1;
0782     else if (module_eta < 2.25f)
0783       return 2;
0784     else if (module_eta < 3.0f)
0785       return 3;
0786     return -1;
0787   }
0788
0789   // Helper function to determine category number for occupancies
0790   ALPAKA_FN_ACC ALPAKA_FN_INLINE int getCategoryNumber(const short module_layers,
0791                                                        const short module_subdets,
0792                                                        const short module_rings) {
0793     if (module_subdets == Barrel) {
0794       return (module_layers <= 3) ? 0 : 1;
0795     } else if (module_subdets == Endcap) {
0796       if (module_layers <= 2) {
0797         return (module_rings >= 11) ? 2 : 3;
0798       } else {
0799         return (module_rings >= 8) ? 2 : 3;
0800       }
0801     }
0802     return -1;
0803   }
0804
0805   struct CreateMDArrayRangesGPU {
0806     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
0807                                   ModulesConst modules,
0808                                   HitsRangesConst hitsRanges,
0809                                   ObjectRanges ranges,
0810                                   const float ptCut) const {
0811       // implementation is 1D with a single block
0812       ALPAKA_ASSERT_ACC((alpaka::getWorkDiv<alpaka::Grid, alpaka::Blocks>(acc)[0] == 1));
0813
0814       // Declare variables in shared memory and set to 0
0815       int& nTotalMDs = alpaka::declareSharedVar<int, __COUNTER__>(acc);
0816       if (cms::alpakatools::once_per_block(acc)) {
0817         nTotalMDs = 0;
0818       }
0819       alpaka::syncBlockThreads(acc);
0820
0821       // Occupancy matrix for 0.8 GeV pT Cut
0822       constexpr int p08_occupancy_matrix[4][4] = {
0823           {49, 42, 37, 41},  // category 0
0824           {100, 100, 0, 0},  // category 1
0825           {0, 16, 19, 0},    // category 2
0826           {0, 14, 20, 25}    // category 3
0827       };
0828
0829       // Occupancy matrix for 0.6 GeV pT Cut, 99.99%
0830       constexpr int p06_occupancy_matrix[4][4] = {
0831           {60, 57, 54, 48},  // category 0
0832           {259, 195, 0, 0},  // category 1
0833           {0, 23, 28, 0},    // category 2
0834           {0, 25, 25, 33}    // category 3
0835       };
0836
0837       // Select the appropriate occupancy matrix based on ptCut
0838       const auto& occupancy_matrix = (ptCut < 0.8f) ? p06_occupancy_matrix : p08_occupancy_matrix;
0839
0840       for (uint16_t i : cms::alpakatools::uniform_elements(acc, modules.nLowerModules())) {
0841         const int nLower = hitsRanges.hitRangesnLower()[i];
0842         const int nUpper = hitsRanges.hitRangesnUpper()[i];
0843         const int dynamicMDs = nLower * nUpper;
0844
0845         // Matrix-based cap
0846         short module_layers = modules.layers()[i];
0847         short module_subdets = modules.subdets()[i];
0848         short module_rings = modules.rings()[i];
0849         float module_eta = alpaka::math::abs(acc, modules.eta()[i]);
0850
0851         int category_number = getCategoryNumber(module_layers, module_subdets, module_rings);
0852         int eta_number = getEtaBin(module_eta);
0853
0854 #ifdef WARNINGS
0855         if (category_number == -1 || eta_number == -1) {
0856           printf("Unhandled case in createMDArrayRangesGPU! Module index = %i\n", i);
0857         }
0858 #endif
0859
0860         int occupancy = (category_number != -1 && eta_number != -1)
0861                             ? alpaka::math::min(acc, dynamicMDs, occupancy_matrix[category_number][eta_number])
0862                             : 0;
0863         unsigned int nTotMDs = alpaka::atomicAdd(acc, &nTotalMDs, occupancy, alpaka::hierarchy::Threads{});
0864
0865         ranges.miniDoubletModuleIndices()[i] = nTotMDs;
0866         ranges.miniDoubletModuleOccupancy()[i] = occupancy;
0867       }
0868
0869       // Wait for all threads to finish before reporting final values
0870       alpaka::syncBlockThreads(acc);
0871       if (cms::alpakatools::once_per_block(acc)) {
0872         ranges.miniDoubletModuleIndices()[modules.nLowerModules()] = nTotalMDs;
0873         ranges.nTotalMDs() = nTotalMDs;
0874       }
0875     }
0876   };
0877
0878   struct AddMiniDoubletRangesToEventExplicit {
0879     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
0880                                   ModulesConst modules,
0881                                   MiniDoubletsOccupancy mdsOccupancy,
0882                                   ObjectRanges ranges,
0883                                   HitsRangesConst hitsRanges) const {
0884       // implementation is 1D with a single block
0885       ALPAKA_ASSERT_ACC((alpaka::getWorkDiv<alpaka::Grid, alpaka::Blocks>(acc)[0] == 1));
0886
0887       for (uint16_t i : cms::alpakatools::uniform_elements(acc, modules.nLowerModules())) {
0888         if (mdsOccupancy.nMDs()[i] == 0 or hitsRanges.hitRanges()[i][0] == -1) {
0889           ranges.mdRanges()[i][0] = -1;
0890           ranges.mdRanges()[i][1] = -1;
0891         } else {
0892           ranges.mdRanges()[i][0] = ranges.miniDoubletModuleIndices()[i];
0893           ranges.mdRanges()[i][1] = ranges.miniDoubletModuleIndices()[i] + mdsOccupancy.nMDs()[i] - 1;
0894         }
0895       }
0896     }
0897   };
0898 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE::lst
0899
0900 #endif