File indexing completed on 2023-10-25 09:59:51
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
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
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
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
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
0218
0219
0220 DEFINE_FWK_MODULE(SubjetFilterJetProducer);