Back to home page

Project CMSSW displayed by LXR

 
 

    


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[/*nCALO*/],
0046                                           TkObj track[/*nTRACK*/],
0047                                           MuObj mu[/*nMU*/]) 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[/*nTRACK*/],
0056                                           PFNeutralObj outne[/*nSELCALO*/],
0057                                           PFChargedObj outmu[/*nMU*/]) 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   // TK-MU Linking
0149   std::vector<int> iMu;
0150   pfalgo_mu_ref(in, out, iMu);
0151 
0152   ////////////////////////////////////////////////////
0153   // TK-HAD Linking
0154 
0155   // initialize sum track pt
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   // initialize good track bit
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   // initialize output
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   // for each track, find the closest calo
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;  // for emulator info
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];  //  + (calo_sumtkErr2[ic] >> 1)); // this multiplies by 1.5 = sqrt(1.5)^2 ~ (1.2)^2
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   // copy out charged hadrons
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], /*isMu=*/(iMu[it] != -1), isEle[it]);
0240       // extra emulator information
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   // copy out neutral hadrons with sorting and cropping
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);  // FIXME
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 }