File indexing completed on 2023-03-17 11:17:54
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021
0022
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
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
0063
0064
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
0088 void HGCalPhotonIDValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0089 using namespace edm;
0090
0091 auto photonsH = iEvent.getHandle(photonsToken_);
0092
0093
0094 for (auto&& kv : maps_) {
0095 kv.second.clear();
0096 kv.second.reserve(photonsH->size());
0097 }
0098
0099
0100 phoIDHelper_->eventInit(iEvent, iSetup);
0101
0102 for (const auto& pho : *photonsH) {
0103 if (pho.isEB()) {
0104
0105 for (auto&& kv : maps_) {
0106 kv.second.push_back(0.);
0107 }
0108 } else {
0109 phoIDHelper_->computeHGCAL(pho, radius_);
0110
0111
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
0125
0126
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
0137
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
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
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
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
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
0175
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
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
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
0193 iEvent.put(std::move(out), kv.first);
0194 }
0195 }
0196
0197
0198 void HGCalPhotonIDValueMapProducer::beginStream(edm::StreamID) {}
0199
0200
0201 void HGCalPhotonIDValueMapProducer::endStream() {}
0202
0203
0204 void HGCalPhotonIDValueMapProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0205
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("hgcalRecHitMapProducer"));
0219 descriptions.add("hgcalPhotonIDValueMap", desc);
0220 }
0221
0222
0223 DEFINE_FWK_MODULE(HGCalPhotonIDValueMapProducer);