Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-12-22 23:30:01 

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