Back to home page

Project CMSSW displayed by LXR

 
 

    


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   // read input collection
0065   edm::Handle<edm::View<reco::Jet>> jets;
0066   iEvent.getByToken(src_token_, jets);
0067 
0068   // Get Weights Collection
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     // prepare room for output
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       // Here, the daughters are the "end" node, so this is a PFJet
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 {  // Otherwise, this is a BasicJet, so you need to descend further.
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       }  // end if basic jet
0131     }  // end if daughter pointer is nonnull and available
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 // ParameterSet description for module
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);