Phase2L1GT/plugins/L1GTCorrelationalCut.h

0001 #ifndef L1Trigger_Phase2L1GT_L1GTCorrelationalCut_h
0002 #define L1Trigger_Phase2L1GT_L1GTCorrelationalCut_h
0003
0004 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "DataFormats/L1Trigger/interface/P2GTCandidate.h"
0007
0008 #include "L1Trigger/Phase2L1GT/interface/L1GTInvariantMassError.h"
0009
0010 #include "L1Trigger/Phase2L1GT/interface/L1GTScales.h"
0011
0012 #include "L1GTSingleInOutLUT.h"
0013 #include "L1GTOptionalParam.h"
0014
0015 #include <optional>
0016 #include <cinttypes>
0017
0018 namespace l1t {
0019
0020   class L1GTCorrelationalCut {
0021   public:
0022     L1GTCorrelationalCut(const edm::ParameterSet& config,
0023                          const edm::ParameterSet& lutConfig,
0024                          const L1GTScales& scales,
0025                          bool enable_sanity_checks = false,
0026                          bool inv_mass_checks = false)
0027         : scales_(scales),
0028           coshEtaLUT_(lutConfig.getParameterSet("cosh_eta_lut")),
0029           coshEtaLUT2_(lutConfig.getParameterSet("cosh_eta_lut2")),
0030           cosPhiLUT_(lutConfig.getParameterSet("cos_phi_lut")),
0031           minDEta_(getOptionalParam<int, double>(
0032               "minDEta", config, [&scales](double value) { return scales.to_hw_eta_floor(value); })),
0033           maxDEta_(getOptionalParam<int, double>(
0034               "maxDEta", config, [&scales](double value) { return scales.to_hw_eta_ceil(value); })),
0035           minDPhi_(getOptionalParam<int, double>(
0036               "minDPhi", config, [&scales](double value) { return scales.to_hw_phi_floor(value); })),
0037           maxDPhi_(getOptionalParam<int, double>(
0038               "maxDPhi", config, [&scales](double value) { return scales.to_hw_phi_ceil(value); })),
0039           minDz_(getOptionalParam<int, double>(
0040               "minDz", config, [&scales](double value) { return scales.to_hw_z0_floor(value); })),
0041           maxDz_(getOptionalParam<int, double>(
0042               "maxDz", config, [&scales](double value) { return scales.to_hw_z0_ceil(value); })),
0043           minDRSquared_(getOptionalParam<int, double>(
0044               "minDR", config, [&scales](double value) { return scales.to_hw_dRSquared_floor(value); })),
0045           maxDRSquared_(getOptionalParam<int, double>(
0046               "maxDR", config, [&scales](double value) { return scales.to_hw_dRSquared_ceil(value); })),
0047           minInvMassSqrDiv2_scale1_(getOptionalParam<int64_t, double>(
0048               "minInvMass",
0049               config,
0050               [&](double value) {
0051                 return std::floor(scales.to_hw_InvMassSqrDiv2(value) * coshEtaLUT_.output_scale());
0052               })),
0053           maxInvMassSqrDiv2_scale1_(getOptionalParam<int64_t, double>(
0054               "maxInvMass",
0055               config,
0056               [&](double value) {
0057                 return std::ceil(scales.to_hw_InvMassSqrDiv2(value) * coshEtaLUT_.output_scale());
0058               })),
0059           minInvMassSqrDiv2_scale2_(getOptionalParam<int64_t, double>(
0060               "minInvMass",
0061               config,
0062               [&](double value) {
0063                 return std::floor(scales.to_hw_InvMassSqrDiv2(value) * coshEtaLUT2_.output_scale());
0064               })),
0065           maxInvMassSqrDiv2_scale2_(getOptionalParam<int64_t, double>(
0066               "maxInvMass",
0067               config,
0068               [&](double value) {
0069                 return std::ceil(scales.to_hw_InvMassSqrDiv2(value) * coshEtaLUT2_.output_scale());
0070               })),
0071           minTransMassSqrDiv2_(getOptionalParam<int64_t, double>(
0072               "minTransMass",
0073               config,
0074               [&](double value) {
0075                 return std::floor(scales.to_hw_TransMassSqrDiv2(value) * cosPhiLUT_.output_scale());
0076               })),
0077           maxTransMassSqrDiv2_(getOptionalParam<int64_t, double>(
0078               "maxTransMass",
0079               config,
0080               [&](double value) {
0081                 return std::ceil(scales.to_hw_TransMassSqrDiv2(value) * cosPhiLUT_.output_scale());
0082               })),
0083           minPTSquared_(getOptionalParam<int64_t, double>(
0084               "minCombPt",
0085               config,
0086               [&](double value) { return std::floor(scales.to_hw_PtSquared(value) * cosPhiLUT_.output_scale()); })),
0087           maxPTSquared_(getOptionalParam<int64_t, double>(
0088               "maxCombPt",
0089               config,
0090               [&](double value) { return std::ceil(scales.to_hw_PtSquared(value) * cosPhiLUT_.output_scale()); })),
0091           minInvMassSqrOver2DRSqr_scale1_(getOptionalParam<double, double>(
0092               "minInvMassOverDR",
0093               config,
0094               [&](double value) {
0095                 return std::floor(scales.to_hw_InvMassSqrDiv2(value) * coshEtaLUT_.output_scale());
0096               })),
0097           maxInvMassSqrOver2DRSqr_scale1_(getOptionalParam<double, double>(
0098               "maxInvMassOverDR",
0099               config,
0100               [&](double value) {
0101                 return std::ceil(scales.to_hw_InvMassSqrDiv2(value) * coshEtaLUT_.output_scale());
0102               })),
0103           minInvMassSqrOver2DRSqr_scale2_(getOptionalParam<double, double>(
0104               "minInvMassOverDR",
0105               config,
0106               [&](double value) {
0107                 return std::floor(scales.to_hw_InvMassSqrDiv2(value) * coshEtaLUT2_.output_scale());
0108               })),
0109           maxInvMassSqrOver2DRSqr_scale2_(getOptionalParam<double, double>(
0110               "maxInvMassOverDR",
0111               config,
0112               [&](double value) {
0113                 return std::ceil(scales.to_hw_InvMassSqrDiv2(value) * coshEtaLUT2_.output_scale());
0114               })),
0115           os_(config.getParameter<bool>("os")),
0116           ss_(config.getParameter<bool>("ss")),
0117           enable_sanity_checks_(enable_sanity_checks),
0118           inv_mass_checks_(inv_mass_checks) {}
0119
0120     bool checkObjects(const P2GTCandidate& obj1,
0121                       const P2GTCandidate& obj2,
0122                       InvariantMassErrorCollection& massErrors) const {
0123       bool res = true;
0124
0125       std::optional<uint32_t> dEta;
0126
0127       if (minDEta_ || maxDEta_ || minDRSquared_ || maxDRSquared_ || minInvMassSqrDiv2_scale1_ ||
0128           maxInvMassSqrDiv2_scale1_ || minInvMassSqrDiv2_scale2_ || maxInvMassSqrDiv2_scale2_ ||
0129           minInvMassSqrOver2DRSqr_scale1_ || maxInvMassSqrOver2DRSqr_scale1_ || minInvMassSqrOver2DRSqr_scale2_ ||
0130           maxInvMassSqrOver2DRSqr_scale2_) {
0131         dEta = std::abs(obj1.hwEta().to_int() - obj2.hwEta().to_int());
0132         res &= minDEta_ ? dEta > minDEta_ : true;
0133         res &= maxDEta_ ? dEta < maxDEta_ : true;
0134       }
0135
0136       constexpr int HW_PI = 1 << (P2GTCandidate::hwPhi_t::width - 1);  // assumes phi in [-pi, pi)
0137
0138       std::optional<uint32_t> dPhi;
0139
0140       if (minDPhi_ || maxDPhi_ || minDRSquared_ || maxDRSquared_ || minInvMassSqrDiv2_scale1_ ||
0141           maxInvMassSqrDiv2_scale1_ || minInvMassSqrDiv2_scale2_ || maxInvMassSqrDiv2_scale2_ || minTransMassSqrDiv2_ ||
0142           maxTransMassSqrDiv2_ || minPTSquared_ || maxPTSquared_ || minInvMassSqrOver2DRSqr_scale1_ ||
0143           maxInvMassSqrOver2DRSqr_scale1_ || minInvMassSqrOver2DRSqr_scale2_ || maxInvMassSqrOver2DRSqr_scale2_) {
0144         // Ensure dPhi is always the smaller angle, i.e. always between [0, pi]
0145         dPhi = HW_PI - std::abs(std::abs(obj1.hwPhi().to_int() - obj2.hwPhi().to_int()) - HW_PI);
0146       }
0147
0148       res &= minDPhi_ ? dPhi > minDPhi_ : true;
0149       res &= maxDPhi_ ? dPhi < maxDPhi_ : true;
0150
0151       if (minDz_ || maxDz_) {
0152         uint32_t dZ = abs(obj1.hwZ0() - obj2.hwZ0());
0153         res &= minDz_ ? dZ > minDz_ : true;
0154         res &= maxDz_ ? dZ < maxDz_ : true;
0155       }
0156
0157       uint32_t dRSquared = 0;
0158       if (minDRSquared_ || maxDRSquared_ || minInvMassSqrOver2DRSqr_scale1_ || maxInvMassSqrOver2DRSqr_scale1_ ||
0159           minInvMassSqrOver2DRSqr_scale2_ || maxInvMassSqrOver2DRSqr_scale2_) {
0160         dRSquared = dEta.value() * dEta.value() + dPhi.value() * dPhi.value();
0161         res &= minDRSquared_ ? dRSquared > minDRSquared_ : true;
0162         res &= maxDRSquared_ ? dRSquared < maxDRSquared_ : true;
0163       }
0164
0165       res &= os_ ? obj1.hwCharge() != obj2.hwCharge() : true;
0166       res &= ss_ ? obj1.hwCharge() == obj2.hwCharge() : true;
0167
0168       int32_t lutCoshDEta = 0;
0169       if (dEta) {
0170         lutCoshDEta = dEta < L1GTSingleInOutLUT::DETA_LUT_SPLIT
0171                           ? coshEtaLUT_[dEta.value()]
0172                           : coshEtaLUT2_[dEta.value() - L1GTSingleInOutLUT::DETA_LUT_SPLIT];
0173       }
0174
0175       // Optimization: (cos(x + pi) = -cos(x), x in [0, pi))
0176       int32_t lutCosDPhi = 0;
0177       if (dPhi) {
0178         lutCosDPhi = dPhi >= HW_PI ? -cosPhiLUT_[dPhi.value()] : cosPhiLUT_[dPhi.value()];
0179       }
0180
0181       if (enable_sanity_checks_ && dEta && dPhi) {
0182         // Check whether the LUT error is smaller or equal than the expected maximum LUT error
0183         double coshEtaLUTMax =
0184             dEta < L1GTSingleInOutLUT::DETA_LUT_SPLIT ? coshEtaLUT_.hwMax_error() : coshEtaLUT2_.hwMax_error();
0185         double etaLUTScale =
0186             dEta < L1GTSingleInOutLUT::DETA_LUT_SPLIT ? coshEtaLUT_.output_scale() : coshEtaLUT2_.output_scale();
0187
0188         if (std::abs(lutCoshDEta - etaLUTScale * std::cosh(dEta.value() * scales_.eta_lsb())) > coshEtaLUTMax) {
0189           edm::LogError("COSH LUT") << "Difference larger than max LUT error: " << coshEtaLUTMax
0190                                     << ", lut: " << lutCoshDEta
0191                                     << ", calc: " << etaLUTScale * std::cosh(dEta.value() * scales_.eta_lsb())
0192                                     << ", dEta: " << dEta.value() << ", scale: " << etaLUTScale;
0193         }
0194
0195         if (std::abs(lutCosDPhi - cosPhiLUT_.output_scale() * std::cos(dPhi.value() * scales_.phi_lsb())) >
0196             cosPhiLUT_.hwMax_error()) {
0197           edm::LogError("COS LUT") << "Difference larger than max LUT error: " << cosPhiLUT_.hwMax_error()
0198                                    << ", lut: " << lutCosDPhi << ", calc: "
0199                                    << cosPhiLUT_.output_scale() * std::cos(dPhi.value() * scales_.phi_lsb());
0200         }
0201       }
0202
0203       int64_t invMassSqrDiv2 = 0;
0204       if (minInvMassSqrDiv2_scale1_ || maxInvMassSqrDiv2_scale1_ || minInvMassSqrDiv2_scale2_ ||
0205           maxInvMassSqrDiv2_scale2_ || minInvMassSqrOver2DRSqr_scale1_ || maxInvMassSqrOver2DRSqr_scale1_ ||
0206           minInvMassSqrOver2DRSqr_scale2_ || maxInvMassSqrOver2DRSqr_scale2_) {
0207         if (dEta < L1GTSingleInOutLUT::DETA_LUT_SPLIT) {
0208           // dEta [0, 2pi)
0209           invMassSqrDiv2 = obj1.hwPT().to_int64() * obj2.hwPT().to_int64() * (lutCoshDEta - lutCosDPhi);
0210           res &= minInvMassSqrDiv2_scale1_ ? invMassSqrDiv2 > minInvMassSqrDiv2_scale1_ : true;
0211           res &= maxInvMassSqrDiv2_scale1_ ? invMassSqrDiv2 < maxInvMassSqrDiv2_scale1_ : true;
0212         } else {
0213           // dEta [2pi, 4pi), ignore cos
0214           invMassSqrDiv2 = obj1.hwPT().to_int64() * obj2.hwPT().to_int64() * lutCoshDEta;
0215           res &= minInvMassSqrDiv2_scale2_ ? invMassSqrDiv2 > minInvMassSqrDiv2_scale2_ : true;
0216           res &= maxInvMassSqrDiv2_scale2_ ? invMassSqrDiv2 < maxInvMassSqrDiv2_scale2_ : true;
0217         }
0218
0219         if (inv_mass_checks_) {
0220           double precInvMass =
0221               scales_.pT_lsb() *
0222               std::sqrt(2 * obj1.hwPT().to_double() * obj2.hwPT().to_double() *
0223                         (std::cosh(dEta.value() * scales_.eta_lsb()) - std::cos(dPhi.value() * scales_.phi_lsb())));
0224
0225           double lutInvMass =
0226               scales_.pT_lsb() * std::sqrt(2 * static_cast<double>(invMassSqrDiv2) /
0227                                            (dEta < L1GTSingleInOutLUT::DETA_LUT_SPLIT ? coshEtaLUT_.output_scale()
0228                                                                                       : coshEtaLUT2_.output_scale()));
0229
0230           double error = std::abs(precInvMass - lutInvMass);
0231           massErrors.emplace_back(InvariantMassError{error, error / precInvMass, precInvMass});
0232         }
0233       }
0234
0235       if (minPTSquared_ || maxPTSquared_) {
0236         int64_t pTSquared = obj1.hwPT().to_int64() * obj1.hwPT().to_int64() *
0237                                 static_cast<int64_t>(std::round(cosPhiLUT_.output_scale())) +
0238                             obj2.hwPT().to_int64() * obj2.hwPT().to_int64() *
0239                                 static_cast<int64_t>(std::round(cosPhiLUT_.output_scale())) +
0240                             2 * obj1.hwPT().to_int64() * obj2.hwPT().to_int64() * lutCosDPhi;
0241         res &= minPTSquared_ ? pTSquared > minPTSquared_.value() : true;
0242         res &= maxPTSquared_ ? pTSquared < maxPTSquared_.value() : true;
0243       }
0244
0245       if (minTransMassSqrDiv2_ || maxTransMassSqrDiv2_) {
0246         int64_t transMassDiv2 = obj1.hwPT().to_int64() * obj2.hwPT().to_int64() *
0247                                 (static_cast<int64_t>(std::round(cosPhiLUT_.output_scale())) - lutCosDPhi);
0248         res &= minTransMassSqrDiv2_ ? transMassDiv2 > minTransMassSqrDiv2_.value() : true;
0249         res &= maxTransMassSqrDiv2_ ? transMassDiv2 < maxTransMassSqrDiv2_.value() : true;
0250       }
0251
0252       if (minInvMassSqrOver2DRSqr_scale1_ || maxInvMassSqrOver2DRSqr_scale1_ || minInvMassSqrOver2DRSqr_scale2_ ||
0253           maxInvMassSqrOver2DRSqr_scale2_) {
0254         ap_uint<96> invMassSqrDiv2Shift = ap_uint<96>(invMassSqrDiv2)
0255                                           << L1GTScales::INV_MASS_SQR_OVER_2_DR_SQR_RESOLUTION;
0256
0257         if (dEta < L1GTSingleInOutLUT::DETA_LUT_SPLIT) {
0258           res &= minInvMassSqrOver2DRSqr_scale1_
0259                      ? invMassSqrDiv2Shift > minInvMassSqrOver2DRSqr_scale1_.value() * dRSquared
0260                      : true;
0261           res &= maxInvMassSqrOver2DRSqr_scale1_
0262                      ? invMassSqrDiv2Shift < maxInvMassSqrOver2DRSqr_scale1_.value() * dRSquared
0263                      : true;
0264         } else {
0265           res &= minInvMassSqrOver2DRSqr_scale2_
0266                      ? invMassSqrDiv2Shift > minInvMassSqrOver2DRSqr_scale2_.value() * dRSquared
0267                      : true;
0268           res &= maxInvMassSqrOver2DRSqr_scale2_
0269                      ? invMassSqrDiv2Shift < maxInvMassSqrOver2DRSqr_scale2_.value() * dRSquared
0270                      : true;
0271         }
0272       }
0273
0274       return res;
0275     }
0276
0277     static void fillLUTDescriptions(edm::ParameterSetDescription& desc) {
0278       edm::ParameterSetDescription coshLUTDesc;
0279       L1GTSingleInOutLUT::fillLUTDescriptions(coshLUTDesc);
0280       desc.add<edm::ParameterSetDescription>("cosh_eta_lut", coshLUTDesc);
0281
0282       edm::ParameterSetDescription coshLUT2Desc;
0283       L1GTSingleInOutLUT::fillLUTDescriptions(coshLUT2Desc);
0284       desc.add<edm::ParameterSetDescription>("cosh_eta_lut2", coshLUT2Desc);
0285
0286       edm::ParameterSetDescription cosLUTDesc;
0287       L1GTSingleInOutLUT::fillLUTDescriptions(cosLUTDesc);
0288       desc.add<edm::ParameterSetDescription>("cos_phi_lut", cosLUTDesc);
0289     }
0290
0291     static void fillPSetDescription(edm::ParameterSetDescription& desc) {
0292       desc.addOptional<double>("minDEta");
0293       desc.addOptional<double>("maxDEta");
0294       desc.addOptional<double>("minDPhi");
0295       desc.addOptional<double>("maxDPhi");
0296       desc.addOptional<double>("minDR");
0297       desc.addOptional<double>("maxDR");
0298       desc.addOptional<double>("minDz");
0299       desc.addOptional<double>("maxDz");
0300       desc.addOptional<double>("minInvMass");
0301       desc.addOptional<double>("maxInvMass");
0302       desc.addOptional<double>("minTransMass");
0303       desc.addOptional<double>("maxTransMass");
0304       desc.addOptional<double>("minCombPt");
0305       desc.addOptional<double>("maxCombPt");
0306       desc.addOptional<double>("minInvMassOverDR");
0307       desc.addOptional<double>("maxInvMassOverDR");
0308       desc.add<bool>("os", false);
0309       desc.add<bool>("ss", false);
0310     }
0311
0312   private:
0313     const L1GTScales& scales_;
0314
0315     const L1GTSingleInOutLUT coshEtaLUT_;   // [0, 2pi)
0316     const L1GTSingleInOutLUT coshEtaLUT2_;  // [2pi, 4pi)
0317     const L1GTSingleInOutLUT cosPhiLUT_;
0318
0319     const std::optional<int> minDEta_;
0320     const std::optional<int> maxDEta_;
0321     const std::optional<int> minDPhi_;
0322     const std::optional<int> maxDPhi_;
0323     const std::optional<int> minDz_;
0324     const std::optional<int> maxDz_;
0325
0326     const std::optional<int> minDRSquared_;
0327     const std::optional<int> maxDRSquared_;
0328
0329     const std::optional<int64_t> minInvMassSqrDiv2_scale1_;
0330     const std::optional<int64_t> maxInvMassSqrDiv2_scale1_;
0331
0332     const std::optional<int64_t> minInvMassSqrDiv2_scale2_;
0333     const std::optional<int64_t> maxInvMassSqrDiv2_scale2_;
0334
0335     const std::optional<int64_t> minTransMassSqrDiv2_;
0336     const std::optional<int64_t> maxTransMassSqrDiv2_;
0337
0338     const std::optional<int64_t> minPTSquared_;
0339     const std::optional<int64_t> maxPTSquared_;
0340
0341     const std::optional<int64_t> minInvMassSqrOver2DRSqr_scale1_;
0342     const std::optional<int64_t> maxInvMassSqrOver2DRSqr_scale1_;
0343
0344     const std::optional<int64_t> minInvMassSqrOver2DRSqr_scale2_;
0345     const std::optional<int64_t> maxInvMassSqrOver2DRSqr_scale2_;
0346
0347     const bool os_;  // Opposite sign
0348     const bool ss_;  // Same sign
0349
0350     const bool enable_sanity_checks_;
0351     const bool inv_mass_checks_;
0352   };
0353
0354 }  // namespace l1t
0355
0356 #endif  // L1Trigger_Phase2L1GT_L1GTCorrelationalCut_h