IsolationAlgos/plugins/CITKPFIsolationSumProducerForPUPPI.cc

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