Back to home page

Project CMSSW displayed by LXR
 
 

    

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 &region,
0219                                               const std::vector<PVObjEmu> &pv,
0220                                               const std::vector<PFChargedObjEmu> &pfch /*[nTrack]*/,
0221                                               std::vector<PuppiObjEmu> &outallch /*[nTrack]*/) 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_) {                           // 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     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;                  // decimal bits the discriminator values
0284   const int alpha_bits = LINPUPPI_alpha_bits;            // decimal bits of the alpha values
0285   const int alphaSlope_bits = LINPUPPI_alphaSlope_bits;  // decimal bits of the alphaSlope values
0286   const int ptSlope_bits = LINPUPPI_ptSlope_bits;        // decimal bits of the ptSlope values
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;  // in pt scale
0292   const int ptZeroPh = ptZeroPh_[ieta] / LINPUPPI_ptLSB;  // in pt scale
0293   const int alphaCrop = alphaCrop_[ieta] * (1 << x2_bits);
0294   const int alphaSlopeNe =
0295       alphaSlope_[ieta] * std::log(2.) *
0296       (1 << alphaSlope_bits);  // we put a log(2) here since we compute alpha as log2(sum) instead of ln(sum)
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   // -- simplest version
0304   //int alpha = sum > 0 ? int(std::log2(float(sum) * LINPUPPI_pt2DR2_scale / (1<<sum_bitShift)) * (1 << alpha_bits) + 0.5) :  0;
0305   // -- re-written bringing terms out of the log
0306   //int alpha = sum > 0 ? int(std::log2(float(sum))*(1 << alpha_bits) + (std::log2(LINPUPPI_pt2DR2_scale) - sum_bitShift)*(1 << alpha_bits) + 0.5 ) :  0;
0307   // -- re-written for a LUT implementation of the log2
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));  // the maximum value of this term is log2lut_bits * (1 << alpha_bits) ~ 10*16 = 160 => fits in ap_uint<4+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   // -- re-written to fit in a single LUT
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     /*if (in <= 3) dbgPrintf("ref [%d]:  x2a(sum = %9lu): logarg = %9lu, sumterm = %9d, table[logarg] = %9d, ret pre-crop = %9d\n", 
0338           in, sum, logarg, 
0339           alphaSlope * int((std::log2(LINPUPPI_pt2DR2_scale) - sum_bitShift)*(1 << alpha_bits) + 0.5) - alphaSlope * alphaZero,
0340           alphaSlope * int(std::log2(float(logarg))*(1 << alpha_bits)), 
0341           x2a_lut); */
0342   } else {
0343     //if (in <= 3) dbgPrintf("ref [%d]:  x2a(sum = %9lu): logarg = %9lu, ret pre-crop = %9d\n",
0344     //        in, sum, logarg, x2a_lut);
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 &region,
0385                                              const std::vector<HadCaloObjEmu> &caloin /*[nIn]*/,
0386                                              std::vector<PuppiObjEmu> &outallne_nocut /*[nIn]*/,
0387                                              std::vector<PuppiObjEmu> &outallne /*[nIn]*/,
0388                                              std::vector<PuppiObjEmu> &outselne /*[nOut]*/) 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;  // 2 ^ sum_bitShift times (int pt^2)/(int dr2)
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);  // if dr is inside puppi cone
0407       if (dr2 <= dR2Max_) {
0408         ap_uint<9> dr2short = (dr2 >= dR2Min_ ? dr2 : dR2Min_) >> 5;  // reduce precision to make divide LUT cheaper
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         //      dr2short >= (dR2Min_ >> 5) = 2
0412         //      num <= (PTMAX2 >> 5) << sum_bitShift = (2^11) << 15 = 2^26
0413         //      ==> term <= 2^25
0414         //dbgPrintf("ref term [%2d,%2d]: dr = %8d  pt2_shift = %8lu  term = %12lu\n", in, it, dr2, std::min<uint64_t>(pt2 >> 5, PTMAX2 >> 5), term);
0415         assert(uint64_t(PTMAX2 << (sum_bitShift - 5)) / (dR2Min_ >> 5) <= (1 << 25));
0416         assert(term < (1 << 25));
0417         sum += term;
0418         //dbgPrintf("    pT cand %5.1f    pT item %5.1f    dR = %.3f   term = %.1f [dbl] = %lu [int]\n",
0419         //            caloin[in].floatPt(), caloin[it].floatPt(), std::sqrt(dr2*LINPUPPI_DR2LSB),
0420         //            double(std::min<uint64_t>(pt2 >> 5, 131071)<<15)/double(std::max<int>(dr2,dR2Min_) >> 5),
0421         //            term);
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 &region,
0436                                           const std::vector<TkObjEmu> &track /*[nTrack]*/,
0437                                           const std::vector<PVObjEmu> &pv, /*[nVtx]*/
0438                                           const std::vector<PFNeutralObjEmu> &pfallne /*[nIn]*/,
0439                                           std::vector<PuppiObjEmu> &outallne_nocut /*[nIn]*/,
0440                                           std::vector<PuppiObjEmu> &outallne /*[nIn]*/,
0441                                           std::vector<PuppiObjEmu> &outselne /*[nOut]*/) 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;  // 2 ^ sum_bitShift times (int pt^2)/(int dr2)
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);  // if dr is inside puppi cone
0472       if (dr2 <= dR2Max_) {
0473         ap_uint<9> dr2short = (dr2 >= dR2Min_ ? dr2 : dR2Min_) >> 5;  // reduce precision to make divide LUT cheaper
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         //      dr2short >= (dR2Min_ >> 5) = 2
0477         //      num <= (PTMAX2 >> 5) << sum_bitShift = (2^11) << 15 = 2^26
0478         //      ==> term <= 2^25
0479         //dbgPrintf("ref term [%2d,%2d]: dr = %8d  pt2_shift = %8lu  term = %12lu\n", in, it, dr2, std::min<uint64_t>(pt2 >> 5, PTMAX2 >> 5), term);
0480         assert(uint64_t(PTMAX2 << (sum_bitShift - 5)) / (dR2Min_ >> 5) <= (1 << 25));
0481         assert(term < (1 << 25));
0482         sum += term;
0483         //dbgPrintf("    pT cand %5.1f    pT item %5.1f    dR = %.3f   term = %.1f [dbl] = %lu [int]\n",
0484         //            pfallne[in].floatPt(), track[it].floatPt(), std::sqrt(dr2*LINPUPPI_DR2LSB),
0485         //            double(std::min<uint64_t>(pt2 >> 5, 131071)<<15)/double(std::max<int>(dr2,dR2Min_) >> 5),
0486         //            term);
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 {  // fakePuppi: keep the full candidate, but set the Puppi weight and some debug info into it
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 &region,
0552                                              const std::vector<HadCaloObjEmu> &caloin /*[nIn]*/,
0553                                              std::vector<PuppiObjEmu> &outallne_nocut /*[nIn]*/,
0554                                              std::vector<PuppiObjEmu> &outallne /*[nIn]*/,
0555                                              std::vector<PuppiObjEmu> &outselne /*[nOut]*/) 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);  // if dr is inside puppi cone
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 &region,
0590                                           const std::vector<TkObjEmu> &track /*[nTrack]*/,
0591                                           const std::vector<PVObjEmu> &pv,
0592                                           const std::vector<PFNeutralObjEmu> &pfallne /*[nIn]*/,
0593                                           std::vector<PuppiObjEmu> &outallne_nocut /*[nIn]*/,
0594                                           std::vector<PuppiObjEmu> &outallne /*[nIn]*/,
0595                                           std::vector<PuppiObjEmu> &outselne /*[nOut]*/) 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);  // if dr is inside puppi cone
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) {  // within tracker
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     // ensure proper sizes of the vectors, to get accurate sorting wrt firmware
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     // trim if needed
0665     while (!out.puppi.empty() && out.puppi.back().hwPt == 0)
0666       out.puppi.pop_back();
0667     out.puppi.shrink_to_fit();
0668   } else {  // forward
0669     std::vector<PuppiObjEmu> outallne_nocut, outallne;
0670     fwdlinpuppi_ref(in.region, in.hadcalo, outallne_nocut, outallne, out.puppi);
0671   }
0672 }
[ 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!