EgammaTools/plugins/HGCalElectronIDValueMapProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:   RecoEgamma/EgammaTools
0004 // Class:    HGCalElectronIDValueMapProducer
0005 //
0006 /**\class HGCalElectronIDValueMapProducer HGCalElectronIDValueMapProducer.cc RecoEgamma/EgammaTools/plugins/HGCalElectronIDValueMapProducer.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/GsfElectron.h"
0035 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0036 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0037
0038 #include "RecoEgamma/EgammaTools/interface/HGCalEgammaIDHelper.h"
0039 #include "RecoEgamma/EgammaTools/interface/LongDeps.h"
0040
0041 class HGCalElectronIDValueMapProducer : public edm::stream::EDProducer<> {
0042 public:
0043   explicit HGCalElectronIDValueMapProducer(const edm::ParameterSet&);
0044   ~HGCalElectronIDValueMapProducer() override;
0045
0046   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0047
0048 private:
0049   void beginStream(edm::StreamID) override;
0050   void produce(edm::Event&, const edm::EventSetup&) override;
0051   void endStream() override;
0052
0053   // ----------member data ---------------------------
0054   edm::EDGetTokenT<edm::View<reco::GsfElectron>> electronsToken_;
0055   float radius_;
0056
0057   static const std::vector<std::string> valuesProduced_;
0058   std::map<const std::string, std::vector<float>> maps_;
0059
0060   std::unique_ptr<HGCalEgammaIDHelper> eIDHelper_;
0061 };
0062
0063 // All the ValueMap names to output are defined in the auto-generated python cfi
0064 // so that potential consumers can configure themselves in a simple manner
0065 // Would be cool to use compile-time validation, but need constexpr strings, e.g. std::string_view in C++17
0066 const std::vector<std::string> HGCalElectronIDValueMapProducer::valuesProduced_ = {
0067     "ecOrigEt",      "ecOrigEnergy",  "ecEt",
0068     "ecEnergy",      "ecEnergyEE",    "ecEnergyFH",
0069     "ecEnergyBH",    "pcaEig1",       "pcaEig2",
0070     "pcaEig3",       "pcaSig1",       "pcaSig2",
0071     "pcaSig3",       "pcaAxisX",      "pcaAxisY",
0072     "pcaAxisZ",      "pcaPositionX",  "pcaPositionY",
0073     "pcaPositionZ",  "sigmaUU",       "sigmaVV",
0074     "sigmaEE",       "sigmaPP",       "nLayers",
0075     "firstLayer",    "lastLayer",     "e4oEtot",
0076     "layerEfrac10",  "layerEfrac90",  "measuredDepth",
0077     "expectedDepth", "expectedSigma", "depthCompatibility",
0078     "caloIsoRing0",  "caloIsoRing1",  "caloIsoRing2",
0079     "caloIsoRing3",  "caloIsoRing4",
0080 };
0081
0082 HGCalElectronIDValueMapProducer::HGCalElectronIDValueMapProducer(const edm::ParameterSet& iConfig)
0083     : electronsToken_(consumes(iConfig.getParameter<edm::InputTag>("electrons"))),
0084       radius_(iConfig.getParameter<double>("pcaRadius")) {
0085   for (const auto& key : valuesProduced_) {
0086     maps_[key] = {};
0087     produces<edm::ValueMap<float>>(key);
0088   }
0089
0090   eIDHelper_ = std::make_unique<HGCalEgammaIDHelper>(iConfig, consumesCollector());
0091 }
0092
0093 HGCalElectronIDValueMapProducer::~HGCalElectronIDValueMapProducer() {}
0094
0095 // ------------ method called to produce the data  ------------
0096 void HGCalElectronIDValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0097   using namespace edm;
0098
0099   auto electronsH = iEvent.getHandle(electronsToken_);
0100
0101   // Clear previous map
0102   for (auto&& kv : maps_) {
0103     kv.second.clear();
0104     kv.second.reserve(electronsH->size());
0105   }
0106
0107   // Set up helper tool
0108   eIDHelper_->eventInit(iEvent, iSetup);
0109
0110   for (const auto& electron : *electronsH) {
0111     if (electron.isEB()) {
0112       // Fill some dummy value
0113       for (auto&& kv : maps_) {
0114         kv.second.push_back(0.);
0115       }
0116     } else {
0117       eIDHelper_->computeHGCAL(electron, radius_);
0118
0119       // check the PCA has worked out
0120       if (eIDHelper_->sigmaUU() == -1) {
0121         for (auto&& kv : maps_) {
0122           kv.second.push_back(0.);
0123         }
0124         continue;
0125       }
0126
0127       hgcal::LongDeps ld(eIDHelper_->energyPerLayer(radius_, true));
0128       float measuredDepth, expectedDepth, expectedSigma;
0129       float depthCompatibility = eIDHelper_->clusterDepthCompatibility(ld, measuredDepth, expectedDepth, expectedSigma);
0130
0131       // Fill here all the ValueMaps from their appropriate functions
0132
0133       // Energies / PT
0134       const auto* eleCluster = electron.electronCluster().get();
0135       const double sinTheta = eleCluster->position().rho() / eleCluster->position().r();
0136       maps_["ecOrigEt"].push_back(eleCluster->energy() * sinTheta);
0137       maps_["ecOrigEnergy"].push_back(eleCluster->energy());
0138
0139       // energies calculated in an cylinder around the axis of the electron cluster
0140       float ec_tot_energy = ld.energyEE() + ld.energyFH() + ld.energyBH();
0141       maps_["ecEt"].push_back(ec_tot_energy * sinTheta);
0142       maps_["ecEnergy"].push_back(ec_tot_energy);
0143       maps_["ecEnergyEE"].push_back(ld.energyEE());
0144       maps_["ecEnergyFH"].push_back(ld.energyFH());
0145       maps_["ecEnergyBH"].push_back(ld.energyBH());
0146
0147       // Cluster shapes
0148       // PCA related
0149       maps_["pcaEig1"].push_back(eIDHelper_->eigenValues()(0));
0150       maps_["pcaEig2"].push_back(eIDHelper_->eigenValues()(1));
0151       maps_["pcaEig3"].push_back(eIDHelper_->eigenValues()(2));
0152       maps_["pcaSig1"].push_back(eIDHelper_->sigmas()(0));
0153       maps_["pcaSig2"].push_back(eIDHelper_->sigmas()(1));
0154       maps_["pcaSig3"].push_back(eIDHelper_->sigmas()(2));
0155       maps_["pcaAxisX"].push_back(eIDHelper_->axis().x());
0156       maps_["pcaAxisY"].push_back(eIDHelper_->axis().y());
0157       maps_["pcaAxisZ"].push_back(eIDHelper_->axis().z());
0158       maps_["pcaPositionX"].push_back(eIDHelper_->barycenter().x());
0159       maps_["pcaPositionY"].push_back(eIDHelper_->barycenter().y());
0160       maps_["pcaPositionZ"].push_back(eIDHelper_->barycenter().z());
0161
0162       // transverse shapes
0163       maps_["sigmaUU"].push_back(eIDHelper_->sigmaUU());
0164       maps_["sigmaVV"].push_back(eIDHelper_->sigmaVV());
0165       maps_["sigmaEE"].push_back(eIDHelper_->sigmaEE());
0166       maps_["sigmaPP"].push_back(eIDHelper_->sigmaPP());
0167
0168       // long profile
0169       maps_["nLayers"].push_back(ld.nLayers());
0170       maps_["firstLayer"].push_back(ld.firstLayer());
0171       maps_["lastLayer"].push_back(ld.lastLayer());
0172       maps_["e4oEtot"].push_back(ld.e4oEtot());
0173       maps_["layerEfrac10"].push_back(ld.layerEfrac10());
0174       maps_["layerEfrac90"].push_back(ld.layerEfrac90());
0175
0176       // depth
0177       maps_["measuredDepth"].push_back(measuredDepth);
0178       maps_["expectedDepth"].push_back(expectedDepth);
0179       maps_["expectedSigma"].push_back(expectedSigma);
0180       maps_["depthCompatibility"].push_back(depthCompatibility);
0181
0182       // Isolation
0183       maps_["caloIsoRing0"].push_back(eIDHelper_->getIsolationRing(0));
0184       maps_["caloIsoRing1"].push_back(eIDHelper_->getIsolationRing(1));
0185       maps_["caloIsoRing2"].push_back(eIDHelper_->getIsolationRing(2));
0186       maps_["caloIsoRing3"].push_back(eIDHelper_->getIsolationRing(3));
0187       maps_["caloIsoRing4"].push_back(eIDHelper_->getIsolationRing(4));
0188     }
0189   }
0190
0191   // Check we didn't make up a new variable and forget it in valuesProduced_
0192   if (maps_.size() != valuesProduced_.size()) {
0193     throw cms::Exception("HGCalElectronIDValueMapProducer")
0194         << "We have a miscoded value map producer, since map size changed";
0195   }
0196
0197   for (auto&& kv : maps_) {
0198     // Check we didn't forget any values
0199     if (kv.second.size() != electronsH->size()) {
0200       throw cms::Exception("HGCalElectronIDValueMapProducer")
0201           << "We have a miscoded value map producer, since the variable " << kv.first << " wasn't filled.";
0202     }
0203     // Do the filling
0204     auto out = std::make_unique<edm::ValueMap<float>>();
0205     edm::ValueMap<float>::Filler filler(*out);
0206     filler.insert(electronsH, kv.second.begin(), kv.second.end());
0207     filler.fill();
0208     // and put it into the event
0209     iEvent.put(std::move(out), kv.first);
0210   }
0211 }
0212
0213 // ------------ method called once each stream before processing any runs, lumis or events  ------------
0214 void HGCalElectronIDValueMapProducer::beginStream(edm::StreamID) {}
0215
0216 // ------------ method called once each stream after processing all runs, lumis and events  ------------
0217 void HGCalElectronIDValueMapProducer::endStream() {}
0218
0219 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0220 void HGCalElectronIDValueMapProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0221   // Auto-generate hgcalElectronIDValueMap_cfi
0222   edm::ParameterSetDescription desc;
0223   desc.add<edm::InputTag>("electrons", edm::InputTag("ecalDrivenGsfElectronsHGC"));
0224   desc.add<double>("pcaRadius", 3.0);
0225   desc.add<std::vector<std::string>>("variables", valuesProduced_);
0226   desc.add<std::vector<double>>("dEdXWeights")
0227       ->setComment("This must be copied from dEdX_weights in RecoLocalCalo.HGCalRecProducers.HGCalRecHit_cfi");
0228   desc.add<unsigned int>("isoNRings", 5);
0229   desc.add<double>("isoDeltaR", 0.15);
0230   desc.add<double>("isoDeltaRmin", 0.0);
0231   desc.add<edm::InputTag>("EERecHits", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0232   desc.add<edm::InputTag>("FHRecHits", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
0233   desc.add<edm::InputTag>("BHRecHits", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
0234   desc.add<edm::InputTag>("hitMapTag", edm::InputTag("recHitMapProducer", "hgcalRecHitMap"));
0235   descriptions.add("hgcalElectronIDValueMap", desc);
0236 }
0237
0238 //define this as a plug-in
0239 DEFINE_FWK_MODULE(HGCalElectronIDValueMapProducer);