EgammaTools/plugins/HGCalPhotonIDValueMapProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:   RecoEgamma/EgammaTools
0004 // Class:    HGCalPhotonIDValueMapProducer
0005 //
0006 /**\class HGCalPhotonIDValueMapProducer HGCalPhotonIDValueMapProducer.cc RecoEgamma/EgammaTools/plugins/HGCalPhotonIDValueMapProducer.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011     [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Nicholas Charles Smith
0015 //      Created:  Wed, 05 Apr 2017 12:17:43 GMT
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/stream/EDProducer.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/Utilities/interface/StreamID.h"
0031
0032 #include "DataFormats/Common/interface/ValueMap.h"
0033
0034 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0035 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
0036
0037 #include "RecoEgamma/EgammaTools/interface/HGCalEgammaIDHelper.h"
0038 #include "RecoEgamma/EgammaTools/interface/LongDeps.h"
0039
0040 class HGCalPhotonIDValueMapProducer : public edm::stream::EDProducer<> {
0041 public:
0042   explicit HGCalPhotonIDValueMapProducer(const edm::ParameterSet&);
0043   ~HGCalPhotonIDValueMapProducer() override;
0044
0045   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0046
0047 private:
0048   void beginStream(edm::StreamID) override;
0049   void produce(edm::Event&, const edm::EventSetup&) override;
0050   void endStream() override;
0051
0052   // ----------member data ---------------------------
0053   edm::EDGetTokenT<edm::View<reco::Photon>> photonsToken_;
0054   float radius_;
0055
0056   static const std::vector<std::string> valuesProduced_;
0057   std::map<const std::string, std::vector<float>> maps_;
0058
0059   std::unique_ptr<HGCalEgammaIDHelper> phoIDHelper_;
0060 };
0061
0062 // All the ValueMap names to output are defined in the auto-generated python cfi
0063 // so that potential consumers can configure themselves in a simple manner
0064 // Would be cool to use compile-time validation, but need constexpr strings, e.g. std::string_view in C++17
0065 const std::vector<std::string> HGCalPhotonIDValueMapProducer::valuesProduced_ = {
0066     "seedEt",       "seedEnergy",    "seedEnergyEE",  "seedEnergyFH",  "seedEnergyBH",
0067     "pcaEig1",      "pcaEig2",       "pcaEig3",       "pcaSig1",       "pcaSig2",
0068     "pcaSig3",      "sigmaUU",       "sigmaVV",       "sigmaEE",       "sigmaPP",
0069     "nLayers",      "firstLayer",    "lastLayer",     "e4oEtot",       "layerEfrac10",
0070     "layerEfrac90", "measuredDepth", "expectedDepth", "expectedSigma", "depthCompatibility",
0071     "caloIsoRing0", "caloIsoRing1",  "caloIsoRing2",  "caloIsoRing3",  "caloIsoRing4",
0072 };
0073
0074 HGCalPhotonIDValueMapProducer::HGCalPhotonIDValueMapProducer(const edm::ParameterSet& iConfig)
0075     : photonsToken_(consumes(iConfig.getParameter<edm::InputTag>("photons"))),
0076       radius_(iConfig.getParameter<double>("pcaRadius")) {
0077   for (const auto& key : valuesProduced_) {
0078     maps_[key] = {};
0079     produces<edm::ValueMap<float>>(key);
0080   }
0081
0082   phoIDHelper_ = std::make_unique<HGCalEgammaIDHelper>(iConfig, consumesCollector());
0083 }
0084
0085 HGCalPhotonIDValueMapProducer::~HGCalPhotonIDValueMapProducer() {}
0086
0087 // ------------ method called to produce the data  ------------
0088 void HGCalPhotonIDValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0089   using namespace edm;
0090
0091   auto photonsH = iEvent.getHandle(photonsToken_);
0092
0093   // Clear previous map
0094   for (auto&& kv : maps_) {
0095     kv.second.clear();
0096     kv.second.reserve(photonsH->size());
0097   }
0098
0099   // Set up helper tool
0100   phoIDHelper_->eventInit(iEvent, iSetup);
0101
0102   for (const auto& pho : *photonsH) {
0103     if (pho.isEB()) {
0104       // Fill some dummy value
0105       for (auto&& kv : maps_) {
0106         kv.second.push_back(0.);
0107       }
0108     } else {
0109       phoIDHelper_->computeHGCAL(pho, radius_);
0110
0111       // check the PCA has worked out
0112       if (phoIDHelper_->sigmaUU() == -1) {
0113         for (auto&& kv : maps_) {
0114           kv.second.push_back(0.);
0115         }
0116         continue;
0117       }
0118
0119       hgcal::LongDeps ld(phoIDHelper_->energyPerLayer(radius_, true));
0120       float measuredDepth, expectedDepth, expectedSigma;
0121       float depthCompatibility =
0122           phoIDHelper_->clusterDepthCompatibility(ld, measuredDepth, expectedDepth, expectedSigma);
0123
0124       // Fill here all the ValueMaps from their appropriate functions
0125
0126       // energies calculated in an cylinder around the axis of the pho cluster
0127       float seed_tot_energy = ld.energyEE() + ld.energyFH() + ld.energyBH();
0128       const double div_cosh_eta =
0129           pho.superCluster()->seed()->position().rho() / pho.superCluster()->seed()->position().r();
0130       maps_["seedEt"].push_back(seed_tot_energy * div_cosh_eta);
0131       maps_["seedEnergy"].push_back(seed_tot_energy);
0132       maps_["seedEnergyEE"].push_back(ld.energyEE());
0133       maps_["seedEnergyFH"].push_back(ld.energyFH());
0134       maps_["seedEnergyBH"].push_back(ld.energyBH());
0135
0136       // Cluster shapes
0137       // PCA related
0138       maps_["pcaEig1"].push_back(phoIDHelper_->eigenValues()(0));
0139       maps_["pcaEig2"].push_back(phoIDHelper_->eigenValues()(1));
0140       maps_["pcaEig3"].push_back(phoIDHelper_->eigenValues()(2));
0141       maps_["pcaSig1"].push_back(phoIDHelper_->sigmas()(0));
0142       maps_["pcaSig2"].push_back(phoIDHelper_->sigmas()(1));
0143       maps_["pcaSig3"].push_back(phoIDHelper_->sigmas()(2));
0144
0145       // transverse shapes
0146       maps_["sigmaUU"].push_back(phoIDHelper_->sigmaUU());
0147       maps_["sigmaVV"].push_back(phoIDHelper_->sigmaVV());
0148       maps_["sigmaEE"].push_back(phoIDHelper_->sigmaEE());
0149       maps_["sigmaPP"].push_back(phoIDHelper_->sigmaPP());
0150
0151       // long profile
0152       maps_["nLayers"].push_back(ld.nLayers());
0153       maps_["firstLayer"].push_back(ld.firstLayer());
0154       maps_["lastLayer"].push_back(ld.lastLayer());
0155       maps_["e4oEtot"].push_back(ld.e4oEtot());
0156       maps_["layerEfrac10"].push_back(ld.layerEfrac10());
0157       maps_["layerEfrac90"].push_back(ld.layerEfrac90());
0158
0159       // depth
0160       maps_["measuredDepth"].push_back(measuredDepth);
0161       maps_["expectedDepth"].push_back(expectedDepth);
0162       maps_["expectedSigma"].push_back(expectedSigma);
0163       maps_["depthCompatibility"].push_back(depthCompatibility);
0164
0165       // Isolation
0166       maps_["caloIsoRing0"].push_back(phoIDHelper_->getIsolationRing(0));
0167       maps_["caloIsoRing1"].push_back(phoIDHelper_->getIsolationRing(1));
0168       maps_["caloIsoRing2"].push_back(phoIDHelper_->getIsolationRing(2));
0169       maps_["caloIsoRing3"].push_back(phoIDHelper_->getIsolationRing(3));
0170       maps_["caloIsoRing4"].push_back(phoIDHelper_->getIsolationRing(4));
0171     }
0172   }
0173
0174   // Check we didn't make up a new variable and forget it in the constructor
0175   // (or some other pathology)
0176   if (maps_.size() != valuesProduced_.size()) {
0177     throw cms::Exception("HGCalPhotonIDValueMapProducer")
0178         << "We have a miscoded value map producer, since map size changed";
0179   }
0180
0181   for (auto&& kv : maps_) {
0182     // Check we didn't forget any values
0183     if (kv.second.size() != photonsH->size()) {
0184       throw cms::Exception("HGCalPhotonIDValueMapProducer")
0185           << "We have a miscoded value map producer, since the variable " << kv.first << " wasn't filled.";
0186     }
0187     // Do the filling
0188     auto out = std::make_unique<edm::ValueMap<float>>();
0189     edm::ValueMap<float>::Filler filler(*out);
0190     filler.insert(photonsH, kv.second.begin(), kv.second.end());
0191     filler.fill();
0192     // and put it into the event
0193     iEvent.put(std::move(out), kv.first);
0194   }
0195 }
0196
0197 // ------------ method called once each stream before processing any runs, lumis or events  ------------
0198 void HGCalPhotonIDValueMapProducer::beginStream(edm::StreamID) {}
0199
0200 // ------------ method called once each stream after processing all runs, lumis and events  ------------
0201 void HGCalPhotonIDValueMapProducer::endStream() {}
0202
0203 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0204 void HGCalPhotonIDValueMapProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0205   // hgcalPhotonIDValueMap
0206   edm::ParameterSetDescription desc;
0207   desc.add<edm::InputTag>("photons", edm::InputTag("photonsHGC"));
0208   desc.add<double>("pcaRadius", 3.0);
0209   desc.add<std::vector<std::string>>("variables", valuesProduced_);
0210   desc.add<std::vector<double>>("dEdXWeights")
0211       ->setComment("This must be copied from dEdX_weights in RecoLocalCalo.HGCalRecProducers.HGCalRecHit_cfi");
0212   desc.add<unsigned int>("isoNRings", 5);
0213   desc.add<double>("isoDeltaR", 0.15);
0214   desc.add<double>("isoDeltaRmin", 0.0);
0215   desc.add<edm::InputTag>("EERecHits", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0216   desc.add<edm::InputTag>("FHRecHits", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
0217   desc.add<edm::InputTag>("BHRecHits", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
0218   desc.add<edm::InputTag>("hitMapTag", edm::InputTag("recHitMapProducer", "hgcalRecHitMap"));
0219   descriptions.add("hgcalPhotonIDValueMap", desc);
0220 }
0221
0222 //define this as a plug-in
0223 DEFINE_FWK_MODULE(HGCalPhotonIDValueMapProducer);