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