Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-06-10 01:53:56 

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