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