Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:21:23

0001 #ifndef L1Trigger_Phase2L1GT_L1GTDeltaCut_h
0002 #define L1Trigger_Phase2L1GT_L1GTDeltaCut_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 "FWCore/Framework/interface/MakerMacros.h"
0011 
0012 #include "L1Trigger/Phase2L1GT/interface/L1GTScales.h"
0013 
0014 #include "L1GTSingleInOutLUT.h"
0015 #include "L1GTOptionalParam.h"
0016 
0017 #include <optional>
0018 
0019 namespace l1t {
0020 
0021   class L1GTDeltaCut {
0022   public:
0023     static constexpr uint32_t DETA_LUT_SPLIT = 1 << 13;  // hw 2pi
0024 
0025     L1GTDeltaCut(const edm::ParameterSet& config,
0026                  const edm::ParameterSet& lutConfig,
0027                  const L1GTScales& scales,
0028                  bool enable_sanity_checks = false,
0029                  bool inv_mass_checks = false)
0030         : scales_(scales),
0031           coshEtaLUT_(lutConfig.getParameterSet("cosh_eta_lut")),
0032           coshEtaLUT2_(lutConfig.getParameterSet("cosh_eta_lut2")),
0033           cosPhiLUT_(lutConfig.getParameterSet("cos_phi_lut")),
0034           minDEta_(getOptionalParam<int, double>(
0035               "minDEta", config, [&scales](double value) { return scales.to_hw_eta(value); })),
0036           maxDEta_(getOptionalParam<int, double>(
0037               "maxDEta", config, [&scales](double value) { return scales.to_hw_eta(value); })),
0038           minDPhi_(getOptionalParam<int, double>(
0039               "minDPhi", config, [&scales](double value) { return scales.to_hw_phi(value); })),
0040           maxDPhi_(getOptionalParam<int, double>(
0041               "maxDPhi", config, [&scales](double value) { return scales.to_hw_phi(value); })),
0042           minDz_(getOptionalParam<int, double>(
0043               "minDz", config, [&scales](double value) { return scales.to_hw_z0(value); })),
0044           maxDz_(getOptionalParam<int, double>(
0045               "maxDz", config, [&scales](double value) { return scales.to_hw_z0(value); })),
0046           minDRSquared_(getOptionalParam<int, double>(
0047               "minDR", config, [&scales](double value) { return scales.to_hw_dRSquared(value); })),
0048           maxDRSquared_(getOptionalParam<int, double>(
0049               "maxDR", config, [&scales](double value) { return scales.to_hw_dRSquared(value); })),
0050           minInvMassSqrDiv2_(getOptionalParam<double, double>(
0051               "minInvMass", config, [&scales](double value) { return scales.to_hw_InvMassSqrDiv2(value); })),
0052           maxInvMassSqrDiv2_(getOptionalParam<double, double>(
0053               "maxInvMass", config, [&scales](double value) { return scales.to_hw_InvMassSqrDiv2(value); })),
0054           minTransMassSqrDiv2_(getOptionalParam<double, double>(
0055               "minTransMass", config, [&scales](double value) { return scales.to_hw_TransMassSqrDiv2(value); })),
0056           maxTransMassSqrDiv2_(getOptionalParam<double, double>(
0057               "maxTransMass", config, [&scales](double value) { return scales.to_hw_TransMassSqrDiv2(value); })),
0058           minPTSquared_(getOptionalParam<double, double>(
0059               "minCombPt", config, [&scales](double value) { return scales.to_hw_PtSquared(value); })),
0060           maxPTSquared_(getOptionalParam<double, double>(
0061               "maxCombPt", config, [&scales](double value) { return scales.to_hw_PtSquared(value); })),
0062           os_(config.getParameter<bool>("os")),
0063           ss_(config.getParameter<bool>("ss")),
0064           enable_sanity_checks_(enable_sanity_checks),
0065           inv_mass_checks_(inv_mass_checks) {}
0066 
0067     bool checkObjects(const P2GTCandidate& obj1,
0068                       const P2GTCandidate& obj2,
0069                       InvariantMassErrorCollection& massErrors) const {
0070       bool res = true;
0071 
0072       std::optional<uint32_t> dEta;
0073 
0074       if (minDEta_ || maxDEta_ || minDRSquared_ || maxDRSquared_ || minInvMassSqrDiv2_ || maxInvMassSqrDiv2_) {
0075         dEta = (obj1.hwEta() > obj2.hwEta()) ? obj1.hwEta().to_int() - obj2.hwEta().to_int()
0076                                              : obj2.hwEta().to_int() - obj1.hwEta().to_int();
0077         res &= minDEta_ ? dEta > minDEta_ : true;
0078         res &= maxDEta_ ? dEta < maxDEta_ : true;
0079       }
0080 
0081       constexpr int HW_PI = 1 << (P2GTCandidate::hwPhi_t::width - 1);  // assumes phi in [-pi, pi)
0082 
0083       // Ensure dPhi is always the smaller angle, i.e. always between [0, pi]
0084       std::optional<uint32_t> dPhi;
0085 
0086       if (minDPhi_ || maxDPhi_ || minDRSquared_ || maxDRSquared_ || minInvMassSqrDiv2_ || maxInvMassSqrDiv2_ ||
0087           minTransMassSqrDiv2_ || maxTransMassSqrDiv2_ || minPTSquared_ || maxPTSquared_) {
0088         dPhi = HW_PI - abs(abs(obj1.hwPhi().to_int() - obj2.hwPhi().to_int()) - HW_PI);
0089       }
0090 
0091       res &= minDPhi_ ? dPhi > minDPhi_ : true;
0092       res &= maxDPhi_ ? dPhi < maxDPhi_ : true;
0093 
0094       if (minDz_ || maxDz_) {
0095         uint32_t dZ = abs(obj1.hwZ0() - obj2.hwZ0());
0096         res &= minDz_ ? dZ > minDz_ : true;
0097         res &= maxDz_ ? dZ < maxDz_ : true;
0098       }
0099 
0100       if (minDRSquared_ || maxDRSquared_) {
0101         uint32_t dRSquared = dEta.value() * dEta.value() + dPhi.value() * dPhi.value();
0102         res &= minDRSquared_ ? dRSquared > minDRSquared_ : true;
0103         res &= maxDRSquared_ ? dRSquared < maxDRSquared_ : true;
0104       }
0105 
0106       res &= os_ ? obj1.hwCharge() != obj2.hwCharge() : true;
0107       res &= ss_ ? obj1.hwCharge() == obj2.hwCharge() : true;
0108 
0109       int32_t lutCoshDEta = 0;
0110       if (dEta) {
0111         lutCoshDEta = dEta < DETA_LUT_SPLIT ? coshEtaLUT_[dEta.value()] : coshEtaLUT2_[dEta.value() - DETA_LUT_SPLIT];
0112       }
0113 
0114       // Optimization: (cos(x + pi) = -cos(x), x in [0, pi))
0115       int32_t lutCosDPhi = 0;
0116       if (dPhi) {
0117         lutCosDPhi = dPhi >= HW_PI ? -cosPhiLUT_[dPhi.value()] : cosPhiLUT_[dPhi.value()];
0118       }
0119 
0120       if (enable_sanity_checks_ && dEta && dPhi) {
0121         // Check whether the LUT error is smaller or equal than the expected maximum LUT error
0122         double coshEtaLUTMax = dEta < DETA_LUT_SPLIT ? coshEtaLUT_.hwMax_error() : coshEtaLUT2_.hwMax_error();
0123         double etaLUTScale = dEta < DETA_LUT_SPLIT ? coshEtaLUT_.output_scale() : coshEtaLUT2_.output_scale();
0124 
0125         if (std::abs(lutCoshDEta - etaLUTScale * std::cosh(dEta.value() * scales_.eta_lsb())) > coshEtaLUTMax) {
0126           edm::LogError("COSH LUT") << "Difference larger than max LUT error: " << coshEtaLUTMax
0127                                     << ", lut: " << lutCoshDEta
0128                                     << ", calc: " << etaLUTScale * std::cosh(dEta.value() * scales_.eta_lsb())
0129                                     << ", dEta: " << dEta.value() << ", scale: " << etaLUTScale;
0130         }
0131 
0132         if (std::abs(lutCosDPhi - cosPhiLUT_.output_scale() * std::cos(dPhi.value() * scales_.phi_lsb())) >
0133             cosPhiLUT_.hwMax_error()) {
0134           edm::LogError("COS LUT") << "Difference larger than max LUT error: " << cosPhiLUT_.hwMax_error()
0135                                    << ", lut: " << lutCosDPhi << ", calc: "
0136                                    << cosPhiLUT_.output_scale() * std::cos(dPhi.value() * scales_.phi_lsb());
0137         }
0138       }
0139 
0140       if (minInvMassSqrDiv2_ || maxInvMassSqrDiv2_) {
0141         int64_t invMassSqrDiv2;
0142         if (dEta < DETA_LUT_SPLIT) {
0143           // dEta [0, 2pi)
0144           invMassSqrDiv2 = obj1.hwPT().to_int64() * obj2.hwPT().to_int64() * (lutCoshDEta - lutCosDPhi);
0145           res &= minInvMassSqrDiv2_
0146                      ? invMassSqrDiv2 > std::round(minInvMassSqrDiv2_.value() * coshEtaLUT_.output_scale())
0147                      : true;
0148           res &= maxInvMassSqrDiv2_
0149                      ? invMassSqrDiv2 < std::round(maxInvMassSqrDiv2_.value() * coshEtaLUT_.output_scale())
0150                      : true;
0151         } else {
0152           // dEta [2pi, 4pi), ignore cos
0153           invMassSqrDiv2 = obj1.hwPT().to_int64() * obj2.hwPT().to_int64() * lutCoshDEta;
0154           res &= minInvMassSqrDiv2_
0155                      ? invMassSqrDiv2 > std::round(minInvMassSqrDiv2_.value() * coshEtaLUT2_.output_scale())
0156                      : true;
0157           res &= maxInvMassSqrDiv2_
0158                      ? invMassSqrDiv2 < std::round(maxInvMassSqrDiv2_.value() * coshEtaLUT2_.output_scale())
0159                      : true;
0160         }
0161 
0162         if (inv_mass_checks_) {
0163           double precInvMass =
0164               scales_.pT_lsb() *
0165               std::sqrt(2 * obj1.hwPT().to_double() * obj2.hwPT().to_double() *
0166                         (std::cosh(dEta.value() * scales_.eta_lsb()) - std::cos(dPhi.value() * scales_.phi_lsb())));
0167 
0168           double lutInvMass = scales_.pT_lsb() * std::sqrt(2 * static_cast<double>(invMassSqrDiv2) /
0169                                                            (dEta < DETA_LUT_SPLIT ? coshEtaLUT_.output_scale()
0170                                                                                   : coshEtaLUT2_.output_scale()));
0171 
0172           double error = std::abs(precInvMass - lutInvMass);
0173           massErrors.emplace_back(InvariantMassError{error, error / precInvMass, precInvMass});
0174         }
0175       }
0176 
0177       if (minPTSquared_ || maxPTSquared_) {
0178         int64_t pTSquared = obj1.hwPT().to_int64() * obj1.hwPT().to_int64() *
0179                                 static_cast<int64_t>(std::round(cosPhiLUT_.output_scale())) +
0180                             obj2.hwPT().to_int64() * obj2.hwPT().to_int64() *
0181                                 static_cast<int64_t>(std::round(cosPhiLUT_.output_scale())) +
0182                             2 * obj1.hwPT().to_int64() * obj2.hwPT().to_int64() * lutCosDPhi;
0183         res &= minPTSquared_ ? pTSquared > std::round(minPTSquared_.value() * cosPhiLUT_.output_scale()) : true;
0184         res &= maxPTSquared_ ? pTSquared < std::round(maxPTSquared_.value() * cosPhiLUT_.output_scale()) : true;
0185       }
0186 
0187       if (minTransMassSqrDiv2_ || maxTransMassSqrDiv2_) {
0188         int64_t transMassDiv2 = obj1.hwPT().to_int64() * obj2.hwPT().to_int64() *
0189                                 (static_cast<int64_t>(coshEtaLUT_.output_scale()) - lutCosDPhi);
0190         res &= minTransMassSqrDiv2_
0191                    ? transMassDiv2 > std::round(minTransMassSqrDiv2_.value() * coshEtaLUT_.output_scale())
0192                    : true;
0193         res &= maxTransMassSqrDiv2_
0194                    ? transMassDiv2 < std::round(maxTransMassSqrDiv2_.value() * coshEtaLUT_.output_scale())
0195                    : true;
0196       }
0197 
0198       return res;
0199     }
0200 
0201     static void fillLUTDescriptions(edm::ParameterSetDescription& desc) {
0202       edm::ParameterSetDescription coshLUTDesc;
0203       L1GTSingleInOutLUT::fillLUTDescriptions(coshLUTDesc);
0204       desc.add<edm::ParameterSetDescription>("cosh_eta_lut", coshLUTDesc);
0205 
0206       edm::ParameterSetDescription coshLUT2Desc;
0207       L1GTSingleInOutLUT::fillLUTDescriptions(coshLUT2Desc);
0208       desc.add<edm::ParameterSetDescription>("cosh_eta_lut2", coshLUT2Desc);
0209 
0210       edm::ParameterSetDescription cosLUTDesc;
0211       L1GTSingleInOutLUT::fillLUTDescriptions(cosLUTDesc);
0212       desc.add<edm::ParameterSetDescription>("cos_phi_lut", cosLUTDesc);
0213     }
0214 
0215     static void fillPSetDescription(edm::ParameterSetDescription& desc) {
0216       desc.addOptional<double>("minDEta");
0217       desc.addOptional<double>("maxDEta");
0218       desc.addOptional<double>("minDPhi");
0219       desc.addOptional<double>("maxDPhi");
0220       desc.addOptional<double>("minDR");
0221       desc.addOptional<double>("maxDR");
0222       desc.addOptional<double>("minDz");
0223       desc.addOptional<double>("maxDz");
0224       desc.addOptional<double>("minInvMass");
0225       desc.addOptional<double>("maxInvMass");
0226       desc.addOptional<double>("minTransMass");
0227       desc.addOptional<double>("maxTransMass");
0228       desc.addOptional<double>("minCombPt");
0229       desc.addOptional<double>("maxCombPt");
0230       desc.add<bool>("os", false);
0231       desc.add<bool>("ss", false);
0232     }
0233 
0234   private:
0235     const L1GTScales& scales_;
0236 
0237     const L1GTSingleInOutLUT coshEtaLUT_;   // [0, 2pi)
0238     const L1GTSingleInOutLUT coshEtaLUT2_;  // [2pi, 4pi)
0239     const L1GTSingleInOutLUT cosPhiLUT_;
0240 
0241     const std::optional<int> minDEta_;
0242     const std::optional<int> maxDEta_;
0243     const std::optional<int> minDPhi_;
0244     const std::optional<int> maxDPhi_;
0245     const std::optional<int> minDz_;
0246     const std::optional<int> maxDz_;
0247 
0248     const std::optional<int> minDRSquared_;
0249     const std::optional<int> maxDRSquared_;
0250 
0251     const std::optional<double> minInvMassSqrDiv2_;
0252     const std::optional<double> maxInvMassSqrDiv2_;
0253     const std::optional<double> minTransMassSqrDiv2_;
0254     const std::optional<double> maxTransMassSqrDiv2_;
0255 
0256     const std::optional<double> minPTSquared_;
0257     const std::optional<double> maxPTSquared_;
0258 
0259     const bool os_;  // Opposite sign
0260     const bool ss_;  // Same sign
0261 
0262     const bool enable_sanity_checks_;
0263     const bool inv_mass_checks_;
0264   };
0265 
0266 }  // namespace l1t
0267 
0268 #endif  // L1Trigger_Phase2L1GT_L1GTDeltaCut_h