File indexing completed on 2024-04-06 12:21:31
0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/pf/pfalgo2hgc_ref.h"
0002 #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h"
0003
0004 #include <cmath>
0005 #include <cstdio>
0006 #include <algorithm>
0007 #include <memory>
0008
0009 #ifdef CMSSW_GIT_HASH
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0012
0013 l1ct::PFAlgo2HGCEmulator::PFAlgo2HGCEmulator(const edm::ParameterSet& iConfig)
0014 : PFAlgoEmulatorBase(iConfig.getParameter<uint32_t>("nTrack"),
0015 iConfig.getParameter<uint32_t>("nCalo"),
0016 iConfig.getParameter<uint32_t>("nMu"),
0017 iConfig.getParameter<uint32_t>("nSelCalo"),
0018 l1ct::Scales::makeDR2FromFloatDR(iConfig.getParameter<double>("trackMuDR")),
0019 l1ct::Scales::makeDR2FromFloatDR(iConfig.getParameter<double>("trackCaloDR")),
0020 l1ct::Scales::makePtFromFloat(iConfig.getParameter<double>("maxInvisiblePt")),
0021 l1ct::Scales::makePtFromFloat(iConfig.getParameter<double>("tightTrackMaxInvisiblePt"))) {
0022 debug_ = iConfig.getUntrackedParameter<bool>("debug", false);
0023 loadPtErrBins(iConfig);
0024 }
0025
0026 edm::ParameterSetDescription l1ct::PFAlgo2HGCEmulator::getParameterSetDescription() {
0027 edm::ParameterSetDescription description;
0028 description.add<unsigned int>("nTrack", 30);
0029 description.add<unsigned int>("nCalo", 20);
0030 description.add<unsigned int>("nMu", 4);
0031 description.add<unsigned int>("nSelCalo", 20);
0032 description.add<double>("trackMuDR", 0.2);
0033 description.add<double>("trackCaloDR", 0.1);
0034 description.add<double>("maxInvisiblePt", 10.0);
0035 description.add<double>("tightTrackMaxInvisiblePt", 20);
0036 addCaloResolutionParameterSetDescription(description);
0037 description.addUntracked<bool>("debug", false);
0038 return description;
0039 }
0040
0041 #endif
0042
0043 void l1ct::PFAlgo2HGCEmulator::toFirmware(const PFInputRegion& in,
0044 PFRegion& region,
0045 HadCaloObj calo[],
0046 TkObj track[],
0047 MuObj mu[]) const {
0048 region = in.region;
0049 l1ct::toFirmware(in.track, nTRACK_, track);
0050 l1ct::toFirmware(in.hadcalo, nCALO_, calo);
0051 l1ct::toFirmware(in.muon, nMU_, mu);
0052 }
0053
0054 void l1ct::PFAlgo2HGCEmulator::toFirmware(const OutputRegion& out,
0055 PFChargedObj outch[],
0056 PFNeutralObj outne[],
0057 PFChargedObj outmu[]) const {
0058 l1ct::toFirmware(out.pfcharged, nTRACK_, outch);
0059 l1ct::toFirmware(out.pfneutral, nSELCALO_, outne);
0060 l1ct::toFirmware(out.pfmuon, nMU_, outmu);
0061 }
0062
0063 void l1ct::PFAlgo2HGCEmulator::run(const PFInputRegion& in, OutputRegion& out) const {
0064 unsigned int nTRACK = std::min<unsigned>(nTRACK_, in.track.size());
0065 unsigned int nCALO = std::min<unsigned>(nCALO_, in.hadcalo.size());
0066 unsigned int nSELCALO = std::min<unsigned>(nSELCALO_, in.hadcalo.size());
0067 unsigned int nMU = std::min<unsigned>(nMU_, in.muon.size());
0068
0069 if (debug_) {
0070 dbgPrintf("FW\nFW \t region eta %+5.2f [ %+5.2f , %+5.2f ], phi %+5.2f [ %+5.2f , %+5.2f ] packed %s\n",
0071 in.region.floatEtaCenter(),
0072 in.region.floatEtaMinExtra(),
0073 in.region.floatEtaMaxExtra(),
0074 in.region.floatPhiCenter(),
0075 in.region.floatPhiCenter() - in.region.floatPhiHalfWidthExtra(),
0076 in.region.floatPhiCenter() + in.region.floatPhiHalfWidthExtra(),
0077 in.region.pack().to_string(16).c_str());
0078
0079 dbgPrintf("FW \t N(track) %3lu N(calo) %3lu N(mu) %3lu\n", in.track.size(), in.hadcalo.size(), in.muon.size());
0080
0081 for (unsigned int i = 0; i < nTRACK; ++i) {
0082 if (in.track[i].hwPt == 0)
0083 continue;
0084 dbgPrintf(
0085 "FW \t track %3d: pt %8.2f [ %8d ] calo eta %+5.2f [ %+5d ] calo phi %+5.2f [ %+5d ] vtx eta %+5.2f "
0086 "vtx phi %+5.2f charge %+2d qual %d fid %d glb eta %+5.2f phi %+5.2f packed %s\n",
0087 i,
0088 in.track[i].floatPt(),
0089 in.track[i].intPt(),
0090 in.track[i].floatEta(),
0091 in.track[i].intEta(),
0092 in.track[i].floatPhi(),
0093 in.track[i].intPhi(),
0094 in.track[i].floatVtxEta(),
0095 in.track[i].floatVtxPhi(),
0096 in.track[i].intCharge(),
0097 int(in.track[i].hwQuality),
0098 int(in.region.isFiducial(in.track[i].hwEta, in.track[i].hwPhi)),
0099 in.region.floatGlbEta(in.track[i].hwVtxEta()),
0100 in.region.floatGlbPhi(in.track[i].hwVtxPhi()),
0101 in.track[i].pack().to_string(16).c_str());
0102 }
0103 for (unsigned int i = 0; i < nCALO; ++i) {
0104 if (in.hadcalo[i].hwPt == 0)
0105 continue;
0106 dbgPrintf(
0107 "FW \t calo %3d: pt %8.2f [ %8d ] calo eta %+5.2f [ %+5d ] calo phi %+5.2f [ %+5d ] calo emPt %8.2f [ "
0108 "%6d ] emID %2d fid %d glb eta %+5.2f phi %+5.2f packed %s \n",
0109 i,
0110 in.hadcalo[i].floatPt(),
0111 in.hadcalo[i].intPt(),
0112 in.hadcalo[i].floatEta(),
0113 in.hadcalo[i].intEta(),
0114 in.hadcalo[i].floatPhi(),
0115 in.hadcalo[i].intPhi(),
0116 in.hadcalo[i].floatEmPt(),
0117 in.hadcalo[i].intEmPt(),
0118 in.hadcalo[i].hwEmID.to_int(),
0119 int(in.region.isFiducial(in.hadcalo[i].hwEta, in.hadcalo[i].hwPhi)),
0120 in.region.floatGlbEtaOf(in.hadcalo[i]),
0121 in.region.floatGlbPhiOf(in.hadcalo[i]),
0122 in.hadcalo[i].pack().to_string(16).c_str());
0123 }
0124 for (unsigned int i = 0; i < nMU; ++i) {
0125 if (in.muon[i].hwPt == 0)
0126 continue;
0127 dbgPrintf(
0128 "FW \t muon %3d: pt %8.2f [ %8d ] calo eta %+5.2f [ %+5d ] calo phi %+5.2f [ %+5d ] "
0129 "vtx eta %+5.2f vtx phi %+5.2f charge %+2d qual %2d glb eta %+5.2f phi %+5.2f packed %s \n",
0130 i,
0131 in.muon[i].floatPt(),
0132 in.muon[i].intPt(),
0133 in.muon[i].floatEta(),
0134 in.muon[i].intEta(),
0135 in.muon[i].floatPhi(),
0136 in.muon[i].intPhi(),
0137 in.muon[i].floatVtxEta(),
0138 in.muon[i].floatVtxPhi(),
0139 in.muon[i].intCharge(),
0140 int(in.muon[i].hwQuality),
0141 in.region.floatGlbEta(in.muon[i].hwVtxEta()),
0142 in.region.floatGlbPhi(in.muon[i].hwVtxPhi()),
0143 in.muon[i].pack().to_string(16).c_str());
0144 }
0145 }
0146
0147
0148
0149 std::vector<int> iMu;
0150 pfalgo_mu_ref(in, out, iMu);
0151
0152
0153
0154
0155
0156 std::vector<pt_t> calo_sumtk(nCALO), calo_subpt(nCALO);
0157 std::vector<pt2_t> calo_sumtkErr2(nCALO);
0158 for (unsigned int ic = 0; ic < nCALO; ++ic) {
0159 calo_sumtk[ic] = 0;
0160 calo_sumtkErr2[ic] = 0;
0161 }
0162
0163
0164 std::vector<bool> track_good(nTRACK, false);
0165 std::vector<bool> isEle(nTRACK, false);
0166 for (unsigned int it = 0; it < nTRACK; ++it) {
0167 if (!in.track[it].isPFLoose())
0168 continue;
0169 pt_t ptInv = in.track[it].isPFTight() ? tk_MAXINVPT_TIGHT_ : tk_MAXINVPT_LOOSE_;
0170 track_good[it] = (in.track[it].hwPt < ptInv) || (iMu[it] != -1);
0171 isEle[it] = false;
0172 }
0173
0174
0175 out.pfcharged.resize(nTRACK);
0176 out.pfneutral.resize(nSELCALO);
0177 for (unsigned int ipf = 0; ipf < nTRACK; ++ipf)
0178 out.pfcharged[ipf].clear();
0179 for (unsigned int ipf = 0; ipf < nSELCALO; ++ipf)
0180 out.pfneutral[ipf].clear();
0181
0182
0183 std::vector<int> tk2calo(nTRACK, -1);
0184 for (unsigned int it = 0; it < nTRACK; ++it) {
0185 if (in.track[it].hwPt > 0 && in.track[it].isPFLoose() && iMu[it] == -1) {
0186 pt_t tkCaloPtErr = ptErr_ref(in.region, in.track[it]);
0187 int ibest = best_match_with_pt_ref(dR2MAX_TK_CALO_, in.hadcalo, in.track[it], tkCaloPtErr);
0188 if (ibest != -1) {
0189 if (debug_)
0190 dbgPrintf("FW \t track %3d pt %8.2f caloPtErr %6.2f matched to calo %3d pt %8.2f\n",
0191 it,
0192 in.track[it].floatPt(),
0193 Scales::floatPt(tkCaloPtErr),
0194 ibest,
0195 in.hadcalo[ibest].floatPt());
0196 track_good[it] = true;
0197 isEle[it] = in.hadcalo[ibest].hwIsEM();
0198 calo_sumtk[ibest] += in.track[it].hwPt;
0199 calo_sumtkErr2[ibest] += tkCaloPtErr * tkCaloPtErr;
0200 }
0201 tk2calo[it] = ibest;
0202 }
0203 }
0204
0205 for (unsigned int ic = 0; ic < nCALO; ++ic) {
0206 if (calo_sumtk[ic] > 0) {
0207 pt_t ptdiff = in.hadcalo[ic].hwPt - calo_sumtk[ic];
0208 pt2_t sigmamult =
0209 calo_sumtkErr2[ic];
0210 if (debug_ && (in.hadcalo[ic].hwPt > 0)) {
0211 dbgPrintf(
0212 "FW \t calo %3d pt %8.2f [ %7d ] eta %+5.2f [ %+5d ] has a sum track pt %8.2f, difference %7.2f +- %.2f "
0213 "\n",
0214 ic,
0215 in.hadcalo[ic].floatPt(),
0216 in.hadcalo[ic].intPt(),
0217 in.hadcalo[ic].floatEta(),
0218 in.hadcalo[ic].intEta(),
0219 Scales::floatPt(calo_sumtk[ic]),
0220 Scales::floatPt(ptdiff),
0221 std::sqrt(Scales::floatPt(calo_sumtkErr2[ic])));
0222 }
0223 if (ptdiff > 0 && ptdiff * ptdiff > sigmamult) {
0224 calo_subpt[ic] = ptdiff;
0225 } else {
0226 calo_subpt[ic] = 0;
0227 }
0228 } else {
0229 calo_subpt[ic] = in.hadcalo[ic].hwPt;
0230 }
0231 if (debug_ && (in.hadcalo[ic].hwPt > 0))
0232 dbgPrintf(
0233 "FW \t calo' %3d pt %8.2f ---> %8.2f \n", ic, in.hadcalo[ic].floatPt(), Scales::floatPt(calo_subpt[ic]));
0234 }
0235
0236
0237 for (unsigned int it = 0; it < nTRACK; ++it) {
0238 if (in.track[it].hwPt > 0 && track_good[it]) {
0239 fillPFCand(in.track[it], out.pfcharged[it], (iMu[it] != -1), isEle[it]);
0240
0241 if (tk2calo[it] != -1)
0242 out.pfcharged[it].srcCluster = in.hadcalo[tk2calo[it]].src;
0243 if (iMu[it] != -1)
0244 out.pfcharged[it].srcMu = in.muon[iMu[it]].src;
0245 }
0246 }
0247
0248
0249 std::vector<PFNeutralObjEmu> outne_all(nCALO);
0250 for (unsigned int ipf = 0; ipf < nCALO; ++ipf)
0251 outne_all[ipf].clear();
0252 for (unsigned int ic = 0; ic < nCALO; ++ic) {
0253 if (calo_subpt[ic] > 0) {
0254 fillPFCand(in.hadcalo[ic], outne_all[ic], in.hadcalo[ic].hwIsEM());
0255 outne_all[ic].hwPt = calo_subpt[ic];
0256 outne_all[ic].hwEmPt = in.hadcalo[ic].hwIsEM() ? calo_subpt[ic] : pt_t(0);
0257 }
0258 }
0259
0260 if (nCALO_ == nSELCALO_) {
0261 std::swap(outne_all, out.pfneutral);
0262 } else {
0263 ptsort_ref(nCALO, nSELCALO, outne_all, out.pfneutral);
0264 }
0265
0266 if (debug_) {
0267 for (unsigned int i = 0; i < nTRACK; ++i) {
0268 if (out.pfcharged[i].hwPt == 0)
0269 continue;
0270 dbgPrintf(
0271 "FW \t outch %3d: pt %8.2f [ %8d ] calo eta %+5.2f [ %+5d ] calo phi %+5.2f [ %+5d ] pid %d packed %s\n",
0272 i,
0273 out.pfcharged[i].floatPt(),
0274 out.pfcharged[i].intPt(),
0275 out.pfcharged[i].floatEta(),
0276 out.pfcharged[i].intEta(),
0277 out.pfcharged[i].floatPhi(),
0278 out.pfcharged[i].intPhi(),
0279 out.pfcharged[i].intId(),
0280 out.pfcharged[i].pack().to_string(16).c_str());
0281 }
0282 for (unsigned int i = 0; i < nSELCALO; ++i) {
0283 if (out.pfneutral[i].hwPt == 0)
0284 continue;
0285 dbgPrintf(
0286 "FW \t outne %3d: pt %8.2f [ %8d ] calo eta %+5.2f [ %+5d ] calo phi %+5.2f [ %+5d ] pid %d packed %s\n",
0287 i,
0288 out.pfneutral[i].floatPt(),
0289 out.pfneutral[i].intPt(),
0290 out.pfneutral[i].floatEta(),
0291 out.pfneutral[i].intEta(),
0292 out.pfneutral[i].floatPhi(),
0293 out.pfneutral[i].intPhi(),
0294 out.pfneutral[i].intId(),
0295 out.pfneutral[i].pack().to_string(16).c_str());
0296 }
0297 }
0298 }