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 ®ion,
0253 const std::vector<PVObjEmu> &pv,
0254 const std::vector<PFChargedObjEmu> &pfch ,
0255 std::vector<PuppiObjEmu> &outallch ) 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_) {
0279 outallch[i].setHwDxy(dxy_t(pv[0].hwZ0));
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 ®ion, 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
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
0336
0337
0338
0339
0340
0341
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
0366
0367
0368
0369
0370
0371
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 ®ion,
0407 const std::vector<HadCaloObjEmu> &caloin ,
0408 std::vector<PuppiObjEmu> &outallne_nocut ,
0409 std::vector<PuppiObjEmu> &outallne ,
0410 std::vector<PuppiObjEmu> &outselne ) 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;
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_) {
0427 unsigned int dr2short = (dr2 >= dR2Min_ ? dr2 : dR2Min_) >> 5;
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 ®ion,
0446 const std::vector<TkObjEmu> &track ,
0447 const std::vector<PVObjEmu> &pv,
0448 const std::vector<PFNeutralObjEmu> &pfallne ,
0449 std::vector<PuppiObjEmu> &outallne_nocut ,
0450 std::vector<PuppiObjEmu> &outallne ,
0451 std::vector<PuppiObjEmu> &outselne ) 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;
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_) {
0486 unsigned int dr2short = (dr2 >= dR2Min_ ? dr2 : dR2Min_) >> 5;
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
0491
0492
0493
0494
0495
0496
0497
0498
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 {
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 ®ion,
0566 const std::vector<HadCaloObjEmu> &caloin ,
0567 std::vector<PuppiObjEmu> &outallne_nocut ,
0568 std::vector<PuppiObjEmu> &outallne ,
0569 std::vector<PuppiObjEmu> &outselne ) 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_) {
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 ®ion,
0602 const std::vector<TkObjEmu> &track ,
0603 const std::vector<PVObjEmu> &pv,
0604 const std::vector<PFNeutralObjEmu> &pfallne ,
0605 std::vector<PuppiObjEmu> &outallne_nocut ,
0606 std::vector<PuppiObjEmu> &outallne ,
0607 std::vector<PuppiObjEmu> &outselne ) 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);
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) {
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
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
0681 while (!out.puppi.empty() && out.puppi.back().hwPt == 0)
0682 out.puppi.pop_back();
0683 out.puppi.shrink_to_fit();
0684 } else {
0685 std::vector<PuppiObjEmu> outallne_nocut, outallne;
0686 fwdlinpuppi_ref(in.region, in.hadcalo, outallne_nocut, outallne, out.puppi);
0687 }
0688 }
