Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-10 01:54:00

0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/PFAlgo3.h"
0002 #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h"
0003 
0004 #include "DataFormats/L1TParticleFlow/interface/PFCandidate.h"
0005 
0006 #include "FWCore/Utilities/interface/Exception.h"
0007 
0008 #include "DataFormats/Math/interface/deltaR.h"
0009 
0010 namespace {
0011   template <typename T1, typename T2>
0012   float floatDR(const T1 &t1, const T2 &t2) {
0013     return deltaR(t1.floatEta(), t1.floatPhi(), t2.floatEta(), t2.floatPhi());
0014   }
0015 }  // namespace
0016 
0017 using namespace l1tpf_impl;
0018 
0019 PFAlgo3::PFAlgo3(const edm::ParameterSet &iConfig) : PFAlgoBase(iConfig) {
0020   debug_ = iConfig.getUntrackedParameter<int>("debugPFAlgo3", iConfig.getUntrackedParameter<int>("debug", 0));
0021   edm::ParameterSet linkcfg = iConfig.getParameter<edm::ParameterSet>("linking");
0022   drMatchMu_ = linkcfg.getParameter<double>("trackMuDR");
0023 
0024   std::string muMatchMode = linkcfg.getParameter<std::string>("trackMuMatch");
0025   if (muMatchMode == "boxBestByPtRatio")
0026     muMatchMode_ = MuMatchMode::BoxBestByPtRatio;
0027   else if (muMatchMode == "drBestByPtRatio")
0028     muMatchMode_ = MuMatchMode::DrBestByPtRatio;
0029   else if (muMatchMode == "drBestByPtDiff")
0030     muMatchMode_ = MuMatchMode::DrBestByPtDiff;
0031   else
0032     throw cms::Exception("Configuration", "bad value for trackMuMatch configurable");
0033 
0034   std::string tkCaloLinkMetric = linkcfg.getParameter<std::string>("trackCaloLinkMetric");
0035   if (tkCaloLinkMetric == "bestByDR")
0036     tkCaloLinkMetric_ = TkCaloLinkMetric::BestByDR;
0037   else if (tkCaloLinkMetric == "bestByDRPt")
0038     tkCaloLinkMetric_ = TkCaloLinkMetric::BestByDRPt;
0039   else if (tkCaloLinkMetric == "bestByDR2Pt2")
0040     tkCaloLinkMetric_ = TkCaloLinkMetric::BestByDR2Pt2;
0041   else
0042     throw cms::Exception("Configuration", "bad value for tkCaloLinkMetric configurable");
0043 
0044   drMatch_ = linkcfg.getParameter<double>("trackCaloDR");
0045   ptMatchLow_ = linkcfg.getParameter<double>("trackCaloNSigmaLow");
0046   ptMatchHigh_ = linkcfg.getParameter<double>("trackCaloNSigmaHigh");
0047   useTrackCaloSigma_ = linkcfg.getParameter<bool>("useTrackCaloSigma");
0048   maxInvisiblePt_ = linkcfg.getParameter<double>("maxInvisiblePt");
0049 
0050   drMatchEm_ = linkcfg.getParameter<double>("trackEmDR");
0051   trackEmUseAlsoTrackSigma_ = linkcfg.getParameter<bool>("trackEmUseAlsoTrackSigma");
0052   trackEmMayUseCaloMomenta_ = linkcfg.getParameter<bool>("trackEmMayUseCaloMomenta");
0053   emCaloUseAlsoCaloSigma_ = linkcfg.getParameter<bool>("emCaloUseAlsoCaloSigma");
0054   ptMinFracMatchEm_ = linkcfg.getParameter<double>("caloEmPtMinFrac");
0055   drMatchEmHad_ = linkcfg.getParameter<double>("emCaloDR");
0056   emHadSubtractionPtSlope_ = linkcfg.getParameter<double>("emCaloSubtractionPtSlope");
0057   caloReLinkStep_ = linkcfg.getParameter<bool>("caloReLink");
0058   caloReLinkDr_ = linkcfg.getParameter<double>("caloReLinkDR");
0059   caloReLinkThreshold_ = linkcfg.getParameter<double>("caloReLinkThreshold");
0060   rescaleTracks_ = linkcfg.getParameter<bool>("rescaleTracks");
0061   caloTrkWeightedAverage_ = linkcfg.getParameter<bool>("useCaloTrkWeightedAverage");
0062   sumTkCaloErr2_ = linkcfg.getParameter<bool>("sumTkCaloErr2");
0063   ecalPriority_ = linkcfg.getParameter<bool>("ecalPriority");
0064   tightTrackMaxInvisiblePt_ = linkcfg.getParameter<double>("tightTrackMaxInvisiblePt");
0065   sortInputs_ = iConfig.getParameter<bool>("sortInputs");
0066 }
0067 
0068 void PFAlgo3::runPF(Region &r) const {
0069   initRegion(r, sortInputs_);
0070 
0071   /// ------------- first step (can all go in parallel) ----------------
0072 
0073   if (debug_) {
0074     dbgPrintf(
0075         "PFAlgo3\nPFAlgo3 region eta [ %+5.2f , %+5.2f ], phi [ %+5.2f , %+5.2f ], fiducial eta [ %+5.2f , %+5.2f ], "
0076         "phi [ %+5.2f , %+5.2f ]\n",
0077         r.etaMin - r.etaExtra,
0078         r.etaMax + r.etaExtra,
0079         r.phiCenter - r.phiHalfWidth - r.phiExtra,
0080         r.phiCenter + r.phiHalfWidth + r.phiExtra,
0081         r.etaMin,
0082         r.etaMax,
0083         r.phiCenter - r.phiHalfWidth,
0084         r.phiCenter + r.phiHalfWidth);
0085     dbgPrintf("PFAlgo3 \t N(track) %3lu   N(em) %3lu   N(calo) %3lu   N(mu) %3lu\n",
0086               r.track.size(),
0087               r.emcalo.size(),
0088               r.calo.size(),
0089               r.muon.size());
0090     for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0091       const auto &tk = r.track[itk];
0092       dbgPrintf(
0093           "PFAlgo3 \t track %3d: pt %7.2f +- %5.2f  vtx eta %+5.2f  vtx phi %+5.2f  calo eta %+5.2f  calo phi %+5.2f  "
0094           "fid %1d  calo ptErr %7.2f stubs %2d chi2 %7.1f quality %d\n",
0095           itk,
0096           tk.floatPt(),
0097           tk.floatPtErr(),
0098           tk.floatVtxEta(),
0099           tk.floatVtxPhi(),
0100           tk.floatEta(),
0101           tk.floatPhi(),
0102           int(r.fiducialLocal(tk.floatEta(), tk.floatPhi())),
0103           tk.floatCaloPtErr(),
0104           int(tk.hwStubs),
0105           tk.hwChi2 * 0.1f,
0106           int(tk.hwFlags));
0107     }
0108     for (int iem = 0, nem = r.emcalo.size(); iem < nem; ++iem) {
0109       const auto &em = r.emcalo[iem];
0110       dbgPrintf(
0111           "PFAlgo3 \t EM    %3d: pt %7.2f +- %5.2f  vtx eta %+5.2f  vtx phi %+5.2f  calo eta %+5.2f  calo phi %+5.2f  "
0112           "fid %1d  calo ptErr %7.2f\n",
0113           iem,
0114           em.floatPt(),
0115           em.floatPtErr(),
0116           em.floatEta(),
0117           em.floatPhi(),
0118           em.floatEta(),
0119           em.floatPhi(),
0120           int(r.fiducialLocal(em.floatEta(), em.floatPhi())),
0121           em.floatPtErr());
0122     }
0123     for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0124       auto &calo = r.calo[ic];
0125       dbgPrintf(
0126           "PFAlgo3 \t calo  %3d: pt %7.2f +- %5.2f  vtx eta %+5.2f  vtx phi %+5.2f  calo eta %+5.2f  calo phi %+5.2f  "
0127           "fid %1d  calo ptErr %7.2f em pt %7.2f \n",
0128           ic,
0129           calo.floatPt(),
0130           calo.floatPtErr(),
0131           calo.floatEta(),
0132           calo.floatPhi(),
0133           calo.floatEta(),
0134           calo.floatPhi(),
0135           int(r.fiducialLocal(calo.floatEta(), calo.floatPhi())),
0136           calo.floatPtErr(),
0137           calo.floatEmPt());
0138     }
0139     for (int im = 0, nm = r.muon.size(); im < nm; ++im) {
0140       auto &mu = r.muon[im];
0141       dbgPrintf(
0142           "PFAlgo3 \t muon  %3d: pt %7.2f           vtx eta %+5.2f  vtx phi %+5.2f  calo eta %+5.2f  calo phi %+5.2f  "
0143           "fid %1d \n",
0144           im,
0145           mu.floatPt(),
0146           mu.floatEta(),
0147           mu.floatPhi(),
0148           mu.floatEta(),
0149           mu.floatPhi(),
0150           int(r.fiducialLocal(mu.floatEta(), mu.floatPhi())));
0151     }
0152   }
0153 
0154   std::vector<int> tk2mu(r.track.size(), -1), mu2tk(r.muon.size(), -1);
0155   link_tk2mu(r, tk2mu, mu2tk);
0156 
0157   // match all tracks to the closest EM cluster
0158   std::vector<int> tk2em(r.track.size(), -1);
0159   link_tk2em(r, tk2em);
0160 
0161   // match all em to the closest had (can happen in parallel to the above)
0162   std::vector<int> em2calo(r.emcalo.size(), -1);
0163   link_em2calo(r, em2calo);
0164 
0165   /// ------------- next step (needs the previous) ----------------
0166   // for each EM cluster, count and add up the pt of all the corresponding tracks (skipping muons)
0167   std::vector<int> em2ntk(r.emcalo.size(), 0);
0168   std::vector<float> em2sumtkpt(r.emcalo.size(), 0);
0169   std::vector<float> em2sumtkpterr(r.emcalo.size(), 0);
0170   sum_tk2em(r, tk2em, em2ntk, em2sumtkpt, em2sumtkpterr);
0171 
0172   /// ------------- next step (needs the previous) ----------------
0173   // process ecal clusters after linking
0174   emcalo_algo(r, em2ntk, em2sumtkpt, em2sumtkpterr);
0175 
0176   /// ------------- next step (needs the previous) ----------------
0177   // promote all flagged tracks to electrons
0178   emtk_algo(r, tk2em, em2ntk, em2sumtkpterr);
0179   sub_em2calo(r, em2calo);
0180 
0181   /// ------------- next step (needs the previous) ----------------
0182   // track to calo matching (first iteration, with a lower bound on the calo pt; there may be another one later)
0183   std::vector<int> tk2calo(r.track.size(), -1);
0184   link_tk2calo(r, tk2calo);
0185 
0186   /// ------------- next step (needs the previous) ----------------
0187   // for each calo, compute the sum of the track pt
0188   std::vector<int> calo2ntk(r.calo.size(), 0);
0189   std::vector<float> calo2sumtkpt(r.calo.size(), 0);
0190   std::vector<float> calo2sumtkpterr(r.calo.size(), 0);
0191   sum_tk2calo(r, tk2calo, calo2ntk, calo2sumtkpt, calo2sumtkpterr);
0192 
0193   // in the meantime, promote unlinked low pt tracks to hadrons
0194   unlinkedtk_algo(r, tk2calo);
0195 
0196   /// ------------- next step (needs the previous) ----------------
0197   /// OPTIONAL STEP: try to recover split hadron showers (v1.0):
0198   //     off by default, as it seems to not do much in jets even if it helps remove tails in single-pion events
0199   if (caloReLinkStep_)
0200     calo_relink(r, calo2ntk, calo2sumtkpt, calo2sumtkpterr);
0201 
0202   /// ------------- next step (needs the previous) ----------------
0203   // process matched calo clusters, compare energy to sum track pt
0204   std::vector<float> calo2alpha(r.calo.size(), 1);
0205   linkedcalo_algo(r, calo2ntk, calo2sumtkpt, calo2sumtkpterr, calo2alpha);
0206 
0207   /// ------------- next step (needs the previous) ----------------
0208   /// process matched tracks, if necessary rescale or average
0209   linkedtk_algo(r, tk2calo, calo2ntk, calo2alpha);
0210   // process unmatched calo clusters
0211   unlinkedcalo_algo(r);
0212   // finally do muons
0213   save_muons(r, tk2mu);
0214 }
0215 
0216 void PFAlgo3::link_tk2mu(Region &r, std::vector<int> &tk2mu, std::vector<int> &mu2tk) const {
0217   // do a rectangular match for the moment; make a box of the same are as a 0.2 cone
0218   int intDrMuonMatchBox = std::ceil(drMatchMu_ * CaloCluster::ETAPHI_SCALE * std::sqrt(M_PI / 4));
0219   for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0220     tk2mu[itk] = false;
0221   }
0222   for (int imu = 0, nmu = r.muon.size(); imu < nmu; ++imu) {
0223     const auto &mu = r.muon[imu];
0224     if (debug_)
0225       dbgPrintf(
0226           "PFAlgo3 \t muon  %3d (pt %7.2f, eta %+5.2f, phi %+5.2f) \n", imu, mu.floatPt(), mu.floatEta(), mu.floatPhi());
0227     float minDistance = 9e9;
0228     switch (muMatchMode_) {
0229       case MuMatchMode::BoxBestByPtRatio:
0230         minDistance = 4.;
0231         break;
0232       case MuMatchMode::DrBestByPtRatio:
0233         minDistance = 4.;
0234         break;
0235       case MuMatchMode::DrBestByPtDiff:
0236         minDistance = 0.5 * mu.floatPt();
0237         break;
0238     }
0239     int imatch = -1;
0240     for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0241       const auto &tk = r.track[itk];
0242       if (!tk.quality(l1tpf_impl::InputTrack::PFLOOSE))
0243         continue;
0244       int deta = std::abs(mu.hwEta - tk.hwEta);
0245       int dphi = std::abs((mu.hwPhi - tk.hwPhi) % CaloCluster::PHI_WRAP);
0246       float dr = floatDR(mu, tk);
0247       float dpt = std::abs(mu.floatPt() - tk.floatPt());
0248       float dptr = (mu.hwPt > tk.hwPt ? mu.floatPt() / tk.floatPt() : tk.floatPt() / mu.floatPt());
0249       bool ok = false;
0250       float distance = 9e9;
0251       switch (muMatchMode_) {
0252         case MuMatchMode::BoxBestByPtRatio:
0253           ok = (deta < intDrMuonMatchBox) && (dphi < intDrMuonMatchBox);
0254           distance = dptr;
0255           break;
0256         case MuMatchMode::DrBestByPtRatio:
0257           ok = (dr < drMatchMu_);
0258           distance = dptr;
0259           break;
0260         case MuMatchMode::DrBestByPtDiff:
0261           ok = (dr < drMatchMu_);
0262           distance = dpt;
0263           break;
0264       }
0265       if (debug_ && dr < 0.4) {
0266         dbgPrintf(
0267             "PFAlgo3 \t\t possible match with track %3d (pt %7.2f, caloeta %+5.2f, calophi %+5.2f, dr %.2f, eta "
0268             "%+5.2f, phi %+5.2f, dr %.2f):  angular %1d, distance %.3f (vs %.3f)\n",
0269             itk,
0270             tk.floatPt(),
0271             tk.floatEta(),
0272             tk.floatPhi(),
0273             dr,
0274             tk.floatVtxEta(),
0275             tk.floatVtxPhi(),
0276             deltaR(mu.floatEta(), mu.floatPhi(), tk.floatVtxEta(), tk.floatVtxPhi()),
0277             (ok ? 1 : 0),
0278             distance,
0279             minDistance);
0280       }
0281       if (!ok)
0282         continue;
0283       // FIXME for the moment, we do the floating point matching in pt
0284       if (distance < minDistance) {
0285         minDistance = distance;
0286         imatch = itk;
0287       }
0288     }
0289     if (debug_ && imatch > -1)
0290       dbgPrintf("PFAlgo3 \t muon  %3d (pt %7.2f) linked to track %3d (pt %7.2f)\n",
0291                 imu,
0292                 mu.floatPt(),
0293                 imatch,
0294                 r.track[imatch].floatPt());
0295     if (debug_ && imatch == -1)
0296       dbgPrintf("PFAlgo3 \t muon  %3d (pt %7.2f) not linked to any track\n", imu, mu.floatPt());
0297     mu2tk[imu] = imatch;
0298     if (imatch > -1) {
0299       tk2mu[imatch] = imu;
0300       r.track[imatch].muonLink = true;
0301     }
0302   }
0303 }
0304 
0305 void PFAlgo3::link_tk2em(Region &r, std::vector<int> &tk2em) const {
0306   // match all tracks to the closest EM cluster
0307   for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0308     const auto &tk = r.track[itk];
0309     if (!tk.quality(l1tpf_impl::InputTrack::PFLOOSE))
0310       continue;
0311     //if (tk.muonLink) continue; // not necessary I think
0312     float drbest = drMatchEm_;
0313     for (int iem = 0, nem = r.emcalo.size(); iem < nem; ++iem) {
0314       const auto &em = r.emcalo[iem];
0315       float dr = floatDR(tk, em);
0316       if (dr < drbest) {
0317         tk2em[itk] = iem;
0318         drbest = dr;
0319       }
0320     }
0321     if (debug_ && tk2em[itk] != -1)
0322       dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) matches to EM   %3d (pt %7.2f) with dr %.3f\n",
0323                 itk,
0324                 tk.floatPt(),
0325                 tk2em[itk],
0326                 tk2em[itk] == -1 ? 0.0 : r.emcalo[tk2em[itk]].floatPt(),
0327                 drbest);
0328   }
0329 }
0330 
0331 void PFAlgo3::link_em2calo(Region &r, std::vector<int> &em2calo) const {
0332   // match all em to the closest had (can happen in parallel to the above)
0333   for (int iem = 0, nem = r.emcalo.size(); iem < nem; ++iem) {
0334     const auto &em = r.emcalo[iem];
0335     float drbest = drMatchEmHad_;
0336     for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0337       const auto &calo = r.calo[ic];
0338       if (calo.floatEmPt() < ptMinFracMatchEm_ * em.floatPt())
0339         continue;
0340       float dr = floatDR(calo, em);
0341       if (dr < drbest) {
0342         em2calo[iem] = ic;
0343         drbest = dr;
0344       }
0345     }
0346     if (debug_ && em2calo[iem] != -1)
0347       dbgPrintf("PFAlgo3 \t EM    %3d (pt %7.2f) matches to calo %3d (pt %7.2f, empt %7.2f) with dr %.3f\n",
0348                 iem,
0349                 em.floatPt(),
0350                 em2calo[iem],
0351                 em2calo[iem] == -1 ? 0.0 : r.calo[em2calo[iem]].floatPt(),
0352                 em2calo[iem] == -1 ? 0.0 : r.calo[em2calo[iem]].floatEmPt(),
0353                 drbest);
0354   }
0355 }
0356 
0357 void PFAlgo3::sum_tk2em(Region &r,
0358                         const std::vector<int> &tk2em,
0359                         std::vector<int> &em2ntk,
0360                         std::vector<float> &em2sumtkpt,
0361                         std::vector<float> &em2sumtkpterr) const {
0362   // for each EM cluster, count and add up the pt of all the corresponding tracks (skipping muons)
0363   for (int iem = 0, nem = r.emcalo.size(); iem < nem; ++iem) {
0364     const auto &em = r.emcalo[iem];
0365     if (r.globalAbsEta(em.floatEta()) > 2.5)
0366       continue;
0367     for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0368       if (tk2em[itk] == iem) {
0369         const auto &tk = r.track[itk];
0370         if (tk.muonLink)
0371           continue;
0372         em2ntk[iem]++;
0373         em2sumtkpt[iem] += tk.floatPt();
0374         em2sumtkpterr[iem] += tk.floatPtErr();
0375       }
0376     }
0377   }
0378 }
0379 
0380 void PFAlgo3::emcalo_algo(Region &r,
0381                           const std::vector<int> &em2ntk,
0382                           const std::vector<float> &em2sumtkpt,
0383                           const std::vector<float> &em2sumtkpterr) const {
0384   // process ecal clusters after linking
0385   for (int iem = 0, nem = r.emcalo.size(); iem < nem; ++iem) {
0386     auto &em = r.emcalo[iem];
0387     em.isEM = false;
0388     em.used = false;
0389     em.hwFlags = 0;
0390     if (r.globalAbsEta(em.floatEta()) > 2.5)
0391       continue;
0392     if (debug_)
0393       dbgPrintf("PFAlgo3 \t EM    %3d (pt %7.2f) has %2d tracks (sumpt %7.2f, sumpterr %7.2f), ptdif %7.2f +- %7.2f\n",
0394                 iem,
0395                 em.floatPt(),
0396                 em2ntk[iem],
0397                 em2sumtkpt[iem],
0398                 em2sumtkpterr[iem],
0399                 em.floatPt() - em2sumtkpt[iem],
0400                 std::max<float>(em2sumtkpterr[iem], em.floatPtErr()));
0401     if (em2ntk[iem] == 0) {  // Photon
0402       em.isEM = true;
0403       addCaloToPF(r, em);
0404       em.used = true;
0405       if (debug_)
0406         dbgPrintf("PFAlgo3 \t EM    %3d (pt %7.2f)    ---> promoted to photon\n", iem, em.floatPt());
0407       continue;
0408     }
0409     float ptdiff = em.floatPt() - em2sumtkpt[iem];
0410     float pterr = trackEmUseAlsoTrackSigma_ ? std::max<float>(em2sumtkpterr[iem], em.floatPtErr()) : em.floatPtErr();
0411     // avoid "pt = inf +- inf" track to become an electron.
0412     if (pterr > 2 * em.floatPt()) {
0413       pterr = 2 * em.floatPt();
0414       if (debug_)
0415         dbgPrintf("PFAlgo3 \t EM    %3d (pt %7.2f)    ---> clamp pterr ---> new ptdiff %7.2f +- %7.2f\n",
0416                   iem,
0417                   em.floatPt(),
0418                   ptdiff,
0419                   pterr);
0420     }
0421 
0422     if (ptdiff > -ptMatchLow_ * pterr) {
0423       em.isEM = true;
0424       em.used = true;
0425       // convert leftover to a photon if significant
0426       if (ptdiff > +ptMatchHigh_ * pterr) {
0427         auto &p = addCaloToPF(r, em);
0428         p.setFloatPt(ptdiff);
0429         if (debug_)
0430           dbgPrintf("PFAlgo3 \t EM    %3d (pt %7.2f)    ---> promoted to electron(s) + photon (pt %7.2f)\n",
0431                     iem,
0432                     em.floatPt(),
0433                     ptdiff);
0434       } else {
0435         em.hwFlags = 1;  // may use calo momentum
0436         if (debug_)
0437           dbgPrintf("PFAlgo3 \t EM    %3d (pt %7.2f)    ---> promoted to electron(s)\n", iem, em.floatPt());
0438       }
0439     } else {
0440       em.isEM = false;
0441       em.used = false;
0442       em.hwFlags = 0;
0443       //discardCalo(r, em, 2);
0444     }
0445   }
0446 }
0447 
0448 void PFAlgo3::emtk_algo(Region &r,
0449                         const std::vector<int> &tk2em,
0450                         const std::vector<int> &em2ntk,
0451                         const std::vector<float> &em2sumtkpterr) const {
0452   // promote all flagged tracks to electrons
0453   for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0454     auto &tk = r.track[itk];
0455     if (tk2em[itk] == -1 || tk.muonLink)
0456       continue;
0457     const auto &em = r.emcalo[tk2em[itk]];
0458     if (em.isEM) {
0459       auto &p = addTrackToPF(r, tk);
0460       p.cluster.src = em.src;
0461       // FIXME to check if this is useful
0462       if (trackEmMayUseCaloMomenta_ && em2ntk[tk2em[itk]] == 1 && em.hwFlags == 1) {
0463         if (em.floatPtErr() < em2sumtkpterr[tk2em[itk]]) {
0464           p.setFloatPt(em.floatPt());
0465         }
0466       }
0467       if (debug_)
0468         dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) matched to EM   %3d (pt %7.2f) promoted to electron with pt %7.2f\n",
0469                   itk,
0470                   tk.floatPt(),
0471                   tk2em[itk],
0472                   em.floatPt(),
0473                   p.floatPt());
0474       p.hwId = l1t::PFCandidate::Electron;
0475       tk.used = true;
0476     }
0477   }
0478 }
0479 
0480 void PFAlgo3::sub_em2calo(Region &r, const std::vector<int> &em2calo) const {
0481   // subtract EM component from Calo clusters for all photons and electrons (within tracker coverage)
0482   // kill clusters that end up below their own uncertainty, or that loose 90% of the energy,
0483   // unless they still have live EM clusters pointing to them
0484   for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0485     auto &calo = r.calo[ic];
0486     float pt0 = calo.floatPt(), ept0 = calo.floatEmPt(), pt = pt0, ept = ept0;
0487     bool keepme = false;
0488     for (int iem = 0, nem = r.emcalo.size(); iem < nem; ++iem) {
0489       if (em2calo[iem] == ic) {
0490         const auto &em = r.emcalo[iem];
0491         if (em.isEM) {
0492           if (debug_)
0493             dbgPrintf(
0494                 "PFAlgo3 \t EM    %3d (pt %7.2f) is  subtracted from calo %3d (pt %7.2f) scaled by %.3f (deltaPt = "
0495                 "%7.2f)\n",
0496                 iem,
0497                 em.floatPt(),
0498                 ic,
0499                 calo.floatPt(),
0500                 emHadSubtractionPtSlope_,
0501                 emHadSubtractionPtSlope_ * em.floatPt());
0502           pt -= emHadSubtractionPtSlope_ * em.floatPt();
0503           ept -= em.floatPt();
0504         } else {
0505           keepme = true;
0506           if (debug_)
0507             dbgPrintf(
0508                 "PFAlgo3 \t EM    %3d (pt %7.2f) not subtracted from calo %3d (pt %7.2f), and calo marked to be kept "
0509                 "after EM subtraction\n",
0510                 iem,
0511                 em.floatPt(),
0512                 ic,
0513                 calo.floatPt());
0514         }
0515       }
0516     }
0517     if (pt < pt0) {
0518       if (debug_)
0519         dbgPrintf(
0520             "PFAlgo3 \t calo  %3d (pt %7.2f +- %7.2f) has a subtracted pt of %7.2f, empt %7.2f -> %7.2f, isem %d\n",
0521             ic,
0522             calo.floatPt(),
0523             calo.floatPtErr(),
0524             pt,
0525             ept0,
0526             ept,
0527             calo.isEM);
0528       calo.setFloatPt(pt);
0529       calo.setFloatEmPt(ept);
0530       if (!keepme &&
0531           ((emCaloUseAlsoCaloSigma_ ? pt < calo.floatPtErr() : false) || pt <= 0.125 * pt0 ||
0532            (calo.isEM && ept <= 0.125 * ept0))) {  // the <= is important since in firmware the pt0/8 can be zero
0533         if (debug_)
0534           dbgPrintf("PFAlgo3 \t calo  %3d (pt %7.2f)    ----> discarded\n", ic, calo.floatPt());
0535         calo.used = true;
0536         calo.setFloatPt(pt0);  //discardCalo(r, calo, 1);  // log this as discarded, for debugging
0537       }
0538     }
0539   }
0540 }
0541 
0542 void PFAlgo3::link_tk2calo(Region &r, std::vector<int> &tk2calo) const {
0543   // track to calo matching (first iteration, with a lower bound on the calo pt; there may be another one later)
0544   for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0545     const auto &tk = r.track[itk];
0546     if (!tk.quality(l1tpf_impl::InputTrack::PFLOOSE))
0547       continue;
0548     if (tk.muonLink || tk.used)
0549       continue;  // not necessary but just a waste of CPU otherwise
0550     float drbest = drMatch_, dptscale = 0;
0551     switch (tkCaloLinkMetric_) {
0552       case TkCaloLinkMetric::BestByDR:
0553         drbest = drMatch_;
0554         break;
0555       case TkCaloLinkMetric::BestByDRPt:
0556         drbest = 1.0;
0557         dptscale = drMatch_ / tk.floatCaloPtErr();
0558         break;
0559       case TkCaloLinkMetric::BestByDR2Pt2:
0560         drbest = 1.0;
0561         dptscale = drMatch_ / tk.floatCaloPtErr();
0562         break;
0563     }
0564     float minCaloPt = tk.floatPt() - ptMatchLow_ * tk.floatCaloPtErr();
0565     if (debug_)
0566       dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) to be matched to calo, min pT %7.2f\n", itk, tk.floatPt(), minCaloPt);
0567     for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0568       auto &calo = r.calo[ic];
0569       if (calo.used || calo.floatPt() <= minCaloPt)
0570         continue;
0571       float dr = floatDR(tk, calo), dq;
0572       switch (tkCaloLinkMetric_) {
0573         case TkCaloLinkMetric::BestByDR:
0574           if (dr < drbest) {
0575             tk2calo[itk] = ic;
0576             drbest = dr;
0577           }
0578           break;
0579         case TkCaloLinkMetric::BestByDRPt:
0580           dq = dr + std::max<float>(tk.floatPt() - calo.floatPt(), 0.) * dptscale;
0581           //if (debug_ && dr < 0.2) dbgPrintf("PFAlgo3 \t\t\t track %3d (pt %7.2f) vs calo %3d (pt %7.2f): dr %.3f, dq %.3f\n", itk, tk.floatPt(), ic, calo.floatPt(), dr, dq);
0582           if (dr < drMatch_ && dq < drbest) {
0583             tk2calo[itk] = ic;
0584             drbest = dq;
0585           }
0586           break;
0587         case TkCaloLinkMetric::BestByDR2Pt2:
0588           dq = hypot(dr, std::max<float>(tk.floatPt() - calo.floatPt(), 0.) * dptscale);
0589           //if (debug_ && dr < 0.2) dbgPrintf("PFAlgo3 \t\t\t track %3d (pt %7.2f) vs calo %3d (pt %7.2f): dr %.3f, dq %.3f\n", itk, tk.floatPt(), ic, calo.floatPt(), dr, dq);
0590           if (dr < drMatch_ && dq < drbest) {
0591             tk2calo[itk] = ic;
0592             drbest = dq;
0593           }
0594           break;
0595       }
0596     }
0597     if (debug_ && tk2calo[itk] != -1)
0598       dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) matches to calo %3d (pt %7.2f) with dist %.3f\n",
0599                 itk,
0600                 tk.floatPt(),
0601                 tk2calo[itk],
0602                 tk2calo[itk] == -1 ? 0.0 : r.calo[tk2calo[itk]].floatPt(),
0603                 drbest);
0604     // now we re-do this for debugging sake, it may be done for real later
0605     if (debug_ && tk2calo[itk] == -1) {
0606       int ibest = -1;
0607       drbest = 0.3;
0608       for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0609         auto &calo = r.calo[ic];
0610         if (calo.used)
0611           continue;
0612         float dr = floatDR(tk, calo);
0613         if (dr < drbest) {
0614           ibest = ic;
0615           drbest = dr;
0616         }
0617       }
0618       if (ibest != -1)
0619         dbgPrintf(
0620             "PFAlgo3 \t track %3d (pt %7.2f) would match to calo %3d (pt %7.2f) with dr %.3f if the pt min and dr "
0621             "requirement had been relaxed\n",
0622             itk,
0623             tk.floatPt(),
0624             ibest,
0625             r.calo[ibest].floatPt(),
0626             drbest);
0627     }
0628   }
0629 }
0630 
0631 void PFAlgo3::sum_tk2calo(Region &r,
0632                           const std::vector<int> &tk2calo,
0633                           std::vector<int> &calo2ntk,
0634                           std::vector<float> &calo2sumtkpt,
0635                           std::vector<float> &calo2sumtkpterr) const {
0636   // for each calo, compute the sum of the track pt
0637   for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0638     const auto &calo = r.calo[ic];
0639     if (r.globalAbsEta(calo.floatEta()) > 2.5)
0640       continue;
0641     for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0642       if (tk2calo[itk] == ic) {
0643         const auto &tk = r.track[itk];
0644         if (tk.muonLink || tk.used)
0645           continue;
0646         calo2ntk[ic]++;
0647         calo2sumtkpt[ic] += tk.floatPt();
0648         calo2sumtkpterr[ic] += std::pow(tk.floatCaloPtErr(), sumTkCaloErr2_ ? 2 : 1);
0649       }
0650     }
0651     if (sumTkCaloErr2_ && calo2sumtkpterr[ic] > 0)
0652       calo2sumtkpterr[ic] = std::sqrt(calo2sumtkpterr[ic]);
0653   }
0654 }
0655 
0656 void PFAlgo3::unlinkedtk_algo(Region &r, const std::vector<int> &tk2calo) const {
0657   // in the meantime, promote unlinked low pt tracks to hadrons
0658   for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0659     auto &tk = r.track[itk];
0660     if (tk2calo[itk] != -1 || tk.muonLink || tk.used || !tk.quality(l1tpf_impl::InputTrack::PFLOOSE))
0661       continue;
0662     float maxPt = tk.quality(l1tpf_impl::InputTrack::PFTIGHT) ? tightTrackMaxInvisiblePt_ : maxInvisiblePt_;
0663     if (tk.floatPt() < maxPt) {
0664       if (debug_)
0665         dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) not matched to calo, kept as charged hadron\n", itk, tk.floatPt());
0666       auto &p = addTrackToPF(r, tk);
0667       p.hwStatus = GoodTK_NoCalo;
0668       tk.used = true;
0669     } else {
0670       if (debug_)
0671         dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) not matched to calo, dropped\n", itk, tk.floatPt());
0672       //discardTrack(r, tk, BadTK_NoCalo); // log this as discarded, for debugging
0673     }
0674   }
0675 }
0676 
0677 void PFAlgo3::calo_relink(Region &r,
0678                           const std::vector<int> &calo2ntk,
0679                           const std::vector<float> &calo2sumtkpt,
0680                           const std::vector<float> &calo2sumtkpterr) const {
0681   /// OPTIONAL STEP: try to recover split hadron showers (v1.0):
0682   //     take hadrons that are not track matched, close by a hadron which has an excess of track pt vs calo pt
0683   //     add this pt to the calo pt of the other cluster
0684   //     off by default, as it seems to not do much in jets even if it helps remove tails in single-pion events
0685   std::vector<float> addtopt(r.calo.size(), 0);
0686   for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0687     auto &calo = r.calo[ic];
0688     if (calo2ntk[ic] != 0 || calo.used || r.globalAbsEta(calo.floatEta()) > 2.5)
0689       continue;
0690     int i2best = -1;
0691     float drbest = caloReLinkDr_;
0692     for (int ic2 = 0; ic2 < nc; ++ic2) {
0693       const auto &calo2 = r.calo[ic2];
0694       if (calo2ntk[ic2] == 0 || calo2.used || r.globalAbsEta(calo2.floatEta()) > 2.5)
0695         continue;
0696       float dr = floatDR(calo, calo2);
0697       //// uncomment below for more verbose debugging
0698       //if (debug_ && dr < 0.5) dbgPrintf("PFAlgo3 \t calo  %3d (pt %7.2f) with no tracks is at dr %.3f from calo %3d with pt %7.2f (sum tk pt %7.2f), track excess %7.2f +- %7.2f\n", ic, calo.floatPt(), dr, ic2, calo2.floatPt(), calo2sumtkpt[ic2], calo2sumtkpt[ic2] - calo2.floatPt(), useTrackCaloSigma_ ? calo2sumtkpterr[ic2] : calo2.floatPtErr());
0699       if (dr < drbest) {
0700         float ptdiff =
0701             calo2sumtkpt[ic2] - calo2.floatPt() + (useTrackCaloSigma_ ? calo2sumtkpterr[ic2] : calo2.floatPtErr());
0702         if (ptdiff >= caloReLinkThreshold_ * calo.floatPt()) {
0703           i2best = ic2;
0704           drbest = dr;
0705         }
0706       }
0707     }
0708     if (i2best != -1) {
0709       const auto &calo2 = r.calo[i2best];
0710       if (debug_)
0711         dbgPrintf(
0712             "PFAlgo3 \t calo  %3d (pt %7.2f) with no tracks matched within dr %.3f with calo %3d with pt %7.2f (sum tk "
0713             "pt %7.2f), track excess %7.2f +- %7.2f\n",
0714             ic,
0715             calo.floatPt(),
0716             drbest,
0717             i2best,
0718             calo2.floatPt(),
0719             calo2sumtkpt[i2best],
0720             calo2sumtkpt[i2best] - calo2.floatPt(),
0721             useTrackCaloSigma_ ? calo2sumtkpterr[i2best] : calo2.floatPtErr());
0722       calo.used = true;
0723       addtopt[i2best] += calo.floatPt();
0724     }
0725   }
0726   // we do this at the end, so that the above loop is parallelizable
0727   for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0728     if (addtopt[ic]) {
0729       auto &calo = r.calo[ic];
0730       if (debug_)
0731         dbgPrintf("PFAlgo3 \t calo  %3d (pt %7.2f, sum tk pt %7.2f) is increased to pt %7.2f after merging\n",
0732                   ic,
0733                   calo.floatPt(),
0734                   calo2sumtkpt[ic],
0735                   calo.floatPt() + addtopt[ic]);
0736       calo.setFloatPt(calo.floatPt() + addtopt[ic]);
0737     }
0738   }
0739 }
0740 
0741 void PFAlgo3::linkedcalo_algo(Region &r,
0742                               const std::vector<int> &calo2ntk,
0743                               const std::vector<float> &calo2sumtkpt,
0744                               const std::vector<float> &calo2sumtkpterr,
0745                               std::vector<float> &calo2alpha) const {
0746   /// ------------- next step (needs the previous) ----------------
0747   // process matched calo clusters, compare energy to sum track pt
0748   for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0749     auto &calo = r.calo[ic];
0750     if (calo2ntk[ic] == 0 || calo.used)
0751       continue;
0752     float ptdiff = calo.floatPt() - calo2sumtkpt[ic];
0753     float pterr = useTrackCaloSigma_ ? calo2sumtkpterr[ic] : calo.floatPtErr();
0754     if (debug_)
0755       dbgPrintf(
0756           "PFAlgo3 \t calo  %3d (pt %7.2f +- %7.2f, empt %7.2f) has %2d tracks (sumpt %7.2f, sumpterr %7.2f), ptdif "
0757           "%7.2f +- %7.2f\n",
0758           ic,
0759           calo.floatPt(),
0760           calo.floatPtErr(),
0761           calo.floatEmPt(),
0762           calo2ntk[ic],
0763           calo2sumtkpt[ic],
0764           calo2sumtkpterr[ic],
0765           ptdiff,
0766           pterr);
0767     if (ptdiff > +ptMatchHigh_ * pterr) {
0768       if (ecalPriority_) {
0769         if (calo.floatEmPt() > 1) {
0770           float emptdiff = std::min(ptdiff, calo.floatEmPt());
0771           if (debug_)
0772             dbgPrintf(
0773                 "PFAlgo3 \t calo  %3d (pt %7.2f, empt %7.2f)    ---> make photon with pt %7.2f, reduce ptdiff to %7.2f "
0774                 "+- %7.2f\n",
0775                 ic,
0776                 calo.floatPt(),
0777                 calo.floatEmPt(),
0778                 emptdiff,
0779                 ptdiff - emptdiff,
0780                 pterr);
0781           auto &p = addCaloToPF(r, calo);
0782           p.setFloatPt(emptdiff);
0783           p.hwId = l1t::PFCandidate::Photon;
0784           ptdiff -= emptdiff;
0785         }
0786         if (ptdiff > 2) {
0787           if (debug_)
0788             dbgPrintf("PFAlgo3 \t calo  %3d (pt %7.2f, empt %7.2f)    ---> make also neutral hadron with pt %7.2f\n",
0789                       ic,
0790                       calo.floatPt(),
0791                       calo.floatEmPt(),
0792                       ptdiff);
0793           auto &p = addCaloToPF(r, calo);
0794           p.setFloatPt(ptdiff);
0795           p.hwId = l1t::PFCandidate::NeutralHadron;
0796         }
0797       } else {
0798         if (debug_)
0799           dbgPrintf("PFAlgo3 \t calo  %3d (pt %7.2f)    ---> promoted to neutral with pt %7.2f\n",
0800                     ic,
0801                     calo.floatPt(),
0802                     ptdiff);
0803         auto &p = addCaloToPF(r, calo);
0804         p.setFloatPt(ptdiff);
0805         calo.hwFlags = 0;
0806       }
0807     } else if (ptdiff > -ptMatchLow_ * pterr) {
0808       // nothing to do (weighted average happens when we process the tracks)
0809       calo.hwFlags = 1;
0810       if (debug_)
0811         dbgPrintf(
0812             "PFAlgo3 \t calo  %3d (pt %7.2f)    ---> to be deleted, will use tracks instead\n", ic, calo.floatPt());
0813       //discardCalo(r, calo, 0); // log this as discarded, for debugging
0814     } else {
0815       // tracks overshoot, rescale to tracks to calo
0816       calo2alpha[ic] = rescaleTracks_ ? calo.floatPt() / calo2sumtkpt[ic] : 1.0;
0817       calo.hwFlags = 2;
0818       if (debug_ && rescaleTracks_)
0819         dbgPrintf("PFAlgo3 \t calo  %3d (pt %7.2f)    ---> tracks overshoot and will be scaled down by %.4f\n",
0820                   ic,
0821                   calo.floatPt(),
0822                   calo2alpha[ic]);
0823       if (debug_ && !rescaleTracks_)
0824         dbgPrintf("PFAlgo3 \t calo  %3d (pt %7.2f)    ---> tracks overshoot by %.4f\n",
0825                   ic,
0826                   calo.floatPt(),
0827                   calo2sumtkpt[ic] / calo.floatPt());
0828     }
0829     calo.used = true;
0830   }
0831 }
0832 
0833 void PFAlgo3::linkedtk_algo(Region &r,
0834                             const std::vector<int> &tk2calo,
0835                             const std::vector<int> &calo2ntk,
0836                             const std::vector<float> &calo2alpha) const {
0837   // process matched tracks, if necessary rescale or average
0838   for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0839     auto &tk = r.track[itk];
0840     if (tk2calo[itk] == -1 || tk.muonLink || tk.used)
0841       continue;
0842     auto &p = addTrackToPF(r, tk);
0843     tk.used = true;
0844     const auto &calo = r.calo[tk2calo[itk]];
0845     p.cluster.src = calo.src;
0846     if (calo.hwFlags == 1) {
0847       // can do weighted average if there's just one track
0848       if (calo2ntk[tk2calo[itk]] == 1 && caloTrkWeightedAverage_) {
0849         p.hwStatus = GoodTK_Calo_TkPt;
0850         float ptavg = tk.floatPt();
0851         if (tk.floatPtErr() > 0) {
0852           float wcalo = 1.0 / std::pow(tk.floatCaloPtErr(), 2);
0853           float wtk = 1.0 / std::pow(tk.floatPtErr(), 2);
0854           ptavg = (calo.floatPt() * wcalo + tk.floatPt() * wtk) / (wcalo + wtk);
0855           p.hwStatus = GoodTK_Calo_TkCaloPt;
0856         }
0857         p.setFloatPt(ptavg);
0858         if (debug_)
0859           dbgPrintf(
0860               "PFAlgo3 \t track %3d (pt %7.2f +- %7.2f) combined with calo %3d (pt %7.2f +- %7.2f (from tk) yielding "
0861               "candidate of pt %7.2f\n",
0862               itk,
0863               tk.floatPt(),
0864               tk.floatPtErr(),
0865               tk2calo[itk],
0866               calo.floatPt(),
0867               tk.floatCaloPtErr(),
0868               ptavg);
0869       } else {
0870         p.hwStatus = GoodTK_Calo_TkPt;
0871         if (debug_)
0872           dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) linked to calo %3d promoted to charged hadron\n",
0873                     itk,
0874                     tk.floatPt(),
0875                     tk2calo[itk]);
0876       }
0877     } else if (calo.hwFlags == 2) {
0878       // must rescale
0879       p.setFloatPt(tk.floatPt() * calo2alpha[tk2calo[itk]]);
0880       p.hwStatus = GoodTk_Calo_CaloPt;
0881       if (debug_)
0882         dbgPrintf(
0883             "PFAlgo3 \t track %3d (pt %7.2f, stubs %2d chi2 %7.1f) linked to calo %3d promoted to charged hadron with "
0884             "pt %7.2f after maybe rescaling\n",
0885             itk,
0886             tk.floatPt(),
0887             int(tk.hwStubs),
0888             tk.hwChi2 * 0.1f,
0889             tk2calo[itk],
0890             p.floatPt());
0891     }
0892   }
0893 }
0894 
0895 void PFAlgo3::unlinkedcalo_algo(Region &r) const {
0896   // process unmatched calo clusters
0897   for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0898     if (!r.calo[ic].used) {
0899       addCaloToPF(r, r.calo[ic]);
0900       if (debug_)
0901         dbgPrintf("PFAlgo3 \t calo  %3d (pt %7.2f) not linked, promoted to neutral\n", ic, r.calo[ic].floatPt());
0902     }
0903   }
0904 }
0905 
0906 void PFAlgo3::save_muons(Region &r, const std::vector<int> &tk2mu) const {
0907   // finally do muons
0908   for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0909     if (r.track[itk].muonLink) {
0910       auto &p = addTrackToPF(r, r.track[itk]);
0911       p.muonsrc = r.muon[tk2mu[itk]].src;
0912       if (debug_)
0913         dbgPrintf("PFAlgo3 \t track %3d (pt %7.2f) promoted to muon.\n", itk, r.track[itk].floatPt());
0914     }
0915   }
0916 }