File indexing completed on 2022-10-01 01:02:14
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
0011 using namespace l1ct;
0012
0013 #ifdef CMSSW_GIT_HASH
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015
0016 l1ct::PFTkEGAlgoEmuConfig::PFTkEGAlgoEmuConfig(const edm::ParameterSet &pset)
0017 : nTRACK(pset.getParameter<uint32_t>("nTRACK")),
0018 nTRACK_EGIN(pset.getParameter<uint32_t>("nTRACK_EGIN")),
0019 nEMCALO_EGIN(pset.getParameter<uint32_t>("nEMCALO_EGIN")),
0020 nEM_EGOUT(pset.getParameter<uint32_t>("nEM_EGOUT")),
0021 filterHwQuality(pset.getParameter<bool>("filterHwQuality")),
0022 doBremRecovery(pset.getParameter<bool>("doBremRecovery")),
0023 writeBeforeBremRecovery(pset.getParameter<bool>("writeBeforeBremRecovery")),
0024 caloHwQual(pset.getParameter<int>("caloHwQual")),
0025 doEndcapHwQual(pset.getParameter<bool>("doEndcapHwQual")),
0026 emClusterPtMin(pset.getParameter<double>("caloEtMin")),
0027 dEtaMaxBrem(pset.getParameter<double>("dEtaMaxBrem")),
0028 dPhiMaxBrem(pset.getParameter<double>("dPhiMaxBrem")),
0029 absEtaBoundaries(pset.getParameter<std::vector<double>>("absEtaBoundaries")),
0030 dEtaValues(pset.getParameter<std::vector<double>>("dEtaValues")),
0031 dPhiValues(pset.getParameter<std::vector<double>>("dPhiValues")),
0032 trkQualityPtMin(pset.getParameter<double>("trkQualityPtMin")),
0033 writeEgSta(pset.getParameter<bool>("writeEGSta")),
0034 tkIsoParams_tkEle(pset.getParameter<edm::ParameterSet>("tkIsoParametersTkEle")),
0035 tkIsoParams_tkEm(pset.getParameter<edm::ParameterSet>("tkIsoParametersTkEm")),
0036 pfIsoParams_tkEle(pset.getParameter<edm::ParameterSet>("pfIsoParametersTkEle")),
0037 pfIsoParams_tkEm(pset.getParameter<edm::ParameterSet>("pfIsoParametersTkEm")),
0038 doTkIso(pset.getParameter<bool>("doTkIso")),
0039 doPfIso(pset.getParameter<bool>("doPfIso")),
0040 hwIsoTypeTkEle(static_cast<EGIsoEleObjEmu::IsoType>(pset.getParameter<uint32_t>("hwIsoTypeTkEle"))),
0041 hwIsoTypeTkEm(static_cast<EGIsoObjEmu::IsoType>(pset.getParameter<uint32_t>("hwIsoTypeTkEm"))),
0042 debug(pset.getUntrackedParameter<uint32_t>("debug", 0)) {}
0043
0044 l1ct::PFTkEGAlgoEmuConfig::IsoParameters::IsoParameters(const edm::ParameterSet &pset)
0045 : IsoParameters(pset.getParameter<double>("tkQualityPtMin"),
0046 pset.getParameter<double>("dZ"),
0047 pset.getParameter<double>("dRMin"),
0048 pset.getParameter<double>("dRMax")) {}
0049
0050 #endif
0051
0052 void PFTkEGAlgoEmulator::toFirmware(const PFInputRegion &in,
0053 PFRegion ®ion,
0054 EmCaloObj emcalo[],
0055 TkObj track[]) const {
0056 region = in.region;
0057 l1ct::toFirmware(in.track, cfg.nTRACK_EGIN, track);
0058 l1ct::toFirmware(in.emcalo, cfg.nEMCALO_EGIN, emcalo);
0059 if (debug_ > 0)
0060 dbgCout() << "# of inpput tracks: " << in.track.size() << " (max: " << cfg.nTRACK_EGIN << ")"
0061 << " emcalo: " << in.emcalo.size() << "(" << cfg.nEMCALO_EGIN << ")" << std::endl;
0062 }
0063
0064 void PFTkEGAlgoEmulator::toFirmware(const OutputRegion &out, EGIsoObj out_egphs[], EGIsoEleObj out_egeles[]) const {
0065 l1ct::toFirmware(out.egphoton, cfg.nEM_EGOUT, out_egphs);
0066 l1ct::toFirmware(out.egelectron, cfg.nEM_EGOUT, out_egeles);
0067 if (debug_ > 0)
0068 dbgCout() << "# output photons: " << out.egphoton.size() << " electrons: " << out.egelectron.size() << std::endl;
0069 }
0070
0071 void PFTkEGAlgoEmulator::toFirmware(
0072 const PFInputRegion &in, const l1ct::PVObjEmu &pvin, PFRegion ®ion, TkObj track[], PVObj &pv) const {
0073 region = in.region;
0074 l1ct::toFirmware(in.track, cfg.nTRACK, track);
0075 pv = pvin;
0076 if (debug_ > 0)
0077 dbgCout() << "# of inpput tracks: " << in.track.size() << " (max: " << cfg.nTRACK << ")" << std::endl;
0078 }
0079
0080 float PFTkEGAlgoEmulator::deltaPhi(float phi1, float phi2) const {
0081
0082 float x = phi1 - phi2;
0083 float o2pi = 1. / (2. * M_PI);
0084 if (std::abs(x) <= float(M_PI))
0085 return x;
0086 float n = std::round(x * o2pi);
0087 return x - n * float(2. * M_PI);
0088 }
0089
0090 void PFTkEGAlgoEmulator::link_emCalo2emCalo(const std::vector<EmCaloObjEmu> &emcalo,
0091 std::vector<int> &emCalo2emCalo) const {
0092
0093 for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0094 auto &calo = emcalo[ic];
0095 if (emCalo2emCalo[ic] != -1)
0096 continue;
0097
0098 for (int jc = ic + 1; jc < nc; ++jc) {
0099 if (emCalo2emCalo[jc] != -1)
0100 continue;
0101
0102 auto &otherCalo = emcalo[jc];
0103
0104 if (fabs(otherCalo.floatEta() - calo.floatEta()) < cfg.dEtaMaxBrem &&
0105 fabs(deltaPhi(otherCalo.floatPhi(), calo.floatPhi())) < cfg.dPhiMaxBrem) {
0106 emCalo2emCalo[jc] = ic;
0107 }
0108 }
0109 }
0110 }
0111
0112 void PFTkEGAlgoEmulator::link_emCalo2tk(const PFRegionEmu &r,
0113 const std::vector<EmCaloObjEmu> &emcalo,
0114 const std::vector<TkObjEmu> &track,
0115 std::vector<int> &emCalo2tk) const {
0116 unsigned int nTrackMax = std::min<unsigned>(track.size(), cfg.nTRACK_EGIN);
0117 for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0118 auto &calo = emcalo[ic];
0119
0120 float dPtMin = 999;
0121 for (unsigned int itk = 0; itk < nTrackMax; ++itk) {
0122 const auto &tk = track[itk];
0123 if (tk.floatPt() < cfg.trkQualityPtMin)
0124 continue;
0125
0126 float d_phi = deltaPhi(tk.floatPhi(), calo.floatPhi());
0127 float d_eta = tk.floatEta() - calo.floatEta();
0128
0129 auto eta_index =
0130 std::distance(cfg.absEtaBoundaries.begin(),
0131 std::lower_bound(
0132 cfg.absEtaBoundaries.begin(), cfg.absEtaBoundaries.end(), abs(r.floatGlbEta(calo.hwEta)))) -
0133 1;
0134
0135 float dEtaMax = cfg.dEtaValues[eta_index];
0136 float dPhiMax = cfg.dPhiValues[eta_index];
0137
0138 if ((((d_phi / dPhiMax) * (d_phi / dPhiMax)) + ((d_eta / dEtaMax) * (d_eta / dEtaMax))) < 1.) {
0139
0140 if (fabs(tk.floatPt() - calo.floatPt()) < dPtMin) {
0141 emCalo2tk[ic] = itk;
0142 dPtMin = fabs(tk.floatPt() - calo.floatPt());
0143 }
0144 }
0145 }
0146 }
0147 }
0148
0149 void PFTkEGAlgoEmulator::sel_emCalo(unsigned int nmax_sel,
0150 const std::vector<EmCaloObjEmu> &emcalo,
0151 std::vector<EmCaloObjEmu> &emcalo_sel) const {
0152 for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0153 const auto &calo = emcalo[ic];
0154 if ((calo.hwPt == 0) || (cfg.filterHwQuality && calo.hwEmID != cfg.caloHwQual) ||
0155 (calo.floatPt() < cfg.emClusterPtMin))
0156 continue;
0157 emcalo_sel.push_back(calo);
0158 if (emcalo_sel.size() >= nmax_sel)
0159 break;
0160 }
0161 }
0162
0163 void PFTkEGAlgoEmulator::run(const PFInputRegion &in, OutputRegion &out) const {
0164 if (debug_ > 1) {
0165 for (int ic = 0, nc = in.emcalo.size(); ic < nc; ++ic) {
0166 const auto &calo = in.emcalo[ic];
0167 if (calo.hwPt > 0)
0168 dbgCout() << "[REF] IN calo[" << ic << "] pt: " << calo.hwPt << " eta: " << calo.hwEta
0169 << " (glb eta: " << in.region.floatGlbEta(calo.hwEta) << ") phi: " << calo.hwPhi
0170 << "(glb phi: " << in.region.floatGlbPhi(calo.hwPhi) << ") qual: " << std::bitset<4>(calo.hwEmID)
0171 << std::endl;
0172 }
0173 }
0174
0175
0176
0177
0178 std::vector<EmCaloObjEmu> emcalo_sel;
0179 sel_emCalo(cfg.nEMCALO_EGIN, in.emcalo, emcalo_sel);
0180
0181 std::vector<int> emCalo2emCalo(emcalo_sel.size(), -1);
0182 if (cfg.doBremRecovery)
0183 link_emCalo2emCalo(emcalo_sel, emCalo2emCalo);
0184
0185 std::vector<int> emCalo2tk(emcalo_sel.size(), -1);
0186 link_emCalo2tk(in.region, emcalo_sel, in.track, emCalo2tk);
0187
0188 out.egsta.clear();
0189 std::vector<EGIsoObjEmu> egobjs;
0190 std::vector<EGIsoEleObjEmu> egeleobjs;
0191 eg_algo(in.region, emcalo_sel, in.track, emCalo2emCalo, emCalo2tk, out.egsta, egobjs, egeleobjs);
0192
0193 unsigned int nEGOut = std::min<unsigned>(cfg.nEM_EGOUT, egobjs.size());
0194 unsigned int nEGEleOut = std::min<unsigned>(cfg.nEM_EGOUT, egeleobjs.size());
0195
0196
0197 out.egphoton.clear();
0198 out.egelectron.clear();
0199 ptsort_ref(egobjs.size(), nEGOut, egobjs, out.egphoton);
0200 ptsort_ref(egeleobjs.size(), nEGEleOut, egeleobjs, out.egelectron);
0201 }
0202
0203 void PFTkEGAlgoEmulator::eg_algo(const PFRegionEmu ®ion,
0204 const std::vector<EmCaloObjEmu> &emcalo,
0205 const std::vector<TkObjEmu> &track,
0206 const std::vector<int> &emCalo2emCalo,
0207 const std::vector<int> &emCalo2tk,
0208 std::vector<EGObjEmu> &egstas,
0209 std::vector<EGIsoObjEmu> &egobjs,
0210 std::vector<EGIsoEleObjEmu> &egeleobjs) const {
0211 for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) {
0212 auto &calo = emcalo[ic];
0213
0214
0215 if (!region.isFiducial(calo))
0216 continue;
0217
0218 if (debug_ > 3)
0219 dbgCout() << "[REF] SEL emcalo with pt: " << calo.hwPt << " qual: " << calo.hwEmID << " eta: " << calo.hwEta
0220 << " phi " << calo.hwPhi << std::endl;
0221
0222 int itk = emCalo2tk[ic];
0223
0224
0225 if (!cfg.doBremRecovery || cfg.writeBeforeBremRecovery) {
0226
0227 unsigned int egQual = calo.hwEmID;
0228
0229
0230
0231 if (cfg.writeBeforeBremRecovery && cfg.doBremRecovery) {
0232 egQual = calo.hwEmID | 0x8;
0233 }
0234
0235 addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, egQual, calo.hwPt, itk);
0236 }
0237
0238 if (!cfg.doBremRecovery)
0239 continue;
0240
0241
0242 if (emCalo2emCalo[ic] != -1)
0243 continue;
0244
0245 pt_t ptBremReco = calo.hwPt;
0246 std::vector<unsigned int> components;
0247
0248 for (int jc = ic; jc < nc; ++jc) {
0249 if (emCalo2emCalo[jc] == ic) {
0250 auto &otherCalo = emcalo[jc];
0251 ptBremReco += otherCalo.hwPt;
0252 components.push_back(jc);
0253 }
0254 }
0255
0256
0257
0258 addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, calo.hwEmID, ptBremReco, itk, components);
0259 }
0260 }
0261
0262 EGObjEmu &PFTkEGAlgoEmulator::addEGStaToPF(std::vector<EGObjEmu> &egobjs,
0263 const EmCaloObjEmu &calo,
0264 const unsigned int hwQual,
0265 const pt_t ptCorr,
0266 const std::vector<unsigned int> &components) const {
0267 EGObjEmu egsta;
0268 egsta.clear();
0269 egsta.hwPt = ptCorr;
0270 egsta.hwEta = calo.hwEta;
0271 egsta.hwPhi = calo.hwPhi;
0272 egsta.hwQual = hwQual;
0273 egobjs.push_back(egsta);
0274
0275 if (debug_ > 2)
0276 dbgCout() << "[REF] EGSta pt: " << egsta.hwPt << " eta: " << egsta.hwEta << " phi: " << egsta.hwPhi
0277 << " qual: " << std::bitset<4>(egsta.hwQual) << " packed: " << egsta.pack().to_string(16) << std::endl;
0278
0279 return egobjs.back();
0280 }
0281
0282 EGIsoObjEmu &PFTkEGAlgoEmulator::addEGIsoToPF(std::vector<EGIsoObjEmu> &egobjs,
0283 const EmCaloObjEmu &calo,
0284 const unsigned int hwQual,
0285 const pt_t ptCorr) const {
0286 EGIsoObjEmu egiso;
0287 egiso.clear();
0288 egiso.hwPt = ptCorr;
0289 egiso.hwEta = calo.hwEta;
0290 egiso.hwPhi = calo.hwPhi;
0291 unsigned int egHwQual = hwQual;
0292 if (cfg.doEndcapHwQual) {
0293
0294
0295 egHwQual = (hwQual & 0x9) | ((hwQual & 0x1) << 2);
0296 }
0297 egiso.hwQual = egHwQual;
0298 egiso.srcCluster = calo.src;
0299 egobjs.push_back(egiso);
0300
0301 if (debug_ > 2)
0302 dbgCout() << "[REF] EGIsoObjEmu pt: " << egiso.hwPt << " eta: " << egiso.hwEta << " phi: " << egiso.hwPhi
0303 << " qual: " << std::bitset<4>(egiso.hwQual) << " packed: " << egiso.pack().to_string(16) << std::endl;
0304
0305 return egobjs.back();
0306 }
0307
0308 EGIsoEleObjEmu &PFTkEGAlgoEmulator::addEGIsoEleToPF(std::vector<EGIsoEleObjEmu> &egobjs,
0309 const EmCaloObjEmu &calo,
0310 const TkObjEmu &track,
0311 const unsigned int hwQual,
0312 const pt_t ptCorr) const {
0313 EGIsoEleObjEmu egiso;
0314 egiso.clear();
0315 egiso.hwPt = ptCorr;
0316 egiso.hwEta = calo.hwEta;
0317 egiso.hwPhi = calo.hwPhi;
0318 unsigned int egHwQual = hwQual;
0319 if (cfg.doEndcapHwQual) {
0320
0321
0322 egHwQual = (hwQual & 0x9) | ((hwQual & 0x1) << 1);
0323 }
0324 egiso.hwQual = egHwQual;
0325 egiso.hwDEta = track.hwVtxEta() - egiso.hwEta;
0326 egiso.hwDPhi = abs(track.hwVtxPhi() - egiso.hwPhi);
0327 egiso.hwZ0 = track.hwZ0;
0328 egiso.hwCharge = track.hwCharge;
0329 egiso.srcCluster = calo.src;
0330 egiso.srcTrack = track.src;
0331 egobjs.push_back(egiso);
0332
0333 if (debug_ > 2)
0334 dbgCout() << "[REF] EGIsoEleObjEmu pt: " << egiso.hwPt << " eta: " << egiso.hwEta << " phi: " << egiso.hwPhi
0335 << " qual: " << std::bitset<4>(egiso.hwQual) << " packed: " << egiso.pack().to_string(16) << std::endl;
0336
0337 return egobjs.back();
0338 }
0339
0340 void PFTkEGAlgoEmulator::addEgObjsToPF(std::vector<EGObjEmu> &egstas,
0341 std::vector<EGIsoObjEmu> &egobjs,
0342 std::vector<EGIsoEleObjEmu> &egeleobjs,
0343 const std::vector<EmCaloObjEmu> &emcalo,
0344 const std::vector<TkObjEmu> &track,
0345 const int calo_idx,
0346 const unsigned int hwQual,
0347 const pt_t ptCorr,
0348 const int tk_idx,
0349 const std::vector<unsigned int> &components) const {
0350 int sta_idx = -1;
0351 if (writeEgSta()) {
0352 addEGStaToPF(egstas, emcalo[calo_idx], hwQual, ptCorr, components);
0353 sta_idx = egstas.size() - 1;
0354 }
0355 EGIsoObjEmu &egobj = addEGIsoToPF(egobjs, emcalo[calo_idx], hwQual, ptCorr);
0356 egobj.sta_idx = sta_idx;
0357 if (tk_idx != -1) {
0358 EGIsoEleObjEmu &eleobj = addEGIsoEleToPF(egeleobjs, emcalo[calo_idx], track[tk_idx], hwQual, ptCorr);
0359 eleobj.sta_idx = sta_idx;
0360 }
0361 }
0362
0363 void PFTkEGAlgoEmulator::runIso(const PFInputRegion &in,
0364 const std::vector<l1ct::PVObjEmu> &pvs,
0365 OutputRegion &out) const {
0366 if (cfg.doTkIso) {
0367 compute_isolation(out.egelectron, in.track, cfg.tkIsoParams_tkEle, pvs[0].hwZ0);
0368 compute_isolation(out.egphoton, in.track, cfg.tkIsoParams_tkEm, pvs[0].hwZ0);
0369 }
0370 if (cfg.doPfIso) {
0371 compute_isolation(out.egelectron, out.pfcharged, out.pfneutral, cfg.pfIsoParams_tkEle, pvs[0].hwZ0);
0372 compute_isolation(out.egphoton, out.pfcharged, out.pfneutral, cfg.pfIsoParams_tkEm, pvs[0].hwZ0);
0373 }
0374
0375 std::for_each(out.egelectron.begin(), out.egelectron.end(), [&](EGIsoEleObjEmu &obj) {
0376 obj.hwIso = obj.hwIsoVar(cfg.hwIsoTypeTkEle);
0377 });
0378 std::for_each(
0379 out.egphoton.begin(), out.egphoton.end(), [&](EGIsoObjEmu &obj) { obj.hwIso = obj.hwIsoVar(cfg.hwIsoTypeTkEm); });
0380 }
0381
0382 void PFTkEGAlgoEmulator::compute_isolation(std::vector<EGIsoObjEmu> &egobjs,
0383 const std::vector<TkObjEmu> &objects,
0384 const PFTkEGAlgoEmuConfig::IsoParameters ¶ms,
0385 z0_t z0) const {
0386 for (int ic = 0, nc = egobjs.size(); ic < nc; ++ic) {
0387 auto &egphoton = egobjs[ic];
0388 iso_t sumPt = 0.;
0389 iso_t sumPtPV = 0.;
0390 compute_sumPt(sumPt, sumPtPV, objects, cfg.nTRACK, egphoton, params, z0);
0391 egphoton.setHwIso(EGIsoObjEmu::IsoType::TkIso, sumPt);
0392 egphoton.setHwIso(EGIsoObjEmu::IsoType::TkIsoPV, sumPtPV);
0393 }
0394 }
0395
0396 void PFTkEGAlgoEmulator::compute_isolation(std::vector<EGIsoEleObjEmu> &egobjs,
0397 const std::vector<TkObjEmu> &objects,
0398 const PFTkEGAlgoEmuConfig::IsoParameters ¶ms,
0399 z0_t z0) const {
0400 for (int ic = 0, nc = egobjs.size(); ic < nc; ++ic) {
0401 auto &egele = egobjs[ic];
0402 iso_t sumPt = 0.;
0403 iso_t sumPtPV = 0.;
0404 compute_sumPt(sumPt, sumPtPV, objects, cfg.nTRACK, egele, params, z0);
0405 egele.setHwIso(EGIsoEleObjEmu::IsoType::TkIso, sumPtPV);
0406 }
0407 }
0408
0409 void PFTkEGAlgoEmulator::compute_isolation(std::vector<EGIsoObjEmu> &egobjs,
0410 const std::vector<PFChargedObjEmu> &charged,
0411 const std::vector<PFNeutralObjEmu> &neutrals,
0412 const PFTkEGAlgoEmuConfig::IsoParameters ¶ms,
0413 z0_t z0) const {
0414 for (int ic = 0, nc = egobjs.size(); ic < nc; ++ic) {
0415 auto &egphoton = egobjs[ic];
0416 iso_t sumPt = 0.;
0417 iso_t sumPtPV = 0.;
0418
0419 compute_sumPt(sumPt, sumPtPV, charged, charged.size(), egphoton, params, z0);
0420 compute_sumPt(sumPt, sumPtPV, neutrals, neutrals.size(), egphoton, params, z0);
0421 egphoton.setHwIso(EGIsoObjEmu::IsoType::PfIso, sumPt);
0422 egphoton.setHwIso(EGIsoObjEmu::IsoType::PfIsoPV, sumPtPV);
0423 }
0424 }
0425
0426 void PFTkEGAlgoEmulator::compute_isolation(std::vector<EGIsoEleObjEmu> &egobjs,
0427 const std::vector<PFChargedObjEmu> &charged,
0428 const std::vector<PFNeutralObjEmu> &neutrals,
0429 const PFTkEGAlgoEmuConfig::IsoParameters ¶ms,
0430 z0_t z0) const {
0431 for (int ic = 0, nc = egobjs.size(); ic < nc; ++ic) {
0432 auto &egele = egobjs[ic];
0433 iso_t sumPt = 0.;
0434 iso_t sumPtPV = 0.;
0435 compute_sumPt(sumPt, sumPtPV, charged, charged.size(), egele, params, z0);
0436 compute_sumPt(sumPt, sumPtPV, neutrals, neutrals.size(), egele, params, z0);
0437 egele.setHwIso(EGIsoEleObjEmu::IsoType::PfIso, sumPtPV);
0438 }
0439 }