Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h"
0002 #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h"
0003 
0004 #include <cmath>
0005 #include <cstdio>
0006 #include <algorithm>
0007 #include <memory>
0008 #include <iostream>
0009 #include <bitset>
0010 #include <vector>
0011 
0012 using namespace l1ct;
0013 
0014 #ifdef CMSSW_GIT_HASH
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0017 
0018 l1ct::PFTkEGAlgoEmuConfig::PFTkEGAlgoEmuConfig(const edm::ParameterSet &pset)
0019     : nTRACK(pset.getParameter<uint32_t>("nTRACK")),
0020       nTRACK_EGIN(pset.getParameter<uint32_t>("nTRACK_EGIN")),
0021       nEMCALO_EGIN(pset.getParameter<uint32_t>("nEMCALO_EGIN")),
0022       nEM_EGOUT(pset.getParameter<uint32_t>("nEM_EGOUT")),
0023       filterHwQuality(pset.getParameter<bool>("filterHwQuality")),
0024       doBremRecovery(pset.getParameter<bool>("doBremRecovery")),
0025       writeBeforeBremRecovery(pset.getParameter<bool>("writeBeforeBremRecovery")),
0026       caloHwQual(pset.getParameter<int>("caloHwQual")),
0027       doEndcapHwQual(pset.getParameter<bool>("doEndcapHwQual")),
0028       emClusterPtMin(pset.getParameter<double>("caloEtMin")),
0029       dEtaMaxBrem(pset.getParameter<double>("dEtaMaxBrem")),
0030       dPhiMaxBrem(pset.getParameter<double>("dPhiMaxBrem")),
0031       absEtaBoundaries(pset.getParameter<std::vector<double>>("absEtaBoundaries")),
0032       dEtaValues(pset.getParameter<std::vector<double>>("dEtaValues")),
0033       dPhiValues(pset.getParameter<std::vector<double>>("dPhiValues")),
0034       trkQualityPtMin(pset.getParameter<double>("trkQualityPtMin")),
0035       doCompositeTkEle(pset.getParameter<bool>("doCompositeTkEle")),
0036       nCompCandPerCluster(pset.getParameter<uint32_t>("nCompCandPerCluster")),
0037       writeEgSta(pset.getParameter<bool>("writeEGSta")),
0038       tkIsoParams_tkEle(pset.getParameter<edm::ParameterSet>("tkIsoParametersTkEle")),
0039       tkIsoParams_tkEm(pset.getParameter<edm::ParameterSet>("tkIsoParametersTkEm")),
0040       pfIsoParams_tkEle(pset.getParameter<edm::ParameterSet>("pfIsoParametersTkEle")),
0041       pfIsoParams_tkEm(pset.getParameter<edm::ParameterSet>("pfIsoParametersTkEm")),
0042       doTkIso(pset.getParameter<bool>("doTkIso")),
0043       doPfIso(pset.getParameter<bool>("doPfIso")),
0044       hwIsoTypeTkEle(static_cast<EGIsoEleObjEmu::IsoType>(pset.getParameter<uint32_t>("hwIsoTypeTkEle"))),
0045       hwIsoTypeTkEm(static_cast<EGIsoObjEmu::IsoType>(pset.getParameter<uint32_t>("hwIsoTypeTkEm"))),
0046       compIDparams(pset.getParameter<edm::ParameterSet>("compositeParametersTkEle")),
0047       debug(pset.getUntrackedParameter<uint32_t>("debug", 0)) {}
0048 
0049 edm::ParameterSetDescription l1ct::PFTkEGAlgoEmuConfig::getParameterSetDescription() {
0050   edm::ParameterSetDescription description;
0051   description.add<unsigned int>("nTRACK");
0052   description.add<unsigned int>("nTRACK_EGIN");
0053   description.add<unsigned int>("nEMCALO_EGIN");
0054   description.add<unsigned int>("nEM_EGOUT");
0055   description.add<bool>("doBremRecovery", false);
0056   description.add<bool>("writeBeforeBremRecovery", false);
0057   description.add<bool>("filterHwQuality", false);
0058   description.add<int>("caloHwQual", 4);
0059   description.add<bool>("doEndcapHwQual", false);
0060   description.add<double>("dEtaMaxBrem", 0.02);
0061   description.add<double>("dPhiMaxBrem", 0.1);
0062   description.add<std::vector<double>>("absEtaBoundaries",
0063                                        {
0064                                            0.0,
0065                                            0.9,
0066                                            1.5,
0067                                        });
0068   description.add<std::vector<double>>("dEtaValues",
0069                                        {
0070                                            0.025,
0071                                            0.015,
0072                                            0.01,
0073                                        });
0074   description.add<std::vector<double>>("dPhiValues",
0075                                        {
0076                                            0.07,
0077                                            0.07,
0078                                            0.07,
0079                                        });
0080   description.add<double>("caloEtMin", 0.0);
0081   description.add<double>("trkQualityPtMin", 10.0);
0082   description.add<bool>("writeEGSta", false);
0083   description.add<edm::ParameterSetDescription>("tkIsoParametersTkEm", IsoParameters::getParameterSetDescription());
0084   description.add<edm::ParameterSetDescription>("tkIsoParametersTkEle", IsoParameters::getParameterSetDescription());
0085   description.add<edm::ParameterSetDescription>("pfIsoParametersTkEm", IsoParameters::getParameterSetDescription());
0086   description.add<edm::ParameterSetDescription>("pfIsoParametersTkEle", IsoParameters::getParameterSetDescription());
0087   description.add<bool>("doTkIso", true);
0088   description.add<bool>("doPfIso", true);
0089   description.add<unsigned int>("hwIsoTypeTkEle", 0);
0090   description.add<unsigned int>("hwIsoTypeTkEm", 2);
0091   description.add<bool>("doCompositeTkEle", false);
0092   description.add<unsigned int>("nCompCandPerCluster", 3);
0093   description.add<edm::ParameterSetDescription>("compositeParametersTkEle",
0094                                                 CompIDParameters::getParameterSetDescription());
0095   return description;
0096 }
0097 
0098 l1ct::PFTkEGAlgoEmuConfig::IsoParameters::IsoParameters(const edm::ParameterSet &pset)
0099     : IsoParameters(pset.getParameter<double>("tkQualityPtMin"),
0100                     pset.getParameter<double>("dZ"),
0101                     pset.getParameter<double>("dRMin"),
0102                     pset.getParameter<double>("dRMax")) {}
0103 
0104 edm::ParameterSetDescription l1ct::PFTkEGAlgoEmuConfig::IsoParameters::getParameterSetDescription() {
0105   edm::ParameterSetDescription description;
0106   description.add<double>("tkQualityPtMin");
0107   description.add<double>("dZ", 0.6);
0108   description.add<double>("dRMin");
0109   description.add<double>("dRMax");
0110   return description;
0111 }
0112 
0113 l1ct::PFTkEGAlgoEmuConfig::CompIDParameters::CompIDParameters(const edm::ParameterSet &pset)
0114     : CompIDParameters(pset.getParameter<double>("loose_wp"),
0115                        pset.getParameter<double>("tight_wp"),
0116                        pset.getParameter<std::string>("model")) {}
0117 
0118 edm::ParameterSetDescription l1ct::PFTkEGAlgoEmuConfig::CompIDParameters::getParameterSetDescription() {
0119   edm::ParameterSetDescription description;
0120   description.add<double>("loose_wp", -0.732422);
0121   description.add<double>("tight_wp", 0.214844);
0122   description.add<std::string>("model", "L1Trigger/Phase2L1ParticleFlow/data/compositeID.json");
0123   return description;
0124 }
0125 #endif
0126 
0127 PFTkEGAlgoEmulator::PFTkEGAlgoEmulator(const PFTkEGAlgoEmuConfig &config)
0128     : cfg(config), composite_bdt_(nullptr), debug_(cfg.debug) {
0129   if (cfg.doCompositeTkEle) {
0130 #ifdef CMSSW_GIT_HASH
0131     auto resolvedFileName = edm::FileInPath(cfg.compIDparams.conifer_model).fullPath();
0132 #else
0133     auto resolvedFileName = cfg.compIDparams.conifer_model;
0134 #endif
0135     composite_bdt_ = new conifer::BDT<bdt_feature_t, bdt_score_t, false>(resolvedFileName);
0136   }
0137 }
0138 
0139 void PFTkEGAlgoEmulator::toFirmware(const PFInputRegion &in,
0140                                     PFRegion &region,
0141                                     EmCaloObj emcalo[/*nCALO*/],
0142                                     TkObj track[/*nTRACK*/]) const {
0143   region = in.region;
0144   l1ct::toFirmware(in.track, cfg.nTRACK_EGIN, track);
0145   l1ct::toFirmware(in.emcalo, cfg.nEMCALO_EGIN, emcalo);
0146   if (debug_ > 0)
0147     dbgCout() << "# of inpput tracks: " << in.track.size() << " (max: " << cfg.nTRACK_EGIN << ")"
0148               << " emcalo: " << in.emcalo.size() << "(" << cfg.nEMCALO_EGIN << ")" << std::endl;
0149 }
0150 
0151 void PFTkEGAlgoEmulator::toFirmware(const OutputRegion &out, EGIsoObj out_egphs[], EGIsoEleObj out_egeles[]) const {
0152   l1ct::toFirmware(out.egphoton, cfg.nEM_EGOUT, out_egphs);
0153   l1ct::toFirmware(out.egelectron, cfg.nEM_EGOUT, out_egeles);
0154   if (debug_ > 0)
0155     dbgCout() << "# output photons: " << out.egphoton.size() << " electrons: " << out.egelectron.size() << std::endl;
0156 }
0157 
0158 void PFTkEGAlgoEmulator::toFirmware(
0159     const PFInputRegion &in, const l1ct::PVObjEmu &pvin, PFRegion &region, TkObj track[/*nTRACK*/], PVObj &pv) const {
0160   region = in.region;
0161   l1ct::toFirmware(in.track, cfg.nTRACK, track);
0162   pv = pvin;
0163   if (debug_ > 0)
0164     dbgCout() << "# of inpput tracks: " << in.track.size() << " (max: " << cfg.nTRACK << ")" << std::endl;
0165 }
0166 
0167 float PFTkEGAlgoEmulator::deltaPhi(float phi1, float phi2) const {
0168   // reduce to [-pi,pi]
0169   float x = phi1 - phi2;
0170   float o2pi = 1. / (2. * M_PI);
0171   if (std::abs(x) <= float(M_PI))
0172     return x;
0173   float n = std::round(x * o2pi);
0174   return x - n * float(2. * M_PI);
0175 }
0176 
0177 void PFTkEGAlgoEmulator::link_emCalo2emCalo(const std::vector<EmCaloObjEmu> &emcalo,
0178                                             std::vector<int> &emCalo2emCalo) const {
0179   // NOTE: we assume the input to be sorted!!!
0180   for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0181     auto &calo = emcalo[ic];
0182     if (emCalo2emCalo[ic] != -1)
0183       continue;
0184 
0185     for (int jc = ic + 1; jc < nc; ++jc) {
0186       if (emCalo2emCalo[jc] != -1)
0187         continue;
0188 
0189       auto &otherCalo = emcalo[jc];
0190 
0191       if (fabs(otherCalo.floatEta() - calo.floatEta()) < cfg.dEtaMaxBrem &&
0192           fabs(deltaPhi(otherCalo.floatPhi(), calo.floatPhi())) < cfg.dPhiMaxBrem) {
0193         emCalo2emCalo[jc] = ic;
0194       }
0195     }
0196   }
0197 }
0198 
0199 void PFTkEGAlgoEmulator::link_emCalo2tk_elliptic(const PFRegionEmu &r,
0200                                                  const std::vector<EmCaloObjEmu> &emcalo,
0201                                                  const std::vector<TkObjEmu> &track,
0202                                                  std::vector<int> &emCalo2tk) const {
0203   unsigned int nTrackMax = std::min<unsigned>(track.size(), cfg.nTRACK_EGIN);
0204   for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0205     auto &calo = emcalo[ic];
0206 
0207     float dPtMin = 999;
0208     for (unsigned int itk = 0; itk < nTrackMax; ++itk) {
0209       const auto &tk = track[itk];
0210       if (tk.floatPt() < cfg.trkQualityPtMin)
0211         continue;
0212 
0213       float d_phi = deltaPhi(tk.floatPhi(), calo.floatPhi());
0214       float d_eta = tk.floatEta() - calo.floatEta();  // We only use it squared
0215 
0216       auto eta_index =
0217           std::distance(cfg.absEtaBoundaries.begin(),
0218                         std::lower_bound(
0219                             cfg.absEtaBoundaries.begin(), cfg.absEtaBoundaries.end(), abs(r.floatGlbEta(calo.hwEta)))) -
0220           1;
0221 
0222       float dEtaMax = cfg.dEtaValues[eta_index];
0223       float dPhiMax = cfg.dPhiValues[eta_index];
0224 
0225       if (debug_ > 2 && calo.hwPt > 0) {
0226         dbgCout() << "[REF] tried to link calo " << ic << " (pt " << calo.intPt() << ", eta " << calo.intEta()
0227                   << ", phi " << calo.intPhi() << ") "
0228                   << " to tk " << itk << " (pt " << tk.intPt() << ", eta " << tk.intEta() << ", phi " << tk.intPhi()
0229                   << "): "
0230                   << " eta_index " << eta_index << ", "
0231                   << " dEta " << d_eta << " (max " << dEtaMax << "), dPhi " << d_phi << " (max " << dPhiMax << ") "
0232                   << " ellipse = "
0233                   << (((d_phi / dPhiMax) * (d_phi / dPhiMax)) + ((d_eta / dEtaMax) * (d_eta / dEtaMax))) << "\n";
0234       }
0235       if ((((d_phi / dPhiMax) * (d_phi / dPhiMax)) + ((d_eta / dEtaMax) * (d_eta / dEtaMax))) < 1.) {
0236         // NOTE: for now we implement only best pt match. This is NOT what is done in the L1TkElectronTrackProducer
0237         if (fabs(tk.floatPt() - calo.floatPt()) < dPtMin) {
0238           emCalo2tk[ic] = itk;
0239           dPtMin = fabs(tk.floatPt() - calo.floatPt());
0240         }
0241       }
0242     }
0243   }
0244 }
0245 
0246 void PFTkEGAlgoEmulator::link_emCalo2tk_composite(const PFRegionEmu &r,
0247                                                   const std::vector<EmCaloObjEmu> &emcalo,
0248                                                   const std::vector<TkObjEmu> &track,
0249                                                   std::vector<int> &emCalo2tk,
0250                                                   std::vector<id_score_t> &emCaloTkBdtScore) const {
0251   unsigned int nTrackMax = std::min<unsigned>(track.size(), cfg.nTRACK_EGIN);
0252   for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0253     auto &calo = emcalo[ic];
0254 
0255     std::vector<CompositeCandidate> candidates;
0256 
0257     for (unsigned int itk = 0; itk < nTrackMax; ++itk) {
0258       const auto &tk = track[itk];
0259       if (tk.floatPt() <= cfg.trkQualityPtMin)
0260         continue;
0261 
0262       float d_phi = deltaPhi(tk.floatPhi(), calo.floatPhi());
0263       float d_eta = tk.floatEta() - calo.floatEta();  // We only use it squared
0264       float dR2 = (d_phi * d_phi) + (d_eta * d_eta);
0265 
0266       if (dR2 < 0.04) {
0267         // Only store indices, dR and dpT for now. The other quantities are computed only for the best nCandPerCluster.
0268         CompositeCandidate cand;
0269         cand.cluster_idx = ic;
0270         cand.track_idx = itk;
0271         cand.dpt = std::abs(tk.floatPt() - calo.floatPt());
0272         candidates.push_back(cand);
0273       }
0274     }
0275     // FIXME: find best sort criteria, for now we use dpt
0276     std::sort(candidates.begin(),
0277               candidates.end(),
0278               [](const CompositeCandidate &a, const CompositeCandidate &b) -> bool { return a.dpt < b.dpt; });
0279     unsigned int nCandPerCluster = std::min<unsigned int>(candidates.size(), cfg.nCompCandPerCluster);
0280     if (nCandPerCluster == 0)
0281       continue;
0282 
0283     id_score_t maxScore = -pow(2, l1ct::id_score_t::iwidth - 1);
0284     int ibest = -1;
0285     for (unsigned int icand = 0; icand < nCandPerCluster; icand++) {
0286       auto &cand = candidates[icand];
0287       const std::vector<EmCaloObjEmu> &emcalo_sel = emcalo;
0288       id_score_t score = compute_composite_score(cand, emcalo_sel, track, cfg.compIDparams);
0289       if ((score > cfg.compIDparams.bdtScore_loose_wp) && (score > maxScore)) {
0290         maxScore = score;
0291         ibest = icand;
0292       }
0293     }
0294     if (ibest != -1) {
0295       emCalo2tk[ic] = candidates[ibest].track_idx;
0296       emCaloTkBdtScore[ic] = maxScore;
0297     }
0298   }
0299 }
0300 
0301 id_score_t PFTkEGAlgoEmulator::compute_composite_score(CompositeCandidate &cand,
0302                                                        const std::vector<EmCaloObjEmu> &emcalo,
0303                                                        const std::vector<TkObjEmu> &track,
0304                                                        const PFTkEGAlgoEmuConfig::CompIDParameters &params) const {
0305   // Get the cluster/track objects that form the composite candidate
0306   const auto &calo = emcalo[cand.cluster_idx];
0307   const auto &tk = track[cand.track_idx];
0308 
0309   // Prepare the input features
0310   bdt_feature_t hoe = calo.hwHoe;
0311   bdt_feature_t tkpt = tk.hwPt;
0312   bdt_feature_t srrtot = calo.hwSrrTot;
0313   bdt_feature_t deta = tk.hwEta - calo.hwEta;
0314   ap_ufixed<18, 0> calo_invPt = l1ct::invert_with_shift<pt_t, ap_ufixed<18, 0>, 1024>(calo.hwPt);
0315   bdt_feature_t dpt = tk.hwPt * calo_invPt;
0316   bdt_feature_t meanz = calo.hwMeanZ;
0317   bdt_feature_t dphi = tk.hwPhi - calo.hwPhi;
0318   bdt_feature_t nstubs = tk.hwStubs;
0319   bdt_feature_t chi2rphi = tk.hwRedChi2RPhi;
0320   bdt_feature_t chi2rz = tk.hwRedChi2RZ;
0321   bdt_feature_t chi2bend = tk.hwRedChi2Bend;
0322 
0323   // Run BDT inference
0324   std::vector<bdt_feature_t> inputs = {tkpt, hoe, srrtot, deta, dphi, dpt, meanz, nstubs, chi2rphi, chi2rz, chi2bend};
0325   std::vector<bdt_score_t> bdt_score = composite_bdt_->decision_function(inputs);
0326 
0327   return bdt_score[0] / 4;
0328 }
0329 
0330 void PFTkEGAlgoEmulator::sel_emCalo(unsigned int nmax_sel,
0331                                     const std::vector<EmCaloObjEmu> &emcalo,
0332                                     std::vector<EmCaloObjEmu> &emcalo_sel) const {
0333   for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0334     const auto &calo = emcalo[ic];
0335     if ((calo.hwPt == 0) || (cfg.filterHwQuality && calo.hwEmID != cfg.caloHwQual) ||
0336         (calo.floatPt() < cfg.emClusterPtMin))
0337       continue;
0338     emcalo_sel.push_back(calo);
0339     if (emcalo_sel.size() >= nmax_sel)
0340       break;
0341   }
0342 }
0343 
0344 void PFTkEGAlgoEmulator::run(const PFInputRegion &in, OutputRegion &out) const {
0345   if (debug_ > 1) {
0346     for (int ic = 0, nc = in.emcalo.size(); ic < nc; ++ic) {
0347       const auto &calo = in.emcalo[ic];
0348       if (calo.hwPt > 0)
0349         dbgCout() << "[REF] IN calo[" << ic << "] pt: " << calo.hwPt << " eta: " << calo.hwEta
0350                   << " (glb eta: " << in.region.floatGlbEta(calo.hwEta) << ") phi: " << calo.hwPhi
0351                   << "(glb phi: " << in.region.floatGlbPhi(calo.hwPhi) << ") qual: " << std::bitset<4>(calo.hwEmID)
0352                   << std::endl;
0353     }
0354   }
0355   // FIXME: can be removed in the endcap since now running with the "interceptor".
0356   // Might still be needed in barrel
0357   // filter and select first N elements of input clusters
0358   std::vector<EmCaloObjEmu> emcalo_sel;
0359   sel_emCalo(cfg.nEMCALO_EGIN, in.emcalo, emcalo_sel);
0360 
0361   std::vector<int> emCalo2emCalo(emcalo_sel.size(), -1);
0362   if (cfg.doBremRecovery)
0363     link_emCalo2emCalo(emcalo_sel, emCalo2emCalo);
0364 
0365   std::vector<int> emCalo2tk(emcalo_sel.size(), -1);
0366   std::vector<id_score_t> emCaloTkBdtScore(emcalo_sel.size(), 0);
0367 
0368   if (cfg.doCompositeTkEle) {
0369     link_emCalo2tk_composite(in.region, emcalo_sel, in.track, emCalo2tk, emCaloTkBdtScore);
0370   } else {
0371     link_emCalo2tk_elliptic(in.region, emcalo_sel, in.track, emCalo2tk);
0372   }
0373 
0374   out.egsta.clear();
0375   std::vector<EGIsoObjEmu> egobjs;
0376   std::vector<EGIsoEleObjEmu> egeleobjs;
0377   eg_algo(in.region, emcalo_sel, in.track, emCalo2emCalo, emCalo2tk, emCaloTkBdtScore, out.egsta, egobjs, egeleobjs);
0378 
0379   unsigned int nEGOut = std::min<unsigned>(cfg.nEM_EGOUT, egobjs.size());
0380   unsigned int nEGEleOut = std::min<unsigned>(cfg.nEM_EGOUT, egeleobjs.size());
0381 
0382   // init output containers
0383   out.egphoton.clear();
0384   out.egelectron.clear();
0385   ptsort_ref(egobjs.size(), nEGOut, egobjs, out.egphoton);
0386   ptsort_ref(egeleobjs.size(), nEGEleOut, egeleobjs, out.egelectron);
0387 }
0388 
0389 void PFTkEGAlgoEmulator::eg_algo(const PFRegionEmu &region,
0390                                  const std::vector<EmCaloObjEmu> &emcalo,
0391                                  const std::vector<TkObjEmu> &track,
0392                                  const std::vector<int> &emCalo2emCalo,
0393                                  const std::vector<int> &emCalo2tk,
0394                                  const std::vector<id_score_t> &emCaloTkBdtScore,
0395                                  std::vector<EGObjEmu> &egstas,
0396                                  std::vector<EGIsoObjEmu> &egobjs,
0397                                  std::vector<EGIsoEleObjEmu> &egeleobjs) const {
0398   for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0399     auto &calo = emcalo[ic];
0400 
0401     // discard immediately EG objects that would not fall in the fiducial eta-phi region
0402     if (!region.isFiducial(calo))
0403       continue;
0404 
0405     if (debug_ > 3)
0406       dbgCout() << "[REF] SEL emcalo with pt: " << calo.hwPt << " qual: " << calo.hwEmID << " eta: " << calo.hwEta
0407                 << " phi " << calo.hwPhi << std::endl;
0408 
0409     int itk = emCalo2tk[ic];
0410     const id_score_t &bdt = emCaloTkBdtScore[ic];
0411 
0412     // check if brem recovery is on
0413     if (!cfg.doBremRecovery || cfg.writeBeforeBremRecovery) {
0414       // 1. create EG objects before brem recovery
0415       unsigned int egQual = calo.hwEmID;
0416       // If we write both objects with and without brem-recovery
0417       // bit 3 is used for the brem-recovery bit: if set = no recovery
0418       // (for consistency with the barrel hwQual where by default the brem recovery is done upstream)
0419       if (cfg.writeBeforeBremRecovery && cfg.doBremRecovery) {
0420         egQual = calo.hwEmID | 0x8;
0421       }
0422 
0423       addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, egQual, calo.hwPt, itk, bdt);
0424     }
0425 
0426     if (!cfg.doBremRecovery)
0427       continue;
0428 
0429     // check if the cluster has already been used in a brem reclustering
0430     if (emCalo2emCalo[ic] != -1)
0431       continue;
0432 
0433     pt_t ptBremReco = calo.hwPt;
0434     std::vector<unsigned int> components;
0435 
0436     for (int jc = ic; jc < nc; ++jc) {
0437       if (emCalo2emCalo[jc] == ic) {
0438         auto &otherCalo = emcalo[jc];
0439         ptBremReco += otherCalo.hwPt;
0440         components.push_back(jc);
0441       }
0442     }
0443 
0444     // 2. create EG objects with brem recovery
0445     addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, calo.hwEmID, ptBremReco, itk, bdt, components);
0446   }
0447 }
0448 
0449 EGObjEmu &PFTkEGAlgoEmulator::addEGStaToPF(std::vector<EGObjEmu> &egobjs,
0450                                            const EmCaloObjEmu &calo,
0451                                            const unsigned int hwQual,
0452                                            const pt_t ptCorr,
0453                                            const std::vector<unsigned int> &components) const {
0454   EGObjEmu egsta;
0455   egsta.clear();
0456   egsta.hwPt = ptCorr;
0457   egsta.hwEta = calo.hwEta;
0458   egsta.hwPhi = calo.hwPhi;
0459   egsta.hwQual = hwQual;
0460   egobjs.push_back(egsta);
0461 
0462   if (debug_ > 2)
0463     dbgCout() << "[REF] EGSta pt: " << egsta.hwPt << " eta: " << egsta.hwEta << " phi: " << egsta.hwPhi
0464               << " qual: " << std::bitset<4>(egsta.hwQual) << " packed: " << egsta.pack().to_string(16) << std::endl;
0465 
0466   return egobjs.back();
0467 }
0468 
0469 EGIsoObjEmu &PFTkEGAlgoEmulator::addEGIsoToPF(std::vector<EGIsoObjEmu> &egobjs,
0470                                               const EmCaloObjEmu &calo,
0471                                               const unsigned int hwQual,
0472                                               const pt_t ptCorr) const {
0473   EGIsoObjEmu egiso;
0474   egiso.clear();
0475   egiso.hwPt = ptCorr;
0476   egiso.hwEta = calo.hwEta;
0477   egiso.hwPhi = calo.hwPhi;
0478   unsigned int egHwQual = hwQual;
0479   if (cfg.doEndcapHwQual) {
0480     // 1. zero-suppress the loose EG-ID (bit 1)
0481     // 2. for now use the standalone tight definition (bit 0) to set the tight point for photons (bit 2)
0482     egHwQual = (hwQual & 0x9) | ((hwQual & 0x1) << 2);
0483   }
0484   egiso.hwQual = egHwQual;
0485   egiso.srcCluster = calo.src;
0486   egobjs.push_back(egiso);
0487 
0488   if (debug_ > 2)
0489     dbgCout() << "[REF] EGIsoObjEmu pt: " << egiso.hwPt << " eta: " << egiso.hwEta << " phi: " << egiso.hwPhi
0490               << " qual: " << std::bitset<4>(egiso.hwQual) << " packed: " << egiso.pack().to_string(16) << std::endl;
0491 
0492   return egobjs.back();
0493 }
0494 
0495 EGIsoEleObjEmu &PFTkEGAlgoEmulator::addEGIsoEleToPF(std::vector<EGIsoEleObjEmu> &egobjs,
0496                                                     const EmCaloObjEmu &calo,
0497                                                     const TkObjEmu &track,
0498                                                     const unsigned int hwQual,
0499                                                     const pt_t ptCorr,
0500                                                     const id_score_t bdtScore) const {
0501   EGIsoEleObjEmu egiso;
0502   egiso.clear();
0503   egiso.hwPt = ptCorr;
0504   egiso.hwEta = calo.hwEta;
0505   egiso.hwPhi = calo.hwPhi;
0506   unsigned int egHwQual = hwQual;
0507   if (cfg.doEndcapHwQual) {
0508     if (cfg.doCompositeTkEle) {
0509       // tight ele WP is set for tight BDT score
0510       egHwQual = (hwQual & 0x9) | ((bdtScore >= cfg.compIDparams.bdtScore_tight_wp) << 1);
0511     } else {
0512       // 1. zero-suppress the loose EG-ID (bit 1)
0513       // 2. for now use the standalone tight definition (bit 0) to set the tight point for eles (bit 1)
0514       egHwQual = (hwQual & 0x9) | ((hwQual & 0x1) << 1);
0515     }
0516   }
0517   egiso.hwQual = egHwQual;
0518   egiso.hwDEta = track.hwVtxEta() - egiso.hwEta;
0519   egiso.hwDPhi = abs(track.hwVtxPhi() - egiso.hwPhi);
0520   egiso.hwZ0 = track.hwZ0;
0521   egiso.hwCharge = track.hwCharge;
0522   egiso.srcCluster = calo.src;
0523   egiso.srcTrack = track.src;
0524   egiso.hwIDScore = bdtScore;
0525   egobjs.push_back(egiso);
0526 
0527   if (debug_ > 2)
0528     dbgCout() << "[REF] EGIsoEleObjEmu pt: " << egiso.hwPt << " eta: " << egiso.hwEta << " phi: " << egiso.hwPhi
0529               << " qual: " << std::bitset<4>(egiso.hwQual) << " packed: " << egiso.pack().to_string(16) << std::endl;
0530 
0531   return egobjs.back();
0532 }
0533 
0534 void PFTkEGAlgoEmulator::addEgObjsToPF(std::vector<EGObjEmu> &egstas,
0535                                        std::vector<EGIsoObjEmu> &egobjs,
0536                                        std::vector<EGIsoEleObjEmu> &egeleobjs,
0537                                        const std::vector<EmCaloObjEmu> &emcalo,
0538                                        const std::vector<TkObjEmu> &track,
0539                                        const int calo_idx,
0540                                        const unsigned int hwQual,
0541                                        const pt_t ptCorr,
0542                                        const int tk_idx,
0543                                        const id_score_t bdtScore,
0544                                        const std::vector<unsigned int> &components) const {
0545   int src_idx = -1;
0546   if (writeEgSta()) {
0547     addEGStaToPF(egstas, emcalo[calo_idx], hwQual, ptCorr, components);
0548     src_idx = egstas.size() - 1;
0549   }
0550   EGIsoObjEmu &egobj = addEGIsoToPF(egobjs, emcalo[calo_idx], hwQual, ptCorr);
0551   egobj.src_idx = src_idx;
0552   if (tk_idx != -1) {
0553     EGIsoEleObjEmu &eleobj = addEGIsoEleToPF(egeleobjs, emcalo[calo_idx], track[tk_idx], hwQual, ptCorr, bdtScore);
0554     eleobj.src_idx = src_idx;
0555   }
0556 }
0557 
0558 void PFTkEGAlgoEmulator::runIso(const PFInputRegion &in,
0559                                 const std::vector<l1ct::PVObjEmu> &pvs,
0560                                 OutputRegion &out) const {
0561   if (cfg.doTkIso) {
0562     compute_isolation(out.egelectron, in.track, cfg.tkIsoParams_tkEle, pvs[0].hwZ0);
0563     compute_isolation(out.egphoton, in.track, cfg.tkIsoParams_tkEm, pvs[0].hwZ0);
0564   }
0565   if (cfg.doPfIso) {
0566     compute_isolation(out.egelectron, out.pfcharged, out.pfneutral, cfg.pfIsoParams_tkEle, pvs[0].hwZ0);
0567     compute_isolation(out.egphoton, out.pfcharged, out.pfneutral, cfg.pfIsoParams_tkEm, pvs[0].hwZ0);
0568   }
0569 
0570   std::for_each(out.egelectron.begin(), out.egelectron.end(), [&](EGIsoEleObjEmu &obj) {
0571     obj.hwIso = obj.hwIsoVar(cfg.hwIsoTypeTkEle);
0572   });
0573   std::for_each(
0574       out.egphoton.begin(), out.egphoton.end(), [&](EGIsoObjEmu &obj) { obj.hwIso = obj.hwIsoVar(cfg.hwIsoTypeTkEm); });
0575 }
0576 
0577 void PFTkEGAlgoEmulator::compute_isolation(std::vector<EGIsoObjEmu> &egobjs,
0578                                            const std::vector<TkObjEmu> &objects,
0579                                            const PFTkEGAlgoEmuConfig::IsoParameters &params,
0580                                            z0_t z0) const {
0581   for (int ic = 0, nc = egobjs.size(); ic < nc; ++ic) {
0582     auto &egphoton = egobjs[ic];
0583     iso_t sumPt = 0.;
0584     iso_t sumPtPV = 0.;
0585     compute_sumPt(sumPt, sumPtPV, objects, cfg.nTRACK, egphoton, params, z0);
0586     egphoton.setHwIso(EGIsoObjEmu::IsoType::TkIso, sumPt);
0587     egphoton.setHwIso(EGIsoObjEmu::IsoType::TkIsoPV, sumPtPV);
0588   }
0589 }
0590 
0591 void PFTkEGAlgoEmulator::compute_isolation(std::vector<EGIsoEleObjEmu> &egobjs,
0592                                            const std::vector<TkObjEmu> &objects,
0593                                            const PFTkEGAlgoEmuConfig::IsoParameters &params,
0594                                            z0_t z0) const {
0595   for (int ic = 0, nc = egobjs.size(); ic < nc; ++ic) {
0596     auto &egele = egobjs[ic];
0597     iso_t sumPt = 0.;
0598     iso_t sumPtPV = 0.;
0599     compute_sumPt(sumPt, sumPtPV, objects, cfg.nTRACK, egele, params, z0);
0600     egele.setHwIso(EGIsoEleObjEmu::IsoType::TkIso, sumPtPV);
0601   }
0602 }
0603 
0604 void PFTkEGAlgoEmulator::compute_isolation(std::vector<EGIsoObjEmu> &egobjs,
0605                                            const std::vector<PFChargedObjEmu> &charged,
0606                                            const std::vector<PFNeutralObjEmu> &neutrals,
0607                                            const PFTkEGAlgoEmuConfig::IsoParameters &params,
0608                                            z0_t z0) const {
0609   for (int ic = 0, nc = egobjs.size(); ic < nc; ++ic) {
0610     auto &egphoton = egobjs[ic];
0611     iso_t sumPt = 0.;
0612     iso_t sumPtPV = 0.;
0613     // FIXME: set max # of PF objects for iso
0614     compute_sumPt(sumPt, sumPtPV, charged, charged.size(), egphoton, params, z0);
0615     compute_sumPt(sumPt, sumPtPV, neutrals, neutrals.size(), egphoton, params, z0);
0616     egphoton.setHwIso(EGIsoObjEmu::IsoType::PfIso, sumPt);
0617     egphoton.setHwIso(EGIsoObjEmu::IsoType::PfIsoPV, sumPtPV);
0618   }
0619 }
0620 
0621 void PFTkEGAlgoEmulator::compute_isolation(std::vector<EGIsoEleObjEmu> &egobjs,
0622                                            const std::vector<PFChargedObjEmu> &charged,
0623                                            const std::vector<PFNeutralObjEmu> &neutrals,
0624                                            const PFTkEGAlgoEmuConfig::IsoParameters &params,
0625                                            z0_t z0) const {
0626   for (int ic = 0, nc = egobjs.size(); ic < nc; ++ic) {
0627     auto &egele = egobjs[ic];
0628     iso_t sumPt = 0.;
0629     iso_t sumPtPV = 0.;
0630     compute_sumPt(sumPt, sumPtPV, charged, charged.size(), egele, params, z0);
0631     compute_sumPt(sumPt, sumPtPV, neutrals, neutrals.size(), egele, params, z0);
0632     egele.setHwIso(EGIsoEleObjEmu::IsoType::PfIso, sumPtPV);
0633   }
0634 }
[ 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!