File indexing completed on 2024-12-20 03:14:05
0001 #include <memory>
0002
0003 #include "RecoJets/JetProducers/interface/ECFAdder.h"
0004 #include "fastjet/PseudoJet.hh"
0005 #include "fastjet/ClusterSequence.hh"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010
0011 ECFAdder::ECFAdder(const edm::ParameterSet& iConfig)
0012 : src_(iConfig.getParameter<edm::InputTag>("src")),
0013 src_token_(consumes<edm::View<reco::Jet>>(src_)),
0014 Njets_(iConfig.getParameter<std::vector<unsigned>>("Njets")),
0015 cuts_(iConfig.getParameter<std::vector<std::string>>("cuts")),
0016 ecftype_(iConfig.getParameter<std::string>("ecftype")),
0017 alpha_(iConfig.getParameter<double>("alpha")),
0018 beta_(iConfig.getParameter<double>("beta")) {
0019 if (cuts_.size() != Njets_.size()) {
0020 throw cms::Exception("ConfigurationError") << "cuts and Njets must be the same size in ECFAdder" << std::endl;
0021 }
0022
0023 edm::InputTag srcWeights = iConfig.getParameter<edm::InputTag>("srcWeights");
0024 if (!srcWeights.label().empty())
0025 input_weights_token_ = consumes<edm::ValueMap<float>>(srcWeights);
0026
0027 for (std::vector<unsigned>::const_iterator n = Njets_.begin(); n != Njets_.end(); ++n) {
0028 std::ostringstream ecfN_str;
0029 std::shared_ptr<fastjet::FunctionOfPseudoJet<double>> pfunc;
0030
0031 if (ecftype_ == "ECF" || ecftype_.empty()) {
0032 ecfN_str << "ecf" << *n;
0033 pfunc = std::make_shared<fastjet::contrib::EnergyCorrelator>(*n, beta_, fastjet::contrib::EnergyCorrelator::pt_R);
0034 } else if (ecftype_ == "C") {
0035 ecfN_str << "ecfC" << *n;
0036 pfunc = std::make_shared<fastjet::contrib::EnergyCorrelatorCseries>(
0037 *n, beta_, fastjet::contrib::EnergyCorrelator::pt_R);
0038 } else if (ecftype_ == "D") {
0039 ecfN_str << "ecfD" << *n;
0040 pfunc = std::make_shared<fastjet::contrib::EnergyCorrelatorGeneralizedD2>(
0041 alpha_, beta_, fastjet::contrib::EnergyCorrelator::pt_R);
0042 } else if (ecftype_ == "N") {
0043 ecfN_str << "ecfN" << *n;
0044 pfunc = std::make_shared<fastjet::contrib::EnergyCorrelatorNseries>(
0045 *n, beta_, fastjet::contrib::EnergyCorrelator::pt_R);
0046 } else if (ecftype_ == "M") {
0047 ecfN_str << "ecfM" << *n;
0048 pfunc = std::make_shared<fastjet::contrib::EnergyCorrelatorMseries>(
0049 *n, beta_, fastjet::contrib::EnergyCorrelator::pt_R);
0050 } else if (ecftype_ == "U") {
0051 ecfN_str << "ecfU" << *n;
0052 pfunc = std::make_shared<fastjet::contrib::EnergyCorrelatorUseries>(
0053 *n, beta_, fastjet::contrib::EnergyCorrelator::pt_R);
0054 }
0055 variables_.push_back(ecfN_str.str());
0056 produces<edm::ValueMap<float>>(ecfN_str.str());
0057 routine_.push_back(pfunc);
0058
0059 selectors_.push_back(StringCutObjectSelector<reco::Jet>(cuts_[n - Njets_.begin()]));
0060 }
0061 }
0062
0063 void ECFAdder::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0064
0065 edm::Handle<edm::View<reco::Jet>> jets;
0066 iEvent.getByToken(src_token_, jets);
0067
0068
0069 if (!input_weights_token_.isUninitialized())
0070 weightsHandle_ = &iEvent.get(input_weights_token_);
0071
0072 unsigned i = 0;
0073 for (std::vector<unsigned>::const_iterator n = Njets_.begin(); n != Njets_.end(); ++n) {
0074
0075 std::vector<float> ecfN;
0076 ecfN.reserve(jets->size());
0077
0078 for (typename edm::View<reco::Jet>::const_iterator jetIt = jets->begin(); jetIt != jets->end(); ++jetIt) {
0079 edm::Ptr<reco::Jet> jetPtr = jets->ptrAt(jetIt - jets->begin());
0080
0081 float t = -1.0;
0082 if (selectors_[n - Njets_.begin()](*jetIt))
0083 t = getECF(i, jetPtr);
0084
0085 ecfN.push_back(t);
0086 }
0087
0088 auto outT = std::make_unique<edm::ValueMap<float>>();
0089 edm::ValueMap<float>::Filler fillerT(*outT);
0090 fillerT.insert(jets, ecfN.begin(), ecfN.end());
0091 fillerT.fill();
0092
0093 iEvent.put(std::move(outT), variables_[i]);
0094 ++i;
0095 }
0096 }
0097
0098 float ECFAdder::getECF(unsigned index, const edm::Ptr<reco::Jet>& object) const {
0099 std::vector<fastjet::PseudoJet> fjParticles;
0100 for (unsigned k = 0; k < object->numberOfDaughters(); ++k) {
0101 const reco::CandidatePtr& dp = object->daughterPtr(k);
0102 if (dp.isNonnull() && dp.isAvailable()) {
0103
0104 if (dp->numberOfDaughters() == 0) {
0105 if (object->isWeighted()) {
0106 if (input_weights_token_.isUninitialized())
0107 throw cms::Exception("MissingConstituentWeight")
0108 << "ECFAdder: No weights (e.g. PUPPI) given for weighted jet collection" << std::endl;
0109 float w = (*weightsHandle_)[dp];
0110 fjParticles.push_back(fastjet::PseudoJet(dp->px() * w, dp->py() * w, dp->pz() * w, dp->energy() * w));
0111 } else
0112 fjParticles.push_back(fastjet::PseudoJet(dp->px(), dp->py(), dp->pz(), dp->energy()));
0113 } else {
0114 auto subjet = dynamic_cast<reco::Jet const*>(dp.get());
0115 for (unsigned l = 0; l < subjet->numberOfDaughters(); ++l) {
0116 if (subjet != nullptr) {
0117 const reco::CandidatePtr& ddp = subjet->daughterPtr(l);
0118 if (subjet->isWeighted()) {
0119 if (input_weights_token_.isUninitialized())
0120 throw cms::Exception("MissingConstituentWeight")
0121 << "ECFAdder: No weights (e.g. PUPPI) given for weighted jet collection" << std::endl;
0122 float w = (*weightsHandle_)[ddp];
0123 fjParticles.push_back(fastjet::PseudoJet(ddp->px() * w, ddp->py() * w, ddp->pz() * w, ddp->energy() * w));
0124 } else
0125 fjParticles.push_back(fastjet::PseudoJet(ddp->px(), ddp->py(), ddp->pz(), ddp->energy()));
0126 } else {
0127 edm::LogWarning("MissingJetConstituent") << "BasicJet constituent required for ECF computation is missing!";
0128 }
0129 }
0130 }
0131 }
0132 else
0133 edm::LogWarning("MissingJetConstituent") << "Jet constituent required for ECF computation is missing!";
0134 }
0135 if (fjParticles.size() > Njets_[index]) {
0136 return routine_[index]->result(join(fjParticles));
0137 } else {
0138 return -1.0;
0139 }
0140 }
0141
0142
0143 void ECFAdder::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0144 edm::ParameterSetDescription iDesc;
0145 iDesc.setComment("Energy Correlation Functions adder");
0146
0147 iDesc.add<edm::InputTag>("src", edm::InputTag("no default"))->setComment("input collection");
0148 iDesc.add<std::vector<unsigned>>("Njets", {1, 2, 3})->setComment("Number of jets to emulate");
0149 iDesc.add<std::vector<std::string>>("cuts", {"", "", ""})
0150 ->setComment("Jet selections for each N value. Size must match Njets.");
0151 iDesc.add<double>("alpha", 1.0)->setComment("alpha factor, only valid for N2");
0152 iDesc.add<double>("beta", 1.0)->setComment("angularity factor");
0153 iDesc.add<std::string>("ecftype", "")->setComment("ECF type: ECF or empty; C; D; N; M; U;");
0154 iDesc.add<edm::InputTag>("srcWeights", edm::InputTag("puppi"));
0155 descriptions.add("ECFAdder", iDesc);
0156 }
0157
0158 DEFINE_FWK_MODULE(ECFAdder);