Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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