File indexing completed on 2024-09-07 04:37:34
0001 #include <memory>
0002
0003 #include "RecoJets/JetProducers/interface/NjettinessAdder.h"
0004
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006
0007 NjettinessAdder::NjettinessAdder(const edm::ParameterSet& iConfig)
0008 : src_(iConfig.getParameter<edm::InputTag>("src")),
0009 src_token_(consumes<edm::View<reco::Jet>>(src_)),
0010 Njets_(iConfig.getParameter<std::vector<unsigned>>("Njets")),
0011 measureDefinition_(iConfig.getParameter<unsigned>("measureDefinition")),
0012 beta_(iConfig.getParameter<double>("beta")),
0013 R0_(iConfig.getParameter<double>("R0")),
0014 Rcutoff_(iConfig.getParameter<double>("Rcutoff")),
0015 axesDefinition_(iConfig.getParameter<unsigned>("axesDefinition")),
0016 nPass_(iConfig.getParameter<int>("nPass")),
0017 akAxesR0_(iConfig.getParameter<double>("akAxesR0")) {
0018 edm::InputTag srcWeights = iConfig.getParameter<edm::InputTag>("srcWeights");
0019 if (!srcWeights.label().empty())
0020 input_weights_token_ = consumes<edm::ValueMap<float>>(srcWeights);
0021
0022 for (std::vector<unsigned>::const_iterator n = Njets_.begin(); n != Njets_.end(); ++n) {
0023 std::ostringstream tauN_str;
0024 tauN_str << "tau" << *n;
0025
0026 produces<edm::ValueMap<float>>(tauN_str.str());
0027 }
0028
0029
0030 fastjet::contrib::NormalizedMeasure normalizedMeasure(beta_, R0_);
0031 fastjet::contrib::UnnormalizedMeasure unnormalizedMeasure(beta_);
0032 fastjet::contrib::OriginalGeometricMeasure geometricMeasure(beta_);
0033 fastjet::contrib::NormalizedCutoffMeasure normalizedCutoffMeasure(beta_, R0_, Rcutoff_);
0034 fastjet::contrib::UnnormalizedCutoffMeasure unnormalizedCutoffMeasure(beta_, Rcutoff_);
0035
0036
0037 fastjet::contrib::MeasureDefinition const* measureDef = nullptr;
0038 switch (measureDefinition_) {
0039 case UnnormalizedMeasure:
0040 measureDef = &unnormalizedMeasure;
0041 break;
0042 case OriginalGeometricMeasure:
0043 measureDef = &geometricMeasure;
0044 break;
0045 case NormalizedCutoffMeasure:
0046 measureDef = &normalizedCutoffMeasure;
0047 break;
0048 case UnnormalizedCutoffMeasure:
0049 measureDef = &unnormalizedCutoffMeasure;
0050 break;
0051
0052 case NormalizedMeasure:
0053 default:
0054 measureDef = &normalizedMeasure;
0055 break;
0056 }
0057
0058
0059 fastjet::contrib::KT_Axes kt_axes;
0060 fastjet::contrib::CA_Axes ca_axes;
0061 fastjet::contrib::AntiKT_Axes antikt_axes(akAxesR0_);
0062 fastjet::contrib::WTA_KT_Axes wta_kt_axes;
0063 fastjet::contrib::WTA_CA_Axes wta_ca_axes;
0064 fastjet::contrib::OnePass_KT_Axes onepass_kt_axes;
0065 fastjet::contrib::OnePass_CA_Axes onepass_ca_axes;
0066 fastjet::contrib::OnePass_AntiKT_Axes onepass_antikt_axes(akAxesR0_);
0067 fastjet::contrib::OnePass_WTA_KT_Axes onepass_wta_kt_axes;
0068 fastjet::contrib::OnePass_WTA_CA_Axes onepass_wta_ca_axes;
0069 fastjet::contrib::MultiPass_Axes multipass_axes(nPass_);
0070
0071 fastjet::contrib::AxesDefinition const* axesDef = nullptr;
0072 switch (axesDefinition_) {
0073 case KT_Axes:
0074 default:
0075 axesDef = &kt_axes;
0076 break;
0077 case CA_Axes:
0078 axesDef = &ca_axes;
0079 break;
0080 case AntiKT_Axes:
0081 axesDef = &antikt_axes;
0082 break;
0083 case WTA_KT_Axes:
0084 axesDef = &wta_kt_axes;
0085 break;
0086 case WTA_CA_Axes:
0087 axesDef = &wta_ca_axes;
0088 break;
0089 case OnePass_KT_Axes:
0090 axesDef = &onepass_kt_axes;
0091 break;
0092 case OnePass_CA_Axes:
0093 axesDef = &onepass_ca_axes;
0094 break;
0095 case OnePass_AntiKT_Axes:
0096 axesDef = &onepass_antikt_axes;
0097 break;
0098 case OnePass_WTA_KT_Axes:
0099 axesDef = &onepass_wta_kt_axes;
0100 break;
0101 case OnePass_WTA_CA_Axes:
0102 axesDef = &onepass_wta_ca_axes;
0103 break;
0104 case MultiPass_Axes:
0105 axesDef = &multipass_axes;
0106 break;
0107 };
0108
0109 routine_ = std::make_unique<fastjet::contrib::Njettiness>(*axesDef, *measureDef);
0110 }
0111
0112 void NjettinessAdder::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0113
0114 edm::Handle<edm::View<reco::Jet>> jets;
0115 iEvent.getByToken(src_token_, jets);
0116
0117
0118 if (!input_weights_token_.isUninitialized())
0119 weightsHandle_ = &iEvent.get(input_weights_token_);
0120
0121 for (std::vector<unsigned>::const_iterator n = Njets_.begin(); n != Njets_.end(); ++n) {
0122 std::ostringstream tauN_str;
0123 tauN_str << "tau" << *n;
0124
0125
0126 std::vector<float> tauN;
0127 tauN.reserve(jets->size());
0128
0129 for (typename edm::View<reco::Jet>::const_iterator jetIt = jets->begin(); jetIt != jets->end(); ++jetIt) {
0130 edm::Ptr<reco::Jet> jetPtr = jets->ptrAt(jetIt - jets->begin());
0131
0132 float t = getTau(*n, jetPtr);
0133
0134 tauN.push_back(t);
0135 }
0136
0137 auto outT = std::make_unique<edm::ValueMap<float>>();
0138 edm::ValueMap<float>::Filler fillerT(*outT);
0139 fillerT.insert(jets, tauN.begin(), tauN.end());
0140 fillerT.fill();
0141
0142 iEvent.put(std::move(outT), tauN_str.str());
0143 }
0144 }
0145
0146 float NjettinessAdder::getTau(unsigned num, const edm::Ptr<reco::Jet>& object) const {
0147 std::vector<fastjet::PseudoJet> fjParticles;
0148 for (unsigned k = 0; k < object->numberOfDaughters(); ++k) {
0149 const reco::CandidatePtr& dp = object->daughterPtr(k);
0150 if (dp.isNonnull() && dp.isAvailable()) {
0151
0152 if (dp->numberOfDaughters() == 0) {
0153 if (object->isWeighted()) {
0154 if (input_weights_token_.isUninitialized())
0155 throw cms::Exception("MissingConstituentWeight")
0156 << "ECFAdder: No weights (e.g. PUPPI) given for weighted jet collection" << std::endl;
0157 float w = (*weightsHandle_)[dp];
0158 fjParticles.push_back(fastjet::PseudoJet(dp->px() * w, dp->py() * w, dp->pz() * w, dp->energy() * w));
0159 } else
0160 fjParticles.push_back(fastjet::PseudoJet(dp->px(), dp->py(), dp->pz(), dp->energy()));
0161 } else {
0162 auto subjet = dynamic_cast<reco::Jet const*>(dp.get());
0163 for (unsigned l = 0; l < subjet->numberOfDaughters(); ++l) {
0164 if (subjet != nullptr) {
0165 const reco::CandidatePtr& ddp = subjet->daughterPtr(l);
0166 if (subjet->isWeighted()) {
0167 if (input_weights_token_.isUninitialized())
0168 throw cms::Exception("MissingConstituentWeight")
0169 << "ECFAdder: No weights (e.g. PUPPI) given for weighted jet collection" << std::endl;
0170 float w = (*weightsHandle_)[ddp];
0171 fjParticles.push_back(fastjet::PseudoJet(ddp->px() * w, ddp->py() * w, ddp->pz() * w, ddp->energy() * w));
0172 } else
0173 fjParticles.push_back(fastjet::PseudoJet(ddp->px(), ddp->py(), ddp->pz(), ddp->energy()));
0174 } else {
0175 edm::LogWarning("MissingJetConstituent")
0176 << "BasicJet constituent required for N-jettiness computation is missing!";
0177 }
0178 }
0179 }
0180 }
0181 else
0182 edm::LogWarning("MissingJetConstituent") << "Jet constituent required for N-jettiness computation is missing!";
0183 }
0184
0185 return routine_->getTau(num, fjParticles);
0186 }
0187
0188 DEFINE_FWK_MODULE(NjettinessAdder);