From 9d00a354a0cfeae450f9210b71ff08a193e621a3 Mon Sep 17 00:00:00 2001 From: Badder Date: Thu, 27 Jan 2022 14:53:31 +0100 Subject: [PATCH] DeepSC graph evaluation move to 12_3_X --- .../interface/PFECALSuperClusterAlgo.h | 46 +- .../src/PFECALSuperClusterAlgo.cc | 243 +++++--- .../particleFlowSuperClusterECAL_cff.py | 1 + .../particleFlowSuperClusterECAL_cfi.py | 4 +- ...particleFlowSuperClusteringSequence_cff.py | 3 +- .../src/PFECALSuperClusterProducer.cc | 209 +++++-- RecoEcal/EgammaCoreTools/BuildFile.xml | 6 + .../interface/CalibratedPFCluster.h | 24 + .../interface/DeepSCGraphEvaluation.h | 73 +++ .../interface/EcalClustersGraph.h | 135 +++++ .../EgammaCoreTools/interface/GraphMatrix.h | 546 ++++++++++++++++++ .../src/DeepSCGraphEvaluation.cc | 213 +++++++ .../EgammaCoreTools/src/EcalClustersGraph.cc | 526 +++++++++++++++++ 13 files changed, 1892 insertions(+), 137 deletions(-) create mode 100644 RecoEcal/EgammaCoreTools/interface/CalibratedPFCluster.h create mode 100644 RecoEcal/EgammaCoreTools/interface/DeepSCGraphEvaluation.h create mode 100644 RecoEcal/EgammaCoreTools/interface/EcalClustersGraph.h create mode 100644 RecoEcal/EgammaCoreTools/interface/GraphMatrix.h create mode 100644 RecoEcal/EgammaCoreTools/src/DeepSCGraphEvaluation.cc create mode 100644 RecoEcal/EgammaCoreTools/src/EcalClustersGraph.cc diff --git a/RecoEcal/EgammaClusterAlgos/interface/PFECALSuperClusterAlgo.h b/RecoEcal/EgammaClusterAlgos/interface/PFECALSuperClusterAlgo.h index 4e3558e89fe90..c15905139f66c 100644 --- a/RecoEcal/EgammaClusterAlgos/interface/PFECALSuperClusterAlgo.h +++ b/RecoEcal/EgammaClusterAlgos/interface/PFECALSuperClusterAlgo.h @@ -13,11 +13,23 @@ #include "DataFormats/EgammaReco/interface/BasicCluster.h" #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h" +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "DataFormats/EcalDetId/interface/EEDetId.h" +#include "DataFormats/CaloRecHit/interface/CaloCluster.h" #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h" +#include "Geometry/CaloTopology/interface/CaloTopology.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" +#include "Geometry/Records/interface/CaloTopologyRecord.h" +#include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h" + #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h" #include "RecoEcal/EgammaClusterAlgos/interface/SCEnergyCorrectorSemiParm.h" +#include "RecoEcal/EgammaCoreTools/interface/GraphMatrix.h" +#include "RecoEcal/EgammaCoreTools/interface/EcalClustersGraph.h" +#include "RecoEcal/EgammaCoreTools/interface/CalibratedPFCluster.h" #include "FWCore/Framework/interface/Event.h" #include "FWCore/Framework/interface/ESHandle.h" @@ -37,6 +49,8 @@ #include "CondFormats/EcalObjects/interface/EcalSCDynamicDPhiParameters.h" #include "CondFormats/DataRecord/interface/EcalSCDynamicDPhiParametersRcd.h" +#include "RecoEcal/EgammaCoreTools/interface/DeepSCGraphEvaluation.h" + #include #include @@ -48,31 +62,18 @@ \date July 2012 */ +// class SCProducerCache; + class PFECALSuperClusterAlgo { public: - enum clustering_type { kBOX = 1, kMustache = 2 }; + enum clustering_type { kBOX = 1, kMustache = 2, kDeepSC = 3 }; enum energy_weight { kRaw, kCalibratedNoPS, kCalibratedTotal }; - // simple class for associating calibrated energies - class CalibratedPFCluster { - public: - CalibratedPFCluster(const edm::Ptr& p) : cluptr(p) {} - - double energy() const { return cluptr->correctedEnergy(); } - double energy_nocalib() const { return cluptr->energy(); } - double eta() const { return cluptr->positionREP().eta(); } - double phi() const { return cluptr->positionREP().phi(); } - - edm::Ptr the_ptr() const { return cluptr; } - - private: - edm::Ptr cluptr; - }; typedef std::shared_ptr CalibratedClusterPtr; typedef std::vector CalibratedClusterPtrVector; /// constructor - PFECALSuperClusterAlgo(); + PFECALSuperClusterAlgo(const SCProducerCache* cache); void setVerbosityLevel(bool verbose) { verbose_ = verbose; } @@ -110,6 +111,9 @@ class PFECALSuperClusterAlgo { void setCrackCorrections(bool applyCrackCorrections) { applyCrackCorrections_ = applyCrackCorrections; } void setTokens(const edm::ParameterSet&, edm::ConsumesCollector&&); + + void setTensorflowObjects(); + void update(const edm::EventSetup&); void updateSCParams(const edm::EventSetup&); @@ -132,6 +136,12 @@ class PFECALSuperClusterAlgo { const reco::BeamSpot* beamSpot_; const ESChannelStatus* channelStatus_; + const CaloGeometry* geometry_; + const CaloSubdetectorGeometry* ebGeom_; + const CaloSubdetectorGeometry* eeGeom_; + const CaloSubdetectorGeometry* esGeom_; + const CaloTopology* topology_; + const EcalMustacheSCParameters* mustacheSCParams_; const EcalSCDynamicDPhiParameters* scDynamicDPhiParams_; @@ -173,6 +183,8 @@ class PFECALSuperClusterAlgo { bool applyCrackCorrections_; bool threshIsET_; + const reco::SCProducerCache* SCProducerCache_; + // OOT photons bool isOOTCollection_; edm::EDGetTokenT inputTagBarrelRecHits_; diff --git a/RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc b/RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc index 92b0ee63e7f6d..2c99090de7796 100644 --- a/RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc +++ b/RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc @@ -1,8 +1,11 @@ #include "RecoEcal/EgammaClusterAlgos/interface/PFECALSuperClusterAlgo.h" +#include "RecoEcal/EgammaCoreTools/interface/EcalClustersGraph.h" #include "CommonTools/ParticleFlow/interface/PFClusterWidthAlgo.h" #include "RecoParticleFlow/PFClusterTools/interface/LinkByRecHit.h" #include "DataFormats/ParticleFlowReco/interface/PFLayer.h" #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h" +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "DataFormats/EcalDetId/interface/EEDetId.h" #include "DataFormats/EcalDetId/interface/ESDetId.h" #include "DataFormats/HcalDetId/interface/HcalDetId.h" #include "RecoEcal/EgammaCoreTools/interface/Mustache.h" @@ -26,8 +29,8 @@ namespace { typedef edm::View PFClusterView; typedef edm::Ptr PFClusterPtr; typedef edm::PtrVector PFClusterPtrVector; - typedef PFECALSuperClusterAlgo::CalibratedClusterPtr CalibClusterPtr; - typedef PFECALSuperClusterAlgo::CalibratedClusterPtrVector CalibClusterPtrVector; + typedef std::shared_ptr CalibClusterPtr; + typedef std::vector CalibClusterPtrVector; typedef std::pair EEPSPair; bool sortByKey(const EEPSPair& a, const EEPSPair& b) { return a.first < b.first; } @@ -81,6 +84,52 @@ namespace { return x_rechits_match / x_rechits_tot > majority; } + std::vector clusterLocalPosition(const CalibClusterPtr& cluster, + const CaloSubdetectorGeometry* ebGeom_, + const CaloSubdetectorGeometry* eeGeom_) { + std::vector position; // ieta,iphi,iz or ix,iy,iz + position.resize(3); + reco::CaloCluster caloBC(*cluster->the_ptr()); + math::XYZPoint caloPos = caloBC.position(); + if (cluster->the_ptr()->layer() == PFLayer::ECAL_BARREL) { + EBDetId id(ebGeom_->getClosestCell(GlobalPoint(caloPos.x(), caloPos.y(), caloPos.z()))); + position[0] = id.ieta(); + position[1] = id.ieta(); + position[2] = 0; + } else if (cluster->the_ptr()->layer() == PFLayer::ECAL_ENDCAP) { + EEDetId id(eeGeom_->getClosestCell(GlobalPoint(caloPos.x(), caloPos.y(), caloPos.z()))); + position[0] = id.ix(); + position[1] = id.iy(); + position[2] = id.zside(); + } + return position; + } + + DetId clusterDetId(const CalibClusterPtr& cluster, + const CaloSubdetectorGeometry* ebGeom_, + const CaloSubdetectorGeometry* eeGeom_) { + DetId clId; + reco::CaloCluster caloBC(*cluster->the_ptr()); + math::XYZPoint caloPos = caloBC.position(); + if (cluster->the_ptr()->layer() == PFLayer::ECAL_BARREL) { + EBDetId id(ebGeom_->getClosestCell(GlobalPoint(caloPos.x(), caloPos.y(), caloPos.z()))); + clId = id; + } else if (cluster->the_ptr()->layer() == PFLayer::ECAL_ENDCAP) { + EEDetId id(eeGeom_->getClosestCell(GlobalPoint(caloPos.x(), caloPos.y(), caloPos.z()))); + clId = id; + } + return clId; + } + + double clusterZside(const CalibClusterPtr& cluster) { + double zSide = 0.; + if (cluster->the_ptr()->layer() == PFLayer::ECAL_ENDCAP && cluster->eta() < 0.) + zSide = -1.; + if (cluster->the_ptr()->layer() == PFLayer::ECAL_ENDCAP && cluster->eta() > 0.) + zSide = +1.; + return zSide; + } + bool isClustered(const CalibClusterPtr& x, const CalibClusterPtr seed, const PFECALSuperClusterAlgo::clustering_type type, @@ -107,14 +156,14 @@ namespace { } // namespace -PFECALSuperClusterAlgo::PFECALSuperClusterAlgo() : beamSpot_(nullptr) {} +PFECALSuperClusterAlgo::PFECALSuperClusterAlgo(const reco::SCProducerCache* cache) : beamSpot_(nullptr),SCProducerCache_(cache) {} void PFECALSuperClusterAlgo::setPFClusterCalibration(const std::shared_ptr& calib) { _pfEnergyCalibration = calib; } void PFECALSuperClusterAlgo::setTokens(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& cc) { - inputTagPFClusters_ = cc.consumes >(iConfig.getParameter("PFClusters")); + inputTagPFClusters_ = cc.consumes>(iConfig.getParameter("PFClusters")); inputTagPFClustersES_ = cc.consumes(iConfig.getParameter("ESAssociation")); inputTagBeamSpot_ = cc.consumes(iConfig.getParameter("BeamSpot")); @@ -137,7 +186,8 @@ void PFECALSuperClusterAlgo::setTokens(const edm::ParameterSet& iConfig, edm::Co regr_->setTokens(regconf, cc); } - if (isOOTCollection_) { // OOT photons only + if (isOOTCollection_ || _clustype == PFECALSuperClusterAlgo::kDeepSC) { // OOT photons or DeepSC + //std::cout << "_clustype:" << _clustype << std::endl; inputTagBarrelRecHits_ = cc.consumes(iConfig.getParameter("barrelRecHits")); inputTagEndcapRecHits_ = cc.consumes(iConfig.getParameter("endcapRecHits")); } @@ -153,6 +203,17 @@ void PFECALSuperClusterAlgo::update(const edm::EventSetup& setup) { edm::ESHandle esChannelStatusHandle_ = setup.getHandle(esChannelStatusToken_); channelStatus_ = esChannelStatusHandle_.product(); + + edm::ESHandle caloGeometryHandle_; + setup.get().get(caloGeometryHandle_); + geometry_ = caloGeometryHandle_.product(); + ebGeom_ = caloGeometryHandle_->getSubdetectorGeometry(DetId::Ecal, EcalBarrel); + eeGeom_ = caloGeometryHandle_->getSubdetectorGeometry(DetId::Ecal, EcalEndcap); + esGeom_ = caloGeometryHandle_->getSubdetectorGeometry(DetId::Ecal, EcalPreshower); + + edm::ESHandle caloTopologyHandle_; + setup.get().get(caloTopologyHandle_); + topology_ = caloTopologyHandle_.product(); } void PFECALSuperClusterAlgo::updateSCParams(const edm::EventSetup& setup) { @@ -167,7 +228,7 @@ void PFECALSuperClusterAlgo::updateSCParams(const edm::EventSetup& setup) { void PFECALSuperClusterAlgo::loadAndSortPFClusters(const edm::Event& iEvent) { //load input collections //Load the pfcluster collections - edm::Handle > pfclustersHandle; + edm::Handle> pfclustersHandle; iEvent.getByToken(inputTagPFClusters_, pfclustersHandle); edm::Handle psAssociationHandle; @@ -196,7 +257,7 @@ void PFECALSuperClusterAlgo::loadAndSortPFClusters(const edm::Event& iEvent) { //Select PF clusters available for the clustering for (size_t i = 0; i < clusters.size(); ++i) { auto cluster = clusters.ptrAt(i); - LogDebug("PFClustering") << "Loading PFCluster i=" << cluster.key() << " energy=" << cluster->energy() << std::endl; + //LogDebug("PFClustering") << "Loading PFCluster i=" << cluster.key() << " energy=" << cluster->energy() << std::endl; // protection for sim clusters if (cluster->caloID().detectors() == 0 && cluster->hitsAndFractions().empty()) @@ -224,13 +285,13 @@ void PFECALSuperClusterAlgo::loadAndSortPFClusters(const edm::Event& iEvent) { std::sort(_clustersEB.begin(), _clustersEB.end(), greaterByEt); std::sort(_clustersEE.begin(), _clustersEE.end(), greaterByEt); - // set recHit collections for OOT photons - if (isOOTCollection_) { + // set recHit collections for OOT photons and DeepSC + if (isOOTCollection_ || _clustype == PFECALSuperClusterAlgo::kDeepSC) { edm::Handle barrelRecHitsHandle; iEvent.getByToken(inputTagBarrelRecHits_, barrelRecHitsHandle); if (!barrelRecHitsHandle.isValid()) { throw cms::Exception("PFECALSuperClusterAlgo") - << "If you use OOT photons, need to specify proper barrel rec hit collection"; + << "If you use OOT photons or DeepSC, need to specify proper barrel rec hit collection"; } barrelRecHits_ = barrelRecHitsHandle.product(); @@ -238,7 +299,7 @@ void PFECALSuperClusterAlgo::loadAndSortPFClusters(const edm::Event& iEvent) { iEvent.getByToken(inputTagEndcapRecHits_, endcapRecHitsHandle); if (!endcapRecHitsHandle.isValid()) { throw cms::Exception("PFECALSuperClusterAlgo") - << "If you use OOT photons, need to specify proper endcap rec hit collection"; + << "If you use OOT photons or DeepSC, need to specify proper endcap rec hit collection"; } endcapRecHits_ = endcapRecHitsHandle.product(); } @@ -253,18 +314,61 @@ void PFECALSuperClusterAlgo::run() { void PFECALSuperClusterAlgo::buildAllSuperClusters(CalibClusterPtrVector& clusters, double seedthresh) { auto seedable = std::bind(isSeed, _1, seedthresh, threshIsET_); - // make sure only seeds appear at the front of the list of clusters - std::stable_partition(clusters.begin(), clusters.end(), seedable); - // in each iteration we are working on a list that is already sorted - // in the cluster energy and remains so through each iteration - // NB: since clusters is sorted in loadClusters any_of has O(1) - // timing until you run out of seeds! - while (std::any_of(clusters.cbegin(), clusters.cend(), seedable)) { - buildSuperCluster(clusters.front(), clusters); + + if (_clustype != PFECALSuperClusterAlgo::kDeepSC) { + // make sure only seeds appear at the front of the list of clusters + std::stable_partition(clusters.begin(), clusters.end(), seedable); + + // in each iteration we are working on a list that is already sorted + // in the cluster energy and remains so through each iteration + // NB: since clusters is sorted in loadClusters any_of has O(1) + // timing until you run out of seeds! + while (std::any_of(clusters.cbegin(), clusters.cend(), seedable)) { + buildSuperCluster(clusters.front(), clusters); + } + + } else { + //TEST EcalClustersGraph + + // make sure only seeds appear at the front of the list of clusters + auto last_seed = std::stable_partition(clusters.begin(), clusters.end(), seedable); + + EcalClustersGraph ecalClusterGraph_ {clusters, + static_cast(std::distance(clusters.begin(), last_seed)), + topology_, + ebGeom_, + eeGeom_, + barrelRecHits_, + endcapRecHits_, + SCProducerCache_}; + + // check which clusters are inside the dynamic windows ('1') and which are out ('0') + ecalClusterGraph_.initWindows(); + + // for each pfCluster inside a window ('1'),i.e. a matrix row, fill the variables needed for evaluating the GraphNet + ecalClusterGraph_.fillVariables(); + + // for each window evaluate the GrpahNet score of any pfCluster inside the windwo ('1') + ecalClusterGraph_.evaluateScores(); + ecalClusterGraph_.printDebugInfo(); + + // first keep all pfClusters with a score greater than a threshold (seed-eta and seed-et dependent), + // then reduce elements and remove duplicates (pfClusters in many windows) + ecalClusterGraph_.setThresholds(); + ecalClusterGraph_.selectClusters(); + + // for each window make a superCluster out of the remaining pfClusters in the window ('1') + std::vector> windows = ecalClusterGraph_.getWindows(); + for (unsigned int iw = 0; iw < windows.size(); iw++) + buildSuperCluster(windows.at(iw).first, windows.at(iw).second); + + ecalClusterGraph_.clearWindows(); } } void PFECALSuperClusterAlgo::buildSuperCluster(CalibClusterPtr& seed, CalibClusterPtrVector& clusters) { + CalibratedClusterPtrVector clustered; + double etawidthSuperCluster = 0.0; double phiwidthSuperCluster = 0.0; bool isEE = false; @@ -286,59 +390,68 @@ void PFECALSuperClusterAlgo::buildSuperCluster(CalibClusterPtr& seed, CalibClust default: break; } - auto isClusteredWithSeed = std::bind(isClustered, - _1, - seed, - _clustype, - mustacheSCParams_, - scDynamicDPhiParams_, - useDynamicDPhi_, - etawidthSuperCluster, - phiwidthSuperCluster); - auto matchesSeedByRecHit = std::bind(isLinkedByRecHit, _1, seed, satelliteThreshold_, fractionForMajority_, 0.1, 0.2); - - // this function shuffles the list of clusters into a list - // where all clustered sub-clusters are at the front - // and returns a pointer to the first unclustered cluster. - // The relative ordering of clusters is preserved - // (i.e. both resulting sub-lists are sorted by energy). - auto not_clustered = std::stable_partition(clusters.begin(), clusters.end(), isClusteredWithSeed); - // satellite cluster merging - // it was found that large clusters can split! - if (doSatelliteClusterMerge_) { - not_clustered = std::stable_partition(not_clustered, clusters.end(), matchesSeedByRecHit); - } - if (verbose_) { - edm::LogInfo("PFClustering") << "Dumping cluster detail"; - edm::LogVerbatim("PFClustering") << "\tPassed seed: e = " << seed->energy_nocalib() << " eta = " << seed->eta() - << " phi = " << seed->phi() << std::endl; - for (auto clus = clusters.cbegin(); clus != not_clustered; ++clus) { - edm::LogVerbatim("PFClustering") << "\t\tClustered cluster: e = " << (*clus)->energy_nocalib() - << " eta = " << (*clus)->eta() << " phi = " << (*clus)->phi() << std::endl; + if (_clustype != PFECALSuperClusterAlgo::kDeepSC) { + auto isClusteredWithSeed = std::bind(isClustered, + _1, + seed, + _clustype, + mustacheSCParams_, + scDynamicDPhiParams_, + useDynamicDPhi_, + etawidthSuperCluster, + phiwidthSuperCluster); + + auto matchesSeedByRecHit = + std::bind(isLinkedByRecHit, _1, seed, satelliteThreshold_, fractionForMajority_, 0.1, 0.2); + + // this function shuffles the list of clusters into a list + // where all clustered sub-clusters are at the front + // and returns a pointer to the first unclustered cluster. + // The relative ordering of clusters is preserved + // (i.e. both resulting sub-lists are sorted by energy). + auto not_clustered = std::stable_partition(clusters.begin(), clusters.end(), isClusteredWithSeed); + // satellite cluster merging + // it was found that large clusters can split! + if (doSatelliteClusterMerge_) { + not_clustered = std::stable_partition(not_clustered, clusters.end(), matchesSeedByRecHit); } - for (auto clus = not_clustered; clus != clusters.end(); ++clus) { - edm::LogVerbatim("PFClustering") << "\tNon-Clustered cluster: e = " << (*clus)->energy_nocalib() - << " eta = " << (*clus)->eta() << " phi = " << (*clus)->phi() << std::endl; + + if (verbose_) { + edm::LogInfo("PFClustering") << "Dumping cluster detail"; + edm::LogVerbatim("PFClustering") << "\tPassed seed: e = " << seed->energy_nocalib() << " eta = " << seed->eta() + << " phi = " << seed->phi() << std::endl; + for (auto clus = clusters.cbegin(); clus != not_clustered; ++clus) { + edm::LogVerbatim("PFClustering") << "\t\tClustered cluster: e = " << (*clus)->energy_nocalib() + << " eta = " << (*clus)->eta() << " phi = " << (*clus)->phi() << std::endl; + } + for (auto clus = not_clustered; clus != clusters.end(); ++clus) { + edm::LogVerbatim("PFClustering") << "\tNon-Clustered cluster: e = " << (*clus)->energy_nocalib() + << " eta = " << (*clus)->eta() << " phi = " << (*clus)->phi() << std::endl; + } } - } - if (not_clustered == clusters.begin()) { - if (dropUnseedable_) { - clusters.erase(clusters.begin()); - return; - } else { - throw cms::Exception("PFECALSuperClusterAlgo::buildSuperCluster") - << "Cluster is not seedable!" << std::endl - << "\tNon-Clustered cluster: e = " << (*not_clustered)->energy_nocalib() - << " eta = " << (*not_clustered)->eta() << " phi = " << (*not_clustered)->phi() << std::endl; + if (not_clustered == clusters.begin()) { + if (dropUnseedable_) { + clusters.erase(clusters.begin()); + return; + } else { + throw cms::Exception("PFECALSuperClusterAlgo::buildSuperCluster") + << "Cluster is not seedable!" << std::endl + << "\tNon-Clustered cluster: e = " << (*not_clustered)->energy_nocalib() + << " eta = " << (*not_clustered)->eta() << " phi = " << (*not_clustered)->phi() << std::endl; + } } + + // move the clustered clusters out of available cluster list + // and into a temporary vector for building the SC + CalibratedClusterPtrVector clustered_tmp(clusters.begin(), not_clustered); + clustered = clustered_tmp; + clusters.erase(clusters.begin(), not_clustered); + } else { + clustered = clusters; } - // move the clustered clusters out of available cluster list - // and into a temporary vector for building the SC - CalibratedClusterPtrVector clustered(clusters.begin(), not_clustered); - clusters.erase(clusters.begin(), not_clustered); // need the vector of raw pointers for a PF width class std::vector bare_ptrs; // calculate necessary parameters and build the SC diff --git a/RecoEcal/EgammaClusterProducers/python/particleFlowSuperClusterECAL_cff.py b/RecoEcal/EgammaClusterProducers/python/particleFlowSuperClusterECAL_cff.py index 966c1add65c99..484d3c06aa405 100644 --- a/RecoEcal/EgammaClusterProducers/python/particleFlowSuperClusterECAL_cff.py +++ b/RecoEcal/EgammaClusterProducers/python/particleFlowSuperClusterECAL_cff.py @@ -1,6 +1,7 @@ import FWCore.ParameterSet.Config as cms from RecoEcal.EgammaClusterProducers.particleFlowSuperClusterECAL_cfi import * +from RecoEcal.EgammaClusterProducers.particleFlowSuperClusterECALDeepSC_cfi import * from RecoEcal.EgammaClusterProducers.particleFlowSuperClusterECALMustache_cfi import * from RecoEcal.EgammaClusterProducers.particleFlowSuperClusterECALBox_cfi import * from RecoEcal.EgammaClusterProducers.particleFlowSuperClusterECALOnly_cfi import * diff --git a/RecoEcal/EgammaClusterProducers/python/particleFlowSuperClusterECAL_cfi.py b/RecoEcal/EgammaClusterProducers/python/particleFlowSuperClusterECAL_cfi.py index 2c81cf48e0e4f..c23964e18fc04 100644 --- a/RecoEcal/EgammaClusterProducers/python/particleFlowSuperClusterECAL_cfi.py +++ b/RecoEcal/EgammaClusterProducers/python/particleFlowSuperClusterECAL_cfi.py @@ -1,9 +1,11 @@ import FWCore.ParameterSet.Config as cms from RecoEcal.EgammaClusterProducers.particleFlowSuperClusterECALMustache_cfi import particleFlowSuperClusterECALMustache as _particleFlowSuperClusterECALMustache +from RecoEcal.EgammaClusterProducers.particleFlowSuperClusterECALDeepSC_cfi import particleFlowSuperClusterECALDeepSC as _particleFlowSuperClusterECALDeepSC -# define the default ECAL clustering (Mustache or Box) +# define the default ECAL clustering (Mustache or Box or DeepSC) particleFlowSuperClusterECAL = _particleFlowSuperClusterECALMustache.clone() +particleFlowDeepSuperClusterECAL = _particleFlowSuperClusterECALDeepSC.clone() from Configuration.ProcessModifiers.pp_on_AA_cff import pp_on_AA pp_on_AA.toModify(particleFlowSuperClusterECAL, useDynamicDPhiWindow = False, diff --git a/RecoEcal/EgammaClusterProducers/python/particleFlowSuperClusteringSequence_cff.py b/RecoEcal/EgammaClusterProducers/python/particleFlowSuperClusteringSequence_cff.py index 4185150bcd19a..231b240bd3498 100644 --- a/RecoEcal/EgammaClusterProducers/python/particleFlowSuperClusteringSequence_cff.py +++ b/RecoEcal/EgammaClusterProducers/python/particleFlowSuperClusteringSequence_cff.py @@ -8,7 +8,8 @@ # Producer for energy corrections #from RecoEcal.EgammaClusterProducers.correctedDynamicHybridSuperClusters_cfi import * # PFECAL super clusters, either hybrid-clustering clone (Box) or mustache. -particleFlowSuperClusteringTask = cms.Task(particleFlowSuperClusterECAL) +#particleFlowSuperClusteringTask = cms.Task(particleFlowSuperClusterECAL) +particleFlowSuperClusteringTask = cms.Task(cms.Task(particleFlowSuperClusterECAL),cms.Task(particleFlowDeepSuperClusterECAL)) particleFlowSuperClusteringSequence = cms.Sequence(particleFlowSuperClusteringTask) particleFlowSuperClusterHGCal = particleFlowSuperClusterECAL.clone() diff --git a/RecoEcal/EgammaClusterProducers/src/PFECALSuperClusterProducer.cc b/RecoEcal/EgammaClusterProducers/src/PFECALSuperClusterProducer.cc index 7dfc9045f177a..32fc2c5bba615 100644 --- a/RecoEcal/EgammaClusterProducers/src/PFECALSuperClusterProducer.cc +++ b/RecoEcal/EgammaClusterProducers/src/PFECALSuperClusterProducer.cc @@ -1,9 +1,3 @@ -/**\class PFECALSuperClusterProducer - -\author Nicolas Chanon -Additional authors for Mustache: Y. Gershtein, R. Patel, L. Gray -\date July 2012 -*/ #include "CondFormats/DataRecord/interface/GBRWrapperRcd.h" #include "CondFormats/GBRForest/interface/GBRForest.h" @@ -35,20 +29,26 @@ Additional authors for Mustache: Y. Gershtein, R. Patel, L. Gray #include "RecoEcal/EgammaClusterAlgos/interface/PFECALSuperClusterAlgo.h" #include "RecoEcal/EgammaClusterAlgos/interface/SCEnergyCorrectorSemiParm.h" #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h" - +#include "RecoEcal/EgammaCoreTools/interface/DeepSCGraphEvaluation.h" #include "TVector2.h" #include #include -class PFECALSuperClusterProducer : public edm::stream::EDProducer<> { +class PFECALSuperClusterProducer : public edm::stream::EDProducer> { public: - explicit PFECALSuperClusterProducer(const edm::ParameterSet&); + explicit PFECALSuperClusterProducer(const edm::ParameterSet&, const reco::SCProducerCache* gcache); ~PFECALSuperClusterProducer() override; void beginLuminosityBlock(const edm::LuminosityBlock&, const edm::EventSetup&) override; void produce(edm::Event&, const edm::EventSetup&) override; + static std::unique_ptr initializeGlobalCache(const edm::ParameterSet& config){ + return std::make_unique(config); + } + + static void globalEndJob(const reco::SCProducerCache*){}; + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); private: @@ -86,13 +86,15 @@ using namespace edm; namespace { const std::string ClusterType__BOX("Box"); const std::string ClusterType__Mustache("Mustache"); + const std::string ClusterType__DeepSC("DeepSC"); const std::string EnergyWeight__Raw("Raw"); const std::string EnergyWeight__CalibratedNoPS("CalibratedNoPS"); const std::string EnergyWeight__CalibratedTotal("CalibratedTotal"); } // namespace -PFECALSuperClusterProducer::PFECALSuperClusterProducer(const edm::ParameterSet& iConfig) { +PFECALSuperClusterProducer::PFECALSuperClusterProducer(const edm::ParameterSet& iConfig, const reco::SCProducerCache* gcache): + superClusterAlgo_(gcache) { verbose_ = iConfig.getUntrackedParameter("verbose", false); superClusterAlgo_.setUseRegression(iConfig.getParameter("useRegression")); @@ -105,14 +107,12 @@ PFECALSuperClusterProducer::PFECALSuperClusterProducer(const edm::ParameterSet& _theclusteringtype = PFECALSuperClusterAlgo::kBOX; } else if (_typename == ClusterType__Mustache) { _theclusteringtype = PFECALSuperClusterAlgo::kMustache; + } else if (_typename == ClusterType__DeepSC) { + _theclusteringtype = PFECALSuperClusterAlgo::kDeepSC; } else { throw cms::Exception("InvalidClusteringType") << "You have not chosen a valid clustering type," - << " please choose from \"Box\" or \"Mustache\"!"; + << " please choose from \"Box\" or \"Mustache\" or \"DeepSC\"!"; } - superClusterAlgo_.setClusteringType(_theclusteringtype); - superClusterAlgo_.setUseDynamicDPhi(iConfig.getParameter("useDynamicDPhiWindow")); - // clusteringType and useDynamicDPhi need to be defined before setting the tokens in order to esConsume only the necessary records - superClusterAlgo_.setTokens(iConfig, consumesCollector()); std::string _weightname = iConfig.getParameter("EnergyWeight"); if (_weightname == EnergyWeight__Raw) { @@ -130,6 +130,8 @@ PFECALSuperClusterProducer::PFECALSuperClusterProducer(const edm::ParameterSet& // parameters for clustering bool seedThresholdIsET = iConfig.getParameter("seedThresholdIsET"); + bool useDynamicDPhi = iConfig.getParameter("useDynamicDPhiWindow"); + double threshPFClusterSeedBarrel = iConfig.getParameter("thresh_PFClusterSeedBarrel"); double threshPFClusterBarrel = iConfig.getParameter("thresh_PFClusterBarrel"); @@ -142,11 +144,19 @@ PFECALSuperClusterProducer::PFECALSuperClusterProducer(const edm::ParameterSet& double phiwidthSuperClusterEndcap = iConfig.getParameter("phiwidth_SuperClusterEndcap"); double etawidthSuperClusterEndcap = iConfig.getParameter("etawidth_SuperClusterEndcap"); + //double threshPFClusterMustacheOutBarrel = iConfig.getParameter("thresh_PFClusterMustacheOutBarrel"); + //double threshPFClusterMustacheOutEndcap = iConfig.getParameter("thresh_PFClusterMustacheOutEndcap"); + double doSatelliteClusterMerge = iConfig.getParameter("doSatelliteClusterMerge"); double satelliteClusterSeedThreshold = iConfig.getParameter("satelliteClusterSeedThreshold"); double satelliteMajorityFraction = iConfig.getParameter("satelliteMajorityFraction"); bool dropUnseedable = iConfig.getParameter("dropUnseedable"); + superClusterAlgo_.setClusteringType(_theclusteringtype); + superClusterAlgo_.setUseDynamicDPhi(useDynamicDPhi); + // clusteringType and useDynamicDPhi need to be defined before setting the tokens in order to esConsume only the necessary records + superClusterAlgo_.setTokens(iConfig, consumesCollector()); + superClusterAlgo_.setVerbosityLevel(verbose_); superClusterAlgo_.setEnergyWeighting(_theenergyweight); superClusterAlgo_.setUseETForSeeding(seedThresholdIsET); @@ -169,6 +179,8 @@ PFECALSuperClusterProducer::PFECALSuperClusterProducer(const edm::ParameterSet& superClusterAlgo_.setSatelliteThreshold(satelliteClusterSeedThreshold); superClusterAlgo_.setMajorityFraction(satelliteMajorityFraction); superClusterAlgo_.setDropUnseedable(dropUnseedable); + //superClusterAlgo_.setThreshPFClusterMustacheOutBarrel( threshPFClusterMustacheOutBarrel ); + //superClusterAlgo_.setThreshPFClusterMustacheOutEndcap( threshPFClusterMustacheOutEndcap ); //Load the ECAL energy calibration thePFEnergyCalibration_ = std::make_shared(); @@ -331,42 +343,133 @@ void PFECALSuperClusterProducer::produce(edm::Event& iEvent, const edm::EventSet } void PFECALSuperClusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { - edm::ParameterSetDescription desc; - desc.add("PFSuperClusterCollectionEndcap", "particleFlowSuperClusterECALEndcap"); - desc.add("doSatelliteClusterMerge", false); - desc.add("thresh_PFClusterBarrel", 0.0); - desc.add("PFBasicClusterCollectionBarrel", "particleFlowBasicClusterECALBarrel"); - desc.add("useRegression", true); - desc.add("satelliteMajorityFraction", 0.5); - desc.add("thresh_PFClusterEndcap", 0.0); - desc.add("ESAssociation", edm::InputTag("particleFlowClusterECAL")); - desc.add("PFBasicClusterCollectionPreshower", "particleFlowBasicClusterECALPreshower"); - desc.add("use_preshower", true); - desc.addUntracked("verbose", false); - desc.add("thresh_SCEt", 4.0); - desc.add("etawidth_SuperClusterEndcap", 0.04); - desc.add("phiwidth_SuperClusterEndcap", 0.6); - desc.add("useDynamicDPhiWindow", true); - desc.add("PFSuperClusterCollectionBarrel", "particleFlowSuperClusterECALBarrel"); - desc.add("regressionConfig", SCEnergyCorrectorSemiParm::makePSetDescription()); - desc.add("applyCrackCorrections", false); - desc.add("satelliteClusterSeedThreshold", 50.0); - desc.add("etawidth_SuperClusterBarrel", 0.04); - desc.add("PFBasicClusterCollectionEndcap", "particleFlowBasicClusterECALEndcap"); - desc.add("PFClusters", edm::InputTag("particleFlowClusterECAL")); - desc.add("thresh_PFClusterSeedBarrel", 1.0); - desc.add("ClusteringType", "Mustache"); - desc.add("EnergyWeight", "Raw"); - desc.add("BeamSpot", edm::InputTag("offlineBeamSpot")); - desc.add("thresh_PFClusterSeedEndcap", 1.0); - desc.add("phiwidth_SuperClusterBarrel", 0.6); - desc.add("thresh_PFClusterES", 0.0); - desc.add("seedThresholdIsET", true); - desc.add("isOOTCollection", false); - desc.add("barrelRecHits", edm::InputTag("ecalRecHit", "EcalRecHitsEB")); - desc.add("endcapRecHits", edm::InputTag("ecalRecHit", "EcalRecHitsEE")); - desc.add("PFSuperClusterCollectionEndcapWithPreshower", - "particleFlowSuperClusterECALEndcapWithPreshower"); - desc.add("dropUnseedable", false); - descriptions.add("particleFlowSuperClusterECALMustache", desc); + //DeepSC + edm::ParameterSetDescription desc2; + desc2.add("PFSuperClusterCollectionEndcap", "particleFlowDeepSuperClusterECALEndcap"); + desc2.add("doSatelliteClusterMerge", false); + desc2.add("thresh_PFClusterBarrel", 0.0); + desc2.add("PFBasicClusterCollectionBarrel", "particleFlowBasicClusterECALBarrel"); + desc2.add("useRegression", true); + desc2.add("satelliteMajorityFraction", 0.5); + desc2.add("thresh_PFClusterEndcap", 0.0); + desc2.add("ESAssociation", edm::InputTag("particleFlowClusterECAL")); + desc2.add("PFBasicClusterCollectionPreshower", "particleFlowBasicClusterECALPreshower"); + desc2.add("use_preshower", true); + desc2.addUntracked("verbose", false); + desc2.add("thresh_SCEt", 4.0); + desc2.add("etawidth_SuperClusterEndcap", 0.04); + desc2.add("phiwidth_SuperClusterEndcap", 0.6); + desc2.add("useDynamicDPhiWindow", true); + desc2.add("PFSuperClusterCollectionBarrel", "particleFlowDeepSuperClusterECALBarrel"); + { + edm::ParameterSetDescription psd0; + psd0.add("isHLT", false); + psd0.add("isPhaseII", false); + psd0.add("applySigmaIetaIphiBug", false); + psd0.add("ecalRecHitsEE", edm::InputTag("ecalRecHit", "EcalRecHitsEE")); + psd0.add("ecalRecHitsEB", edm::InputTag("ecalRecHit", "EcalRecHitsEB")); + psd0.add("regressionKeyEB", "pfscecal_EBCorrection_offline_v2"); + psd0.add("regressionKeyEE", "pfscecal_EECorrection_offline_v2"); + psd0.add("uncertaintyKeyEB", "pfscecal_EBUncertainty_offline_v2"); + psd0.add("uncertaintyKeyEE", "pfscecal_EEUncertainty_offline_v2"); + psd0.add("regressionMinEB", 0.2); + psd0.add("regressionMaxEB", 2.); + psd0.add("regressionMinEE", 0.2); + psd0.add("regressionMaxEE", 2.); + psd0.add("uncertaintyMinEB", 0.0002); + psd0.add("uncertaintyMaxEB", 0.5); + psd0.add("uncertaintyMinEE", 0.0002); + psd0.add("uncertaintyMaxEE", 0.5); + psd0.add("vertexCollection", edm::InputTag("offlinePrimaryVertices")); + psd0.add("eRecHitThreshold", 1.); + psd0.add("hgcalRecHits", edm::InputTag("")); + psd0.add("hgcalCylinderR", 2.7999999523162842); + desc2.add("regressionConfig", psd0); + } + desc2.add("applyCrackCorrections", false); + desc2.add("satelliteClusterSeedThreshold", 50.0); + desc2.add("etawidth_SuperClusterBarrel", 0.04); + desc2.add("PFBasicClusterCollectionEndcap", "particleFlowBasicClusterECALEndcap"); + desc2.add("PFClusters", edm::InputTag("particleFlowClusterECAL")); + desc2.add("thresh_PFClusterSeedBarrel", 1.0); + desc2.add("ClusteringType", "DeepSC"); + desc2.add("EnergyWeight", "Raw"); + desc2.add("BeamSpot", edm::InputTag("offlineBeamSpot")); + desc2.add("thresh_PFClusterSeedEndcap", 1.0); + desc2.add("phiwidth_SuperClusterBarrel", 0.6); + desc2.add("thresh_PFClusterES", 0.0); + desc2.add("seedThresholdIsET", true); + desc2.add("isOOTCollection", false); + desc2.add("barrelRecHits", edm::InputTag("ecalRecHit", "EcalRecHitsEB")); + desc2.add("endcapRecHits", edm::InputTag("ecalRecHit", "EcalRecHitsEE")); + desc2.add("PFSuperClusterCollectionEndcapWithPreshower", + "particleFlowDeepSuperClusterECALEndcapWithPreshower"); + desc2.add("dropUnseedable", false); + { + edm::ParameterSetDescription psd1; + psd1.add("modelFile","RecoEcal/EgammaClusterProducers/data/DeepSCGraph/model/model.pb"); + psd1.add("scalerFileClusterFeatures", "RecoEcal/EgammaClusterProducers/data/DeepSCGraph/model/scaler_clusters.txt"); + psd1.add("scalerFileWindowFeatures", "RecoEcal/EgammaClusterProducers/data/DeepSCGraph//model/scaler_window.txt"); + psd1.add("nClusterFeatures", 12 ); + psd1.add("nWindowFeatures", 18); + psd1.add("maxNClusters", 45); + psd1.add("maxNRechits", 40); + desc2.add("deepSuperClusterGraphConfig", psd1); + } + descriptions.add("particleFlowSuperClusterECALDeepSC", desc2); + + //Mustache + edm::ParameterSetDescription desc1; + desc1.add("PFSuperClusterCollectionEndcap", "particleFlowSuperClusterECALEndcap"); + desc1.add("doSatelliteClusterMerge", false); + desc1.add("thresh_PFClusterBarrel", 0.0); + desc1.add("PFBasicClusterCollectionBarrel", "particleFlowBasicClusterECALBarrel"); + desc1.add("useRegression", true); + desc1.add("satelliteMajorityFraction", 0.5); + desc1.add("thresh_PFClusterEndcap", 0.0); + desc1.add("ESAssociation", edm::InputTag("particleFlowClusterECAL")); + desc1.add("PFBasicClusterCollectionPreshower", "particleFlowBasicClusterECALPreshower"); + desc1.add("use_preshower", true); + desc1.addUntracked("verbose", false); + desc1.add("thresh_SCEt", 4.0); + desc1.add("etawidth_SuperClusterEndcap", 0.04); + desc1.add("phiwidth_SuperClusterEndcap", 0.6); + desc1.add("useDynamicDPhiWindow", true); + desc1.add("PFSuperClusterCollectionBarrel", "particleFlowSuperClusterECALBarrel"); + { + edm::ParameterSetDescription psd0; + psd0.add("isHLT", false); + psd0.add("applySigmaIetaIphiBug", false); + psd0.add("ecalRecHitsEE", edm::InputTag("ecalRecHit", "EcalRecHitsEE")); + psd0.add("ecalRecHitsEB", edm::InputTag("ecalRecHit", "EcalRecHitsEB")); + psd0.add("regressionKeyEB", "pfscecal_EBCorrection_offline_v2"); + psd0.add("regressionKeyEE", "pfscecal_EECorrection_offline_v2"); + psd0.add("uncertaintyKeyEB", "pfscecal_EBUncertainty_offline_v2"); + psd0.add("uncertaintyKeyEE", "pfscecal_EEUncertainty_offline_v2"); + psd0.add("vertexCollection", edm::InputTag("offlinePrimaryVertices")); + psd0.add("eRecHitThreshold", 1.); + desc1.add("regressionConfig", psd0); + } + desc1.add("applyCrackCorrections", false); + desc1.add("satelliteClusterSeedThreshold", 50.0); + desc1.add("etawidth_SuperClusterBarrel", 0.04); + desc1.add("PFBasicClusterCollectionEndcap", "particleFlowBasicClusterECALEndcap"); + desc1.add("PFClusters", edm::InputTag("particleFlowClusterECAL")); + desc1.add("thresh_PFClusterSeedBarrel", 1.0); + desc1.add("ClusteringType", "Mustache"); + desc1.add("EnergyWeight", "Raw"); + desc1.add("BeamSpot", edm::InputTag("offlineBeamSpot")); + desc1.add("thresh_PFClusterSeedEndcap", 1.0); + desc1.add("phiwidth_SuperClusterBarrel", 0.6); + desc1.add("thresh_PFClusterES", 0.0); + desc1.add("seedThresholdIsET", true); + desc1.add("isOOTCollection", false); + desc1.add("barrelRecHits", edm::InputTag("ecalRecHit", "EcalRecHitsEB")); + desc1.add("endcapRecHits", edm::InputTag("ecalRecHit", "EcalRecHitsEE")); + desc1.add("PFSuperClusterCollectionEndcapWithPreshower", + "particleFlowSuperClusterECALEndcapWithPreshower"); + desc1.add("dropUnseedable", false); + descriptions.add("particleFlowSuperClusterECALMustache", desc1); + } + diff --git a/RecoEcal/EgammaCoreTools/BuildFile.xml b/RecoEcal/EgammaCoreTools/BuildFile.xml index 82cb929f52e52..fc92b4d53b001 100644 --- a/RecoEcal/EgammaCoreTools/BuildFile.xml +++ b/RecoEcal/EgammaCoreTools/BuildFile.xml @@ -1,16 +1,22 @@ + + + + + + diff --git a/RecoEcal/EgammaCoreTools/interface/CalibratedPFCluster.h b/RecoEcal/EgammaCoreTools/interface/CalibratedPFCluster.h new file mode 100644 index 0000000000000..b7f77d544eb4f --- /dev/null +++ b/RecoEcal/EgammaCoreTools/interface/CalibratedPFCluster.h @@ -0,0 +1,24 @@ +#ifndef RecoEcal_EgammaCoreTools_CalibratedPFCluster_h +#define RecoEcal_EgammaCoreTools_CalibratedPFCluster_h + +#include "DataFormats/ParticleFlowReco/interface/PFCluster.h" +#include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h" +#include "DataFormats/CaloRecHit/interface/CaloCluster.h" + +// simple class for associating calibrated energies +class CalibratedPFCluster { +public: + CalibratedPFCluster(const edm::Ptr& p) : cluptr(p) {} + + double energy() const { return cluptr->correctedEnergy(); } + double energy_nocalib() const { return cluptr->energy(); } + double eta() const { return cluptr->positionREP().eta(); } + double phi() const { return cluptr->positionREP().phi(); } + + edm::Ptr the_ptr() const { return cluptr; } + +private: + edm::Ptr cluptr; +}; + +#endif diff --git a/RecoEcal/EgammaCoreTools/interface/DeepSCGraphEvaluation.h b/RecoEcal/EgammaCoreTools/interface/DeepSCGraphEvaluation.h new file mode 100644 index 0000000000000..475a02feefe57 --- /dev/null +++ b/RecoEcal/EgammaCoreTools/interface/DeepSCGraphEvaluation.h @@ -0,0 +1,73 @@ +#ifndef RecoEcal_EgammaCoreTools_DeepSCGraphEvaluation_h +#define RecoEcal_EgammaCoreTools_DeepSCGraphEvaluation_h + +#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include +#include +#include +#include +#include + +//author: Davide Valsecchi +//description: +// Handles Tensorflow DNN graphs and variables scaler configuration. +// To be used for DeepSC. + +namespace reco { + + struct DeepSCConfiguration { + std::string modelFile; + std::string scalerFileClusterFeatures; + std::string scalerFileWindowFeatures; + uint nClusterFeatures; + uint nWindowFeatures; + static constexpr uint nRechitsFeatures = 4; + uint maxNClusters; + uint maxNRechits; + }; + + struct DeepSCInputs { + uint batchSize; + std::vector>> clustersX; + std::vector>>> hitsX; + std::vector> windowX; + std::vector> isSeed; + }; + + class DeepSCGraphEvaluation { + public: + DeepSCGraphEvaluation(const DeepSCConfiguration&); + ~DeepSCGraphEvaluation(); + + std::vector scaleClusterFeatures(const std::vector& input) const; + std::vector scaleWindowFeatures(const std::vector& inputs) const; + + std::vector> evaluate(const DeepSCInputs& inputs) const; + + private: + void initTensorFlowGraphAndSession(); + uint readScalerConfig(std::string file, std::vector>& scalingParams); + + void prepareTensorflowInput(const DeepSCInputs& inputs) const; + + const DeepSCConfiguration cfg_; + std::unique_ptr graphDef_; + tensorflow::Session * session_; + + std::vector> scalerParamsClusters_; + std::vector> scalerParamsWindows_; + + }; + + + + class SCProducerCache { + public: + SCProducerCache(const edm::ParameterSet& conf); + std::unique_ptr deepSCEvaluator; + }; + +}; // namespace reco + +#endif diff --git a/RecoEcal/EgammaCoreTools/interface/EcalClustersGraph.h b/RecoEcal/EgammaCoreTools/interface/EcalClustersGraph.h new file mode 100644 index 0000000000000..33de6b64e9e59 --- /dev/null +++ b/RecoEcal/EgammaCoreTools/interface/EcalClustersGraph.h @@ -0,0 +1,135 @@ +#ifndef RecoEcal_EgammaCoreTools_EcalClustersGraph_h +#define RecoEcal_EgammaCoreTools_EcalClustersGraph_h + +/** + \file + Tools for manipulating ECAL Clusters as graphs + \author Davide Valsecchi, Badder Marzocchi + \date 05 October 2020 +*/ + +#include +#include +#include +#include +#include +#include +#include +#include "TRandom.h" + +#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" +#include "FWCore/Utilities/interface/isFinite.h" + +#include "DataFormats/CaloRecHit/interface/CaloCluster.h" +#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h" +#include "DataFormats/EgammaReco/interface/SuperCluster.h" +#include "DataFormats/ParticleFlowReco/interface/PFCluster.h" + +#include "DataFormats/ParticleFlowReco/interface/PFLayer.h" +#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h" +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "DataFormats/EcalDetId/interface/EEDetId.h" + +#include "Geometry/CaloTopology/interface/CaloTopology.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" +#include "Geometry/Records/interface/CaloTopologyRecord.h" +#include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h" +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h" +#include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h" + +#include "RecoEcal/EgammaCoreTools/interface/CalibratedPFCluster.h" +#include "RecoEcal/EgammaCoreTools/interface/GraphMatrix.h" +#include "RecoEcal/EgammaCoreTools/interface/DeepSCGraphEvaluation.h" + +using namespace std; +using namespace reco; +namespace ublas = boost::numeric::ublas; + +namespace reco { + + class EcalClustersGraph { + typedef std::shared_ptr CalibratedClusterPtr; + typedef std::vector CalibratedClusterPtrVector; + + private: + CalibratedClusterPtrVector clusters_; + uint nSeeds_; + uint nCls_; + + // Adjacency matrix defining which clusters are inside the seeds windows. + // row: seeds (Et ordered), column: clusters (Et ordered) + GraphMatrix inWindows_; + // Adjacency matrix defining how much each cluster is linked to the seed + // row: seeds (Et ordered), column: clusters (Et ordered) + GraphMatrix scoreMatrix_; + GraphMatrix clusterMatrix_; + + //To compute the input variables + const CaloTopology* topology_; + const CaloSubdetectorGeometry* ebGeom_; + const CaloSubdetectorGeometry* eeGeom_; + const EcalRecHitCollection* recHitsEB_; + const EcalRecHitCollection* recHitsEE_; + const SCProducerCache* SCProducerCache_; + + std::array locCov_; + std::pair widths_; + std::vector thresholds_; + DeepSCInputs inputs_; + TRandom* Rnd; + + public: + EcalClustersGraph(CalibratedClusterPtrVector clusters, + int nSeeds, + const CaloTopology* topology, + const CaloSubdetectorGeometry* ebGeom, + const CaloSubdetectorGeometry* eeGeom, + const EcalRecHitCollection* recHitsEB, + const EcalRecHitCollection* recHitsEE, + const SCProducerCache* cache); + + std::vector clusterPosition(const CaloCluster* cluster); + + double deltaPhi(double seed_phi, double cluster_phi) { + double dphi = seed_phi - cluster_phi; + if (dphi > TMath::Pi()) + dphi -= 2 * TMath::Pi(); + if (dphi < -TMath::Pi()) + dphi += 2 * TMath::Pi(); + return dphi; + } + + double deltaEta(double seed_eta, double cluster_eta) { + double deta = 0.; + if (seed_eta > 0.) + deta = cluster_eta - seed_eta; + if (seed_eta <= 0.) + deta = seed_eta - cluster_eta; + return deta; + } + std::vector dynamicWindow(double seedEta); + + std::pair computeCovariances(const CaloCluster* cluster); + std::vector computeShowerShapes(const CaloCluster* cluster, bool full5x5); + std::vector computeVariables(const CaloCluster* seed, const CaloCluster* cluster); + std::vector> fillHits(const CaloCluster* cluster); + std::vector computeWindowVariables(const std::vector>& clusters); + + void fillVariables(); + + double scoreThreshold(const CaloCluster* cluster); + void initWindows(); + void clearWindows(); + + void setThresholds(); + void evaluateScores(); + void selectClusters(); + + void printDebugInfo(); + std::vector> getWindows(); + }; + +} // namespace reco +#endif diff --git a/RecoEcal/EgammaCoreTools/interface/GraphMatrix.h b/RecoEcal/EgammaCoreTools/interface/GraphMatrix.h new file mode 100644 index 0000000000000..9a62a28f78c49 --- /dev/null +++ b/RecoEcal/EgammaCoreTools/interface/GraphMatrix.h @@ -0,0 +1,546 @@ +#ifndef RecoEcal_EgammaCoreTools_GraphMatrix_h +#define RecoEcal_EgammaCoreTools_GraphMatrix_h + +/** + \file + Tools for manipulating ECAL Clusters as graphs + \author Davide Valsecchi, Badder Marzocchi + \date 05 October 2020 +*/ + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; +namespace ublas = boost::numeric::ublas; +typedef size_t size_type; + +template +class GraphMatrix; + +template +std::ostream& operator<<(ostream& s, const GraphMatrix& m); + +template +class GraphMatrix { + friend std::ostream& operator<<<>(ostream& s, const GraphMatrix& m); + typedef typename ublas::matrix::const_iterator1 const_iterator1; + typedef typename ublas::matrix::const_iterator2 const_iterator2; + typedef typename ublas::matrix::iterator1 iterator1; + typedef typename ublas::matrix::iterator2 iterator2; + typedef typename ublas::matrix::const_reverse_iterator1 const_reverse_iterator1; + typedef typename ublas::matrix::const_reverse_iterator2 const_reverse_iterator2; + typedef typename ublas::matrix::reverse_iterator1 reverse_iterator1; + typedef typename ublas::matrix::reverse_iterator2 reverse_iterator2; + +public: + explicit GraphMatrix(); + explicit GraphMatrix(const size_type r, const size_type c); + explicit GraphMatrix(const size_type r, const size_type c, const std::vector> elements, bool isRow); + explicit GraphMatrix(const size_type r, const size_type c, const std::vector>* elements, bool isRow); + explicit GraphMatrix(const ublas::matrix& inMatrix); + ~GraphMatrix(); + + //methods + void Clear(); + void Resize(const size_type r, const size_type c); + size_type nRows() const { return nRows_; }; + size_type nColumns() const { return nColumns_; }; + T Get(const size_type r, const size_type c) const { return matrix_(r, c); }; + ublas::matrix Get() const { return matrix_; }; + std::vector GetRow(const size_type r); + std::vector GetColumn(const size_type c); + void SetZero(const size_type r, const size_type c) { matrix_.erase_element(r, c); }; + void SetRowZero(size_type r); + void SetColumnZero(size_type c); + void Set(const size_type r, const size_type c, const T val) { matrix_.insert_element(r, c, val); }; + void SetRow(size_type r, std::vector row); + void SetRow(size_type r, std::vector* row); + void SetColumn(size_type c, std::vector column); + void SetColumn(size_type c, std::vector* column); + void SetRows(const size_type r, const size_type c, std::vector> rows); + void SetRows(const size_type r, const size_type c, std::vector>* rows); + void SetColumns(const size_type r, const size_type c, std::vector> columns); + void SetColumns(const size_type r, const size_type c, std::vector>* columns); + int nZeros(int i, bool isRow); + int nZeros(std::vector elements); + int nZeros(std::vector* elements); + bool AllZeros(std::vector elements); + bool AllZeros(std::vector* elements); + GraphMatrix Unit(const size_type s) { return GraphMatrix(ublas::identity_matrix(s)); }; + GraphMatrix Zero(const size_type r, const size_type c) { return GraphMatrix(ublas::zero_matrix(r, c)); }; + GraphMatrix Transpose() { return GraphMatrix(trans(matrix_)); }; + GraphMatrix ReduceElements(T maxVal, T defaultGood, std::vector& threshold, bool isRow); + GraphMatrix RemoveDuplicates(T val, bool isRow); + + //iterators + const_iterator1 begin1() const { return matrix_.begin1(); }; + const_iterator2 begin2() const { return matrix_.begin2(); }; + const_iterator1 end1() const { return matrix_.end1(); }; + const_iterator2 end2() const { return matrix_.end2(); }; + iterator1 begin1() { return matrix_.begin1(); }; + iterator2 begin2() { return matrix_.begin2(); }; + iterator1 end1() { return matrix_.end1(); }; + iterator2 end2() { return matrix_.end2(); }; + const_reverse_iterator1 rbegin1() const { return matrix_.rbegin1(); }; + const_reverse_iterator2 rbegin2() const { return matrix_.rbegin2(); }; + const_reverse_iterator1 rend1() const { return matrix_.rend1(); }; + const_reverse_iterator2 rend2() const { return matrix_.rend2(); }; + reverse_iterator1 rbegin1() { return matrix_.rbegin1(); }; + reverse_iterator2 rbegin2() { return matrix_.rbegin2(); }; + reverse_iterator1 rend1() { return matrix_.rend1(); }; + reverse_iterator2 rend2() { return matrix_.rend2(); }; + + //operators + GraphMatrix& operator=(const GraphMatrix& other); + bool operator==(const GraphMatrix& other); + bool operator!=(const GraphMatrix& other); + GraphMatrix operator+(const GraphMatrix& other) { return GraphMatrix(matrix_ + other.Get()); }; + GraphMatrix operator-(const GraphMatrix& other) { return GraphMatrix(matrix_ - other.Get()); }; + GraphMatrix operator*(const GraphMatrix& other) { return GraphMatrix(prod(matrix_, other.Get())); }; + GraphMatrix operator*(const T& other) { return GraphMatrix(matrix_ *= other); }; + GraphMatrix operator/(const T& other) { return GraphMatrix(matrix_ /= other); }; + std::vector operator*(const std::vector& other) { return ToVector(prod(matrix_, ToVector(other))); }; + +private: + size_type nRows_; + size_type nColumns_; + ublas::matrix matrix_; + + //methods + ublas::vector ToVector(const std::vector vector); + std::vector ToVector(const ublas::vector vector); + std::vector ToVector(const ublas::matrix_row> row); + std::vector ToVector(const ublas::matrix_column> column); +}; + +template +GraphMatrix::GraphMatrix() { + nRows_ = 0; + nColumns_ = 0; + matrix_ = ublas::zero_matrix(); +} + +template +GraphMatrix::GraphMatrix(const size_type r, const size_type c) { + nRows_ = r; + nColumns_ = c; + matrix_ = ublas::zero_matrix(nRows_, nColumns_); +} + +template +GraphMatrix::GraphMatrix(const size_type r, + const size_type c, + const std::vector> elements, + bool isRow) { + matrix_.resize(r, c); + nRows_ = r; + nColumns_ = c; + + if (isRow) { + std::vector> rows_ = elements; + rows_.resize(nRows_); + for (typename std::vector>::iterator it = rows_.begin(); it != rows_.end(); ++it) + (*it).resize(nColumns_); + + typename std::vector>::iterator it; + typename std::vector::iterator jt; + for (it = rows_.begin(); it != rows_.end(); ++it) + for (jt = (*it).begin(); jt != (*it).end(); jt++) + matrix_(it - rows_.begin(), jt - (*it).begin()) = *jt; + } else { + std::vector> columns_ = elements; + columns_.resize(nColumns_); + for (typename std::vector>::iterator it = columns_.begin(); it != columns_.end(); ++it) + (*it).resize(nRows_); + + typename std::vector>::iterator it; + typename std::vector::iterator jt; + for (it = columns_.begin(); it != columns_.end(); ++it) + for (jt = (*it).begin(); jt != (*it).end(); jt++) + matrix_(jt - (*it).begin(), it - columns_.begin()) = *jt; + } +} + +template +GraphMatrix::GraphMatrix(const size_type r, + const size_type c, + const std::vector>* elements, + bool isRow) { + matrix_.resize(r, c); + nRows_ = r; + nColumns_ = c; + + if (isRow) { + std::vector> rows_ = *elements; + rows_.resize(nRows_); + for (typename std::vector>::iterator it = rows_.begin(); it != rows_.end(); ++it) + (*it).resize(nColumns_); + + typename std::vector>::iterator it; + typename std::vector::iterator jt; + for (it = rows_.begin(); it != rows_.end(); ++it) + for (jt = (*it).begin(); jt != (*it).end(); jt++) + matrix_(it - rows_.begin(), jt - (*it).begin()) = *jt; + } else { + std::vector> columns_ = *elements; + columns_.resize(nColumns_); + for (typename std::vector>::iterator it = columns_.begin(); it != columns_.end(); ++it) + (*it).resize(nRows_); + + typename std::vector>::iterator it; + typename std::vector::iterator jt; + for (it = columns_.begin(); it != columns_.end(); ++it) + for (jt = (*it).begin(); jt != (*it).end(); jt++) + matrix_(jt - (*it).begin(), it - columns_.begin()) = *jt; + } +} + +template +GraphMatrix::GraphMatrix(const ublas::matrix& inMatrix) { + nRows_ = inMatrix.size1(); + nColumns_ = inMatrix.size2(); + matrix_ = inMatrix; +} + +template +GraphMatrix::~GraphMatrix() { + nRows_ = 0; + nColumns_ = 0; + matrix_.clear(); +} + +//operators +template +std::ostream& operator<<(std::ostream& s, const GraphMatrix& m) { + s << "GraphMatrix: " << std::endl; + for (auto is = m.begin1(); is != m.end1(); is++) { + for (auto ic = is.begin(); ic != is.end(); ic++) { + s << std::setprecision(3) << *ic << " "; + } + s << endl; + } + s << endl; + return s; +} + +template +GraphMatrix& GraphMatrix::operator=(const GraphMatrix& other) { + if (this == &other) + return *this; + this->Resize(other.nRows_, other.nColumns_); + matrix_ = other.Get(); + return *this; +} + +template +bool GraphMatrix::operator==(const GraphMatrix& other) { + if (matrix_.size1() != other.nRows() || matrix_.size2() != other.nColumns()) + return false; + else { + for (auto is = matrix_.begin1(); is != matrix_.end1(); ++is) + for (auto ic = is.begin(); ic != is.end(); ic++) + if (matrix_(is - matrix_.begin1(), ic - is.begin()) != other.Get(is - matrix_.begin1(), ic - is.begin())) + return false; + } + return true; +} + +template +bool GraphMatrix::operator!=(const GraphMatrix& other) { + if (matrix_.size1() != other.nRows() || matrix_.size2() != other.nColumns()) + return true; + else { + for (auto is = matrix_.begin1(); is != matrix_.end1(); ++is) + for (auto ic = is.begin(); ic != is.end(); ic++) + if (matrix_(is - matrix_.begin1(), ic - is.begin()) != other.Get(is - matrix_.begin1(), ic - is.begin())) + return true; + } + return false; +} + +//other methods +template +void GraphMatrix::Clear() { + nRows_ = 0; + nColumns_ = 0; + matrix_.clear(); +} + +template +void GraphMatrix::Resize(const size_type r, const size_type c) { + nRows_ = r; + nColumns_ = c; + matrix_.resize(nRows_, nColumns_); +} + +template +std::vector GraphMatrix::GetRow(const size_type r) { + size_type row_ = r; + if (row_ >= nRows_) + row_ = nRows_ - 1; + return ToVector(row(matrix_, row_)); +} + +template +std::vector GraphMatrix::GetColumn(const size_type c) { + size_type column_ = c; + if (column_ >= nColumns_) + column_ = nColumns_ - 1; + return ToVector(column(matrix_, column_)); +}; + +template +void GraphMatrix::SetRowZero(size_type r) { + std::vector row_; + row_.resize(nColumns_); + if (r >= nRows_) { + Resize(nRows_ + 1, nColumns_); + r = nRows_ - 1; + } + + for (typename std::vector::iterator it = row_.begin(); it != row_.end(); ++it) + matrix_.erase_element(r, it - row_.begin()); + + nRows_ = matrix_.size1(); + nColumns_ = matrix_.size2(); +} + +template +void GraphMatrix::SetRow(size_type r, std::vector row) { + std::vector row_ = row; + row_.resize(nColumns_); + if (r >= nRows_) { + Resize(nRows_ + 1, nColumns_); + r = nRows_ - 1; + } + + for (typename std::vector::iterator it = row_.begin(); it != row_.end(); ++it) + matrix_.insert_element(r, it - row_.begin(), *it); + + nRows_ = matrix_.size1(); + nColumns_ = matrix_.size2(); +} + +template +void GraphMatrix::SetRow(size_type r, std::vector* row) { + std::vector row_ = *row; + row_.resize(nColumns_); + if (r >= nRows_) { + Resize(nRows_ + 1, nColumns_); + r = nRows_ - 1; + } + + for (typename std::vector::iterator it = row_.begin(); it != row_.end(); ++it) + matrix_.insert_element(r, it - row_.begin(), *it); + + nRows_ = matrix_.size1(); + nColumns_ = matrix_.size2(); +} + +template +void GraphMatrix::SetColumnZero(size_type c) { + std::vector column_; + column_.resize(nRows_); + if (c >= nColumns_) { + Resize(nRows_, nColumns_ + 1); + c = nColumns_ - 1; + } + + for (typename std::vector::iterator it = column_.begin(); it != column_.end(); ++it) + matrix_.erase_element(it - column_.begin(), c); + + nRows_ = matrix_.size1(); + nColumns_ = matrix_.size2(); +} + +template +void GraphMatrix::SetColumn(size_type c, std::vector column) { + std::vector column_ = column; + column_.resize(nRows_); + if (c >= nColumns_) { + Resize(nRows_, nColumns_ + 1); + c = nColumns_ - 1; + } + + for (typename std::vector::iterator it = column_.begin(); it != column_.end(); ++it) + matrix_.insert_element(it - column_.begin(), c, *it); + + nRows_ = matrix_.size1(); + nColumns_ = matrix_.size2(); +} + +template +void GraphMatrix::SetColumn(size_type c, std::vector* column) { + std::vector column_ = *column; + column_.resize(nRows_); + if (c >= nColumns_) { + Resize(nRows_, nColumns_ + 1); + c = nColumns_ - 1; + } + + for (typename std::vector::iterator it = column_.begin(); it != column_.end(); ++it) + matrix_.insert_element(it - column_.begin(), c, *it); + + nRows_ = matrix_.size1(); + nColumns_ = matrix_.size2(); +} + +template +int GraphMatrix::nZeros(int i, bool isRow) { + std::vector vector_; + if (isRow) + vector_ = this->GetRow(i); + else + vector_ = this->GetColumn(i); + + int n = std::count(vector_.begin(), vector_.end(), T(0)); + return n; +} + +template +int GraphMatrix::nZeros(std::vector elements) { + int n = std::count(elements.begin(), elements.end(), T(0)); + return n; +} + +template +int GraphMatrix::nZeros(std::vector* elements) { + int n = std::count(elements->begin(), elements->end(), T(0)); + return n; +} + +template +bool GraphMatrix::AllZeros(std::vector elements) { + bool isZeros = false; + if (nZeros(&elements) == (int)elements.size()) + isZeros = true; + return isZeros; +} + +template +bool GraphMatrix::AllZeros(std::vector* elements) { + bool isZeros = false; + if (nZeros(elements) == (int)elements->size()) + isZeros = true; + return isZeros; +} + +template +GraphMatrix GraphMatrix::ReduceElements(T maxVal, T defaultGood, std::vector& thresholds, bool isRow) { + GraphMatrix m_(nRows_, nColumns_); + + size_type n; + if (isRow) + n = nRows_; + else + n = nColumns_; + thresholds.resize(n); + + for (size_type i = 0; i < n; i++) { + std::vector vector_; + if (isRow) + vector_ = this->GetRow(i); + else + vector_ = this->GetColumn(i); + + T threshold = thresholds[i]; + if (AllZeros(&vector_)) { + if (isRow) + m_.SetRow(i, vector_); + else + m_.SetColumn(i, vector_); + } else { + std::replace_if( + vector_.begin(), vector_.end(), [threshold](T& v) { return v <= threshold; }, T(0)); + std::replace_if( + vector_.begin(), vector_.end(), [threshold](T& v) { return v > threshold; }, defaultGood); + + if (isRow) + m_.SetRow(i, vector_); + else + m_.SetColumn(i, vector_); + } + } + + return m_; +} + +template +GraphMatrix GraphMatrix::RemoveDuplicates(T val, bool isRow) { + GraphMatrix m_(nRows_, nColumns_); + + size_type n; + if (isRow) + n = nRows_; + else + n = nColumns_; + + for (size_type i = 0; i < n; i++) { + std::vector vector_; + if (isRow) + vector_ = this->GetRow(i); + else + vector_ = this->GetColumn(i); + + if (AllZeros(&vector_)) { + if (isRow) + m_.SetRow(i, vector_); + else + m_.SetColumn(i, vector_); + } else { + typename std::vector::iterator it = std::find(vector_.begin(), vector_.end(), val); + + T val = *it; + int index = it - vector_.begin(); + + vector_ = std::vector(vector_.size(), T(0)); + if (it != vector_.end()) + vector_[index] = val; + + if (isRow) + m_.SetRow(i, vector_); + else + m_.SetColumn(i, vector_); + } + } + + return m_; +} + +//Private methods +template +ublas::vector GraphMatrix::ToVector(const std::vector vector) { + ublas::vector v(vector.size()); + std::copy(vector.begin(), vector.end(), v.begin()); + return v; +} + +template +std::vector GraphMatrix::ToVector(const ublas::vector vector) { + std::vector v(vector.size()); + std::copy(vector.begin(), vector.end(), v.begin()); + return v; +} + +template +std::vector GraphMatrix::ToVector(const ublas::matrix_row> row) { + std::vector v(row.size()); + std::copy(row.begin(), row.end(), v.begin()); + return v; +} + +template +std::vector GraphMatrix::ToVector(const ublas::matrix_column> column) { + std::vector v(column.size()); + std::copy(column.begin(), column.end(), v.begin()); + return v; +} + +#endif diff --git a/RecoEcal/EgammaCoreTools/src/DeepSCGraphEvaluation.cc b/RecoEcal/EgammaCoreTools/src/DeepSCGraphEvaluation.cc new file mode 100644 index 0000000000000..2d047f0cdf00b --- /dev/null +++ b/RecoEcal/EgammaCoreTools/src/DeepSCGraphEvaluation.cc @@ -0,0 +1,213 @@ +#include "RecoEcal/EgammaCoreTools/interface/DeepSCGraphEvaluation.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/FileInPath.h" +#include "TMath.h" +#include +#include +using namespace reco; + +DeepSCGraphEvaluation::DeepSCGraphEvaluation(const DeepSCConfiguration& cfg) : cfg_(cfg) { + tensorflow::setLogging("0"); + // Init TF graph and session objects + initTensorFlowGraphAndSession(); + // Init scaler configs + uint nClFeat = readScalerConfig(cfg_.scalerFileClusterFeatures, scalerParamsClusters_); + if (nClFeat != cfg_.nClusterFeatures) { + throw cms::Exception("WrongConfiguration") << "Mismatch between number of input features for Clusters and " + << "parameters in the scaler file."; + } + uint nClWind = readScalerConfig(cfg_.scalerFileWindowFeatures, scalerParamsWindows_); + if (nClWind != cfg_.nWindowFeatures) { + throw cms::Exception("WrongConfiguration") << "Mismatch between number of input features for Clusters and " + << "parameters in the scaler file."; + } +} + +DeepSCGraphEvaluation::~DeepSCGraphEvaluation(){ + if(session_ != nullptr) tensorflow::closeSession(session_); +} + +void DeepSCGraphEvaluation::initTensorFlowGraphAndSession() { + // load the graph definition + LogDebug("DeepSCGraphEvaluation") << "Loading graph"; + graphDef_ = + std::unique_ptr(tensorflow::loadGraphDef(edm::FileInPath(cfg_.modelFile).fullPath())); + LogDebug("DeepSCGraphEvaluation") << "Starting TF sessions"; + session_ = tensorflow::createSession(graphDef_.get()); + LogDebug("DeepSCGraphEvaluation") << "TF ready"; +} + +uint DeepSCGraphEvaluation::readScalerConfig(std::string file, std::vector>& scalingParams) { + LogDebug("DeepSCGraphEvaluation") << "Reading scaler file: "<< edm::FileInPath(file).fullPath(); + std::ifstream inputfile{edm::FileInPath(file).fullPath()}; + int ninputs = 0; + if (inputfile.fail()) { + throw cms::Exception("MissingFile") << "Scaler file not found: " << file; + } else { + // Now read mean, scale factors for each variable + float par1, par2; + while (inputfile >> par1 >> par2) { + scalingParams.push_back(std::make_pair(par1, par2)); + ninputs += 1; + } + } + return ninputs; +} + +std::vector DeepSCGraphEvaluation::scaleClusterFeatures( + const std::vector& input) const { + std::vector out(input.size()); + for (size_t i = 0; i < input.size(); i++) { + const auto& [par1, par2] = scalerParamsClusters_[i]; + out[i] = (input[i] - par1) / par2; + } + return out; +} + +std::vector DeepSCGraphEvaluation::scaleWindowFeatures(const std::vector& input) const { + std::vector out(input.size()); + for (size_t i = 0; i < input.size(); i++) { + const auto& [par1, par2] = scalerParamsWindows_[i]; + out[i] = (input[i] - par1) / par2; + } + return out; +} + +std::vector> DeepSCGraphEvaluation::evaluate(const DeepSCInputs& inputs) const { + /* + Evaluate the DeepSC model + */ + LogDebug("DeepSCGraphEvaluation") << "Starting evaluation"; + // Input tensors initialization + tensorflow::Tensor clsX {tensorflow::DT_FLOAT, {inputs.batchSize, cfg_.maxNClusters, cfg_.nClusterFeatures}}; + tensorflow::Tensor windX {tensorflow::DT_FLOAT, {inputs.batchSize, cfg_.nWindowFeatures}}; + tensorflow::Tensor hitsX {tensorflow::DT_FLOAT, {inputs.batchSize, cfg_.maxNClusters, cfg_.maxNRechits, cfg_.nRechitsFeatures }}; + tensorflow::Tensor isSeedX {tensorflow::DT_FLOAT, {inputs.batchSize, cfg_.maxNClusters, 1}}; + tensorflow::Tensor nClsSize {tensorflow::DT_FLOAT, {inputs.batchSize}}; + + float * C = clsX.flat().data(); + // Look on batch dim + for (const auto & cls_data : inputs.clustersX ){ + // Loop on clusters + for (size_t k = 0; k < cfg_.maxNClusters; k++){ + // Loop on features + for (size_t z=0; z < cfg_.nClusterFeatures; z++, C++){//--> note the double loop on the tensor pointer + if (k < cls_data.size()){ + *C = float(cls_data[k][z]); + }else{ + *C = 0.; + } + } + } + } + + float * W = windX.flat().data(); + // Look on batch dim + for (const auto & wind_features : inputs.windowX ){ + // Loop on features + for (size_t k = 0; k < cfg_.nWindowFeatures; k++, W++){ //--> note the double loop on the tensor pointer + *W = float(wind_features[k]); + } + } + + float * H = hitsX.flat().data(); + size_t iW = -1; + // Look on batch dim + for (const auto & hits_data : inputs.hitsX ){ + iW++; + size_t ncls_in_window = hits_data.size(); + // Loop on clusters + for (size_t k = 0; k < cfg_.maxNClusters; k++){ //--> note the triple loop on the tensor pointer + // Check padding + size_t nhits_in_cluster; + if (k < ncls_in_window) nhits_in_cluster = hits_data[k].size(); + else nhits_in_cluster = 0; + + // Loop on hits + for (size_t j=0; j < cfg_.maxNRechits; j++){//--> note the triple loop on the tensor pointer + // Check the number of clusters and hits for padding + bool ok = j < nhits_in_cluster; + // Loop on rechits features + for (size_t z=0; z< cfg_.nRechitsFeatures; z++, H++){//--> note the triple loop on the tensor pointe + if (ok) *H = float( hits_data[k][j][z]); + else *H = 0.; + } + } + } + } + + float * S = isSeedX.flat().data(); + // Look on batch dim + for (const auto & isSeed_data : inputs.isSeed ){ + // Loop on clusters + for (size_t k = 0; k < cfg_.maxNClusters; k++, S++){ //--> note the double loop on the tensor pointer + if (k < isSeed_data.size()){ + *S = float(isSeed_data[k]); + }else{ + *S = 0.; + } + } + } + + float * M = nClsSize.flat().data(); + for (size_t k = 0; k < inputs.batchSize; k++, M++){ + *M = float(inputs.clustersX[k].size()); + } + + std::vector> feed_dict = { + { "input_1", clsX }, + { "input_2", windX}, + { "input_3", hitsX}, + { "input_4", isSeedX}, + { "input_5", nClsSize} + }; + + // prepare tensorflow outputs + std::vector outputs_tf; + // // Define the output and run + std::vector> outputs_clustering; + // // Run the models + LogDebug("DeepSCGraphEvaluation") << "Run model"; + tensorflow::run(session_, feed_dict, {"Identity", "Identity_1","Identity_2","Identity_3"}, &outputs_tf); + + // Reading the 1st output: clustering probability + // const auto& r = outputs_tf[0].tensor(); + + float * y_cl = outputs_tf[0].flat().data(); + // Iterate on the clusters for each window + for (size_t b = 0; b< inputs.batchSize; b++) { + uint ncls = inputs.clustersX[b].size(); + std::vector cl_output(ncls); + for (size_t c = 0; c < ncls; c++){ + float y = y_cl[b*cfg_.maxNClusters + c]; + cl_output[c] = 1 / (1 + TMath::Exp(- y)); + } + std::cout << b << ") "; + std::for_each(cl_output.begin(), cl_output.end(), [](float x){std::cout <("ClusteringType"); + const auto& pset_dnn = conf.getParameter("deepSuperClusterGraphConfig"); + + if (clustering_type == "DeepSC") { + config.modelFile = pset_dnn.getParameter("modelFile"); + config.scalerFileClusterFeatures = pset_dnn.getParameter("scalerFileClusterFeatures"); + config.scalerFileWindowFeatures = pset_dnn.getParameter("scalerFileWindowFeatures"); + config.nClusterFeatures = pset_dnn.getParameter("nClusterFeatures"); + config.nWindowFeatures = pset_dnn.getParameter("nWindowFeatures"); + config.maxNClusters = pset_dnn.getParameter("maxNClusters"); + config.maxNRechits = pset_dnn.getParameter("maxNRechits"); + deepSCEvaluator = std::make_unique(config); + } +} diff --git a/RecoEcal/EgammaCoreTools/src/EcalClustersGraph.cc b/RecoEcal/EgammaCoreTools/src/EcalClustersGraph.cc new file mode 100644 index 0000000000000..f5096466c19c4 --- /dev/null +++ b/RecoEcal/EgammaCoreTools/src/EcalClustersGraph.cc @@ -0,0 +1,526 @@ +#include "RecoEcal/EgammaCoreTools/interface/EcalClustersGraph.h" +#include +#include +#include "TVector2.h" +#include "TMath.h" +#include + +using namespace std; +using namespace reco; + +typedef std::shared_ptr CalibratedClusterPtr; +typedef std::vector CalibratedClusterPtrVector; + +EcalClustersGraph::EcalClustersGraph(CalibratedClusterPtrVector clusters, + int nSeeds, + const CaloTopology* topology, + const CaloSubdetectorGeometry* ebGeom, + const CaloSubdetectorGeometry* eeGeom, + const EcalRecHitCollection* recHitsEB, + const EcalRecHitCollection* recHitsEE, + const SCProducerCache* cache) + : clusters_(clusters), + nSeeds_(nSeeds), + topology_(topology), + ebGeom_(ebGeom), + eeGeom_(eeGeom), + recHitsEB_(recHitsEB), + recHitsEE_(recHitsEE), + SCProducerCache_(cache) { + nCls_ = clusters_.size(); + inWindows_ = GraphMatrix(nSeeds_, nCls_); + scoreMatrix_ = GraphMatrix(nSeeds_, nCls_); + clusterMatrix_ = GraphMatrix(nSeeds_, nCls_); + Rnd = new TRandom(); + + // Prepare the batch size of the tensor inputs == number of windows + inputs_.clustersX.resize(nSeeds_); + inputs_.windowX.resize(nSeeds_); + inputs_.hitsX.resize(nSeeds_); + inputs_.isSeed.resize(nSeeds_); + + LogDebug("EcalClustersGraph") << "EcalClustersGraph created. nSeeds " << nSeeds_ << ", nClusters " << nCls_ << endl; +} + +std::vector EcalClustersGraph::clusterPosition(const CaloCluster* cluster) { + std::vector coordinates; + coordinates.resize(3); + int ieta = -999; + int iphi = -999; + int iz = -99; + + math::XYZPoint caloPos = cluster->position(); + if (PFLayer::fromCaloID(cluster->caloID()) == PFLayer::ECAL_BARREL) { + EBDetId eb_id(ebGeom_->getClosestCell(GlobalPoint(caloPos.x(), caloPos.y(), caloPos.z()))); + ieta = eb_id.ieta(); + iphi = eb_id.iphi(); + iz = 0; + } else if (PFLayer::fromCaloID(cluster->caloID()) == PFLayer::ECAL_ENDCAP) { + EEDetId ee_id(eeGeom_->getClosestCell(GlobalPoint(caloPos.x(), caloPos.y(), caloPos.z()))); + ieta = ee_id.ix(); + iphi = ee_id.iy(); + if (ee_id.zside() < 0) + iz = -1; + if (ee_id.zside() > 0) + iz = 1; + } + + coordinates[0] = ieta; + coordinates[1] = iphi; + coordinates[2] = iz; + return coordinates; +} + +std::vector EcalClustersGraph::dynamicWindow(double seedEta) { + std::vector window; + window.resize(3); + + double eta = fabs(seedEta); + double deta_down = 0.; + double deta_up = 0.; + double dphi = 0.; + + //deta_down + if (eta < 2.1) + deta_down = -0.075; + else if (eta >= 2.1 && eta < 2.5) + deta_down = -0.1875 * eta + 0.31875; + else if (eta >= 2.5) + deta_down = -0.15; + + //deta_up + if (eta >= 0 && eta < 0.1) + deta_up = 0.075; + else if (eta >= 0.1 && eta < 1.3) + deta_up = 0.0758929 - 0.0178571 * eta + 0.0892857 * (eta * eta); + else if (eta >= 1.3 && eta < 1.7) + deta_up = 0.2; + else if (eta >= 1.7 && eta < 1.9) + deta_up = 0.625 - 0.25 * eta; + else if (eta >= 1.9) + deta_up = 0.15; + + //dphi + if (eta < 1.9) + dphi = 0.6; + else if (eta >= 1.9 && eta < 2.7) + dphi = 1.075 - 0.25 * eta; + else if (eta >= 2.7) + dphi = 0.4; + + window[0] = deta_down; + window[1] = deta_up; + window[2] = dphi; + + return window; +} + +void EcalClustersGraph::initWindows() { + for (uint is = 0; is < nSeeds_; is++) { + std::vector seedLocal = clusterPosition((*clusters_.at(is)).the_ptr().get()); + double seed_eta = clusters_.at(is)->eta(); + double seed_phi = clusters_.at(is)->phi(); + inWindows_.Set(is, is, 1); + std::vector width = dynamicWindow(seed_eta); + + for (uint icl = is + 1; icl < nCls_; icl++) { + std::vector clusterLocal = clusterPosition((*clusters_.at(icl)).the_ptr().get()); + double cl_eta = clusters_.at(icl)->eta(); + double cl_phi = clusters_.at(icl)->phi(); + double dphi = deltaPhi(seed_phi, cl_phi); + double deta = deltaEta(seed_eta, cl_eta); + + int isIn = 0; + if (seedLocal[2] == clusterLocal[2] && deta >= width[0] && deta <= width[1] && fabs(dphi) <= width[2]) + isIn = 1; + + inWindows_.Set(is, icl, isIn); + //Save also symmetric part of the adj matrix + if (icl < nSeeds_) + inWindows_.Set(icl, is, isIn); + } + } +} + +void EcalClustersGraph::clearWindows() { + inWindows_.Clear(); + scoreMatrix_.Clear(); + clusterMatrix_.Clear(); +} + +std::pair EcalClustersGraph::computeCovariances(const CaloCluster* cluster) +{ + + double etaWidth = 0.; + double phiWidth = 0.; + double numeratorEtaWidth = 0; + double numeratorPhiWidth = 0; + + double clEnergy = cluster->energy(); + double denominator = clEnergy; + + double clEta = cluster->position().eta(); + double clPhi = cluster->position().phi(); + + std::shared_ptr this_cell; + EcalRecHitCollection::const_iterator rHit; + + const std::vector >& detId = cluster->hitsAndFractions(); + // Loop over recHits associated with the given SuperCluster + for (std::vector >::const_iterator hit = detId.begin(); hit != detId.end(); ++hit) { + if(PFLayer::fromCaloID(cluster->caloID()) == PFLayer::ECAL_BARREL){ + rHit = recHitsEB_->find((*hit).first); + //FIXME: THIS IS JUST A WORKAROUND A FIX SHOULD BE APPLIED + if (rHit == recHitsEB_->end()) { + continue; + } + }else if(PFLayer::fromCaloID(cluster->caloID()) == PFLayer::ECAL_ENDCAP){ + rHit = recHitsEE_->find((*hit).first); + //FIXME: THIS IS JUST A WORKAROUND A FIX SHOULD BE APPLIED + if (rHit == recHitsEE_->end()) { + continue; + } + } + + if(PFLayer::fromCaloID(cluster->caloID()) == PFLayer::ECAL_BARREL){ + this_cell = ebGeom_->getGeometry(rHit->id()); + }else if(PFLayer::fromCaloID(cluster->caloID()) == PFLayer::ECAL_ENDCAP){ + this_cell = eeGeom_->getGeometry(rHit->id()); + } + if (this_cell == nullptr) { + //edm::LogInfo("SuperClusterShapeAlgo") << "pointer to the cell in Calculate_Covariances is NULL!"; + continue; + } + + GlobalPoint position = this_cell->getPosition(); + //take into account energy fractions + double energyHit = rHit->energy() * hit->second; + + //form differences + double dPhi = position.phi() - clPhi; + if (dPhi > +Geom::pi()) { + dPhi = Geom::twoPi() - dPhi; + } + if (dPhi < -Geom::pi()) { + dPhi = Geom::twoPi() + dPhi; + } + + double dEta = position.eta() - clEta; + + if (energyHit > 0) { + numeratorEtaWidth += energyHit * dEta * dEta; + numeratorPhiWidth += energyHit * dPhi * dPhi; + } + + etaWidth = sqrt(numeratorEtaWidth / denominator); + phiWidth = sqrt(numeratorPhiWidth / denominator); + } + + return std::make_pair(etaWidth,phiWidth); +} + +std::vector EcalClustersGraph::computeShowerShapes(const CaloCluster* cluster, bool full5x5=false) +{ + std::vector showerVars_; + showerVars_.resize(8); + float e1=1.; + float e4=0.; + + if(full5x5){ + if(PFLayer::fromCaloID(cluster->caloID()) == PFLayer::ECAL_BARREL){ + locCov_ = noZS::EcalClusterTools::localCovariances(*cluster, recHitsEB_, topology_); + widths_ = computeCovariances(cluster); + e1 = noZS::EcalClusterTools::eMax(*cluster, recHitsEB_); + e4 = noZS::EcalClusterTools::eTop(*cluster, recHitsEB_, topology_) + + noZS::EcalClusterTools::eRight(*cluster, recHitsEB_, topology_) + + noZS::EcalClusterTools::eBottom(*cluster, recHitsEB_, topology_) + + noZS::EcalClusterTools::eLeft(*cluster, recHitsEB_, topology_); + showerVars_[0] = noZS::EcalClusterTools::e3x3(*cluster, recHitsEB_, topology_)/cluster->energy(); //r9 + showerVars_[1] = sqrt(locCov_[0]); //sigmaietaieta + showerVars_[2] = locCov_[1]; //sigmaietaiphi + showerVars_[3] = (!edm::isFinite(locCov_[2])) ? 0. : sqrt(locCov_[2]); //sigmaiphiiphi + showerVars_[4] = (e1!=0.) ? 1.-e4/e1 : -999.; //swiss_cross + showerVars_[5] = cluster->hitsAndFractions().size(); //nXtals + showerVars_[6] = widths_.first; //etaWidth + showerVars_[7] = widths_.second; //phiWidth + }else if(PFLayer::fromCaloID(cluster->caloID()) == PFLayer::ECAL_ENDCAP){ + locCov_ = noZS::EcalClusterTools::localCovariances(*cluster, recHitsEE_, topology_); + widths_ = computeCovariances(cluster); + e1 = noZS::EcalClusterTools::eMax(*cluster, recHitsEE_); + e4 = noZS::EcalClusterTools::eTop(*cluster, recHitsEE_, topology_) + + noZS::EcalClusterTools::eRight(*cluster, recHitsEE_, topology_) + + noZS::EcalClusterTools::eBottom(*cluster, recHitsEE_, topology_) + + noZS::EcalClusterTools::eLeft(*cluster, recHitsEE_, topology_); + showerVars_[0] = noZS::EcalClusterTools::e3x3(*cluster, recHitsEE_, topology_)/cluster->energy(); //r9 + showerVars_[1] = sqrt(locCov_[0]); //sigmaietaieta + showerVars_[2] = locCov_[1]; //sigmaietaiphi + showerVars_[3] = (!edm::isFinite(locCov_[2])) ? 0. : sqrt(locCov_[2]); //sigmaiphiiphi + showerVars_[4] = (e1!=0.) ? 1.-e4/e1 : -999.; //swiss_cross + showerVars_[5] = cluster->hitsAndFractions().size(); //nXtals + showerVars_[6] = widths_.first; //etaWidth + showerVars_[7] = widths_.second; //phiWidth + } + }else{ + if(PFLayer::fromCaloID(cluster->caloID()) == PFLayer::ECAL_BARREL){ + locCov_ = EcalClusterTools::localCovariances(*cluster, recHitsEB_, topology_); + widths_ = computeCovariances(cluster); + e1 = EcalClusterTools::eMax(*cluster, recHitsEB_); + e4 = EcalClusterTools::eTop(*cluster, recHitsEB_, topology_) + + EcalClusterTools::eRight(*cluster, recHitsEB_, topology_) + + EcalClusterTools::eBottom(*cluster, recHitsEB_, topology_) + + EcalClusterTools::eLeft(*cluster, recHitsEB_, topology_); + showerVars_[0] = EcalClusterTools::e3x3(*cluster, recHitsEB_, topology_)/cluster->energy(); //r9 + showerVars_[1] = sqrt(locCov_[0]); //sigmaietaieta + showerVars_[2] = locCov_[1]; //sigmaietaiphi + showerVars_[3] = (!edm::isFinite(locCov_[2])) ? 0. : sqrt(locCov_[2]); //sigmaiphiiphi + showerVars_[4] = (e1!=0.) ? 1.-e4/e1 : -999.; //swiss_cross + showerVars_[5] = cluster->hitsAndFractions().size(); //nXtals + showerVars_[6] = widths_.first; //etaWidth + showerVars_[7] = widths_.second; //phiWidth + }else if(PFLayer::fromCaloID(cluster->caloID()) == PFLayer::ECAL_ENDCAP){ + locCov_ = EcalClusterTools::localCovariances(*cluster, recHitsEE_, topology_); + widths_ = computeCovariances(cluster); + e1 = EcalClusterTools::eMax(*cluster, recHitsEE_); + e4 = EcalClusterTools::eTop(*cluster, recHitsEE_, topology_) + + EcalClusterTools::eRight(*cluster, recHitsEE_, topology_) + + EcalClusterTools::eBottom(*cluster, recHitsEE_, topology_) + + EcalClusterTools::eLeft(*cluster, recHitsEE_, topology_); + showerVars_[0] = EcalClusterTools::e3x3(*cluster, recHitsEE_, topology_)/cluster->energy(); //r9 + showerVars_[1] = sqrt(locCov_[0]); //sigmaietaieta + showerVars_[2] = locCov_[1]; //sigmaietaiphi + showerVars_[3] = (!edm::isFinite(locCov_[2])) ? 0. : sqrt(locCov_[2]); //sigmaiphiiphi + showerVars_[4] = (e1!=0.) ? 1.-e4/e1 : -999.; //swiss_cross + showerVars_[5] = cluster->hitsAndFractions().size(); //nXtals + showerVars_[6] = widths_.first; //etaWidth + showerVars_[7] = widths_.second; //phiWidth + } + } + + return showerVars_; +} + +std::vector> EcalClustersGraph::fillHits(const CaloCluster* cluster) { + const std::vector>& hitsAndFractions = cluster->hitsAndFractions(); + std::vector> out (hitsAndFractions.size()); + if (hitsAndFractions.size()==0){ + edm::LogError("EcalClustersGraph") << "No hits in cluster!!"; + } + for (unsigned int i = 0; i < hitsAndFractions.size(); i++) { + std::vector rechit (DeepSCConfiguration::nRechitsFeatures); + if (hitsAndFractions[i].first.subdetId() == EcalBarrel) { + double energy = (*recHitsEB_->find(hitsAndFractions[i].first)).energy(); + EBDetId eb_id(hitsAndFractions[i].first); + rechit[0] = eb_id.ieta(); //ieta + rechit[1] = eb_id.iphi(); //iphi + rechit[2] = 0.; //iz + // rechit[3] = energy; //energy + rechit[3] = energy * hitsAndFractions[i].second; //energy * fraction + // rechit[5] = hitsAndFractions[i].second; //fraction + } else if (hitsAndFractions[i].first.subdetId() == EcalEndcap) { + double energy = (*recHitsEE_->find(hitsAndFractions[i].first)).energy(); + EEDetId ee_id(hitsAndFractions[i].first); + rechit[0] = ee_id.ix(); //ix + rechit[1] = ee_id.iy(); //iy + if (ee_id.zside() < 0) + rechit[2] = -1.; //iz + if (ee_id.zside() > 0) + rechit[2] = +1.; //iz + rechit[3] = energy * hitsAndFractions[i].second; //energy * fraction + // rechit[3] = energy; //energy + // rechit[5] = hitsAndFractions[i].second; //fraction + }else{ + edm::LogError("EcalClustersGraph") << "Rechit is not either EB or EE!!"; + } + out[i] = rechit; + } + return out; +} + +std::vector EcalClustersGraph::computeVariables(const CaloCluster* seed, const CaloCluster* cluster) { + std::vector cl_vars(12); //TODO == PUT dynamic configuration + //showerShapes_ = computeShowerShapes(cluster,false); + std::vector clusterLocal = clusterPosition(cluster); + + cl_vars[0] = cluster->energy(); //cl_energy + cl_vars[1] = cluster->energy() / TMath::CosH(cluster->eta()); //cl_et + cl_vars[2] = cluster->eta(); //cl_eta + cl_vars[3] = cluster->phi(); //cl_phi + cl_vars[4] = clusterLocal[0]; //cl_ieta/ix + cl_vars[5] = clusterLocal[1]; //cl_iphi/iy + cl_vars[6] = clusterLocal[2]; //cl_iz + cl_vars[7] = deltaEta(seed->eta(), cluster->eta()); //cl_dEta + cl_vars[8] = deltaPhi(seed->phi(), cluster->phi()); //cl_dPhi + cl_vars[9] = seed->energy() - cluster->energy(); //cl_dEnergy + cl_vars[10] = + (seed->energy() / TMath::CosH(seed->eta())) - (cluster->energy() / TMath::CosH(cluster->eta())); //cl_dEt + cl_vars[11] = cluster->hitsAndFractions().size(); // nxtals + // cl_vars[12] = showerShapes_[0]; //cl_r9 + // cl_vars[13] = showerShapes_[1]; //cl_sigmaietaieta + // cl_vars[14] = showerShapes_[2]; //cl_sigmaietaiphi + // cl_vars[15] = showerShapes_[3]; //cl_sigmaiphiiphi + // cl_vars[16] = showerShapes_[4]; //cl_swiss_cross + // cl_vars[17] = showerShapes_[5]; //cl_nXtals + // cl_vars[18] = showerShapes_[6]; //cl_etaWidth + // cl_vars[19] = showerShapes_[7]; //cl_phiWidth + return cl_vars; +} + +std::vector EcalClustersGraph::computeWindowVariables(const std::vector>& clusters) { + size_t nCls = clusters.size(); + size_t nFeatures = clusters[0].size(); + std::vector min(nFeatures); + std::vector max(nFeatures); + std::vector sum(nFeatures); + for (const auto& vec : clusters) { + for (size_t i = 0; i < nFeatures; i++) { + const auto& x = vec[i]; + sum[i] += x; + if (x < min[i]) + min[i] = x; + if (x > max[i]) + max[i] = x; + } + } + std::vector out(18); + out[0] = max[0]; // max_en_cluster + out[1] = max[1]; // max_et_cluster + out[2] = max[7]; // max_deta_cluster + out[3] = max[8]; // max_dphi_cluster + out[4] = max[9]; // max_den + out[5] = max[10]; // max_det + out[6] = min[0]; // min_en_cluster + out[7] = min[1]; // min_et_cluster + out[8] = min[7]; // min_deta + out[9] = min[8]; // min_dphi + out[10] = min[9]; // min_den + out[11] = min[10]; // min_det + out[12] = sum[0] / nCls; // mean_en_cluster + out[13] = sum[1] / nCls; // mean_et_cluster + out[14] = sum[7] / nCls; // mean_deta + out[15] = sum[8] / nCls; // mean_dphi + out[16] = sum[9] / nCls; // mean_den + out[17] = sum[10] / nCls; // mean_det + return out; +} + +void EcalClustersGraph::fillVariables() { + + LogDebug("EcalClustersGraph") << "Preparing variables for all windows"; + + //Looping on all the seeds (window) + for (uint is = 0; is < nSeeds_; is++) { + uint nClsInWindow = 0; + const auto seedPointer = (*clusters_.at(is)).the_ptr().get(); + std::vector> unscaledClusterFeatures; + // Loop on all the clusters + for (uint ic = 0; ic < nCls_; ic++) { + if (inWindows_.Get(is, ic) == 1) { + const auto clPointer = (*clusters_.at(ic)).the_ptr().get(); + const auto & rawClX = computeVariables(seedPointer, clPointer); + unscaledClusterFeatures.push_back(rawClX); + inputs_.clustersX[is].push_back(SCProducerCache_->deepSCEvaluator->scaleClusterFeatures(rawClX)); + inputs_.hitsX[is].push_back(fillHits(clPointer)); + inputs_.isSeed[is].push_back(ic == is); + nClsInWindow++; + } + } + inputs_.windowX[is] = SCProducerCache_->deepSCEvaluator->scaleWindowFeatures(computeWindowVariables(unscaledClusterFeatures)); + + } + + inputs_.batchSize = nSeeds_; + + LogDebug("EcalClustersGraph") << "N. Windows: "<< inputs_.clustersX.size(); + + // LogDebug("EcalClustersGraph") << "Check hits: Seed | Cluster | Hits"; + // for (uint i = 0; i< nSeeds_;i++){ + // const size_t ncls = inputs_.hitsX[i].size(); + // for (size_t j = 0; jdeepSCEvaluator->evaluate(inputs_); + for (uint i = 0; i < nSeeds_; ++i){ + uint k = 0; + for (uint j = 0; j < nCls_; ++j) { + if (inWindows_.Get(i, j) == 1){ + scoreMatrix_.Set(i, j, scores[i][k]); + k++; + } + else{ + scoreMatrix_.Set(i, j, 0.); + } + } + } +} + +void EcalClustersGraph::printDebugInfo(){ + LogDebug("EcalClustersGraph") << "In window matrix:"; + for (uint i = 0; i < nSeeds_; ++i){ + for (uint j = 0; j < nCls_; ++j) { + std::cout << inWindows_.Get(i, j) << ","; + } + std::cout << std::endl; + } + LogDebug("EcalClustersGraph") << "Score matrix:"; + for (uint i = 0; i < nSeeds_; ++i){ + for (uint j = 0; j < nCls_; ++j) { + std::cout << scoreMatrix_.Get(i, j) << ","; + } + std::cout << std::endl; + } + LogDebug("EcalClustersGraph") << "Clusters ieta,iphi,iz,en"; + for (uint j = 0; j < nCls_; ++j) { + const auto cluster = (*clusters_.at(j)).the_ptr().get(); + std::vector clusterLocal = clusterPosition(cluster); + std::cout << clusterLocal[0] << "," << clusterLocal[1]<< "," << clusterLocal[2] << "," << cluster->energy() << std::endl; + } + +} + +void EcalClustersGraph::setThresholds() { + //test: place holder code + thresholds_ = std::vector(nSeeds_, 0.5); +} + +void EcalClustersGraph::selectClusters() { + //test + clusterMatrix_ = scoreMatrix_.ReduceElements(1, 1, thresholds_, false); + GraphMatrix clusterMatrixNoDuplicate_ = clusterMatrix_.RemoveDuplicates(1., false); + for (size_type r = 0; r < clusterMatrixNoDuplicate_.nRows(); r++) { + std::vector row = clusterMatrixNoDuplicate_.GetRow(r); + std::vector subRow(row.begin(), row.begin() + clusterMatrixNoDuplicate_.nRows()); + if (GraphMatrix().AllZeros(&subRow)) + clusterMatrix_.SetRowZero(r); + } + clusterMatrix_ = clusterMatrix_.RemoveDuplicates(1., false); + //std::cout << "clusterMatrix: " << clusterMatrix_ << std::endl; +} + +std::vector> EcalClustersGraph::getWindows() { + std::vector> windows; + for (size_type ir = 0; ir < clusterMatrix_.nRows(); ir++) { + if (GraphMatrix().AllZeros(clusterMatrix_.GetRow(ir))) + continue; + + CalibratedClusterPtr seed = clusters_[ir]; + CalibratedClusterPtrVector clusters_inWindow; + for (size_type ic = 0; ic < clusterMatrix_.nColumns(); ic++) + if (clusterMatrix_.Get(ir, ic) != 0.) + clusters_inWindow.push_back(clusters_[ic]); + windows.push_back(std::make_pair(seed, clusters_inWindow)); + } + return windows; +}