Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!