Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-06-06 01:33:38 

0001 #ifndef RecoTracker_LSTCore_src_alpaka_Quintuplet_h
0002 #define RecoTracker_LSTCore_src_alpaka_Quintuplet_h
0003 
0004 #include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
0005 #include "FWCore/Utilities/interface/isFinite.h"
0006 #include "FWCore/Utilities/interface/CMSUnrollLoop.h"
0007 
0008 #include "RecoTracker/LSTCore/interface/ObjectRangesSoA.h"
0009 #include "RecoTracker/LSTCore/interface/MiniDoubletsSoA.h"
0010 #include "RecoTracker/LSTCore/interface/SegmentsSoA.h"
0011 #include "RecoTracker/LSTCore/interface/TripletsSoA.h"
0012 #include "RecoTracker/LSTCore/interface/QuintupletsSoA.h"
0013 #include "RecoTracker/LSTCore/interface/alpaka/Common.h"
0014 #include "RecoTracker/LSTCore/interface/ModulesSoA.h"
0015 #include "RecoTracker/LSTCore/interface/EndcapGeometry.h"
0016 #include "RecoTracker/LSTCore/interface/ObjectRangesSoA.h"
0017 #include "RecoTracker/LSTCore/interface/Circle.h"
0018 
0019 #include "NeuralNetwork.h"
0020 
0021 namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
0022   ALPAKA_FN_ACC ALPAKA_FN_INLINE void addQuintupletToMemory(TripletsConst triplets,
0023                                                             Quintuplets quintuplets,
0024                                                             unsigned int innerTripletIndex,
0025                                                             unsigned int outerTripletIndex,
0026                                                             uint16_t lowerModule1,
0027                                                             uint16_t lowerModule2,
0028                                                             uint16_t lowerModule3,
0029                                                             uint16_t lowerModule4,
0030                                                             uint16_t lowerModule5,
0031                                                             float innerRadius,
0032                                                             float bridgeRadius,
0033                                                             float outerRadius,
0034                                                             float regressionCenterX,
0035                                                             float regressionCenterY,
0036                                                             float regressionRadius,
0037                                                             float rzChiSquared,
0038                                                             float rPhiChiSquared,
0039                                                             float nonAnchorChiSquared,
0040                                                             float dBeta1,
0041                                                             float dBeta2,
0042                                                             float pt,
0043                                                             float eta,
0044                                                             float phi,
0045                                                             float scores,
0046                                                             uint8_t layer,
0047                                                             unsigned int quintupletIndex,
0048                                                             const float (&t5Embed)[Params_T5::kEmbed],
0049                                                             bool tightCutFlag) {
0050     quintuplets.tripletIndices()[quintupletIndex][0] = innerTripletIndex;
0051     quintuplets.tripletIndices()[quintupletIndex][1] = outerTripletIndex;
0052 
0053     quintuplets.lowerModuleIndices()[quintupletIndex][0] = lowerModule1;
0054     quintuplets.lowerModuleIndices()[quintupletIndex][1] = lowerModule2;
0055     quintuplets.lowerModuleIndices()[quintupletIndex][2] = lowerModule3;
0056     quintuplets.lowerModuleIndices()[quintupletIndex][3] = lowerModule4;
0057     quintuplets.lowerModuleIndices()[quintupletIndex][4] = lowerModule5;
0058     quintuplets.innerRadius()[quintupletIndex] = __F2H(innerRadius);
0059     quintuplets.outerRadius()[quintupletIndex] = __F2H(outerRadius);
0060     quintuplets.pt()[quintupletIndex] = __F2H(pt);
0061     quintuplets.eta()[quintupletIndex] = __F2H(eta);
0062     quintuplets.phi()[quintupletIndex] = __F2H(phi);
0063     quintuplets.score_rphisum()[quintupletIndex] = __F2H(scores);
0064     quintuplets.isDup()[quintupletIndex] = 0;
0065     quintuplets.tightCutFlag()[quintupletIndex] = tightCutFlag;
0066     quintuplets.regressionRadius()[quintupletIndex] = regressionRadius;
0067     quintuplets.regressionCenterX()[quintupletIndex] = regressionCenterX;
0068     quintuplets.regressionCenterY()[quintupletIndex] = regressionCenterY;
0069     quintuplets.logicalLayers()[quintupletIndex][0] = triplets.logicalLayers()[innerTripletIndex][0];
0070     quintuplets.logicalLayers()[quintupletIndex][1] = triplets.logicalLayers()[innerTripletIndex][1];
0071     quintuplets.logicalLayers()[quintupletIndex][2] = triplets.logicalLayers()[innerTripletIndex][2];
0072     quintuplets.logicalLayers()[quintupletIndex][3] = triplets.logicalLayers()[outerTripletIndex][1];
0073     quintuplets.logicalLayers()[quintupletIndex][4] = triplets.logicalLayers()[outerTripletIndex][2];
0074 
0075     quintuplets.hitIndices()[quintupletIndex][0] = triplets.hitIndices()[innerTripletIndex][0];
0076     quintuplets.hitIndices()[quintupletIndex][1] = triplets.hitIndices()[innerTripletIndex][1];
0077     quintuplets.hitIndices()[quintupletIndex][2] = triplets.hitIndices()[innerTripletIndex][2];
0078     quintuplets.hitIndices()[quintupletIndex][3] = triplets.hitIndices()[innerTripletIndex][3];
0079     quintuplets.hitIndices()[quintupletIndex][4] = triplets.hitIndices()[innerTripletIndex][4];
0080     quintuplets.hitIndices()[quintupletIndex][5] = triplets.hitIndices()[innerTripletIndex][5];
0081     quintuplets.hitIndices()[quintupletIndex][6] = triplets.hitIndices()[outerTripletIndex][2];
0082     quintuplets.hitIndices()[quintupletIndex][7] = triplets.hitIndices()[outerTripletIndex][3];
0083     quintuplets.hitIndices()[quintupletIndex][8] = triplets.hitIndices()[outerTripletIndex][4];
0084     quintuplets.hitIndices()[quintupletIndex][9] = triplets.hitIndices()[outerTripletIndex][5];
0085     quintuplets.bridgeRadius()[quintupletIndex] = bridgeRadius;
0086     quintuplets.rzChiSquared()[quintupletIndex] = rzChiSquared;
0087     quintuplets.chiSquared()[quintupletIndex] = rPhiChiSquared;
0088     quintuplets.nonAnchorChiSquared()[quintupletIndex] = nonAnchorChiSquared;
0089     quintuplets.dBeta1()[quintupletIndex] = dBeta1;
0090     quintuplets.dBeta2()[quintupletIndex] = dBeta2;
0091 
0092     CMS_UNROLL_LOOP
0093     for (unsigned int i = 0; i < Params_T5::kEmbed; ++i) {
0094       quintuplets.t5Embed()[quintupletIndex][i] = t5Embed[i];
0095     }
0096   }
0097 
0098   //bounds can be found at http://uaf-10.t2.ucsd.edu/~bsathian/SDL/T5_RZFix/t5_rz_thresholds.txt
0099   template <typename TAcc>
0100   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passT5RZConstraint(TAcc const& acc,
0101                                                          ModulesConst modules,
0102                                                          MiniDoubletsConst mds,
0103                                                          unsigned int firstMDIndex,
0104                                                          unsigned int secondMDIndex,
0105                                                          unsigned int thirdMDIndex,
0106                                                          unsigned int fourthMDIndex,
0107                                                          unsigned int fifthMDIndex,
0108                                                          uint16_t lowerModuleIndex1,
0109                                                          uint16_t lowerModuleIndex2,
0110                                                          uint16_t lowerModuleIndex3,
0111                                                          uint16_t lowerModuleIndex4,
0112                                                          uint16_t lowerModuleIndex5,
0113                                                          float& rzChiSquared,
0114                                                          float inner_pt,
0115                                                          float innerRadius,
0116                                                          float g,
0117                                                          float f,
0118                                                          bool& tightCutFlag) {
0119     //(g,f) is the center of the circle fitted by the innermost 3 points on x,y coordinates
0120     const float rt1 = mds.anchorRt()[firstMDIndex] / 100;  //in the unit of m instead of cm
0121     const float rt2 = mds.anchorRt()[secondMDIndex] / 100;
0122     const float rt3 = mds.anchorRt()[thirdMDIndex] / 100;
0123     const float rt4 = mds.anchorRt()[fourthMDIndex] / 100;
0124     const float rt5 = mds.anchorRt()[fifthMDIndex] / 100;
0125 
0126     const float z1 = mds.anchorZ()[firstMDIndex] / 100;
0127     const float z2 = mds.anchorZ()[secondMDIndex] / 100;
0128     const float z3 = mds.anchorZ()[thirdMDIndex] / 100;
0129     const float z4 = mds.anchorZ()[fourthMDIndex] / 100;
0130     const float z5 = mds.anchorZ()[fifthMDIndex] / 100;
0131 
0132     // Using lst_layer numbering convention defined in ModuleMethods.h
0133     const int layer1 = modules.lstLayers()[lowerModuleIndex1];
0134     const int layer2 = modules.lstLayers()[lowerModuleIndex2];
0135     const int layer3 = modules.lstLayers()[lowerModuleIndex3];
0136     const int layer4 = modules.lstLayers()[lowerModuleIndex4];
0137     const int layer5 = modules.lstLayers()[lowerModuleIndex5];
0138 
0139     //slope computed using the internal T3s
0140     const int moduleType1 = modules.moduleType()[lowerModuleIndex1];  //0 is ps, 1 is 2s
0141     const int moduleType2 = modules.moduleType()[lowerModuleIndex2];
0142     const int moduleType3 = modules.moduleType()[lowerModuleIndex3];
0143     const int moduleType4 = modules.moduleType()[lowerModuleIndex4];
0144     const int moduleType5 = modules.moduleType()[lowerModuleIndex5];
0145 
0146     const float x1 = mds.anchorX()[firstMDIndex] / 100;
0147     const float x2 = mds.anchorX()[secondMDIndex] / 100;
0148     const float x3 = mds.anchorX()[thirdMDIndex] / 100;
0149     const float x4 = mds.anchorX()[fourthMDIndex] / 100;
0150     const float y1 = mds.anchorY()[firstMDIndex] / 100;
0151     const float y2 = mds.anchorY()[secondMDIndex] / 100;
0152     const float y3 = mds.anchorY()[thirdMDIndex] / 100;
0153     const float y4 = mds.anchorY()[fourthMDIndex] / 100;
0154 
0155     float residual = 0;
0156     float error2 = 0;
0157     float x_center = g / 100, y_center = f / 100;
0158     float x_init = mds.anchorX()[thirdMDIndex] / 100;
0159     float y_init = mds.anchorY()[thirdMDIndex] / 100;
0160     float z_init = mds.anchorZ()[thirdMDIndex] / 100;
0161     float rt_init = mds.anchorRt()[thirdMDIndex] / 100;  //use the second MD as initial point
0162 
0163     if (moduleType3 == 1)  // 1: if MD3 is in 2s layer
0164     {
0165       x_init = mds.anchorX()[secondMDIndex] / 100;
0166       y_init = mds.anchorY()[secondMDIndex] / 100;
0167       z_init = mds.anchorZ()[secondMDIndex] / 100;
0168       rt_init = mds.anchorRt()[secondMDIndex] / 100;
0169     }
0170 
0171     // start from a circle of inner T3.
0172     // to determine the charge
0173     int charge = 0;
0174     float slope3c = (y3 - y_center) / (x3 - x_center);
0175     float slope1c = (y1 - y_center) / (x1 - x_center);
0176     // these 4 "if"s basically separate the x-y plane into 4 quarters. It determines geometrically how a circle and line slope goes and their positions, and we can get the charges correspondingly.
0177     if ((y3 - y_center) > 0 && (y1 - y_center) > 0) {
0178       if (slope1c > 0 && slope3c < 0)
0179         charge = -1;  // on x axis of a quarter, 3 hits go anti-clockwise
0180       else if (slope1c < 0 && slope3c > 0)
0181         charge = 1;  // on x axis of a quarter, 3 hits go clockwise
0182       else if (slope3c > slope1c)
0183         charge = -1;
0184       else if (slope3c < slope1c)
0185         charge = 1;
0186     } else if ((y3 - y_center) < 0 && (y1 - y_center) < 0) {
0187       if (slope1c < 0 && slope3c > 0)
0188         charge = 1;
0189       else if (slope1c > 0 && slope3c < 0)
0190         charge = -1;
0191       else if (slope3c > slope1c)
0192         charge = -1;
0193       else if (slope3c < slope1c)
0194         charge = 1;
0195     } else if ((y3 - y_center) < 0 && (y1 - y_center) > 0) {
0196       if ((x3 - x_center) > 0 && (x1 - x_center) > 0)
0197         charge = 1;
0198       else if ((x3 - x_center) < 0 && (x1 - x_center) < 0)
0199         charge = -1;
0200     } else if ((y3 - y_center) > 0 && (y1 - y_center) < 0) {
0201       if ((x3 - x_center) > 0 && (x1 - x_center) > 0)
0202         charge = -1;
0203       else if ((x3 - x_center) < 0 && (x1 - x_center) < 0)
0204         charge = 1;
0205     }
0206 
0207     float pseudo_phi = alpaka::math::atan(
0208         acc, (y_init - y_center) / (x_init - x_center));  //actually represent pi/2-phi, wrt helix axis z
0209     float Pt = inner_pt, Px = Pt * alpaka::math::abs(acc, alpaka::math::sin(acc, pseudo_phi)),
0210           Py = Pt * alpaka::math::abs(acc, cos(pseudo_phi));
0211 
0212     // Above line only gives you the correct value of Px and Py, but signs of Px and Py calculated below.
0213     // We look at if the circle is clockwise or anti-clock wise, to make it simpler, we separate the x-y plane into 4 quarters.
0214     if (x_init > x_center && y_init > y_center)  //1st quad
0215     {
0216       if (charge == 1)
0217         Py = -Py;
0218       if (charge == -1)
0219         Px = -Px;
0220     }
0221     if (x_init < x_center && y_init > y_center)  //2nd quad
0222     {
0223       if (charge == -1) {
0224         Px = -Px;
0225         Py = -Py;
0226       }
0227     }
0228     if (x_init < x_center && y_init < y_center)  //3rd quad
0229     {
0230       if (charge == 1)
0231         Px = -Px;
0232       if (charge == -1)
0233         Py = -Py;
0234     }
0235     if (x_init > x_center && y_init < y_center)  //4th quad
0236     {
0237       if (charge == 1) {
0238         Px = -Px;
0239         Py = -Py;
0240       }
0241     }
0242 
0243     // But if the initial T5 curve goes across quarters(i.e. cross axis to separate the quarters), need special redeclaration of Px,Py signs on these to avoid errors
0244     if (moduleType3 == 0) {  // 0 is ps
0245       if (x4 < x3 && x3 < x2)
0246         Px = -alpaka::math::abs(acc, Px);
0247       else if (x4 > x3 && x3 > x2)
0248         Px = alpaka::math::abs(acc, Px);
0249       if (y4 < y3 && y3 < y2)
0250         Py = -alpaka::math::abs(acc, Py);
0251       else if (y4 > y3 && y3 > y2)
0252         Py = alpaka::math::abs(acc, Py);
0253     } else if (moduleType3 == 1)  // 1 is 2s
0254     {
0255       if (x3 < x2 && x2 < x1)
0256         Px = -alpaka::math::abs(acc, Px);
0257       else if (x3 > x2 && x2 > x1)
0258         Px = alpaka::math::abs(acc, Px);
0259       if (y3 < y2 && y2 < y1)
0260         Py = -alpaka::math::abs(acc, Py);
0261       else if (y3 > y2 && y2 > y1)
0262         Py = alpaka::math::abs(acc, Py);
0263     }
0264 
0265     //to get Pz, we use pt/pz=ds/dz, ds is the arclength between MD1 and MD3.
0266     float AO = alpaka::math::sqrt(acc, (x1 - x_center) * (x1 - x_center) + (y1 - y_center) * (y1 - y_center));
0267     float BO =
0268         alpaka::math::sqrt(acc, (x_init - x_center) * (x_init - x_center) + (y_init - y_center) * (y_init - y_center));
0269     float AB2 = (x1 - x_init) * (x1 - x_init) + (y1 - y_init) * (y1 - y_init);
0270     float dPhi = alpaka::math::acos(acc, (AO * AO + BO * BO - AB2) / (2 * AO * BO));
0271     float ds = innerRadius / 100 * dPhi;
0272 
0273     float Pz = (z_init - z1) / ds * Pt;
0274     float p = alpaka::math::sqrt(acc, Px * Px + Py * Py + Pz * Pz);
0275 
0276     float a = -2.f * k2Rinv1GeVf * 100 * charge;  // multiply by 100 to make the correct length units
0277 
0278     float zsi, rtsi;
0279     int layeri, moduleTypei;
0280     rzChiSquared = 0;
0281     for (size_t i = 2; i < 6; i++) {
0282       if (i == 2) {
0283         zsi = z2;
0284         rtsi = rt2;
0285         layeri = layer2;
0286         moduleTypei = moduleType2;
0287       } else if (i == 3) {
0288         zsi = z3;
0289         rtsi = rt3;
0290         layeri = layer3;
0291         moduleTypei = moduleType3;
0292       } else if (i == 4) {
0293         zsi = z4;
0294         rtsi = rt4;
0295         layeri = layer4;
0296         moduleTypei = moduleType4;
0297       } else if (i == 5) {
0298         zsi = z5;
0299         rtsi = rt5;
0300         layeri = layer5;
0301         moduleTypei = moduleType5;
0302       }
0303 
0304       if (moduleType3 == 0) {  //0: ps
0305         if (i == 3)
0306           continue;
0307       } else {
0308         if (i == 2)
0309           continue;
0310       }
0311 
0312       // calculation is copied from PixelTriplet.h computePT3RZChiSquared
0313       float diffr = 0, diffz = 0;
0314 
0315       float rou = a / p;
0316       // for endcap
0317       float s = (zsi - z_init) * p / Pz;
0318       float x = x_init + Px / a * alpaka::math::sin(acc, rou * s) - Py / a * (1 - alpaka::math::cos(acc, rou * s));
0319       float y = y_init + Py / a * alpaka::math::sin(acc, rou * s) + Px / a * (1 - alpaka::math::cos(acc, rou * s));
0320       diffr = (rtsi - alpaka::math::sqrt(acc, x * x + y * y)) * 100;
0321 
0322       // for barrel
0323       if (layeri <= 6) {
0324         float paraA =
0325             rt_init * rt_init + 2 * (Px * Px + Py * Py) / (a * a) + 2 * (y_init * Px - x_init * Py) / a - rtsi * rtsi;
0326         float paraB = 2 * (x_init * Px + y_init * Py) / a;
0327         float paraC = 2 * (y_init * Px - x_init * Py) / a + 2 * (Px * Px + Py * Py) / (a * a);
0328         float A = paraB * paraB + paraC * paraC;
0329         float B = 2 * paraA * paraB;
0330         float C = paraA * paraA - paraC * paraC;
0331         float sol1 = (-B + alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
0332         float sol2 = (-B - alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
0333         float solz1 = alpaka::math::asin(acc, sol1) / rou * Pz / p + z_init;
0334         float solz2 = alpaka::math::asin(acc, sol2) / rou * Pz / p + z_init;
0335         float diffz1 = (solz1 - zsi) * 100;
0336         float diffz2 = (solz2 - zsi) * 100;
0337         if (edm::isNotFinite(diffz1))
0338           diffz = diffz2;
0339         else if (edm::isNotFinite(diffz2))
0340           diffz = diffz1;
0341         else {
0342           diffz = (alpaka::math::abs(acc, diffz1) < alpaka::math::abs(acc, diffz2)) ? diffz1 : diffz2;
0343         }
0344       }
0345       residual = (layeri > 6) ? diffr : diffz;
0346 
0347       //PS Modules
0348       if (moduleTypei == 0) {
0349         error2 = kPixelPSZpitch * kPixelPSZpitch;
0350       } else  //2S modules
0351       {
0352         error2 = kStrip2SZpitch * kStrip2SZpitch;
0353       }
0354 
0355       //check the tilted module, side: PosZ, NegZ, Center(for not tilted)
0356       float drdz;
0357       short side, subdets;
0358       if (i == 2) {
0359         drdz = alpaka::math::abs(acc, modules.drdzs()[lowerModuleIndex2]);
0360         side = modules.sides()[lowerModuleIndex2];
0361         subdets = modules.subdets()[lowerModuleIndex2];
0362       }
0363       if (i == 3) {
0364         drdz = alpaka::math::abs(acc, modules.drdzs()[lowerModuleIndex3]);
0365         side = modules.sides()[lowerModuleIndex3];
0366         subdets = modules.subdets()[lowerModuleIndex3];
0367       }
0368       if (i == 2 || i == 3) {
0369         residual = (layeri <= 6 && ((side == Center) or (drdz < 1))) ? diffz : diffr;
0370         float projection_missing2 = 1.f;
0371         if (drdz < 1)
0372           projection_missing2 =
0373               ((subdets == Endcap) or (side == Center)) ? 1.f : 1.f / (1 + drdz * drdz);  // cos(atan(drdz)), if dr/dz<1
0374         if (drdz > 1)
0375           projection_missing2 = ((subdets == Endcap) or (side == Center))
0376                                     ? 1.f
0377                                     : (drdz * drdz) / (1 + drdz * drdz);  //sin(atan(drdz)), if dr/dz>1
0378         error2 = error2 * projection_missing2;
0379       }
0380       rzChiSquared += 12 * (residual * residual) / error2;
0381     }
0382     // for set rzchi2 cut
0383     // if the 5 points are linear, helix calculation gives nan
0384     if (inner_pt > 100 || edm::isNotFinite(rzChiSquared)) {
0385       float slope;
0386       if (moduleType1 == 0 and moduleType2 == 0 and moduleType3 == 1)  //PSPS2S
0387       {
0388         slope = (z2 - z1) / (rt2 - rt1);
0389       } else {
0390         slope = (z3 - z1) / (rt3 - rt1);
0391       }
0392       float residual4_linear = (layer4 <= 6) ? ((z4 - z1) - slope * (rt4 - rt1)) : ((rt4 - rt1) - (z4 - z1) / slope);
0393       float residual5_linear = (layer4 <= 6) ? ((z5 - z1) - slope * (rt5 - rt1)) : ((rt5 - rt1) - (z5 - z1) / slope);
0394 
0395       // creating a chi squared type quantity
0396       // 0-> PS, 1->2S
0397       residual4_linear = (moduleType4 == 0) ? residual4_linear / kPixelPSZpitch : residual4_linear / kStrip2SZpitch;
0398       residual5_linear = (moduleType5 == 0) ? residual5_linear / kPixelPSZpitch : residual5_linear / kStrip2SZpitch;
0399       residual4_linear = residual4_linear * 100;
0400       residual5_linear = residual5_linear * 100;
0401 
0402       rzChiSquared = 12 * (residual4_linear * residual4_linear + residual5_linear * residual5_linear);
0403       return rzChiSquared < 4.677f;
0404     }
0405 
0406     // when building T5, apply 99% chi2 cuts as default, and add to pT5 collection. But when adding T5 to TC collections, apply 95% cut to reduce the fake rate
0407     tightCutFlag = false;
0408     // The category numbers are related to module regions and layers, decoding of the region numbers can be found here in slide 2 table. https://github.com/SegmentLinking/TrackLooper/files/11420927/part.2.pdf
0409     // The commented numbers after each case is the region code, and can look it up from the table to see which category it belongs to. For example, //0 means T5 built with Endcap 1,2,3,4,5 ps modules
0410     if (layer1 == 7 and layer2 == 8 and layer3 == 9 and layer4 == 10 and layer5 == 11)  //0
0411     {
0412       if (rzChiSquared < 94.470f)
0413         tightCutFlag = true;
0414       return true;
0415     } else if (layer1 == 7 and layer2 == 8 and layer3 == 9 and layer4 == 10 and layer5 == 16)  //1
0416     {
0417       if (rzChiSquared < 22.099f)
0418         tightCutFlag = true;
0419       return rzChiSquared < 37.956f;
0420     } else if (layer1 == 7 and layer2 == 8 and layer3 == 9 and layer4 == 15 and layer5 == 16)  //2
0421     {
0422       if (rzChiSquared < 7.992f)
0423         tightCutFlag = true;
0424       return rzChiSquared < 11.622f;
0425     } else if (layer1 == 1 and layer2 == 7 and layer3 == 8 and layer4 == 9) {
0426       if (layer5 == 10)  //3
0427       {
0428         if (rzChiSquared < 111.390f)
0429           tightCutFlag = true;
0430         return true;
0431       }
0432       if (layer5 == 15)  //4
0433       {
0434         if (rzChiSquared < 18.351f)
0435           tightCutFlag = true;
0436         return rzChiSquared < 37.941f;
0437       }
0438     } else if (layer1 == 1 and layer2 == 2 and layer3 == 7) {
0439       if (layer4 == 8 and layer5 == 9)  //5
0440       {
0441         if (rzChiSquared < 116.148f)
0442           tightCutFlag = true;
0443         return true;
0444       }
0445       if (layer4 == 8 and layer5 == 14)  //6
0446       {
0447         if (rzChiSquared < 19.352f)
0448           tightCutFlag = true;
0449         return rzChiSquared < 52.561f;
0450       } else if (layer4 == 13 and layer5 == 14)  //7
0451       {
0452         if (rzChiSquared < 10.392f)
0453           tightCutFlag = true;
0454         return rzChiSquared < 13.76f;
0455       }
0456     } else if (layer1 == 1 and layer2 == 2 and layer3 == 3) {
0457       if (layer4 == 7 and layer5 == 8)  //8
0458       {
0459         if (rzChiSquared < 27.824f)
0460           tightCutFlag = true;
0461         return rzChiSquared < 44.247f;
0462       } else if (layer4 == 7 and layer5 == 13)  //9
0463       {
0464         if (rzChiSquared < 18.145f)
0465           tightCutFlag = true;
0466         return rzChiSquared < 33.752f;
0467       } else if (layer4 == 12 and layer5 == 13)  //10
0468       {
0469         if (rzChiSquared < 13.308f)
0470           tightCutFlag = true;
0471         return rzChiSquared < 21.213f;
0472       } else if (layer4 == 4 and layer5 == 5)  //11
0473       {
0474         if (rzChiSquared < 15.627f)
0475           tightCutFlag = true;
0476         return rzChiSquared < 29.035f;
0477       } else if (layer4 == 4 and layer5 == 12)  //12
0478       {
0479         if (rzChiSquared < 14.64f)
0480           tightCutFlag = true;
0481         return rzChiSquared < 23.037f;
0482       }
0483     } else if (layer1 == 2 and layer2 == 7 and layer3 == 8) {
0484       if (layer4 == 9 and layer5 == 15)  //14
0485       {
0486         if (rzChiSquared < 24.662f)
0487           tightCutFlag = true;
0488         return rzChiSquared < 41.036f;
0489       } else if (layer4 == 14 and layer5 == 15)  //15
0490       {
0491         if (rzChiSquared < 8.866f)
0492           tightCutFlag = true;
0493         return rzChiSquared < 14.092f;
0494       }
0495     } else if (layer1 == 2 and layer2 == 3 and layer3 == 7) {
0496       if (layer4 == 8 and layer5 == 14)  //16
0497       {
0498         if (rzChiSquared < 23.730f)
0499           tightCutFlag = true;
0500         return rzChiSquared < 23.748f;
0501       }
0502       if (layer4 == 13 and layer5 == 14)  //17
0503       {
0504         if (rzChiSquared < 10.772f)
0505           tightCutFlag = true;
0506         return rzChiSquared < 17.945f;
0507       }
0508     } else if (layer1 == 2 and layer2 == 3 and layer3 == 4) {
0509       if (layer4 == 5 and layer5 == 6)  //18
0510       {
0511         if (rzChiSquared < 6.065f)
0512           tightCutFlag = true;
0513         return rzChiSquared < 8.803f;
0514       } else if (layer4 == 5 and layer5 == 12)  //19
0515       {
0516         if (rzChiSquared < 5.693f)
0517           tightCutFlag = true;
0518         return rzChiSquared < 7.930f;
0519       }
0520 
0521       else if (layer4 == 12 and layer5 == 13)  //20
0522       {
0523         if (rzChiSquared < 5.473f)
0524           tightCutFlag = true;
0525         return rzChiSquared < 7.626f;
0526       }
0527     }
0528     return true;
0529   }
0530 
0531   template <typename TAcc>
0532   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool T5HasCommonMiniDoublet(TripletsConst triplets,
0533                                                              SegmentsConst segments,
0534                                                              unsigned int innerTripletIndex,
0535                                                              unsigned int outerTripletIndex) {
0536     unsigned int innerOuterSegmentIndex = triplets.segmentIndices()[innerTripletIndex][1];
0537     unsigned int outerInnerSegmentIndex = triplets.segmentIndices()[outerTripletIndex][0];
0538     unsigned int innerOuterOuterMiniDoubletIndex =
0539         segments.mdIndices()[innerOuterSegmentIndex][1];  //inner triplet outer segment outer MD index
0540     unsigned int outerInnerInnerMiniDoubletIndex =
0541         segments.mdIndices()[outerInnerSegmentIndex][0];  //outer triplet inner segment inner MD index
0542 
0543     return (innerOuterOuterMiniDoubletIndex == outerInnerInnerMiniDoubletIndex);
0544   }
0545 
0546   template <typename TAcc>
0547   ALPAKA_FN_ACC ALPAKA_FN_INLINE void computeSigmasForRegression(TAcc const& acc,
0548                                                                  ModulesConst modules,
0549                                                                  const uint16_t* lowerModuleIndices,
0550                                                                  float* delta1,
0551                                                                  float* delta2,
0552                                                                  float* slopes,
0553                                                                  bool* isFlat,
0554                                                                  unsigned int nPoints = 5,
0555                                                                  bool anchorHits = true) {
0556     /*
0557     Bool anchorHits required to deal with a weird edge case wherein 
0558     the hits ultimately used in the regression are anchor hits, but the
0559     lower modules need not all be Pixel Modules (in case of PS). Similarly,
0560     when we compute the chi squared for the non-anchor hits, the "partner module"
0561     need not always be a PS strip module, but all non-anchor hits sit on strip 
0562     modules.
0563     */
0564 
0565     ModuleType moduleType;
0566     short moduleSubdet, moduleSide;
0567     float inv1 = kWidthPS / kWidth2S;
0568     float inv2 = kPixelPSZpitch / kWidth2S;
0569     float inv3 = kStripPSZpitch / kWidth2S;
0570     for (size_t i = 0; i < nPoints; i++) {
0571       moduleType = modules.moduleType()[lowerModuleIndices[i]];
0572       moduleSubdet = modules.subdets()[lowerModuleIndices[i]];
0573       moduleSide = modules.sides()[lowerModuleIndices[i]];
0574       const float& drdz = modules.drdzs()[lowerModuleIndices[i]];
0575       slopes[i] = modules.dxdys()[lowerModuleIndices[i]];
0576       //category 1 - barrel PS flat
0577       if (moduleSubdet == Barrel and moduleType == PS and moduleSide == Center) {
0578         delta1[i] = inv1;
0579         delta2[i] = inv1;
0580         slopes[i] = -999.f;
0581         isFlat[i] = true;
0582       }
0583       //category 2 - barrel 2S
0584       else if (moduleSubdet == Barrel and moduleType == TwoS) {
0585         delta1[i] = 1.f;
0586         delta2[i] = 1.f;
0587         slopes[i] = -999.f;
0588         isFlat[i] = true;
0589       }
0590       //category 3 - barrel PS tilted
0591       else if (moduleSubdet == Barrel and moduleType == PS and moduleSide != Center) {
0592         delta1[i] = inv1;
0593         isFlat[i] = false;
0594 
0595         if (anchorHits) {
0596           delta2[i] = (inv2 * drdz / alpaka::math::sqrt(acc, 1 + drdz * drdz));
0597         } else {
0598           delta2[i] = (inv3 * drdz / alpaka::math::sqrt(acc, 1 + drdz * drdz));
0599         }
0600       }
0601       //category 4 - endcap PS
0602       else if (moduleSubdet == Endcap and moduleType == PS) {
0603         delta1[i] = inv1;
0604         isFlat[i] = false;
0605 
0606         /*
0607         despite the type of the module layer of the lower module index,
0608         all anchor hits are on the pixel side and all non-anchor hits are
0609         on the strip side!
0610         */
0611         if (anchorHits) {
0612           delta2[i] = inv2;
0613         } else {
0614           delta2[i] = inv3;
0615         }
0616       }
0617       //category 5 - endcap 2S
0618       else if (moduleSubdet == Endcap and moduleType == TwoS) {
0619         delta1[i] = 1.f;
0620         delta2[i] = 500.f * inv1;
0621         isFlat[i] = false;
0622       } else {
0623 #ifdef WARNINGS
0624         printf("ERROR!!!!! I SHOULDN'T BE HERE!!!! subdet = %d, type = %d, side = %d\n",
0625                moduleSubdet,
0626                moduleType,
0627                moduleSide);
0628 #endif
0629       }
0630     }
0631   }
0632 
0633   template <typename TAcc>
0634   ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeRadiusUsingRegression(TAcc const& acc,
0635                                                                     unsigned int nPoints,
0636                                                                     float* xs,
0637                                                                     float* ys,
0638                                                                     float* delta1,
0639                                                                     float* delta2,
0640                                                                     float* slopes,
0641                                                                     bool* isFlat,
0642                                                                     float& g,
0643                                                                     float& f,
0644                                                                     float* sigmas2,
0645                                                                     float& chiSquared) {
0646     float radius = 0.f;
0647 
0648     // Some extra variables
0649     // the two variables will be called x1 and x2, and y (which is x^2 + y^2)
0650 
0651     float sigmaX1Squared = 0.f;
0652     float sigmaX2Squared = 0.f;
0653     float sigmaX1X2 = 0.f;
0654     float sigmaX1y = 0.f;
0655     float sigmaX2y = 0.f;
0656     float sigmaY = 0.f;
0657     float sigmaX1 = 0.f;
0658     float sigmaX2 = 0.f;
0659     float sigmaOne = 0.f;
0660 
0661     float xPrime, yPrime, absArctanSlope, angleM;
0662     for (size_t i = 0; i < nPoints; i++) {
0663       // Computing sigmas is a very tricky affair
0664       // if the module is tilted or endcap, we need to use the slopes properly!
0665 
0666       absArctanSlope = ((slopes[i] != kVerticalModuleSlope && edm::isFinite(slopes[i]))
0667                             ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i]))
0668                             : kPi / 2.f);
0669 
0670       if (xs[i] > 0 and ys[i] > 0) {
0671         angleM = kPi / 2.f - absArctanSlope;
0672       } else if (xs[i] < 0 and ys[i] > 0) {
0673         angleM = absArctanSlope + kPi / 2.f;
0674       } else if (xs[i] < 0 and ys[i] < 0) {
0675         angleM = -(absArctanSlope + kPi / 2.f);
0676       } else if (xs[i] > 0 and ys[i] < 0) {
0677         angleM = -(kPi / 2.f - absArctanSlope);
0678       } else {
0679         angleM = 0;
0680       }
0681 
0682       if (not isFlat[i]) {
0683         xPrime = xs[i] * alpaka::math::cos(acc, angleM) + ys[i] * alpaka::math::sin(acc, angleM);
0684         yPrime = ys[i] * alpaka::math::cos(acc, angleM) - xs[i] * alpaka::math::sin(acc, angleM);
0685       } else {
0686         xPrime = xs[i];
0687         yPrime = ys[i];
0688       }
0689       sigmas2[i] = 4 * ((xPrime * delta1[i]) * (xPrime * delta1[i]) + (yPrime * delta2[i]) * (yPrime * delta2[i]));
0690 
0691       sigmaX1Squared += (xs[i] * xs[i]) / sigmas2[i];
0692       sigmaX2Squared += (ys[i] * ys[i]) / sigmas2[i];
0693       sigmaX1X2 += (xs[i] * ys[i]) / sigmas2[i];
0694       sigmaX1y += (xs[i] * (xs[i] * xs[i] + ys[i] * ys[i])) / sigmas2[i];
0695       sigmaX2y += (ys[i] * (xs[i] * xs[i] + ys[i] * ys[i])) / sigmas2[i];
0696       sigmaY += (xs[i] * xs[i] + ys[i] * ys[i]) / sigmas2[i];
0697       sigmaX1 += xs[i] / sigmas2[i];
0698       sigmaX2 += ys[i] / sigmas2[i];
0699       sigmaOne += 1.0f / sigmas2[i];
0700     }
0701     float denominator = (sigmaX1X2 - sigmaX1 * sigmaX2) * (sigmaX1X2 - sigmaX1 * sigmaX2) -
0702                         (sigmaX1Squared - sigmaX1 * sigmaX1) * (sigmaX2Squared - sigmaX2 * sigmaX2);
0703 
0704     float twoG = ((sigmaX2y - sigmaX2 * sigmaY) * (sigmaX1X2 - sigmaX1 * sigmaX2) -
0705                   (sigmaX1y - sigmaX1 * sigmaY) * (sigmaX2Squared - sigmaX2 * sigmaX2)) /
0706                  denominator;
0707     float twoF = ((sigmaX1y - sigmaX1 * sigmaY) * (sigmaX1X2 - sigmaX1 * sigmaX2) -
0708                   (sigmaX2y - sigmaX2 * sigmaY) * (sigmaX1Squared - sigmaX1 * sigmaX1)) /
0709                  denominator;
0710 
0711     float c = -(sigmaY - twoG * sigmaX1 - twoF * sigmaX2) / sigmaOne;
0712     g = 0.5f * twoG;
0713     f = 0.5f * twoF;
0714     if (g * g + f * f - c < 0) {
0715 #ifdef WARNINGS
0716       printf("FATAL! r^2 < 0!\n");
0717 #endif
0718       chiSquared = -1;
0719       return -1;
0720     }
0721 
0722     radius = alpaka::math::sqrt(acc, g * g + f * f - c);
0723     // compute chi squared
0724     chiSquared = 0.f;
0725     for (size_t i = 0; i < nPoints; i++) {
0726       chiSquared += (xs[i] * xs[i] + ys[i] * ys[i] - twoG * xs[i] - twoF * ys[i] + c) *
0727                     (xs[i] * xs[i] + ys[i] * ys[i] - twoG * xs[i] - twoF * ys[i] + c) / sigmas2[i];
0728     }
0729     return radius;
0730   }
0731 
0732   template <typename TAcc>
0733   ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeChiSquared(TAcc const& acc,
0734                                                          unsigned int nPoints,
0735                                                          float* xs,
0736                                                          float* ys,
0737                                                          float* delta1,
0738                                                          float* delta2,
0739                                                          float* slopes,
0740                                                          bool* isFlat,
0741                                                          float g,
0742                                                          float f,
0743                                                          float radius) {
0744     // given values of (g, f, radius) and a set of points (and its uncertainties)
0745     // compute chi squared
0746     float c = g * g + f * f - radius * radius;
0747     float chiSquared = 0.f;
0748     float absArctanSlope, angleM, xPrime, yPrime, sigma2;
0749     for (size_t i = 0; i < nPoints; i++) {
0750       absArctanSlope = ((slopes[i] != kVerticalModuleSlope && edm::isFinite(slopes[i]))
0751                             ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i]))
0752                             : kPi / 2.f);
0753       if (xs[i] > 0 and ys[i] > 0) {
0754         angleM = kPi / 2.f - absArctanSlope;
0755       } else if (xs[i] < 0 and ys[i] > 0) {
0756         angleM = absArctanSlope + kPi / 2.f;
0757       } else if (xs[i] < 0 and ys[i] < 0) {
0758         angleM = -(absArctanSlope + kPi / 2.f);
0759       } else if (xs[i] > 0 and ys[i] < 0) {
0760         angleM = -(kPi / 2.f - absArctanSlope);
0761       } else {
0762         angleM = 0;
0763       }
0764 
0765       if (not isFlat[i]) {
0766         xPrime = xs[i] * alpaka::math::cos(acc, angleM) + ys[i] * alpaka::math::sin(acc, angleM);
0767         yPrime = ys[i] * alpaka::math::cos(acc, angleM) - xs[i] * alpaka::math::sin(acc, angleM);
0768       } else {
0769         xPrime = xs[i];
0770         yPrime = ys[i];
0771       }
0772       sigma2 = 4 * ((xPrime * delta1[i]) * (xPrime * delta1[i]) + (yPrime * delta2[i]) * (yPrime * delta2[i]));
0773       chiSquared += (xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) *
0774                     (xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) / sigma2;
0775     }
0776     return chiSquared;
0777   }
0778 
0779   template <typename TAcc>
0780   ALPAKA_FN_ACC ALPAKA_FN_INLINE void runDeltaBetaIterations(TAcc const& acc,
0781                                                              float& betaIn,
0782                                                              float& betaOut,
0783                                                              float betaAv,
0784                                                              float& pt_beta,
0785                                                              float sdIn_dr,
0786                                                              float sdOut_dr,
0787                                                              float dr,
0788                                                              float lIn) {
0789     if (lIn == 0) {
0790       betaOut += alpaka::math::copysign(
0791           acc,
0792           alpaka::math::asin(
0793               acc, alpaka::math::min(acc, sdOut_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), kSinAlphaMax)),
0794           betaOut);
0795       return;
0796     }
0797 
0798     if (betaIn * betaOut > 0.f and
0799         (alpaka::math::abs(acc, pt_beta) < 4.f * kPt_betaMax or
0800          (lIn >= 11 and alpaka::math::abs(acc, pt_beta) <
0801                             8.f * kPt_betaMax)))  //and the pt_beta is well-defined; less strict for endcap-endcap
0802     {
0803       const float betaInUpd =
0804           betaIn +
0805           alpaka::math::copysign(
0806               acc,
0807               alpaka::math::asin(
0808                   acc, alpaka::math::min(acc, sdIn_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), kSinAlphaMax)),
0809               betaIn);  //FIXME: need a faster version
0810       const float betaOutUpd =
0811           betaOut +
0812           alpaka::math::copysign(
0813               acc,
0814               alpaka::math::asin(
0815                   acc, alpaka::math::min(acc, sdOut_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), kSinAlphaMax)),
0816               betaOut);  //FIXME: need a faster version
0817       betaAv = 0.5f * (betaInUpd + betaOutUpd);
0818 
0819       //1st update
0820       const float pt_beta_inv =
0821           1.f / alpaka::math::abs(acc, dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv));  //get a better pt estimate
0822 
0823       betaIn += alpaka::math::copysign(
0824           acc,
0825           alpaka::math::asin(acc, alpaka::math::min(acc, sdIn_dr * k2Rinv1GeVf * pt_beta_inv, kSinAlphaMax)),
0826           betaIn);  //FIXME: need a faster version
0827       betaOut += alpaka::math::copysign(
0828           acc,
0829           alpaka::math::asin(acc, alpaka::math::min(acc, sdOut_dr * k2Rinv1GeVf * pt_beta_inv, kSinAlphaMax)),
0830           betaOut);  //FIXME: need a faster version
0831       //update the av and pt
0832       betaAv = 0.5f * (betaIn + betaOut);
0833       //2nd update
0834       pt_beta = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv);  //get a better pt estimate
0835     } else if (lIn < 11 && alpaka::math::abs(acc, betaOut) < 0.2f * alpaka::math::abs(acc, betaIn) &&
0836                alpaka::math::abs(acc, pt_beta) < 12.f * kPt_betaMax)  //use betaIn sign as ref
0837     {
0838       const float pt_betaIn = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaIn);
0839 
0840       const float betaInUpd =
0841           betaIn +
0842           alpaka::math::copysign(
0843               acc,
0844               alpaka::math::asin(
0845                   acc, alpaka::math::min(acc, sdIn_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_betaIn), kSinAlphaMax)),
0846               betaIn);  //FIXME: need a faster version
0847       const float betaOutUpd =
0848           betaOut +
0849           alpaka::math::copysign(
0850               acc,
0851               alpaka::math::asin(
0852                   acc,
0853                   alpaka::math::min(acc, sdOut_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_betaIn), kSinAlphaMax)),
0854               betaIn);  //FIXME: need a faster version
0855       betaAv = (alpaka::math::abs(acc, betaOut) > 0.2f * alpaka::math::abs(acc, betaIn))
0856                    ? (0.5f * (betaInUpd + betaOutUpd))
0857                    : betaInUpd;
0858 
0859       //1st update
0860       pt_beta = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv);  //get a better pt estimate
0861       betaIn += alpaka::math::copysign(
0862           acc,
0863           alpaka::math::asin(
0864               acc, alpaka::math::min(acc, sdIn_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), kSinAlphaMax)),
0865           betaIn);  //FIXME: need a faster version
0866       betaOut += alpaka::math::copysign(
0867           acc,
0868           alpaka::math::asin(
0869               acc, alpaka::math::min(acc, sdOut_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), kSinAlphaMax)),
0870           betaIn);  //FIXME: need a faster version
0871       //update the av and pt
0872       betaAv = 0.5f * (betaIn + betaOut);
0873       //2nd update
0874       pt_beta = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv);  //get a better pt estimate
0875     }
0876   }
0877 
0878   template <typename TAcc>
0879   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletdBetaCutBBBB(TAcc const& acc,
0880                                                                 ModulesConst modules,
0881                                                                 MiniDoubletsConst mds,
0882                                                                 SegmentsConst segments,
0883                                                                 uint16_t innerInnerLowerModuleIndex,
0884                                                                 uint16_t innerOuterLowerModuleIndex,
0885                                                                 uint16_t outerInnerLowerModuleIndex,
0886                                                                 uint16_t outerOuterLowerModuleIndex,
0887                                                                 unsigned int innerSegmentIndex,
0888                                                                 unsigned int outerSegmentIndex,
0889                                                                 unsigned int firstMDIndex,
0890                                                                 unsigned int secondMDIndex,
0891                                                                 unsigned int thirdMDIndex,
0892                                                                 unsigned int fourthMDIndex,
0893                                                                 float& dBeta,
0894                                                                 const float ptCut) {
0895     float rt_InLo = mds.anchorRt()[firstMDIndex];
0896     float rt_InOut = mds.anchorRt()[secondMDIndex];
0897     float rt_OutLo = mds.anchorRt()[thirdMDIndex];
0898 
0899     float z_InLo = mds.anchorZ()[firstMDIndex];
0900     float z_OutLo = mds.anchorZ()[thirdMDIndex];
0901 
0902     float r3_InLo = alpaka::math::sqrt(acc, z_InLo * z_InLo + rt_InLo * rt_InLo);
0903     float drt_InSeg = rt_InOut - rt_InLo;
0904 
0905     float thetaMuls2 = (kMulsInGeV * kMulsInGeV) * (0.1f + 0.2f * (rt_OutLo - rt_InLo) / 50.f) * (r3_InLo / rt_InLo);
0906 
0907     float midPointX = 0.5f * (mds.anchorX()[firstMDIndex] + mds.anchorX()[thirdMDIndex]);
0908     float midPointY = 0.5f * (mds.anchorY()[firstMDIndex] + mds.anchorY()[thirdMDIndex]);
0909     float diffX = mds.anchorX()[thirdMDIndex] - mds.anchorX()[firstMDIndex];
0910     float diffY = mds.anchorY()[thirdMDIndex] - mds.anchorY()[firstMDIndex];
0911 
0912     float dPhi = cms::alpakatools::deltaPhi(acc, midPointX, midPointY, diffX, diffY);
0913 
0914     // First obtaining the raw betaIn and betaOut values without any correction and just purely based on the mini-doublet hit positions
0915     float alpha_InLo = __H2F(segments.dPhiChanges()[innerSegmentIndex]);
0916     float alpha_OutLo = __H2F(segments.dPhiChanges()[outerSegmentIndex]);
0917 
0918     bool isEC_lastLayer = modules.subdets()[outerOuterLowerModuleIndex] == Endcap and
0919                           modules.moduleType()[outerOuterLowerModuleIndex] == TwoS;
0920 
0921     float alpha_OutUp, alpha_OutUp_highEdge, alpha_OutUp_lowEdge;
0922 
0923     alpha_OutUp = cms::alpakatools::reducePhiRange(
0924         acc,
0925         cms::alpakatools::phi(acc,
0926                               mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex],
0927                               mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]) -
0928             mds.anchorPhi()[fourthMDIndex]);
0929 
0930     alpha_OutUp_highEdge = alpha_OutUp;
0931     alpha_OutUp_lowEdge = alpha_OutUp;
0932 
0933     float tl_axis_x = mds.anchorX()[fourthMDIndex] - mds.anchorX()[firstMDIndex];
0934     float tl_axis_y = mds.anchorY()[fourthMDIndex] - mds.anchorY()[firstMDIndex];
0935     float tl_axis_highEdge_x = tl_axis_x;
0936     float tl_axis_highEdge_y = tl_axis_y;
0937     float tl_axis_lowEdge_x = tl_axis_x;
0938     float tl_axis_lowEdge_y = tl_axis_y;
0939 
0940     float betaIn =
0941         alpha_InLo - cms::alpakatools::reducePhiRange(
0942                          acc, cms::alpakatools::phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[firstMDIndex]);
0943 
0944     float betaInRHmin = betaIn;
0945     float betaInRHmax = betaIn;
0946     float betaOut =
0947         -alpha_OutUp + cms::alpakatools::reducePhiRange(
0948                            acc, cms::alpakatools::phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[fourthMDIndex]);
0949 
0950     float betaOutRHmin = betaOut;
0951     float betaOutRHmax = betaOut;
0952 
0953     if (isEC_lastLayer) {
0954       alpha_OutUp_highEdge = cms::alpakatools::reducePhiRange(
0955           acc,
0956           cms::alpakatools::phi(acc,
0957                                 mds.anchorHighEdgeX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex],
0958                                 mds.anchorHighEdgeY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]) -
0959               mds.anchorHighEdgePhi()[fourthMDIndex]);
0960       alpha_OutUp_lowEdge = cms::alpakatools::reducePhiRange(
0961           acc,
0962           cms::alpakatools::phi(acc,
0963                                 mds.anchorLowEdgeX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex],
0964                                 mds.anchorLowEdgeY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]) -
0965               mds.anchorLowEdgePhi()[fourthMDIndex]);
0966 
0967       tl_axis_highEdge_x = mds.anchorHighEdgeX()[fourthMDIndex] - mds.anchorX()[firstMDIndex];
0968       tl_axis_highEdge_y = mds.anchorHighEdgeY()[fourthMDIndex] - mds.anchorY()[firstMDIndex];
0969       tl_axis_lowEdge_x = mds.anchorLowEdgeX()[fourthMDIndex] - mds.anchorX()[firstMDIndex];
0970       tl_axis_lowEdge_y = mds.anchorLowEdgeY()[fourthMDIndex] - mds.anchorY()[firstMDIndex];
0971 
0972       betaOutRHmin = -alpha_OutUp_highEdge + cms::alpakatools::reducePhiRange(
0973                                                  acc,
0974                                                  cms::alpakatools::phi(acc, tl_axis_highEdge_x, tl_axis_highEdge_y) -
0975                                                      mds.anchorHighEdgePhi()[fourthMDIndex]);
0976       betaOutRHmax = -alpha_OutUp_lowEdge +
0977                      cms::alpakatools::reducePhiRange(acc,
0978                                                       cms::alpakatools::phi(acc, tl_axis_lowEdge_x, tl_axis_lowEdge_y) -
0979                                                           mds.anchorLowEdgePhi()[fourthMDIndex]);
0980     }
0981 
0982     //beta computation
0983     float drt_tl_axis = alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y);
0984 
0985     //innerOuterAnchor - innerInnerAnchor
0986     const float rt_InSeg = alpaka::math::sqrt(acc,
0987                                               (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) *
0988                                                       (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) +
0989                                                   (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]) *
0990                                                       (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]));
0991 
0992     float betaAv = 0.5f * (betaIn + betaOut);
0993     float pt_beta = drt_tl_axis * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv);
0994     int lIn = 5;
0995     int lOut = isEC_lastLayer ? 11 : 5;
0996     float sdOut_dr = alpaka::math::sqrt(acc,
0997                                         (mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex]) *
0998                                                 (mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex]) +
0999                                             (mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]) *
1000                                                 (mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]));
1001     float sdOut_d = mds.anchorRt()[fourthMDIndex] - mds.anchorRt()[thirdMDIndex];
1002 
1003     runDeltaBetaIterations(acc, betaIn, betaOut, betaAv, pt_beta, rt_InSeg, sdOut_dr, drt_tl_axis, lIn);
1004 
1005     const float betaInMMSF = (alpaka::math::abs(acc, betaInRHmin + betaInRHmax) > 0)
1006                                  ? (2.f * betaIn / alpaka::math::abs(acc, betaInRHmin + betaInRHmax))
1007                                  : 0.f;  //mean value of min,max is the old betaIn
1008     const float betaOutMMSF = (alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax) > 0)
1009                                   ? (2.f * betaOut / alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax))
1010                                   : 0.f;
1011     betaInRHmin *= betaInMMSF;
1012     betaInRHmax *= betaInMMSF;
1013     betaOutRHmin *= betaOutMMSF;
1014     betaOutRHmax *= betaOutMMSF;
1015 
1016     float min_ptBeta_maxPtBeta = alpaka::math::min(
1017         acc, alpaka::math::abs(acc, pt_beta), kPt_betaMax);  //need to confimm the range-out value of 7 GeV
1018     const float dBetaMuls2 = thetaMuls2 * 16.f / (min_ptBeta_maxPtBeta * min_ptBeta_maxPtBeta);
1019 
1020     const float alphaInAbsReg =
1021         alpaka::math::max(acc,
1022                           alpaka::math::abs(acc, alpha_InLo),
1023                           alpaka::math::asin(acc, alpaka::math::min(acc, rt_InLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1024     const float alphaOutAbsReg =
1025         alpaka::math::max(acc,
1026                           alpaka::math::abs(acc, alpha_OutLo),
1027                           alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1028     const float dBetaInLum = lIn < 11 ? 0.0f : alpaka::math::abs(acc, alphaInAbsReg * kDeltaZLum / z_InLo);
1029     const float dBetaOutLum = lOut < 11 ? 0.0f : alpaka::math::abs(acc, alphaOutAbsReg * kDeltaZLum / z_OutLo);
1030     const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum);
1031     const float sinDPhi = alpaka::math::sin(acc, dPhi);
1032 
1033     float dBetaROut = 0;
1034     if (isEC_lastLayer) {
1035       dBetaROut = (alpaka::math::sqrt(acc,
1036                                       mds.anchorHighEdgeX()[fourthMDIndex] * mds.anchorHighEdgeX()[fourthMDIndex] +
1037                                           mds.anchorHighEdgeY()[fourthMDIndex] * mds.anchorHighEdgeY()[fourthMDIndex]) -
1038                    alpaka::math::sqrt(acc,
1039                                       mds.anchorLowEdgeX()[fourthMDIndex] * mds.anchorLowEdgeX()[fourthMDIndex] +
1040                                           mds.anchorLowEdgeY()[fourthMDIndex] * mds.anchorLowEdgeY()[fourthMDIndex])) *
1041                   sinDPhi / drt_tl_axis;
1042     }
1043 
1044     const float dBetaROut2 = dBetaROut * dBetaROut;
1045 
1046     float dBetaRes = 0.02f / alpaka::math::min(acc, sdOut_d, drt_InSeg);
1047     float dBetaCut2 =
1048         (dBetaRes * dBetaRes * 2.0f + dBetaMuls2 + dBetaLum2 + dBetaROut2 +
1049          0.25f *
1050              (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)) *
1051              (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)));
1052 
1053     dBeta = betaIn - betaOut;
1054     return dBeta * dBeta <= dBetaCut2;
1055   }
1056 
1057   template <typename TAcc>
1058   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletdBetaCutBBEE(TAcc const& acc,
1059                                                                 ModulesConst modules,
1060                                                                 MiniDoubletsConst mds,
1061                                                                 SegmentsConst segments,
1062                                                                 uint16_t innerInnerLowerModuleIndex,
1063                                                                 uint16_t innerOuterLowerModuleIndex,
1064                                                                 uint16_t outerInnerLowerModuleIndex,
1065                                                                 uint16_t outerOuterLowerModuleIndex,
1066                                                                 unsigned int innerSegmentIndex,
1067                                                                 unsigned int outerSegmentIndex,
1068                                                                 unsigned int firstMDIndex,
1069                                                                 unsigned int secondMDIndex,
1070                                                                 unsigned int thirdMDIndex,
1071                                                                 unsigned int fourthMDIndex,
1072                                                                 float& dBeta,
1073                                                                 const float ptCut) {
1074     float rt_InLo = mds.anchorRt()[firstMDIndex];
1075     float rt_InOut = mds.anchorRt()[secondMDIndex];
1076     float rt_OutLo = mds.anchorRt()[thirdMDIndex];
1077 
1078     float z_InLo = mds.anchorZ()[firstMDIndex];
1079     float z_OutLo = mds.anchorZ()[thirdMDIndex];
1080 
1081     float rIn = alpaka::math::sqrt(acc, z_InLo * z_InLo + rt_InLo * rt_InLo);
1082     const float thetaMuls2 = (kMulsInGeV * kMulsInGeV) * (0.1f + 0.2f * (rt_OutLo - rt_InLo) / 50.f) * (rIn / rt_InLo);
1083 
1084     float midPointX = 0.5f * (mds.anchorX()[firstMDIndex] + mds.anchorX()[thirdMDIndex]);
1085     float midPointY = 0.5f * (mds.anchorY()[firstMDIndex] + mds.anchorY()[thirdMDIndex]);
1086     float diffX = mds.anchorX()[thirdMDIndex] - mds.anchorX()[firstMDIndex];
1087     float diffY = mds.anchorY()[thirdMDIndex] - mds.anchorY()[firstMDIndex];
1088 
1089     float dPhi = cms::alpakatools::deltaPhi(acc, midPointX, midPointY, diffX, diffY);
1090 
1091     float sdIn_alpha = __H2F(segments.dPhiChanges()[innerSegmentIndex]);
1092     float sdIn_alpha_min = __H2F(segments.dPhiChangeMins()[innerSegmentIndex]);
1093     float sdIn_alpha_max = __H2F(segments.dPhiChangeMaxs()[innerSegmentIndex]);
1094     float sdOut_alpha = sdIn_alpha;
1095 
1096     float sdOut_dPhiPos =
1097         cms::alpakatools::reducePhiRange(acc, mds.anchorPhi()[fourthMDIndex] - mds.anchorPhi()[thirdMDIndex]);
1098 
1099     float sdOut_dPhiChange = __H2F(segments.dPhiChanges()[outerSegmentIndex]);
1100     float sdOut_dPhiChange_min = __H2F(segments.dPhiChangeMins()[outerSegmentIndex]);
1101     float sdOut_dPhiChange_max = __H2F(segments.dPhiChangeMaxs()[outerSegmentIndex]);
1102 
1103     float sdOut_alphaOutRHmin = cms::alpakatools::reducePhiRange(acc, sdOut_dPhiChange_min - sdOut_dPhiPos);
1104     float sdOut_alphaOutRHmax = cms::alpakatools::reducePhiRange(acc, sdOut_dPhiChange_max - sdOut_dPhiPos);
1105     float sdOut_alphaOut = cms::alpakatools::reducePhiRange(acc, sdOut_dPhiChange - sdOut_dPhiPos);
1106 
1107     float tl_axis_x = mds.anchorX()[fourthMDIndex] - mds.anchorX()[firstMDIndex];
1108     float tl_axis_y = mds.anchorY()[fourthMDIndex] - mds.anchorY()[firstMDIndex];
1109 
1110     float betaIn =
1111         sdIn_alpha - cms::alpakatools::reducePhiRange(
1112                          acc, cms::alpakatools::phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[firstMDIndex]);
1113 
1114     float betaInRHmin = betaIn;
1115     float betaInRHmax = betaIn;
1116     float betaOut =
1117         -sdOut_alphaOut + cms::alpakatools::reducePhiRange(
1118                               acc, cms::alpakatools::phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[fourthMDIndex]);
1119 
1120     float betaOutRHmin = betaOut;
1121     float betaOutRHmax = betaOut;
1122 
1123     bool isEC_secondLayer = (modules.subdets()[innerOuterLowerModuleIndex] == Endcap) and
1124                             (modules.moduleType()[innerOuterLowerModuleIndex] == TwoS);
1125 
1126     if (isEC_secondLayer) {
1127       betaInRHmin = betaIn - sdIn_alpha_min + sdIn_alpha;
1128       betaInRHmax = betaIn - sdIn_alpha_max + sdIn_alpha;
1129     }
1130 
1131     betaOutRHmin = betaOut - sdOut_alphaOutRHmin + sdOut_alphaOut;
1132     betaOutRHmax = betaOut - sdOut_alphaOutRHmax + sdOut_alphaOut;
1133 
1134     float swapTemp;
1135     if (alpaka::math::abs(acc, betaOutRHmin) > alpaka::math::abs(acc, betaOutRHmax)) {
1136       swapTemp = betaOutRHmin;
1137       betaOutRHmin = betaOutRHmax;
1138       betaOutRHmax = swapTemp;
1139     }
1140 
1141     if (alpaka::math::abs(acc, betaInRHmin) > alpaka::math::abs(acc, betaInRHmax)) {
1142       swapTemp = betaInRHmin;
1143       betaInRHmin = betaInRHmax;
1144       betaInRHmax = swapTemp;
1145     }
1146 
1147     float sdIn_dr = alpaka::math::sqrt(acc,
1148                                        (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) *
1149                                                (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) +
1150                                            (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]) *
1151                                                (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]));
1152     float sdIn_d = rt_InOut - rt_InLo;
1153 
1154     float dr = alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y);
1155 
1156     float betaAv = 0.5f * (betaIn + betaOut);
1157     float pt_beta = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv);
1158 
1159     float lIn = 5;
1160     float lOut = 11;
1161 
1162     float sdOut_dr = alpaka::math::sqrt(acc,
1163                                         (mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex]) *
1164                                                 (mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex]) +
1165                                             (mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]) *
1166                                                 (mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]));
1167     float sdOut_d = mds.anchorRt()[fourthMDIndex] - mds.anchorRt()[thirdMDIndex];
1168 
1169     runDeltaBetaIterations(acc, betaIn, betaOut, betaAv, pt_beta, sdIn_dr, sdOut_dr, dr, lIn);
1170 
1171     const float betaInMMSF = (alpaka::math::abs(acc, betaInRHmin + betaInRHmax) > 0)
1172                                  ? (2.f * betaIn / alpaka::math::abs(acc, betaInRHmin + betaInRHmax))
1173                                  : 0.;  //mean value of min,max is the old betaIn
1174     const float betaOutMMSF = (alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax) > 0)
1175                                   ? (2.f * betaOut / alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax))
1176                                   : 0.;
1177     betaInRHmin *= betaInMMSF;
1178     betaInRHmax *= betaInMMSF;
1179     betaOutRHmin *= betaOutMMSF;
1180     betaOutRHmax *= betaOutMMSF;
1181 
1182     float min_ptBeta_maxPtBeta = alpaka::math::min(
1183         acc, alpaka::math::abs(acc, pt_beta), kPt_betaMax);  //need to confirm the range-out value of 7 GeV
1184     const float dBetaMuls2 = thetaMuls2 * 16.f / (min_ptBeta_maxPtBeta * min_ptBeta_maxPtBeta);
1185 
1186     const float alphaInAbsReg =
1187         alpaka::math::max(acc,
1188                           alpaka::math::abs(acc, sdIn_alpha),
1189                           alpaka::math::asin(acc, alpaka::math::min(acc, rt_InLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1190     const float alphaOutAbsReg =
1191         alpaka::math::max(acc,
1192                           alpaka::math::abs(acc, sdOut_alpha),
1193                           alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1194     const float dBetaInLum = lIn < 11 ? 0.0f : alpaka::math::abs(acc, alphaInAbsReg * kDeltaZLum / z_InLo);
1195     const float dBetaOutLum = lOut < 11 ? 0.0f : alpaka::math::abs(acc, alphaOutAbsReg * kDeltaZLum / z_OutLo);
1196     const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum);
1197     const float sinDPhi = alpaka::math::sin(acc, dPhi);
1198 
1199     const float dBetaRIn2 = 0;  // TODO-RH
1200     float dBetaROut = 0;
1201     if (modules.moduleType()[outerOuterLowerModuleIndex] == TwoS) {
1202       dBetaROut = (alpaka::math::sqrt(acc,
1203                                       mds.anchorHighEdgeX()[fourthMDIndex] * mds.anchorHighEdgeX()[fourthMDIndex] +
1204                                           mds.anchorHighEdgeY()[fourthMDIndex] * mds.anchorHighEdgeY()[fourthMDIndex]) -
1205                    alpaka::math::sqrt(acc,
1206                                       mds.anchorLowEdgeX()[fourthMDIndex] * mds.anchorLowEdgeX()[fourthMDIndex] +
1207                                           mds.anchorLowEdgeY()[fourthMDIndex] * mds.anchorLowEdgeY()[fourthMDIndex])) *
1208                   sinDPhi / dr;
1209     }
1210 
1211     const float dBetaROut2 = dBetaROut * dBetaROut;
1212 
1213     float dBetaRes = 0.02f / alpaka::math::min(acc, sdOut_d, sdIn_d);
1214     float dBetaCut2 =
1215         (dBetaRes * dBetaRes * 2.0f + dBetaMuls2 + dBetaLum2 + dBetaRIn2 + dBetaROut2 +
1216          0.25f *
1217              (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)) *
1218              (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)));
1219     dBeta = betaIn - betaOut;
1220     return dBeta * dBeta <= dBetaCut2;
1221   }
1222 
1223   template <typename TAcc>
1224   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletdBetaCutEEEE(TAcc const& acc,
1225                                                                 ModulesConst modules,
1226                                                                 MiniDoubletsConst mds,
1227                                                                 SegmentsConst segments,
1228                                                                 uint16_t innerInnerLowerModuleIndex,
1229                                                                 uint16_t innerOuterLowerModuleIndex,
1230                                                                 uint16_t outerInnerLowerModuleIndex,
1231                                                                 uint16_t outerOuterLowerModuleIndex,
1232                                                                 unsigned int innerSegmentIndex,
1233                                                                 unsigned int outerSegmentIndex,
1234                                                                 unsigned int firstMDIndex,
1235                                                                 unsigned int secondMDIndex,
1236                                                                 unsigned int thirdMDIndex,
1237                                                                 unsigned int fourthMDIndex,
1238                                                                 float& dBeta,
1239                                                                 const float ptCut) {
1240     float rt_InLo = mds.anchorRt()[firstMDIndex];
1241     float rt_InOut = mds.anchorRt()[secondMDIndex];
1242     float rt_OutLo = mds.anchorRt()[thirdMDIndex];
1243 
1244     float z_InLo = mds.anchorZ()[firstMDIndex];
1245     float z_OutLo = mds.anchorZ()[thirdMDIndex];
1246 
1247     float thetaMuls2 = (kMulsInGeV * kMulsInGeV) * (0.1f + 0.2f * (rt_OutLo - rt_InLo) / 50.f);
1248     float sdIn_alpha = __H2F(segments.dPhiChanges()[innerSegmentIndex]);
1249     float sdOut_alpha = sdIn_alpha;  //weird
1250     float sdOut_dPhiPos =
1251         cms::alpakatools::reducePhiRange(acc, mds.anchorPhi()[fourthMDIndex] - mds.anchorPhi()[thirdMDIndex]);
1252 
1253     float sdOut_dPhiChange = __H2F(segments.dPhiChanges()[outerSegmentIndex]);
1254     float sdOut_dPhiChange_min = __H2F(segments.dPhiChangeMins()[outerSegmentIndex]);
1255     float sdOut_dPhiChange_max = __H2F(segments.dPhiChangeMaxs()[outerSegmentIndex]);
1256 
1257     float sdOut_alphaOutRHmin = cms::alpakatools::reducePhiRange(acc, sdOut_dPhiChange_min - sdOut_dPhiPos);
1258     float sdOut_alphaOutRHmax = cms::alpakatools::reducePhiRange(acc, sdOut_dPhiChange_max - sdOut_dPhiPos);
1259     float sdOut_alphaOut = cms::alpakatools::reducePhiRange(acc, sdOut_dPhiChange - sdOut_dPhiPos);
1260 
1261     float tl_axis_x = mds.anchorX()[fourthMDIndex] - mds.anchorX()[firstMDIndex];
1262     float tl_axis_y = mds.anchorY()[fourthMDIndex] - mds.anchorY()[firstMDIndex];
1263 
1264     float betaIn =
1265         sdIn_alpha - cms::alpakatools::reducePhiRange(
1266                          acc, cms::alpakatools::phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[firstMDIndex]);
1267 
1268     float sdIn_alphaRHmin = __H2F(segments.dPhiChangeMins()[innerSegmentIndex]);
1269     float sdIn_alphaRHmax = __H2F(segments.dPhiChangeMaxs()[innerSegmentIndex]);
1270     float betaInRHmin = betaIn + sdIn_alphaRHmin - sdIn_alpha;
1271     float betaInRHmax = betaIn + sdIn_alphaRHmax - sdIn_alpha;
1272 
1273     float betaOut =
1274         -sdOut_alphaOut + cms::alpakatools::reducePhiRange(
1275                               acc, cms::alpakatools::phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[fourthMDIndex]);
1276 
1277     float betaOutRHmin = betaOut - sdOut_alphaOutRHmin + sdOut_alphaOut;
1278     float betaOutRHmax = betaOut - sdOut_alphaOutRHmax + sdOut_alphaOut;
1279 
1280     float swapTemp;
1281     if (alpaka::math::abs(acc, betaOutRHmin) > alpaka::math::abs(acc, betaOutRHmax)) {
1282       swapTemp = betaOutRHmin;
1283       betaOutRHmin = betaOutRHmax;
1284       betaOutRHmax = swapTemp;
1285     }
1286 
1287     if (alpaka::math::abs(acc, betaInRHmin) > alpaka::math::abs(acc, betaInRHmax)) {
1288       swapTemp = betaInRHmin;
1289       betaInRHmin = betaInRHmax;
1290       betaInRHmax = swapTemp;
1291     }
1292     float sdIn_dr = alpaka::math::sqrt(acc,
1293                                        (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) *
1294                                                (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) +
1295                                            (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]) *
1296                                                (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]));
1297     float sdIn_d = rt_InOut - rt_InLo;
1298 
1299     float dr = alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y);
1300 
1301     float betaAv = 0.5f * (betaIn + betaOut);
1302     float pt_beta = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv);
1303 
1304     int lIn = 11;   //endcap
1305     int lOut = 13;  //endcap
1306 
1307     float sdOut_dr = alpaka::math::sqrt(acc,
1308                                         (mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex]) *
1309                                                 (mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex]) +
1310                                             (mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]) *
1311                                                 (mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]));
1312     float sdOut_d = mds.anchorRt()[fourthMDIndex] - mds.anchorRt()[thirdMDIndex];
1313 
1314     runDeltaBetaIterations(acc, betaIn, betaOut, betaAv, pt_beta, sdIn_dr, sdOut_dr, dr, lIn);
1315 
1316     const float betaInMMSF = (alpaka::math::abs(acc, betaInRHmin + betaInRHmax) > 0)
1317                                  ? (2.f * betaIn / alpaka::math::abs(acc, betaInRHmin + betaInRHmax))
1318                                  : 0.;  //mean value of min,max is the old betaIn
1319     const float betaOutMMSF = (alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax) > 0)
1320                                   ? (2.f * betaOut / alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax))
1321                                   : 0.;
1322     betaInRHmin *= betaInMMSF;
1323     betaInRHmax *= betaInMMSF;
1324     betaOutRHmin *= betaOutMMSF;
1325     betaOutRHmax *= betaOutMMSF;
1326 
1327     float min_ptBeta_maxPtBeta = alpaka::math::min(
1328         acc, alpaka::math::abs(acc, pt_beta), kPt_betaMax);  //need to confirm the range-out value of 7 GeV
1329     const float dBetaMuls2 = thetaMuls2 * 16.f / (min_ptBeta_maxPtBeta * min_ptBeta_maxPtBeta);
1330 
1331     const float alphaInAbsReg =
1332         alpaka::math::max(acc,
1333                           alpaka::math::abs(acc, sdIn_alpha),
1334                           alpaka::math::asin(acc, alpaka::math::min(acc, rt_InLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1335     const float alphaOutAbsReg =
1336         alpaka::math::max(acc,
1337                           alpaka::math::abs(acc, sdOut_alpha),
1338                           alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1339     const float dBetaInLum = lIn < 11 ? 0.0f : alpaka::math::abs(acc, alphaInAbsReg * kDeltaZLum / z_InLo);
1340     const float dBetaOutLum = lOut < 11 ? 0.0f : alpaka::math::abs(acc, alphaOutAbsReg * kDeltaZLum / z_OutLo);
1341     const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum);
1342 
1343     float dBetaRes = 0.02f / alpaka::math::min(acc, sdOut_d, sdIn_d);
1344     float dBetaCut2 =
1345         (dBetaRes * dBetaRes * 2.0f + dBetaMuls2 + dBetaLum2 +
1346          0.25f *
1347              (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)) *
1348              (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)));
1349     dBeta = betaIn - betaOut;
1350     return dBeta * dBeta <= dBetaCut2;
1351   }
1352 
1353   template <typename TAcc>
1354   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletdBetaAlgoSelector(TAcc const& acc,
1355                                                                      ModulesConst modules,
1356                                                                      MiniDoubletsConst mds,
1357                                                                      SegmentsConst segments,
1358                                                                      uint16_t innerInnerLowerModuleIndex,
1359                                                                      uint16_t innerOuterLowerModuleIndex,
1360                                                                      uint16_t outerInnerLowerModuleIndex,
1361                                                                      uint16_t outerOuterLowerModuleIndex,
1362                                                                      unsigned int innerSegmentIndex,
1363                                                                      unsigned int outerSegmentIndex,
1364                                                                      unsigned int firstMDIndex,
1365                                                                      unsigned int secondMDIndex,
1366                                                                      unsigned int thirdMDIndex,
1367                                                                      unsigned int fourthMDIndex,
1368                                                                      float& dBeta,
1369                                                                      const float ptCut) {
1370     short innerInnerLowerModuleSubdet = modules.subdets()[innerInnerLowerModuleIndex];
1371     short innerOuterLowerModuleSubdet = modules.subdets()[innerOuterLowerModuleIndex];
1372     short outerInnerLowerModuleSubdet = modules.subdets()[outerInnerLowerModuleIndex];
1373     short outerOuterLowerModuleSubdet = modules.subdets()[outerOuterLowerModuleIndex];
1374 
1375     if (innerInnerLowerModuleSubdet == Barrel and innerOuterLowerModuleSubdet == Barrel and
1376         outerInnerLowerModuleSubdet == Barrel and outerOuterLowerModuleSubdet == Barrel) {
1377       return runQuintupletdBetaCutBBBB(acc,
1378                                        modules,
1379                                        mds,
1380                                        segments,
1381                                        innerInnerLowerModuleIndex,
1382                                        innerOuterLowerModuleIndex,
1383                                        outerInnerLowerModuleIndex,
1384                                        outerOuterLowerModuleIndex,
1385                                        innerSegmentIndex,
1386                                        outerSegmentIndex,
1387                                        firstMDIndex,
1388                                        secondMDIndex,
1389                                        thirdMDIndex,
1390                                        fourthMDIndex,
1391                                        dBeta,
1392                                        ptCut);
1393     } else if (innerInnerLowerModuleSubdet == Barrel and innerOuterLowerModuleSubdet == Barrel and
1394                outerInnerLowerModuleSubdet == Endcap and outerOuterLowerModuleSubdet == Endcap) {
1395       return runQuintupletdBetaCutBBEE(acc,
1396                                        modules,
1397                                        mds,
1398                                        segments,
1399                                        innerInnerLowerModuleIndex,
1400                                        innerOuterLowerModuleIndex,
1401                                        outerInnerLowerModuleIndex,
1402                                        outerOuterLowerModuleIndex,
1403                                        innerSegmentIndex,
1404                                        outerSegmentIndex,
1405                                        firstMDIndex,
1406                                        secondMDIndex,
1407                                        thirdMDIndex,
1408                                        fourthMDIndex,
1409                                        dBeta,
1410                                        ptCut);
1411     } else if (innerInnerLowerModuleSubdet == Barrel and innerOuterLowerModuleSubdet == Barrel and
1412                outerInnerLowerModuleSubdet == Barrel and outerOuterLowerModuleSubdet == Endcap) {
1413       return runQuintupletdBetaCutBBBB(acc,
1414                                        modules,
1415                                        mds,
1416                                        segments,
1417                                        innerInnerLowerModuleIndex,
1418                                        innerOuterLowerModuleIndex,
1419                                        outerInnerLowerModuleIndex,
1420                                        outerOuterLowerModuleIndex,
1421                                        innerSegmentIndex,
1422                                        outerSegmentIndex,
1423                                        firstMDIndex,
1424                                        secondMDIndex,
1425                                        thirdMDIndex,
1426                                        fourthMDIndex,
1427                                        dBeta,
1428                                        ptCut);
1429     } else if (innerInnerLowerModuleSubdet == Barrel and innerOuterLowerModuleSubdet == Endcap and
1430                outerInnerLowerModuleSubdet == Endcap and outerOuterLowerModuleSubdet == Endcap) {
1431       return runQuintupletdBetaCutBBEE(acc,
1432                                        modules,
1433                                        mds,
1434                                        segments,
1435                                        innerInnerLowerModuleIndex,
1436                                        innerOuterLowerModuleIndex,
1437                                        outerInnerLowerModuleIndex,
1438                                        outerOuterLowerModuleIndex,
1439                                        innerSegmentIndex,
1440                                        outerSegmentIndex,
1441                                        firstMDIndex,
1442                                        secondMDIndex,
1443                                        thirdMDIndex,
1444                                        fourthMDIndex,
1445                                        dBeta,
1446                                        ptCut);
1447     } else if (innerInnerLowerModuleSubdet == Endcap and innerOuterLowerModuleSubdet == Endcap and
1448                outerInnerLowerModuleSubdet == Endcap and outerOuterLowerModuleSubdet == Endcap) {
1449       return runQuintupletdBetaCutEEEE(acc,
1450                                        modules,
1451                                        mds,
1452                                        segments,
1453                                        innerInnerLowerModuleIndex,
1454                                        innerOuterLowerModuleIndex,
1455                                        outerInnerLowerModuleIndex,
1456                                        outerOuterLowerModuleIndex,
1457                                        innerSegmentIndex,
1458                                        outerSegmentIndex,
1459                                        firstMDIndex,
1460                                        secondMDIndex,
1461                                        thirdMDIndex,
1462                                        fourthMDIndex,
1463                                        dBeta,
1464                                        ptCut);
1465     }
1466 
1467     return false;
1468   }
1469 
1470   template <typename TAcc>
1471   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletDefaultAlgo(TAcc const& acc,
1472                                                                ModulesConst modules,
1473                                                                MiniDoubletsConst mds,
1474                                                                SegmentsConst segments,
1475                                                                TripletsConst triplets,
1476                                                                uint16_t lowerModuleIndex1,
1477                                                                uint16_t lowerModuleIndex2,
1478                                                                uint16_t lowerModuleIndex3,
1479                                                                uint16_t lowerModuleIndex4,
1480                                                                uint16_t lowerModuleIndex5,
1481                                                                unsigned int innerTripletIndex,
1482                                                                unsigned int outerTripletIndex,
1483                                                                float& innerRadius,
1484                                                                float& outerRadius,
1485                                                                float& bridgeRadius,
1486                                                                float& regressionCenterX,
1487                                                                float& regressionCenterY,
1488                                                                float& regressionRadius,
1489                                                                float& rzChiSquared,
1490                                                                float& chiSquared,
1491                                                                float& nonAnchorChiSquared,
1492                                                                float& dBeta1,
1493                                                                float& dBeta2,
1494                                                                bool& tightCutFlag,
1495                                                                float (&t5Embed)[Params_T5::kEmbed],
1496                                                                const float ptCut) {
1497     unsigned int firstSegmentIndex = triplets.segmentIndices()[innerTripletIndex][0];
1498     unsigned int secondSegmentIndex = triplets.segmentIndices()[innerTripletIndex][1];
1499     unsigned int thirdSegmentIndex = triplets.segmentIndices()[outerTripletIndex][0];
1500     unsigned int fourthSegmentIndex = triplets.segmentIndices()[outerTripletIndex][1];
1501 
1502     unsigned int innerOuterOuterMiniDoubletIndex =
1503         segments.mdIndices()[secondSegmentIndex][1];  //inner triplet outer segment outer MD index
1504     unsigned int outerInnerInnerMiniDoubletIndex =
1505         segments.mdIndices()[thirdSegmentIndex][0];  //outer triplet inner segment inner MD index
1506 
1507     //this cut reduces the number of candidates by a factor of 3, i.e., 2 out of 3 warps can end right here!
1508     if (innerOuterOuterMiniDoubletIndex != outerInnerInnerMiniDoubletIndex)
1509       return false;
1510 
1511     unsigned int firstMDIndex = segments.mdIndices()[firstSegmentIndex][0];
1512     unsigned int secondMDIndex = segments.mdIndices()[secondSegmentIndex][0];
1513     unsigned int thirdMDIndex = segments.mdIndices()[secondSegmentIndex][1];
1514     unsigned int fourthMDIndex = segments.mdIndices()[thirdSegmentIndex][1];
1515     unsigned int fifthMDIndex = segments.mdIndices()[fourthSegmentIndex][1];
1516 
1517     float x1 = mds.anchorX()[firstMDIndex];
1518     float x2 = mds.anchorX()[secondMDIndex];
1519     float x3 = mds.anchorX()[thirdMDIndex];
1520     float x4 = mds.anchorX()[fourthMDIndex];
1521     float x5 = mds.anchorX()[fifthMDIndex];
1522 
1523     float y1 = mds.anchorY()[firstMDIndex];
1524     float y2 = mds.anchorY()[secondMDIndex];
1525     float y3 = mds.anchorY()[thirdMDIndex];
1526     float y4 = mds.anchorY()[fourthMDIndex];
1527     float y5 = mds.anchorY()[fifthMDIndex];
1528 
1529     float g, f;
1530     outerRadius = triplets.radius()[outerTripletIndex];
1531     std::tie(bridgeRadius, g, f) = computeRadiusFromThreeAnchorHits(acc, x2, y2, x3, y3, x4, y4);
1532     innerRadius = triplets.radius()[innerTripletIndex];
1533 
1534     bool inference = lst::t5dnn::runInference(acc,
1535                                               mds,
1536                                               firstMDIndex,
1537                                               secondMDIndex,
1538                                               thirdMDIndex,
1539                                               fourthMDIndex,
1540                                               fifthMDIndex,
1541                                               innerRadius,
1542                                               outerRadius,
1543                                               bridgeRadius);
1544     tightCutFlag = tightCutFlag and inference;  // T5-in-TC cut
1545     if (!inference)                             // T5-building cut
1546       return false;
1547 
1548     if (not runQuintupletdBetaAlgoSelector(acc,
1549                                            modules,
1550                                            mds,
1551                                            segments,
1552                                            lowerModuleIndex1,
1553                                            lowerModuleIndex2,
1554                                            lowerModuleIndex3,
1555                                            lowerModuleIndex4,
1556                                            firstSegmentIndex,
1557                                            thirdSegmentIndex,
1558                                            firstMDIndex,
1559                                            secondMDIndex,
1560                                            thirdMDIndex,
1561                                            fourthMDIndex,
1562                                            dBeta1,
1563                                            ptCut))
1564       return false;
1565 
1566     if (not runQuintupletdBetaAlgoSelector(acc,
1567                                            modules,
1568                                            mds,
1569                                            segments,
1570                                            lowerModuleIndex1,
1571                                            lowerModuleIndex2,
1572                                            lowerModuleIndex4,
1573                                            lowerModuleIndex5,
1574                                            firstSegmentIndex,
1575                                            fourthSegmentIndex,
1576                                            firstMDIndex,
1577                                            secondMDIndex,
1578                                            fourthMDIndex,
1579                                            fifthMDIndex,
1580                                            dBeta2,
1581                                            ptCut))
1582       return false;
1583 
1584     g = triplets.centerX()[innerTripletIndex];
1585     f = triplets.centerY()[innerTripletIndex];
1586 
1587     float inner_pt = 2 * k2Rinv1GeVf * innerRadius;
1588 
1589     if (not passT5RZConstraint(acc,
1590                                modules,
1591                                mds,
1592                                firstMDIndex,
1593                                secondMDIndex,
1594                                thirdMDIndex,
1595                                fourthMDIndex,
1596                                fifthMDIndex,
1597                                lowerModuleIndex1,
1598                                lowerModuleIndex2,
1599                                lowerModuleIndex3,
1600                                lowerModuleIndex4,
1601                                lowerModuleIndex5,
1602                                rzChiSquared,
1603                                inner_pt,
1604                                innerRadius,
1605                                g,
1606                                f,
1607                                tightCutFlag))
1608       return false;
1609 
1610     lst::t5embdnn::runEmbed(acc,
1611                             mds,
1612                             firstMDIndex,
1613                             secondMDIndex,
1614                             thirdMDIndex,
1615                             fourthMDIndex,
1616                             fifthMDIndex,
1617                             innerRadius,
1618                             outerRadius,
1619                             bridgeRadius,
1620                             triplets.fakeScore()[innerTripletIndex],
1621                             triplets.promptScore()[innerTripletIndex],
1622                             triplets.displacedScore()[innerTripletIndex],
1623                             triplets.fakeScore()[outerTripletIndex],
1624                             triplets.promptScore()[outerTripletIndex],
1625                             triplets.displacedScore()[outerTripletIndex],
1626                             t5Embed);
1627 
1628     // 5 categories for sigmas
1629     float sigmas2[5], delta1[5], delta2[5], slopes[5];
1630     bool isFlat[5];
1631 
1632     float xVec[] = {x1, x2, x3, x4, x5};
1633     float yVec[] = {y1, y2, y3, y4, y5};
1634     const uint16_t lowerModuleIndices[] = {
1635         lowerModuleIndex1, lowerModuleIndex2, lowerModuleIndex3, lowerModuleIndex4, lowerModuleIndex5};
1636 
1637     computeSigmasForRegression(acc, modules, lowerModuleIndices, delta1, delta2, slopes, isFlat);
1638     regressionRadius = computeRadiusUsingRegression(acc,
1639                                                     Params_T5::kLayers,
1640                                                     xVec,
1641                                                     yVec,
1642                                                     delta1,
1643                                                     delta2,
1644                                                     slopes,
1645                                                     isFlat,
1646                                                     regressionCenterX,
1647                                                     regressionCenterY,
1648                                                     sigmas2,
1649                                                     chiSquared);
1650 
1651     //compute the other chisquared
1652     //non anchor is always shifted for tilted and endcap!
1653     float nonAnchorDelta1[Params_T5::kLayers], nonAnchorDelta2[Params_T5::kLayers], nonAnchorSlopes[Params_T5::kLayers];
1654     float nonAnchorxs[] = {mds.outerX()[firstMDIndex],
1655                            mds.outerX()[secondMDIndex],
1656                            mds.outerX()[thirdMDIndex],
1657                            mds.outerX()[fourthMDIndex],
1658                            mds.outerX()[fifthMDIndex]};
1659     float nonAnchorys[] = {mds.outerY()[firstMDIndex],
1660                            mds.outerY()[secondMDIndex],
1661                            mds.outerY()[thirdMDIndex],
1662                            mds.outerY()[fourthMDIndex],
1663                            mds.outerY()[fifthMDIndex]};
1664 
1665     computeSigmasForRegression(acc,
1666                                modules,
1667                                lowerModuleIndices,
1668                                nonAnchorDelta1,
1669                                nonAnchorDelta2,
1670                                nonAnchorSlopes,
1671                                isFlat,
1672                                Params_T5::kLayers,
1673                                false);
1674     nonAnchorChiSquared = computeChiSquared(acc,
1675                                             Params_T5::kLayers,
1676                                             nonAnchorxs,
1677                                             nonAnchorys,
1678                                             nonAnchorDelta1,
1679                                             nonAnchorDelta2,
1680                                             nonAnchorSlopes,
1681                                             isFlat,
1682                                             regressionCenterX,
1683                                             regressionCenterY,
1684                                             regressionRadius);
1685     return true;
1686   }
1687 
1688   struct CreateQuintuplets {
1689     ALPAKA_FN_ACC void operator()(Acc3D const& acc,
1690                                   ModulesConst modules,
1691                                   MiniDoubletsConst mds,
1692                                   SegmentsConst segments,
1693                                   Triplets triplets,
1694                                   TripletsOccupancyConst tripletsOccupancy,
1695                                   Quintuplets quintuplets,
1696                                   QuintupletsOccupancy quintupletsOccupancy,
1697                                   ObjectRangesConst ranges,
1698                                   uint16_t nEligibleT5Modules,
1699                                   const float ptCut) const {
1700       for (int iter : cms::alpakatools::uniform_elements_z(acc, nEligibleT5Modules)) {
1701         uint16_t lowerModule1 = ranges.indicesOfEligibleT5Modules()[iter];
1702         short layer2_adjustment;
1703         int layer = modules.layers()[lowerModule1];
1704         if (layer == 1) {
1705           layer2_adjustment = 1;
1706         }  // get upper segment to be in second layer
1707         else if (layer == 2) {
1708           layer2_adjustment = 0;
1709         }  // get lower segment to be in second layer
1710         else {
1711           continue;
1712         }
1713         unsigned int nInnerTriplets = tripletsOccupancy.nTriplets()[lowerModule1];
1714         for (unsigned int innerTripletArrayIndex : cms::alpakatools::uniform_elements_y(acc, nInnerTriplets)) {
1715           unsigned int innerTripletIndex = ranges.tripletModuleIndices()[lowerModule1] + innerTripletArrayIndex;
1716           uint16_t lowerModule2 = triplets.lowerModuleIndices()[innerTripletIndex][1];
1717           uint16_t lowerModule3 = triplets.lowerModuleIndices()[innerTripletIndex][2];
1718           unsigned int nOuterTriplets = tripletsOccupancy.nTriplets()[lowerModule3];
1719           for (unsigned int outerTripletArrayIndex : cms::alpakatools::uniform_elements_x(acc, nOuterTriplets)) {
1720             unsigned int outerTripletIndex = ranges.tripletModuleIndices()[lowerModule3] + outerTripletArrayIndex;
1721             uint16_t lowerModule4 = triplets.lowerModuleIndices()[outerTripletIndex][1];
1722             uint16_t lowerModule5 = triplets.lowerModuleIndices()[outerTripletIndex][2];
1723 
1724             float innerRadius, outerRadius, bridgeRadius, regressionCenterX, regressionCenterY, regressionRadius,
1725                 rzChiSquared, chiSquared, nonAnchorChiSquared, dBeta1, dBeta2;  //required for making distributions
1726 
1727             float t5Embed[Params_T5::kEmbed] = {0.f};
1728 
1729             bool tightCutFlag = false;
1730             bool success = runQuintupletDefaultAlgo(acc,
1731                                                     modules,
1732                                                     mds,
1733                                                     segments,
1734                                                     triplets,
1735                                                     lowerModule1,
1736                                                     lowerModule2,
1737                                                     lowerModule3,
1738                                                     lowerModule4,
1739                                                     lowerModule5,
1740                                                     innerTripletIndex,
1741                                                     outerTripletIndex,
1742                                                     innerRadius,
1743                                                     outerRadius,
1744                                                     bridgeRadius,
1745                                                     regressionCenterX,
1746                                                     regressionCenterY,
1747                                                     regressionRadius,
1748                                                     rzChiSquared,
1749                                                     chiSquared,
1750                                                     nonAnchorChiSquared,
1751                                                     dBeta1,
1752                                                     dBeta2,
1753                                                     tightCutFlag,
1754                                                     t5Embed,
1755                                                     ptCut);
1756 
1757             if (success) {
1758               int totOccupancyQuintuplets = alpaka::atomicAdd(
1759                   acc, &quintupletsOccupancy.totOccupancyQuintuplets()[lowerModule1], 1u, alpaka::hierarchy::Threads{});
1760               if (totOccupancyQuintuplets >= ranges.quintupletModuleOccupancy()[lowerModule1]) {
1761 #ifdef WARNINGS
1762                 printf("Quintuplet excess alert! Module index = %d, Occupancy = %d\n",
1763                        lowerModule1,
1764                        totOccupancyQuintuplets);
1765 #endif
1766               } else {
1767                 int quintupletModuleIndex = alpaka::atomicAdd(
1768                     acc, &quintupletsOccupancy.nQuintuplets()[lowerModule1], 1u, alpaka::hierarchy::Threads{});
1769                 //this if statement should never get executed!
1770                 if (ranges.quintupletModuleIndices()[lowerModule1] == -1) {
1771 #ifdef WARNINGS
1772                   printf("Quintuplets : no memory for module at module index = %d\n", lowerModule1);
1773 #endif
1774                 } else {
1775                   unsigned int quintupletIndex = ranges.quintupletModuleIndices()[lowerModule1] + quintupletModuleIndex;
1776                   float phi = mds.anchorPhi()[segments.mdIndices()[triplets.segmentIndices()[innerTripletIndex][0]]
1777                                                                   [layer2_adjustment]];
1778                   float eta = mds.anchorEta()[segments.mdIndices()[triplets.segmentIndices()[innerTripletIndex][0]]
1779                                                                   [layer2_adjustment]];
1780                   float pt = (innerRadius + outerRadius) * k2Rinv1GeVf;
1781                   float scores = chiSquared + nonAnchorChiSquared;
1782                   addQuintupletToMemory(triplets,
1783                                         quintuplets,
1784                                         innerTripletIndex,
1785                                         outerTripletIndex,
1786                                         lowerModule1,
1787                                         lowerModule2,
1788                                         lowerModule3,
1789                                         lowerModule4,
1790                                         lowerModule5,
1791                                         innerRadius,
1792                                         bridgeRadius,
1793                                         outerRadius,
1794                                         regressionCenterX,
1795                                         regressionCenterY,
1796                                         regressionRadius,
1797                                         rzChiSquared,
1798                                         chiSquared,
1799                                         nonAnchorChiSquared,
1800                                         dBeta1,
1801                                         dBeta2,
1802                                         pt,
1803                                         eta,
1804                                         phi,
1805                                         scores,
1806                                         layer,
1807                                         quintupletIndex,
1808                                         t5Embed,
1809                                         tightCutFlag);
1810 
1811                   triplets.partOfT5()[quintuplets.tripletIndices()[quintupletIndex][0]] = true;
1812                   triplets.partOfT5()[quintuplets.tripletIndices()[quintupletIndex][1]] = true;
1813                 }
1814               }
1815             }
1816           }
1817         }
1818       }
1819     }
1820   };
1821 
1822   struct CreateEligibleModulesListForQuintuplets {
1823     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
1824                                   ModulesConst modules,
1825                                   TripletsOccupancyConst tripletsOccupancy,
1826                                   ObjectRanges ranges,
1827                                   Triplets triplets,
1828                                   const float ptCut) const {
1829       // implementation is 1D with a single block
1830       ALPAKA_ASSERT_ACC((alpaka::getWorkDiv<alpaka::Grid, alpaka::Blocks>(acc)[0] == 1));
1831 
1832       // Initialize variables in shared memory and set to 0
1833       int& nEligibleT5Modulesx = alpaka::declareSharedVar<int, __COUNTER__>(acc);
1834       int& nTotalQuintupletsx = alpaka::declareSharedVar<int, __COUNTER__>(acc);
1835       if (cms::alpakatools::once_per_block(acc)) {
1836         nTotalQuintupletsx = 0;
1837         nEligibleT5Modulesx = 0;
1838       }
1839       alpaka::syncBlockThreads(acc);
1840 
1841       // Occupancy matrix for 0.8 GeV pT Cut
1842       constexpr int p08_occupancy_matrix[4][4] = {
1843           {336, 414, 231, 146},  // category 0
1844           {0, 0, 0, 0},          // category 1
1845           {0, 0, 0, 0},          // category 2
1846           {0, 0, 191, 106}       // category 3
1847       };
1848 
1849       // Occupancy matrix for 0.6 GeV pT Cut, 99.99%
1850       constexpr int p06_occupancy_matrix[4][4] = {
1851           {325, 237, 217, 176},  // category 0
1852           {0, 0, 0, 0},          // category 1
1853           {0, 0, 0, 0},          // category 2
1854           {0, 0, 129, 180}       // category 3
1855       };
1856 
1857       // Select the appropriate occupancy matrix based on ptCut
1858       const auto& occupancy_matrix = (ptCut < 0.8f) ? p06_occupancy_matrix : p08_occupancy_matrix;
1859 
1860       for (int i : cms::alpakatools::uniform_elements(acc, modules.nLowerModules())) {
1861         // Condition for a quintuple to exist for a module
1862         // TCs don't exist for layers 5 and 6 barrel, and layers 2,3,4,5 endcap
1863         short module_rings = modules.rings()[i];
1864         short module_layers = modules.layers()[i];
1865         short module_subdets = modules.subdets()[i];
1866         float module_eta = alpaka::math::abs(acc, modules.eta()[i]);
1867 
1868         if (tripletsOccupancy.nTriplets()[i] == 0)
1869           continue;
1870         if (module_subdets == Barrel && module_layers >= 3)
1871           continue;
1872         if (module_subdets == Endcap && module_layers > 1)
1873           continue;
1874 
1875         int dynamic_count = 0;
1876 
1877         // How many triplets are in module i?
1878         int nTriplets_i = tripletsOccupancy.nTriplets()[i];
1879         int firstTripletIdx = ranges.tripletModuleIndices()[i];
1880 
1881         // Loop over all triplets that live in module i
1882         for (int t = 0; t < nTriplets_i; t++) {
1883           int tripletIndex = firstTripletIdx + t;
1884           uint16_t outerModule = triplets.lowerModuleIndices()[tripletIndex][2];
1885           dynamic_count += tripletsOccupancy.nTriplets()[outerModule];
1886         }
1887 
1888         int category_number = getCategoryNumber(module_layers, module_subdets, module_rings);
1889         int eta_number = getEtaBin(module_eta);
1890 
1891 #ifdef WARNINGS
1892         if (category_number == -1 || eta_number == -1) {
1893           printf("Unhandled case in createEligibleModulesListForQuintupletsGPU! Module index = %i\n", i);
1894         }
1895 #endif
1896 
1897         // Get matrix-based cap (use dynamic_count as fallback)
1898         int matrix_cap =
1899             (category_number != -1 && eta_number != -1) ? occupancy_matrix[category_number][eta_number] : 0;
1900 
1901         // Cap occupancy at minimum of dynamic count and matrix value
1902         int occupancy = alpaka::math::min(acc, dynamic_count, matrix_cap);
1903 
1904         int nEligibleT5Modules = alpaka::atomicAdd(acc, &nEligibleT5Modulesx, 1, alpaka::hierarchy::Threads{});
1905         int nTotQ = alpaka::atomicAdd(acc, &nTotalQuintupletsx, occupancy, alpaka::hierarchy::Threads{});
1906 
1907         ranges.quintupletModuleIndices()[i] = nTotQ;
1908         ranges.indicesOfEligibleT5Modules()[nEligibleT5Modules] = i;
1909         ranges.quintupletModuleOccupancy()[i] = occupancy;
1910       }
1911 
1912       // Wait for all threads to finish before reporting final values
1913       alpaka::syncBlockThreads(acc);
1914       if (cms::alpakatools::once_per_block(acc)) {
1915         ranges.nEligibleT5Modules() = static_cast<uint16_t>(nEligibleT5Modulesx);
1916         ranges.nTotalQuints() = static_cast<unsigned int>(nTotalQuintupletsx);
1917       }
1918     }
1919   };
1920 
1921   struct AddQuintupletRangesToEventExplicit {
1922     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
1923                                   ModulesConst modules,
1924                                   QuintupletsOccupancyConst quintupletsOccupancy,
1925                                   ObjectRanges ranges) const {
1926       // implementation is 1D with a single block
1927       ALPAKA_ASSERT_ACC((alpaka::getWorkDiv<alpaka::Grid, alpaka::Blocks>(acc)[0] == 1));
1928 
1929       for (uint16_t i : cms::alpakatools::uniform_elements(acc, modules.nLowerModules())) {
1930         if (quintupletsOccupancy.nQuintuplets()[i] == 0 or ranges.quintupletModuleIndices()[i] == -1) {
1931           ranges.quintupletRanges()[i][0] = -1;
1932           ranges.quintupletRanges()[i][1] = -1;
1933         } else {
1934           ranges.quintupletRanges()[i][0] = ranges.quintupletModuleIndices()[i];
1935           ranges.quintupletRanges()[i][1] =
1936               ranges.quintupletModuleIndices()[i] + quintupletsOccupancy.nQuintuplets()[i] - 1;
1937         }
1938       }
1939     }
1940   };
1941 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE::lst
1942 #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!