JetProducers/plugins/SubjetFilterJetProducer.cc

0001 ////////////////////////////////////////////////////////////////////////////////
0002 //
0003 // SubjetFilterJetProducer
0004 // -----------------------
0005 //
0006 // For a description of the Subjet/Filter algorithm, see e.g.
0007 // http://arXiv.org/abs/0802.2470
0008 //
0009 // This implementation is largely based on fastjet_boosted_higgs.cc provided
0010 // with the fastjet package.
0011 //
0012 // If you did, you know that the algorithm produces a 'fatjet', two
0013 // associated 'subjets', and two or three associated 'filterjets'. This producer
0014 // will therefore write three corresponding jet collections. The association
0015 // between a fatjet and its subjets and filterjets is done by making all (!)
0016 // of them daughters of the fatjet, such that the first two daughter jets always
0017 // correspond to the subjets, while all remaining correspond to the filterjets.
0018 //
0019 // The real work is done in RecoJets/JetAlgorithms/src/SubjetFilterAlgorithm.cc
0020 //
0021 // see: https://twiki.cern.ch/twiki/bin/view/CMS/SWGuideSubjetFilterJetProducer
0022 //
0023 //                       David Lopes-Pegna                 <david.lopes@cern.ch>
0024 //            25/11/2009 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
0025 ////////////////////////////////////////////////////////////////////////////////
0026
0027 #include "SubjetFilterJetProducer.h"
0028 #include "RecoJets/JetProducers/interface/JetSpecific.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030
0031 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0032 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0033 using namespace std;
0034
0035 ////////////////////////////////////////////////////////////////////////////////
0036 // construction / destruction
0037 ////////////////////////////////////////////////////////////////////////////////
0038
0039 //______________________________________________________________________________
0040 SubjetFilterJetProducer::SubjetFilterJetProducer(const edm::ParameterSet& iConfig)
0041     : VirtualJetProducer(iConfig),
0042       alg_(iConfig.getParameter<string>("@module_label"),
0043            iConfig.getParameter<string>("jetAlgorithm"),
0044            iConfig.getParameter<unsigned>("nFatMax"),
0045            iConfig.getParameter<double>("rParam"),
0046            iConfig.getParameter<double>("rFilt"),
0047            iConfig.getParameter<double>("jetPtMin"),
0048            iConfig.getParameter<double>("massDropCut"),
0049            iConfig.getParameter<double>("asymmCut"),
0050            iConfig.getParameter<bool>("asymmCutLater"),
0051            iConfig.getParameter<bool>("doAreaFastjet"),
0052            iConfig.getParameter<double>("Ghost_EtaMax"),
0053            iConfig.getParameter<int>("Active_Area_Repeats"),
0054            iConfig.getParameter<double>("GhostArea"),
0055            iConfig.getUntrackedParameter<bool>("verbose", false)) {
0056   produces<reco::BasicJetCollection>("fat");
0057   makeProduces(moduleLabel_, "sub");
0058   makeProduces(moduleLabel_, "filter");
0059 }
0060
0061 //______________________________________________________________________________
0062 SubjetFilterJetProducer::~SubjetFilterJetProducer() {}
0063
0064 ////////////////////////////////////////////////////////////////////////////////
0065 // implementation of member functions
0066 ////////////////////////////////////////////////////////////////////////////////
0067
0068 //______________________________________________________________________________
0069 void SubjetFilterJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0070   VirtualJetProducer::produce(iEvent, iSetup);
0071 }
0072
0073 //______________________________________________________________________________
0074 void SubjetFilterJetProducer::endJob() { cout << alg_.summary() << endl; }
0075
0076 //______________________________________________________________________________
0077 void SubjetFilterJetProducer::runAlgorithm(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0078   alg_.run(fjInputs_, fjCompoundJets_, iSetup);
0079 }
0080
0081 //______________________________________________________________________________
0082 void SubjetFilterJetProducer::inputTowers() {
0083   fjCompoundJets_.clear();
0084   VirtualJetProducer::inputTowers();
0085 }
0086
0087 //______________________________________________________________________________
0088 void SubjetFilterJetProducer::output(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0089   // Write jets and constitutents. Will use fjCompoundJets_.
0090   switch (jetTypeE) {
0091     case JetType::CaloJet:
0092       writeCompoundJets<reco::CaloJet>(iEvent, iSetup);
0093       break;
0094     case JetType::PFJet:
0095       writeCompoundJets<reco::PFJet>(iEvent, iSetup);
0096       break;
0097     case JetType::GenJet:
0098       writeCompoundJets<reco::GenJet>(iEvent, iSetup);
0099       break;
0100     case JetType::BasicJet:
0101       writeCompoundJets<reco::BasicJet>(iEvent, iSetup);
0102       break;
0103     default:
0104       edm::LogError("InvalidInput") << " invalid jet type in SubjetFilterJetProducer\n";
0105       break;
0106   };
0107 }
0108
0109 //______________________________________________________________________________
0110 template <class T>
0111 void SubjetFilterJetProducer::writeCompoundJets(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0112   auto fatJets = std::make_unique<reco::BasicJetCollection>();
0113   auto subJets = std::make_unique<vector<T>>();
0114   auto filterJets = std::make_unique<vector<T>>();
0115
0116   edm::OrphanHandle<vector<T>> subJetsAfterPut;
0117   edm::OrphanHandle<vector<T>> filterJetsAfterPut;
0118
0119   vector<vector<int>> subIndices(fjCompoundJets_.size());
0120   vector<vector<int>> filterIndices(fjCompoundJets_.size());
0121
0122   vector<math::XYZTLorentzVector> p4FatJets;
0123   vector<double> areaFatJets;
0124
0125   vector<CompoundPseudoJet>::const_iterator itBegin(fjCompoundJets_.begin());
0126   vector<CompoundPseudoJet>::const_iterator itEnd(fjCompoundJets_.end());
0127   vector<CompoundPseudoJet>::const_iterator it(itBegin);
0128
0129   [[maybe_unused]] const CaloGeometry* pGeometry = nullptr;
0130   [[maybe_unused]] const HcalTopology* pTopology = nullptr;
0131   if constexpr (std::is_same_v<T, reco::CaloJet>) {
0132     pGeometry = &getGeometry(iSetup);
0133     pTopology = &getTopology(iSetup);
0134   }
0135
0136   for (; it != itEnd; ++it) {
0137     int jetIndex = it - itBegin;
0138     fastjet::PseudoJet fatJet = it->hardJet();
0139     p4FatJets.push_back(math::XYZTLorentzVector(fatJet.px(), fatJet.py(), fatJet.pz(), fatJet.e()));
0140     areaFatJets.push_back(it->hardJetArea());
0141
0142     vector<CompoundPseudoSubJet>::const_iterator itSubBegin(it->subjets().begin());
0143     vector<CompoundPseudoSubJet>::const_iterator itSubEnd(it->subjets().end());
0144     vector<CompoundPseudoSubJet>::const_iterator itSub(itSubBegin);
0145
0146     for (; itSub != itSubEnd; ++itSub) {
0147       int subJetIndex = itSub - itSubBegin;
0148       fastjet::PseudoJet fjSubJet = itSub->subjet();
0149       math::XYZTLorentzVector p4SubJet(fjSubJet.px(), fjSubJet.py(), fjSubJet.pz(), fjSubJet.e());
0150       reco::Particle::Point point(0, 0, 0);
0151
0152       vector<reco::CandidatePtr> subJetConstituents;
0153       const vector<int>& subJetConstituentIndices = itSub->constituents();
0154       vector<int>::const_iterator itIndexBegin(subJetConstituentIndices.begin());
0155       vector<int>::const_iterator itIndexEnd(subJetConstituentIndices.end());
0156       vector<int>::const_iterator itIndex(itIndexBegin);
0157       for (; itIndex != itIndexEnd; ++itIndex)
0158         if ((*itIndex) < static_cast<int>(inputs_.size()))
0159           subJetConstituents.push_back(inputs_[*itIndex]);
0160
0161       T subJet;
0162       if constexpr (std::is_same_v<T, reco::CaloJet>) {
0163         reco::writeSpecific(subJet, p4SubJet, point, subJetConstituents, *pGeometry, *pTopology);
0164       } else {
0165         reco::writeSpecific(subJet, p4SubJet, point, subJetConstituents);
0166       }
0167       subJet.setJetArea(itSub->subjetArea());
0168
0169       if (subJetIndex < 2) {
0170         subIndices[jetIndex].push_back(subJets->size());
0171         subJets->push_back(subJet);
0172       } else {
0173         filterIndices[jetIndex].push_back(filterJets->size());
0174         filterJets->push_back(subJet);
0175       }
0176     }
0177   }
0178
0179   subJetsAfterPut = iEvent.put(std::move(subJets), "sub");
0180   filterJetsAfterPut = iEvent.put(std::move(filterJets), "filter");
0181
0182   vector<math::XYZTLorentzVector>::const_iterator itP4Begin(p4FatJets.begin());
0183   vector<math::XYZTLorentzVector>::const_iterator itP4End(p4FatJets.end());
0184   vector<math::XYZTLorentzVector>::const_iterator itP4(itP4Begin);
0185   for (; itP4 != itP4End; ++itP4) {
0186     int fatIndex = itP4 - itP4Begin;
0187     vector<int>& fatToSub = subIndices[fatIndex];
0188     vector<int>& fatToFilter = filterIndices[fatIndex];
0189
0190     vector<reco::CandidatePtr> i_fatJetConstituents;
0191
0192     vector<int>::const_iterator itSubBegin(fatToSub.begin());
0193     vector<int>::const_iterator itSubEnd(fatToSub.end());
0194     vector<int>::const_iterator itSub(itSubBegin);
0195     for (; itSub != itSubEnd; ++itSub) {
0196       reco::CandidatePtr candPtr(subJetsAfterPut, (*itSub), false);
0197       i_fatJetConstituents.push_back(candPtr);
0198     }
0199
0200     vector<int>::const_iterator itFilterBegin(fatToFilter.begin());
0201     vector<int>::const_iterator itFilterEnd(fatToFilter.end());
0202     vector<int>::const_iterator itFilter(itFilterBegin);
0203     for (; itFilter != itFilterEnd; ++itFilter) {
0204       reco::CandidatePtr candPtr(filterJetsAfterPut, (*itFilter), false);
0205       i_fatJetConstituents.push_back(candPtr);
0206     }
0207
0208     reco::Particle::Point point(0, 0, 0);
0209     fatJets->push_back(reco::BasicJet((*itP4), point, i_fatJetConstituents));
0210     fatJets->back().setJetArea(areaFatJets[fatIndex]);
0211   }
0212
0213   iEvent.put(std::move(fatJets), "fat");
0214 }
0215
0216 ////////////////////////////////////////////////////////////////////////////////
0217 // DEFINE THIS AS A CMSSW FWK PLUGIN
0218 ////////////////////////////////////////////////////////////////////////////////
0219
0220 DEFINE_FWK_MODULE(SubjetFilterJetProducer);