Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-06-20 01:53:35 

0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h"
0002 #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h"
0003 #ifdef CMSSW_GIT_HASH
0004 #include "DataFormats/L1TParticleFlow/interface/PFTrack.h"
0005 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloCrystalCluster.h"
0006 #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
0007 #endif
0008 #include <cmath>
0009 #include <cstdio>
0010 #include <algorithm>
0011 #include <memory>
0012 #include <iostream>
0013 #include <bitset>
0014 #include <vector>
0015 
0016 using namespace l1ct;
0017 
0018 #ifdef CMSSW_GIT_HASH
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0021 
0022 l1ct::PFTkEGAlgoEmuConfig::PFTkEGAlgoEmuConfig(const edm::ParameterSet &pset)
0023     : PFTkEGAlgoEmuConfig(pset.getParameter<uint32_t>("nTRACK"),
0024                           pset.getParameter<uint32_t>("nTRACK_EGIN"),
0025                           pset.getParameter<uint32_t>("nEMCALO_EGIN"),
0026                           pset.getParameter<uint32_t>("nEM_EGOUT"),
0027                           pset.getParameter<bool>("filterHwQuality"),
0028                           pset.getParameter<bool>("doBremRecovery"),
0029                           pset.getParameter<bool>("writeBeforeBremRecovery"),
0030                           pset.getParameter<int>("caloHwQual"),
0031                           pset.getParameter<bool>("doEndcapHwQual"),
0032                           pset.getParameter<double>("caloEtMin"),
0033                           pset.getParameter<double>("dEtaMaxBrem"),
0034                           pset.getParameter<double>("dPhiMaxBrem"),
0035                           pset.getParameter<std::vector<double>>("absEtaBoundaries"),
0036                           pset.getParameter<std::vector<double>>("dEtaValues"),
0037                           pset.getParameter<std::vector<double>>("dPhiValues"),
0038                           pset.getParameter<double>("trkQualityPtMin"),
0039                           pset.getParameter<uint32_t>("algorithm"),
0040                           pset.getParameter<uint32_t>("nCompCandPerCluster"),
0041                           pset.getParameter<bool>("writeEGSta"),
0042                           IsoParameters(pset.getParameter<edm::ParameterSet>("tkIsoParametersTkEle")),
0043                           IsoParameters(pset.getParameter<edm::ParameterSet>("tkIsoParametersTkEm")),
0044                           IsoParameters(pset.getParameter<edm::ParameterSet>("pfIsoParametersTkEle")),
0045                           IsoParameters(pset.getParameter<edm::ParameterSet>("pfIsoParametersTkEm")),
0046                           pset.getParameter<bool>("doTkIso"),
0047                           pset.getParameter<bool>("doPfIso"),
0048                           static_cast<EGIsoEleObjEmu::IsoType>(pset.getParameter<uint32_t>("hwIsoTypeTkEle")),
0049                           static_cast<EGIsoObjEmu::IsoType>(pset.getParameter<uint32_t>("hwIsoTypeTkEm")),
0050                           pset.getParameter<std::vector<edm::ParameterSet>>("compositeParametersTkEle")
0051                               .at(pset.getParameter<uint32_t>("algorithm")),
0052                           pset.getUntrackedParameter<uint32_t>("debug", 0)) {}
0053 
0054 edm::ParameterSetDescription l1ct::PFTkEGAlgoEmuConfig::getParameterSetDescription() {
0055   edm::ParameterSetDescription description;
0056   description.addUntracked<unsigned int>("debug", 0);
0057   description.add<unsigned int>("nTRACK");
0058   description.add<unsigned int>("nTRACK_EGIN");
0059   description.add<unsigned int>("nEMCALO_EGIN");
0060   description.add<unsigned int>("nEM_EGOUT");
0061   description.add<bool>("doBremRecovery", false);
0062   description.add<bool>("writeBeforeBremRecovery", false);
0063   description.add<bool>("filterHwQuality", false);
0064   description.add<int>("caloHwQual", 4);
0065   description.add<bool>("doEndcapHwQual", false);
0066   description.add<double>("dEtaMaxBrem", 0.02);
0067   description.add<double>("dPhiMaxBrem", 0.1);
0068   description.add<std::vector<double>>("absEtaBoundaries",
0069                                        {
0070                                            0.0,
0071                                            0.9,
0072                                            1.5,
0073                                        });
0074   description.add<std::vector<double>>("dEtaValues",
0075                                        {
0076                                            0.025,
0077                                            0.015,
0078                                            0.01,
0079                                        });
0080   description.add<std::vector<double>>("dPhiValues",
0081                                        {
0082                                            0.07,
0083                                            0.07,
0084                                            0.07,
0085                                        });
0086   description.add<double>("caloEtMin", 0.0);
0087   description.add<double>("trkQualityPtMin", 10.0);
0088   description.add<bool>("writeEGSta", false);
0089   description.add<edm::ParameterSetDescription>("tkIsoParametersTkEm", IsoParameters::getParameterSetDescription());
0090   description.add<edm::ParameterSetDescription>("tkIsoParametersTkEle", IsoParameters::getParameterSetDescription());
0091   description.add<edm::ParameterSetDescription>("pfIsoParametersTkEm", IsoParameters::getParameterSetDescription());
0092   description.add<edm::ParameterSetDescription>("pfIsoParametersTkEle", IsoParameters::getParameterSetDescription());
0093   description.add<bool>("doTkIso", true);
0094   description.add<bool>("doPfIso", true);
0095   description.add<unsigned int>("hwIsoTypeTkEle", 0);
0096   description.add<unsigned int>("hwIsoTypeTkEm", 2);
0097   description.add<unsigned int>("algorithm", 0);
0098   description.add<unsigned int>("nCompCandPerCluster", 3);
0099 
0100   description.addVPSet("compositeParametersTkEle", CompIDParameters::getParameterSetDescription());
0101 
0102   return description;
0103 }
0104 
0105 l1ct::PFTkEGAlgoEmuConfig::IsoParameters::IsoParameters(const edm::ParameterSet &pset)
0106     : IsoParameters(pset.getParameter<double>("tkQualityPtMin"),
0107                     pset.getParameter<double>("dZ"),
0108                     pset.getParameter<double>("dRMin"),
0109                     pset.getParameter<double>("dRMax")) {}
0110 
0111 edm::ParameterSetDescription l1ct::PFTkEGAlgoEmuConfig::IsoParameters::getParameterSetDescription() {
0112   edm::ParameterSetDescription description;
0113   description.add<double>("tkQualityPtMin");
0114   description.add<double>("dZ", 0.6);
0115   description.add<double>("dRMin");
0116   description.add<double>("dRMax");
0117   return description;
0118 }
0119 
0120 l1ct::PFTkEGAlgoEmuConfig::CompIDParameters::CompIDParameters(const edm::ParameterSet &pset)
0121     : CompIDParameters(pset.getParameter<edm::ParameterSet>("loose_wp").getParameter<std::vector<double>>("bins"),
0122                        pset.getParameter<edm::ParameterSet>("loose_wp").getParameter<std::vector<double>>("values"),
0123                        pset.getParameter<edm::ParameterSet>("tight_wp").getParameter<std::vector<double>>("bins"),
0124                        pset.getParameter<edm::ParameterSet>("tight_wp").getParameter<std::vector<double>>("values"),
0125                        pset.getParameter<std::string>("model"),
0126                        pset.getParameter<double>("dPhi_max"),
0127                        pset.getParameter<double>("dEta_max")) {}
0128 
0129 edm::ParameterSetDescription l1ct::PFTkEGAlgoEmuConfig::CompIDParameters::getParameterSetDescription() {
0130   edm::ParameterSetDescription wp_description;
0131   wp_description.addOptional<std::vector<double>>("bins");
0132   wp_description.addOptional<std::vector<double>>("values");
0133 
0134   edm::ParameterSetDescription description;
0135   description.addOptional<edm::ParameterSetDescription>("loose_wp", wp_description);
0136   description.addOptional<edm::ParameterSetDescription>("tight_wp", wp_description);
0137   description.addOptional<std::string>("model");
0138   description.add<double>("dPhi_max", 0.2);
0139   description.add<double>("dEta_max", 0.2);
0140   return description;
0141 }
0142 #endif
0143 
0144 //Constructor to be used with createWP factory methods
0145 l1ct::PFTkEGAlgoEmuConfig::CompIDParameters::CompIDParameters(const std::vector<double> &loose_wp_bins,
0146                                                               const std::vector<double> &loose_wp,
0147                                                               const std::vector<double> &tight_wp_bins,
0148                                                               const std::vector<double> &tight_wp,
0149                                                               const std::string &model,
0150                                                               double dphi_max,
0151                                                               double deta_max)
0152     : loose_wp_bins_(loose_wp_bins),
0153       loose_wp_(loose_wp),
0154       tight_wp_bins_(tight_wp_bins),
0155       tight_wp_(tight_wp),
0156       conifer_model_(model),
0157       dPhi_max_(dphi_max),
0158       dEta_max_(deta_max) {}
0159 
0160 l1ct::TkEGEleAssociationModel::TkEGEleAssociationModel(const l1ct::PFTkEGAlgoEmuConfig::CompIDParameters &params,
0161                                                        int debug)
0162     : loose_wp_(createWP(params.loose_wp_bins_, params.loose_wp_)),
0163       tight_wp_(createWP(params.tight_wp_bins_, params.tight_wp_)),
0164       dphi2_max_(params.dPhi_max_ * params.dPhi_max_),
0165       deta2_max_(params.dEta_max_ * params.dEta_max_),
0166       debug_(debug) {}
0167 
0168 bool l1ct::TkEGEleAssociationModel::geometric_match(const EmCaloObjEmu &calo, const TkObjEmu &tk) const {
0169   float d_phi = PFTkEGAlgoEmulator::deltaPhi(tk.floatPhi(), calo.floatPhi());
0170   float d_eta = tk.floatEta() - calo.floatEta();  // We only use it squared
0171   return ((d_phi * d_phi / dphi2_max_) + (d_eta * d_eta / deta2_max_) < 1.);
0172 }
0173 
0174 l1ct::TkEgCID_EE_v0::TkEgCID_EE_v0(const l1ct::PFTkEGAlgoEmuConfig::CompIDParameters &params, int debug)
0175     : TkEGEleAssociationModel(params, debug) {
0176 #ifdef CMSSW_GIT_HASH
0177   auto resolvedFileName = edm::FileInPath(params.conifer_model_).fullPath();
0178 #else
0179   auto resolvedFileName = params.conifer_model_;
0180 #endif
0181   model_ = std::make_unique<conifer::BDT<bdt_feature_t, bdt_score_t, false>>(resolvedFileName);
0182 }
0183 
0184 id_score_t l1ct::TkEgCID_EE_v0::compute_score(const CompositeCandidate &cand,
0185                                               const std::vector<EmCaloObjEmu> &emcalo,
0186                                               const std::vector<TkObjEmu> &track,
0187                                               const std::vector<float> additional_vars) const {
0188   // Get the cluster/track objects that form the composite candidate
0189   const auto &calo = emcalo[cand.cluster_idx];
0190   const auto &tk = track[cand.track_idx];
0191 
0192   // Prepare the input features
0193   bdt_feature_t hoe = calo.hwHoe;
0194   bdt_feature_t tkpt = tk.hwPt;
0195   bdt_feature_t srrtot = calo.hwSrrTot;
0196   bdt_feature_t deta = tk.hwEta - calo.hwEta;
0197   ap_ufixed<18, 0> calo_invPt = l1ct::invert_with_shift<pt_t, ap_ufixed<18, 0>, 1024>(calo.hwPt);
0198   bdt_feature_t dpt = tk.hwPt * calo_invPt;
0199   bdt_feature_t meanz = calo.hwMeanZ;
0200   bdt_feature_t dphi = tk.hwPhi - calo.hwPhi;
0201   bdt_feature_t nstubs = tk.hwStubs;
0202   bdt_feature_t chi2rphi = tk.hwRedChi2RPhi;
0203   bdt_feature_t chi2rz = tk.hwRedChi2RZ;
0204   bdt_feature_t chi2bend = tk.hwRedChi2Bend;
0205 
0206   // Run BDT inference
0207   std::vector<bdt_feature_t> inputs = {tkpt, hoe, srrtot, deta, dphi, dpt, meanz, nstubs, chi2rphi, chi2rz, chi2bend};
0208   std::vector<bdt_score_t> bdt_score = model_->decision_function(inputs);
0209   constexpr unsigned int MAX_SCORE = (1 << (bdt_score_t::iwidth - 1));
0210   return bdt_score[0] / MAX_SCORE;
0211 }
0212 
0213 l1ct::TkEgCID_EE_v1::TkEgCID_EE_v1(const l1ct::PFTkEGAlgoEmuConfig::CompIDParameters &params, int debug)
0214     : TkEGEleAssociationModel(params, debug) {
0215 #ifdef CMSSW_GIT_HASH
0216   auto resolvedFileName = edm::FileInPath(params.conifer_model_).fullPath();
0217 #else
0218   auto resolvedFileName = params.conifer_model_;
0219 #endif
0220   model_ = std::make_unique<conifer::BDT<bdt_feature_t, bdt_score_t, false>>(resolvedFileName);
0221 }
0222 
0223 id_score_t l1ct::TkEgCID_EE_v1::compute_score(const CompositeCandidate &cand,
0224                                               const std::vector<EmCaloObjEmu> &emcalo,
0225                                               const std::vector<TkObjEmu> &track,
0226                                               const std::vector<float> additional_vars) const {
0227   float sumTkPt = additional_vars[1];
0228 #ifdef CMSSW_GIT_HASH
0229   // NOTE: this is not yet ready for emulation!
0230   // Get the cluster/track objects that form the composite candidate
0231   const auto &calo = emcalo[cand.cluster_idx];
0232   const auto &tk = track[cand.track_idx];
0233   const l1t::PFTrack *pftk = tk.src;
0234   const l1t::HGCalMulticluster *cl3d = dynamic_cast<const l1t::HGCalMulticluster *>(calo.src);
0235 
0236   // Prepare the input features
0237   bdt_feature_t cl_coreshowerlength = cl3d->coreShowerLength();
0238   bdt_feature_t cl_meanz = std::fabs(cl3d->zBarycenter());
0239   bdt_feature_t cl_spptot = cl3d->sigmaPhiPhiTot();
0240   bdt_feature_t cl_seetot = cl3d->sigmaEtaEtaTot();
0241   bdt_feature_t cl_szz = cl3d->sigmaZZ();
0242   bdt_feature_t cl_multiClassPionIdScore = calo.floatPiProb();
0243   bdt_feature_t cl_multiClassEmIdScore = calo.floatEmProb();
0244   bdt_feature_t tk_ptFrac = pftk->pt() / sumTkPt;
0245   bdt_feature_t cltk_ptRatio = calo.floatPt() / pftk->pt();
0246   bdt_feature_t cltk_absDeta = fabs(cl3d->eta() - pftk->caloEta());
0247   bdt_feature_t cltk_absDphi = fabs(cl3d->phi() - pftk->caloPhi());
0248 
0249   // Run BDT inference
0250   std::vector<bdt_feature_t> inputs = {cl_coreshowerlength,
0251                                        cl_meanz,
0252                                        cl_spptot,
0253                                        cl_seetot,
0254                                        cl_szz,
0255                                        cl_multiClassPionIdScore,
0256                                        cl_multiClassEmIdScore,
0257                                        tk_ptFrac,
0258                                        cltk_ptRatio,
0259                                        cltk_absDeta,
0260                                        cltk_absDphi};
0261   std::vector<bdt_score_t> bdt_score = model_->decision_function(inputs);
0262   // std::cout << "  out BDT score: " << bdt_score[0] << std::endl;
0263   constexpr unsigned int MAX_SCORE = 1 << (bdt_score_t::iwidth - 1);
0264   return bdt_score[0] / MAX_SCORE;
0265 #else
0266   return 0;
0267 #endif
0268 }
0269 
0270 l1ct::TkEgCID_EB_v0::TkEgCID_EB_v0(const l1ct::PFTkEGAlgoEmuConfig::CompIDParameters &params, int debug)
0271     : TkEGEleAssociationModel(params, debug) {
0272 #ifdef CMSSW_GIT_HASH
0273   auto resolvedFileName = edm::FileInPath(params.conifer_model_).fullPath();
0274 #else
0275   auto resolvedFileName = params.conifer_model_;
0276 #endif
0277   model_ = std::make_unique<conifer::BDT<bdt_feature_t, bdt_score_t, false>>(resolvedFileName);
0278 }
0279 
0280 id_score_t l1ct::TkEgCID_EB_v0::compute_score(const CompositeCandidate &cand,
0281                                               const std::vector<EmCaloObjEmu> &emcalo,
0282                                               const std::vector<TkObjEmu> &track,
0283                                               const std::vector<float> additional_vars) const {
0284   unsigned int nTkMatch = (unsigned int)(additional_vars[0]);
0285   float sumTkPt = additional_vars[1];
0286 
0287 #ifdef CMSSW_GIT_HASH
0288   // NOTE: not yet ready for HLS testbench
0289   // Get the cluster/track objects that form the composite candidate
0290   const auto &calo = emcalo[cand.cluster_idx];
0291   const auto &tk = track[cand.track_idx];
0292   const l1tp2::CaloCrystalCluster *crycl = dynamic_cast<const l1tp2::CaloCrystalCluster *>(calo.src);
0293 
0294   // Prepare the input features
0295   // NOTE: 16 bit estimate for the inversion is approximate
0296   ap_ufixed<16, 0> calo_invPt = l1ct::invert_with_shift<pt_t, ap_ufixed<16, 0>, 1024>(calo.hwPt);
0297   // NOTE: this could be computed once per cluster and passed directly to the function
0298   ap_ufixed<16, 0> sumTk_invPt = l1ct::invert_with_shift<pt_t, ap_ufixed<16, 0>, 1024>(pt_t(sumTkPt));
0299   ap_ufixed<16, 0> tk_invPt = l1ct::invert_with_shift<pt_t, ap_ufixed<16, 0>, 1024>(tk.hwPt);
0300 
0301   constexpr std::array<float, 1 << l1ct::redChi2Bin_t::width> chi2RPhiBins = {
0302       {0.0, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 10.0, 15.0, 20.0, 35.0, 60.0, 200.0}};
0303 
0304   bdt_feature_t cl_pt = calo.floatPt();
0305   bdt_feature_t cl_ss = crycl->e2x5() / crycl->e5x5();
0306   bdt_feature_t cl_relIso = iso_t(crycl->isolation()) * calo_invPt;
0307   bdt_feature_t cl_staWP = calo.hwEmID & 0x1;
0308   bdt_feature_t cl_looseTkWP = calo.hwEmID & 0x2;
0309   bdt_feature_t tk_chi2RPhi = chi2RPhiBins[tk.hwRedChi2RPhi.to_int()];
0310   bdt_feature_t tk_ptFrac = tk.hwPt * sumTk_invPt;
0311   bdt_feature_t cltk_ptRatio = calo.hwPt * tk_invPt;
0312   bdt_feature_t cltk_nTkMatch = nTkMatch;
0313   bdt_feature_t cltk_absDeta = fabs(tk.floatEta() - calo.floatEta());
0314   bdt_feature_t cltk_absDphi = fabs(tk.floatPhi() - calo.floatPhi());
0315 
0316   // Run BDT inference
0317   std::vector<bdt_feature_t> inputs = {cl_pt,
0318                                        cl_ss,
0319                                        cl_relIso,
0320                                        cl_staWP,
0321                                        cl_looseTkWP,
0322                                        tk_chi2RPhi,
0323                                        tk_ptFrac,
0324                                        cltk_ptRatio,
0325                                        cltk_nTkMatch,
0326                                        cltk_absDeta,
0327                                        cltk_absDphi};
0328   std::vector<bdt_score_t> bdt_score = model_->decision_function(inputs);
0329   // std::cout << "  out BDT score: " << bdt_score[0] << std::endl;
0330   constexpr unsigned int MAX_SCORE = 1 << (bdt_score_t::iwidth - 1);
0331   return bdt_score[0] / MAX_SCORE;  // normalize to [-1,1]
0332 #else
0333   return 0;
0334 #endif
0335 }
0336 
0337 l1ct::TkEgCID_EB_v1::TkEgCID_EB_v1(const l1ct::PFTkEGAlgoEmuConfig::CompIDParameters &params, int debug)
0338     : TkEGEleAssociationModel(params, debug) {
0339 #ifdef CMSSW_GIT_HASH
0340   auto resolvedFileName = edm::FileInPath(params.conifer_model_).fullPath();
0341 #else
0342   auto resolvedFileName = params.conifer_model_;
0343 #endif
0344   model_ = std::make_unique<conifer::BDT<bdt_feature_t, bdt_score_t, false>>(resolvedFileName);
0345 }
0346 
0347 id_score_t l1ct::TkEgCID_EB_v1::compute_score(const CompositeCandidate &cand,
0348                                               const std::vector<EmCaloObjEmu> &emcalo,
0349                                               const std::vector<TkObjEmu> &track,
0350                                               const std::vector<float> additional_vars) const {
0351   unsigned int nTkMatch = (unsigned int)(additional_vars[0]);
0352   float sumTkPt = additional_vars[1];
0353 
0354   // NOTE: not yet ready for HLS testbench
0355   // Get the cluster/track objects that form the composite candidate
0356   const auto &calo = emcalo[cand.cluster_idx];
0357   const auto &tk = track[cand.track_idx];
0358 
0359   // Prepare the input features
0360   // NOTE: this could be computed once per cluster and passed directly to the function
0361   ap_ufixed<16, 0> tk_invPt = l1ct::invert_with_shift<pt_t, ap_ufixed<16, 0>, 1024>(tk.hwPt);
0362 
0363   constexpr std::array<float, 1 << l1ct::redChi2Bin_t::width> chi2RPhiBins = {
0364       {0.0, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 10.0, 15.0, 20.0, 35.0, 60.0, 200.0}};
0365 
0366   float cl_pt = calo.floatPt();
0367   //This two ratios will be computed in the calotrigger and passed to the CTL1 in 6 bits
0368   float cl_ss = emcalo[cand.cluster_idx].hwShowerShape.to_float();
0369   float cl_relIso = emcalo[cand.cluster_idx].hwRelIso.to_float();
0370   float cl_staWP = calo.hwEmID & 0x1;
0371   float cl_looseTkWP = (calo.hwEmID & 0x2) == 0x2;
0372   float tk_chi2RPhi = chi2RPhiBins[tk.hwRedChi2RPhi.to_int()];
0373   float tk_ptFrac = sumTkPt * tk_invPt.to_float();
0374   float cltk_ptRatio = calo.hwPt * tk_invPt;
0375   float cltk_nTkMatch = nTkMatch;
0376   float cltk_absDeta = fabs(tk.hwEta.to_int() - calo.hwEta.to_int());
0377   float cltk_absDphi = fabs(tk.hwPhi.to_int() - calo.hwPhi.to_int());
0378 
0379   // Scaling
0380   bdt_feature_t scaled_cl_pt = scale(cl_pt, 1.5, 5);
0381   bdt_feature_t scaled_cl_ss = scale(cl_ss, 0.1875, -1);
0382   bdt_feature_t scaled_cl_relIso = scale(cl_relIso, 0.0, -1);
0383   bdt_feature_t scaled_cl_staWP = scale(cl_staWP, 0.0, 0);
0384   bdt_feature_t scaled_cl_looseTkWP = scale(cl_looseTkWP, 0.0, 0);
0385   bdt_feature_t scaled_tk_chi2RPhi = scale(tk_chi2RPhi, 0.0, 3);
0386   bdt_feature_t scaled_tk_ptFrac = scale(tk_ptFrac, 1.0, 5);
0387   bdt_feature_t scaled_cltk_ptRatio = scale(cltk_ptRatio, 0.0003669276, 4);
0388   bdt_feature_t scaled_cltk_nTkMatch = scale(cltk_nTkMatch, 1.0, 3);
0389   bdt_feature_t scaled_cltk_absDeta = scale(cltk_absDeta, 0.0, 2);
0390   bdt_feature_t scaled_cltk_absDphi = scale(cltk_absDphi, 0.0, 5);
0391 
0392   // Run BDT inference
0393   std::vector<bdt_feature_t> inputs = {scaled_cl_pt,
0394                                        scaled_cl_ss,
0395                                        scaled_cl_relIso,
0396                                        scaled_cl_staWP,
0397                                        scaled_cl_looseTkWP,
0398                                        scaled_tk_chi2RPhi,
0399                                        scaled_tk_ptFrac,
0400                                        scaled_cltk_ptRatio,
0401                                        scaled_cltk_nTkMatch,
0402                                        scaled_cltk_absDeta,
0403                                        scaled_cltk_absDphi};
0404   std::vector<bdt_score_t> bdt_score = model_->decision_function(inputs);
0405   if (debug_ > 3) {
0406     dbgCout() << "[REF] EM calo pt: " << calo.hwPt << " tk pt " << tk.hwPt << std::endl;
0407     if (debug_ > 5) {
0408       dbgCout() << " .  [0] cl_pt: " << cl_pt << std::endl;
0409       dbgCout() << " .  [0] scaled cl_pt: " << scaled_cl_pt << std::endl;
0410       dbgCout() << " .  [1] scaled cl_ss: " << scaled_cl_ss << std::endl;
0411       dbgCout() << " .  [2] scaled cl_relIso: " << scaled_cl_relIso << std::endl;
0412       dbgCout() << " .  [3] scaled cl_staWP: " << scaled_cl_staWP << std::endl;
0413       dbgCout() << " .  [4] scaled cl_looseTkWP: " << scaled_cl_looseTkWP << std::endl;
0414       dbgCout() << " .  [5] scaled tk_chi2RPhi: " << scaled_tk_chi2RPhi << std::endl;
0415       dbgCout() << " .  [6] scaled tk_ptFrac: " << scaled_tk_ptFrac << std::endl;
0416       dbgCout() << " .  [7] scaled cltk_ptRatio: " << scaled_cltk_ptRatio << std::endl;
0417       dbgCout() << " .  [8] scaled cltk_nTkMatch: " << scaled_cltk_nTkMatch << std::endl;
0418       dbgCout() << " .  [9] scaled cltk_absDeta: " << scaled_cltk_absDeta << std::endl;
0419       dbgCout() << " .  [10] scaled cltk_absDphi: " << scaled_cltk_absDphi << std::endl;
0420     }
0421     dbgCout() << "  out BDT score: " << bdt_score[0] / 8 << std::endl;
0422   }
0423 
0424   return bdt_score[0] / 8;  // normalize to [-1,1]
0425 }
0426 
0427 PFTkEGAlgoEmulator::PFTkEGAlgoEmulator(const PFTkEGAlgoEmuConfig &config)
0428     : cfg(config), tkEleModel_(nullptr), debug_(cfg.debug) {
0429   if (cfg.algorithm == PFTkEGAlgoEmuConfig::Algo::compositeEE_v0) {
0430     tkEleModel_ = std::make_unique<TkEgCID_EE_v0>(cfg.compIDparams, cfg.debug);
0431   } else if (cfg.algorithm == PFTkEGAlgoEmuConfig::Algo::compositeEB_v0) {
0432     tkEleModel_ = std::make_unique<TkEgCID_EB_v0>(cfg.compIDparams, cfg.debug);
0433   } else if (cfg.algorithm == PFTkEGAlgoEmuConfig::Algo::compositeEE_v1) {
0434     tkEleModel_ = std::make_unique<TkEgCID_EE_v1>(cfg.compIDparams, cfg.debug);
0435   } else if (cfg.algorithm == PFTkEGAlgoEmuConfig::Algo::compositeEB_v1) {
0436     tkEleModel_ = std::make_unique<TkEgCID_EB_v1>(cfg.compIDparams, cfg.debug);
0437   }
0438 }
0439 
0440 void PFTkEGAlgoEmulator::toFirmware(const PFInputRegion &in,
0441                                     PFRegion &region,
0442                                     EmCaloObj emcalo[/*nCALO*/],
0443                                     TkObj track[/*nTRACK*/]) const {
0444   region = in.region;
0445   l1ct::toFirmware(in.track, cfg.nTRACK_EGIN, track);
0446   l1ct::toFirmware(in.emcalo, cfg.nEMCALO_EGIN, emcalo);
0447   if (debug_ > 0)
0448     dbgCout() << "# of inpput tracks: " << in.track.size() << " (max: " << cfg.nTRACK_EGIN << ")"
0449               << " emcalo: " << in.emcalo.size() << "(" << cfg.nEMCALO_EGIN << ")" << std::endl;
0450 }
0451 
0452 void PFTkEGAlgoEmulator::toFirmware(const OutputRegion &out, EGIsoObj out_egphs[], EGIsoEleObj out_egeles[]) const {
0453   l1ct::toFirmware(out.egphoton, cfg.nEM_EGOUT, out_egphs);
0454   l1ct::toFirmware(out.egelectron, cfg.nEM_EGOUT, out_egeles);
0455   if (debug_ > 0)
0456     dbgCout() << "# output photons: " << out.egphoton.size() << " electrons: " << out.egelectron.size() << std::endl;
0457 }
0458 
0459 void PFTkEGAlgoEmulator::toFirmware(
0460     const PFInputRegion &in, const l1ct::PVObjEmu &pvin, PFRegion &region, TkObj track[/*nTRACK*/], PVObj &pv) const {
0461   region = in.region;
0462   l1ct::toFirmware(in.track, cfg.nTRACK, track);
0463   pv = pvin;
0464   if (debug_ > 0)
0465     dbgCout() << "# of inpput tracks: " << in.track.size() << " (max: " << cfg.nTRACK << ")" << std::endl;
0466 }
0467 
0468 float PFTkEGAlgoEmulator::deltaPhi(float phi1, float phi2) {
0469   // reduce to [-pi,pi]
0470   float x = phi1 - phi2;
0471   float o2pi = 1. / (2. * M_PI);
0472   if (std::abs(x) <= float(M_PI))
0473     return x;
0474   float n = std::round(x * o2pi);
0475   return x - n * float(2. * M_PI);
0476 }
0477 
0478 void PFTkEGAlgoEmulator::link_emCalo2emCalo(const std::vector<EmCaloObjEmu> &emcalo,
0479                                             std::vector<int> &emCalo2emCalo) const {
0480   // NOTE: we assume the input to be sorted!!!
0481   for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0482     auto &calo = emcalo[ic];
0483     if (emCalo2emCalo[ic] != -1)
0484       continue;
0485 
0486     for (int jc = ic + 1; jc < nc; ++jc) {
0487       if (emCalo2emCalo[jc] != -1)
0488         continue;
0489 
0490       auto &otherCalo = emcalo[jc];
0491 
0492       if (fabs(otherCalo.floatEta() - calo.floatEta()) < cfg.dEtaMaxBrem &&
0493           fabs(deltaPhi(otherCalo.floatPhi(), calo.floatPhi())) < cfg.dPhiMaxBrem) {
0494         emCalo2emCalo[jc] = ic;
0495       }
0496     }
0497   }
0498 }
0499 
0500 void PFTkEGAlgoEmulator::link_emCalo2tk_elliptic(const PFRegionEmu &r,
0501                                                  const std::vector<EmCaloObjEmu> &emcalo,
0502                                                  const std::vector<TkObjEmu> &track,
0503                                                  std::vector<int> &emCalo2tk) const {
0504   unsigned int nTrackMax = std::min<unsigned>(track.size(), cfg.nTRACK_EGIN);
0505   for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0506     auto &calo = emcalo[ic];
0507 
0508     float dPtMin = 999;
0509     for (unsigned int itk = 0; itk < nTrackMax; ++itk) {
0510       const auto &tk = track[itk];
0511       if (tk.floatPt() < cfg.trkQualityPtMin)
0512         continue;
0513 
0514       float d_phi = deltaPhi(tk.floatPhi(), calo.floatPhi());
0515       float d_eta = tk.floatEta() - calo.floatEta();  // We only use it squared
0516 
0517       auto eta_index =
0518           std::distance(cfg.absEtaBoundaries.begin(),
0519                         std::lower_bound(
0520                             cfg.absEtaBoundaries.begin(), cfg.absEtaBoundaries.end(), abs(r.floatGlbEta(calo.hwEta)))) -
0521           1;
0522 
0523       float dEtaMax = cfg.dEtaValues[eta_index];
0524       float dPhiMax = cfg.dPhiValues[eta_index];
0525 
0526       if (debug_ > 2 && calo.hwPt > 0) {
0527         dbgCout() << "[REF] tried to link calo " << ic << " (pt " << calo.intPt() << ", eta " << calo.intEta()
0528                   << ", phi " << calo.intPhi() << ") "
0529                   << " to tk " << itk << " (pt " << tk.intPt() << ", eta " << tk.intEta() << ", phi " << tk.intPhi()
0530                   << "): "
0531                   << " eta_index " << eta_index << ", "
0532                   << " dEta " << d_eta << " (max " << dEtaMax << "), dPhi " << d_phi << " (max " << dPhiMax << ") "
0533                   << " ellipse = "
0534                   << (((d_phi / dPhiMax) * (d_phi / dPhiMax)) + ((d_eta / dEtaMax) * (d_eta / dEtaMax))) << "\n";
0535       }
0536       if ((((d_phi / dPhiMax) * (d_phi / dPhiMax)) + ((d_eta / dEtaMax) * (d_eta / dEtaMax))) < 1.) {
0537         // NOTE: for now we implement only best pt match. This is NOT what is done in the L1TkElectronTrackProducer
0538         if (fabs(tk.floatPt() - calo.floatPt()) < dPtMin) {
0539           emCalo2tk[ic] = itk;
0540           dPtMin = fabs(tk.floatPt() - calo.floatPt());
0541         }
0542       }
0543     }
0544   }
0545 }
0546 
0547 void PFTkEGAlgoEmulator::link_emCalo2tk_composite_eb_ee(const PFRegionEmu &r,
0548                                                         const std::vector<EmCaloObjEmu> &emcalo,
0549                                                         const std::vector<TkObjEmu> &track,
0550                                                         std::vector<int> &emCalo2tk,
0551                                                         std::vector<id_score_t> &emCaloTkBdtScore) const {
0552   unsigned int nTrackMax = std::min<unsigned>(track.size(), cfg.nTRACK_EGIN);
0553   for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0554     auto &calo = emcalo[ic];
0555 
0556     std::vector<CompositeCandidate> candidates;
0557     float sumTkPt = 0.;
0558     unsigned int nTkMatch = 0;
0559     for (unsigned int itk = 0; itk < nTrackMax; ++itk) {
0560       const auto &tk = track[itk];
0561       if (tk.floatPt() <= cfg.trkQualityPtMin)
0562         continue;
0563       bool keep = tkEleModel_->geometric_match(calo, tk);
0564       if (keep) {
0565         // Only store indices, dR and dpT for now. The other quantities are computed only for the best nCandPerCluster.
0566         CompositeCandidate cand;
0567         cand.cluster_idx = ic;
0568         cand.track_idx = itk;
0569         cand.dpt = std::abs(tk.floatPt() - calo.floatPt());
0570         candidates.push_back(cand);
0571         sumTkPt += tk.floatPt();
0572         nTkMatch++;
0573       }
0574       if (debug_ > 3)
0575         dbgCout() << "[REF] tried to link calo " << ic << " (pt " << calo.floatPt() << ", eta " << calo.intEta()
0576                   << ", phi " << calo.intPhi() << ") "
0577                   << " to tk " << itk << " (pt " << tk.floatPt() << ", eta " << tk.intEta() << ", phi " << tk.intPhi()
0578                   << ") keep: " << keep << std::endl;
0579     }
0580     // we use dpt as sort criteria
0581     std::sort(candidates.begin(),
0582               candidates.end(),
0583               [](const CompositeCandidate &a, const CompositeCandidate &b) -> bool { return a.dpt < b.dpt; });
0584     unsigned int nCandPerCluster = std::min<unsigned int>(candidates.size(), cfg.nCompCandPerCluster);
0585     if (nCandPerCluster == 0)
0586       continue;
0587 
0588     id_score_t maxScore = -(1 << (l1ct::id_score_t::iwidth - 1));
0589     int ibest = -1;
0590     for (unsigned int icand = 0; icand < nCandPerCluster; icand++) {
0591       auto &cand = candidates[icand];
0592       const std::vector<EmCaloObjEmu> &emcalo_sel = emcalo;
0593       id_score_t score = tkEleModel_->compute_score(cand, emcalo_sel, track, {float(nTkMatch), sumTkPt});
0594       if ((tkEleModel_->apply_wp_loose(score, emcalo_sel[cand.cluster_idx].floatPt())) && (score > maxScore)) {
0595         maxScore = score;
0596         ibest = icand;
0597       }
0598     }
0599     if (ibest != -1) {
0600       emCalo2tk[ic] = candidates[ibest].track_idx;
0601       emCaloTkBdtScore[ic] = maxScore;
0602     }
0603   }
0604 }
0605 
0606 void PFTkEGAlgoEmulator::sel_emCalo(unsigned int nmax_sel,
0607                                     const std::vector<EmCaloObjEmu> &emcalo,
0608                                     std::vector<EmCaloObjEmu> &emcalo_sel) const {
0609   for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0610     const auto &calo = emcalo[ic];
0611     if ((calo.hwPt == 0) || (cfg.filterHwQuality && calo.hwEmID != cfg.caloHwQual) ||
0612         (calo.floatPt() < cfg.emClusterPtMin))
0613       continue;
0614     emcalo_sel.push_back(calo);
0615     if (emcalo_sel.size() >= nmax_sel)
0616       break;
0617   }
0618 }
0619 
0620 void PFTkEGAlgoEmulator::run(const PFInputRegion &in, OutputRegion &out) const {
0621   if (debug_ > 1) {
0622     for (int ic = 0, nc = in.emcalo.size(); ic < nc; ++ic) {
0623       const auto &calo = in.emcalo[ic];
0624       if (calo.hwPt > 0)
0625         dbgCout() << "[REF] IN calo[" << ic << "] pt: " << calo.hwPt << " eta: " << calo.hwEta
0626                   << " (glb eta: " << in.region.floatGlbEta(calo.hwEta) << ") phi: " << calo.hwPhi
0627                   << "(glb phi: " << in.region.floatGlbPhi(calo.hwPhi) << ") qual: " << std::bitset<4>(calo.hwEmID)
0628                   << std::endl;
0629     }
0630   }
0631   // FIXME: can be removed in the endcap since now running with the "interceptor".
0632   // Might still be needed in barrel
0633   // filter and select first N elements of input clusters
0634   std::vector<EmCaloObjEmu> emcalo_sel;
0635   sel_emCalo(cfg.nEMCALO_EGIN, in.emcalo, emcalo_sel);
0636 
0637   std::vector<int> emCalo2emCalo(emcalo_sel.size(), -1);
0638   if (cfg.doBremRecovery)
0639     link_emCalo2emCalo(emcalo_sel, emCalo2emCalo);
0640 
0641   std::vector<int> emCalo2tk(emcalo_sel.size(), -1);
0642   std::vector<id_score_t> emCaloTkBdtScore(emcalo_sel.size(), 0);
0643 
0644   if (cfg.algorithm == PFTkEGAlgoEmuConfig::Algo::elliptic) {
0645     link_emCalo2tk_elliptic(in.region, emcalo_sel, in.track, emCalo2tk);
0646   } else {
0647     link_emCalo2tk_composite_eb_ee(in.region, emcalo_sel, in.track, emCalo2tk, emCaloTkBdtScore);
0648   }
0649 
0650   out.egsta.clear();
0651   std::vector<EGIsoObjEmu> egobjs;
0652   std::vector<EGIsoEleObjEmu> egeleobjs;
0653   eg_algo(in.region, emcalo_sel, in.track, emCalo2emCalo, emCalo2tk, emCaloTkBdtScore, out.egsta, egobjs, egeleobjs);
0654 
0655   unsigned int nEGOut = std::min<unsigned>(cfg.nEM_EGOUT, egobjs.size());
0656   unsigned int nEGEleOut = std::min<unsigned>(cfg.nEM_EGOUT, egeleobjs.size());
0657 
0658   // init output containers
0659   out.egphoton.clear();
0660   out.egelectron.clear();
0661   ptsort_ref(egobjs.size(), nEGOut, egobjs, out.egphoton);
0662   ptsort_ref(egeleobjs.size(), nEGEleOut, egeleobjs, out.egelectron);
0663 }
0664 
0665 void PFTkEGAlgoEmulator::eg_algo(const PFRegionEmu &region,
0666                                  const std::vector<EmCaloObjEmu> &emcalo,
0667                                  const std::vector<TkObjEmu> &track,
0668                                  const std::vector<int> &emCalo2emCalo,
0669                                  const std::vector<int> &emCalo2tk,
0670                                  const std::vector<id_score_t> &emCaloTkBdtScore,
0671                                  std::vector<EGObjEmu> &egstas,
0672                                  std::vector<EGIsoObjEmu> &egobjs,
0673                                  std::vector<EGIsoEleObjEmu> &egeleobjs) const {
0674   for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0675     auto &calo = emcalo[ic];
0676 
0677     // discard immediately EG objects that would not fall in the fiducial eta-phi region
0678     if (!region.isFiducial(calo))
0679       continue;
0680 
0681     if (debug_ > 3)
0682       dbgCout() << "[REF] SEL emcalo with pt: " << calo.hwPt << " qual: " << calo.hwEmID << " eta: " << calo.hwEta
0683                 << " phi " << calo.hwPhi << std::endl;
0684 
0685     int itk = emCalo2tk[ic];
0686     const id_score_t &bdt = emCaloTkBdtScore[ic];
0687 
0688     // check if brem recovery is on
0689     if (!cfg.doBremRecovery || cfg.writeBeforeBremRecovery) {
0690       // 1. create EG objects before brem recovery
0691       unsigned int egQual = calo.hwEmID;
0692       // If we write both objects with and without brem-recovery
0693       // bit 3 is used for the brem-recovery bit: if set = no recovery
0694       // (for consistency with the barrel hwQual where by default the brem recovery is done upstream)
0695       if (cfg.writeBeforeBremRecovery && cfg.doBremRecovery) {
0696         egQual = calo.hwEmID | 0x8;
0697       }
0698 
0699       addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, egQual, calo.hwPt, itk, bdt);
0700     }
0701 
0702     if (!cfg.doBremRecovery)
0703       continue;
0704 
0705     // check if the cluster has already been used in a brem reclustering
0706     if (emCalo2emCalo[ic] != -1)
0707       continue;
0708 
0709     pt_t ptBremReco = calo.hwPt;
0710     std::vector<unsigned int> components;
0711 
0712     for (int jc = ic; jc < nc; ++jc) {
0713       if (emCalo2emCalo[jc] == ic) {
0714         auto &otherCalo = emcalo[jc];
0715         ptBremReco += otherCalo.hwPt;
0716         components.push_back(jc);
0717       }
0718     }
0719 
0720     // 2. create EG objects with brem recovery
0721     addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, calo.hwEmID, ptBremReco, itk, bdt, components);
0722   }
0723 }
0724 
0725 EGObjEmu &PFTkEGAlgoEmulator::addEGStaToPF(std::vector<EGObjEmu> &egobjs,
0726                                            const EmCaloObjEmu &calo,
0727                                            const unsigned int hwQual,
0728                                            const pt_t ptCorr,
0729                                            const std::vector<unsigned int> &components) const {
0730   EGObjEmu egsta;
0731   egsta.clear();
0732   egsta.hwPt = ptCorr;
0733   egsta.hwEta = calo.hwEta;
0734   egsta.hwPhi = calo.hwPhi;
0735   egsta.hwQual = hwQual;
0736   egobjs.push_back(egsta);
0737 
0738   if (debug_ > 2)
0739     dbgCout() << "[REF] EGSta pt: " << egsta.hwPt << " eta: " << egsta.hwEta << " phi: " << egsta.hwPhi
0740               << " qual: " << std::bitset<4>(egsta.hwQual) << " packed: " << egsta.pack().to_string(16) << std::endl;
0741 
0742   return egobjs.back();
0743 }
0744 
0745 EGIsoObjEmu &PFTkEGAlgoEmulator::addEGIsoToPF(std::vector<EGIsoObjEmu> &egobjs,
0746                                               const EmCaloObjEmu &calo,
0747                                               const unsigned int hwQual,
0748                                               const pt_t ptCorr) const {
0749   EGIsoObjEmu egiso;
0750   egiso.clear();
0751   egiso.hwPt = ptCorr;
0752   egiso.hwEta = calo.hwEta;
0753   egiso.hwPhi = calo.hwPhi;
0754   unsigned int egHwQual = hwQual;
0755   if (cfg.doEndcapHwQual) {
0756     // 1. zero-suppress the loose EG-ID (bit 1)
0757     // 2. for now use the standalone tight definition (bit 0) to set the tight point for photons (bit 2)
0758     egHwQual = (hwQual & 0x9) | ((hwQual & 0x1) << 2);
0759   }
0760   egiso.hwQual = egHwQual;
0761   egiso.srcCluster = calo.src;
0762   egobjs.push_back(egiso);
0763 
0764   if (debug_ > 2)
0765     dbgCout() << "[REF] EGIsoObjEmu pt: " << egiso.hwPt << " eta: " << egiso.hwEta << " phi: " << egiso.hwPhi
0766               << " qual: " << std::bitset<4>(egiso.hwQual) << " packed: " << egiso.pack().to_string(16) << std::endl;
0767 
0768   return egobjs.back();
0769 }
0770 
0771 EGIsoEleObjEmu &PFTkEGAlgoEmulator::addEGIsoEleToPF(std::vector<EGIsoEleObjEmu> &egobjs,
0772                                                     const EmCaloObjEmu &calo,
0773                                                     const TkObjEmu &track,
0774                                                     const unsigned int hwQual,
0775                                                     const pt_t ptCorr,
0776                                                     const id_score_t bdtScore) const {
0777   EGIsoEleObjEmu egiso;
0778   egiso.clear();
0779   egiso.hwPt = ptCorr;
0780   egiso.hwEta = calo.hwEta;
0781   egiso.hwPhi = calo.hwPhi;
0782   unsigned int egHwQual = hwQual;
0783   if (cfg.algorithm == PFTkEGAlgoEmuConfig::Algo::compositeEE_v0 ||
0784       cfg.algorithm == PFTkEGAlgoEmuConfig::Algo::compositeEB_v0 ||
0785       cfg.algorithm == PFTkEGAlgoEmuConfig::Algo::compositeEE_v1 ||
0786       cfg.algorithm == PFTkEGAlgoEmuConfig::Algo::compositeEB_v1) {
0787     //Set tight WP
0788     bool is_tight = tkEleModel_->apply_wp_tight(bdtScore, egiso.floatPt());
0789     egHwQual = ((hwQual & 0x9) | (is_tight << 1));
0790 
0791   } else if (cfg.doEndcapHwQual) {
0792     // 1. zero-suppress the loose EG-ID (bit 1)
0793     // 2. for now use the standalone tight definition (bit 0) to set the tight point for eles (bit 1)
0794     egHwQual = (hwQual & 0x9) | ((hwQual & 0x1) << 1);
0795   }
0796 
0797   egiso.hwQual = egHwQual;
0798   egiso.hwDEta = track.hwVtxEta() - egiso.hwEta;
0799   egiso.hwDPhi = abs(track.hwVtxPhi() - egiso.hwPhi);
0800   egiso.hwZ0 = track.hwZ0;
0801   egiso.hwCharge = track.hwCharge;
0802   egiso.srcCluster = calo.src;
0803   egiso.srcTrack = track.src;
0804   egiso.hwIDScore = bdtScore;
0805   egobjs.push_back(egiso);
0806 
0807   if (debug_ > 2)
0808     dbgCout() << "[REF] EGIsoEleObjEmu pt: " << egiso.hwPt << " eta: " << egiso.hwEta << " phi: " << egiso.hwPhi
0809               << " qual: " << std::bitset<4>(egiso.hwQual) << " packed: " << egiso.pack().to_string(16) << std::endl;
0810 
0811   return egobjs.back();
0812 }
0813 
0814 void PFTkEGAlgoEmulator::addEgObjsToPF(std::vector<EGObjEmu> &egstas,
0815                                        std::vector<EGIsoObjEmu> &egobjs,
0816                                        std::vector<EGIsoEleObjEmu> &egeleobjs,
0817                                        const std::vector<EmCaloObjEmu> &emcalo,
0818                                        const std::vector<TkObjEmu> &track,
0819                                        const int calo_idx,
0820                                        const unsigned int hwQual,
0821                                        const pt_t ptCorr,
0822                                        const int tk_idx,
0823                                        const id_score_t bdtScore,
0824                                        const std::vector<unsigned int> &components) const {
0825   int src_idx = -1;
0826   if (writeEgSta()) {
0827     addEGStaToPF(egstas, emcalo[calo_idx], hwQual, ptCorr, components);
0828     src_idx = egstas.size() - 1;
0829   }
0830   EGIsoObjEmu &egobj = addEGIsoToPF(egobjs, emcalo[calo_idx], hwQual, ptCorr);
0831   egobj.src_idx = src_idx;
0832   if (tk_idx != -1) {
0833     EGIsoEleObjEmu &eleobj = addEGIsoEleToPF(egeleobjs, emcalo[calo_idx], track[tk_idx], hwQual, ptCorr, bdtScore);
0834     eleobj.src_idx = src_idx;
0835   }
0836 }
0837 
0838 void PFTkEGAlgoEmulator::runIso(const PFInputRegion &in,
0839                                 const std::vector<l1ct::PVObjEmu> &pvs,
0840                                 OutputRegion &out) const {
0841   if (cfg.doTkIso) {
0842     compute_isolation(out.egelectron, in.track, cfg.tkIsoParams_tkEle, pvs[0].hwZ0);
0843     compute_isolation(out.egphoton, in.track, cfg.tkIsoParams_tkEm, pvs[0].hwZ0);
0844   }
0845   if (cfg.doPfIso) {
0846     compute_isolation(out.egelectron, out.pfcharged, out.pfneutral, cfg.pfIsoParams_tkEle, pvs[0].hwZ0);
0847     compute_isolation(out.egphoton, out.pfcharged, out.pfneutral, cfg.pfIsoParams_tkEm, pvs[0].hwZ0);
0848   }
0849 
0850   std::for_each(out.egelectron.begin(), out.egelectron.end(), [&](EGIsoEleObjEmu &obj) {
0851     obj.hwIso = obj.hwIsoVar(cfg.hwIsoTypeTkEle);
0852   });
0853   std::for_each(
0854       out.egphoton.begin(), out.egphoton.end(), [&](EGIsoObjEmu &obj) { obj.hwIso = obj.hwIsoVar(cfg.hwIsoTypeTkEm); });
0855 }
0856 
0857 void PFTkEGAlgoEmulator::compute_isolation(std::vector<EGIsoObjEmu> &egobjs,
0858                                            const std::vector<TkObjEmu> &objects,
0859                                            const PFTkEGAlgoEmuConfig::IsoParameters &params,
0860                                            z0_t z0) const {
0861   for (int ic = 0, nc = egobjs.size(); ic < nc; ++ic) {
0862     auto &egphoton = egobjs[ic];
0863     iso_t sumPt = 0.;
0864     iso_t sumPtPV = 0.;
0865     compute_sumPt(sumPt, sumPtPV, objects, cfg.nTRACK, egphoton, params, z0);
0866     egphoton.setHwIso(EGIsoObjEmu::IsoType::TkIso, sumPt);
0867     egphoton.setHwIso(EGIsoObjEmu::IsoType::TkIsoPV, sumPtPV);
0868   }
0869 }
0870 
0871 void PFTkEGAlgoEmulator::compute_isolation(std::vector<EGIsoEleObjEmu> &egobjs,
0872                                            const std::vector<TkObjEmu> &objects,
0873                                            const PFTkEGAlgoEmuConfig::IsoParameters &params,
0874                                            z0_t z0) const {
0875   for (int ic = 0, nc = egobjs.size(); ic < nc; ++ic) {
0876     auto &egele = egobjs[ic];
0877     iso_t sumPt = 0.;
0878     iso_t sumPtPV = 0.;
0879     compute_sumPt(sumPt, sumPtPV, objects, cfg.nTRACK, egele, params, z0);
0880     egele.setHwIso(EGIsoEleObjEmu::IsoType::TkIso, sumPtPV);
0881   }
0882 }
0883 
0884 void PFTkEGAlgoEmulator::compute_isolation(std::vector<EGIsoObjEmu> &egobjs,
0885                                            const std::vector<PFChargedObjEmu> &charged,
0886                                            const std::vector<PFNeutralObjEmu> &neutrals,
0887                                            const PFTkEGAlgoEmuConfig::IsoParameters &params,
0888                                            z0_t z0) const {
0889   for (int ic = 0, nc = egobjs.size(); ic < nc; ++ic) {
0890     auto &egphoton = egobjs[ic];
0891     iso_t sumPt = 0.;
0892     iso_t sumPtPV = 0.;
0893     // FIXME: set max # of PF objects for iso
0894     compute_sumPt(sumPt, sumPtPV, charged, charged.size(), egphoton, params, z0);
0895     compute_sumPt(sumPt, sumPtPV, neutrals, neutrals.size(), egphoton, params, z0);
0896     egphoton.setHwIso(EGIsoObjEmu::IsoType::PfIso, sumPt);
0897     egphoton.setHwIso(EGIsoObjEmu::IsoType::PfIsoPV, sumPtPV);
0898   }
0899 }
0900 
0901 void PFTkEGAlgoEmulator::compute_isolation(std::vector<EGIsoEleObjEmu> &egobjs,
0902                                            const std::vector<PFChargedObjEmu> &charged,
0903                                            const std::vector<PFNeutralObjEmu> &neutrals,
0904                                            const PFTkEGAlgoEmuConfig::IsoParameters &params,
0905                                            z0_t z0) const {
0906   for (int ic = 0, nc = egobjs.size(); ic < nc; ++ic) {
0907     auto &egele = egobjs[ic];
0908     iso_t sumPt = 0.;
0909     iso_t sumPtPV = 0.;
0910     compute_sumPt(sumPt, sumPtPV, charged, charged.size(), egele, params, z0);
0911     compute_sumPt(sumPt, sumPtPV, neutrals, neutrals.size(), egele, params, z0);
0912     egele.setHwIso(EGIsoEleObjEmu::IsoType::PfIso, sumPtPV);
0913   }
0914 }
[ 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!