Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-05-09 22:40:11 

0001 
0002 #include <vector>
0003 #include <numeric>
0004 #include <string>
0005 
0006 ////////////////////
0007 // FRAMEWORK HEADERS
0008 #include "FWCore/Framework/interface/global/EDProducer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Utilities/interface/InputTag.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 
0014 #include "DataFormats/L1TParticleFlow/interface/PFCandidate.h"
0015 #include "DataFormats/L1TParticleFlow/interface/PFJet.h"
0016 #include "L1Trigger/Phase2L1ParticleFlow/interface/corrector.h"
0017 #include "DataFormats/Math/interface/deltaR.h"
0018 
0019 // bitwise emulation headers
0020 #include "L1Trigger/Phase2L1ParticleFlow/interface/jetmet/L1SeedConePFJetEmulator.h"
0021 #include "DataFormats/L1TParticleFlow/interface/gt_datatypes.h"
0022 
0023 class L1SeedConePFJetProducer : public edm::global::EDProducer<> {
0024 public:
0025   explicit L1SeedConePFJetProducer(const edm::ParameterSet&);
0026   ~L1SeedConePFJetProducer() override;
0027   static void fillDescriptions(edm::ConfigurationDescriptions& description);
0028 
0029 private:
0030   /// ///////////////// ///
0031   /// MANDATORY METHODS ///
0032   void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
0033   /// ///////////////// ///
0034 
0035   const float coneSize;
0036   const bool wideConeJet;
0037   const unsigned nJets;
0038   const bool HW;
0039   const bool debug;
0040   const bool doCorrections;
0041   L1SCJetEmu emulator;
0042   edm::EDGetTokenT<std::vector<l1t::PFCandidate>> l1PFToken;
0043   l1tpf::corrector corrector;
0044 
0045   std::vector<l1t::PFJet> processEvent_SW(std::vector<edm::Ptr<l1t::PFCandidate>>& parts) const;
0046   std::vector<l1t::PFJet> processEvent_HW(std::vector<edm::Ptr<l1t::PFCandidate>>& parts) const;
0047 
0048   l1t::PFJet makeJet_SW(const std::vector<edm::Ptr<l1t::PFCandidate>>& parts,
0049                         const edm::Ptr<l1t::PFCandidate>& seed) const;
0050 
0051   std::pair<std::vector<L1SCJetEmu::Particle>, std::unordered_map<const l1t::PFCandidate*, edm::Ptr<l1t::PFCandidate>>>
0052   convertEDMToHW(std::vector<edm::Ptr<l1t::PFCandidate>>& edmParticles) const;
0053 
0054   std::vector<l1t::PFJet> convertHWToEDM(
0055       std::vector<L1SCJetEmu::Jet> hwJets,
0056       std::unordered_map<const l1t::PFCandidate*, edm::Ptr<l1t::PFCandidate>> constituentMap) const;
0057 };
0058 
0059 L1SeedConePFJetProducer::L1SeedConePFJetProducer(const edm::ParameterSet& cfg)
0060     : coneSize(cfg.getParameter<double>("coneSize")),
0061       wideConeJet(cfg.getParameter<bool>("wideConeJet")),
0062       nJets(cfg.getParameter<unsigned>("nJets")),
0063       HW(cfg.getParameter<bool>("HW")),
0064       debug(cfg.getParameter<bool>("debug")),
0065       doCorrections(cfg.getParameter<bool>("doCorrections")),
0066       emulator(L1SCJetEmu(debug, coneSize, nJets)),
0067       l1PFToken(consumes<std::vector<l1t::PFCandidate>>(cfg.getParameter<edm::InputTag>("L1PFObjects"))) {
0068   produces<l1t::PFJetCollection>();
0069   if (doCorrections) {
0070     corrector = l1tpf::corrector(
0071         cfg.getParameter<std::string>("correctorFile"), cfg.getParameter<std::string>("correctorDir"), -1., debug, HW);
0072   }
0073 }
0074 
0075 void L1SeedConePFJetProducer::produce(edm::StreamID /*unused*/,
0076                                       edm::Event& iEvent,
0077                                       const edm::EventSetup& iSetup) const {
0078   std::unique_ptr<l1t::PFJetCollection> newPFJetCollection(new l1t::PFJetCollection);
0079 
0080   edm::Handle<l1t::PFCandidateCollection> l1PFCandidates;
0081   iEvent.getByToken(l1PFToken, l1PFCandidates);
0082   std::vector<edm::Ptr<l1t::PFCandidate>> particles;
0083   for (unsigned i = 0; i < (*l1PFCandidates).size(); i++) {
0084     particles.push_back(edm::Ptr<l1t::PFCandidate>(l1PFCandidates, i));
0085   }
0086 
0087   std::vector<l1t::PFJet> jets;
0088   if (HW) {
0089     jets = processEvent_HW(particles);
0090   } else {
0091     jets = processEvent_SW(particles);
0092   }
0093   std::sort(jets.begin(), jets.end(), [](l1t::PFJet i, l1t::PFJet j) { return (i.pt() > j.pt()); });
0094   newPFJetCollection->swap(jets);
0095   iEvent.put(std::move(newPFJetCollection));  // Add jets to the event
0096 }
0097 
0098 /////////////
0099 // DESTRUCTOR
0100 L1SeedConePFJetProducer::~L1SeedConePFJetProducer() {}
0101 
0102 l1t::PFJet L1SeedConePFJetProducer::makeJet_SW(const std::vector<edm::Ptr<l1t::PFCandidate>>& parts,
0103                                                const edm::Ptr<l1t::PFCandidate>& seed) const {
0104   auto sumpt = [](float a, const edm::Ptr<l1t::PFCandidate>& b) { return a + b->pt(); };
0105 
0106   // Sum the pt
0107   float pt = std::accumulate(parts.begin(), parts.end(), 0., sumpt);
0108 
0109   // pt weighted d eta
0110   std::vector<float> pt_deta;
0111   pt_deta.resize(parts.size());
0112   std::transform(parts.begin(), parts.end(), pt_deta.begin(), [&seed, &pt](const edm::Ptr<l1t::PFCandidate>& part) {
0113     return (part->pt() / pt) * (part->eta() - seed->eta());  // may have to derefernce seed
0114   });
0115   // Accumulate the pt weighted etas. Init to the seed eta, start accumulating at begin()+1 to skip seed
0116   float eta = std::accumulate(pt_deta.begin() + 1, pt_deta.end(), seed->eta());
0117 
0118   // pt weighted d phi
0119   std::vector<float> pt_dphi;
0120   pt_dphi.resize(parts.size());
0121   std::transform(parts.begin(), parts.end(), pt_dphi.begin(), [&seed, &pt](const edm::Ptr<l1t::PFCandidate>& part) {
0122     return (part->pt() / pt) * reco::deltaPhi(part->phi(), seed->phi());
0123   });
0124   // Accumulate the pt weighted phis. Init to the seed phi, start accumulating at begin()+1 to skip seed
0125   float phi = std::accumulate(pt_dphi.begin() + 1, pt_dphi.end(), seed->phi());
0126 
0127   // Calculate the mass
0128   std::vector<float> en;
0129   en.resize(parts.size());
0130   std::transform(parts.begin(), parts.end(), en.begin(), [](const edm::Ptr<l1t::PFCandidate>& part) {
0131     return std::pow(std::pow((part->pt() * std::cosh(part->eta())), 2) + std::pow(part->mass(), 2), 0.5);
0132   });
0133   float en_tot = std::accumulate(en.begin(), en.end(), 0.0);
0134 
0135   auto sumpx = [](float a, const edm::Ptr<l1t::PFCandidate>& b) { return a + (b->pt() * std::cos(b->phi())); };
0136   float px_tot = std::accumulate(parts.begin(), parts.end(), 0.0, sumpx);
0137 
0138   auto sumpy = [](float a, const edm::Ptr<l1t::PFCandidate>& b) { return a + (b->pt() * std::sin(b->phi())); };
0139   float py_tot = std::accumulate(parts.begin(), parts.end(), 0.0, sumpy);
0140 
0141   auto sumpz = [](float a, const edm::Ptr<l1t::PFCandidate>& b) { return a + (b->pt() * std::sinh(b->eta())); };
0142   float pz_tot = std::accumulate(parts.begin(), parts.end(), 0.0, sumpz);
0143 
0144   float mass = std::sqrt((en_tot * en_tot) - (px_tot * px_tot) - (py_tot * py_tot) - (pz_tot * pz_tot));
0145 
0146   l1t::PFJet jet(pt, eta, phi, mass);
0147 
0148   for (auto it = parts.begin(); it != parts.end(); it++) {
0149     jet.addConstituent(*it);
0150   }
0151 
0152   if (doCorrections) {
0153     jet.calibratePt(corrector.correctedPt(jet.pt(), jet.eta()));
0154   }
0155 
0156   return jet;
0157 }
0158 
0159 std::vector<l1t::PFJet> L1SeedConePFJetProducer::processEvent_SW(std::vector<edm::Ptr<l1t::PFCandidate>>& parts) const {
0160   // The floating point algorithm simulation
0161   std::stable_sort(parts.begin(), parts.end(), [](edm::Ptr<l1t::PFCandidate> i, edm::Ptr<l1t::PFCandidate> j) {
0162     return (i->pt() > j->pt());  // this sorts the candidates by pT
0163   });
0164   std::vector<l1t::PFJet> jets;  // make vector of jets
0165   jets.reserve(nJets);           // reserve enough entries for nJets
0166 
0167   while (!parts.empty() &&
0168          jets.size() < nJets) {  // whilst theres candidates in the array and nJets havent yet been found
0169     edm::Ptr<l1t::PFCandidate> seed =
0170         parts.at(0);  // If use external seeds true, use external seeds, else use highest pt cand
0171 
0172     // Get the particles within a coneSize of the seed
0173     std::vector<edm::Ptr<l1t::PFCandidate>> particlesInCone;
0174     std::copy_if(
0175         parts.begin(), parts.end(), std::back_inserter(particlesInCone), [&](const edm::Ptr<l1t::PFCandidate>& part) {
0176           return reco::deltaR<l1t::PFCandidate, l1t::PFCandidate>(*seed, *part) <= coneSize;
0177         });
0178 
0179     jets.push_back(makeJet_SW(particlesInCone, seed));
0180     // remove the clustered particles
0181     parts.erase(std::remove_if(parts.begin(),
0182                                parts.end(),
0183                                [&](const edm::Ptr<l1t::PFCandidate>& part) {
0184                                  return reco::deltaR<l1t::PFCandidate, l1t::PFCandidate>(*seed, *part) <= coneSize;
0185                                }),
0186                 parts.end());
0187   }
0188 
0189   return jets;
0190 }
0191 
0192 std::vector<l1t::PFJet> L1SeedConePFJetProducer::processEvent_HW(std::vector<edm::Ptr<l1t::PFCandidate>>& parts) const {
0193   // The fixed point emulator
0194   // Convert the EDM format to the hardware format, and call the standalone emulator
0195   std::pair<std::vector<L1SCJetEmu::Particle>, std::unordered_map<const l1t::PFCandidate*, edm::Ptr<l1t::PFCandidate>>>
0196       particles = convertEDMToHW(parts);
0197   std::vector<L1SCJetEmu::Jet> jets = emulator.emulateEvent(particles.first);
0198   return convertHWToEDM(jets, particles.second);
0199 }
0200 
0201 std::pair<std::vector<L1SCJetEmu::Particle>, std::unordered_map<const l1t::PFCandidate*, edm::Ptr<l1t::PFCandidate>>>
0202 L1SeedConePFJetProducer::convertEDMToHW(std::vector<edm::Ptr<l1t::PFCandidate>>& edmParticles) const {
0203   std::vector<l1ct::PuppiObjEmu> hwParticles;
0204   std::unordered_map<const l1t::PFCandidate*, edm::Ptr<l1t::PFCandidate>> candidateMap;
0205   std::for_each(edmParticles.begin(), edmParticles.end(), [&](edm::Ptr<l1t::PFCandidate>& edmParticle) {
0206     l1ct::PuppiObjEmu particle;
0207     particle.initFromBits(edmParticle->encodedPuppi64());
0208     particle.srcCand = edmParticle.get();
0209     candidateMap.insert(std::make_pair(edmParticle.get(), edmParticle));
0210     hwParticles.push_back(particle);
0211   });
0212   return std::make_pair(hwParticles, candidateMap);
0213 }
0214 
0215 std::vector<l1t::PFJet> L1SeedConePFJetProducer::convertHWToEDM(
0216     std::vector<L1SCJetEmu::Jet> hwJets,
0217     std::unordered_map<const l1t::PFCandidate*, edm::Ptr<l1t::PFCandidate>> constituentMap) const {
0218   std::vector<l1t::PFJet> edmJets;
0219   std::for_each(hwJets.begin(), hwJets.end(), [&](L1SCJetEmu::Jet jet) {
0220     if (doCorrections) {
0221       float correctedPt = corrector.correctedPt(jet.floatPt(), jet.floatEta());
0222       jet.hwPt = correctedPt;
0223     }
0224     l1gt::Jet gtJet = jet.toGT();
0225     l1t::PFJet edmJet(l1gt::Scales::floatPt(gtJet.v3.pt),
0226                       l1gt::Scales::floatEta(gtJet.v3.eta),
0227                       l1gt::Scales::floatPhi(gtJet.v3.phi),
0228                       0,  // dont set the mass by default
0229                       gtJet.v3.pt.V,
0230                       gtJet.v3.eta.V,
0231                       gtJet.v3.phi.V);
0232     edmJet.setEncodedJet(l1t::PFJet::HWEncoding::CT, jet.pack());
0233 
0234     if (wideConeJet) {
0235       edmJet.setEncodedJet(l1t::PFJet::HWEncoding::GTWide, jet.toGTWide().pack());
0236       edmJet.setMass(std::sqrt(l1gt::Scales::floatMassSq(jet.toGTWide().hwMassSq)));
0237     } else {
0238       edmJet.setEncodedJet(l1t::PFJet::HWEncoding::GT, jet.toGT().pack());
0239     }
0240 
0241     // get back the references to the constituents
0242     std::vector<edm::Ptr<l1t::PFCandidate>> constituents;
0243     std::for_each(jet.constituents.begin(), jet.constituents.end(), [&](auto constituent) {
0244       edmJet.addConstituent(constituentMap[constituent.srcCand]);
0245     });
0246     edmJets.push_back(edmJet);
0247   });
0248   return edmJets;
0249 }
0250 
0251 void L1SeedConePFJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0252   edm::ParameterSetDescription desc;
0253   desc.add<edm::InputTag>("L1PFObjects", edm::InputTag("l1tLayer1", "Puppi"));
0254   desc.add<uint32_t>("nJets", 16);
0255   desc.add<double>("coneSize", 0.4);
0256   desc.add<bool>("wideConeJet", false);
0257   desc.add<bool>("HW", false);
0258   desc.add<bool>("debug", false);
0259   desc.add<bool>("doCorrections", false);
0260   desc.add<std::string>("correctorFile", "");
0261   desc.add<std::string>("correctorDir", "");
0262   descriptions.addWithDefaultLabel(desc);
0263 }
0264 
0265 DEFINE_FWK_MODULE(L1SeedConePFJetProducer);
[ 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!