Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-05-09 22:40:21 

0001 #ifndef RecoTracker_LSTCore_src_alpaka_PixelTriplet_h
0002 #define RecoTracker_LSTCore_src_alpaka_PixelTriplet_h
0003 
0004 #include "RecoTracker/LSTCore/interface/alpaka/Common.h"
0005 #include "RecoTracker/LSTCore/interface/ModulesSoA.h"
0006 #include "RecoTracker/LSTCore/interface/ObjectRangesSoA.h"
0007 #include "RecoTracker/LSTCore/interface/MiniDoubletsSoA.h"
0008 #include "RecoTracker/LSTCore/interface/PixelTripletsSoA.h"
0009 #include "RecoTracker/LSTCore/interface/PixelSegmentsSoA.h"
0010 #include "RecoTracker/LSTCore/interface/QuintupletsSoA.h"
0011 #include "RecoTracker/LSTCore/interface/SegmentsSoA.h"
0012 #include "RecoTracker/LSTCore/interface/TripletsSoA.h"
0013 
0014 #include "Quintuplet.h"
0015 
0016 namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
0017 
0018   template <typename TAcc>
0019   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runTripletDefaultAlgoPPBB(TAcc const& acc,
0020                                                                 ModulesConst modules,
0021                                                                 ObjectRangesConst ranges,
0022                                                                 MiniDoubletsConst mds,
0023                                                                 SegmentsConst segments,
0024                                                                 PixelSeedsConst pixelSeeds,
0025                                                                 uint16_t pixelModuleIndex,
0026                                                                 uint16_t outerInnerLowerModuleIndex,
0027                                                                 uint16_t outerOuterLowerModuleIndex,
0028                                                                 unsigned int innerSegmentIndex,
0029                                                                 unsigned int outerSegmentIndex,
0030                                                                 unsigned int firstMDIndex,
0031                                                                 unsigned int secondMDIndex,
0032                                                                 unsigned int thirdMDIndex,
0033                                                                 unsigned int fourthMDIndex,
0034                                                                 const float ptCut);
0035   template <typename TAcc>
0036   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runTripletDefaultAlgoPPEE(TAcc const& acc,
0037                                                                 ModulesConst modules,
0038                                                                 ObjectRangesConst ranges,
0039                                                                 MiniDoubletsConst mds,
0040                                                                 SegmentsConst segments,
0041                                                                 PixelSeedsConst pixelSeeds,
0042                                                                 uint16_t pixelModuleIndex,
0043                                                                 uint16_t outerInnerLowerModuleIndex,
0044                                                                 uint16_t outerOuterLowerModuleIndex,
0045                                                                 unsigned int innerSegmentIndex,
0046                                                                 unsigned int outerSegmentIndex,
0047                                                                 unsigned int firstMDIndex,
0048                                                                 unsigned int secondMDIndex,
0049                                                                 unsigned int thirdMDIndex,
0050                                                                 unsigned int fourthMDIndex,
0051                                                                 const float ptCut);
0052 
0053   ALPAKA_FN_ACC ALPAKA_FN_INLINE void addPixelTripletToMemory(MiniDoubletsConst mds,
0054                                                               SegmentsConst segments,
0055                                                               TripletsConst triplets,
0056                                                               PixelTriplets pixelTriplets,
0057                                                               unsigned int pixelSegmentIndex,
0058                                                               unsigned int tripletIndex,
0059                                                               float pixelRadius,
0060                                                               float tripletRadius,
0061                                                               float centerX,
0062                                                               float centerY,
0063                                                               float rPhiChiSquared,
0064                                                               float rPhiChiSquaredInwards,
0065                                                               float rzChiSquared,
0066                                                               unsigned int pixelTripletIndex,
0067                                                               float pt,
0068                                                               float eta,
0069                                                               float phi,
0070                                                               float eta_pix,
0071                                                               float phi_pix,
0072                                                               float pixelRadiusError,
0073                                                               float score) {
0074     pixelTriplets.pixelSegmentIndices()[pixelTripletIndex] = pixelSegmentIndex;
0075     pixelTriplets.tripletIndices()[pixelTripletIndex] = tripletIndex;
0076     pixelTriplets.pixelRadius()[pixelTripletIndex] = __F2H(pixelRadius);
0077     pixelTriplets.pixelRadiusError()[pixelTripletIndex] = __F2H(pixelRadiusError);
0078     pixelTriplets.tripletRadius()[pixelTripletIndex] = __F2H(tripletRadius);
0079     pixelTriplets.pt()[pixelTripletIndex] = __F2H(pt);
0080     pixelTriplets.eta()[pixelTripletIndex] = __F2H(eta);
0081     pixelTriplets.phi()[pixelTripletIndex] = __F2H(phi);
0082     pixelTriplets.eta_pix()[pixelTripletIndex] = __F2H(eta_pix);
0083     pixelTriplets.phi_pix()[pixelTripletIndex] = __F2H(phi_pix);
0084     pixelTriplets.isDup()[pixelTripletIndex] = false;
0085     pixelTriplets.score()[pixelTripletIndex] = __F2H(score);
0086 
0087     pixelTriplets.centerX()[pixelTripletIndex] = __F2H(centerX);
0088     pixelTriplets.centerY()[pixelTripletIndex] = __F2H(centerY);
0089     pixelTriplets.logicalLayers()[pixelTripletIndex][0] = 0;
0090     pixelTriplets.logicalLayers()[pixelTripletIndex][1] = 0;
0091     pixelTriplets.logicalLayers()[pixelTripletIndex][2] = triplets.logicalLayers()[tripletIndex][0];
0092     pixelTriplets.logicalLayers()[pixelTripletIndex][3] = triplets.logicalLayers()[tripletIndex][1];
0093     pixelTriplets.logicalLayers()[pixelTripletIndex][4] = triplets.logicalLayers()[tripletIndex][2];
0094 
0095     pixelTriplets.lowerModuleIndices()[pixelTripletIndex][0] = segments.innerLowerModuleIndices()[pixelSegmentIndex];
0096     pixelTriplets.lowerModuleIndices()[pixelTripletIndex][1] = segments.outerLowerModuleIndices()[pixelSegmentIndex];
0097     pixelTriplets.lowerModuleIndices()[pixelTripletIndex][2] = triplets.lowerModuleIndices()[tripletIndex][0];
0098     pixelTriplets.lowerModuleIndices()[pixelTripletIndex][3] = triplets.lowerModuleIndices()[tripletIndex][1];
0099     pixelTriplets.lowerModuleIndices()[pixelTripletIndex][4] = triplets.lowerModuleIndices()[tripletIndex][2];
0100 
0101     unsigned int pixelInnerMD = segments.mdIndices()[pixelSegmentIndex][0];
0102     unsigned int pixelOuterMD = segments.mdIndices()[pixelSegmentIndex][1];
0103 
0104     pixelTriplets.hitIndices()[pixelTripletIndex][0] = mds.anchorHitIndices()[pixelInnerMD];
0105     pixelTriplets.hitIndices()[pixelTripletIndex][1] = mds.outerHitIndices()[pixelInnerMD];
0106     pixelTriplets.hitIndices()[pixelTripletIndex][2] = mds.anchorHitIndices()[pixelOuterMD];
0107     pixelTriplets.hitIndices()[pixelTripletIndex][3] = mds.outerHitIndices()[pixelOuterMD];
0108 
0109     pixelTriplets.hitIndices()[pixelTripletIndex][4] = triplets.hitIndices()[tripletIndex][0];
0110     pixelTriplets.hitIndices()[pixelTripletIndex][5] = triplets.hitIndices()[tripletIndex][1];
0111     pixelTriplets.hitIndices()[pixelTripletIndex][6] = triplets.hitIndices()[tripletIndex][2];
0112     pixelTriplets.hitIndices()[pixelTripletIndex][7] = triplets.hitIndices()[tripletIndex][3];
0113     pixelTriplets.hitIndices()[pixelTripletIndex][8] = triplets.hitIndices()[tripletIndex][4];
0114     pixelTriplets.hitIndices()[pixelTripletIndex][9] = triplets.hitIndices()[tripletIndex][5];
0115     pixelTriplets.rPhiChiSquared()[pixelTripletIndex] = rPhiChiSquared;
0116     pixelTriplets.rPhiChiSquaredInwards()[pixelTripletIndex] = rPhiChiSquaredInwards;
0117     pixelTriplets.rzChiSquared()[pixelTripletIndex] = rzChiSquared;
0118   }
0119 
0120   template <typename TAcc>
0121   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runPixelTrackletDefaultAlgopT3(TAcc const& acc,
0122                                                                      ModulesConst modules,
0123                                                                      ObjectRangesConst ranges,
0124                                                                      MiniDoubletsConst mds,
0125                                                                      SegmentsConst segments,
0126                                                                      PixelSeedsConst pixelSeeds,
0127                                                                      uint16_t pixelLowerModuleIndex,
0128                                                                      uint16_t outerInnerLowerModuleIndex,
0129                                                                      uint16_t outerOuterLowerModuleIndex,
0130                                                                      unsigned int innerSegmentIndex,
0131                                                                      unsigned int outerSegmentIndex,
0132                                                                      const float ptCut) {
0133     short outerInnerLowerModuleSubdet = modules.subdets()[outerInnerLowerModuleIndex];
0134     short outerOuterLowerModuleSubdet = modules.subdets()[outerOuterLowerModuleIndex];
0135 
0136     unsigned int firstMDIndex = segments.mdIndices()[innerSegmentIndex][0];
0137     unsigned int secondMDIndex = segments.mdIndices()[innerSegmentIndex][1];
0138 
0139     unsigned int thirdMDIndex = segments.mdIndices()[outerSegmentIndex][0];
0140     unsigned int fourthMDIndex = segments.mdIndices()[outerSegmentIndex][1];
0141 
0142     if (outerInnerLowerModuleSubdet == Barrel and
0143         (outerOuterLowerModuleSubdet == Barrel or outerOuterLowerModuleSubdet == Endcap)) {
0144       return runTripletDefaultAlgoPPBB(acc,
0145                                        modules,
0146                                        ranges,
0147                                        mds,
0148                                        segments,
0149                                        pixelSeeds,
0150                                        pixelLowerModuleIndex,
0151                                        outerInnerLowerModuleIndex,
0152                                        outerOuterLowerModuleIndex,
0153                                        innerSegmentIndex,
0154                                        outerSegmentIndex,
0155                                        firstMDIndex,
0156                                        secondMDIndex,
0157                                        thirdMDIndex,
0158                                        fourthMDIndex,
0159                                        ptCut);
0160     } else if (outerInnerLowerModuleSubdet == Endcap and outerOuterLowerModuleSubdet == Endcap) {
0161       return runTripletDefaultAlgoPPEE(acc,
0162                                        modules,
0163                                        ranges,
0164                                        mds,
0165                                        segments,
0166                                        pixelSeeds,
0167                                        pixelLowerModuleIndex,
0168                                        outerInnerLowerModuleIndex,
0169                                        outerOuterLowerModuleIndex,
0170                                        innerSegmentIndex,
0171                                        outerSegmentIndex,
0172                                        firstMDIndex,
0173                                        secondMDIndex,
0174                                        thirdMDIndex,
0175                                        fourthMDIndex,
0176                                        ptCut);
0177     }
0178     return false;
0179   }
0180 
0181   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passPT3RZChiSquaredCuts(ModulesConst modules,
0182                                                               uint16_t lowerModuleIndex1,
0183                                                               uint16_t lowerModuleIndex2,
0184                                                               uint16_t lowerModuleIndex3,
0185                                                               float rzChiSquared) {
0186     const int layer1 =
0187         modules.layers()[lowerModuleIndex1] + 6 * (modules.subdets()[lowerModuleIndex1] == Endcap) +
0188         5 * (modules.subdets()[lowerModuleIndex1] == Endcap and modules.moduleType()[lowerModuleIndex1] == TwoS);
0189     const int layer2 =
0190         modules.layers()[lowerModuleIndex2] + 6 * (modules.subdets()[lowerModuleIndex2] == Endcap) +
0191         5 * (modules.subdets()[lowerModuleIndex2] == Endcap and modules.moduleType()[lowerModuleIndex2] == TwoS);
0192     const int layer3 =
0193         modules.layers()[lowerModuleIndex3] + 6 * (modules.subdets()[lowerModuleIndex3] == Endcap) +
0194         5 * (modules.subdets()[lowerModuleIndex3] == Endcap and modules.moduleType()[lowerModuleIndex3] == TwoS);
0195 
0196     if (layer1 == 8 and layer2 == 9 and layer3 == 10) {
0197       return rzChiSquared < 13.6067f;
0198     } else if (layer1 == 8 and layer2 == 9 and layer3 == 15) {
0199       return rzChiSquared < 5.5953f;
0200     } else if (layer1 == 7 and layer2 == 8 and layer3 == 9) {
0201       return rzChiSquared < 3.9263f;
0202     }
0203     /*
0204     else if(layer1 == 7 and layer2 == 8 and layer3 == 14)
0205     {
0206         // PS+PS+2S in endcap layers 1+2+3, which is not really feasible in the current geometry,
0207         // without skipping barrel layers 1 and 2 (not allowed by algorithm logic).
0208     }
0209     */
0210     else if (layer1 == 1 and layer2 == 2 and layer3 == 3) {
0211       return rzChiSquared < 9.4377f;
0212     } else if (layer1 == 1 and layer2 == 2 and layer3 == 7) {
0213       return rzChiSquared < 9.9975f;
0214     } else if (layer1 == 1 and layer2 == 7 and layer3 == 8) {
0215       return rzChiSquared < 8.6369f;
0216     } else if (layer1 == 2 and layer2 == 3 and layer3 == 7) {
0217       return rzChiSquared < 37.945f;
0218     } else if (layer1 == 2 and layer2 == 3 and layer3 == 12) {
0219       return rzChiSquared < 43.0167f;
0220     } else if (layer1 == 2 and layer2 == 3 and layer3 == 4) {
0221       return rzChiSquared < 8.6923f;
0222     } else if (layer1 == 2 and layer2 == 7 and layer3 == 8) {
0223       return rzChiSquared < 11.9672f;
0224     } else if (layer1 == 2 and layer2 == 7 and layer3 == 13) {
0225       return rzChiSquared < 16.2133f;
0226     }
0227 
0228     //default - category not found!
0229     return true;
0230   }
0231 
0232   //TODO: merge this one and the pT5 function later into a single function
0233   template <typename TAcc>
0234   ALPAKA_FN_ACC ALPAKA_FN_INLINE float computePT3RPhiChiSquared(TAcc const& acc,
0235                                                                 ModulesConst modules,
0236                                                                 uint16_t* lowerModuleIndices,
0237                                                                 float g,
0238                                                                 float f,
0239                                                                 float radius,
0240                                                                 float* xs,
0241                                                                 float* ys) {
0242     float delta1[3]{}, delta2[3]{}, slopes[3]{};
0243     bool isFlat[3]{};
0244     float chiSquared = 0;
0245     float inv1 = kWidthPS / kWidth2S;
0246     float inv2 = kPixelPSZpitch / kWidth2S;
0247     for (size_t i = 0; i < 3; i++) {
0248       ModuleType moduleType = modules.moduleType()[lowerModuleIndices[i]];
0249       short moduleSubdet = modules.subdets()[lowerModuleIndices[i]];
0250       short moduleSide = modules.sides()[lowerModuleIndices[i]];
0251       float drdz = modules.drdzs()[lowerModuleIndices[i]];
0252       slopes[i] = modules.dxdys()[lowerModuleIndices[i]];
0253       //category 1 - barrel PS flat
0254       if (moduleSubdet == Barrel and moduleType == PS and moduleSide == Center) {
0255         delta1[i] = inv1;
0256         delta2[i] = inv1;
0257         slopes[i] = -999;
0258         isFlat[i] = true;
0259       }
0260       //category 2 - barrel 2S
0261       else if (moduleSubdet == Barrel and moduleType == TwoS) {
0262         delta1[i] = 1;
0263         delta2[i] = 1;
0264         slopes[i] = -999;
0265         isFlat[i] = true;
0266       }
0267       //category 3 - barrel PS tilted
0268       else if (moduleSubdet == Barrel and moduleType == PS and moduleSide != Center) {
0269         delta1[i] = inv1;
0270         isFlat[i] = false;
0271         delta2[i] = (inv2 * drdz / alpaka::math::sqrt(acc, 1 + drdz * drdz));
0272       }
0273       //category 4 - endcap PS
0274       else if (moduleSubdet == Endcap and moduleType == PS) {
0275         delta1[i] = inv1;
0276         isFlat[i] = false;
0277 
0278         /*
0279         despite the type of the module layer of the lower module index, all anchor
0280         hits are on the pixel side and all non-anchor hits are on the strip side!
0281         */
0282         delta2[i] = inv2;
0283       }
0284       //category 5 - endcap 2S
0285       else if (moduleSubdet == Endcap and moduleType == TwoS) {
0286         delta1[i] = 1;
0287         delta2[i] = 500 * inv1;
0288         isFlat[i] = false;
0289       }
0290 #ifdef WARNINGS
0291       else {
0292         printf("ERROR!!!!! I SHOULDN'T BE HERE!!!! subdet = %d, type = %d, side = %d\n",
0293                moduleSubdet,
0294                moduleType,
0295                moduleSide);
0296       }
0297 #endif
0298     }
0299     chiSquared = computeChiSquared(acc, 3, xs, ys, delta1, delta2, slopes, isFlat, g, f, radius);
0300 
0301     return chiSquared;
0302   }
0303 
0304   ALPAKA_FN_ACC ALPAKA_FN_INLINE float computePT3RPhiChiSquaredInwards(
0305       float g, float f, float r, float* xPix, float* yPix) {
0306     float residual = (xPix[0] - g) * (xPix[0] - g) + (yPix[0] - f) * (yPix[0] - f) - r * r;
0307     float chiSquared = residual * residual;
0308     residual = (xPix[1] - g) * (xPix[1] - g) + (yPix[1] - f) * (yPix[1] - f) - r * r;
0309     chiSquared += residual * residual;
0310 
0311     chiSquared *= 0.5f;
0312     return chiSquared;
0313   }
0314 
0315   //90pc threshold
0316   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passPT3RPhiChiSquaredCuts(ModulesConst modules,
0317                                                                 uint16_t lowerModuleIndex1,
0318                                                                 uint16_t lowerModuleIndex2,
0319                                                                 uint16_t lowerModuleIndex3,
0320                                                                 float chiSquared) {
0321     const int layer1 =
0322         modules.layers()[lowerModuleIndex1] + 6 * (modules.subdets()[lowerModuleIndex1] == Endcap) +
0323         5 * (modules.subdets()[lowerModuleIndex1] == Endcap and modules.moduleType()[lowerModuleIndex1] == TwoS);
0324     const int layer2 =
0325         modules.layers()[lowerModuleIndex2] + 6 * (modules.subdets()[lowerModuleIndex2] == Endcap) +
0326         5 * (modules.subdets()[lowerModuleIndex2] == Endcap and modules.moduleType()[lowerModuleIndex2] == TwoS);
0327     const int layer3 =
0328         modules.layers()[lowerModuleIndex3] + 6 * (modules.subdets()[lowerModuleIndex3] == Endcap) +
0329         5 * (modules.subdets()[lowerModuleIndex3] == Endcap and modules.moduleType()[lowerModuleIndex3] == TwoS);
0330 
0331     if (layer1 == 8 and layer2 == 9 and layer3 == 10) {
0332       return chiSquared < 7.003f;
0333     } else if (layer1 == 8 and layer2 == 9 and layer3 == 15) {
0334       return chiSquared < 0.5f;
0335     } else if (layer1 == 7 and layer2 == 8 and layer3 == 9) {
0336       return chiSquared < 8.046f;
0337     } else if (layer1 == 7 and layer2 == 8 and layer3 == 14) {
0338       return chiSquared < 0.575f;
0339     } else if (layer1 == 1 and layer2 == 2 and layer3 == 7) {
0340       return chiSquared < 5.304f;
0341     } else if (layer1 == 1 and layer2 == 2 and layer3 == 3) {
0342       return chiSquared < 10.6211f;
0343     } else if (layer1 == 1 and layer2 == 7 and layer3 == 8) {
0344       return chiSquared < 4.617f;
0345     } else if (layer1 == 2 and layer2 == 7 and layer3 == 8) {
0346       return chiSquared < 8.046f;
0347     } else if (layer1 == 2 and layer2 == 7 and layer3 == 13) {
0348       return chiSquared < 0.435f;
0349     } else if (layer1 == 2 and layer2 == 3 and layer3 == 7) {
0350       return chiSquared < 9.244f;
0351     } else if (layer1 == 2 and layer2 == 3 and layer3 == 12) {
0352       return chiSquared < 0.287f;
0353     } else if (layer1 == 2 and layer2 == 3 and layer3 == 4) {
0354       return chiSquared < 18.509f;
0355     }
0356 
0357     return true;
0358   }
0359 
0360   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passPT3RPhiChiSquaredInwardsCuts(ModulesConst modules,
0361                                                                        uint16_t lowerModuleIndex1,
0362                                                                        uint16_t lowerModuleIndex2,
0363                                                                        uint16_t lowerModuleIndex3,
0364                                                                        float chiSquared) {
0365     const int layer1 =
0366         modules.layers()[lowerModuleIndex1] + 6 * (modules.subdets()[lowerModuleIndex1] == Endcap) +
0367         5 * (modules.subdets()[lowerModuleIndex1] == Endcap and modules.moduleType()[lowerModuleIndex1] == TwoS);
0368     const int layer2 =
0369         modules.layers()[lowerModuleIndex2] + 6 * (modules.subdets()[lowerModuleIndex2] == Endcap) +
0370         5 * (modules.subdets()[lowerModuleIndex2] == Endcap and modules.moduleType()[lowerModuleIndex2] == TwoS);
0371     const int layer3 =
0372         modules.layers()[lowerModuleIndex3] + 6 * (modules.subdets()[lowerModuleIndex3] == Endcap) +
0373         5 * (modules.subdets()[lowerModuleIndex3] == Endcap and modules.moduleType()[lowerModuleIndex3] == TwoS);
0374 
0375     if (layer1 == 7 and layer2 == 8 and layer3 == 9)  // endcap layer 1,2,3, ps
0376     {
0377       return chiSquared < 22016.8055f;
0378     } else if (layer1 == 7 and layer2 == 8 and layer3 == 14)  // endcap layer 1,2,3 layer3->2s
0379     {
0380       return chiSquared < 935179.56807f;
0381     } else if (layer1 == 8 and layer2 == 9 and layer3 == 10)  // endcap layer 2,3,4
0382     {
0383       return chiSquared < 29064.12959f;
0384     } else if (layer1 == 8 and layer2 == 9 and layer3 == 15)  // endcap layer 2,3,4, layer3->2s
0385     {
0386       return chiSquared < 935179.5681f;
0387     } else if (layer1 == 1 and layer2 == 2 and layer3 == 3)  // barrel 1,2,3
0388     {
0389       return chiSquared < 1370.0113195101474f;
0390     } else if (layer1 == 1 and layer2 == 2 and layer3 == 7)  // barrel 1,2 endcap 1
0391     {
0392       return chiSquared < 5492.110048314815f;
0393     } else if (layer1 == 2 and layer2 == 3 and layer3 == 4)  // barrel 2,3,4
0394     {
0395       return chiSquared < 4160.410806470067f;
0396     } else if (layer1 == 1 and layer2 == 7 and layer3 == 8)  // barrel 1, endcap 1,2
0397     {
0398       return chiSquared < 29064.129591225726f;
0399     } else if (layer1 == 2 and layer2 == 3 and layer3 == 7)  // barrel 2,3 endcap 1
0400     {
0401       return chiSquared < 12634.215376250893f;
0402     } else if (layer1 == 2 and layer2 == 3 and layer3 == 12)  // barrel 2,3, endcap 1->2s
0403     {
0404       return chiSquared < 353821.69361145404f;
0405     } else if (layer1 == 2 and layer2 == 7 and layer3 == 8)  // barrel2, endcap 1,2
0406     {
0407       return chiSquared < 33393.26076341235f;
0408     } else if (layer1 == 2 and layer2 == 7 and layer3 == 13)  //barrel 2, endcap 1, endcap2->2s
0409     {
0410       return chiSquared < 935179.5680742573f;
0411     }
0412 
0413     return true;
0414   }
0415 
0416   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool checkIntervalOverlappT3(float firstMin,
0417                                                               float firstMax,
0418                                                               float secondMin,
0419                                                               float secondMax) {
0420     return ((firstMin <= secondMin) && (secondMin < firstMax)) || ((secondMin < firstMin) && (firstMin < secondMax));
0421   }
0422 
0423   /*bounds for high Pt taken from : http://uaf-10.t2.ucsd.edu/~bsathian/SDL/T5_efficiency/efficiencies/new_efficiencies/efficiencies_20210513_T5_recovering_high_Pt_efficiencies/highE_radius_matching/highE_bounds.txt */
0424   template <typename TAcc>
0425   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passRadiusCriterionBBB(TAcc const& acc,
0426                                                              float pixelRadius,
0427                                                              float pixelRadiusError,
0428                                                              float tripletRadius) {
0429     float tripletInvRadiusErrorBound = 0.15624f;
0430     float pixelInvRadiusErrorBound = 0.17235f;
0431 
0432     if (pixelRadius > 2.0f * kR1GeVf) {
0433       pixelInvRadiusErrorBound = 0.6375f;
0434       tripletInvRadiusErrorBound = 0.6588f;
0435     }
0436 
0437     float tripletRadiusInvMax = (1 + tripletInvRadiusErrorBound) / tripletRadius;
0438     float tripletRadiusInvMin = alpaka::math::max(acc, (1 - tripletInvRadiusErrorBound) / tripletRadius, 0.0f);
0439 
0440     float pixelRadiusInvMax =
0441         alpaka::math::max(acc, (1 + pixelInvRadiusErrorBound) / pixelRadius, 1.f / (pixelRadius - pixelRadiusError));
0442     float pixelRadiusInvMin =
0443         alpaka::math::min(acc, (1 - pixelInvRadiusErrorBound) / pixelRadius, 1.f / (pixelRadius + pixelRadiusError));
0444 
0445     return checkIntervalOverlappT3(tripletRadiusInvMin, tripletRadiusInvMax, pixelRadiusInvMin, pixelRadiusInvMax);
0446   }
0447 
0448   template <typename TAcc>
0449   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passRadiusCriterionBBE(TAcc const& acc,
0450                                                              float pixelRadius,
0451                                                              float pixelRadiusError,
0452                                                              float tripletRadius) {
0453     float tripletInvRadiusErrorBound = 0.45972f;
0454     float pixelInvRadiusErrorBound = 0.19644f;
0455 
0456     if (pixelRadius > 2.0f * kR1GeVf) {
0457       pixelInvRadiusErrorBound = 0.6805f;
0458       tripletInvRadiusErrorBound = 0.8557f;
0459     }
0460 
0461     float tripletRadiusInvMax = (1 + tripletInvRadiusErrorBound) / tripletRadius;
0462     float tripletRadiusInvMin = alpaka::math::max(acc, (1 - tripletInvRadiusErrorBound) / tripletRadius, 0.0f);
0463 
0464     float pixelRadiusInvMax =
0465         alpaka::math::max(acc, (1 + pixelInvRadiusErrorBound) / pixelRadius, 1.f / (pixelRadius - pixelRadiusError));
0466     float pixelRadiusInvMin =
0467         alpaka::math::min(acc, (1 - pixelInvRadiusErrorBound) / pixelRadius, 1.f / (pixelRadius + pixelRadiusError));
0468 
0469     return checkIntervalOverlappT3(tripletRadiusInvMin, tripletRadiusInvMax, pixelRadiusInvMin, pixelRadiusInvMax);
0470   }
0471 
0472   template <typename TAcc>
0473   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passRadiusCriterionBEE(TAcc const& acc,
0474                                                              float pixelRadius,
0475                                                              float pixelRadiusError,
0476                                                              float tripletRadius) {
0477     float tripletInvRadiusErrorBound = 1.59294f;
0478     float pixelInvRadiusErrorBound = 0.255181f;
0479 
0480     if (pixelRadius > 2.0f * kR1GeVf)  //as good as not having selections
0481     {
0482       pixelInvRadiusErrorBound = 2.2091f;
0483       tripletInvRadiusErrorBound = 2.3548f;
0484     }
0485 
0486     float tripletRadiusInvMax = (1 + tripletInvRadiusErrorBound) / tripletRadius;
0487     float tripletRadiusInvMin = alpaka::math::max(acc, (1 - tripletInvRadiusErrorBound) / tripletRadius, 0.0f);
0488 
0489     float pixelRadiusInvMax =
0490         alpaka::math::max(acc, (1 + pixelInvRadiusErrorBound) / pixelRadius, 1.f / (pixelRadius - pixelRadiusError));
0491     float pixelRadiusInvMin =
0492         alpaka::math::min(acc, (1 - pixelInvRadiusErrorBound) / pixelRadius, 1.f / (pixelRadius + pixelRadiusError));
0493     pixelRadiusInvMin = alpaka::math::max(acc, pixelRadiusInvMin, 0.0f);
0494 
0495     return checkIntervalOverlappT3(tripletRadiusInvMin, tripletRadiusInvMax, pixelRadiusInvMin, pixelRadiusInvMax);
0496   }
0497 
0498   template <typename TAcc>
0499   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passRadiusCriterionEEE(TAcc const& acc,
0500                                                              float pixelRadius,
0501                                                              float pixelRadiusError,
0502                                                              float tripletRadius) {
0503     float tripletInvRadiusErrorBound = 1.7006f;
0504     float pixelInvRadiusErrorBound = 0.26367f;
0505 
0506     if (pixelRadius > 2.0f * kR1GeVf)  //as good as not having selections
0507     {
0508       pixelInvRadiusErrorBound = 2.286f;
0509       tripletInvRadiusErrorBound = 2.436f;
0510     }
0511 
0512     float tripletRadiusInvMax = (1 + tripletInvRadiusErrorBound) / tripletRadius;
0513     float tripletRadiusInvMin = alpaka::math::max(acc, (1 - tripletInvRadiusErrorBound) / tripletRadius, 0.0f);
0514 
0515     float pixelRadiusInvMax =
0516         alpaka::math::max(acc, (1 + pixelInvRadiusErrorBound) / pixelRadius, 1.f / (pixelRadius - pixelRadiusError));
0517     float pixelRadiusInvMin =
0518         alpaka::math::min(acc, (1 - pixelInvRadiusErrorBound) / pixelRadius, 1.f / (pixelRadius + pixelRadiusError));
0519     pixelRadiusInvMin = alpaka::math::max(acc, 0.0f, pixelRadiusInvMin);
0520 
0521     return checkIntervalOverlappT3(tripletRadiusInvMin, tripletRadiusInvMax, pixelRadiusInvMin, pixelRadiusInvMax);
0522   }
0523 
0524   template <typename TAcc>
0525   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passRadiusCriterion(TAcc const& acc,
0526                                                           ModulesConst modules,
0527                                                           float pixelRadius,
0528                                                           float pixelRadiusError,
0529                                                           float tripletRadius,
0530                                                           int16_t lowerModuleIndex,
0531                                                           uint16_t middleModuleIndex,
0532                                                           uint16_t upperModuleIndex) {
0533     if (modules.subdets()[lowerModuleIndex] == Endcap) {
0534       return passRadiusCriterionEEE(acc, pixelRadius, pixelRadiusError, tripletRadius);
0535     } else if (modules.subdets()[middleModuleIndex] == Endcap) {
0536       return passRadiusCriterionBEE(acc, pixelRadius, pixelRadiusError, tripletRadius);
0537     } else if (modules.subdets()[upperModuleIndex] == Endcap) {
0538       return passRadiusCriterionBBE(acc, pixelRadius, pixelRadiusError, tripletRadius);
0539     } else {
0540       return passRadiusCriterionBBB(acc, pixelRadius, pixelRadiusError, tripletRadius);
0541     }
0542   }
0543 
0544   template <typename TAcc>
0545   ALPAKA_FN_ACC ALPAKA_FN_INLINE float computePT3RZChiSquared(TAcc const& acc,
0546                                                               ModulesConst modules,
0547                                                               const uint16_t* lowerModuleIndices,
0548                                                               const float* rtPix,
0549                                                               const float* xPix,
0550                                                               const float* yPix,
0551                                                               const float* zPix,
0552                                                               const float* rts,
0553                                                               const float* xs,
0554                                                               const float* ys,
0555                                                               const float* zs,
0556                                                               float pixelSegmentPt,
0557                                                               float pixelSegmentPx,
0558                                                               float pixelSegmentPy,
0559                                                               float pixelSegmentPz,
0560                                                               int pixelSegmentCharge) {
0561     float residual = 0;
0562     float error2 = 0;
0563     float RMSE = 0;
0564 
0565     float Px = pixelSegmentPx, Py = pixelSegmentPy, Pz = pixelSegmentPz;
0566     int charge = pixelSegmentCharge;
0567     float x1 = xPix[1] / 100;
0568     float y1 = yPix[1] / 100;
0569     float z1 = zPix[1] / 100;
0570     float r1 = rtPix[1] / 100;
0571 
0572     float a = -2.f * k2Rinv1GeVf * 100 * charge;  // multiply by 100 to make the correct length units
0573 
0574     for (size_t i = 0; i < Params_T3::kLayers; i++) {
0575       float zsi = zs[i] / 100;
0576       float rtsi = rts[i] / 100;
0577       uint16_t lowerModuleIndex = lowerModuleIndices[i];
0578       const int moduleType = modules.moduleType()[lowerModuleIndex];
0579       const int moduleSide = modules.sides()[lowerModuleIndex];
0580       const int moduleSubdet = modules.subdets()[lowerModuleIndex];
0581 
0582       // calculation is detailed documented here https://indico.cern.ch/event/1185895/contributions/4982756/attachments/2526561/4345805/helix%20pT3%20summarize.pdf
0583       float diffr, diffz;
0584       float p = alpaka::math::sqrt(acc, Px * Px + Py * Py + Pz * Pz);
0585 
0586       float rou = a / p;
0587       if (moduleSubdet == Endcap) {
0588         float s = (zsi - z1) * p / Pz;
0589         float x = x1 + Px / a * alpaka::math::sin(acc, rou * s) - Py / a * (1 - alpaka::math::cos(acc, rou * s));
0590         float y = y1 + Py / a * alpaka::math::sin(acc, rou * s) + Px / a * (1 - alpaka::math::cos(acc, rou * s));
0591         diffr = alpaka::math::abs(acc, rtsi - alpaka::math::sqrt(acc, x * x + y * y)) * 100;
0592       }
0593 
0594       if (moduleSubdet == Barrel) {
0595         float paraA = r1 * r1 + 2 * (Px * Px + Py * Py) / (a * a) + 2 * (y1 * Px - x1 * Py) / a - rtsi * rtsi;
0596         float paraB = 2 * (x1 * Px + y1 * Py) / a;
0597         float paraC = 2 * (y1 * Px - x1 * Py) / a + 2 * (Px * Px + Py * Py) / (a * a);
0598         float A = paraB * paraB + paraC * paraC;
0599         float B = 2 * paraA * paraB;
0600         float C = paraA * paraA - paraC * paraC;
0601         float sol1 = (-B + alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
0602         float sol2 = (-B - alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
0603         float solz1 = alpaka::math::asin(acc, sol1) / rou * Pz / p + z1;
0604         float solz2 = alpaka::math::asin(acc, sol2) / rou * Pz / p + z1;
0605         float diffz1 = alpaka::math::abs(acc, solz1 - zsi) * 100;
0606         float diffz2 = alpaka::math::abs(acc, solz2 - zsi) * 100;
0607         diffz = alpaka::math::min(acc, diffz1, diffz2);
0608       }
0609 
0610       residual = moduleSubdet == Barrel ? diffz : diffr;
0611 
0612       //PS Modules
0613       if (moduleType == 0) {
0614         error2 = kPixelPSZpitch * kPixelPSZpitch;
0615       } else  //2S modules
0616       {
0617         error2 = kStrip2SZpitch * kStrip2SZpitch;
0618       }
0619 
0620       //special dispensation to tilted PS modules!
0621       if (moduleType == 0 and moduleSubdet == Barrel and moduleSide != Center) {
0622         float drdz = modules.drdzs()[lowerModuleIndex];
0623         error2 /= (1 + drdz * drdz);
0624       }
0625       RMSE += (residual * residual) / error2;
0626     }
0627 
0628     RMSE = alpaka::math::sqrt(acc, 0.2f * RMSE);  // Divided by the degree of freedom 5.
0629 
0630     return RMSE;
0631   }
0632 
0633   template <typename TAcc>
0634   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runPixelTripletDefaultAlgo(TAcc const& acc,
0635                                                                  ModulesConst modules,
0636                                                                  ObjectRangesConst ranges,
0637                                                                  MiniDoubletsConst mds,
0638                                                                  SegmentsConst segments,
0639                                                                  PixelSeedsConst pixelSeeds,
0640                                                                  PixelSegmentsConst pixelSegments,
0641                                                                  TripletsConst triplets,
0642                                                                  unsigned int pixelSegmentIndex,
0643                                                                  unsigned int tripletIndex,
0644                                                                  float& pixelRadius,
0645                                                                  float& tripletRadius,
0646                                                                  float& centerX,
0647                                                                  float& centerY,
0648                                                                  float& rzChiSquared,
0649                                                                  float& rPhiChiSquared,
0650                                                                  float& rPhiChiSquaredInwards,
0651                                                                  float& pixelRadiusError,
0652                                                                  const float ptCut,
0653                                                                  bool runDNN = true,
0654                                                                  bool runChiSquaredCuts = true) {
0655     //run pT4 compatibility between the pixel segment and inner segment, and between the pixel and outer segment of the triplet
0656     uint16_t pixelModuleIndex = segments.innerLowerModuleIndices()[pixelSegmentIndex];
0657 
0658     uint16_t lowerModuleIndex = triplets.lowerModuleIndices()[tripletIndex][0];
0659     uint16_t middleModuleIndex = triplets.lowerModuleIndices()[tripletIndex][1];
0660     uint16_t upperModuleIndex = triplets.lowerModuleIndices()[tripletIndex][2];
0661 
0662     {
0663       // pixel segment vs inner segment of the triplet
0664       if (not runPixelTrackletDefaultAlgopT3(acc,
0665                                              modules,
0666                                              ranges,
0667                                              mds,
0668                                              segments,
0669                                              pixelSeeds,
0670                                              pixelModuleIndex,
0671                                              lowerModuleIndex,
0672                                              middleModuleIndex,
0673                                              pixelSegmentIndex,
0674                                              triplets.segmentIndices()[tripletIndex][0],
0675                                              ptCut))
0676         return false;
0677 
0678       //pixel segment vs outer segment of triplet
0679       if (not runPixelTrackletDefaultAlgopT3(acc,
0680                                              modules,
0681                                              ranges,
0682                                              mds,
0683                                              segments,
0684                                              pixelSeeds,
0685                                              pixelModuleIndex,
0686                                              middleModuleIndex,
0687                                              upperModuleIndex,
0688                                              pixelSegmentIndex,
0689                                              triplets.segmentIndices()[tripletIndex][1],
0690                                              ptCut))
0691         return false;
0692     }
0693 
0694     //pt matching between the pixel ptin and the triplet circle pt
0695     unsigned int pixelSegmentArrayIndex = pixelSegmentIndex - ranges.segmentModuleIndices()[pixelModuleIndex];
0696     float pixelSegmentPt = pixelSeeds.ptIn()[pixelSegmentArrayIndex];
0697     float pixelSegmentPtError = pixelSeeds.ptErr()[pixelSegmentArrayIndex];
0698     float pixelSegmentPx = pixelSeeds.px()[pixelSegmentArrayIndex];
0699     float pixelSegmentPy = pixelSeeds.py()[pixelSegmentArrayIndex];
0700     float pixelSegmentPz = pixelSeeds.pz()[pixelSegmentArrayIndex];
0701     int pixelSegmentCharge = pixelSeeds.charge()[pixelSegmentArrayIndex];
0702 
0703     float pixelG = pixelSegments.circleCenterX()[pixelSegmentArrayIndex];
0704     float pixelF = pixelSegments.circleCenterY()[pixelSegmentArrayIndex];
0705     float pixelRadiusPCA = pixelSegments.circleRadius()[pixelSegmentArrayIndex];
0706 
0707     unsigned int pixelInnerMDIndex = segments.mdIndices()[pixelSegmentIndex][0];
0708     unsigned int pixelOuterMDIndex = segments.mdIndices()[pixelSegmentIndex][1];
0709 
0710     pixelRadius = pixelSegmentPt * kR1GeVf;
0711     pixelRadiusError = pixelSegmentPtError * kR1GeVf;
0712     unsigned int tripletInnerSegmentIndex = triplets.segmentIndices()[tripletIndex][0];
0713     unsigned int tripletOuterSegmentIndex = triplets.segmentIndices()[tripletIndex][1];
0714 
0715     unsigned int firstMDIndex = segments.mdIndices()[tripletInnerSegmentIndex][0];
0716     unsigned int secondMDIndex = segments.mdIndices()[tripletInnerSegmentIndex][1];
0717     unsigned int thirdMDIndex = segments.mdIndices()[tripletOuterSegmentIndex][1];
0718 
0719     float xs[Params_T3::kLayers] = {
0720         mds.anchorX()[firstMDIndex], mds.anchorX()[secondMDIndex], mds.anchorX()[thirdMDIndex]};
0721     float ys[Params_T3::kLayers] = {
0722         mds.anchorY()[firstMDIndex], mds.anchorY()[secondMDIndex], mds.anchorY()[thirdMDIndex]};
0723 
0724     float g, f;
0725     tripletRadius = triplets.radius()[tripletIndex];
0726     g = triplets.centerX()[tripletIndex];
0727     f = triplets.centerY()[tripletIndex];
0728 
0729     if (not passRadiusCriterion(acc,
0730                                 modules,
0731                                 pixelRadius,
0732                                 pixelRadiusError,
0733                                 tripletRadius,
0734                                 lowerModuleIndex,
0735                                 middleModuleIndex,
0736                                 upperModuleIndex))
0737       return false;
0738 
0739     uint16_t lowerModuleIndices[Params_T3::kLayers] = {lowerModuleIndex, middleModuleIndex, upperModuleIndex};
0740 
0741     if (runDNN || runChiSquaredCuts) {
0742       float rts[Params_T3::kLayers] = {
0743           mds.anchorRt()[firstMDIndex], mds.anchorRt()[secondMDIndex], mds.anchorRt()[thirdMDIndex]};
0744       float zs[Params_T3::kLayers] = {
0745           mds.anchorZ()[firstMDIndex], mds.anchorZ()[secondMDIndex], mds.anchorZ()[thirdMDIndex]};
0746       float rtPix[Params_pLS::kLayers] = {mds.anchorRt()[pixelInnerMDIndex], mds.anchorRt()[pixelOuterMDIndex]};
0747       float xPix[Params_pLS::kLayers] = {mds.anchorX()[pixelInnerMDIndex], mds.anchorX()[pixelOuterMDIndex]};
0748       float yPix[Params_pLS::kLayers] = {mds.anchorY()[pixelInnerMDIndex], mds.anchorY()[pixelOuterMDIndex]};
0749       float zPix[Params_pLS::kLayers] = {mds.anchorZ()[pixelInnerMDIndex], mds.anchorZ()[pixelOuterMDIndex]};
0750 
0751       rzChiSquared = computePT3RZChiSquared(acc,
0752                                             modules,
0753                                             lowerModuleIndices,
0754                                             rtPix,
0755                                             xPix,
0756                                             yPix,
0757                                             zPix,
0758                                             rts,
0759                                             xs,
0760                                             ys,
0761                                             zs,
0762                                             pixelSegmentPt,
0763                                             pixelSegmentPx,
0764                                             pixelSegmentPy,
0765                                             pixelSegmentPz,
0766                                             pixelSegmentCharge);
0767       if (runChiSquaredCuts && pixelSegmentPt < 5.0f) {
0768         if (!passPT3RZChiSquaredCuts(modules, lowerModuleIndex, middleModuleIndex, upperModuleIndex, rzChiSquared))
0769           return false;
0770       }
0771 
0772       rPhiChiSquared =
0773           computePT3RPhiChiSquared(acc, modules, lowerModuleIndices, pixelG, pixelF, pixelRadiusPCA, xs, ys);
0774       if (runChiSquaredCuts && pixelSegmentPt < 5.0f) {
0775         if (!passPT3RPhiChiSquaredCuts(modules, lowerModuleIndex, middleModuleIndex, upperModuleIndex, rPhiChiSquared))
0776           return false;
0777       }
0778 
0779       rPhiChiSquaredInwards = computePT3RPhiChiSquaredInwards(g, f, tripletRadius, xPix, yPix);
0780       if (runChiSquaredCuts && pixelSegmentPt < 5.0f) {
0781         if (!passPT3RPhiChiSquaredInwardsCuts(
0782                 modules, lowerModuleIndex, middleModuleIndex, upperModuleIndex, rPhiChiSquaredInwards))
0783           return false;
0784       }
0785     }
0786 
0787     centerX = 0;
0788     centerY = 0;
0789 
0790     if (runDNN and !lst::pt3dnn::runInference(acc,
0791                                               rPhiChiSquared,
0792                                               tripletRadius,
0793                                               pixelRadius,
0794                                               pixelRadiusError,
0795                                               rzChiSquared,
0796                                               pixelSeeds.eta()[pixelSegmentArrayIndex],
0797                                               pixelSegmentPt)) {
0798       return false;
0799     }
0800 
0801     return true;
0802   }
0803 
0804   struct CreatePixelTripletsFromMap {
0805     ALPAKA_FN_ACC void operator()(Acc3D const& acc,
0806                                   ModulesConst modules,
0807                                   ModulesPixelConst modulesPixel,
0808                                   ObjectRangesConst ranges,
0809                                   MiniDoubletsConst mds,
0810                                   SegmentsConst segments,
0811                                   PixelSeedsConst pixelSeeds,
0812                                   PixelSegmentsConst pixelSegments,
0813                                   Triplets triplets,
0814                                   TripletsOccupancyConst tripletsOccupancy,
0815                                   PixelTriplets pixelTriplets,
0816                                   unsigned int* connectedPixelSize,
0817                                   unsigned int* connectedPixelIndex,
0818                                   unsigned int nPixelSegments,
0819                                   const float ptCut) const {
0820       for (unsigned int i_pLS : cms::alpakatools::uniform_elements_z(acc, nPixelSegments)) {
0821         auto iLSModule_max = connectedPixelIndex[i_pLS] + connectedPixelSize[i_pLS];
0822 
0823         for (unsigned int iLSModule :
0824              cms::alpakatools::uniform_elements_y(acc, connectedPixelIndex[i_pLS], iLSModule_max)) {
0825           uint16_t tripletLowerModuleIndex =
0826               modulesPixel.connectedPixels()
0827                   [iLSModule];  //connected pixels will have the appropriate lower module index by default!
0828 #ifdef WARNINGS
0829           if (tripletLowerModuleIndex >= modules.nLowerModules()) {
0830             printf("tripletLowerModuleIndex %d >= modules.nLowerModules %d \n",
0831                    tripletLowerModuleIndex,
0832                    modules.nLowerModules());
0833             continue;  //sanity check
0834           }
0835 #endif
0836           //Removes 2S-2S :FIXME: filter these out in the pixel map
0837           if (modules.moduleType()[tripletLowerModuleIndex] == TwoS)
0838             continue;
0839 
0840           uint16_t pixelModuleIndex = modules.nLowerModules();
0841           unsigned int nOuterTriplets = tripletsOccupancy.nTriplets()[tripletLowerModuleIndex];
0842           if (nOuterTriplets == 0)
0843             continue;
0844 
0845           unsigned int pixelSegmentIndex = ranges.segmentModuleIndices()[pixelModuleIndex] + i_pLS;
0846 
0847           if (pixelSegments.isDup()[i_pLS])
0848             continue;
0849           if (pixelSegments.partOfPT5()[i_pLS])
0850             continue;  //don't make pT3s for those pixels that are part of pT5
0851 
0852           short layer2_adjustment;
0853           if (modules.layers()[tripletLowerModuleIndex] == 1) {
0854             layer2_adjustment = 1;
0855           }  //get upper segment to be in second layer
0856           else if (modules.layers()[tripletLowerModuleIndex] == 2) {
0857             layer2_adjustment = 0;
0858           }  // get lower segment to be in second layer
0859           else {
0860             continue;
0861           }
0862 
0863           //fetch the triplet
0864           for (unsigned int outerTripletArrayIndex : cms::alpakatools::uniform_elements_x(acc, nOuterTriplets)) {
0865             unsigned int outerTripletIndex =
0866                 ranges.tripletModuleIndices()[tripletLowerModuleIndex] + outerTripletArrayIndex;
0867             if (modules.moduleType()[triplets.lowerModuleIndices()[outerTripletIndex][1]] == TwoS)
0868               continue;  //REMOVES PS-2S
0869 
0870             if (triplets.partOfPT5()[outerTripletIndex])
0871               continue;  //don't create pT3s for T3s accounted in pT5s
0872 
0873             float pixelRadius, tripletRadius, rPhiChiSquared, rzChiSquared, rPhiChiSquaredInwards, centerX, centerY,
0874                 pixelRadiusError;
0875             bool success = runPixelTripletDefaultAlgo(acc,
0876                                                       modules,
0877                                                       ranges,
0878                                                       mds,
0879                                                       segments,
0880                                                       pixelSeeds,
0881                                                       pixelSegments,
0882                                                       triplets,
0883                                                       pixelSegmentIndex,
0884                                                       outerTripletIndex,
0885                                                       pixelRadius,
0886                                                       tripletRadius,
0887                                                       centerX,
0888                                                       centerY,
0889                                                       rzChiSquared,
0890                                                       rPhiChiSquared,
0891                                                       rPhiChiSquaredInwards,
0892                                                       pixelRadiusError,
0893                                                       ptCut);
0894 
0895             if (success) {
0896               float phi =
0897                   mds.anchorPhi()[segments
0898                                       .mdIndices()[triplets.segmentIndices()[outerTripletIndex][0]][layer2_adjustment]];
0899               float eta =
0900                   mds.anchorEta()[segments
0901                                       .mdIndices()[triplets.segmentIndices()[outerTripletIndex][0]][layer2_adjustment]];
0902               float eta_pix = pixelSeeds.eta()[i_pLS];
0903               float phi_pix = pixelSeeds.phi()[i_pLS];
0904               float pt = pixelSeeds.ptIn()[i_pLS];
0905               float score = rPhiChiSquared + rPhiChiSquaredInwards;
0906               unsigned int totOccupancyPixelTriplets =
0907                   alpaka::atomicAdd(acc, &pixelTriplets.totOccupancyPixelTriplets(), 1u, alpaka::hierarchy::Threads{});
0908               if (totOccupancyPixelTriplets >= n_max_pixel_triplets) {
0909 #ifdef WARNINGS
0910                 printf("Pixel Triplet excess alert!\n");
0911 #endif
0912               } else {
0913                 unsigned int pixelTripletIndex =
0914                     alpaka::atomicAdd(acc, &pixelTriplets.nPixelTriplets(), 1u, alpaka::hierarchy::Threads{});
0915                 addPixelTripletToMemory(mds,
0916                                         segments,
0917                                         triplets,
0918                                         pixelTriplets,
0919                                         pixelSegmentIndex,
0920                                         outerTripletIndex,
0921                                         pixelRadius,
0922                                         tripletRadius,
0923                                         centerX,
0924                                         centerY,
0925                                         rPhiChiSquared,
0926                                         rPhiChiSquaredInwards,
0927                                         rzChiSquared,
0928                                         pixelTripletIndex,
0929                                         pt,
0930                                         eta,
0931                                         phi,
0932                                         eta_pix,
0933                                         phi_pix,
0934                                         pixelRadiusError,
0935                                         score);
0936                 triplets.partOfPT3()[outerTripletIndex] = true;
0937               }
0938             }
0939           }  // for outerTripletArrayIndex
0940         }  // for iLSModule < iLSModule_max
0941       }  // for i_pLS
0942     }
0943   };
0944 
0945   template <typename TAcc>
0946   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runTripletDefaultAlgoPPBB(TAcc const& acc,
0947                                                                 ModulesConst modules,
0948                                                                 ObjectRangesConst ranges,
0949                                                                 MiniDoubletsConst mds,
0950                                                                 SegmentsConst segments,
0951                                                                 PixelSeedsConst pixelSeeds,
0952                                                                 uint16_t pixelModuleIndex,
0953                                                                 uint16_t outerInnerLowerModuleIndex,
0954                                                                 uint16_t outerOuterLowerModuleIndex,
0955                                                                 unsigned int innerSegmentIndex,
0956                                                                 unsigned int outerSegmentIndex,
0957                                                                 unsigned int firstMDIndex,
0958                                                                 unsigned int secondMDIndex,
0959                                                                 unsigned int thirdMDIndex,
0960                                                                 unsigned int fourthMDIndex,
0961                                                                 const float ptCut) {
0962     float dPhi, betaIn, betaOut, pt_beta, zLo, zHi, zLoPointed, zHiPointed, dPhiCut, betaOutCut;
0963 
0964     bool isPS_OutLo = (modules.moduleType()[outerInnerLowerModuleIndex] == PS);
0965 
0966     float rt_InLo = mds.anchorRt()[firstMDIndex];
0967     float rt_InUp = mds.anchorRt()[secondMDIndex];
0968     float rt_OutLo = mds.anchorRt()[thirdMDIndex];
0969     float rt_OutUp = mds.anchorRt()[fourthMDIndex];
0970 
0971     float z_InUp = mds.anchorZ()[secondMDIndex];
0972     float z_OutLo = mds.anchorZ()[thirdMDIndex];
0973 
0974     float x_InLo = mds.anchorX()[firstMDIndex];
0975     float x_InUp = mds.anchorX()[secondMDIndex];
0976     float x_OutLo = mds.anchorX()[thirdMDIndex];
0977     float x_OutUp = mds.anchorX()[fourthMDIndex];
0978 
0979     float y_InLo = mds.anchorY()[firstMDIndex];
0980     float y_InUp = mds.anchorY()[secondMDIndex];
0981     float y_OutLo = mds.anchorY()[thirdMDIndex];
0982     float y_OutUp = mds.anchorY()[fourthMDIndex];
0983 
0984     float rt_InOut = rt_InUp;
0985 
0986     if (alpaka::math::abs(acc, cms::alpakatools::deltaPhi(acc, x_InUp, y_InUp, x_OutLo, y_OutLo)) > kPi / 2.f)
0987       return false;
0988 
0989     unsigned int pixelSegmentArrayIndex = innerSegmentIndex - ranges.segmentModuleIndices()[pixelModuleIndex];
0990     float ptIn = pixelSeeds.ptIn()[pixelSegmentArrayIndex];
0991     float ptSLo = ptIn;
0992     float px = pixelSeeds.px()[pixelSegmentArrayIndex];
0993     float py = pixelSeeds.py()[pixelSegmentArrayIndex];
0994     float pz = pixelSeeds.pz()[pixelSegmentArrayIndex];
0995     float ptErr = pixelSeeds.ptErr()[pixelSegmentArrayIndex];
0996     float etaErr = pixelSeeds.etaErr()[pixelSegmentArrayIndex];
0997     ptSLo = alpaka::math::max(acc, ptCut, ptSLo - 10.0f * alpaka::math::max(acc, ptErr, 0.005f * ptSLo));
0998     ptSLo = alpaka::math::min(acc, 10.0f, ptSLo);
0999 
1000     float alpha1GeV_OutLo =
1001         alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / ptCut, kSinAlphaMax));
1002     const float rtRatio_OutLoInOut =
1003         rt_OutLo / rt_InOut;  // Outer segment beginning rt divided by inner segment beginning rt;
1004 
1005     float dzDrtScale =
1006         alpaka::math::tan(acc, alpha1GeV_OutLo) / alpha1GeV_OutLo;  // The track can bend in r-z plane slightly
1007     const float zpitch_InLo = 0.05f;
1008     const float zpitch_InOut = 0.05f;
1009     float zpitch_OutLo = (isPS_OutLo ? kPixelPSZpitch : kStrip2SZpitch);
1010     float zGeom = zpitch_InLo + zpitch_OutLo;
1011     zHi = z_InUp + (z_InUp + kDeltaZLum) * (rtRatio_OutLoInOut - 1.f) * (z_InUp < 0.f ? 1.f : dzDrtScale) +
1012           (zpitch_InOut + zpitch_OutLo);
1013     zLo = z_InUp + (z_InUp - kDeltaZLum) * (rtRatio_OutLoInOut - 1.f) * (z_InUp > 0.f ? 1.f : dzDrtScale) -
1014           (zpitch_InOut + zpitch_OutLo);  //slope-correction only on outer end
1015 
1016     if ((z_OutLo < zLo) || (z_OutLo > zHi))
1017       return false;
1018 
1019     const float cosh2Eta = 1.f + (pz * pz) / (ptIn * ptIn);
1020 
1021     const float drt_OutLo_InUp = (rt_OutLo - rt_InUp);
1022 
1023     const float r3_InUp = alpaka::math::sqrt(acc, z_InUp * z_InUp + rt_InUp * rt_InUp);
1024 
1025     float drt_InSeg = rt_InOut - rt_InLo;
1026 
1027     const float thetaMuls2 =
1028         (kMulsInGeV * kMulsInGeV) * (0.1f + 0.2f * (rt_OutLo - rt_InUp) / 50.f) * (r3_InUp / rt_InUp);
1029     const float muls2 = thetaMuls2 * 9.f / (ptCut * ptCut) * 16.f;
1030 
1031     float dzErr = (drt_OutLo_InUp * drt_OutLo_InUp) * (etaErr * etaErr) * cosh2Eta;
1032     dzErr += 0.03f * 0.03f;  // Approximately account for IT module size
1033     dzErr *= 9.f;            // 3 sigma
1034     dzErr += muls2 * (drt_OutLo_InUp * drt_OutLo_InUp) / 3.f * cosh2Eta;
1035     dzErr += zGeom * zGeom;
1036     dzErr = alpaka::math::sqrt(acc, dzErr);
1037 
1038     const float dzDrIn = pz / ptIn;
1039     const float zWindow = dzErr / drt_InSeg * drt_OutLo_InUp + zGeom;
1040     const float dzMean = dzDrIn * drt_OutLo_InUp *
1041                          (1.f + drt_OutLo_InUp * drt_OutLo_InUp * 4 * k2Rinv1GeVf * k2Rinv1GeVf / ptIn / ptIn /
1042                                     24.f);  // with curved path correction
1043     // Constructing upper and lower bound
1044     zLoPointed = z_InUp + dzMean - zWindow;
1045     zHiPointed = z_InUp + dzMean + zWindow;
1046 
1047     if ((z_OutLo < zLoPointed) || (z_OutLo > zHiPointed))
1048       return false;
1049 
1050     const float pvOffset = 0.1f / rt_OutLo;
1051     dPhiCut = alpha1GeV_OutLo + alpaka::math::sqrt(acc, muls2 + pvOffset * pvOffset);
1052 
1053     //no dphipos cut
1054     float midPointX = 0.5f * (x_InLo + x_OutLo);
1055     float midPointY = 0.5f * (y_InLo + y_OutLo);
1056 
1057     float diffX = x_OutLo - x_InLo;
1058     float diffY = y_OutLo - y_InLo;
1059 
1060     dPhi = cms::alpakatools::deltaPhi(acc, midPointX, midPointY, diffX, diffY);
1061 
1062     if (alpaka::math::abs(acc, dPhi) > dPhiCut)
1063       return false;
1064 
1065     //lots of array accesses below this...
1066 
1067     float alpha_InLo = __H2F(segments.dPhiChanges()[innerSegmentIndex]);
1068     float alpha_OutLo = __H2F(segments.dPhiChanges()[outerSegmentIndex]);
1069 
1070     bool isEC_lastLayer = modules.subdets()[outerOuterLowerModuleIndex] == Endcap and
1071                           modules.moduleType()[outerOuterLowerModuleIndex] == TwoS;
1072 
1073     float alpha_OutUp, alpha_OutUp_highEdge, alpha_OutUp_lowEdge;
1074     alpha_OutUp = cms::alpakatools::deltaPhi(acc, x_OutUp, y_OutUp, x_OutUp - x_OutLo, y_OutUp - y_OutLo);
1075 
1076     alpha_OutUp_highEdge = alpha_OutUp;
1077     alpha_OutUp_lowEdge = alpha_OutUp;
1078 
1079     float tl_axis_x = x_OutUp - x_InUp;
1080     float tl_axis_y = y_OutUp - y_InUp;
1081 
1082     float tl_axis_highEdge_x = tl_axis_x;
1083     float tl_axis_highEdge_y = tl_axis_y;
1084 
1085     float tl_axis_lowEdge_x = tl_axis_x;
1086     float tl_axis_lowEdge_y = tl_axis_y;
1087 
1088     betaIn = -cms::alpakatools::deltaPhi(acc, px, py, tl_axis_x, tl_axis_y);
1089     float betaInRHmin = betaIn;
1090     float betaInRHmax = betaIn;
1091 
1092     betaOut = -alpha_OutUp + cms::alpakatools::deltaPhi(acc, x_OutUp, y_OutUp, tl_axis_x, tl_axis_y);
1093 
1094     float betaOutRHmin = betaOut;
1095     float betaOutRHmax = betaOut;
1096 
1097     if (isEC_lastLayer) {
1098       alpha_OutUp_highEdge = cms::alpakatools::deltaPhi(acc,
1099                                                         mds.anchorHighEdgeX()[fourthMDIndex],
1100                                                         mds.anchorHighEdgeY()[fourthMDIndex],
1101                                                         mds.anchorHighEdgeX()[fourthMDIndex] - x_OutLo,
1102                                                         mds.anchorHighEdgeY()[fourthMDIndex] - y_OutLo);
1103       alpha_OutUp_lowEdge = cms::alpakatools::deltaPhi(acc,
1104                                                        mds.anchorLowEdgeX()[fourthMDIndex],
1105                                                        mds.anchorLowEdgeY()[fourthMDIndex],
1106                                                        mds.anchorLowEdgeX()[fourthMDIndex] - x_OutLo,
1107                                                        mds.anchorLowEdgeY()[fourthMDIndex] - y_OutLo);
1108 
1109       tl_axis_highEdge_x = mds.anchorHighEdgeX()[fourthMDIndex] - x_InUp;
1110       tl_axis_highEdge_y = mds.anchorHighEdgeY()[fourthMDIndex] - y_InUp;
1111       tl_axis_lowEdge_x = mds.anchorLowEdgeX()[fourthMDIndex] - x_InUp;
1112       tl_axis_lowEdge_y = mds.anchorLowEdgeY()[fourthMDIndex] - y_InUp;
1113 
1114       betaOutRHmin = -alpha_OutUp_highEdge + cms::alpakatools::deltaPhi(acc,
1115                                                                         mds.anchorHighEdgeX()[fourthMDIndex],
1116                                                                         mds.anchorHighEdgeY()[fourthMDIndex],
1117                                                                         tl_axis_highEdge_x,
1118                                                                         tl_axis_highEdge_y);
1119       betaOutRHmax = -alpha_OutUp_lowEdge + cms::alpakatools::deltaPhi(acc,
1120                                                                        mds.anchorLowEdgeX()[fourthMDIndex],
1121                                                                        mds.anchorLowEdgeY()[fourthMDIndex],
1122                                                                        tl_axis_lowEdge_x,
1123                                                                        tl_axis_lowEdge_y);
1124     }
1125 
1126     //beta computation
1127     float drt_tl_axis = alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y);
1128 
1129     //innerOuterAnchor - innerInnerAnchor
1130     const float rt_InSeg =
1131         alpaka::math::sqrt(acc, (x_InUp - x_InLo) * (x_InUp - x_InLo) + (y_InUp - y_InLo) * (y_InUp - y_InLo));
1132 
1133     //no betaIn cut for the pixels
1134     float betaAv = 0.5f * (betaIn + betaOut);
1135     pt_beta = ptIn;
1136 
1137     int lIn = 0;
1138     int lOut = isEC_lastLayer ? 11 : 5;
1139     float sdOut_dr =
1140         alpaka::math::sqrt(acc, (x_OutUp - x_OutLo) * (x_OutUp - x_OutLo) + (y_OutUp - y_OutLo) * (y_OutUp - y_OutLo));
1141     float sdOut_d = rt_OutUp - rt_OutLo;
1142 
1143     runDeltaBetaIterations(acc, betaIn, betaOut, betaAv, pt_beta, rt_InSeg, sdOut_dr, drt_tl_axis, lIn);
1144 
1145     const float betaInMMSF = (alpaka::math::abs(acc, betaInRHmin + betaInRHmax) > 0)
1146                                  ? (2.f * betaIn / alpaka::math::abs(acc, betaInRHmin + betaInRHmax))
1147                                  : 0.;  //mean value of min,max is the old betaIn
1148     const float betaOutMMSF = (alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax) > 0)
1149                                   ? (2.f * betaOut / alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax))
1150                                   : 0.;
1151     betaInRHmin *= betaInMMSF;
1152     betaInRHmax *= betaInMMSF;
1153     betaOutRHmin *= betaOutMMSF;
1154     betaOutRHmax *= betaOutMMSF;
1155 
1156     float min_ptBeta_ptBetaMax = alpaka::math::min(
1157         acc, alpaka::math::abs(acc, pt_beta), kPt_betaMax);  //need to confirm the range-out value of 7 GeV
1158     const float dBetaMuls2 = thetaMuls2 * 16.f / (min_ptBeta_ptBetaMax * min_ptBeta_ptBetaMax);
1159     const float alphaInAbsReg =
1160         alpaka::math::max(acc,
1161                           alpaka::math::abs(acc, alpha_InLo),
1162                           alpaka::math::asin(acc, alpaka::math::min(acc, rt_InUp * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1163     const float alphaOutAbsReg =
1164         alpaka::math::max(acc,
1165                           alpaka::math::abs(acc, alpha_OutLo),
1166                           alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1167     const float dBetaInLum = lIn < 11 ? 0.0f : alpaka::math::abs(acc, alphaInAbsReg * kDeltaZLum / z_InUp);
1168     const float dBetaOutLum = lOut < 11 ? 0.0f : alpaka::math::abs(acc, alphaOutAbsReg * kDeltaZLum / z_OutLo);
1169     const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum);
1170 
1171     const float sinDPhi = alpaka::math::sin(acc, dPhi);
1172     const float dBetaRIn2 = 0;  // TODO-RH
1173 
1174     float dBetaROut = 0;
1175     if (isEC_lastLayer) {
1176       dBetaROut = (alpaka::math::sqrt(acc,
1177                                       mds.anchorHighEdgeX()[fourthMDIndex] * mds.anchorHighEdgeX()[fourthMDIndex] +
1178                                           mds.anchorHighEdgeY()[fourthMDIndex] * mds.anchorHighEdgeY()[fourthMDIndex]) -
1179                    alpaka::math::sqrt(acc,
1180                                       mds.anchorLowEdgeX()[fourthMDIndex] * mds.anchorLowEdgeX()[fourthMDIndex] +
1181                                           mds.anchorLowEdgeY()[fourthMDIndex] * mds.anchorLowEdgeY()[fourthMDIndex])) *
1182                   sinDPhi / drt_tl_axis;
1183     }
1184 
1185     const float dBetaROut2 = dBetaROut * dBetaROut;
1186 
1187     //FIXME: need faster version
1188     betaOutCut = alpaka::math::asin(acc, alpaka::math::min(acc, drt_tl_axis * k2Rinv1GeVf / ptCut, kSinAlphaMax)) +
1189                  (0.02f / sdOut_d) + alpaka::math::sqrt(acc, dBetaLum2 + dBetaMuls2);
1190 
1191     //Cut #6: The real beta cut
1192     if (alpaka::math::abs(acc, betaOut) >= betaOutCut)
1193       return false;
1194     const float dBetaRes = 0.02f / alpaka::math::min(acc, sdOut_d, drt_InSeg);
1195     const float dBetaCut2 =
1196         (dBetaRes * dBetaRes * 2.0f + dBetaMuls2 + dBetaLum2 + dBetaRIn2 + dBetaROut2 +
1197          0.25f *
1198              (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)) *
1199              (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)));
1200     float dBeta = betaIn - betaOut;
1201     return dBeta * dBeta <= dBetaCut2;
1202   }
1203 
1204   template <typename TAcc>
1205   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runTripletDefaultAlgoPPEE(TAcc const& acc,
1206                                                                 ModulesConst modules,
1207                                                                 ObjectRangesConst ranges,
1208                                                                 MiniDoubletsConst mds,
1209                                                                 SegmentsConst segments,
1210                                                                 PixelSeedsConst pixelSeeds,
1211                                                                 uint16_t pixelModuleIndex,
1212                                                                 uint16_t outerInnerLowerModuleIndex,
1213                                                                 uint16_t outerOuterLowerModuleIndex,
1214                                                                 unsigned int innerSegmentIndex,
1215                                                                 unsigned int outerSegmentIndex,
1216                                                                 unsigned int firstMDIndex,
1217                                                                 unsigned int secondMDIndex,
1218                                                                 unsigned int thirdMDIndex,
1219                                                                 unsigned int fourthMDIndex,
1220                                                                 const float ptCut) {
1221     float dPhi, betaIn, betaOut, pt_beta, rtLo, rtHi, dPhiCut, betaOutCut;
1222 
1223     bool isPS_OutLo = (modules.moduleType()[outerInnerLowerModuleIndex] == PS);
1224 
1225     float z_InUp = mds.anchorZ()[secondMDIndex];
1226     float z_OutLo = mds.anchorZ()[thirdMDIndex];
1227 
1228     if (z_InUp * z_OutLo <= 0)
1229       return false;
1230 
1231     float rt_InLo = mds.anchorRt()[firstMDIndex];
1232     float rt_InUp = mds.anchorRt()[secondMDIndex];
1233     float rt_OutLo = mds.anchorRt()[thirdMDIndex];
1234     float rt_OutUp = mds.anchorRt()[fourthMDIndex];
1235 
1236     float x_InLo = mds.anchorX()[firstMDIndex];
1237     float x_InUp = mds.anchorX()[secondMDIndex];
1238     float x_OutLo = mds.anchorX()[thirdMDIndex];
1239     float x_OutUp = mds.anchorX()[fourthMDIndex];
1240 
1241     float y_InLo = mds.anchorY()[firstMDIndex];
1242     float y_InUp = mds.anchorY()[secondMDIndex];
1243     float y_OutLo = mds.anchorY()[thirdMDIndex];
1244     float y_OutUp = mds.anchorY()[fourthMDIndex];
1245 
1246     unsigned int pixelSegmentArrayIndex = innerSegmentIndex - ranges.segmentModuleIndices()[pixelModuleIndex];
1247 
1248     float ptIn = pixelSeeds.ptIn()[pixelSegmentArrayIndex];
1249     float ptSLo = ptIn;
1250     float px = pixelSeeds.px()[pixelSegmentArrayIndex];
1251     float py = pixelSeeds.py()[pixelSegmentArrayIndex];
1252     float pz = pixelSeeds.pz()[pixelSegmentArrayIndex];
1253     float ptErr = pixelSeeds.ptErr()[pixelSegmentArrayIndex];
1254     float etaErr = pixelSeeds.etaErr()[pixelSegmentArrayIndex];
1255 
1256     ptSLo = alpaka::math::max(acc, ptCut, ptSLo - 10.0f * alpaka::math::max(acc, ptErr, 0.005f * ptSLo));
1257     ptSLo = alpaka::math::min(acc, 10.0f, ptSLo);
1258 
1259     const float zpitch_InLo = 0.05f;
1260     float zpitch_OutLo = (isPS_OutLo ? kPixelPSZpitch : kStrip2SZpitch);
1261     float zGeom = zpitch_InLo + zpitch_OutLo;
1262 
1263     const float slope = alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / ptCut, kSinAlphaMax));
1264     const float dzDrtScale = alpaka::math::tan(acc, slope) / slope;  //FIXME: need approximate value
1265 
1266     const float dLum = alpaka::math::copysign(acc, kDeltaZLum, z_InUp);
1267     bool isOutSgInnerMDPS = modules.moduleType()[outerInnerLowerModuleIndex] == PS;
1268 
1269     const float rtGeom1 = isOutSgInnerMDPS
1270                               ? kPixelPSZpitch
1271                               : kStrip2SZpitch;  //FIXME: make this chosen by configuration for lay11,12 full PS
1272     const float zGeom1 = alpaka::math::copysign(acc, zGeom, z_InUp);  //used in B-E region
1273     rtLo = rt_InUp * (1.f + (z_OutLo - z_InUp - zGeom1) / (z_InUp + zGeom1 + dLum) / dzDrtScale) -
1274            rtGeom1;  //slope correction only on the lower end
1275 
1276     float zInForHi = z_InUp - zGeom1 - dLum;
1277     if (zInForHi * z_InUp < 0)
1278       zInForHi = alpaka::math::copysign(acc, 0.1f, z_InUp);
1279     rtHi = rt_InUp * (1.f + (z_OutLo - z_InUp + zGeom1) / zInForHi) + rtGeom1;
1280 
1281     // Cut #2: rt condition
1282     if ((rt_OutLo < rtLo) || (rt_OutLo > rtHi))
1283       return false;
1284 
1285     const float dzOutInAbs = alpaka::math::abs(acc, z_OutLo - z_InUp);
1286     const float cosh2Eta = 1.f + (pz * pz) / (ptIn * ptIn);
1287     const float multDzDr2 = (dzOutInAbs * dzOutInAbs) * cosh2Eta / ((cosh2Eta - 1.f) * (cosh2Eta - 1.f));
1288     const float r3_InUp = alpaka::math::sqrt(acc, z_InUp * z_InUp + rt_InUp * rt_InUp);
1289     const float thetaMuls2 =
1290         (kMulsInGeV * kMulsInGeV) * (0.1f + 0.2f * (rt_OutLo - rt_InUp) / 50.f) * (r3_InUp / rt_InUp);
1291     const float muls2 = thetaMuls2 * 9.f / (ptCut * ptCut) * 16.f;
1292 
1293     float drtErr = (etaErr * etaErr) * multDzDr2;
1294     drtErr += 0.03f * 0.03f;  // Approximately account for IT module size
1295     drtErr *= 9.f;            // 3 sigma
1296     drtErr += muls2 * multDzDr2 / 3.f * cosh2Eta;
1297     drtErr = alpaka::math::sqrt(acc, drtErr);
1298     const float drtDzIn = alpaka::math::abs(acc, ptIn / pz);
1299 
1300     const float drt_OutLo_InUp = (rt_OutLo - rt_InUp);  // drOutIn
1301 
1302     const float rtWindow = drtErr + rtGeom1;
1303     const float drtMean = drtDzIn * dzOutInAbs *
1304                           (1.f - drt_OutLo_InUp * drt_OutLo_InUp * 4 * k2Rinv1GeVf * k2Rinv1GeVf / ptIn / ptIn /
1305                                      24.f);  // with curved path correction
1306     const float rtLo_point = rt_InUp + drtMean - rtWindow;
1307     const float rtHi_point = rt_InUp + drtMean + rtWindow;
1308 
1309     // Cut #3: rt-z pointed
1310     if ((rt_OutLo < rtLo_point) || (rt_OutLo > rtHi_point))
1311       return false;
1312 
1313     const float alpha1GeV_OutLo =
1314         alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / ptCut, kSinAlphaMax));
1315     const float pvOffset = 0.1f / rt_OutLo;
1316     dPhiCut = alpha1GeV_OutLo + alpaka::math::sqrt(acc, muls2 + pvOffset * pvOffset);
1317 
1318     float midPointX = 0.5f * (x_InLo + x_OutLo);
1319     float midPointY = 0.5f * (y_InLo + y_OutLo);
1320 
1321     float diffX = x_OutLo - x_InLo;
1322     float diffY = y_OutLo - y_InLo;
1323 
1324     dPhi = cms::alpakatools::deltaPhi(acc, midPointX, midPointY, diffX, diffY);
1325 
1326     // Cut #5: deltaPhiChange
1327     if (alpaka::math::abs(acc, dPhi) > dPhiCut)
1328       return false;
1329 
1330     float alpha_InLo = __H2F(segments.dPhiChanges()[innerSegmentIndex]);
1331     float alpha_OutLo = __H2F(segments.dPhiChanges()[outerSegmentIndex]);
1332 
1333     bool isEC_lastLayer = modules.subdets()[outerOuterLowerModuleIndex] == Endcap and
1334                           modules.moduleType()[outerOuterLowerModuleIndex] == TwoS;
1335 
1336     float alpha_OutUp, alpha_OutUp_highEdge, alpha_OutUp_lowEdge;
1337 
1338     alpha_OutUp = cms::alpakatools::deltaPhi(acc, x_OutUp, y_OutUp, x_OutUp - x_OutLo, y_OutUp - y_OutLo);
1339     alpha_OutUp_highEdge = alpha_OutUp;
1340     alpha_OutUp_lowEdge = alpha_OutUp;
1341 
1342     float tl_axis_x = x_OutUp - x_InUp;
1343     float tl_axis_y = y_OutUp - y_InUp;
1344 
1345     float tl_axis_highEdge_x = tl_axis_x;
1346     float tl_axis_highEdge_y = tl_axis_y;
1347 
1348     float tl_axis_lowEdge_x = tl_axis_x;
1349     float tl_axis_lowEdge_y = tl_axis_y;
1350 
1351     betaIn = -cms::alpakatools::deltaPhi(acc, px, py, tl_axis_x, tl_axis_y);
1352     float betaInRHmin = betaIn;
1353     float betaInRHmax = betaIn;
1354 
1355     betaOut = -alpha_OutUp + cms::alpakatools::deltaPhi(acc, x_OutUp, y_OutUp, tl_axis_x, tl_axis_y);
1356     float betaOutRHmin = betaOut;
1357     float betaOutRHmax = betaOut;
1358 
1359     if (isEC_lastLayer) {
1360       alpha_OutUp_highEdge = cms::alpakatools::deltaPhi(acc,
1361                                                         mds.anchorHighEdgeX()[fourthMDIndex],
1362                                                         mds.anchorHighEdgeY()[fourthMDIndex],
1363                                                         mds.anchorHighEdgeX()[fourthMDIndex] - x_OutLo,
1364                                                         mds.anchorHighEdgeY()[fourthMDIndex] - y_OutLo);
1365       alpha_OutUp_lowEdge = cms::alpakatools::deltaPhi(acc,
1366                                                        mds.anchorLowEdgeX()[fourthMDIndex],
1367                                                        mds.anchorLowEdgeY()[fourthMDIndex],
1368                                                        mds.anchorLowEdgeX()[fourthMDIndex] - x_OutLo,
1369                                                        mds.anchorLowEdgeY()[fourthMDIndex] - y_OutLo);
1370 
1371       tl_axis_highEdge_x = mds.anchorHighEdgeX()[fourthMDIndex] - x_InUp;
1372       tl_axis_highEdge_y = mds.anchorHighEdgeY()[fourthMDIndex] - y_InUp;
1373       tl_axis_lowEdge_x = mds.anchorLowEdgeX()[fourthMDIndex] - x_InUp;
1374       tl_axis_lowEdge_y = mds.anchorLowEdgeY()[fourthMDIndex] - y_InUp;
1375 
1376       betaOutRHmin = -alpha_OutUp_highEdge + cms::alpakatools::deltaPhi(acc,
1377                                                                         mds.anchorHighEdgeX()[fourthMDIndex],
1378                                                                         mds.anchorHighEdgeY()[fourthMDIndex],
1379                                                                         tl_axis_highEdge_x,
1380                                                                         tl_axis_highEdge_y);
1381       betaOutRHmax = -alpha_OutUp_lowEdge + cms::alpakatools::deltaPhi(acc,
1382                                                                        mds.anchorLowEdgeX()[fourthMDIndex],
1383                                                                        mds.anchorLowEdgeY()[fourthMDIndex],
1384                                                                        tl_axis_lowEdge_x,
1385                                                                        tl_axis_lowEdge_y);
1386     }
1387 
1388     //beta computation
1389     float drt_tl_axis = alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y);
1390     //no betaIn cut for the pixels
1391     const float rt_InSeg =
1392         alpaka::math::sqrt(acc, (x_InUp - x_InLo) * (x_InUp - x_InLo) + (y_InUp - y_InLo) * (y_InUp - y_InLo));
1393 
1394     float betaAv = 0.5f * (betaIn + betaOut);
1395     pt_beta = ptIn;
1396 
1397     int lIn = 0;
1398     int lOut = isEC_lastLayer ? 11 : 5;
1399     float sdOut_dr =
1400         alpaka::math::sqrt(acc, (x_OutUp - x_OutLo) * (x_OutUp - x_OutLo) + (y_OutUp - y_OutLo) * (y_OutUp - y_OutLo));
1401     float sdOut_d = rt_OutUp - rt_OutLo;
1402 
1403     runDeltaBetaIterations(acc, betaIn, betaOut, betaAv, pt_beta, rt_InSeg, sdOut_dr, drt_tl_axis, lIn);
1404 
1405     const float betaInMMSF = (alpaka::math::abs(acc, betaInRHmin + betaInRHmax) > 0)
1406                                  ? (2.f * betaIn / alpaka::math::abs(acc, betaInRHmin + betaInRHmax))
1407                                  : 0.;  //mean value of min,max is the old betaIn
1408     const float betaOutMMSF = (alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax) > 0)
1409                                   ? (2.f * betaOut / alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax))
1410                                   : 0.;
1411     betaInRHmin *= betaInMMSF;
1412     betaInRHmax *= betaInMMSF;
1413     betaOutRHmin *= betaOutMMSF;
1414     betaOutRHmax *= betaOutMMSF;
1415 
1416     float min_ptBeta_ptBetaMax = alpaka::math::min(
1417         acc, alpaka::math::abs(acc, pt_beta), kPt_betaMax);  //need to confirm the range-out value of 7 GeV
1418     const float dBetaMuls2 = thetaMuls2 * 16.f / (min_ptBeta_ptBetaMax * min_ptBeta_ptBetaMax);
1419 
1420     const float alphaInAbsReg =
1421         alpaka::math::max(acc,
1422                           alpaka::math::abs(acc, alpha_InLo),
1423                           alpaka::math::asin(acc, alpaka::math::min(acc, rt_InUp * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1424     const float alphaOutAbsReg =
1425         alpaka::math::max(acc,
1426                           alpaka::math::abs(acc, alpha_OutLo),
1427                           alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1428     const float dBetaInLum = lIn < 11 ? 0.0f : alpaka::math::abs(acc, alphaInAbsReg * kDeltaZLum / z_InUp);
1429     const float dBetaOutLum = lOut < 11 ? 0.0f : alpaka::math::abs(acc, alphaOutAbsReg * kDeltaZLum / z_OutLo);
1430     const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum);
1431 
1432     const float sinDPhi = alpaka::math::sin(acc, dPhi);
1433     const float dBetaRIn2 = 0;  // TODO-RH
1434 
1435     float dBetaROut = 0;
1436     if (isEC_lastLayer) {
1437       dBetaROut = (alpaka::math::sqrt(acc,
1438                                       mds.anchorHighEdgeX()[fourthMDIndex] * mds.anchorHighEdgeX()[fourthMDIndex] +
1439                                           mds.anchorHighEdgeY()[fourthMDIndex] * mds.anchorHighEdgeY()[fourthMDIndex]) -
1440                    alpaka::math::sqrt(acc,
1441                                       mds.anchorLowEdgeX()[fourthMDIndex] * mds.anchorLowEdgeX()[fourthMDIndex] +
1442                                           mds.anchorLowEdgeY()[fourthMDIndex] * mds.anchorLowEdgeY()[fourthMDIndex])) *
1443                   sinDPhi / drt_tl_axis;
1444     }
1445 
1446     const float dBetaROut2 = dBetaROut * dBetaROut;
1447 
1448     betaOutCut =
1449         alpaka::math::asin(
1450             acc, alpaka::math::min(acc, drt_tl_axis * k2Rinv1GeVf / ptCut, kSinAlphaMax))  //FIXME: need faster version
1451         + (0.02f / sdOut_d) + alpaka::math::sqrt(acc, dBetaLum2 + dBetaMuls2);
1452 
1453     //Cut #6: The real beta cut
1454     if (alpaka::math::abs(acc, betaOut) >= betaOutCut)
1455       return false;
1456 
1457     float drt_InSeg = rt_InUp - rt_InLo;
1458 
1459     const float dBetaRes = 0.02f / alpaka::math::min(acc, sdOut_d, drt_InSeg);
1460     const float dBetaCut2 =
1461         (dBetaRes * dBetaRes * 2.0f + dBetaMuls2 + dBetaLum2 + dBetaRIn2 + dBetaROut2 +
1462          0.25f *
1463              (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)) *
1464              (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)));
1465     float dBeta = betaIn - betaOut;
1466     return dBeta * dBeta <= dBetaCut2;
1467   }
1468 
1469 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE::lst
1470 #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!