Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/pf/pfalgo3_ref.h"
0002 #include "L1Trigger/Phase2L1ParticleFlow/interface/dbgPrintf.h"
0003 
0004 #include <cmath>
0005 #include <cstdio>
0006 #include <algorithm>
0007 #include <vector>
0008 #include <memory>
0009 
0010 #ifdef CMSSW_GIT_HASH
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 
0013 l1ct::PFAlgo3Emulator::PFAlgo3Emulator(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       nEMCALO_(iConfig.getParameter<uint32_t>("nEmCalo")),
0023       nPHOTON_(iConfig.getParameter<uint32_t>("nPhoton")),
0024       nALLNEUTRAL_(iConfig.getParameter<uint32_t>("nAllNeutral")),
0025       dR2MAX_TK_EM_(l1ct::Scales::makeDR2FromFloatDR(iConfig.getParameter<double>("trackEmDR"))),
0026       dR2MAX_EM_CALO_(l1ct::Scales::makeDR2FromFloatDR(iConfig.getParameter<double>("emCaloDR"))) {
0027   debug_ = iConfig.getUntrackedParameter<bool>("debug", false);
0028   loadPtErrBins(iConfig);
0029 }
0030 #endif
0031 
0032 void l1ct::PFAlgo3Emulator::toFirmware(const PFInputRegion& in,
0033                                        PFRegion& region,
0034                                        HadCaloObj calo[/*nCALO*/],
0035                                        EmCaloObj emcalo[/*nEMCALO*/],
0036                                        TkObj track[/*nTRACK*/],
0037                                        MuObj mu[/*nMU*/]) const {
0038   region = in.region;
0039   l1ct::toFirmware(in.track, nTRACK_, track);
0040   l1ct::toFirmware(in.emcalo, nEMCALO_, emcalo);
0041   l1ct::toFirmware(in.hadcalo, nCALO_, calo);
0042   l1ct::toFirmware(in.muon, nMU_, mu);
0043 }
0044 
0045 void l1ct::PFAlgo3Emulator::toFirmware(const OutputRegion& out,
0046                                        PFChargedObj outch[/*nTRACK*/],
0047                                        PFNeutralObj outpho[/*nPHOTON*/],
0048                                        PFNeutralObj outne[/*nSELCALO*/],
0049                                        PFChargedObj outmu[/*nMU*/]) const {
0050   l1ct::toFirmware(out.pfcharged, nTRACK_, outch);
0051   l1ct::toFirmware(out.pfphoton, nPHOTON_, outpho);
0052   l1ct::toFirmware(out.pfneutral, nSELCALO_, outne);
0053   l1ct::toFirmware(out.pfmuon, nMU_, outmu);
0054 }
0055 
0056 int l1ct::PFAlgo3Emulator::tk_best_match_ref(unsigned int dR2MAX,
0057                                              const std::vector<l1ct::EmCaloObjEmu>& calo,
0058                                              const l1ct::TkObjEmu& track) const {
0059   int drmin = dR2MAX, ibest = -1;
0060   for (unsigned int ic = 0, nCAL = std::min<unsigned>(nEMCALO_, calo.size()); ic < nCAL; ++ic) {
0061     if (calo[ic].hwPt <= 0)
0062       continue;
0063     int dr = dr2_int(track.hwEta, track.hwPhi, calo[ic].hwEta, calo[ic].hwPhi);
0064     if (dr < drmin) {
0065       drmin = dr;
0066       ibest = ic;
0067     }
0068   }
0069   return ibest;
0070 }
0071 int l1ct::PFAlgo3Emulator::em_best_match_ref(unsigned int dR2MAX,
0072                                              const std::vector<l1ct::HadCaloObjEmu>& calo,
0073                                              const l1ct::EmCaloObjEmu& em) const {
0074   pt_t emPtMin = em.hwPt >> 1;
0075   int drmin = dR2MAX, ibest = -1;
0076   for (unsigned int ic = 0, nCAL = calo.size(); ic < nCAL; ++ic) {
0077     if (calo[ic].hwEmPt <= emPtMin)
0078       continue;
0079     int dr = dr2_int(em.hwEta, em.hwPhi, calo[ic].hwEta, calo[ic].hwPhi);
0080     if (dr < drmin) {
0081       drmin = dr;
0082       ibest = ic;
0083     }
0084   }
0085   return ibest;
0086 }
0087 
0088 void l1ct::PFAlgo3Emulator::pfalgo3_em_ref(const PFInputRegion& in,
0089                                            const std::vector<int>& iMu /*[nTRACK]*/,
0090                                            std::vector<int>& iEle /*[nTRACK]*/,
0091                                            OutputRegion& out,
0092                                            std::vector<HadCaloObjEmu>& hadcalo_out /*[nCALO]*/) const {
0093   // constants
0094   unsigned int nTRACK = std::min<unsigned>(nTRACK_, in.track.size());
0095   unsigned int nEMCALO = std::min<unsigned>(nEMCALO_, in.emcalo.size());
0096   unsigned int nPHOTON = std::min<unsigned>(nPHOTON_, in.emcalo.size());
0097   unsigned int nCALO = std::min<unsigned>(nCALO_, in.hadcalo.size());
0098 
0099   // initialize sum track pt
0100   std::vector<pt_t> calo_sumtk(nEMCALO);
0101   for (unsigned int ic = 0; ic < nEMCALO; ++ic) {
0102     calo_sumtk[ic] = 0;
0103   }
0104   std::vector<int> tk2em(nTRACK);
0105   std::vector<bool> isEM(nEMCALO);
0106   // for each track, find the closest calo
0107   for (unsigned int it = 0; it < nTRACK; ++it) {
0108     if (in.track[it].hwPt > 0 && in.track[it].isPFLoose() && iMu[it] == -1) {
0109       tk2em[it] = tk_best_match_ref(dR2MAX_TK_EM_, in.emcalo, in.track[it]);
0110       if (tk2em[it] != -1) {
0111         if (debug_)
0112           dbgPrintf(
0113               "FW  \t track  %3d pt %8.2f matched to em calo %3d pt %8.2f (int deltaR2 %d)\n",
0114               it,
0115               in.track[it].floatPt(),
0116               tk2em[it],
0117               in.emcalo[tk2em[it]].floatPt(),
0118               dr2_int(in.track[it].hwEta, in.track[it].hwPhi, in.emcalo[tk2em[it]].hwEta, in.emcalo[tk2em[it]].hwPhi));
0119         calo_sumtk[tk2em[it]] += in.track[it].hwPt;
0120       }
0121     } else {
0122       tk2em[it] = -1;
0123     }
0124   }
0125 
0126   if (debug_) {
0127     for (unsigned int ic = 0; ic < nEMCALO; ++ic) {
0128       if (in.emcalo[ic].hwPt > 0)
0129         dbgPrintf("FW  \t emcalo %3d pt %8.2f has sumtk %8.2f\n",
0130                   ic,
0131                   in.emcalo[ic].floatPt(),
0132                   Scales::floatPt(calo_sumtk[ic]));
0133     }
0134   }
0135 
0136   assert(nEMCALO == nPHOTON);  // code doesn't work otherwise
0137   out.pfphoton.resize(nPHOTON);
0138   for (unsigned int ic = 0; ic < nEMCALO; ++ic) {
0139     pt_t photonPt;
0140     if (calo_sumtk[ic] > 0) {
0141       dpt_t ptdiff = dpt_t(in.emcalo[ic].hwPt) - dpt_t(calo_sumtk[ic]);
0142       pt2_t sigma2 = in.emcalo[ic].hwPtErr * in.emcalo[ic].hwPtErr;
0143       pt2_t sigma2Lo = 4 * sigma2;
0144       const pt2_t& sigma2Hi = sigma2;  // + (sigma2>>1); // cut at 1 sigma instead of old cut at sqrt(1.5) sigma's
0145       pt2_t ptdiff2 = ptdiff * ptdiff;
0146       if ((ptdiff >= 0 && ptdiff2 <= sigma2Hi) || (ptdiff < 0 && ptdiff2 < sigma2Lo)) {
0147         // electron
0148         photonPt = 0;
0149         isEM[ic] = true;
0150         if (debug_)
0151           dbgPrintf("FW  \t emcalo %3d pt %8.2f ptdiff %8.2f [match window: -%.2f / +%.2f] flagged as electron\n",
0152                     ic,
0153                     in.emcalo[ic].floatPt(),
0154                     Scales::floatPt(ptdiff),
0155                     std::sqrt(Scales::floatPt(sigma2Lo)),
0156                     std::sqrt(float(sigma2Hi)));
0157       } else if (ptdiff > 0) {
0158         // electron + photon
0159         photonPt = ptdiff;
0160         isEM[ic] = true;
0161         if (debug_)
0162           dbgPrintf(
0163               "FW  \t emcalo %3d pt %8.2f ptdiff %8.2f [match window: -%.2f / +%.2f] flagged as electron + photon of "
0164               "pt %8.2f\n",
0165               ic,
0166               in.emcalo[ic].floatPt(),
0167               Scales::floatPt(ptdiff),
0168               std::sqrt(Scales::floatPt(sigma2Lo)),
0169               std::sqrt(float(sigma2Hi)),
0170               Scales::floatPt(photonPt));
0171       } else {
0172         // pion
0173         photonPt = 0;
0174         isEM[ic] = false;
0175         if (debug_)
0176           dbgPrintf("FW  \t emcalo %3d pt %8.2f ptdiff %8.2f [match window: -%.2f / +%.2f] flagged as pion\n",
0177                     ic,
0178                     in.emcalo[ic].floatPt(),
0179                     Scales::floatPt(ptdiff),
0180                     std::sqrt(Scales::floatPt(sigma2Lo)),
0181                     std::sqrt(Scales::floatPt(sigma2Hi)));
0182       }
0183     } else {
0184       // photon
0185       isEM[ic] = true;
0186       photonPt = in.emcalo[ic].hwPt;
0187       if (debug_ && in.emcalo[ic].hwPt > 0)
0188         dbgPrintf("FW  \t emcalo %3d pt %8.2f flagged as photon\n", ic, in.emcalo[ic].floatPt());
0189     }
0190     if (photonPt) {
0191       fillPFCand(in.emcalo[ic], out.pfphoton[ic]);
0192       out.pfphoton[ic].hwPt = photonPt;
0193       out.pfphoton[ic].hwEmPt = photonPt;
0194     } else {
0195       out.pfphoton[ic].clear();
0196     }
0197   }
0198 
0199   iEle.resize(nTRACK);
0200   for (unsigned int it = 0; it < nTRACK; ++it) {
0201     iEle[it] = ((tk2em[it] != -1) && isEM[tk2em[it]]) ? tk2em[it] : -1;
0202     if (debug_ && (iEle[it] != -1))
0203       dbgPrintf(
0204           "FW  \t track  %3d pt %8.2f flagged as electron (emcluster %d).\n", it, in.track[it].floatPt(), iEle[it]);
0205   }
0206 
0207   std::vector<int> em2calo(nEMCALO);
0208   for (unsigned int ic = 0; ic < nEMCALO; ++ic) {
0209     em2calo[ic] = em_best_match_ref(dR2MAX_EM_CALO_, in.hadcalo, in.emcalo[ic]);
0210     if (debug_ && (in.emcalo[ic].hwPt > 0)) {
0211       dbgPrintf("FW  \t emcalo %3d pt %8.2f isEM %d matched to hadcalo %3d pt %8.2f emPt %8.2f isEM %d\n",
0212                 ic,
0213                 in.emcalo[ic].floatPt(),
0214                 int(isEM[ic]),
0215                 em2calo[ic],
0216                 (em2calo[ic] >= 0 ? in.hadcalo[em2calo[ic]].floatPt() : -1),
0217                 (em2calo[ic] >= 0 ? in.hadcalo[em2calo[ic]].floatEmPt() : -1),
0218                 (em2calo[ic] >= 0 ? int(in.hadcalo[em2calo[ic]].hwIsEM()) : 0));
0219     }
0220   }
0221 
0222   hadcalo_out.resize(nCALO);
0223   for (unsigned int ih = 0; ih < nCALO; ++ih) {
0224     hadcalo_out[ih] = in.hadcalo[ih];
0225     dpt_t sub = 0;
0226     bool keep = false;
0227     for (unsigned int ic = 0; ic < nEMCALO; ++ic) {
0228       if (em2calo[ic] == int(ih)) {
0229         if (isEM[ic])
0230           sub += in.emcalo[ic].hwPt;
0231         else
0232           keep = true;
0233       }
0234     }
0235     dpt_t emdiff = dpt_t(in.hadcalo[ih].hwEmPt) - sub;  // ok to saturate at zero here
0236     dpt_t alldiff = dpt_t(in.hadcalo[ih].hwPt) - sub;
0237     if (debug_ && (in.hadcalo[ih].hwPt > 0)) {
0238       dbgPrintf("FW  \t calo   %3d pt %8.2f has a subtracted pt of %8.2f, empt %8.2f -> %8.2f   isem %d mustkeep %d \n",
0239                 ih,
0240                 in.hadcalo[ih].floatPt(),
0241                 Scales::floatPt(alldiff),
0242                 in.hadcalo[ih].floatEmPt(),
0243                 Scales::floatPt(emdiff),
0244                 int(in.hadcalo[ih].hwIsEM()),
0245                 keep);
0246     }
0247     if (alldiff <= (in.hadcalo[ih].hwPt >> 4)) {
0248       hadcalo_out[ih].hwPt = 0;    // kill
0249       hadcalo_out[ih].hwEmPt = 0;  // kill
0250       if (debug_ && (in.hadcalo[ih].hwPt > 0))
0251         dbgPrintf("FW  \t calo   %3d pt %8.2f --> discarded (zero pt)\n", ih, in.hadcalo[ih].floatPt());
0252     } else if ((in.hadcalo[ih].hwIsEM() && emdiff <= (in.hadcalo[ih].hwEmPt >> 3)) && !keep) {
0253       hadcalo_out[ih].hwPt = 0;    // kill
0254       hadcalo_out[ih].hwEmPt = 0;  // kill
0255       if (debug_ && (in.hadcalo[ih].hwPt > 0))
0256         dbgPrintf("FW  \t calo   %3d pt %8.2f --> discarded (zero em)\n", ih, in.hadcalo[ih].floatPt());
0257     } else {
0258       hadcalo_out[ih].hwPt = alldiff;
0259       hadcalo_out[ih].hwEmPt = (emdiff > 0 ? pt_t(emdiff) : pt_t(0));
0260     }
0261   }
0262 }
0263 
0264 void l1ct::PFAlgo3Emulator::run(const PFInputRegion& in, OutputRegion& out) const {
0265   // constants
0266   unsigned int nTRACK = std::min<unsigned>(nTRACK_, in.track.size());
0267   unsigned int nEMCALO = std::min<unsigned>(nEMCALO_, in.emcalo.size());
0268   unsigned int nPHOTON = std::min<unsigned>(nPHOTON_, in.emcalo.size());
0269   unsigned int nCALO = std::min<unsigned>(nCALO_, in.hadcalo.size());
0270   unsigned int nSELCALO = std::min<unsigned>(nSELCALO_, in.hadcalo.size());
0271   unsigned int nMU = std::min<unsigned>(nMU_, in.muon.size());
0272 
0273   if (debug_) {
0274     dbgPrintf("FW\nFW  \t region eta %+5.2f [ %+5.2f , %+5.2f ], phi %+5.2f [ %+5.2f , %+5.2f ]   packed %s\n",
0275               in.region.floatEtaCenter(),
0276               in.region.floatEtaMinExtra(),
0277               in.region.floatEtaMaxExtra(),
0278               in.region.floatPhiCenter(),
0279               in.region.floatPhiCenter() - in.region.floatPhiHalfWidthExtra(),
0280               in.region.floatPhiCenter() + in.region.floatPhiHalfWidthExtra(),
0281               in.region.pack().to_string(16).c_str());
0282 
0283     dbgPrintf("FW  \t N(track) %3lu   N(em) %3lu   N(calo) %3lu   N(mu) %3lu\n",
0284               in.track.size(),
0285               in.emcalo.size(),
0286               in.hadcalo.size(),
0287               in.muon.size());
0288 
0289     for (unsigned int i = 0; i < nTRACK; ++i) {
0290       if (in.track[i].hwPt == 0)
0291         continue;
0292       dbgPrintf(
0293           "FW  \t track %3d: pt %8.2f [ %8d ]  calo eta %+5.2f [ %+5d ]  calo phi %+5.2f [ %+5d ]  vtx eta %+5.2f  "
0294           "vtx phi %+5.2f   charge %+2d  qual %2d  fid %d  glb eta %+5.2f phi %+5.2f  packed %s\n",
0295           i,
0296           in.track[i].floatPt(),
0297           in.track[i].intPt(),
0298           in.track[i].floatEta(),
0299           in.track[i].intEta(),
0300           in.track[i].floatPhi(),
0301           in.track[i].intPhi(),
0302           in.track[i].floatVtxEta(),
0303           in.track[i].floatVtxPhi(),
0304           in.track[i].intCharge(),
0305           int(in.track[i].hwQuality),
0306           int(in.region.isFiducial(in.track[i].hwEta, in.track[i].hwPhi)),
0307           in.region.floatGlbEta(in.track[i].hwVtxEta()),
0308           in.region.floatGlbPhi(in.track[i].hwVtxPhi()),
0309           in.track[i].pack().to_string(16).c_str());
0310     }
0311     for (unsigned int i = 0; i < nEMCALO; ++i) {
0312       if (in.emcalo[i].hwPt == 0)
0313         continue;
0314       dbgPrintf(
0315           "FW  \t EM    %3d: pt %8.2f [ %8d ]  calo eta %+5.2f [ %+5d ]  calo phi %+5.2f [ %+5d ]  calo ptErr %8.2f [ "
0316           "%6d ]  emID %2d  fid %d  glb eta %+5.2f phi %+5.2f  packed %s \n",
0317           i,
0318           in.emcalo[i].floatPt(),
0319           in.emcalo[i].intPt(),
0320           in.emcalo[i].floatEta(),
0321           in.emcalo[i].intEta(),
0322           in.emcalo[i].floatPhi(),
0323           in.emcalo[i].intPhi(),
0324           in.emcalo[i].floatPtErr(),
0325           in.emcalo[i].intPtErr(),
0326           in.emcalo[i].hwEmID.to_int(),
0327           int(in.region.isFiducial(in.emcalo[i].hwEta, in.emcalo[i].hwPhi)),
0328           in.region.floatGlbEtaOf(in.emcalo[i]),
0329           in.region.floatGlbPhiOf(in.emcalo[i]),
0330           in.emcalo[i].pack().to_string(16).c_str());
0331     }
0332     for (unsigned int i = 0; i < nCALO; ++i) {
0333       if (in.hadcalo[i].hwPt == 0)
0334         continue;
0335       dbgPrintf(
0336           "FW  \t calo  %3d: pt %8.2f [ %8d ]  calo eta %+5.2f [ %+5d ]  calo phi %+5.2f [ %+5d ]  calo emPt  %8.2f [ "
0337           "%6d ]   emID %2d  fid %d  glb eta %+5.2f phi %+5.2f  packed %s \n",
0338           i,
0339           in.hadcalo[i].floatPt(),
0340           in.hadcalo[i].intPt(),
0341           in.hadcalo[i].floatEta(),
0342           in.hadcalo[i].intEta(),
0343           in.hadcalo[i].floatPhi(),
0344           in.hadcalo[i].intPhi(),
0345           in.hadcalo[i].floatEmPt(),
0346           in.hadcalo[i].intEmPt(),
0347           in.hadcalo[i].hwEmID.to_int(),
0348           int(in.region.isFiducial(in.hadcalo[i].hwEta, in.hadcalo[i].hwPhi)),
0349           in.region.floatGlbEtaOf(in.hadcalo[i]),
0350           in.region.floatGlbPhiOf(in.hadcalo[i]),
0351           in.hadcalo[i].pack().to_string(16).c_str());
0352     }
0353     for (unsigned int i = 0; i < nMU; ++i) {
0354       if (in.muon[i].hwPt == 0)
0355         continue;
0356       dbgPrintf(
0357           "FW  \t muon  %3d: pt %8.2f [ %8d ]  calo eta %+5.2f [ %+5d ]  calo phi %+5.2f [ %+5d ]  "
0358           "vtx eta %+5.2f  vtx phi %+5.2f  charge %+2d  qual %2d  glb eta %+5.2f phi %+5.2f  packed %s \n",
0359           i,
0360           in.muon[i].floatPt(),
0361           in.muon[i].intPt(),
0362           in.muon[i].floatEta(),
0363           in.muon[i].intEta(),
0364           in.muon[i].floatPhi(),
0365           in.muon[i].intPhi(),
0366           in.muon[i].floatVtxEta(),
0367           in.muon[i].floatVtxPhi(),
0368           in.muon[i].intCharge(),
0369           int(in.muon[i].hwQuality),
0370           in.region.floatGlbEta(in.muon[i].hwVtxEta()),
0371           in.region.floatGlbPhi(in.muon[i].hwVtxPhi()),
0372           in.muon[i].pack().to_string(16).c_str());
0373     }
0374     dbgPrintf("%s", "FW\n");
0375   }
0376 
0377   ////////////////////////////////////////////////////
0378   // TK-MU Linking
0379   std::vector<int> iMu;
0380   pfalgo_mu_ref(in, out, iMu);
0381 
0382   ////////////////////////////////////////////////////
0383   // TK-EM Linking
0384   std::vector<int> iEle;
0385   std::vector<HadCaloObjEmu> hadcalo_subem(nCALO);
0386   pfalgo3_em_ref(in, iMu, iEle, out, hadcalo_subem);
0387 
0388   ////////////////////////////////////////////////////
0389   // TK-HAD Linking
0390 
0391   // initialize sum track pt
0392   std::vector<pt_t> calo_sumtk(nCALO), calo_subpt(nCALO);
0393   std::vector<pt2_t> calo_sumtkErr2(nCALO);
0394   for (unsigned int ic = 0; ic < nCALO; ++ic) {
0395     calo_sumtk[ic] = 0;
0396     calo_sumtkErr2[ic] = 0;
0397   }
0398 
0399   // initialize good track bit
0400   std::vector<bool> track_good(nTRACK, false);
0401   for (unsigned int it = 0; it < nTRACK; ++it) {
0402     if (!in.track[it].isPFLoose())
0403       continue;
0404     pt_t ptInv = in.track[it].isPFTight() ? tk_MAXINVPT_TIGHT_ : tk_MAXINVPT_LOOSE_;
0405     track_good[it] = (in.track[it].hwPt < ptInv) || (iEle[it] != -1) || (iMu[it] != -1);
0406   }
0407 
0408   // initialize output
0409   out.pfcharged.resize(nTRACK);
0410   out.pfneutral.resize(nSELCALO);
0411   for (unsigned int ipf = 0; ipf < nTRACK; ++ipf)
0412     out.pfcharged[ipf].clear();
0413   for (unsigned int ipf = 0; ipf < nSELCALO; ++ipf)
0414     out.pfneutral[ipf].clear();
0415 
0416   // for each track, find the closest calo
0417   std::vector<int> tk2calo(nTRACK, -1);
0418   for (unsigned int it = 0; it < nTRACK; ++it) {
0419     if (in.track[it].hwPt > 0 && in.track[it].isPFLoose() && (iEle[it] == -1) && (iMu[it] == -1)) {
0420       pt_t tkCaloPtErr = ptErr_ref(in.region, in.track[it]);
0421       int ibest = best_match_with_pt_ref(dR2MAX_TK_CALO_, hadcalo_subem, in.track[it], tkCaloPtErr);
0422       if (ibest != -1) {
0423         if (debug_)
0424           dbgPrintf(
0425               "FW  \t track  %3d pt %8.2f matched to calo %3d pt %8.2f (int deltaR2 %d)\n",
0426               it,
0427               in.track[it].floatPt(),
0428               ibest,
0429               hadcalo_subem[ibest].floatPt(),
0430               dr2_int(in.track[it].hwEta, in.track[it].hwPhi, hadcalo_subem[ibest].hwEta, hadcalo_subem[ibest].hwPhi));
0431         track_good[it] = true;
0432         calo_sumtk[ibest] += in.track[it].hwPt;
0433         calo_sumtkErr2[ibest] += tkCaloPtErr * tkCaloPtErr;
0434       }
0435       tk2calo[it] = ibest;  // for emulator info
0436     }
0437   }
0438 
0439   for (unsigned int ic = 0; ic < nCALO; ++ic) {
0440     if (calo_sumtk[ic] > 0) {
0441       dpt_t ptdiff = dpt_t(hadcalo_subem[ic].hwPt) - dpt_t(calo_sumtk[ic]);
0442       pt2_t sigmamult = calo_sumtkErr2
0443           [ic];  // before we did (calo_sumtkErr2[ic] + (calo_sumtkErr2[ic] >> 1)); to multiply by 1.5 = sqrt(1.5)^2 ~ (1.2)^2
0444       if (debug_ && (hadcalo_subem[ic].hwPt > 0)) {
0445         dbgPrintf(
0446             "FW  \t calo  %3d pt %8.2f [ %7d ] eta %+5.2f [ %+5d ] has a sum track pt %8.2f, difference %7.2f +- %.2f "
0447             "\n",
0448             ic,
0449             hadcalo_subem[ic].floatPt(),
0450             hadcalo_subem[ic].intPt(),
0451             hadcalo_subem[ic].floatEta(),
0452             hadcalo_subem[ic].intEta(),
0453             Scales::floatPt(calo_sumtk[ic]),
0454             Scales::floatPt(ptdiff),
0455             std::sqrt(Scales::floatPt(calo_sumtkErr2[ic])));
0456       }
0457       if (ptdiff > 0 && ptdiff * ptdiff > sigmamult) {
0458         calo_subpt[ic] = pt_t(ptdiff);
0459       } else {
0460         calo_subpt[ic] = 0;
0461       }
0462     } else {
0463       calo_subpt[ic] = hadcalo_subem[ic].hwPt;
0464     }
0465     if (debug_ && (hadcalo_subem[ic].hwPt > 0))
0466       dbgPrintf(
0467           "FW  \t calo  %3d pt %8.2f ---> %8.2f \n", ic, hadcalo_subem[ic].floatPt(), Scales::floatPt(calo_subpt[ic]));
0468   }
0469 
0470   // copy out charged hadrons
0471   for (unsigned int it = 0; it < nTRACK; ++it) {
0472     if (track_good[it]) {
0473       fillPFCand(in.track[it], out.pfcharged[it], iMu[it] != -1, iEle[it] != -1);
0474       // extra emulator information
0475       if (iEle[it] != -1)
0476         out.pfcharged[it].srcCluster = in.emcalo[iEle[it]].src;
0477       if (iMu[it] != -1)
0478         out.pfcharged[it].srcMu = in.muon[iMu[it]].src;
0479     }
0480   }
0481 
0482   // copy out neutral hadrons
0483   std::vector<PFNeutralObjEmu> outne_all(nCALO);
0484   for (unsigned int ipf = 0; ipf < nCALO; ++ipf)
0485     outne_all[ipf].clear();
0486   for (unsigned int ic = 0; ic < nCALO; ++ic) {
0487     if (calo_subpt[ic] > 0) {
0488       fillPFCand(hadcalo_subem[ic], outne_all[ic]);
0489       outne_all[ic].hwPt = calo_subpt[ic];
0490       outne_all[ic].hwEmPt = hadcalo_subem[ic].hwIsEM() ? calo_subpt[ic] : pt_t(0);  // FIXME
0491     }
0492   }
0493 
0494   if (nCALO_ == nSELCALO_) {
0495     std::swap(outne_all, out.pfneutral);
0496   } else {
0497     ptsort_ref(nCALO, nSELCALO, outne_all, out.pfneutral);
0498   }
0499 
0500   if (debug_) {
0501     dbgPrintf("%s", "FW\n");
0502     for (unsigned int i = 0; i < nTRACK; ++i) {
0503       if (out.pfcharged[i].hwPt == 0)
0504         continue;
0505       dbgPrintf(
0506           "FW  \t outch %3d: pt %8.2f [ %8d ]  calo eta %+5.2f [ %+5d ]  calo phi %+5.2f [ %+5d ]  pid %d  packed %s\n",
0507           i,
0508           out.pfcharged[i].floatPt(),
0509           out.pfcharged[i].intPt(),
0510           out.pfcharged[i].floatEta(),
0511           out.pfcharged[i].intEta(),
0512           out.pfcharged[i].floatPhi(),
0513           out.pfcharged[i].intPhi(),
0514           out.pfcharged[i].intId(),
0515           out.pfcharged[i].pack().to_string(16).c_str());
0516     }
0517     for (unsigned int i = 0; i < nPHOTON; ++i) {
0518       if (out.pfphoton[i].hwPt == 0)
0519         continue;
0520       dbgPrintf(
0521           "FW  \t outph %3d: pt %8.2f [ %8d ]  calo eta %+5.2f [ %+5d ]  calo phi %+5.2f [ %+5d ]  pid %d  packed %s\n",
0522           i,
0523           out.pfphoton[i].floatPt(),
0524           out.pfphoton[i].intPt(),
0525           out.pfphoton[i].floatEta(),
0526           out.pfphoton[i].intEta(),
0527           out.pfphoton[i].floatPhi(),
0528           out.pfphoton[i].intPhi(),
0529           out.pfphoton[i].intId(),
0530           out.pfphoton[i].pack().to_string(16).c_str());
0531     }
0532     for (unsigned int i = 0; i < nSELCALO; ++i) {
0533       if (out.pfneutral[i].hwPt == 0)
0534         continue;
0535       dbgPrintf(
0536           "FW  \t outne %3d: pt %8.2f [ %8d ]  calo eta %+5.2f [ %+5d ]  calo phi %+5.2f [ %+5d ]  pid %d  packed %s\n",
0537           i,
0538           out.pfneutral[i].floatPt(),
0539           out.pfneutral[i].intPt(),
0540           out.pfneutral[i].floatEta(),
0541           out.pfneutral[i].intEta(),
0542           out.pfneutral[i].floatPhi(),
0543           out.pfneutral[i].intPhi(),
0544           out.pfneutral[i].intId(),
0545           out.pfneutral[i].pack().to_string(16).c_str());
0546     }
0547     dbgPrintf("%s", "FW\n");
0548   }
0549 }
0550 
0551 void l1ct::PFAlgo3Emulator::mergeNeutrals(OutputRegion& out) const {
0552   out.pfphoton.reserve(out.pfphoton.size() + out.pfneutral.size());
0553   out.pfphoton.insert(out.pfphoton.end(), out.pfneutral.begin(), out.pfneutral.end());
0554   out.pfphoton.swap(out.pfneutral);
0555   out.pfphoton.clear();
0556 }