IsolationAlgos/plugins/CITKPFIsolationSumProducer.cc

0001 #ifndef PhysicsTools_IsolationAlgos_CITKIsolationSumProducer_H
0002 #define PhysicsTools_IsolationAlgos_CITKIsolationSumProducer_H
0003
0004 #include "FWCore/Framework/interface/stream/EDProducer.h"
0005 #include "FWCore/Framework/interface/ConsumesCollector.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007
0008 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0009 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0010 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
0011 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
0012 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0013 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0014 #include "PhysicsTools/IsolationAlgos/interface/EventDependentAbsVeto.h"
0015
0016 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0017 #include "DataFormats/Candidate/interface/Candidate.h"
0018 #include "PhysicsTools/IsolationAlgos/interface/CITKIsolationConeDefinitionBase.h"
0019 #include "DataFormats/Common/interface/OwnVector.h"
0020
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 #include "FWCore/Framework/interface/Event.h"
0023 #include "FWCore/Framework/interface/EventSetup.h"
0024
0025 #include <string>
0026 #include <unordered_map>
0027
0028 namespace edm {
0029   class Event;
0030 }
0031 namespace edm {
0032   class EventSetup;
0033 }
0034
0035 namespace citk {
0036   class PFIsolationSumProducer : public edm::stream::EDProducer<> {
0037   public:
0038     PFIsolationSumProducer(const edm::ParameterSet&);
0039
0040     ~PFIsolationSumProducer() override {}
0041
0042     void beginLuminosityBlock(const edm::LuminosityBlock&, const edm::EventSetup&) final;
0043
0044     void produce(edm::Event&, const edm::EventSetup&) final;
0045
0046     static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0047
0048   private:
0049     // datamembers
0050     static constexpr unsigned kNPFTypes = 8;
0051     typedef std::unordered_map<std::string, int> TypeMap;
0052     typedef std::vector<std::unique_ptr<IsolationConeDefinitionBase>> IsoTypes;
0053     typedef edm::View<reco::Candidate> CandView;
0054     const TypeMap _typeMap;
0055     edm::EDGetTokenT<CandView> _to_isolate, _isolate_with;
0056     // indexed by pf candidate type
0057     std::array<IsoTypes, kNPFTypes> _isolation_types;
0058     std::array<std::vector<std::string>, kNPFTypes> _product_names;
0059   };
0060 }  // namespace citk
0061
0062 typedef citk::PFIsolationSumProducer CITKPFIsolationSumProducer;
0063
0064 DEFINE_FWK_MODULE(CITKPFIsolationSumProducer);
0065
0066 namespace citk {
0067   PFIsolationSumProducer::PFIsolationSumProducer(const edm::ParameterSet& c)
0068       : _typeMap({{"h+", 1}, {"h0", 5}, {"gamma", 4}, {"electron", 2}, {"muon", 3}, {"HFh", 6}, {"HFgamma", 7}}) {
0069     _to_isolate = consumes<CandView>(c.getParameter<edm::InputTag>("srcToIsolate"));
0070     _isolate_with = consumes<CandView>(c.getParameter<edm::InputTag>("srcForIsolationCone"));
0071     const std::vector<edm::ParameterSet>& isoDefs = c.getParameterSetVector("isolationConeDefinitions");
0072     for (const auto& isodef : isoDefs) {
0073       const std::string& name = isodef.getParameter<std::string>("isolationAlgo");
0074       const float coneSize = isodef.getParameter<double>("coneSize");
0075       char buf[50];
0076       sprintf(buf, "DR%.2f", coneSize);
0077       std::string coneName(buf);
0078       auto decimal = coneName.find('.');
0079       if (decimal != std::string::npos)
0080         coneName.erase(decimal, 1);
0081       const std::string& isotype = isodef.getParameter<std::string>("isolateAgainst");
0082       auto theisolator = CITKIsolationConeDefinitionFactory::get()->create(name, isodef);
0083       theisolator->setConsumes(consumesCollector());
0084       const auto thetype = _typeMap.find(isotype);
0085       if (thetype == _typeMap.end()) {
0086         throw cms::Exception("InvalidIsolationType") << "Isolation type: " << isotype << " is not available in the "
0087                                                      << "list of allowed isolations!.";
0088       }
0089       const std::string dash("-");
0090       std::string pname = isotype + dash + coneName + dash + theisolator->additionalCode();
0091       _product_names[thetype->second].emplace_back(pname);
0092       produces<edm::ValueMap<float>>(pname);
0093       _isolation_types[thetype->second].emplace_back(std::move(theisolator));
0094     }
0095   }
0096
0097   void PFIsolationSumProducer::beginLuminosityBlock(const edm::LuminosityBlock&, const edm::EventSetup& es) {
0098     for (const auto& isolators_for_type : _isolation_types) {
0099       for (const auto& isolator : isolators_for_type) {
0100         isolator->getEventSetupInfo(es);
0101       }
0102     }
0103   }
0104
0105   void PFIsolationSumProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
0106     typedef std::unique_ptr<edm::ValueMap<float>> product_type;
0107     typedef std::vector<float> product_values;
0108     edm::Handle<CandView> to_isolate;
0109     edm::Handle<CandView> isolate_with;
0110     ev.getByToken(_to_isolate, to_isolate);
0111     ev.getByToken(_isolate_with, isolate_with);
0112     // the list of value vectors indexed as "to_isolate"
0113     std::array<std::vector<product_values>, kNPFTypes> the_values;
0114     // get extra event info and setup value cache
0115     unsigned i = 0;
0116     for (const auto& isolators_for_type : _isolation_types) {
0117       the_values[i++].resize(isolators_for_type.size());
0118       for (const auto& isolator : isolators_for_type) {
0119         isolator->getEventInfo(ev);
0120       }
0121     }
0122     reco::PFCandidate helper;  // to translate pdg id to type
0123     // loop over the candidates we are isolating and fill the values
0124     for (size_t c = 0; c < to_isolate->size(); ++c) {
0125       auto cand_to_isolate = to_isolate->ptrAt(c);
0126       std::array<std::vector<float>, kNPFTypes> cand_values;
0127       unsigned k = 0;
0128       for (const auto& isolators_for_type : _isolation_types) {
0129         cand_values[k].resize(isolators_for_type.size());
0130         for (auto& value : cand_values[k])
0131           value = 0.0;
0132         ++k;
0133       }
0134       for (size_t ic = 0; ic < isolate_with->size(); ++ic) {
0135         auto isocand = isolate_with->ptrAt(ic);
0136         auto isotype = helper.translatePdgIdToType(isocand->pdgId());
0137         const auto& isolations = _isolation_types[isotype];
0138         for (unsigned i = 0; i < isolations.size(); ++i) {
0139           if (isolations[i]->isInIsolationCone(cand_to_isolate, isocand)) {
0140             cand_values[isotype][i] += isocand->pt();
0141           }
0142         }
0143       }
0144       // add this candidate to isolation value list
0145       for (unsigned i = 0; i < kNPFTypes; ++i) {
0146         for (unsigned j = 0; j < cand_values[i].size(); ++j) {
0147           the_values[i][j].push_back(cand_values[i][j]);
0148         }
0149       }
0150     }
0151     // fill and put all products
0152     for (unsigned i = 0; i < kNPFTypes; ++i) {
0153       for (unsigned j = 0; j < the_values[i].size(); ++j) {
0154         product_type the_product(new edm::ValueMap<float>);
0155         edm::ValueMap<float>::Filler fillerprod(*the_product);
0156         fillerprod.insert(to_isolate, the_values[i][j].begin(), the_values[i][j].end());
0157         fillerprod.fill();
0158         ev.put(std::move(the_product), _product_names[i][j]);
0159       }
0160     }
0161   }
0162
0163   // ParameterSet description for module
0164   void PFIsolationSumProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0165     edm::ParameterSetDescription iDesc;
0166     iDesc.setComment("isolation sum producer");
0167
0168     iDesc.add<edm::InputTag>("srcToIsolate", edm::InputTag("no default"))
0169         ->setComment("calculate isolation for this collection");
0170     iDesc.add<edm::InputTag>("srcForIsolationCone", edm::InputTag("no default"))
0171         ->setComment("collection for the isolation calculation: like particleFlow ");
0172
0173     edm::ParameterSetDescription descIsoConeDefinitions;
0174     descIsoConeDefinitions.add<std::string>("isolationAlgo", "no default");
0175     descIsoConeDefinitions.add<double>("coneSize", 0.3);
0176     descIsoConeDefinitions.add<std::string>("isolateAgainst", "no default");
0177     descIsoConeDefinitions.add<std::vector<unsigned>>("miniAODVertexCodes", {2, 3});
0178     descIsoConeDefinitions.addOptional<double>("VetoConeSizeBarrel", 0.0);
0179     descIsoConeDefinitions.addOptional<double>("VetoConeSizeEndcaps", 0.0);
0180     descIsoConeDefinitions.addOptional<int>("vertexIndex", 0);
0181     descIsoConeDefinitions.addOptional<edm::InputTag>("particleBasedIsolation", edm::InputTag("no default"))
0182         ->setComment("map for footprint removal that is used for photons");
0183
0184     std::vector<edm::ParameterSet> isolationConeDefinitions;
0185     edm::ParameterSet chargedHadrons, neutralHadrons, photons;
0186     isolationConeDefinitions.push_back(chargedHadrons);
0187     isolationConeDefinitions.push_back(neutralHadrons);
0188     isolationConeDefinitions.push_back(photons);
0189     iDesc.addVPSet("isolationConeDefinitions", descIsoConeDefinitions, isolationConeDefinitions);
0190
0191     descriptions.add("CITKPFIsolationSumProducer", iDesc);
0192   }
0193 }  // namespace citk
0194
0195 #endif