Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-10-01 01:02:14 

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