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