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