Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/PFAlgo2HGC.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 PFAlgo2HGC::PFAlgo2HGC(const edm::ParameterSet &iConfig) : PFAlgoBase(iConfig) {
0020   debug_ = iConfig.getUntrackedParameter<int>("debugPFAlgo2HGC", 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   caloReLinkStep_ = linkcfg.getParameter<bool>("caloReLink");
0051   caloReLinkDr_ = linkcfg.getParameter<double>("caloReLinkDR");
0052   caloReLinkThreshold_ = linkcfg.getParameter<double>("caloReLinkThreshold");
0053   rescaleTracks_ = linkcfg.getParameter<bool>("rescaleTracks");
0054   caloTrkWeightedAverage_ = linkcfg.getParameter<bool>("useCaloTrkWeightedAverage");
0055   sumTkCaloErr2_ = linkcfg.getParameter<bool>("sumTkCaloErr2");
0056   ecalPriority_ = linkcfg.getParameter<bool>("ecalPriority");
0057   tightTrackMaxInvisiblePt_ = linkcfg.getParameter<double>("tightTrackMaxInvisiblePt");
0058   sortInputs_ = iConfig.getParameter<bool>("sortInputs");
0059 }
0060 
0061 void PFAlgo2HGC::runPF(Region &r) const {
0062   initRegion(r, sortInputs_);
0063 
0064   /// ------------- first step (can all go in parallel) ----------------
0065 
0066   if (debug_) {
0067     dbgPrintf(
0068         "PFAlgo2HGC\nPFAlgo2HGC region eta [ %+5.2f , %+5.2f ], phi [ %+5.2f , %+5.2f ], fiducial eta [ %+5.2f , "
0069         "%+5.2f ], phi [ %+5.2f , %+5.2f ]\n",
0070         r.etaMin - r.etaExtra,
0071         r.etaMax + r.etaExtra,
0072         r.phiCenter - r.phiHalfWidth - r.phiExtra,
0073         r.phiCenter + r.phiHalfWidth + r.phiExtra,
0074         r.etaMin,
0075         r.etaMax,
0076         r.phiCenter - r.phiHalfWidth,
0077         r.phiCenter + r.phiHalfWidth);
0078     dbgPrintf(
0079         "PFAlgo2HGC \t N(track) %3lu   N(calo) %3lu   N(mu) %3lu\n", r.track.size(), r.calo.size(), r.muon.size());
0080     for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0081       const auto &tk = r.track[itk];
0082       dbgPrintf(
0083           "PFAlgo2HGC \t track %3d: pt %7.2f +- %5.2f  vtx eta %+5.2f  vtx phi %+5.2f  calo eta %+5.2f  calo phi "
0084           "%+5.2f  fid %1d  calo ptErr %7.2f stubs %2d chi2 %7.1f quality %d\n",
0085           itk,
0086           tk.floatPt(),
0087           tk.floatPtErr(),
0088           tk.floatVtxEta(),
0089           tk.floatVtxPhi(),
0090           tk.floatEta(),
0091           tk.floatPhi(),
0092           int(r.fiducialLocal(tk.floatEta(), tk.floatPhi())),
0093           tk.floatCaloPtErr(),
0094           int(tk.hwStubs),
0095           tk.hwChi2 * 0.1f,
0096           int(tk.hwFlags));
0097     }
0098     for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0099       auto &calo = r.calo[ic];
0100       dbgPrintf(
0101           "PFAlgo2HGC \t calo  %3d: pt %7.2f +- %5.2f  vtx eta %+5.2f  vtx phi %+5.2f  calo eta %+5.2f  calo phi "
0102           "%+5.2f  fid %1d  calo ptErr %7.2f em pt %7.2f isEM %1d \n",
0103           ic,
0104           calo.floatPt(),
0105           calo.floatPtErr(),
0106           calo.floatEta(),
0107           calo.floatPhi(),
0108           calo.floatEta(),
0109           calo.floatPhi(),
0110           int(r.fiducialLocal(calo.floatEta(), calo.floatPhi())),
0111           calo.floatPtErr(),
0112           calo.floatEmPt(),
0113           calo.isEM);
0114     }
0115     for (int im = 0, nm = r.muon.size(); im < nm; ++im) {
0116       auto &mu = r.muon[im];
0117       dbgPrintf(
0118           "PFAlgo2HGC \t muon  %3d: pt %7.2f           vtx eta %+5.2f  vtx phi %+5.2f  calo eta %+5.2f  calo phi "
0119           "%+5.2f  fid %1d\n",
0120           im,
0121           mu.floatPt(),
0122           mu.floatEta(),
0123           mu.floatPhi(),
0124           mu.floatEta(),
0125           mu.floatPhi(),
0126           int(r.fiducialLocal(mu.floatEta(), mu.floatPhi())));
0127     }
0128   }
0129 
0130   std::vector<int> tk2mu(r.track.size(), -1), mu2tk(r.muon.size(), -1);
0131   link_tk2mu(r, tk2mu, mu2tk);
0132 
0133   // track to calo matching (first iteration, with a lower bound on the calo pt; there may be another one later)
0134   std::vector<int> tk2calo(r.track.size(), -1);
0135   link_tk2calo(r, tk2calo);
0136 
0137   /// ------------- next step (needs the previous) ----------------
0138   // for each calo, compute the sum of the track pt
0139   std::vector<int> calo2ntk(r.calo.size(), 0);
0140   std::vector<float> calo2sumtkpt(r.calo.size(), 0);
0141   std::vector<float> calo2sumtkpterr(r.calo.size(), 0);
0142   sum_tk2calo(r, tk2calo, calo2ntk, calo2sumtkpt, calo2sumtkpterr);
0143 
0144   // in the meantime, promote unlinked low pt tracks to hadrons
0145   unlinkedtk_algo(r, tk2calo);
0146 
0147   /// ------------- next step (needs the previous) ----------------
0148   /// OPTIONAL STEP: try to recover split hadron showers (v1.0):
0149   //     off by default, as it seems to not do much in jets even if it helps remove tails in single-pion events
0150   if (caloReLinkStep_)
0151     calo_relink(r, calo2ntk, calo2sumtkpt, calo2sumtkpterr);
0152 
0153   /// ------------- next step (needs the previous) ----------------
0154   // process matched calo clusters, compare energy to sum track pt
0155   std::vector<float> calo2alpha(r.calo.size(), 1);
0156   linkedcalo_algo(r, calo2ntk, calo2sumtkpt, calo2sumtkpterr, calo2alpha);
0157 
0158   /// ------------- next step (needs the previous) ----------------
0159   /// process matched tracks, if necessary rescale or average
0160   linkedtk_algo(r, tk2calo, calo2ntk, calo2alpha);
0161   // process unmatched calo clusters
0162   unlinkedcalo_algo(r);
0163   // finally do muons
0164   save_muons(r, tk2mu);
0165 }
0166 
0167 void PFAlgo2HGC::link_tk2mu(Region &r, std::vector<int> &tk2mu, std::vector<int> &mu2tk) const {
0168   // do a rectangular match for the moment; make a box of the same are as a 0.2 cone
0169   int intDrMuonMatchBox = std::ceil(drMatchMu_ * CaloCluster::ETAPHI_SCALE * std::sqrt(M_PI / 4));
0170   for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0171     tk2mu[itk] = false;
0172   }
0173   for (int imu = 0, nmu = r.muon.size(); imu < nmu; ++imu) {
0174     const auto &mu = r.muon[imu];
0175     if (debug_)
0176       dbgPrintf("PFAlgo2HGC \t muon  %3d (pt %7.2f, eta %+5.2f, phi %+5.2f) \n",
0177                 imu,
0178                 mu.floatPt(),
0179                 mu.floatEta(),
0180                 mu.floatPhi());
0181     float minDistance = 9e9;
0182     switch (muMatchMode_) {
0183       case MuMatchMode::BoxBestByPtRatio:
0184         minDistance = 4.;
0185         break;
0186       case MuMatchMode::DrBestByPtRatio:
0187         minDistance = 4.;
0188         break;
0189       case MuMatchMode::DrBestByPtDiff:
0190         minDistance = 0.5 * mu.floatPt();
0191         break;
0192     }
0193     int imatch = -1;
0194     for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0195       const auto &tk = r.track[itk];
0196       if (!tk.quality(l1tpf_impl::InputTrack::PFLOOSE))
0197         continue;
0198       int deta = std::abs(mu.hwEta - tk.hwEta);
0199       int dphi = std::abs((mu.hwPhi - tk.hwPhi) % CaloCluster::PHI_WRAP);
0200       float dr = floatDR(mu, tk);
0201       float dpt = std::abs(mu.floatPt() - tk.floatPt());
0202       float dptr = (mu.hwPt > tk.hwPt ? mu.floatPt() / tk.floatPt() : tk.floatPt() / mu.floatPt());
0203       bool ok = false;
0204       float distance = 9e9;
0205       switch (muMatchMode_) {
0206         case MuMatchMode::BoxBestByPtRatio:
0207           ok = (deta < intDrMuonMatchBox) && (dphi < intDrMuonMatchBox);
0208           distance = dptr;
0209           break;
0210         case MuMatchMode::DrBestByPtRatio:
0211           ok = (dr < drMatchMu_);
0212           distance = dptr;
0213           break;
0214         case MuMatchMode::DrBestByPtDiff:
0215           ok = (dr < drMatchMu_);
0216           distance = dpt;
0217           break;
0218       }
0219       if (debug_ && dr < 0.4) {
0220         dbgPrintf(
0221             "PFAlgo2HGC \t\t possible match with track %3d (pt %7.2f, caloeta %+5.2f, calophi %+5.2f, dr %.2f, eta "
0222             "%+5.2f, phi %+5.2f, dr %.2f):  angular %1d, distance %.3f (vs %.3f)\n",
0223             itk,
0224             tk.floatPt(),
0225             tk.floatEta(),
0226             tk.floatPhi(),
0227             dr,
0228             tk.floatVtxEta(),
0229             tk.floatVtxPhi(),
0230             deltaR(mu.floatEta(), mu.floatPhi(), tk.floatVtxEta(), tk.floatVtxPhi()),
0231             (ok ? 1 : 0),
0232             distance,
0233             minDistance);
0234       }
0235       if (!ok)
0236         continue;
0237       // FIXME for the moment, we do the floating point matching in pt
0238       if (distance < minDistance) {
0239         minDistance = distance;
0240         imatch = itk;
0241       }
0242     }
0243     if (debug_ && imatch > -1)
0244       dbgPrintf("PFAlgo2HGC \t muon  %3d (pt %7.2f) linked to track %3d (pt %7.2f)\n",
0245                 imu,
0246                 mu.floatPt(),
0247                 imatch,
0248                 r.track[imatch].floatPt());
0249     if (debug_ && imatch == -1)
0250       dbgPrintf("PFAlgo2HGC \t muon  %3d (pt %7.2f) not linked to any track\n", imu, mu.floatPt());
0251     mu2tk[imu] = imatch;
0252     if (imatch > -1) {
0253       tk2mu[imatch] = imu;
0254       r.track[imatch].muonLink = true;
0255     }
0256   }
0257 }
0258 
0259 void PFAlgo2HGC::link_tk2calo(Region &r, std::vector<int> &tk2calo) const {
0260   // track to calo matching (first iteration, with a lower bound on the calo pt; there may be another one later)
0261   for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0262     const auto &tk = r.track[itk];
0263     if (!tk.quality(l1tpf_impl::InputTrack::PFLOOSE))
0264       continue;
0265     if (tk.muonLink || tk.used)
0266       continue;  // not necessary but just a waste of CPU otherwise
0267     float drbest = drMatch_, dptscale = 0;
0268     switch (tkCaloLinkMetric_) {
0269       case TkCaloLinkMetric::BestByDR:
0270         drbest = drMatch_;
0271         break;
0272       case TkCaloLinkMetric::BestByDRPt:
0273         drbest = 1.0;
0274         dptscale = drMatch_ / tk.floatCaloPtErr();
0275         break;
0276       case TkCaloLinkMetric::BestByDR2Pt2:
0277         drbest = 1.0;
0278         dptscale = drMatch_ / tk.floatCaloPtErr();
0279         break;
0280     }
0281     float minCaloPt = tk.floatPt() - ptMatchLow_ * tk.floatCaloPtErr();
0282     if (debug_)
0283       dbgPrintf(
0284           "PFAlgo2HGC \t track %3d (pt %7.2f) to be matched to calo, min pT %7.2f\n", itk, tk.floatPt(), minCaloPt);
0285     for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0286       auto &calo = r.calo[ic];
0287       if (calo.used || calo.floatPt() <= minCaloPt)
0288         continue;
0289       float dr = floatDR(tk, calo), dq;
0290       switch (tkCaloLinkMetric_) {
0291         case TkCaloLinkMetric::BestByDR:
0292           if (dr < drbest) {
0293             tk2calo[itk] = ic;
0294             drbest = dr;
0295           }
0296           break;
0297         case TkCaloLinkMetric::BestByDRPt:
0298           dq = dr + std::max<float>(tk.floatPt() - calo.floatPt(), 0.) * dptscale;
0299           if (debug_ > 2 && dr < 0.3)
0300             dbgPrintf("PFAlgo2HGC \t\t\t track %3d (pt %7.2f) vs calo %3d (pt %7.2f): dr %.3f, dq %.3f\n",
0301                       itk,
0302                       tk.floatPt(),
0303                       ic,
0304                       calo.floatPt(),
0305                       dr,
0306                       dq);
0307           if (dr < drMatch_ && dq < drbest) {
0308             tk2calo[itk] = ic;
0309             drbest = dq;
0310           }
0311           break;
0312         case TkCaloLinkMetric::BestByDR2Pt2:
0313           dq = hypot(dr, std::max<float>(tk.floatPt() - calo.floatPt(), 0.) * dptscale);
0314           if (debug_ > 2 && dr < 0.3)
0315             dbgPrintf("PFAlgo2HGC \t\t\t track %3d (pt %7.2f) vs calo %3d (pt %7.2f): dr %.3f, dq %.3f\n",
0316                       itk,
0317                       tk.floatPt(),
0318                       ic,
0319                       calo.floatPt(),
0320                       dr,
0321                       dq);
0322           if (dr < drMatch_ && dq < drbest) {
0323             tk2calo[itk] = ic;
0324             drbest = dq;
0325           }
0326           break;
0327       }
0328     }
0329     if (debug_ && tk2calo[itk] != -1)
0330       dbgPrintf("PFAlgo2HGC \t track %3d (pt %7.2f) matches to calo %3d (pt %7.2f) with dist %.3f (dr %.3f)\n",
0331                 itk,
0332                 tk.floatPt(),
0333                 tk2calo[itk],
0334                 r.calo[tk2calo[itk]].floatPt(),
0335                 drbest,
0336                 floatDR(tk, r.calo[tk2calo[itk]]));
0337     // now we re-do this for debugging sake, it may be done for real later
0338     if (debug_ && tk2calo[itk] == -1) {
0339       int ibest = -1;
0340       drbest = 0.3;
0341       for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0342         auto &calo = r.calo[ic];
0343         if (calo.used)
0344           continue;
0345         float dr = floatDR(tk, calo);
0346         if (dr < drbest) {
0347           ibest = ic;
0348           drbest = dr;
0349         }
0350       }
0351       if (ibest != -1)
0352         dbgPrintf(
0353             "PFAlgo2HGC \t track %3d (pt %7.2f) would match to calo %3d (pt %7.2f) with dr %.3f if the pt min and dr "
0354             "requirement had been relaxed\n",
0355             itk,
0356             tk.floatPt(),
0357             ibest,
0358             r.calo[ibest].floatPt(),
0359             drbest);
0360     }
0361   }
0362 }
0363 
0364 void PFAlgo2HGC::sum_tk2calo(Region &r,
0365                              const std::vector<int> &tk2calo,
0366                              std::vector<int> &calo2ntk,
0367                              std::vector<float> &calo2sumtkpt,
0368                              std::vector<float> &calo2sumtkpterr) const {
0369   // for each calo, compute the sum of the track pt
0370   for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0371     const auto &calo = r.calo[ic];
0372     if (r.globalAbsEta(calo.floatEta()) > 2.5)
0373       continue;
0374     for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0375       if (tk2calo[itk] == ic) {
0376         const auto &tk = r.track[itk];
0377         if (tk.muonLink || tk.used)
0378           continue;
0379         calo2ntk[ic]++;
0380         calo2sumtkpt[ic] += tk.floatPt();
0381         calo2sumtkpterr[ic] += std::pow(tk.floatCaloPtErr(), sumTkCaloErr2_ ? 2 : 1);
0382       }
0383     }
0384     if (sumTkCaloErr2_ && calo2sumtkpterr[ic] > 0)
0385       calo2sumtkpterr[ic] = std::sqrt(calo2sumtkpterr[ic]);
0386   }
0387 }
0388 
0389 void PFAlgo2HGC::unlinkedtk_algo(Region &r, const std::vector<int> &tk2calo) const {
0390   // in the meantime, promote unlinked low pt tracks to hadrons
0391   for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0392     auto &tk = r.track[itk];
0393     if (tk2calo[itk] != -1 || tk.muonLink || tk.used || !tk.quality(l1tpf_impl::InputTrack::PFLOOSE))
0394       continue;
0395     float maxPt = tk.quality(l1tpf_impl::InputTrack::PFTIGHT) ? tightTrackMaxInvisiblePt_ : maxInvisiblePt_;
0396     if (tk.floatPt() < maxPt) {
0397       if (debug_)
0398         dbgPrintf(
0399             "PFAlgo2HGC \t track %3d (pt %7.2f) not matched to calo, kept as charged hadron\n", itk, tk.floatPt());
0400       auto &p = addTrackToPF(r, tk);
0401       p.hwStatus = GoodTK_NoCalo;
0402       tk.used = true;
0403     } else {
0404       if (debug_)
0405         dbgPrintf("PFAlgo2HGC \t track %3d (pt %7.2f) not matched to calo, dropped\n", itk, tk.floatPt());
0406     }
0407   }
0408 }
0409 
0410 void PFAlgo2HGC::calo_relink(Region &r,
0411                              const std::vector<int> &calo2ntk,
0412                              const std::vector<float> &calo2sumtkpt,
0413                              const std::vector<float> &calo2sumtkpterr) const {
0414   /// OPTIONAL STEP: try to recover split hadron showers (v1.0):
0415   //     take hadrons that are not track matched, close by a hadron which has an excess of track pt vs calo pt
0416   //     add this pt to the calo pt of the other cluster
0417   //     off by default, as it seems to not do much in jets even if it helps remove tails in single-pion events
0418   std::vector<float> addtopt(r.calo.size(), 0);
0419   for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0420     auto &calo = r.calo[ic];
0421     if (calo2ntk[ic] != 0 || calo.used || r.globalAbsEta(calo.floatEta()) > 2.5)
0422       continue;
0423     int i2best = -1;
0424     float drbest = caloReLinkDr_;
0425     for (int ic2 = 0; ic2 < nc; ++ic2) {
0426       const auto &calo2 = r.calo[ic2];
0427       if (calo2ntk[ic2] == 0 || calo2.used || r.globalAbsEta(calo2.floatEta()) > 2.5)
0428         continue;
0429       float dr = floatDR(calo, calo2);
0430       //// uncomment below for more verbose debugging
0431       //if (debug_ && dr < 0.5) dbgPrintf("PFAlgo2HGC \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());
0432       if (dr < drbest) {
0433         float ptdiff =
0434             calo2sumtkpt[ic2] - calo2.floatPt() + (useTrackCaloSigma_ ? calo2sumtkpterr[ic2] : calo2.floatPtErr());
0435         if (ptdiff >= caloReLinkThreshold_ * calo.floatPt()) {
0436           i2best = ic2;
0437           drbest = dr;
0438         }
0439       }
0440     }
0441     if (i2best != -1) {
0442       const auto &calo2 = r.calo[i2best];
0443       if (debug_)
0444         dbgPrintf(
0445             "PFAlgo2HGC \t calo  %3d (pt %7.2f) with no tracks matched within dr %.3f with calo %3d with pt %7.2f (sum "
0446             "tk pt %7.2f), track excess %7.2f +- %7.2f\n",
0447             ic,
0448             calo.floatPt(),
0449             drbest,
0450             i2best,
0451             calo2.floatPt(),
0452             calo2sumtkpt[i2best],
0453             calo2sumtkpt[i2best] - calo2.floatPt(),
0454             useTrackCaloSigma_ ? calo2sumtkpterr[i2best] : calo2.floatPtErr());
0455       calo.used = true;
0456       addtopt[i2best] += calo.floatPt();
0457     }
0458   }
0459   // we do this at the end, so that the above loop is parallelizable
0460   for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0461     if (addtopt[ic]) {
0462       auto &calo = r.calo[ic];
0463       if (debug_)
0464         dbgPrintf("PFAlgo2HGC \t calo  %3d (pt %7.2f, sum tk pt %7.2f) is increased to pt %7.2f after merging\n",
0465                   ic,
0466                   calo.floatPt(),
0467                   calo2sumtkpt[ic],
0468                   calo.floatPt() + addtopt[ic]);
0469       calo.setFloatPt(calo.floatPt() + addtopt[ic]);
0470     }
0471   }
0472 }
0473 
0474 void PFAlgo2HGC::linkedcalo_algo(Region &r,
0475                                  const std::vector<int> &calo2ntk,
0476                                  const std::vector<float> &calo2sumtkpt,
0477                                  const std::vector<float> &calo2sumtkpterr,
0478                                  std::vector<float> &calo2alpha) const {
0479   /// ------------- next step (needs the previous) ----------------
0480   // process matched calo clusters, compare energy to sum track pt
0481   for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0482     auto &calo = r.calo[ic];
0483     if (calo2ntk[ic] == 0 || calo.used)
0484       continue;
0485     float ptdiff = calo.floatPt() - calo2sumtkpt[ic];
0486     float pterr = useTrackCaloSigma_ ? calo2sumtkpterr[ic] : calo.floatPtErr();
0487     if (debug_)
0488       dbgPrintf(
0489           "PFAlgo2HGC \t calo  %3d (pt %7.2f +- %7.2f, empt %7.2f) has %2d tracks (sumpt %7.2f, sumpterr %7.2f), ptdif "
0490           "%7.2f +- %7.2f\n",
0491           ic,
0492           calo.floatPt(),
0493           calo.floatPtErr(),
0494           calo.floatEmPt(),
0495           calo2ntk[ic],
0496           calo2sumtkpt[ic],
0497           calo2sumtkpterr[ic],
0498           ptdiff,
0499           pterr);
0500     if (ptdiff > +ptMatchHigh_ * pterr) {
0501       if (ecalPriority_) {
0502         if (calo.floatEmPt() > 1) {
0503           float emptdiff = std::min(ptdiff, calo.floatEmPt());
0504           if (debug_)
0505             dbgPrintf(
0506                 "PFAlgo2HGC \t calo  %3d (pt %7.2f, empt %7.2f)    ---> make photon with pt %7.2f, reduce ptdiff to "
0507                 "%7.2f +- %7.2f\n",
0508                 ic,
0509                 calo.floatPt(),
0510                 calo.floatEmPt(),
0511                 emptdiff,
0512                 ptdiff - emptdiff,
0513                 pterr);
0514           auto &p = addCaloToPF(r, calo);
0515           p.setFloatPt(emptdiff);
0516           p.hwId = l1t::PFCandidate::Photon;
0517           ptdiff -= emptdiff;
0518         }
0519         if (ptdiff > 2) {
0520           if (debug_)
0521             dbgPrintf("PFAlgo2HGC \t calo  %3d (pt %7.2f, empt %7.2f)    ---> make also neutral hadron with pt %7.2f\n",
0522                       ic,
0523                       calo.floatPt(),
0524                       calo.floatEmPt(),
0525                       ptdiff);
0526           auto &p = addCaloToPF(r, calo);
0527           p.setFloatPt(ptdiff);
0528           p.hwId = l1t::PFCandidate::NeutralHadron;
0529         }
0530       } else {
0531         if (debug_)
0532           dbgPrintf("PFAlgo2HGC \t calo  %3d (pt %7.2f)    ---> promoted to neutral with pt %7.2f\n",
0533                     ic,
0534                     calo.floatPt(),
0535                     ptdiff);
0536         auto &p = addCaloToPF(r, calo);
0537         p.setFloatPt(ptdiff);
0538         calo.hwFlags = 0;
0539       }
0540     } else if (ptdiff > -ptMatchLow_ * pterr) {
0541       // nothing to do (weighted average happens when we process the tracks)
0542       calo.hwFlags = 1;
0543       if (debug_)
0544         dbgPrintf(
0545             "PFAlgo2HGC \t calo  %3d (pt %7.2f)    ---> to be deleted, will use tracks instead\n", ic, calo.floatPt());
0546     } else {
0547       // tracks overshoot, rescale to tracks to calo
0548       calo2alpha[ic] = rescaleTracks_ ? calo.floatPt() / calo2sumtkpt[ic] : 1.0;
0549       calo.hwFlags = 2;
0550       if (debug_ && rescaleTracks_)
0551         dbgPrintf("PFAlgo2HGC \t calo  %3d (pt %7.2f)    ---> tracks overshoot and will be scaled down by %.4f\n",
0552                   ic,
0553                   calo.floatPt(),
0554                   calo2alpha[ic]);
0555       if (debug_ && !rescaleTracks_)
0556         dbgPrintf("PFAlgo2HGC \t calo  %3d (pt %7.2f)    ---> tracks overshoot by %.4f\n",
0557                   ic,
0558                   calo.floatPt(),
0559                   calo2sumtkpt[ic] / calo.floatPt());
0560     }
0561     calo.used = true;
0562   }
0563 }
0564 
0565 void PFAlgo2HGC::linkedtk_algo(Region &r,
0566                                const std::vector<int> &tk2calo,
0567                                const std::vector<int> &calo2ntk,
0568                                const std::vector<float> &calo2alpha) const {
0569   // process matched tracks, if necessary rescale or average
0570   for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0571     auto &tk = r.track[itk];
0572     if (tk2calo[itk] == -1 || tk.muonLink || tk.used)
0573       continue;
0574     auto &p = addTrackToPF(r, tk);
0575     tk.used = true;
0576     const auto &calo = r.calo[tk2calo[itk]];
0577     if (calo.isEM)
0578       p.hwId = l1t::PFCandidate::Electron;
0579     p.cluster.src = calo.src;
0580     if (calo.hwFlags == 1) {
0581       // can do weighted average if there's just one track
0582       if (calo2ntk[tk2calo[itk]] == 1 && caloTrkWeightedAverage_) {
0583         p.hwStatus = GoodTK_Calo_TkPt;
0584         float ptavg = tk.floatPt();
0585         if (tk.floatPtErr() > 0) {
0586           float wcalo = 1.0 / std::pow(tk.floatCaloPtErr(), 2);
0587           float wtk = 1.0 / std::pow(tk.floatPtErr(), 2);
0588           ptavg = (calo.floatPt() * wcalo + tk.floatPt() * wtk) / (wcalo + wtk);
0589           p.hwStatus = GoodTK_Calo_TkCaloPt;
0590         }
0591         p.setFloatPt(ptavg);
0592         if (debug_)
0593           dbgPrintf(
0594               "PFAlgo2HGC \t track %3d (pt %7.2f +- %7.2f) combined with calo %3d (pt %7.2f +- %7.2f (from tk) "
0595               "yielding candidate of pt %7.2f\n",
0596               itk,
0597               tk.floatPt(),
0598               tk.floatPtErr(),
0599               tk2calo[itk],
0600               calo.floatPt(),
0601               tk.floatCaloPtErr(),
0602               ptavg);
0603       } else {
0604         p.hwStatus = GoodTK_Calo_TkPt;
0605         if (debug_)
0606           dbgPrintf("PFAlgo2HGC \t track %3d (pt %7.2f) linked to calo %3d promoted to %s\n",
0607                     itk,
0608                     tk.floatPt(),
0609                     tk2calo[itk],
0610                     (p.hwId == l1t::PFCandidate::Electron ? "electron" : "charged hadron"));
0611       }
0612     } else if (calo.hwFlags == 2) {
0613       // must rescale
0614       p.setFloatPt(tk.floatPt() * calo2alpha[tk2calo[itk]]);
0615       p.hwStatus = GoodTk_Calo_CaloPt;
0616       if (debug_)
0617         dbgPrintf(
0618             "PFAlgo2HGC \t track %3d (pt %7.2f, stubs %2d chi2 %7.1f) linked to calo %3d promoted to %s with pt %7.2f "
0619             "after maybe rescaling\n",
0620             itk,
0621             tk.floatPt(),
0622             int(tk.hwStubs),
0623             tk.hwChi2 * 0.1f,
0624             tk2calo[itk],
0625             (p.hwId == l1t::PFCandidate::Electron ? "electron" : "charged hadron"),
0626             p.floatPt());
0627     }
0628   }
0629 }
0630 
0631 void PFAlgo2HGC::unlinkedcalo_algo(Region &r) const {
0632   // process unmatched calo clusters
0633   for (int ic = 0, nc = r.calo.size(); ic < nc; ++ic) {
0634     if (!r.calo[ic].used) {
0635       addCaloToPF(r, r.calo[ic]);
0636       if (debug_)
0637         dbgPrintf("PFAlgo2HGC \t calo  %3d (pt %7.2f) not linked, promoted to neutral\n", ic, r.calo[ic].floatPt());
0638     }
0639   }
0640 }
0641 
0642 void PFAlgo2HGC::save_muons(Region &r, const std::vector<int> &tk2mu) const {
0643   // finally do muons
0644   for (int itk = 0, ntk = r.track.size(); itk < ntk; ++itk) {
0645     if (r.track[itk].muonLink) {
0646       auto &p = addTrackToPF(r, r.track[itk]);
0647       p.muonsrc = r.muon[tk2mu[itk]].src;
0648       if (debug_)
0649         dbgPrintf("PFAlgo2HGC \t track %3d (pt %7.2f) promoted to muon.\n", itk, r.track[itk].floatPt());
0650     }
0651   }
0652 }