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