Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 #ifndef RecoTracker_LSTCore_src_alpaka_Triplet_h
0002 #define RecoTracker_LSTCore_src_alpaka_Triplet_h
0003 
0004 #include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
0005 #include "FWCore/Utilities/interface/isFinite.h"
0006 
0007 #include "RecoTracker/LSTCore/interface/alpaka/Common.h"
0008 #include "RecoTracker/LSTCore/interface/ModulesSoA.h"
0009 #include "RecoTracker/LSTCore/interface/ObjectRangesSoA.h"
0010 #include "RecoTracker/LSTCore/interface/TripletsSoA.h"
0011 
0012 #include "Segment.h"
0013 #include "MiniDoublet.h"
0014 #include "Hit.h"
0015 
0016 namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
0017 
0018   ALPAKA_FN_ACC ALPAKA_FN_INLINE void addTripletToMemory(ModulesConst modules,
0019                                                          MiniDoubletsConst mds,
0020                                                          SegmentsConst segments,
0021                                                          Triplets& triplets,
0022                                                          unsigned int innerSegmentIndex,
0023                                                          unsigned int outerSegmentIndex,
0024                                                          uint16_t innerInnerLowerModuleIndex,
0025                                                          uint16_t middleLowerModuleIndex,
0026                                                          uint16_t outerOuterLowerModuleIndex,
0027 #ifdef CUT_VALUE_DEBUG
0028                                                          float zOut,
0029                                                          float rtOut,
0030 #endif
0031                                                          float betaIn,
0032                                                          float betaInCut,
0033                                                          float circleRadius,
0034                                                          float circleCenterX,
0035                                                          float circleCenterY,
0036                                                          unsigned int tripletIndex) {
0037     triplets.segmentIndices()[tripletIndex][0] = innerSegmentIndex;
0038     triplets.segmentIndices()[tripletIndex][1] = outerSegmentIndex;
0039     triplets.lowerModuleIndices()[tripletIndex][0] = innerInnerLowerModuleIndex;
0040     triplets.lowerModuleIndices()[tripletIndex][1] = middleLowerModuleIndex;
0041     triplets.lowerModuleIndices()[tripletIndex][2] = outerOuterLowerModuleIndex;
0042 
0043     triplets.betaIn()[tripletIndex] = __F2H(betaIn);
0044     triplets.radius()[tripletIndex] = circleRadius;
0045     triplets.centerX()[tripletIndex] = circleCenterX;
0046     triplets.centerY()[tripletIndex] = circleCenterY;
0047     triplets.logicalLayers()[tripletIndex][0] =
0048         modules.layers()[innerInnerLowerModuleIndex] + (modules.subdets()[innerInnerLowerModuleIndex] == 4) * 6;
0049     triplets.logicalLayers()[tripletIndex][1] =
0050         modules.layers()[middleLowerModuleIndex] + (modules.subdets()[middleLowerModuleIndex] == 4) * 6;
0051     triplets.logicalLayers()[tripletIndex][2] =
0052         modules.layers()[outerOuterLowerModuleIndex] + (modules.subdets()[outerOuterLowerModuleIndex] == 4) * 6;
0053     //get the hits
0054     unsigned int firstMDIndex = segments.mdIndices()[innerSegmentIndex][0];
0055     unsigned int secondMDIndex = segments.mdIndices()[innerSegmentIndex][1];
0056     unsigned int thirdMDIndex = segments.mdIndices()[outerSegmentIndex][1];
0057 
0058     triplets.hitIndices()[tripletIndex][0] = mds.anchorHitIndices()[firstMDIndex];
0059     triplets.hitIndices()[tripletIndex][1] = mds.outerHitIndices()[firstMDIndex];
0060     triplets.hitIndices()[tripletIndex][2] = mds.anchorHitIndices()[secondMDIndex];
0061     triplets.hitIndices()[tripletIndex][3] = mds.outerHitIndices()[secondMDIndex];
0062     triplets.hitIndices()[tripletIndex][4] = mds.anchorHitIndices()[thirdMDIndex];
0063     triplets.hitIndices()[tripletIndex][5] = mds.outerHitIndices()[thirdMDIndex];
0064 #ifdef CUT_VALUE_DEBUG
0065     triplets.zOut()[tripletIndex] = zOut;
0066     triplets.rtOut()[tripletIndex] = rtOut;
0067     triplets.betaInCut()[tripletIndex] = betaInCut;
0068 #endif
0069   }
0070 
0071   template <typename TAcc>
0072   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passRZConstraint(TAcc const& acc,
0073                                                        ModulesConst modules,
0074                                                        MiniDoubletsConst mds,
0075                                                        uint16_t innerInnerLowerModuleIndex,
0076                                                        uint16_t middleLowerModuleIndex,
0077                                                        uint16_t outerOuterLowerModuleIndex,
0078                                                        unsigned int firstMDIndex,
0079                                                        unsigned int secondMDIndex,
0080                                                        unsigned int thirdMDIndex,
0081                                                        float circleRadius,
0082                                                        float circleCenterX,
0083                                                        float circleCenterY) {
0084     // Using lst_layer numbering convention defined in ModuleMethods.h
0085     const int layer1 = modules.lstLayers()[innerInnerLowerModuleIndex];
0086     const int layer2 = modules.lstLayers()[middleLowerModuleIndex];
0087     const int layer3 = modules.lstLayers()[outerOuterLowerModuleIndex];
0088 
0089     //all the values are stored in the unit of cm, in the calculation below we need to be cautious if we want to use the meter unit
0090     //get r and z
0091     const float r1 = mds.anchorRt()[firstMDIndex] / 100;
0092     const float r2 = mds.anchorRt()[secondMDIndex] / 100;
0093     const float r3 = mds.anchorRt()[thirdMDIndex] / 100;
0094 
0095     const float z1 = mds.anchorZ()[firstMDIndex] / 100;
0096     const float z2 = mds.anchorZ()[secondMDIndex] / 100;
0097     const float z3 = mds.anchorZ()[thirdMDIndex] / 100;
0098 
0099     //use linear approximation for regions 9 and 20-24 because it works better (see https://github.com/SegmentLinking/cmssw/pull/92)
0100     float residual = alpaka::math::abs(acc, z2 - ((z3 - z1) / (r3 - r1) * (r2 - r1) + z1));
0101 
0102     //region definitions: https://github.com/user-attachments/assets/2b3c1425-66eb-4524-83de-deb6f3b31f71
0103     if (layer1 == 1 && layer2 == 7) {
0104       return residual < 0.01f;  // Region 9
0105     } else if (layer1 == 3 && layer2 == 4) {
0106       if (layer3 == 5) {
0107         return residual < 0.037127972f;  // Region 20
0108       } else if (layer3 == 12) {
0109         return residual < 0.05f;  // Region 21
0110       }
0111     } else if (layer1 == 4) {
0112       if (layer2 == 12) {
0113         return residual < 0.063831687f;  // Region 22
0114       } else if (layer2 == 5) {
0115         if (layer3 == 6) {
0116           return residual < 0.04362525f;  // Region 23
0117         } else if (layer3 == 12) {
0118           return residual < 0.05f;  // Region 24
0119         }
0120       }
0121     }
0122 
0123     //get the type of module: 0 is ps, 1 is 2s
0124     const int moduleType3 = modules.moduleType()[outerOuterLowerModuleIndex];
0125 
0126     //get the x,y position of each MD
0127     const float x1 = mds.anchorX()[firstMDIndex] / 100;
0128     const float x2 = mds.anchorX()[secondMDIndex] / 100;
0129     const float x3 = mds.anchorX()[thirdMDIndex] / 100;
0130 
0131     const float y1 = mds.anchorY()[firstMDIndex] / 100;
0132     const float y2 = mds.anchorY()[secondMDIndex] / 100;
0133     const float y3 = mds.anchorY()[thirdMDIndex] / 100;
0134 
0135     //set initial and target points
0136     float x_init = x2;
0137     float y_init = y2;
0138     float z_init = z2;
0139     float r_init = r2;
0140 
0141     float z_target = z3;
0142     float r_target = r3;
0143 
0144     float x_other = x1;
0145     float y_other = y1;
0146 
0147     float dz = z2 - z1;
0148 
0149     //use MD2 for regions 5 and 19 because it works better (see https://github.com/SegmentLinking/cmssw/pull/92)
0150     if ((layer1 == 8 && layer2 == 14 && layer3 == 15) || (layer1 == 3 && layer2 == 12 && layer3 == 13)) {
0151       x_init = x1;
0152       y_init = y1;
0153       z_init = z1;
0154       r_init = r1;
0155 
0156       z_target = z2;
0157       r_target = r2;
0158 
0159       x_other = x3;
0160       y_other = y3;
0161 
0162       dz = z3 - z1;
0163     }
0164 
0165     //use the 3 MDs to fit a circle. This is the circle parameters, for circle centers and circle radius
0166     float x_center = circleCenterX / 100;
0167     float y_center = circleCenterY / 100;
0168     float pt = 2 * k2Rinv1GeVf * circleRadius;  //k2Rinv1GeVf is already in cm^(-1)
0169 
0170     //determine the charge
0171     int charge = 0;
0172     if ((x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1) > 0)
0173       charge = -1;
0174     else
0175       charge = 1;
0176 
0177     //get the absolute value of px and py at the initial point
0178     float px = 2 * k2Rinv1GeVf * alpaka::math::abs(acc, (y_init - y_center)) * 100;
0179     float py = 2 * k2Rinv1GeVf * alpaka::math::abs(acc, (x_init - x_center)) * 100;
0180 
0181     //Above line only gives you the correct value of px and py, but signs of px and py calculated below.
0182     //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.
0183     if (x_init > x_center && y_init > y_center)  //1st quad
0184     {
0185       if (charge == 1)
0186         py = -py;
0187       if (charge == -1)
0188         px = -px;
0189     }
0190     if (x_init < x_center && y_init > y_center)  //2nd quad
0191     {
0192       if (charge == -1) {
0193         px = -px;
0194         py = -py;
0195       }
0196     }
0197     if (x_init < x_center && y_init < y_center)  //3rd quad
0198     {
0199       if (charge == 1)
0200         px = -px;
0201       if (charge == -1)
0202         py = -py;
0203     }
0204     if (x_init > x_center && y_init < y_center)  //4th quad
0205     {
0206       if (charge == 1) {
0207         px = -px;
0208         py = -py;
0209       }
0210     }
0211 
0212     //But if the initial T3 curve goes across quarters(i.e. cross axis to separate the quarters), need special redeclaration of px,py signs on these to avoid errors
0213     if (x3 < x2 && x2 < x1)
0214       px = -alpaka::math::abs(acc, px);
0215     else if (x3 > x2 && x2 > x1)
0216       px = alpaka::math::abs(acc, px);
0217     if (y3 < y2 && y2 < y1)
0218       py = -alpaka::math::abs(acc, py);
0219     else if (y3 > y2 && y2 > y1)
0220       py = alpaka::math::abs(acc, py);
0221 
0222     float AO = alpaka::math::sqrt(
0223         acc, (x_other - x_center) * (x_other - x_center) + (y_other - y_center) * (y_other - y_center));
0224     float BO =
0225         alpaka::math::sqrt(acc, (x_init - x_center) * (x_init - x_center) + (y_init - y_center) * (y_init - y_center));
0226     float AB2 = (x_other - x_init) * (x_other - x_init) + (y_other - y_init) * (y_other - y_init);
0227     float dPhi = alpaka::math::acos(acc, (AO * AO + BO * BO - AB2) / (2 * AO * BO));  //Law of Cosines
0228     float ds = circleRadius / 100 * dPhi;
0229     float pz = dz / ds * pt;
0230 
0231     float p = alpaka::math::sqrt(acc, px * px + py * py + pz * pz);
0232     float a = -2.f * k2Rinv1GeVf * 100 * charge;
0233     float rou = a / p;
0234 
0235     float rzChiSquared = 0;
0236     float error = 0;
0237 
0238     //check the tilted module, side: PosZ, NegZ, Center(for not tilted)
0239     float drdz = alpaka::math::abs(acc, modules.drdzs()[outerOuterLowerModuleIndex]);
0240     short side = modules.sides()[outerOuterLowerModuleIndex];
0241     short subdets = modules.subdets()[outerOuterLowerModuleIndex];
0242 
0243     //calculate residual
0244     if (layer3 <= 6 && ((side == lst::Center) or (drdz < 1))) {  // for barrel
0245       float paraA = r_init * r_init + 2 * (px * px + py * py) / (a * a) + 2 * (y_init * px - x_init * py) / a -
0246                     r_target * r_target;
0247       float paraB = 2 * (x_init * px + y_init * py) / a;
0248       float paraC = 2 * (y_init * px - x_init * py) / a + 2 * (px * px + py * py) / (a * a);
0249       float A = paraB * paraB + paraC * paraC;
0250       float B = 2 * paraA * paraB;
0251       float C = paraA * paraA - paraC * paraC;
0252       float sol1 = (-B + alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
0253       float sol2 = (-B - alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
0254       float solz1 = alpaka::math::asin(acc, sol1) / rou * pz / p + z_init;
0255       float solz2 = alpaka::math::asin(acc, sol2) / rou * pz / p + z_init;
0256       float diffz1 = (solz1 - z_target) * 100;
0257       float diffz2 = (solz2 - z_target) * 100;
0258       if (edm::isNotFinite(diffz1))
0259         residual = diffz2;
0260       else if (edm::isNotFinite(diffz2))
0261         residual = diffz1;
0262       else {
0263         residual = (alpaka::math::abs(acc, diffz1) < alpaka::math::abs(acc, diffz2)) ? diffz1 : diffz2;
0264       }
0265     } else {  // for endcap
0266       float s = (z_target - z_init) * p / pz;
0267       float x = x_init + px / a * alpaka::math::sin(acc, rou * s) - py / a * (1 - alpaka::math::cos(acc, rou * s));
0268       float y = y_init + py / a * alpaka::math::sin(acc, rou * s) + px / a * (1 - alpaka::math::cos(acc, rou * s));
0269       residual = (r_target - alpaka::math::sqrt(acc, x * x + y * y)) * 100;
0270     }
0271 
0272     // error
0273     if (moduleType3 == 0) {
0274       error = 0.15f;  //PS
0275     } else {
0276       error = 5.0f;  //2S
0277     }
0278 
0279     float projection_missing2 = 1;
0280     if (drdz < 1)
0281       projection_missing2 = ((subdets == lst::Endcap) or (side == lst::Center))
0282                                 ? 1.f
0283                                 : 1 / (1 + drdz * drdz);  // cos(atan(drdz)), if dr/dz<1
0284     if (drdz > 1)
0285       projection_missing2 = ((subdets == lst::Endcap) or (side == lst::Center))
0286                                 ? 1.f
0287                                 : drdz * drdz / (1 + drdz * drdz);  //sin(atan(drdz)), if dr/dz>1
0288 
0289     rzChiSquared = 12 * (residual * residual) / (error * error * projection_missing2);
0290 
0291     //helix calculation returns NaN, use linear approximation
0292     if (edm::isNotFinite(rzChiSquared) || circleRadius < 0) {
0293       float slope = (z3 - z1) / (r3 - r1);
0294 
0295       residual = (layer3 <= 6) ? ((z3 - z1) - slope * (r3 - r1)) : ((r3 - r1) - (z3 - z1) / slope);
0296       residual = (moduleType3 == 0) ? residual / 0.15f : residual / 5.0f;
0297 
0298       rzChiSquared = 12 * residual * residual;
0299       return rzChiSquared < 2.8e-4;
0300     }
0301 
0302     //cuts for different regions
0303     //region definitions: https://github.com/user-attachments/assets/2b3c1425-66eb-4524-83de-deb6f3b31f71
0304     //for the logic behind the cuts, see https://github.com/SegmentLinking/cmssw/pull/92
0305     if (layer1 == 7) {
0306       if (layer2 == 8) {
0307         if (layer3 == 9) {
0308           return rzChiSquared < 65.47191f;  // Region 0
0309         } else if (layer3 == 14) {
0310           return rzChiSquared < 3.3200853f;  // Region 1
0311         }
0312       } else if (layer2 == 13) {
0313         return rzChiSquared < 17.194584f;  // Region 2
0314       }
0315     } else if (layer1 == 8) {
0316       if (layer2 == 9) {
0317         if (layer3 == 10) {
0318           return rzChiSquared < 114.91959f;  // Region 3
0319         } else if (layer3 == 15) {
0320           return rzChiSquared < 3.4359624f;  // Region 4
0321         }
0322       } else if (layer2 == 14) {
0323         return rzChiSquared < 4.6487956f;  // Region 5
0324       }
0325     } else if (layer1 == 9) {
0326       if (layer2 == 10) {
0327         if (layer3 == 11) {
0328           return rzChiSquared < 97.34339f;  // Region 6
0329         } else if (layer3 == 16) {
0330           return rzChiSquared < 3.095819f;  // Region 7
0331         }
0332       } else if (layer2 == 15) {
0333         return rzChiSquared < 11.477617f;  // Region 8
0334       }
0335     } else if (layer1 == 1) {
0336       if (layer3 == 7) {
0337         return rzChiSquared < 96.949936f;  // Region 10
0338       } else if (layer3 == 3) {
0339         return rzChiSquared < 458.43982f;  // Region 11
0340       }
0341     } else if (layer1 == 2) {
0342       if (layer2 == 7) {
0343         if (layer3 == 8) {
0344           return rzChiSquared < 218.82303f;  // Region 12
0345         } else if (layer3 == 13) {
0346           return rzChiSquared < 3.155554f;  // Region 13
0347         }
0348       } else if (layer2 == 3) {
0349         if (layer3 == 7) {
0350           return rzChiSquared < 235.5005f;  // Region 14
0351         } else if (layer3 == 12) {
0352           return rzChiSquared < 3.8522234f;  // Region 15
0353         } else if (layer3 == 4) {
0354           return rzChiSquared < 3.5852437f;  // Region 16
0355         }
0356       }
0357     } else if (layer1 == 3) {
0358       if (layer2 == 7) {
0359         if (layer3 == 8) {
0360           return rzChiSquared < 42.68f;  // Region 17
0361         } else if (layer3 == 13) {
0362           return rzChiSquared < 3.853796f;  // Region 18
0363         }
0364       } else if (layer2 == 12) {
0365         return rzChiSquared < 6.2774787f;  // Region 19
0366       }
0367     }
0368     return false;
0369   }
0370 
0371   template <typename TAcc>
0372   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passPointingConstraintBBB(TAcc const& acc,
0373                                                                 ModulesConst modules,
0374                                                                 MiniDoubletsConst mds,
0375                                                                 SegmentsConst segments,
0376                                                                 uint16_t innerInnerLowerModuleIndex,
0377                                                                 uint16_t middleLowerModuleIndex,
0378                                                                 uint16_t outerOuterLowerModuleIndex,
0379                                                                 unsigned int firstMDIndex,
0380                                                                 unsigned int secondMDIndex,
0381                                                                 unsigned int thirdMDIndex,
0382                                                                 float& zOut,
0383                                                                 float& rtOut,
0384                                                                 unsigned int innerSegmentIndex,
0385                                                                 float& betaIn,
0386                                                                 float& betaInCut,
0387                                                                 const float ptCut) {
0388     float rtIn = mds.anchorRt()[firstMDIndex];
0389     float rtMid = mds.anchorRt()[secondMDIndex];
0390     float drt_InSeg = rtMid - rtIn;
0391 
0392     // raw betaIn value without any correction, based on the mini-doublet hit positions
0393     float alpha_InLo = __H2F(segments.dPhiChanges()[innerSegmentIndex]);
0394     float tl_axis_x = mds.anchorX()[thirdMDIndex] - mds.anchorX()[firstMDIndex];
0395     float tl_axis_y = mds.anchorY()[thirdMDIndex] - mds.anchorY()[firstMDIndex];
0396     betaIn = alpha_InLo - phi_mpi_pi(acc, phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[firstMDIndex]);
0397 
0398     //beta computation
0399     float drt_tl_axis = alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y);
0400 
0401     //innerOuterAnchor - innerInnerAnchor
0402     const float rt_InSeg = alpaka::math::sqrt(acc,
0403                                               (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) *
0404                                                       (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) +
0405                                                   (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]) *
0406                                                       (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]));
0407     betaInCut =
0408         alpaka::math::asin(acc, alpaka::math::min(acc, (-rt_InSeg + drt_tl_axis) * k2Rinv1GeVf / ptCut, kSinAlphaMax)) +
0409         (0.02f / drt_InSeg);
0410 
0411     //Beta cut
0412     return alpaka::math::abs(acc, betaIn) < betaInCut;
0413   }
0414 
0415   template <typename TAcc>
0416   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passPointingConstraintBBE(TAcc const& acc,
0417                                                                 ModulesConst modules,
0418                                                                 MiniDoubletsConst mds,
0419                                                                 SegmentsConst segments,
0420                                                                 uint16_t innerInnerLowerModuleIndex,
0421                                                                 uint16_t middleLowerModuleIndex,
0422                                                                 uint16_t outerOuterLowerModuleIndex,
0423                                                                 unsigned int firstMDIndex,
0424                                                                 unsigned int secondMDIndex,
0425                                                                 unsigned int thirdMDIndex,
0426                                                                 float& zOut,
0427                                                                 float& rtOut,
0428                                                                 uint16_t innerOuterLowerModuleIndex,
0429                                                                 unsigned int innerSegmentIndex,
0430                                                                 unsigned int outerSegmentIndex,
0431                                                                 float& betaIn,
0432                                                                 float& betaInCut,
0433                                                                 const float ptCut) {
0434     float rt_InLo = mds.anchorRt()[firstMDIndex];
0435     float rt_InOut = mds.anchorRt()[secondMDIndex];
0436 
0437     float sdIn_alpha = __H2F(segments.dPhiChanges()[innerSegmentIndex]);
0438 
0439     float tl_axis_x = mds.anchorX()[thirdMDIndex] - mds.anchorX()[firstMDIndex];
0440     float tl_axis_y = mds.anchorY()[thirdMDIndex] - mds.anchorY()[firstMDIndex];
0441 
0442     betaIn = sdIn_alpha - phi_mpi_pi(acc, phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[firstMDIndex]);
0443 
0444     float betaInRHmin = betaIn;
0445     float betaInRHmax = betaIn;
0446 
0447     float swapTemp;
0448 
0449     if (alpaka::math::abs(acc, betaInRHmin) > alpaka::math::abs(acc, betaInRHmax)) {
0450       swapTemp = betaInRHmin;
0451       betaInRHmin = betaInRHmax;
0452       betaInRHmax = swapTemp;
0453     }
0454 
0455     float sdIn_dr = alpaka::math::sqrt(acc,
0456                                        (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) *
0457                                                (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) +
0458                                            (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]) *
0459                                                (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]));
0460     float sdIn_d = rt_InOut - rt_InLo;
0461 
0462     float dr = alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y);
0463     betaInCut = alpaka::math::asin(acc, alpaka::math::min(acc, (-sdIn_dr + dr) * k2Rinv1GeVf / ptCut, kSinAlphaMax)) +
0464                 (0.02f / sdIn_d);
0465 
0466     //Beta cut
0467     return alpaka::math::abs(acc, betaInRHmin) < betaInCut;
0468   }
0469 
0470   template <typename TAcc>
0471   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passPointingConstraintEEE(TAcc const& acc,
0472                                                                 ModulesConst modules,
0473                                                                 MiniDoubletsConst mds,
0474                                                                 SegmentsConst segments,
0475                                                                 uint16_t innerInnerLowerModuleIndex,
0476                                                                 uint16_t middleLowerModuleIndex,
0477                                                                 uint16_t outerOuterLowerModuleIndex,
0478                                                                 unsigned int firstMDIndex,
0479                                                                 unsigned int secondMDIndex,
0480                                                                 unsigned int thirdMDIndex,
0481                                                                 float& zOut,
0482                                                                 float& rtOut,
0483                                                                 unsigned int innerSegmentIndex,
0484                                                                 unsigned int outerSegmentIndex,
0485                                                                 float& betaIn,
0486                                                                 float& betaInCut,
0487                                                                 const float ptCut) {
0488     float rt_InLo = mds.anchorRt()[firstMDIndex];
0489     float rt_InOut = mds.anchorRt()[secondMDIndex];
0490     float sdIn_alpha = __H2F(segments.dPhiChanges()[innerSegmentIndex]);
0491 
0492     float tl_axis_x = mds.anchorX()[thirdMDIndex] - mds.anchorX()[firstMDIndex];
0493     float tl_axis_y = mds.anchorY()[thirdMDIndex] - mds.anchorY()[firstMDIndex];
0494 
0495     betaIn = sdIn_alpha - phi_mpi_pi(acc, phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[firstMDIndex]);
0496 
0497     float sdIn_alphaRHmin = __H2F(segments.dPhiChangeMins()[innerSegmentIndex]);
0498     float sdIn_alphaRHmax = __H2F(segments.dPhiChangeMaxs()[innerSegmentIndex]);
0499     float betaInRHmin = betaIn + sdIn_alphaRHmin - sdIn_alpha;
0500     float betaInRHmax = betaIn + sdIn_alphaRHmax - sdIn_alpha;
0501 
0502     float swapTemp;
0503 
0504     if (alpaka::math::abs(acc, betaInRHmin) > alpaka::math::abs(acc, betaInRHmax)) {
0505       swapTemp = betaInRHmin;
0506       betaInRHmin = betaInRHmax;
0507       betaInRHmax = swapTemp;
0508     }
0509     float sdIn_dr = alpaka::math::sqrt(acc,
0510                                        (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) *
0511                                                (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) +
0512                                            (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]) *
0513                                                (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]));
0514     float sdIn_d = rt_InOut - rt_InLo;
0515 
0516     float dr = alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y);
0517     betaInCut = alpaka::math::asin(acc, alpaka::math::min(acc, (-sdIn_dr + dr) * k2Rinv1GeVf / ptCut, kSinAlphaMax)) +
0518                 (0.02f / sdIn_d);
0519 
0520     //Beta cut
0521     return alpaka::math::abs(acc, betaInRHmin) < betaInCut;
0522   }
0523 
0524   template <typename TAcc>
0525   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passPointingConstraint(TAcc const& acc,
0526                                                              ModulesConst modules,
0527                                                              MiniDoubletsConst mds,
0528                                                              SegmentsConst segments,
0529                                                              uint16_t innerInnerLowerModuleIndex,
0530                                                              uint16_t middleLowerModuleIndex,
0531                                                              uint16_t outerOuterLowerModuleIndex,
0532                                                              unsigned int firstMDIndex,
0533                                                              unsigned int secondMDIndex,
0534                                                              unsigned int thirdMDIndex,
0535                                                              float& zOut,
0536                                                              float& rtOut,
0537                                                              uint16_t innerOuterLowerModuleIndex,
0538                                                              unsigned int innerSegmentIndex,
0539                                                              unsigned int outerSegmentIndex,
0540                                                              float& betaIn,
0541                                                              float& betaInCut,
0542                                                              const float ptCut) {
0543     short innerInnerLowerModuleSubdet = modules.subdets()[innerInnerLowerModuleIndex];
0544     short middleLowerModuleSubdet = modules.subdets()[middleLowerModuleIndex];
0545     short outerOuterLowerModuleSubdet = modules.subdets()[outerOuterLowerModuleIndex];
0546 
0547     if (innerInnerLowerModuleSubdet == Barrel and middleLowerModuleSubdet == Barrel and
0548         outerOuterLowerModuleSubdet == Barrel) {
0549       return passPointingConstraintBBB(acc,
0550                                        modules,
0551                                        mds,
0552                                        segments,
0553                                        innerInnerLowerModuleIndex,
0554                                        middleLowerModuleIndex,
0555                                        outerOuterLowerModuleIndex,
0556                                        firstMDIndex,
0557                                        secondMDIndex,
0558                                        thirdMDIndex,
0559                                        zOut,
0560                                        rtOut,
0561                                        innerSegmentIndex,
0562                                        betaIn,
0563                                        betaInCut,
0564                                        ptCut);
0565     } else if (innerInnerLowerModuleSubdet == Barrel and middleLowerModuleSubdet == Barrel and
0566                outerOuterLowerModuleSubdet == Endcap) {
0567       return passPointingConstraintBBE(acc,
0568                                        modules,
0569                                        mds,
0570                                        segments,
0571                                        innerInnerLowerModuleIndex,
0572                                        middleLowerModuleIndex,
0573                                        outerOuterLowerModuleIndex,
0574                                        firstMDIndex,
0575                                        secondMDIndex,
0576                                        thirdMDIndex,
0577                                        zOut,
0578                                        rtOut,
0579                                        innerOuterLowerModuleIndex,
0580                                        innerSegmentIndex,
0581                                        outerSegmentIndex,
0582                                        betaIn,
0583                                        betaInCut,
0584                                        ptCut);
0585     } else if (innerInnerLowerModuleSubdet == Barrel and middleLowerModuleSubdet == Endcap and
0586                outerOuterLowerModuleSubdet == Endcap) {
0587       return passPointingConstraintBBE(acc,
0588                                        modules,
0589                                        mds,
0590                                        segments,
0591                                        innerInnerLowerModuleIndex,
0592                                        middleLowerModuleIndex,
0593                                        outerOuterLowerModuleIndex,
0594                                        firstMDIndex,
0595                                        secondMDIndex,
0596                                        thirdMDIndex,
0597                                        zOut,
0598                                        rtOut,
0599                                        innerOuterLowerModuleIndex,
0600                                        innerSegmentIndex,
0601                                        outerSegmentIndex,
0602                                        betaIn,
0603                                        betaInCut,
0604                                        ptCut);
0605 
0606     }
0607 
0608     else if (innerInnerLowerModuleSubdet == Endcap and middleLowerModuleSubdet == Endcap and
0609              outerOuterLowerModuleSubdet == Endcap) {
0610       return passPointingConstraintEEE(acc,
0611                                        modules,
0612                                        mds,
0613                                        segments,
0614                                        innerInnerLowerModuleIndex,
0615                                        middleLowerModuleIndex,
0616                                        outerOuterLowerModuleIndex,
0617                                        firstMDIndex,
0618                                        secondMDIndex,
0619                                        thirdMDIndex,
0620                                        zOut,
0621                                        rtOut,
0622                                        innerSegmentIndex,
0623                                        outerSegmentIndex,
0624                                        betaIn,
0625                                        betaInCut,
0626                                        ptCut);
0627     }
0628     return false;  // failsafe
0629   }
0630 
0631   template <typename TAcc>
0632   ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeRadiusFromThreeAnchorHits(
0633       TAcc const& acc, float x1, float y1, float x2, float y2, float x3, float y3, float& g, float& f) {
0634     float radius = 0.f;
0635 
0636     //(g,f) -> center
0637     //first anchor hit - (x1,y1), second anchor hit - (x2,y2), third anchor hit - (x3, y3)
0638 
0639     float denomInv = 1.0f / ((y1 - y3) * (x2 - x3) - (x1 - x3) * (y2 - y3));
0640 
0641     float xy1sqr = x1 * x1 + y1 * y1;
0642 
0643     float xy2sqr = x2 * x2 + y2 * y2;
0644 
0645     float xy3sqr = x3 * x3 + y3 * y3;
0646 
0647     g = 0.5f * ((y3 - y2) * xy1sqr + (y1 - y3) * xy2sqr + (y2 - y1) * xy3sqr) * denomInv;
0648 
0649     f = 0.5f * ((x2 - x3) * xy1sqr + (x3 - x1) * xy2sqr + (x1 - x2) * xy3sqr) * denomInv;
0650 
0651     float c = ((x2 * y3 - x3 * y2) * xy1sqr + (x3 * y1 - x1 * y3) * xy2sqr + (x1 * y2 - x2 * y1) * xy3sqr) * denomInv;
0652 
0653     if (((y1 - y3) * (x2 - x3) - (x1 - x3) * (y2 - y3) == 0) || (g * g + f * f - c < 0)) {
0654 #ifdef WARNINGS
0655       printf("three collinear points or FATAL! r^2 < 0!\n");
0656 #endif
0657       radius = -1.f;
0658     } else
0659       radius = alpaka::math::sqrt(acc, g * g + f * f - c);
0660 
0661     return radius;
0662   }
0663 
0664   template <typename TAcc>
0665   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runTripletConstraintsAndAlgo(TAcc const& acc,
0666                                                                    ModulesConst modules,
0667                                                                    MiniDoubletsConst mds,
0668                                                                    SegmentsConst segments,
0669                                                                    uint16_t innerInnerLowerModuleIndex,
0670                                                                    uint16_t middleLowerModuleIndex,
0671                                                                    uint16_t outerOuterLowerModuleIndex,
0672                                                                    unsigned int innerSegmentIndex,
0673                                                                    unsigned int outerSegmentIndex,
0674                                                                    float& zOut,
0675                                                                    float& rtOut,
0676                                                                    float& betaIn,
0677                                                                    float& betaInCut,
0678                                                                    float& circleRadius,
0679                                                                    float& circleCenterX,
0680                                                                    float& circleCenterY,
0681                                                                    const float ptCut) {
0682     //this cut reduces the number of candidates by a factor of 4, i.e., 3 out of 4 warps can end right here!
0683     if (segments.mdIndices()[innerSegmentIndex][1] != segments.mdIndices()[outerSegmentIndex][0])
0684       return false;
0685 
0686     unsigned int firstMDIndex = segments.mdIndices()[innerSegmentIndex][0];
0687     unsigned int secondMDIndex = segments.mdIndices()[outerSegmentIndex][0];
0688     unsigned int thirdMDIndex = segments.mdIndices()[outerSegmentIndex][1];
0689 
0690     float x1 = mds.anchorX()[firstMDIndex];
0691     float x2 = mds.anchorX()[secondMDIndex];
0692     float x3 = mds.anchorX()[thirdMDIndex];
0693     float y1 = mds.anchorY()[firstMDIndex];
0694     float y2 = mds.anchorY()[secondMDIndex];
0695     float y3 = mds.anchorY()[thirdMDIndex];
0696 
0697     circleRadius = computeRadiusFromThreeAnchorHits(acc, x1, y1, x2, y2, x3, y3, circleCenterX, circleCenterY);
0698 
0699     if (not passRZConstraint(acc,
0700                              modules,
0701                              mds,
0702                              innerInnerLowerModuleIndex,
0703                              middleLowerModuleIndex,
0704                              outerOuterLowerModuleIndex,
0705                              firstMDIndex,
0706                              secondMDIndex,
0707                              thirdMDIndex,
0708                              circleRadius,
0709                              circleCenterX,
0710                              circleCenterY))
0711       return false;
0712 
0713     if (not passPointingConstraint(acc,
0714                                    modules,
0715                                    mds,
0716                                    segments,
0717                                    innerInnerLowerModuleIndex,
0718                                    middleLowerModuleIndex,
0719                                    outerOuterLowerModuleIndex,
0720                                    firstMDIndex,
0721                                    secondMDIndex,
0722                                    thirdMDIndex,
0723                                    zOut,
0724                                    rtOut,
0725                                    middleLowerModuleIndex,
0726                                    innerSegmentIndex,
0727                                    outerSegmentIndex,
0728                                    betaIn,
0729                                    betaInCut,
0730                                    ptCut))
0731       return false;
0732 
0733     return true;
0734   }
0735 
0736   struct CreateTriplets {
0737     template <typename TAcc>
0738     ALPAKA_FN_ACC void operator()(TAcc const& acc,
0739                                   ModulesConst modules,
0740                                   MiniDoubletsConst mds,
0741                                   SegmentsConst segments,
0742                                   SegmentsOccupancyConst segmentsOccupancy,
0743                                   Triplets triplets,
0744                                   TripletsOccupancy tripletsOccupancy,
0745                                   ObjectRangesConst ranges,
0746                                   uint16_t* index_gpu,
0747                                   uint16_t nonZeroModules,
0748                                   const float ptCut) const {
0749       auto const globalThreadIdx = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc);
0750       auto const gridThreadExtent = alpaka::getWorkDiv<alpaka::Grid, alpaka::Threads>(acc);
0751 
0752       for (uint16_t innerLowerModuleArrayIdx = globalThreadIdx[0]; innerLowerModuleArrayIdx < nonZeroModules;
0753            innerLowerModuleArrayIdx += gridThreadExtent[0]) {
0754         uint16_t innerInnerLowerModuleIndex = index_gpu[innerLowerModuleArrayIdx];
0755         if (innerInnerLowerModuleIndex >= modules.nLowerModules())
0756           continue;
0757 
0758         uint16_t nConnectedModules = modules.nConnectedModules()[innerInnerLowerModuleIndex];
0759         if (nConnectedModules == 0)
0760           continue;
0761 
0762         unsigned int nInnerSegments = segmentsOccupancy.nSegments()[innerInnerLowerModuleIndex];
0763         for (unsigned int innerSegmentArrayIndex = globalThreadIdx[1]; innerSegmentArrayIndex < nInnerSegments;
0764              innerSegmentArrayIndex += gridThreadExtent[1]) {
0765           unsigned int innerSegmentIndex =
0766               ranges.segmentRanges()[innerInnerLowerModuleIndex][0] + innerSegmentArrayIndex;
0767 
0768           // middle lower module - outer lower module of inner segment
0769           uint16_t middleLowerModuleIndex = segments.outerLowerModuleIndices()[innerSegmentIndex];
0770 
0771           unsigned int nOuterSegments = segmentsOccupancy.nSegments()[middleLowerModuleIndex];
0772           for (unsigned int outerSegmentArrayIndex = globalThreadIdx[2]; outerSegmentArrayIndex < nOuterSegments;
0773                outerSegmentArrayIndex += gridThreadExtent[2]) {
0774             unsigned int outerSegmentIndex = ranges.segmentRanges()[middleLowerModuleIndex][0] + outerSegmentArrayIndex;
0775 
0776             uint16_t outerOuterLowerModuleIndex = segments.outerLowerModuleIndices()[outerSegmentIndex];
0777 
0778             float zOut, rtOut, betaIn, betaInCut, circleRadius, circleCenterX, circleCenterY;
0779 
0780             bool success = runTripletConstraintsAndAlgo(acc,
0781                                                         modules,
0782                                                         mds,
0783                                                         segments,
0784                                                         innerInnerLowerModuleIndex,
0785                                                         middleLowerModuleIndex,
0786                                                         outerOuterLowerModuleIndex,
0787                                                         innerSegmentIndex,
0788                                                         outerSegmentIndex,
0789                                                         zOut,
0790                                                         rtOut,
0791                                                         betaIn,
0792                                                         betaInCut,
0793                                                         circleRadius,
0794                                                         circleCenterX,
0795                                                         circleCenterY,
0796                                                         ptCut);
0797 
0798             if (success) {
0799               unsigned int totOccupancyTriplets =
0800                   alpaka::atomicAdd(acc,
0801                                     &tripletsOccupancy.totOccupancyTriplets()[innerInnerLowerModuleIndex],
0802                                     1u,
0803                                     alpaka::hierarchy::Threads{});
0804               if (static_cast<int>(totOccupancyTriplets) >=
0805                   ranges.tripletModuleOccupancy()[innerInnerLowerModuleIndex]) {
0806 #ifdef WARNINGS
0807                 printf("Triplet excess alert! Module index = %d, Occupancy = %d\n",
0808                        innerInnerLowerModuleIndex,
0809                        totOccupancyTriplets);
0810 #endif
0811               } else {
0812                 unsigned int tripletModuleIndex = alpaka::atomicAdd(
0813                     acc, &tripletsOccupancy.nTriplets()[innerInnerLowerModuleIndex], 1u, alpaka::hierarchy::Threads{});
0814                 unsigned int tripletIndex =
0815                     ranges.tripletModuleIndices()[innerInnerLowerModuleIndex] + tripletModuleIndex;
0816                 addTripletToMemory(modules,
0817                                    mds,
0818                                    segments,
0819                                    triplets,
0820                                    innerSegmentIndex,
0821                                    outerSegmentIndex,
0822                                    innerInnerLowerModuleIndex,
0823                                    middleLowerModuleIndex,
0824                                    outerOuterLowerModuleIndex,
0825 #ifdef CUT_VALUE_DEBUG
0826                                    zOut,
0827                                    rtOut,
0828 #endif
0829                                    betaIn,
0830                                    betaInCut,
0831                                    circleRadius,
0832                                    circleCenterX,
0833                                    circleCenterY,
0834                                    tripletIndex);
0835               }
0836             }
0837           }
0838         }
0839       }
0840     }
0841   };
0842 
0843   struct CreateTripletArrayRanges {
0844     template <typename TAcc>
0845     ALPAKA_FN_ACC void operator()(TAcc const& acc,
0846                                   ModulesConst modules,
0847                                   ObjectRanges ranges,
0848                                   SegmentsOccupancyConst segmentsOccupancy,
0849                                   const float ptCut) const {
0850       // implementation is 1D with a single block
0851       static_assert(std::is_same_v<TAcc, ALPAKA_ACCELERATOR_NAMESPACE::Acc1D>, "Should be Acc1D");
0852       ALPAKA_ASSERT_ACC((alpaka::getWorkDiv<alpaka::Grid, alpaka::Blocks>(acc)[0] == 1));
0853 
0854       auto const globalThreadIdx = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc);
0855       auto const gridThreadExtent = alpaka::getWorkDiv<alpaka::Grid, alpaka::Threads>(acc);
0856 
0857       // Initialize variables in shared memory and set to 0
0858       int& nTotalTriplets = alpaka::declareSharedVar<int, __COUNTER__>(acc);
0859       if (cms::alpakatools::once_per_block(acc)) {
0860         nTotalTriplets = 0;
0861       }
0862       alpaka::syncBlockThreads(acc);
0863 
0864       // Occupancy matrix for 0.8 GeV pT Cut
0865       constexpr int p08_occupancy_matrix[4][4] = {
0866           {543, 235, 88, 46},  // category 0
0867           {755, 347, 0, 0},    // category 1
0868           {0, 0, 0, 0},        // category 2
0869           {0, 38, 46, 39}      // category 3
0870       };
0871 
0872       // Occupancy matrix for 0.6 GeV pT Cut, 99.9%
0873       constexpr int p06_occupancy_matrix[4][4] = {
0874           {1146, 544, 216, 83},  // category 0
0875           {1032, 275, 0, 0},     // category 1
0876           {0, 0, 0, 0},          // category 2
0877           {0, 115, 110, 76}      // category 3
0878       };
0879 
0880       // Select the appropriate occupancy matrix based on ptCut
0881       const auto& occupancy_matrix = (ptCut < 0.8f) ? p06_occupancy_matrix : p08_occupancy_matrix;
0882 
0883       for (uint16_t i = globalThreadIdx[0]; i < modules.nLowerModules(); i += gridThreadExtent[0]) {
0884         if (segmentsOccupancy.nSegments()[i] == 0) {
0885           ranges.tripletModuleIndices()[i] = nTotalTriplets;
0886           ranges.tripletModuleOccupancy()[i] = 0;
0887           continue;
0888         }
0889 
0890         short module_rings = modules.rings()[i];
0891         short module_layers = modules.layers()[i];
0892         short module_subdets = modules.subdets()[i];
0893         float module_eta = alpaka::math::abs(acc, modules.eta()[i]);
0894 
0895         int category_number = getCategoryNumber(module_layers, module_subdets, module_rings);
0896         int eta_number = getEtaBin(module_eta);
0897 
0898         int occupancy = 0;
0899         if (category_number != -1 && eta_number != -1) {
0900           occupancy = occupancy_matrix[category_number][eta_number];
0901         }
0902 #ifdef WARNINGS
0903         else {
0904           printf("Unhandled case in createTripletArrayRanges! Module index = %i\n", i);
0905         }
0906 #endif
0907 
0908         ranges.tripletModuleOccupancy()[i] = occupancy;
0909         unsigned int nTotT = alpaka::atomicAdd(acc, &nTotalTriplets, occupancy, alpaka::hierarchy::Threads{});
0910         ranges.tripletModuleIndices()[i] = nTotT;
0911       }
0912 
0913       // Wait for all threads to finish before reporting final values
0914       alpaka::syncBlockThreads(acc);
0915       if (cms::alpakatools::once_per_block(acc)) {
0916         ranges.nTotalTrips() = nTotalTriplets;
0917       }
0918     }
0919   };
0920 
0921   struct AddTripletRangesToEventExplicit {
0922     template <typename TAcc>
0923     ALPAKA_FN_ACC void operator()(TAcc const& acc,
0924                                   ModulesConst modules,
0925                                   TripletsOccupancyConst tripletsOccupancy,
0926                                   ObjectRanges ranges) const {
0927       // implementation is 1D with a single block
0928       static_assert(std::is_same_v<TAcc, ALPAKA_ACCELERATOR_NAMESPACE::Acc1D>, "Should be Acc1D");
0929       ALPAKA_ASSERT_ACC((alpaka::getWorkDiv<alpaka::Grid, alpaka::Blocks>(acc)[0] == 1));
0930 
0931       auto const globalThreadIdx = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc);
0932       auto const gridThreadExtent = alpaka::getWorkDiv<alpaka::Grid, alpaka::Threads>(acc);
0933 
0934       for (uint16_t i = globalThreadIdx[0]; i < modules.nLowerModules(); i += gridThreadExtent[0]) {
0935         if (tripletsOccupancy.nTriplets()[i] == 0) {
0936           ranges.tripletRanges()[i][0] = -1;
0937           ranges.tripletRanges()[i][1] = -1;
0938         } else {
0939           ranges.tripletRanges()[i][0] = ranges.tripletModuleIndices()[i];
0940           ranges.tripletRanges()[i][1] = ranges.tripletModuleIndices()[i] + tripletsOccupancy.nTriplets()[i] - 1;
0941         }
0942       }
0943     }
0944   };
0945 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE::lst
0946 #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!