Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:21:32 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!