Phase2L1GT/plugins/L1GT3BodyCut.h

0001 #ifndef L1Trigger_Phase2L1GT_L1GT3BodyCut_h
0002 #define L1Trigger_Phase2L1GT_L1GT3BodyCut_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 L1GT3BodyCut {
0021   public:
0022     L1GT3BodyCut(const edm::ParameterSet& config,
0023                  const edm::ParameterSet& lutConfig,
0024                  const L1GTScales& scales,
0025                  bool inv_mass_checks = false)
0026         : scales_(scales),
0027           coshEtaLUT_(lutConfig.getParameterSet("cosh_eta_lut")),
0028           coshEtaLUT2_(lutConfig.getParameterSet("cosh_eta_lut2")),
0029           cosPhiLUT_(lutConfig.getParameterSet("cos_phi_lut")),
0030           minInvMassSqrDiv2_(getOptionalParam<int64_t, double>("minInvMass",
0031                                                                config,
0032                                                                [&](double value) {
0033                                                                  return std::round(scales.to_hw_InvMassSqrDiv2(value) *
0034                                                                                    cosPhiLUT_.output_scale());
0035                                                                })),
0036           maxInvMassSqrDiv2_(getOptionalParam<int64_t, double>("maxInvMass",
0037                                                                config,
0038                                                                [&](double value) {
0039                                                                  return std::round(scales.to_hw_InvMassSqrDiv2(value) *
0040                                                                                    cosPhiLUT_.output_scale());
0041                                                                })),
0042           minTransMassSqrDiv2_(getOptionalParam<int64_t, double>(
0043               "minTransMass",
0044               config,
0045               [&](double value) {
0046                 return std::round(scales.to_hw_TransMassSqrDiv2(value) * cosPhiLUT_.output_scale());
0047               })),
0048           maxTransMassSqrDiv2_(getOptionalParam<int64_t, double>(
0049               "maxTransMass",
0050               config,
0051               [&](double value) {
0052                 return std::round(scales.to_hw_TransMassSqrDiv2(value) * cosPhiLUT_.output_scale());
0053               })),
0054           scaleNormalShift_(std::round(std::log2(std::ceil(coshEtaLUT_.output_scale() / coshEtaLUT2_.output_scale())))),
0055           invMassResolutionReduceShift_([&]() {
0056             if (minInvMassSqrDiv2_) {
0057               return std::max<int>(
0058                   std::ceil(std::log2(minInvMassSqrDiv2_.value() * cosPhiLUT_.output_scale() + 1.0)) - 16, 0);
0059             } else if (maxInvMassSqrDiv2_) {
0060               return std::max<int>(std::ceil(std::log2(maxInvMassSqrDiv2_.value() * cosPhiLUT_.output_scale())) - 16,
0061                                    0);
0062             } else {
0063               return 0;
0064             }
0065           }()),
0066           transMassResolutionReduceShift_([&]() {
0067             if (minTransMassSqrDiv2_) {
0068               return std::max<int>(
0069                   std::ceil(std::log2(minTransMassSqrDiv2_.value() * cosPhiLUT_.output_scale() + 1.0)) - 16, 0);
0070             } else if (maxTransMassSqrDiv2_) {
0071               return std::max<int>(std::ceil(std::log2(maxTransMassSqrDiv2_.value() * cosPhiLUT_.output_scale())) - 16,
0072                                    0);
0073             } else {
0074               return 0;
0075             }
0076           }()),
0077           inv_mass_checks_(inv_mass_checks) {}
0078
0079     bool checkObjects(const P2GTCandidate& obj1,
0080                       const P2GTCandidate& obj2,
0081                       const P2GTCandidate& obj3,
0082                       InvariantMassErrorCollection& massErrors) const {
0083       bool res = true;
0084
0085       if (minInvMassSqrDiv2_ || maxInvMassSqrDiv2_) {
0086         int64_t invMassSqrDiv2 = (calc2BodyInvMass(obj1, obj2, massErrors) >> invMassResolutionReduceShift_) +
0087                                  (calc2BodyInvMass(obj1, obj3, massErrors) >> invMassResolutionReduceShift_) +
0088                                  (calc2BodyInvMass(obj2, obj3, massErrors) >> invMassResolutionReduceShift_);
0089
0090         res &= minInvMassSqrDiv2_ ? invMassSqrDiv2 > minInvMassSqrDiv2_.value() >> invMassResolutionReduceShift_ : true;
0091         res &= maxInvMassSqrDiv2_ ? invMassSqrDiv2 < maxInvMassSqrDiv2_.value() >> invMassResolutionReduceShift_ : true;
0092       }
0093
0094       if (minTransMassSqrDiv2_ || maxTransMassSqrDiv2_) {
0095         int64_t transMassDiv2 = (calc2BodyTransMass(obj1, obj2) >> transMassResolutionReduceShift_) +
0096                                 (calc2BodyTransMass(obj1, obj3) >> transMassResolutionReduceShift_) +
0097                                 (calc2BodyTransMass(obj2, obj3) >> transMassResolutionReduceShift_);
0098
0099         res &= minTransMassSqrDiv2_ ? transMassDiv2 > minTransMassSqrDiv2_.value() >> transMassResolutionReduceShift_
0100                                     : true;
0101         res &= maxTransMassSqrDiv2_ ? transMassDiv2 < maxTransMassSqrDiv2_.value() >> transMassResolutionReduceShift_
0102                                     : true;
0103       }
0104
0105       return res;
0106     }
0107
0108     static void fillPSetDescription(edm::ParameterSetDescription& desc) {
0109       desc.addOptional<double>("minInvMass");
0110       desc.addOptional<double>("maxInvMass");
0111       desc.addOptional<double>("minTransMass");
0112       desc.addOptional<double>("maxTransMass");
0113     }
0114
0115   private:
0116     static constexpr int HW_PI = 1 << (P2GTCandidate::hwPhi_t::width - 1);  // assumes phi in [-pi, pi)
0117
0118     int64_t calc2BodyInvMass(const P2GTCandidate& obj1,
0119                              const P2GTCandidate& obj2,
0120                              InvariantMassErrorCollection& massErrors) const {
0121       uint32_t dEta = (obj1.hwEta() > obj2.hwEta()) ? obj1.hwEta().to_int() - obj2.hwEta().to_int()
0122                                                     : obj2.hwEta().to_int() - obj1.hwEta().to_int();
0123       int32_t lutCoshDEta = dEta < L1GTSingleInOutLUT::DETA_LUT_SPLIT
0124                                 ? coshEtaLUT_[dEta]
0125                                 : coshEtaLUT2_[dEta - L1GTSingleInOutLUT::DETA_LUT_SPLIT];
0126
0127       // Ensure dPhi is always the smaller angle, i.e. always between [0, pi]
0128       uint32_t dPhi = HW_PI - abs(abs(obj1.hwPhi().to_int() - obj2.hwPhi().to_int()) - HW_PI);
0129
0130       // Optimization: (cos(x + pi) = -cos(x), x in [0, pi))
0131       int32_t lutCosDPhi = dPhi >= HW_PI ? -cosPhiLUT_[dPhi] : cosPhiLUT_[dPhi];
0132
0133       int64_t invMassSqrDiv2;
0134
0135       if (dEta < L1GTSingleInOutLUT::DETA_LUT_SPLIT) {
0136         // dEta [0, 2pi)
0137         invMassSqrDiv2 = obj1.hwPT().to_int64() * obj2.hwPT().to_int64() * (lutCoshDEta - lutCosDPhi);
0138       } else {
0139         // dEta [2pi, 4pi), ignore cos
0140         invMassSqrDiv2 = obj1.hwPT().to_int64() * obj2.hwPT().to_int64() * lutCoshDEta;
0141       }
0142
0143       if (inv_mass_checks_) {
0144         double precInvMass =
0145             scales_.pT_lsb() * std::sqrt(2 * obj1.hwPT().to_double() * obj2.hwPT().to_double() *
0146                                          (std::cosh(dEta * scales_.eta_lsb()) - std::cos(dPhi * scales_.phi_lsb())));
0147
0148         double lutInvMass =
0149             scales_.pT_lsb() * std::sqrt(2 * static_cast<double>(invMassSqrDiv2) /
0150                                          (dEta < L1GTSingleInOutLUT::DETA_LUT_SPLIT ? coshEtaLUT_.output_scale()
0151                                                                                     : coshEtaLUT2_.output_scale()));
0152
0153         double error = std::abs(precInvMass - lutInvMass);
0154         massErrors.emplace_back(InvariantMassError{error, error / precInvMass, precInvMass});
0155       }
0156
0157       // Normalize scales required due to LUT split in dEta with different scale factors.
0158       return dEta < L1GTSingleInOutLUT::DETA_LUT_SPLIT ? invMassSqrDiv2 : invMassSqrDiv2 << scaleNormalShift_;
0159     }
0160
0161     int64_t calc2BodyTransMass(const P2GTCandidate& obj1, const P2GTCandidate& obj2) const {
0162       // Ensure dPhi is always the smaller angle, i.e. always between [0, pi]
0163       uint32_t dPhi = HW_PI - abs(abs(obj1.hwPhi().to_int() - obj2.hwPhi().to_int()) - HW_PI);
0164
0165       // Optimization: (cos(x + pi) = -cos(x), x in [0, pi))
0166       int32_t lutCosDPhi = dPhi >= HW_PI ? -cosPhiLUT_[dPhi] : cosPhiLUT_[dPhi];
0167
0168       return obj1.hwPT().to_int64() * obj2.hwPT().to_int64() *
0169              (static_cast<int64_t>(std::round(cosPhiLUT_.output_scale())) - lutCosDPhi);
0170     }
0171
0172     const L1GTScales& scales_;
0173
0174     const L1GTSingleInOutLUT coshEtaLUT_;   // [0, 2pi)
0175     const L1GTSingleInOutLUT coshEtaLUT2_;  // [2pi, 4pi)
0176     const L1GTSingleInOutLUT cosPhiLUT_;
0177
0178     const std::optional<int64_t> minInvMassSqrDiv2_;
0179     const std::optional<int64_t> maxInvMassSqrDiv2_;
0180     const std::optional<int64_t> minTransMassSqrDiv2_;
0181     const std::optional<int64_t> maxTransMassSqrDiv2_;
0182
0183     const int scaleNormalShift_;
0184     const int invMassResolutionReduceShift_;
0185     const int transMassResolutionReduceShift_;
0186
0187     const bool inv_mass_checks_;
0188   };
0189
0190 }  // namespace l1t
0191
0192 #endif  // L1Trigger_Phase2L1GT_L1GT3BodyCut_h