Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-06-11 03:00:12

0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/puppi/linpuppi_ref.h"
0002 #include "L1Trigger/Phase2L1ParticleFlow/interface/puppi/linpuppi_bits.h"
0003 #include <cmath>
0004 #include <algorithm>
0005 
0006 #include "L1Trigger/Phase2L1ParticleFlow/interface/common/bitonic_hybrid_sort_ref.h"
0007 #include "L1Trigger/Phase2L1ParticleFlow/interface/common/bitonic_sort_ref.h"
0008 #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h"
0009 
0010 #ifdef CMSSW_GIT_HASH
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0013 #include "FWCore/ParameterSet/interface/allowedValues.h"
0014 #include "FWCore/Utilities/interface/transform.h"
0015 #include "FWCore/Utilities/interface/Exception.h"
0016 #include "FWCore/ParameterSet/interface/FileInPath.h"
0017 #endif
0018 
0019 #include "L1Trigger/Phase2L1ParticleFlow/interface/NNVtxAssoc.h"
0020 
0021 using namespace l1ct;
0022 using namespace linpuppi;
0023 
0024 l1ct::LinPuppiEmulator::LinPuppiEmulator(unsigned int nTrack,
0025                                          unsigned int nIn,
0026                                          unsigned int nOut,
0027                                          unsigned int nVtx,
0028                                          unsigned int dR2Min,
0029                                          unsigned int dR2Max,
0030                                          unsigned int iptMax,
0031                                          unsigned int dzCut,
0032                                          glbeta_t etaCut,
0033                                          double ptSlopeNe_0,
0034                                          double ptSlopeNe_1,
0035                                          double ptSlopePh_0,
0036                                          double ptSlopePh_1,
0037                                          double ptZeroNe_0,
0038                                          double ptZeroNe_1,
0039                                          double ptZeroPh_0,
0040                                          double ptZeroPh_1,
0041                                          double alphaSlope_0,
0042                                          double alphaSlope_1,
0043                                          double alphaZero_0,
0044                                          double alphaZero_1,
0045                                          double alphaCrop_0,
0046                                          double alphaCrop_1,
0047                                          double priorNe_0,
0048                                          double priorNe_1,
0049                                          double priorPh_0,
0050                                          double priorPh_1,
0051                                          pt_t ptCut_0,
0052                                          pt_t ptCut_1,
0053                                          bool useMLAssociation,
0054                                          const double associationThreshold,
0055                                          std::string associationGraphPath,
0056                                          std::vector<double> associationNetworkZ0binning,
0057                                          std::vector<double> associationNetworkEtaBounds,
0058                                          std::vector<double> associationNetworkZ0ResBins,
0059                                          unsigned int nFinalSort,
0060                                          SortAlgo finalSortAlgo)
0061     : nTrack_(nTrack),
0062       nIn_(nIn),
0063       nOut_(nOut),
0064       nVtx_(nVtx),
0065       dR2Min_(dR2Min),
0066       dR2Max_(dR2Max),
0067       iptMax_(iptMax),
0068       dzCut_(dzCut),
0069       absEtaBins_(1, etaCut),
0070       ptSlopeNe_(2),
0071       ptSlopePh_(2),
0072       ptZeroNe_(2),
0073       ptZeroPh_(2),
0074       alphaSlope_(2),
0075       alphaZero_(2),
0076       alphaCrop_(2),
0077       priorNe_(2),
0078       priorPh_(2),
0079       ptCut_(2),
0080       useMLAssociation_(useMLAssociation),
0081       nFinalSort_(nFinalSort ? nFinalSort : nOut),
0082       finalSortAlgo_(finalSortAlgo),
0083       debug_(false),
0084       fakePuppi_(false) {
0085   ptSlopeNe_[0] = ptSlopeNe_0;
0086   ptSlopeNe_[1] = ptSlopeNe_1;
0087   ptSlopePh_[0] = ptSlopePh_0;
0088   ptSlopePh_[1] = ptSlopePh_1;
0089   ptZeroNe_[0] = ptZeroNe_0;
0090   ptZeroNe_[1] = ptZeroNe_1;
0091   ptZeroPh_[0] = ptZeroPh_0;
0092   ptZeroPh_[1] = ptZeroPh_1;
0093   alphaSlope_[0] = alphaSlope_0;
0094   alphaSlope_[1] = alphaSlope_1;
0095   alphaZero_[0] = alphaZero_0;
0096   alphaZero_[1] = alphaZero_1;
0097   alphaCrop_[0] = alphaCrop_0;
0098   alphaCrop_[1] = alphaCrop_1;
0099   priorNe_[0] = priorNe_0;
0100   priorNe_[1] = priorNe_1;
0101   priorPh_[0] = priorPh_0;
0102   priorPh_[1] = priorPh_1;
0103   ptCut_[0] = ptCut_0;
0104   ptCut_[1] = ptCut_1;
0105 
0106   if (useMLAssociation_ and withinCMSSW_) {
0107     nnVtxAssoc_ = std::make_unique<NNVtxAssoc>(NNVtxAssoc(associationGraphPath,
0108                                                           associationThreshold,
0109                                                           associationNetworkZ0binning,
0110                                                           associationNetworkEtaBounds,
0111                                                           associationNetworkZ0ResBins));
0112   }
0113 }
0114 
0115 #ifdef CMSSW_GIT_HASH
0116 l1ct::LinPuppiEmulator::LinPuppiEmulator(const edm::ParameterSet &iConfig)
0117     : nTrack_(iConfig.getParameter<uint32_t>("nTrack")),
0118       nIn_(iConfig.getParameter<uint32_t>("nIn")),
0119       nOut_(iConfig.getParameter<uint32_t>("nOut")),
0120       nVtx_(iConfig.getParameter<uint32_t>("nVtx")),
0121       dR2Min_(l1ct::Scales::makeDR2FromFloatDR(iConfig.getParameter<double>("drMin"))),
0122       dR2Max_(l1ct::Scales::makeDR2FromFloatDR(iConfig.getParameter<double>("dr"))),
0123       iptMax_(l1ct::Scales::intPt(l1ct::Scales::makePtFromFloat(iConfig.getParameter<double>("ptMax")))),
0124       dzCut_(l1ct::Scales::makeZ0(iConfig.getParameter<double>("dZ"))),
0125       absEtaBins_(
0126           edm::vector_transform(iConfig.getParameter<std::vector<double>>("absEtaCuts"), l1ct::Scales::makeGlbEta)),
0127       ptSlopeNe_(iConfig.getParameter<std::vector<double>>("ptSlopes")),
0128       ptSlopePh_(iConfig.getParameter<std::vector<double>>("ptSlopesPhoton")),
0129       ptZeroNe_(iConfig.getParameter<std::vector<double>>("ptZeros")),
0130       ptZeroPh_(iConfig.getParameter<std::vector<double>>("ptZerosPhoton")),
0131       alphaSlope_(iConfig.getParameter<std::vector<double>>("alphaSlopes")),
0132       alphaZero_(iConfig.getParameter<std::vector<double>>("alphaZeros")),
0133       alphaCrop_(iConfig.getParameter<std::vector<double>>("alphaCrop")),
0134       priorNe_(iConfig.getParameter<std::vector<double>>("priors")),
0135       priorPh_(iConfig.getParameter<std::vector<double>>("priorsPhoton")),
0136       ptCut_(edm::vector_transform(iConfig.getParameter<std::vector<double>>("ptCut"), l1ct::Scales::makePtFromFloat)),
0137       useMLAssociation_(iConfig.getParameter<bool>("useMLAssociation")),
0138       nFinalSort_(iConfig.getParameter<uint32_t>("nFinalSort")),
0139       debug_(iConfig.getUntrackedParameter<bool>("debug", false)),
0140       fakePuppi_(iConfig.getParameter<bool>("fakePuppi")) {
0141   if (absEtaBins_.size() + 1 != ptSlopeNe_.size())
0142     throw cms::Exception("Configuration", "size mismatch for ptSlopes parameter");
0143   if (absEtaBins_.size() + 1 != ptSlopePh_.size())
0144     throw cms::Exception("Configuration", "size mismatch for ptSlopesPhoton parameter");
0145   if (absEtaBins_.size() + 1 != ptZeroPh_.size())
0146     throw cms::Exception("Configuration", "size mismatch for ptZeros parameter");
0147   if (absEtaBins_.size() + 1 != ptZeroNe_.size())
0148     throw cms::Exception("Configuration", "size mismatch for ptZerosPhotons parameter");
0149   if (absEtaBins_.size() + 1 != priorPh_.size())
0150     throw cms::Exception("Configuration", "size mismatch for priors parameter");
0151   if (absEtaBins_.size() + 1 != priorNe_.size())
0152     throw cms::Exception("Configuration", "size mismatch for priorsPhotons parameter");
0153   if (absEtaBins_.size() + 1 != alphaSlope_.size())
0154     throw cms::Exception("Configuration", "size mismatch for alphaSlope parameter");
0155   if (absEtaBins_.size() + 1 != alphaZero_.size())
0156     throw cms::Exception("Configuration", "size mismatch for alphaZero parameter");
0157   if (absEtaBins_.size() + 1 != alphaCrop_.size())
0158     throw cms::Exception("Configuration", "size mismatch for alphaCrop parameter");
0159   if (absEtaBins_.size() + 1 != ptCut_.size())
0160     throw cms::Exception("Configuration", "size mismatch for ptCut parameter");
0161   if (useMLAssociation_) {
0162     edm::ParameterSet nnVtxAssocPSet_ = iConfig.getParameter<edm::ParameterSet>("NNVtxAssociation");
0163     edm::FileInPath associationGraphPathFIP =
0164         edm::FileInPath(nnVtxAssocPSet_.getParameter<std::string>("associationGraph"));
0165 
0166     nnVtxAssoc_ = std::make_unique<NNVtxAssoc>(
0167         NNVtxAssoc(associationGraphPathFIP.fullPath(),
0168                    nnVtxAssocPSet_.getParameter<double>("associationThreshold"),
0169                    nnVtxAssocPSet_.getParameter<std::vector<double>>("associationNetworkZ0binning"),
0170                    nnVtxAssocPSet_.getParameter<std::vector<double>>("associationNetworkEtaBounds"),
0171                    nnVtxAssocPSet_.getParameter<std::vector<double>>("associationNetworkZ0ResBins")));
0172   }
0173   const std::string &sortAlgo = iConfig.getParameter<std::string>("finalSortAlgo");
0174   if (sortAlgo == "Insertion")
0175     finalSortAlgo_ = SortAlgo::Insertion;
0176   else if (sortAlgo == "BitonicRUFL")
0177     finalSortAlgo_ = SortAlgo::BitonicRUFL;
0178   else if (sortAlgo == "BitonicHLS")
0179     finalSortAlgo_ = SortAlgo::BitonicHLS;
0180   else if (sortAlgo == "Hybrid")
0181     finalSortAlgo_ = SortAlgo::Hybrid;
0182   else if (sortAlgo == "FoldedHybrid")
0183     finalSortAlgo_ = SortAlgo::FoldedHybrid;
0184   else
0185     throw cms::Exception("Configuration", "unsupported finalSortAlgo '" + sortAlgo + "'");
0186 }
0187 
0188 edm::ParameterSetDescription l1ct::LinPuppiEmulator::getParameterSetDescription() {
0189   edm::ParameterSetDescription description;
0190   description.add<uint32_t>("nTrack");
0191   description.add<uint32_t>("nIn");
0192   description.add<uint32_t>("nOut");
0193   description.add<uint32_t>("nVtx", 1);
0194   description.add<double>("dZ");
0195   description.add<double>("dr");
0196   description.add<double>("drMin");
0197   description.add<double>("ptMax");
0198   description.add<std::vector<double>>("absEtaCuts");
0199   description.add<std::vector<double>>("ptCut");
0200   description.add<bool>("useMLAssociation");
0201   description.add<edm::ParameterSetDescription>("NNVtxAssociation", NNVtxAssoc::getParameterSetDescription());
0202   description.add<std::vector<double>>("ptSlopes");
0203   description.add<std::vector<double>>("ptSlopesPhoton");
0204   description.add<std::vector<double>>("ptZeros");
0205   description.add<std::vector<double>>("ptZerosPhoton");
0206   description.add<std::vector<double>>("alphaSlopes");
0207   description.add<std::vector<double>>("alphaZeros");
0208   description.add<std::vector<double>>("alphaCrop");
0209   description.add<std::vector<double>>("priors");
0210   description.add<std::vector<double>>("priorsPhoton");
0211   description.add<uint32_t>("nFinalSort");
0212   description.ifValue(
0213       edm::ParameterDescription<std::string>("finalSortAlgo", "Insertion", true),
0214       edm::allowedValues<std::string>("Insertion", "BitonicRUFL", "BitonicHLS", "Hybrid", "FoldedHybrid"));
0215   description.add<bool>("fakePuppi", false);
0216   description.addUntracked<bool>("debug", false);
0217   return description;
0218 }
0219 #endif
0220 
0221 void l1ct::LinPuppiEmulator::puppisort_and_crop_ref(unsigned int nOutMax,
0222                                                     const std::vector<l1ct::PuppiObjEmu> &in,
0223                                                     std::vector<l1ct::PuppiObjEmu> &out,
0224                                                     SortAlgo sortAlgo) {
0225   const unsigned int nOut = std::min<unsigned int>(nOutMax, in.size());
0226   out.resize(nOut);
0227   for (unsigned int iout = 0; iout < nOut; ++iout) {
0228     out[iout].clear();
0229   }
0230 
0231   if (sortAlgo == SortAlgo::Insertion) {
0232     for (unsigned int it = 0, nIn = in.size(); it < nIn; ++it) {
0233       for (int iout = int(nOut) - 1; iout >= 0; --iout) {
0234         if (out[iout].hwPt <= in[it].hwPt) {
0235           if (iout == 0 || out[iout - 1].hwPt > in[it].hwPt) {
0236             out[iout] = in[it];
0237           } else {
0238             out[iout] = out[iout - 1];
0239           }
0240         }
0241       }
0242     }
0243   } else if (sortAlgo == SortAlgo::BitonicRUFL) {
0244     bitonic_sort_and_crop_ref(in.size(), nOut, &in[0], &out[0]);
0245   } else if (sortAlgo == SortAlgo::BitonicHLS || sortAlgo == SortAlgo::Hybrid) {
0246     hybrid_bitonic_sort_and_crop_ref(in.size(), nOut, &in[0], &out[0], sortAlgo == SortAlgo::Hybrid);
0247   } else if (sortAlgo == SortAlgo::FoldedHybrid) {
0248     folded_hybrid_bitonic_sort_and_crop_ref(in.size(), nOut, &in[0], &out[0], true);
0249   }
0250 }
0251 
0252 void l1ct::LinPuppiEmulator::linpuppi_chs_ref(const PFRegionEmu &region,
0253                                               const std::vector<PVObjEmu> &pv,
0254                                               const std::vector<PFChargedObjEmu> &pfch /*[nTrack]*/,
0255                                               std::vector<PuppiObjEmu> &outallch /*[nTrack]*/) const {
0256   const unsigned int nTrack = std::min<unsigned int>(nTrack_, pfch.size());
0257   const unsigned int nVtx = std::min<unsigned int>(nVtx_, pv.size());
0258   outallch.resize(nTrack);
0259   for (unsigned int i = 0; i < nTrack; ++i) {
0260     int pZ0 = pfch[i].hwZ0;
0261     int z0diff = -99999;
0262     bool pass_network = false;
0263     for (unsigned int j = 0; j < nVtx; ++j) {
0264       int pZ0Diff = pZ0 - pv[j].hwZ0;
0265       if (std::abs(z0diff) > std::abs(pZ0Diff))
0266         z0diff = pZ0Diff;
0267       if (useMLAssociation_ and withinCMSSW_ &&
0268           nnVtxAssoc_->TTTrackNetworkSelector<const l1ct::PFChargedObjEmu>(region, pfch[i], pv[j]) == 1)
0269         pass_network = true;
0270     }
0271     bool accept = pfch[i].hwPt != 0;
0272     if (!fakePuppi_ && !useMLAssociation_)
0273       accept = accept && region.isFiducial(pfch[i]) && (std::abs(z0diff) <= int(dzCut_) || pfch[i].hwId.isMuon());
0274     if (!fakePuppi_ && useMLAssociation_)
0275       accept = accept && region.isFiducial(pfch[i]) && (pass_network || pfch[i].hwId.isMuon());
0276     if (accept) {
0277       outallch[i].fill(region, pfch[i]);
0278       if (fakePuppi_) {                           // overwrite Dxy & TkQuality with debug information
0279         outallch[i].setHwDxy(dxy_t(pv[0].hwZ0));  ///hack to get this to work
0280         outallch[i].setHwTkQuality(region.isFiducial(pfch[i]) ? 1 : 0);
0281       }
0282       if (debug_ && pfch[i].hwPt > 0)
0283         dbgPrintf("ref candidate %02u pt %7.2f pid %1d   vz %+6d  dz %+6d (cut %5d), fid %1d -> pass, packed %s\n",
0284                   i,
0285                   pfch[i].floatPt(),
0286                   pfch[i].intId(),
0287                   int(pfch[i].hwZ0),
0288                   z0diff,
0289                   dzCut_,
0290                   region.isFiducial(pfch[i]),
0291                   outallch[i].pack().to_string(16).c_str());
0292     } else {
0293       outallch[i].clear();
0294       if (debug_ && pfch[i].hwPt > 0)
0295         dbgPrintf("ref candidate %02u pt %7.2f pid %1d   vz %+6d  dz %+6d (cut %5d), fid %1d -> fail\n",
0296                   i,
0297                   pfch[i].floatPt(),
0298                   pfch[i].intId(),
0299                   int(pfch[i].hwZ0),
0300                   z0diff,
0301                   dzCut_,
0302                   region.isFiducial(pfch[i]));
0303     }
0304   }
0305 }
0306 
0307 unsigned int l1ct::LinPuppiEmulator::find_ieta(const PFRegionEmu &region, eta_t eta) const {
0308   int n = absEtaBins_.size();
0309   glbeta_t abseta = region.hwGlbEta(eta);
0310   if (abseta < 0)
0311     abseta = -abseta;
0312   for (int i = 0; i < n; ++i) {
0313     if (abseta <= absEtaBins_[i])
0314       return i;
0315   }
0316   return n;
0317 }
0318 
0319 std::pair<pt_t, puppiWgt_t> l1ct::LinPuppiEmulator::sum2puppiPt_ref(
0320     sumTerm_t sum, pt_t pt, unsigned int ieta, bool isEM, int icand) const {
0321   const alpha_t logOffset = std::log2(linpuppi::PT2DR2_LSB) - SUM_BITSHIFT;
0322   const ptSlope_t ptSlopeNe = ptSlopeNe_[ieta];
0323   const ptSlope_t ptSlopePh = ptSlopePh_[ieta];
0324   const l1ct::pt_t ptZeroNe = ptZeroNe_[ieta];
0325   const l1ct::pt_t ptZeroPh = ptZeroPh_[ieta];
0326   const x2_t alphaCrop = alphaCrop_[ieta];
0327   // we put a log(2) in all alphaSlopes here since we compute alpha as log2(sum) instead of ln(sum)
0328   const alphaSlope_t alphaSlopeNe = alphaSlope_[ieta] * std::log(2.);
0329   const alphaSlope_t alphaSlopePh = alphaSlope_[ieta] * std::log(2.);
0330   const alpha_t alphaZeroNe = alphaZero_[ieta] / std::log(2.);
0331   const alpha_t alphaZeroPh = alphaZero_[ieta] / std::log(2.);
0332   const x2_t priorNe = priorNe_[ieta];
0333   const x2_t priorPh = priorPh_[ieta];
0334 
0335   // emulate computing
0336   //    alpha = log2(sum)
0337   //    x2a = alphaSlope*(alpha - alphaZero)
0338   // we use a 10-bit LUT for the log2(sum), and to save firmware resources we
0339   // also pack the computation of x2a in the same LUT.
0340   // In this emulator, however, we also compute the separately alpha, as it is
0341   // useful for debugging and comparison to the floating-point emulator.
0342   const int log2lut_bits = 10;
0343   alpha_t alpha = 0;
0344   uint64_t logarg = sum.bits_to_uint64();
0345   if (logarg > 0) {
0346     alpha = logOffset;
0347     while (logarg >= (1 << log2lut_bits)) {
0348       logarg = logarg >> 1;
0349       alpha += 1;
0350     }
0351     alpha += alpha_t(std::log2(float(logarg)));
0352   }
0353   alpha_t alphaZero = (isEM ? alphaZeroPh : alphaZeroNe);
0354   alphaSlope_t alphaSlope = (isEM ? alphaSlopePh : alphaSlopeNe);
0355   x2_t x2a = -alphaSlope * alphaZero;
0356   logarg = sum.bits_to_uint64();
0357   if (logarg > 0) {
0358     x2a += alphaSlope * logOffset;
0359     while (logarg >= (1 << log2lut_bits)) {
0360       logarg = logarg >> 1;
0361       x2a += alphaSlope;
0362     }
0363     x2a += alphaSlope * alpha_t(std::log2(float(logarg)));
0364   }
0365   /* // debug printout, can be useful to study ranges and precisions for LUTs and coefficients
0366     dbgPrintf("ref:  x2a(sum = %10.5f, raw = %9lu): logarg = %9lu, sumterm = %.6f, table[logarg] = %d = %.6f, ret pre-crop = %d = %.6f\n", 
0367           sum.to_float(), sum.bits_to_uint64(), logarg, 
0368           (alphaSlope * (logOffset - alphaZero)).to_float(),
0369           (alphaSlope * alpha_t(std::log2(float(logarg)))).V.to_int(), 
0370           (alphaSlope * alpha_t(std::log2(float(logarg)))).to_float(), 
0371           x2a.V.to_int(), x2a.to_float());
0372   */
0373   x2a = std::min(std::max<x2_t>(x2a, -alphaCrop), alphaCrop);
0374 
0375   l1ct::pt_t ptZero = (isEM ? ptZeroPh : ptZeroNe);
0376   ptSlope_t ptSlope = (isEM ? ptSlopePh : ptSlopeNe);
0377   x2_t x2pt = ptSlope * (pt - ptZero);
0378 
0379   x2_t prior = (isEM ? priorPh : priorNe);
0380 
0381   x2_t x2 = x2a + x2pt - prior;
0382 
0383   puppiWgt_t weight = 1.0 / (1.0 + std::exp(-x2.to_float()));
0384   typedef ap_fixed<pt_t::width, pt_t::iwidth, AP_RND, AP_SAT> pt_rounding_t;
0385   pt_t ptPuppi = pt_rounding_t(pt * weight);
0386 
0387   if (debug_)
0388     dbgPrintf(
0389         "ref candidate %02d pt %7.2f  em %1d  ieta %1d: sum %10.4f alpha %+7.2f   x2a %+7.3f  x2pt %+7.3f   x2 %+7.3f  "
0390         "--> weight %.4f  puppi pt %7.2f\n",
0391         icand,
0392         Scales::floatPt(pt),
0393         int(isEM),
0394         ieta,
0395         sum.to_float(),
0396         alpha.to_float() * std::log(2.),
0397         x2a.to_float(),
0398         x2pt.to_float(),
0399         x2.to_float(),
0400         weight.to_float(),
0401         Scales::floatPt(ptPuppi));
0402 
0403   return std::make_pair(ptPuppi, weight);
0404 }
0405 
0406 void l1ct::LinPuppiEmulator::fwdlinpuppi_ref(const PFRegionEmu &region,
0407                                              const std::vector<HadCaloObjEmu> &caloin /*[nIn]*/,
0408                                              std::vector<PuppiObjEmu> &outallne_nocut /*[nIn]*/,
0409                                              std::vector<PuppiObjEmu> &outallne /*[nIn]*/,
0410                                              std::vector<PuppiObjEmu> &outselne /*[nOut]*/) const {
0411   const unsigned int nIn = std::min<unsigned int>(nIn_, caloin.size());
0412   const unsigned int PTMAX2 = (iptMax_ * iptMax_);
0413 
0414   outallne_nocut.resize(nIn);
0415   outallne.resize(nIn);
0416   for (unsigned int in = 0; in < nIn; ++in) {
0417     outallne_nocut[in].clear();
0418     outallne[in].clear();
0419     if (caloin[in].hwPt == 0)
0420       continue;
0421     sumTerm_t sum = 0;  // (pt^2)/(dr2)
0422     for (unsigned int it = 0; it < nIn; ++it) {
0423       if (it == in || caloin[it].hwPt == 0)
0424         continue;
0425       unsigned int dr2 = dr2_int(caloin[it].hwEta, caloin[it].hwPhi, caloin[in].hwEta, caloin[in].hwPhi);
0426       if (dr2 <= dR2Max_) {                                             // if dr is inside puppi cone
0427         unsigned int dr2short = (dr2 >= dR2Min_ ? dr2 : dR2Min_) >> 5;  // reduce precision to make divide LUT cheaper
0428         unsigned int pt = caloin[it].intPt(), pt2 = pt * pt;
0429         dr2inv_t dr2inv = 1.0f / float(dr2short);
0430         sumTerm_t term = std::min(pt2 >> 5, PTMAX2 >> 5) * dr2inv;
0431         sum += term;
0432       }
0433     }
0434     unsigned int ieta = find_ieta(region, caloin[in].hwEta);
0435     std::pair<pt_t, puppiWgt_t> ptAndW = sum2puppiPt_ref(sum, caloin[in].hwPt, ieta, caloin[in].hwIsEM(), in);
0436 
0437     outallne_nocut[in].fill(region, caloin[in], ptAndW.first, ptAndW.second);
0438     if (region.isFiducial(caloin[in]) && outallne_nocut[in].hwPt >= ptCut_[ieta]) {
0439       outallne[in] = outallne_nocut[in];
0440     }
0441   }
0442   puppisort_and_crop_ref(nOut_, outallne, outselne);
0443 }
0444 
0445 void l1ct::LinPuppiEmulator::linpuppi_ref(const PFRegionEmu &region,
0446                                           const std::vector<TkObjEmu> &track /*[nTrack]*/,
0447                                           const std::vector<PVObjEmu> &pv, /*[nVtx]*/
0448                                           const std::vector<PFNeutralObjEmu> &pfallne /*[nIn]*/,
0449                                           std::vector<PuppiObjEmu> &outallne_nocut /*[nIn]*/,
0450                                           std::vector<PuppiObjEmu> &outallne /*[nIn]*/,
0451                                           std::vector<PuppiObjEmu> &outselne /*[nOut]*/) const {
0452   const unsigned int nIn = std::min<unsigned>(nIn_, pfallne.size());
0453   const unsigned int nTrack = std::min<unsigned int>(nTrack_, track.size());
0454   const unsigned int PTMAX2 = (iptMax_ * iptMax_);
0455 
0456   outallne_nocut.resize(nIn);
0457   outallne.resize(nIn);
0458   for (unsigned int in = 0; in < nIn; ++in) {
0459     outallne_nocut[in].clear();
0460     outallne[in].clear();
0461     if (pfallne[in].hwPt == 0)
0462       continue;
0463     sumTerm_t sum = 0;  // (pt^2)/(dr2)
0464     for (unsigned int it = 0; it < nTrack; ++it) {
0465       if (track[it].hwPt == 0)
0466         continue;
0467 
0468       int pZMin = 99999;
0469       bool pass_network = false;
0470       for (unsigned int v = 0; v < nVtx_; ++v) {
0471         if (v < pv.size()) {
0472           int ppZMin = std::abs(int(track[it].hwZ0 - pv[v].hwZ0));
0473           if (pZMin > ppZMin)
0474             pZMin = ppZMin;
0475           if (useMLAssociation_ and withinCMSSW_ &&
0476               nnVtxAssoc_->TTTrackNetworkSelector<const l1ct::TkObjEmu>(region, track[it], pv[v]) == 1)
0477             pass_network = true;
0478         }
0479       }
0480       if (useMLAssociation_ && !pass_network)
0481         continue;
0482       if (!useMLAssociation_ && std::abs(pZMin) > int(dzCut_))
0483         continue;
0484       unsigned int dr2 = dr2_int(pfallne[in].hwEta, pfallne[in].hwPhi, track[it].hwEta, track[it].hwPhi);
0485       if (dr2 <= dR2Max_) {                                             // if dr is inside puppi cone
0486         unsigned int dr2short = (dr2 >= dR2Min_ ? dr2 : dR2Min_) >> 5;  // reduce precision to make divide LUT cheaper
0487         unsigned int pt = track[it].intPt(), pt2 = pt * pt;
0488         dr2inv_t dr2inv = 1.0f / float(dr2short);
0489         sumTerm_t term = std::min(pt2 >> 5, PTMAX2 >> 5) * dr2inv;
0490         /* // printout useful for comparing internals steps of computation to the ones done in the firmware code
0491         dbgPrintf("cand pT %5.1f eta %+6.2f ref term %2d %2d: dr = %8d  pt2_shift = %8lu  term = %12lu  apterm = %12.6f\n",
0492                   pfallne[it].floatPt(),
0493                   region.floatGlbEtaOf(pfallne[it]),
0494                   in,
0495                   it,
0496                   dr2,
0497                   std::min<uint64_t>(pt2 >> 5, PTMAX2 >> 5),
0498                   term.to_float());
0499         */
0500         sum += term;
0501       }
0502     }
0503 
0504     unsigned int ieta = find_ieta(region, pfallne[in].hwEta);
0505     bool isEM = (pfallne[in].hwId.isPhoton());
0506     std::pair<pt_t, puppiWgt_t> ptAndW = sum2puppiPt_ref(sum, pfallne[in].hwPt, ieta, isEM, in);
0507     if (!fakePuppi_) {
0508       outallne_nocut[in].fill(region, pfallne[in], ptAndW.first, ptAndW.second);
0509       if (region.isFiducial(pfallne[in]) && outallne_nocut[in].hwPt >= ptCut_[ieta]) {
0510         outallne[in] = outallne_nocut[in];
0511       }
0512     } else {  // fakePuppi: keep the full candidate, but set the Puppi weight and some debug info into it
0513       outallne_nocut[in].fill(region, pfallne[in], pfallne[in].hwPt, ptAndW.second);
0514       outallne_nocut[in].hwData[9] = region.isFiducial(pfallne[in]);
0515       outallne_nocut[in].hwData(20, 10) = ptAndW.first(10, 0);
0516       outallne[in] = outallne_nocut[in];
0517     }
0518     if (debug_ && pfallne[in].hwPt > 0 && outallne_nocut[in].hwPt > 0) {
0519       dbgPrintf("ref candidate %02u pt %7.2f  -> puppi pt %7.2f, fiducial %1d, packed %s\n",
0520                 in,
0521                 pfallne[in].floatPt(),
0522                 outallne_nocut[in].floatPt(),
0523                 int(region.isFiducial(pfallne[in])),
0524                 outallne_nocut[in].pack().to_string(16).c_str());
0525     }
0526   }
0527   puppisort_and_crop_ref(nOut_, outallne, outselne);
0528 }
0529 
0530 std::pair<float, float> l1ct::LinPuppiEmulator::sum2puppiPt_flt(
0531     float sum, float pt, unsigned int ieta, bool isEM, int icand) const {
0532   float alphaZero = alphaZero_[ieta], alphaSlope = alphaSlope_[ieta], alphaCrop = alphaCrop_[ieta];
0533   float alpha = sum > 0 ? std::log(sum) : -9e9;
0534   float x2a = std::min(std::max(alphaSlope * (alpha - alphaZero), -alphaCrop), alphaCrop);
0535 
0536   float ptZero = (isEM ? ptZeroPh_[ieta] : ptZeroNe_[ieta]);
0537   float ptSlope = (isEM ? ptSlopePh_[ieta] : ptSlopeNe_[ieta]);
0538   float x2pt = ptSlope * (pt - ptZero);
0539 
0540   float prior = (isEM ? priorPh_[ieta] : priorNe_[ieta]);
0541 
0542   float x2 = x2a + x2pt - prior;
0543 
0544   float weight = 1.0 / (1.0 + std::exp(-x2));
0545 
0546   float puppiPt = pt * weight;
0547   if (debug_)
0548     dbgPrintf(
0549         "flt candidate %02d pt %7.2f  em %1d  ieta %1d: alpha %+7.2f   x2a         %+7.3f  x2pt         %+7.3f   x2    "
0550         "     %+7.3f  --> weight        %.4f  puppi pt %7.2f\n",
0551         icand,
0552         pt,
0553         int(isEM),
0554         ieta,
0555         std::max(alpha, -99.99f),
0556         x2a,
0557         x2pt,
0558         x2,
0559         weight,
0560         puppiPt);
0561 
0562   return std::make_pair(puppiPt, weight);
0563 }
0564 
0565 void l1ct::LinPuppiEmulator::fwdlinpuppi_flt(const PFRegionEmu &region,
0566                                              const std::vector<HadCaloObjEmu> &caloin /*[nIn]*/,
0567                                              std::vector<PuppiObjEmu> &outallne_nocut /*[nIn]*/,
0568                                              std::vector<PuppiObjEmu> &outallne /*[nIn]*/,
0569                                              std::vector<PuppiObjEmu> &outselne /*[nOut]*/) const {
0570   const unsigned int nIn = std::min<unsigned int>(nIn_, caloin.size());
0571   const float f_ptMax = Scales::floatPt(Scales::makePt(iptMax_));
0572 
0573   outallne_nocut.resize(nIn);
0574   outallne.resize(nIn);
0575   for (unsigned int in = 0; in < nIn; ++in) {
0576     outallne_nocut[in].clear();
0577     outallne[in].clear();
0578     if (caloin[in].hwPt == 0)
0579       continue;
0580     float sum = 0;
0581     for (unsigned int it = 0; it < nIn; ++it) {
0582       if (it == in || caloin[it].hwPt == 0)
0583         continue;
0584       unsigned int dr2 = dr2_int(caloin[it].hwEta, caloin[it].hwPhi, caloin[in].hwEta, caloin[in].hwPhi);
0585       if (dr2 <= dR2Max_) {  // if dr is inside puppi cone
0586         sum += std::pow(std::min(caloin[it].floatPt(), f_ptMax), 2) / (std::max(dr2, dR2Min_) * DR2_LSB);
0587       }
0588     }
0589 
0590     unsigned int ieta = find_ieta(region, caloin[in].hwEta);
0591     std::pair<float, float> ptAndW = sum2puppiPt_flt(sum, caloin[in].floatPt(), ieta, caloin[in].hwIsEM(), in);
0592     outallne_nocut[in].fill(region, caloin[in], Scales::makePtFromFloat(ptAndW.first), l1ct::puppiWgt_t(ptAndW.second));
0593     if (region.isFiducial(caloin[in]) && outallne_nocut[in].hwPt >= ptCut_[ieta]) {
0594       outallne[in] = outallne_nocut[in];
0595     }
0596   }
0597 
0598   puppisort_and_crop_ref(nOut_, outallne, outselne);
0599 }
0600 
0601 void l1ct::LinPuppiEmulator::linpuppi_flt(const PFRegionEmu &region,
0602                                           const std::vector<TkObjEmu> &track /*[nTrack]*/,
0603                                           const std::vector<PVObjEmu> &pv,
0604                                           const std::vector<PFNeutralObjEmu> &pfallne /*[nIn]*/,
0605                                           std::vector<PuppiObjEmu> &outallne_nocut /*[nIn]*/,
0606                                           std::vector<PuppiObjEmu> &outallne /*[nIn]*/,
0607                                           std::vector<PuppiObjEmu> &outselne /*[nOut]*/) const {
0608   const unsigned int nIn = std::min<unsigned>(nIn_, pfallne.size());
0609   const unsigned int nTrack = std::min<unsigned int>(nTrack_, track.size());
0610   const float f_ptMax = Scales::floatPt(Scales::makePt(iptMax_));
0611 
0612   outallne_nocut.resize(nIn);
0613   outallne.resize(nIn);
0614   for (unsigned int in = 0; in < nIn; ++in) {
0615     outallne_nocut[in].clear();
0616     outallne[in].clear();
0617     if (pfallne[in].hwPt == 0)
0618       continue;
0619     float sum = 0;
0620     for (unsigned int it = 0; it < nTrack; ++it) {
0621       if (track[it].hwPt == 0)
0622         continue;
0623 
0624       int pZMin = 99999;
0625       bool pass_network = false;
0626       for (unsigned int v = 0, nVtx = std::min<unsigned int>(nVtx_, pv.size()); v < nVtx; ++v) {
0627         int ppZMin = std::abs(int(track[it].hwZ0 - pv[v].hwZ0));
0628         if (pZMin > ppZMin)
0629           pZMin = ppZMin;
0630         if (useMLAssociation_ and withinCMSSW_ &&
0631             nnVtxAssoc_->TTTrackNetworkSelector<const l1ct::TkObjEmu>(region, track[it], pv[v]) == 1)
0632           pass_network = true;
0633       }
0634       if (useMLAssociation_ && !pass_network)
0635         continue;
0636       if (!useMLAssociation_ && std::abs(pZMin) > int(dzCut_))
0637         continue;
0638       unsigned int dr2 = dr2_int(
0639           pfallne[in].hwEta, pfallne[in].hwPhi, track[it].hwEta, track[it].hwPhi);  // if dr is inside puppi cone
0640       if (dr2 <= dR2Max_) {
0641         sum += std::pow(std::min(track[it].floatPt(), f_ptMax), 2) / (std::max(dr2, dR2Min_) * DR2_LSB);
0642       }
0643     }
0644     unsigned int ieta = find_ieta(region, pfallne[in].hwEta);
0645     bool isEM = pfallne[in].hwId.isPhoton();
0646     std::pair<float, float> ptAndW = sum2puppiPt_flt(sum, pfallne[in].floatPt(), ieta, isEM, in);
0647     outallne_nocut[in].fill(
0648         region, pfallne[in], Scales::makePtFromFloat(ptAndW.first), l1ct::puppiWgt_t(ptAndW.second));
0649     if (region.isFiducial(pfallne[in]) && outallne_nocut[in].hwPt >= ptCut_[ieta]) {
0650       outallne[in] = outallne_nocut[in];
0651     }
0652   }
0653   puppisort_and_crop_ref(nOut_, outallne, outselne);
0654 }
0655 
0656 void l1ct::LinPuppiEmulator::run(const PFInputRegion &in,
0657                                  const std::vector<l1ct::PVObjEmu> &pvs,
0658                                  OutputRegion &out) const {
0659   if (debug_) {
0660     dbgPrintf("\nWill run LinPuppi in region eta %+5.2f, phi %+5.2f, pv0 int Z %+d\n",
0661               in.region.floatEtaCenter(),
0662               in.region.floatPhiCenter(),
0663               pvs.front().hwZ0.to_int());
0664   }
0665   if (std::abs(in.region.floatEtaCenter()) < 2.5) {  // within tracker
0666     std::vector<PuppiObjEmu> outallch, outallne_nocut, outallne, outselne;
0667     linpuppi_chs_ref(in.region, pvs, out.pfcharged, outallch);
0668     linpuppi_ref(in.region, in.track, pvs, out.pfneutral, outallne_nocut, outallne, outselne);
0669     // ensure proper sizes of the vectors, to get accurate sorting wrt firmware
0670     const std::vector<PuppiObjEmu> &ne = (nOut_ == nIn_ ? outallne : outselne);
0671     unsigned int nch = outallch.size(), nne = ne.size(), i;
0672     outallch.resize(nTrack_ + nOut_);
0673     for (i = nch; i < nTrack_; ++i)
0674       outallch[i].clear();
0675     for (unsigned int j = 0; j < nne; ++i, ++j)
0676       outallch[i] = ne[j];
0677     for (; i < nTrack_ + nOut_; ++i)
0678       outallch[i].clear();
0679     puppisort_and_crop_ref(nFinalSort_, outallch, out.puppi, finalSortAlgo_);
0680     // trim if needed
0681     while (!out.puppi.empty() && out.puppi.back().hwPt == 0)
0682       out.puppi.pop_back();
0683     out.puppi.shrink_to_fit();
0684   } else {  // forward
0685     std::vector<PuppiObjEmu> outallne_nocut, outallne;
0686     fwdlinpuppi_ref(in.region, in.hadcalo, outallne_nocut, outallne, out.puppi);
0687   }
0688 }