diff --git a/CommonTools/RecoAlgos/interface/MultiVectorManager.h b/CommonTools/RecoAlgos/interface/MultiVectorManager.h index 97e8a7285039c..90985db042446 100644 --- a/CommonTools/RecoAlgos/interface/MultiVectorManager.h +++ b/CommonTools/RecoAlgos/interface/MultiVectorManager.h @@ -53,6 +53,12 @@ class MultiVectorManager { class Iterator { public: + using iterator_category = std::forward_iterator_tag; + using difference_type = std::ptrdiff_t; + using value_type = T; + using pointer = T*; + using reference = T&; + Iterator(const MultiVectorManager& manager, size_t index) : manager(manager), currentIndex(index) {} bool operator!=(const Iterator& other) const { return currentIndex != other.currentIndex; } diff --git a/Configuration/ProcessModifiers/python/ticl_superclustering_dnn_cff.py b/Configuration/ProcessModifiers/python/ticl_superclustering_dnn_cff.py new file mode 100644 index 0000000000000..40cd46626a8b3 --- /dev/null +++ b/Configuration/ProcessModifiers/python/ticl_superclustering_dnn_cff.py @@ -0,0 +1,10 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5 +from Configuration.ProcessModifiers.ticl_superclustering_mustache_pf_cff import ticl_superclustering_mustache_pf +from Configuration.ProcessModifiers.ticl_superclustering_mustache_ticl_cff import ticl_superclustering_mustache_ticl + +# Modifier to run superclustering with DNN +# It is not to be used directly in a Process since it is currently the default for ticl_v5 +# It only exists as convenience in configuration files +ticl_superclustering_dnn = ticl_v5 & (~ticl_superclustering_mustache_pf) & (~ticl_superclustering_mustache_ticl) diff --git a/Configuration/ProcessModifiers/python/ticl_superclustering_mustache_pf_cff.py b/Configuration/ProcessModifiers/python/ticl_superclustering_mustache_pf_cff.py new file mode 100644 index 0000000000000..3487edbac62a9 --- /dev/null +++ b/Configuration/ProcessModifiers/python/ticl_superclustering_mustache_pf_cff.py @@ -0,0 +1,4 @@ +import FWCore.ParameterSet.Config as cms + +# Modifier (to be applied on top of ticl_v5 modifier) to use the pre-TICLv5 superclustering code (that translates tracksters to PFClusters before running Mustache) +ticl_superclustering_mustache_pf = cms.Modifier() diff --git a/Configuration/ProcessModifiers/python/ticl_superclustering_mustache_ticl_cff.py b/Configuration/ProcessModifiers/python/ticl_superclustering_mustache_ticl_cff.py new file mode 100644 index 0000000000000..8b2755f583007 --- /dev/null +++ b/Configuration/ProcessModifiers/python/ticl_superclustering_mustache_ticl_cff.py @@ -0,0 +1,4 @@ +import FWCore.ParameterSet.Config as cms + +# Modifier (to be applied on top of ticl_v5 modifier) to run superclustering using Mustache inside TICL (with tracksters as inputs) +ticl_superclustering_mustache_ticl = cms.Modifier() diff --git a/Configuration/PyReleaseValidation/README.md b/Configuration/PyReleaseValidation/README.md index cb04a44c6764f..3591fbb7283e9 100644 --- a/Configuration/PyReleaseValidation/README.md +++ b/Configuration/PyReleaseValidation/README.md @@ -96,6 +96,9 @@ The offsets currently in use are: * 0.103: Phase-2 aging, 3000fb-1 * 0.201: HGCAL special TICL Pattern recognition Workflows: clue3D * 0.202: HGCAL special TICL Pattern recognition Workflows: FastJet +* 0.203: HGCAL TICLv5 +* 0.204: HGCAL superclustering : using Mustache in TICLv5 +* 0.205: HGCAL superclustering : using old PFCluster-based Mustache algorithm with TICLv5 * 0.302: FastSim Run-3 trackingOnly validation * 0.303: FastSim Run-3 MB for mixing * 0.9001: Sonic Triton diff --git a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py index a7c41bf6b5f5e..becf6d18fd9a1 100644 --- a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py +++ b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py @@ -716,6 +716,44 @@ def condition(self, fragment, stepList, key, hasHarvest): upgradeWFs['ticl_v5'].step3 = {'--procModifiers': 'ticl_v5'} upgradeWFs['ticl_v5'].step4 = {'--procModifiers': 'ticl_v5'} +class UpgradeWorkflow_ticl_v5_superclustering(UpgradeWorkflow): + def setup_(self, step, stepName, stepDict, k, properties): + if 'RecoGlobal' in step: + stepDict[stepName][k] = merge([self.step3, stepDict[step][k]]) + if 'HARVESTGlobal' in step: + stepDict[stepName][k] = merge([self.step4, stepDict[step][k]]) + def condition(self, fragment, stepList, key, hasHarvest): + return (fragment=="ZEE_14" or 'Eta1p7_2p7' in fragment) and '2026' in key +upgradeWFs['ticl_v5_superclustering_mustache_ticl'] = UpgradeWorkflow_ticl_v5_superclustering( + steps = [ + 'RecoGlobal', + 'HARVESTGlobal' + ], + PU = [ + 'RecoGlobal', + 'HARVESTGlobal' + ], + suffix = '_ticl_v5_mustache', + offset = 0.204, +) +upgradeWFs['ticl_v5_superclustering_mustache_ticl'].step3 = {'--procModifiers': 'ticl_v5,ticl_superclustering_mustache_ticl'} +upgradeWFs['ticl_v5_superclustering_mustache_ticl'].step4 = {'--procModifiers': 'ticl_v5,ticl_superclustering_mustache_ticl'} + +upgradeWFs['ticl_v5_superclustering_mustache_pf'] = UpgradeWorkflow_ticl_v5_superclustering( + steps = [ + 'RecoGlobal', + 'HARVESTGlobal' + ], + PU = [ + 'RecoGlobal', + 'HARVESTGlobal' + ], + suffix = '_ticl_v5_mustache_pf', + offset = 0.205, +) +upgradeWFs['ticl_v5_superclustering_mustache_pf'].step3 = {'--procModifiers': 'ticl_v5,ticl_superclustering_mustache_pf'} +upgradeWFs['ticl_v5_superclustering_mustache_pf'].step4 = {'--procModifiers': 'ticl_v5,ticl_superclustering_mustache_pf'} + # Track DNN workflows class UpgradeWorkflow_trackdnn(UpgradeWorkflow): def setup_(self, step, stepName, stepDict, k, properties): diff --git a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py index 4952089002906..2ae2f9ffd8d9f 100644 --- a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py +++ b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py @@ -36,6 +36,7 @@ 'keep *_ticlTrackstersCLUE3DEM_*_*', 'keep *_ticlTrackstersCLUE3DHAD_*_*', 'keep *_ticlTracksterLinks_*_*', + 'keep *_ticlTracksterLinksSuperclustering_*_*', 'keep *_ticlCandidate_*_*', ] ) @@ -60,6 +61,8 @@ 'keep *_ticlSimTracksters_*_*', 'keep *_ticlSimTICLCandidates_*_*', 'keep *_ticlSimTrackstersFromCP_*_*', + 'keep *_tracksterSimTracksterAssociationLinkingSuperclustering_*_*', + 'keep *_tracksterSimTracksterAssociationPRSuperclustering_*_*', ) ) diff --git a/RecoHGCal/TICL/interface/SuperclusteringDNNInputs.h b/RecoHGCal/TICL/interface/SuperclusteringDNNInputs.h new file mode 100644 index 0000000000000..90ed787066022 --- /dev/null +++ b/RecoHGCal/TICL/interface/SuperclusteringDNNInputs.h @@ -0,0 +1,93 @@ +/** Computation of input features for superclustering DNN. Used by plugins/TracksterLinkingBySuperClustering.cc and plugins/SuperclusteringSampleDumper.cc */ +// Author: Theo Cuisset - theo.cuisset@cern.ch +// Date: 11/2023 + +#ifndef __RecoHGCal_TICL_SuperclusteringDNNInputs_H__ +#define __RecoHGCal_TICL_SuperclusteringDNNInputs_H__ + +#include +#include +#include + +namespace ticl { + class Trackster; + + // Abstract base class for DNN input preparation. + class AbstractSuperclusteringDNNInput { + public: + virtual ~AbstractSuperclusteringDNNInput() = default; + + virtual unsigned int featureCount() const { return featureNames().size(); }; + + /** Get name of features. Used for SuperclusteringSampleDumper branch names (inference does not use the names, only the indices) + * The default implementation is meant to be overriden by inheriting classes + */ + virtual std::vector featureNames() const { + std::vector defaultNames; + defaultNames.reserve(featureCount()); + for (unsigned int i = 1; i <= featureCount(); i++) { + defaultNames.push_back(std::string("nb_") + std::to_string(i)); + } + return defaultNames; + } + + /** Compute feature for seed and candidate pair */ + virtual std::vector computeVector(ticl::Trackster const& ts_base, ticl::Trackster const& ts_toCluster) = 0; + }; + + /* First version of DNN by Alessandro Tarabini. Meant as a DNN equivalent of Mustache algorithm (superclustering algo in ECAL) + Uses features : ['DeltaEta', 'DeltaPhi', 'multi_en', 'multi_eta', 'multi_pt', 'seedEta','seedPhi','seedEn', 'seedPt'] + */ + class SuperclusteringDNNInputV1 : public AbstractSuperclusteringDNNInput { + public: + unsigned int featureCount() const override { return 9; } + + std::vector computeVector(ticl::Trackster const& ts_base, ticl::Trackster const& ts_toCluster) override; + + std::vector featureNames() const override { + return {"DeltaEtaBaryc", + "DeltaPhiBaryc", + "multi_en", + "multi_eta", + "multi_pt", + "seedEta", + "seedPhi", + "seedEn", + "seedPt"}; + } + }; + + /* Second version of DNN by Alessandro Tarabini, making use of HGCAL-specific features. + Uses features : ['DeltaEta', 'DeltaPhi', 'multi_en', 'multi_eta', 'multi_pt', 'seedEta','seedPhi','seedEn', 'seedPt', 'theta', 'theta_xz_seedFrame', 'theta_yz_seedFrame', 'theta_xy_cmsFrame', 'theta_yz_cmsFrame', 'theta_xz_cmsFrame', 'explVar', 'explVarRatio'] + */ + class SuperclusteringDNNInputV2 : public AbstractSuperclusteringDNNInput { + public: + unsigned int featureCount() const override { return 17; } + + std::vector computeVector(ticl::Trackster const& ts_base, ticl::Trackster const& ts_toCluster) override; + + std::vector featureNames() const override { + return {"DeltaEtaBaryc", + "DeltaPhiBaryc", + "multi_en", + "multi_eta", + "multi_pt", + "seedEta", + "seedPhi", + "seedEn", + "seedPt", + "theta", + "theta_xz_seedFrame", + "theta_yz_seedFrame", + "theta_xy_cmsFrame", + "theta_yz_cmsFrame", + "theta_xz_cmsFrame", + "explVar", + "explVarRatio"}; + } + }; + + std::unique_ptr makeSuperclusteringDNNInputFromString(std::string dnnVersion); +} // namespace ticl + +#endif \ No newline at end of file diff --git a/RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h b/RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h index 2500cf54fa737..63da223cee6f6 100644 --- a/RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h +++ b/RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h @@ -30,11 +30,22 @@ namespace edm { class EventSetup; } // namespace edm +namespace cms { + namespace Ort { + class ONNXRuntime; + } +} // namespace cms + namespace ticl { class TracksterLinkingAlgoBase { public: - TracksterLinkingAlgoBase(const edm::ParameterSet& conf, edm::ConsumesCollector) - : algo_verbosity_(conf.getParameter("algo_verbosity")) {} + /** \param conf the configuration of the plugin + * \param onnxRuntime the ONNXRuntime, if onnxModelPath was provided in plugin configuration (nullptr otherwise) + */ + TracksterLinkingAlgoBase(const edm::ParameterSet& conf, + edm::ConsumesCollector, + cms::Ort::ONNXRuntime const* onnxRuntime = nullptr) + : algo_verbosity_(conf.getParameter("algo_verbosity")), onnxRuntime_(onnxRuntime) {} virtual ~TracksterLinkingAlgoBase(){}; struct Inputs { @@ -62,10 +73,14 @@ namespace ticl { const edm::ESHandle bfieldH, const edm::ESHandle propH) = 0; + // To be called by TracksterLinksProducer at the start of TracksterLinksProducer::produce. Subclasses can use this to store Event and EventSetup + virtual void setEvent(edm::Event& iEvent, edm::EventSetup const& iEventSetup){}; + static void fillPSetDescription(edm::ParameterSetDescription& desc) { desc.add("algo_verbosity", 0); }; protected: int algo_verbosity_; + cms::Ort::ONNXRuntime const* onnxRuntime_; }; } // namespace ticl diff --git a/RecoHGCal/TICL/plugins/BuildFile.xml b/RecoHGCal/TICL/plugins/BuildFile.xml index d48af46aa1572..7bd12c3249e67 100644 --- a/RecoHGCal/TICL/plugins/BuildFile.xml +++ b/RecoHGCal/TICL/plugins/BuildFile.xml @@ -23,6 +23,8 @@ + + diff --git a/RecoHGCal/TICL/plugins/EGammaSuperclusterProducer.cc b/RecoHGCal/TICL/plugins/EGammaSuperclusterProducer.cc new file mode 100644 index 0000000000000..d42bbf9fee24b --- /dev/null +++ b/RecoHGCal/TICL/plugins/EGammaSuperclusterProducer.cc @@ -0,0 +1,200 @@ +// Authors : Theo Cuisset , Shamik Ghosh +// Date : 01/2024 +/* +Translates TICL superclusters to ECAL supercluster dataformats (reco::SuperCluster and reco::CaloCluster). +Performs similar task as RecoEcal/EgammaCLusterProducers/PFECALSuperClusterProducer +Note that all tracksters are translated to reco::CaloCluster, even those that are not in any SuperCluster +*/ + +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Utilities/interface/FileInPath.h" + +#include +#include +#include +#include + +#include "PhysicsTools/ONNXRuntime/interface/ONNXRuntime.h" + +#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h" +#include "DataFormats/EgammaReco/interface/SuperCluster.h" +#include "DataFormats/EgammaReco/interface/ElectronSeed.h" +#include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" +#include "DataFormats/Math/interface/deltaPhi.h" + +using cms::Ort::ONNXRuntime; + +class EGammaSuperclusterProducer : public edm::stream::EDProducer> { +public: + EGammaSuperclusterProducer(const edm::ParameterSet&, const ONNXRuntime*); + + void produce(edm::Event&, const edm::EventSetup&) override; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + static std::unique_ptr initializeGlobalCache(const edm::ParameterSet& iConfig); + static void globalEndJob(const ONNXRuntime*); + +private: + edm::EDGetTokenT ticlSuperClustersToken_; + edm::EDGetTokenT>> superClusterLinksToken_; + edm::EDGetTokenT ticlTrackstersEMToken_; + edm::EDGetTokenT layerClustersToken_; + float superclusterEtThreshold_; + bool enableRegression_; +}; + +EGammaSuperclusterProducer::EGammaSuperclusterProducer(const edm::ParameterSet& ps, const ONNXRuntime*) + : ticlSuperClustersToken_(consumes(ps.getParameter("ticlSuperClusters"))), + superClusterLinksToken_(consumes>>( + edm::InputTag(ps.getParameter("ticlSuperClusters").label(), + "linkedTracksterIdToInputTracksterId", + ps.getParameter("ticlSuperClusters").process()))), + ticlTrackstersEMToken_(consumes(ps.getParameter("ticlTrackstersEM"))), + layerClustersToken_(consumes(ps.getParameter("layerClusters"))), + superclusterEtThreshold_(ps.getParameter("superclusterEtThreshold")), + enableRegression_(ps.getParameter("enableRegression")) { + produces(); + produces(); // The CaloCluster corresponding to each EM trackster +} + +void EGammaSuperclusterProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + auto const& ticlSuperclusters = iEvent.get(ticlSuperClustersToken_); + auto const& ticlSuperclusterLinks = iEvent.get(superClusterLinksToken_); + auto emTracksters_h = iEvent.getHandle(ticlTrackstersEMToken_); + auto const& emTracksters = *emTracksters_h; + auto const& layerClusters = iEvent.get(layerClustersToken_); + // Output collections : + auto egammaSuperclusters = std::make_unique(); + auto caloClustersEM = std::make_unique(); + + // Fill reco::CaloCluster collection (1-1 mapping to TICL EM tracksters) + for (ticl::Trackster const& emTrackster : emTracksters) { + std::vector> hitsAndFractions; + int iLC = 0; + std::for_each(std::begin(emTrackster.vertices()), std::end(emTrackster.vertices()), [&](unsigned int lcId) { + const auto fraction = 1.f / emTrackster.vertex_multiplicity(iLC++); + for (const auto& cell : layerClusters[lcId].hitsAndFractions()) { + hitsAndFractions.emplace_back(cell.first, cell.second * fraction); + } + }); + + reco::CaloCluster& caloCluster = caloClustersEM->emplace_back( + emTrackster.raw_energy(), // energy + math::XYZPoint(emTrackster.barycenter()), // position + reco::CaloID(reco::CaloID::DET_HGCAL_ENDCAP), + hitsAndFractions, + reco::CaloCluster::particleFlow, // algoID (copying from output of PFECALCSuperClusterProducer) + hitsAndFractions.at(0) + .first // seedId (this may need to be updated once a common definition of the seed of a cluster is adopted for Phase-2) + ); + caloCluster.setCorrectedEnergy(emTrackster.raw_energy()); // Needs to be updated with new supercluster regression + } + + edm::OrphanHandle caloClustersEM_h = iEvent.put(std::move(caloClustersEM)); + + // Fill reco::SuperCluster collection and prepare regression inputs + assert(ticlSuperclusters.size() == ticlSuperclusterLinks.size()); + const unsigned int regressionFeatureCount = 8; + std::vector regressionInputs; + regressionInputs.reserve(ticlSuperclusters.size() * regressionFeatureCount); + unsigned int superclustersPassingSelectionsCount = 0; + for (std::size_t sc_i = 0; sc_i < ticlSuperclusters.size(); sc_i++) { + ticl::Trackster const& ticlSupercluster = ticlSuperclusters[sc_i]; + if (ticlSupercluster.raw_pt() < superclusterEtThreshold_) + continue; + std::vector const& superclusterLink = ticlSuperclusterLinks[sc_i]; + assert(!superclusterLink.empty()); + + reco::CaloClusterPtrVector trackstersEMInSupercluster; + float max_eta = std::numeric_limits::min(); + float max_phi = std::numeric_limits::min(); + float min_eta = std::numeric_limits::max(); + float min_phi = std::numeric_limits::max(); + for (unsigned int tsInSc_id : superclusterLink) { + trackstersEMInSupercluster.push_back(reco::CaloClusterPtr(caloClustersEM_h, tsInSc_id)); + ticl::Trackster const& constituentTrackster = emTracksters[tsInSc_id]; + max_eta = std::max(constituentTrackster.barycenter().eta(), max_eta); + max_phi = std::max(constituentTrackster.barycenter().phi(), max_phi); + min_eta = std::min(constituentTrackster.barycenter().eta(), min_eta); + min_phi = std::min(constituentTrackster.barycenter().phi(), min_phi); + } + egammaSuperclusters->emplace_back( + ticlSupercluster.raw_energy(), + reco::SuperCluster::Point(ticlSupercluster.barycenter()), + reco::CaloClusterPtr(caloClustersEM_h, + superclusterLink[0]), // seed (first trackster in superclusterLink is the seed) + trackstersEMInSupercluster, // clusters + 0., // Epreshower (not relevant for HGCAL) + -1, // phiwidth (not implemented yet) + -1 // etawidth (not implemented yet) + ); + superclustersPassingSelectionsCount++; + + if (enableRegression_) { + regressionInputs.insert( + regressionInputs.end(), + {ticlSupercluster.barycenter().eta(), + ticlSupercluster.barycenter().phi(), + ticlSupercluster.raw_energy(), + std::abs(max_eta - min_eta), + max_phi - min_phi > M_PI ? 2 * static_cast(M_PI) - (max_phi - min_phi) : max_phi - min_phi, + emTracksters[superclusterLink[0]].raw_energy() - + (superclusterLink.size() >= 2 ? emTracksters[superclusterLink.back()].raw_energy() : 0.f), + emTracksters[superclusterLink[0]].raw_energy() / ticlSupercluster.raw_energy(), + static_cast(superclusterLink.size())}); + } + } + + if (enableRegression_ && superclustersPassingSelectionsCount > 0) { + // Run the regression + // ONNXRuntime takes std::vector>& as input (non-const reference) so we have to make a new vector + std::vector> inputs_for_onnx{{std::move(regressionInputs)}}; + std::vector outputs = + globalCache()->run({"input"}, inputs_for_onnx, {}, {}, superclustersPassingSelectionsCount)[0]; + + assert(egammaSuperclusters->size() == outputs.size()); + for (std::size_t sc_i = 0; sc_i < egammaSuperclusters->size(); sc_i++) { + (*egammaSuperclusters)[sc_i].setCorrectedEnergy(outputs[sc_i]); + // correctedEnergyUncertainty is left at its default value + } + } + + iEvent.put(std::move(egammaSuperclusters)); +} + +std::unique_ptr EGammaSuperclusterProducer::initializeGlobalCache(const edm::ParameterSet& iConfig) { + if (iConfig.getParameter("enableRegression")) + return std::make_unique(iConfig.getParameter("regressionModelPath").fullPath()); + else + return std::unique_ptr(nullptr); +} + +void EGammaSuperclusterProducer::globalEndJob(const ONNXRuntime*) {} + +void EGammaSuperclusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("ticlSuperClusters", edm::InputTag("ticlTracksterLinksSuperclusteringDNN")); + desc.add("ticlTrackstersEM", edm::InputTag("ticlTrackstersCLUE3DHigh")) + ->setComment("The trackster collection used before superclustering, ie CLUE3D EM tracksters"); + desc.add("layerClusters", edm::InputTag("hgcalMergeLayerClusters")) + ->setComment("The layer cluster collection that goes with ticlTrackstersEM"); + desc.add("superclusterEtThreshold", 4.)->setComment("Minimum supercluster transverse energy."); + desc.add("enableRegression", true)->setComment("Enable supercluster energy regression"); + desc.add("regressionModelPath", + edm::FileInPath("RecoHGCal/TICL/data/superclustering/regression_v1.onnx")) + ->setComment("Path to regression network (as ONNX model)"); + + descriptions.add("ticlEGammaSuperClusterProducer", desc); +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(EGammaSuperclusterProducer); diff --git a/RecoHGCal/TICL/plugins/SuperclusteringSampleDumper.cc b/RecoHGCal/TICL/plugins/SuperclusteringSampleDumper.cc new file mode 100644 index 0000000000000..7208a3c923a86 --- /dev/null +++ b/RecoHGCal/TICL/plugins/SuperclusteringSampleDumper.cc @@ -0,0 +1,298 @@ +// Original Author: Theo Cuisset +// Created: Nov 2023 +/** Produce samples for electron superclustering DNN training in TICL + + Pairs of seed-candidate tracksters (in compatible geometric windows) are iterated over, in similar manner as in TracksterLinkingBySuperclustering. + For each of these pairs, the DNN features are computed and saved to a TTree. + Also saved is the best (=lowest) association score of the seed trackster with CaloParticles. The association score of the candidate trackster + with the same CaloParticle is also saved. +*/ +#include +#include +#include +#include + +#include + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/ParameterSet/interface/allowedValues.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/ServiceRegistry/interface/Service.h" + +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +#include "DataFormats/Math/interface/deltaPhi.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" +#include "SimDataFormats/Associations/interface/TracksterToSimTracksterAssociator.h" + +#include "RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringDNN.h" +#include "RecoHGCal/TICL/interface/SuperclusteringDNNInputs.h" + +using namespace ticl; + +class SuperclusteringSampleDumper : public edm::one::EDAnalyzer { +public: + explicit SuperclusteringSampleDumper(const edm::ParameterSet&); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void beginJob() override; + void analyze(const edm::Event&, const edm::EventSetup&) override; + bool checkExplainedVarianceRatioCut(ticl::Trackster const& ts) const; + + const edm::EDGetTokenT> tracksters_clue3d_token_; + const edm::EDGetTokenT tsRecoToSimCP_token_; + float deltaEtaWindow_; + float deltaPhiWindow_; + float seedPtThreshold_; + float candidateEnergyThreshold_; + float explVarRatioCut_energyBoundary_; // Boundary energy between low and high energy explVarRatio cut threshold + float explVarRatioMinimum_lowEnergy_; // Cut on explained variance ratio of tracksters to be considered as candidate, for trackster raw_energy < explVarRatioCut_energyBoundary + float explVarRatioMinimum_highEnergy_; // Cut on explained variance ratio of tracksters to be considered as candidate, for trackster raw_energy > explVarRatioCut_energyBoundary + + TTree* output_tree_; + edm::EventID eventId_; + std::unique_ptr dnnInput_; + std::vector> + features_; // Outer index : feature number (split into branches), inner index : inference pair index + std::vector seedTracksterIdx_; // ID of seed trackster used for inference pair + std::vector candidateTracksterIdx_; // ID of candidate trackster used for inference pair + + std::vector + seedTracksterBestAssociationScore_; // Best association score of seed trackster (seedTracksterIdx) with CaloParticle + std::vector + seedTracksterBestAssociation_simTsIdx_; // Index of SimTrackster that has the best association score to the seedTrackster + std::vector seedTracksterBestAssociation_caloParticleEnergy_; // Energy of best associated CaloParticle to seed + + std::vector + candidateTracksterBestAssociationScore_; // Best association score of candidate trackster (seedTracksterIdx) with CaloParticle + std::vector + candidateTracksterBestAssociation_simTsIdx_; // Index of SimTrackster that has the best association score to the candidate + + std::vector + candidateTracksterAssociationWithSeed_score_; // Association score of candidate trackster with the CaloParticle of seedTracksterBestAssociation_simTsIdx_ +}; + +SuperclusteringSampleDumper::SuperclusteringSampleDumper(const edm::ParameterSet& ps) + : tracksters_clue3d_token_(consumes>(ps.getParameter("tracksters"))), + tsRecoToSimCP_token_( + consumes(ps.getParameter("recoToSimAssociatorCP"))), + deltaEtaWindow_(ps.getParameter("deltaEtaWindow")), + deltaPhiWindow_(ps.getParameter("deltaPhiWindow")), + seedPtThreshold_(ps.getParameter("seedPtThreshold")), + candidateEnergyThreshold_(ps.getParameter("candidateEnergyThreshold")), + explVarRatioCut_energyBoundary_(ps.getParameter("candidateEnergyThreshold")), + explVarRatioMinimum_lowEnergy_(ps.getParameter("explVarRatioMinimum_lowEnergy")), + explVarRatioMinimum_highEnergy_(ps.getParameter("explVarRatioMinimum_highEnergy")), + eventId_(), + dnnInput_(makeSuperclusteringDNNInputFromString(ps.getParameter("dnnInputsVersion"))), + features_(dnnInput_->featureCount()) { + usesResource("TFileService"); +} + +void SuperclusteringSampleDumper::beginJob() { + edm::Service fs; + output_tree_ = fs->make("superclusteringTraining", "Superclustering training samples"); + output_tree_->Branch("Event", &eventId_); + output_tree_->Branch("seedTracksterIdx", &seedTracksterIdx_); + output_tree_->Branch("candidateTracksterIdx", &candidateTracksterIdx_); + output_tree_->Branch("seedTracksterBestAssociationScore", &seedTracksterBestAssociationScore_); + output_tree_->Branch("seedTracksterBestAssociation_simTsIdx", &seedTracksterBestAssociation_simTsIdx_); + output_tree_->Branch("seedTracksterBestAssociation_caloParticleEnergy", + &seedTracksterBestAssociation_caloParticleEnergy_); + output_tree_->Branch("candidateTracksterBestAssociationScore", &candidateTracksterBestAssociationScore_); + output_tree_->Branch("candidateTracksterBestAssociation_simTsIdx", &candidateTracksterBestAssociation_simTsIdx_); + output_tree_->Branch("candidateTracksterAssociationWithSeed_score", &candidateTracksterAssociationWithSeed_score_); + std::vector featureNames = dnnInput_->featureNames(); + assert(featureNames.size() == dnnInput_->featureCount()); + for (unsigned int i = 0; i < dnnInput_->featureCount(); i++) { + output_tree_->Branch(("feature_" + featureNames[i]).c_str(), &features_[i]); + } +} + +/** + * Check if trackster passes cut on explained variance ratio. The DNN is trained only on pairs where both seed and candidate pass this cut + * Explained variance ratio is (largest PCA eigenvalue) / (sum of PCA eigenvalues) +*/ +bool SuperclusteringSampleDumper::checkExplainedVarianceRatioCut(ticl::Trackster const& ts) const { + float explVar_denominator = + std::accumulate(std::begin(ts.eigenvalues()), std::end(ts.eigenvalues()), 0.f, std::plus()); + if (explVar_denominator != 0.) { + float explVarRatio = ts.eigenvalues()[0] / explVar_denominator; + if (ts.raw_energy() > explVarRatioCut_energyBoundary_) + return explVarRatio >= explVarRatioMinimum_highEnergy_; + else + return explVarRatio >= explVarRatioMinimum_lowEnergy_; + } else + return false; +} + +void SuperclusteringSampleDumper::analyze(const edm::Event& evt, const edm::EventSetup& iSetup) { + eventId_ = evt.id(); + + edm::Handle> inputTracksters; + evt.getByToken(tracksters_clue3d_token_, inputTracksters); + + edm::Handle assoc_CP_recoToSim; + evt.getByToken(tsRecoToSimCP_token_, assoc_CP_recoToSim); + + /* Sorting tracksters by decreasing order of pT (out-of-place sort). + inputTracksters[trackstersIndicesPt[0]], ..., inputTracksters[trackstersIndicesPt[N]] makes a list of tracksters sorted by decreasing pT + Indices into this pT sorted collection will have the suffix _pt. Thus inputTracksters[index] and inputTracksters[trackstersIndicesPt[index_pt]] are correct + */ + std::vector trackstersIndicesPt(inputTracksters->size()); + std::iota(trackstersIndicesPt.begin(), trackstersIndicesPt.end(), 0); + std::stable_sort( + trackstersIndicesPt.begin(), trackstersIndicesPt.end(), [&inputTracksters](unsigned int i1, unsigned int i2) { + return (*inputTracksters)[i1].raw_pt() > (*inputTracksters)[i2].raw_pt(); + }); + + // Order of loops are reversed compared to SuperclusteringProducer (here outer is seed, inner is candidate), for performance reasons. + // The same pairs seed-candidate should be present, just in a different order + // First loop on seed tracksters + for (unsigned int ts_seed_idx_pt = 0; ts_seed_idx_pt < inputTracksters->size(); ts_seed_idx_pt++) { + const unsigned int ts_seed_idx_input = + trackstersIndicesPt[ts_seed_idx_pt]; // Index of seed trackster in input collection (not in pT sorted collection) + Trackster const& ts_seed = (*inputTracksters)[ts_seed_idx_input]; + + if (ts_seed.raw_pt() < seedPtThreshold_) + break; // All further seeds will have lower pT than threshold (due to pT sorting) + + if (!checkExplainedVarianceRatioCut(ts_seed)) + continue; + + // Find best associated CaloParticle to the seed + auto seed_assocs = assoc_CP_recoToSim->find({inputTracksters, ts_seed_idx_input}); + if (seed_assocs == assoc_CP_recoToSim->end()) + continue; // No CaloParticle associations for the current trackster (extremly unlikely) + ticl::RecoToSimCollectionSimTracksters::data_type const& seed_assocWithBestScore = + *std::min_element(seed_assocs->val.begin(), + seed_assocs->val.end(), + [](ticl::RecoToSimCollectionSimTracksters::data_type const& assoc_1, + ticl::RecoToSimCollectionSimTracksters::data_type const& assoc_2) { + // assoc_* is of type : std::pair> + // Best score is smallest score + return assoc_1.second.second < assoc_2.second.second; + }); + + // Second loop on superclustering candidates tracksters + // Look only at candidate tracksters with lower pT than the seed (so all pairs are only looked at once) + for (unsigned int ts_cand_idx_pt = ts_seed_idx_pt + 1; ts_cand_idx_pt < inputTracksters->size(); ts_cand_idx_pt++) { + Trackster const& ts_cand = (*inputTracksters)[trackstersIndicesPt[ts_cand_idx_pt]]; + // Check that the two tracksters are geometrically compatible for superclustering (using deltaEta, deltaPhi window) + // There is no need to run training or inference for tracksters very far apart + if (!(std::abs(ts_seed.barycenter().Eta() - ts_cand.barycenter().Eta()) < deltaEtaWindow_ && + std::abs(deltaPhi(ts_seed.barycenter().Phi(), ts_cand.barycenter().Phi())) < deltaPhiWindow_ && + ts_cand.raw_energy() >= candidateEnergyThreshold_ && checkExplainedVarianceRatioCut(ts_cand))) + continue; + + std::vector features = dnnInput_->computeVector(ts_seed, ts_cand); + assert(features.size() == features_.size()); + for (unsigned int feature_idx = 0; feature_idx < features_.size(); feature_idx++) { + features_[feature_idx].push_back(features[feature_idx]); + } + seedTracksterIdx_.push_back(trackstersIndicesPt[ts_seed_idx_pt]); + candidateTracksterIdx_.push_back(trackstersIndicesPt[ts_cand_idx_pt]); + + float candidateTracksterBestAssociationScore = 1.; // Best association score of candidate with any CaloParticle + long candidateTracksterBestAssociation_simTsIdx = -1; // Corresponding CaloParticle simTrackster index + float candidateTracksterAssociationWithSeed_score = + 1.; // Association score of candidate with CaloParticle best associated with seed + + // First find associated CaloParticles with candidate + auto cand_assocCP = assoc_CP_recoToSim->find( + edm::Ref(inputTracksters, trackstersIndicesPt[ts_cand_idx_pt])); + if (cand_assocCP != assoc_CP_recoToSim->end()) { + // find the association with best score + ticl::RecoToSimCollectionSimTracksters::data_type const& cand_assocWithBestScore = + *std::min_element(cand_assocCP->val.begin(), + cand_assocCP->val.end(), + [](ticl::RecoToSimCollectionSimTracksters::data_type const& assoc_1, + ticl::RecoToSimCollectionSimTracksters::data_type const& assoc_2) { + // assoc_* is of type : std::pair> + return assoc_1.second.second < assoc_2.second.second; + }); + candidateTracksterBestAssociationScore = cand_assocWithBestScore.second.second; + candidateTracksterBestAssociation_simTsIdx = cand_assocWithBestScore.first.key(); + + // find the association score with the same CaloParticle as the seed + auto cand_assocWithSeedCP = + std::find_if(cand_assocCP->val.begin(), + cand_assocCP->val.end(), + [&seed_assocWithBestScore](ticl::RecoToSimCollectionSimTracksters::data_type const& assoc) { + // assoc is of type : std::pair> + return assoc.first == seed_assocWithBestScore.first; + }); + if (cand_assocWithSeedCP != cand_assocCP->val.end()) { + candidateTracksterAssociationWithSeed_score = cand_assocWithSeedCP->second.second; + } + } + + seedTracksterBestAssociationScore_.push_back(seed_assocWithBestScore.second.second); + seedTracksterBestAssociation_simTsIdx_.push_back(seed_assocWithBestScore.first.key()); + seedTracksterBestAssociation_caloParticleEnergy_.push_back(seed_assocWithBestScore.first->regressed_energy()); + + candidateTracksterBestAssociationScore_.push_back(candidateTracksterBestAssociationScore); + candidateTracksterBestAssociation_simTsIdx_.push_back(candidateTracksterBestAssociation_simTsIdx); + + candidateTracksterAssociationWithSeed_score_.push_back(candidateTracksterAssociationWithSeed_score); + } + } + + output_tree_->Fill(); + for (auto& feats : features_) + feats.clear(); + seedTracksterIdx_.clear(); + candidateTracksterIdx_.clear(); + seedTracksterBestAssociationScore_.clear(); + seedTracksterBestAssociation_simTsIdx_.clear(); + seedTracksterBestAssociation_caloParticleEnergy_.clear(); + candidateTracksterBestAssociationScore_.clear(); + candidateTracksterBestAssociation_simTsIdx_.clear(); + candidateTracksterAssociationWithSeed_score_.clear(); +} + +void SuperclusteringSampleDumper::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("tracksters", edm::InputTag("ticlTrackstersCLUE3DHigh")) + ->setComment("Input trackster collection, same as what is used for superclustering inference."); + desc.add("recoToSimAssociatorCP", + edm::InputTag("tracksterSimTracksterAssociationLinkingbyCLUE3D", "recoToSim")); + desc.ifValue(edm::ParameterDescription("dnnInputsVersion", "v2", true), + edm::allowedValues("v1", "v2")) + ->setComment( + "DNN inputs version tag. Defines which set of features is fed to the DNN. Must match with the actual DNN."); + // Cuts are intentionally looser than those used for inference in TracksterLinkingBySuperClustering.cpp + desc.add("deltaEtaWindow", 0.2) + ->setComment( + "Size of delta eta window to consider for superclustering. Seed-candidate pairs outside this window " + "are not considered for DNN inference."); + desc.add("deltaPhiWindow", 0.7) + ->setComment( + "Size of delta phi window to consider for superclustering. Seed-candidate pairs outside this window " + "are not considered for DNN inference."); + desc.add("seedPtThreshold", 3.) + ->setComment("Minimum transverse momentum of trackster to be considered as seed of a supercluster"); + desc.add("candidateEnergyThreshold", 1.5) + ->setComment("Minimum energy of trackster to be considered as candidate for superclustering"); + desc.add("explVarRatioCut_energyBoundary", 50.) + ->setComment("Boundary energy between low and high energy explVarRatio cut threshold"); + desc.add("explVarRatioMinimum_lowEnergy", 0.85) + ->setComment( + "Cut on explained variance ratio of tracksters to be considered as candidate, " + "for trackster raw_energy < explVarRatioCut_energyBoundary"); + desc.add("explVarRatioMinimum_highEnergy", 0.9) + ->setComment( + "Cut on explained variance ratio of tracksters to be considered as candidate, " + "for trackster raw_energy > explVarRatioCut_energyBoundary"); + descriptions.add("superclusteringSampleDumper", desc); +} + +DEFINE_FWK_MODULE(SuperclusteringSampleDumper); \ No newline at end of file diff --git a/RecoHGCal/TICL/plugins/TICLCandidateProducer.cc b/RecoHGCal/TICL/plugins/TICLCandidateProducer.cc index 67c382988cd15..669e4c8e0a96a 100644 --- a/RecoHGCal/TICL/plugins/TICLCandidateProducer.cc +++ b/RecoHGCal/TICL/plugins/TICLCandidateProducer.cc @@ -307,8 +307,12 @@ void TICLCandidateProducer::produce(edm::Event &evt, const edm::EventSetup &es) generalInterpretationAlgo_->makeCandidates(input, inputTiming_h, *resultTracksters, trackstersInTrackIndices); - assignPCAtoTracksters( - *resultTracksters, layerClusters, layerClustersTimes, rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z(), rhtools_, true); + assignPCAtoTracksters(*resultTracksters, + layerClusters, + layerClustersTimes, + rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z(), + rhtools_, + true); std::vector maskTracksters(resultTracksters->size(), true); edm::OrphanHandle> resultTracksters_h = evt.put(std::move(resultTracksters)); diff --git a/RecoHGCal/TICL/plugins/TICLDumper.cc b/RecoHGCal/TICL/plugins/TICLDumper.cc index e4d2a05b068f2..36a82346ed524 100644 --- a/RecoHGCal/TICL/plugins/TICLDumper.cc +++ b/RecoHGCal/TICL/plugins/TICLDumper.cc @@ -11,6 +11,7 @@ #include // unique_ptr #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" #include "FWCore/ParameterSet/interface/PluginDescription.h" +#include "FWCore/ParameterSet/interface/allowedValues.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/Utilities/interface/ESGetToken.h" #include "FWCore/Framework/interface/ESHandle.h" @@ -21,6 +22,7 @@ #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/ServiceRegistry/interface/Service.h" +#include "DataFormats/Provenance/interface/EventID.h" #include "DataFormats/CaloRecHit/interface/CaloCluster.h" #include "DataFormats/HGCalReco/interface/Trackster.h" #include "DataFormats/HGCalReco/interface/TICLCandidate.h" @@ -33,6 +35,8 @@ #include "DataFormats/HGCalReco/interface/Common.h" #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" +#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h" +#include "DataFormats/EgammaReco/interface/SuperCluster.h" #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h" @@ -57,6 +61,505 @@ #include "FWCore/ServiceRegistry/interface/Service.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" +// Helper class for geometry, magnetic field, etc +class DetectorTools { +public: + DetectorTools(const HGCalDDDConstants& hgcons, + const CaloGeometry& geom, + const MagneticField& bfieldH, + const Propagator& propH) + : hgcons(hgcons), rhtools(), bfield(bfieldH), propagator(propH) { + rhtools.setGeometry(geom); + + // build disks at HGCal front & EM-Had interface for track propagation + float zVal = hgcons.waferZ(1, true); + std::pair rMinMax = hgcons.rangeR(zVal, true); + + float zVal_interface = rhtools.getPositionLayer(rhtools.lastLayerEE()).z(); + std::pair rMinMax_interface = hgcons.rangeR(zVal_interface, true); + + for (int iSide = 0; iSide < 2; ++iSide) { + float zSide = (iSide == 0) ? (-1. * zVal) : zVal; + firstDisk_[iSide] = std::make_unique( + Disk::build(Disk::PositionType(0, 0, zSide), + Disk::RotationType(), + SimpleDiskBounds(rMinMax.first, rMinMax.second, zSide - 0.5, zSide + 0.5)) + .get()); + + zSide = (iSide == 0) ? (-1. * zVal_interface) : zVal_interface; + interfaceDisk_[iSide] = std::make_unique( + Disk::build(Disk::PositionType(0, 0, zSide), + Disk::RotationType(), + SimpleDiskBounds(rMinMax_interface.first, rMinMax_interface.second, zSide - 0.5, zSide + 0.5)) + .get()); + } + } + + const HGCalDDDConstants& hgcons; + std::unique_ptr firstDisk_[2]; + std::unique_ptr interfaceDisk_[2]; + hgcal::RecHitTools rhtools; + const MagneticField& bfield; + const Propagator& propagator; +}; + +// Helper class that dumps a single trackster collection (either tracksters or simTracksters) +class TracksterDumperHelper { +public: + enum class TracksterType { + Trackster, ///< Regular trackster (from RECO) + SimTracksterCP, ///< SimTrackster from CaloParticle + SimTracksterSC ///< SimTrackster from SimCluster + }; + + static TracksterType tracksterTypeFromString(std::string str) { + if (str == "Trackster") + return TracksterType::Trackster; + if (str == "SimTracksterCP") + return TracksterType::SimTracksterCP; + if (str == "SimTracksterSC") + return TracksterType::SimTracksterSC; + throw std::runtime_error("TICLDumper : TracksterDumperHelper : Invalid trackster type " + str); + } + + /** tracksterType : dtermines additional information that will be saved (calo truth information, track information) */ + TracksterDumperHelper(TracksterType tracksterType = TracksterType::Trackster) : tracksterType_(tracksterType) {} + + /** + * To be called once after tree creation. eventId_ should be a pointer to the EventID. + * *Do not copy/move or resize vector holding object after calling this function* + */ + void initTree(TTree* trackster_tree_, edm::EventID* eventId_) { + trackster_tree_->Branch("event", eventId_); + trackster_tree_->Branch("NTracksters", &nTracksters); + trackster_tree_->Branch("NClusters", &nClusters); + trackster_tree_->Branch("time", &trackster_time); + trackster_tree_->Branch("timeError", &trackster_timeError); + trackster_tree_->Branch("regressed_energy", &trackster_regressed_energy); + trackster_tree_->Branch("raw_energy", &trackster_raw_energy); + trackster_tree_->Branch("raw_em_energy", &trackster_raw_em_energy); + trackster_tree_->Branch("raw_pt", &trackster_raw_pt); + trackster_tree_->Branch("raw_em_pt", &trackster_raw_em_pt); + trackster_tree_->Branch("barycenter_x", &trackster_barycenter_x); + trackster_tree_->Branch("barycenter_y", &trackster_barycenter_y); + trackster_tree_->Branch("barycenter_z", &trackster_barycenter_z); + trackster_tree_->Branch("barycenter_eta", &trackster_barycenter_eta); + trackster_tree_->Branch("barycenter_phi", &trackster_barycenter_phi); + trackster_tree_->Branch("EV1", &trackster_EV1); + trackster_tree_->Branch("EV2", &trackster_EV2); + trackster_tree_->Branch("EV3", &trackster_EV3); + trackster_tree_->Branch("eVector0_x", &trackster_eVector0_x); + trackster_tree_->Branch("eVector0_y", &trackster_eVector0_y); + trackster_tree_->Branch("eVector0_z", &trackster_eVector0_z); + trackster_tree_->Branch("sigmaPCA1", &trackster_sigmaPCA1); + trackster_tree_->Branch("sigmaPCA2", &trackster_sigmaPCA2); + trackster_tree_->Branch("sigmaPCA3", &trackster_sigmaPCA3); + if (tracksterType_ != TracksterType::Trackster) { + trackster_tree_->Branch("regressed_pt", &simtrackster_regressed_pt); + trackster_tree_->Branch("pdgID", &simtrackster_pdgID); + trackster_tree_->Branch("trackIdx", &simtrackster_trackIdx); + trackster_tree_->Branch("trackTime", &simtrackster_trackTime); + trackster_tree_->Branch("timeBoundary", &simtrackster_timeBoundary); + trackster_tree_->Branch("boundaryX", &simtrackster_boundaryX); + trackster_tree_->Branch("boundaryY", &simtrackster_boundaryY); + trackster_tree_->Branch("boundaryZ", &simtrackster_boundaryZ); + trackster_tree_->Branch("boundaryEta", &simtrackster_boundaryEta); + trackster_tree_->Branch("boundaryPhi", &simtrackster_boundaryPhi); + trackster_tree_->Branch("boundaryPx", &simtrackster_boundaryPx); + trackster_tree_->Branch("boundaryPy", &simtrackster_boundaryPy); + trackster_tree_->Branch("boundaryPz", &simtrackster_boundaryPz); + trackster_tree_->Branch("track_boundaryX", &simtrackster_track_boundaryX); + trackster_tree_->Branch("track_boundaryY", &simtrackster_track_boundaryY); + trackster_tree_->Branch("track_boundaryZ", &simtrackster_track_boundaryZ); + trackster_tree_->Branch("track_boundaryEta", &simtrackster_track_boundaryEta); + trackster_tree_->Branch("track_boundaryPhi", &simtrackster_track_boundaryPhi); + trackster_tree_->Branch("track_boundaryPx", &simtrackster_track_boundaryPx); + trackster_tree_->Branch("track_boundaryPy", &simtrackster_track_boundaryPy); + trackster_tree_->Branch("track_boundaryPz", &simtrackster_track_boundaryPz); + } + trackster_tree_->Branch("id_probabilities", &trackster_id_probabilities); + trackster_tree_->Branch("vertices_indexes", &trackster_vertices_indexes); + trackster_tree_->Branch("vertices_x", &trackster_vertices_x); + trackster_tree_->Branch("vertices_y", &trackster_vertices_y); + trackster_tree_->Branch("vertices_z", &trackster_vertices_z); + trackster_tree_->Branch("vertices_time", &trackster_vertices_time); + trackster_tree_->Branch("vertices_timeErr", &trackster_vertices_timeErr); + trackster_tree_->Branch("vertices_energy", &trackster_vertices_energy); + trackster_tree_->Branch("vertices_correctedEnergy", &trackster_vertices_correctedEnergy); + trackster_tree_->Branch("vertices_correctedEnergyUncertainty", &trackster_vertices_correctedEnergyUncertainty); + trackster_tree_->Branch("vertices_multiplicity", &trackster_vertices_multiplicity); + } + + void clearVariables() { + nTracksters = 0; + nClusters = 0; + trackster_time.clear(); + trackster_timeError.clear(); + trackster_regressed_energy.clear(); + trackster_raw_energy.clear(); + trackster_raw_em_energy.clear(); + trackster_raw_pt.clear(); + trackster_raw_em_pt.clear(); + trackster_barycenter_x.clear(); + trackster_barycenter_y.clear(); + trackster_barycenter_z.clear(); + trackster_EV1.clear(); + trackster_EV2.clear(); + trackster_EV3.clear(); + trackster_eVector0_x.clear(); + trackster_eVector0_y.clear(); + trackster_eVector0_z.clear(); + trackster_sigmaPCA1.clear(); + trackster_sigmaPCA2.clear(); + trackster_sigmaPCA3.clear(); + + simtrackster_regressed_pt.clear(); + simtrackster_pdgID.clear(); + simtrackster_trackIdx.clear(); + simtrackster_trackTime.clear(); + simtrackster_timeBoundary.clear(); + simtrackster_boundaryX.clear(); + simtrackster_boundaryY.clear(); + simtrackster_boundaryZ.clear(); + simtrackster_boundaryEta.clear(); + simtrackster_boundaryPhi.clear(); + simtrackster_boundaryPx.clear(); + simtrackster_boundaryPy.clear(); + simtrackster_boundaryPz.clear(); + simtrackster_track_boundaryX.clear(); + simtrackster_track_boundaryY.clear(); + simtrackster_track_boundaryZ.clear(); + simtrackster_track_boundaryEta.clear(); + simtrackster_track_boundaryPhi.clear(); + simtrackster_track_boundaryPx.clear(); + simtrackster_track_boundaryPy.clear(); + simtrackster_track_boundaryPz.clear(); + + trackster_barycenter_eta.clear(); + trackster_barycenter_phi.clear(); + trackster_id_probabilities.clear(); + trackster_vertices_indexes.clear(); + trackster_vertices_x.clear(); + trackster_vertices_y.clear(); + trackster_vertices_z.clear(); + trackster_vertices_time.clear(); + trackster_vertices_timeErr.clear(); + trackster_vertices_energy.clear(); + trackster_vertices_correctedEnergy.clear(); + trackster_vertices_correctedEnergyUncertainty.clear(); + trackster_vertices_multiplicity.clear(); + } + + void fillFromEvent(std::vector const& tracksters, + std::vector const& clusters, + edm::ValueMap> const& layerClustersTimes, + DetectorTools const& detectorTools, + edm::Handle> simClusters_h, + edm::Handle> caloparticles_h, + std::vector const& tracks) { + nTracksters = tracksters.size(); + nClusters = clusters.size(); + for (auto trackster_iterator = tracksters.begin(); trackster_iterator != tracksters.end(); ++trackster_iterator) { + //per-trackster analysis + trackster_time.push_back(trackster_iterator->time()); + trackster_timeError.push_back(trackster_iterator->timeError()); + trackster_regressed_energy.push_back(trackster_iterator->regressed_energy()); + trackster_raw_energy.push_back(trackster_iterator->raw_energy()); + trackster_raw_em_energy.push_back(trackster_iterator->raw_em_energy()); + trackster_raw_pt.push_back(trackster_iterator->raw_pt()); + trackster_raw_em_pt.push_back(trackster_iterator->raw_em_pt()); + trackster_barycenter_x.push_back(trackster_iterator->barycenter().x()); + trackster_barycenter_y.push_back(trackster_iterator->barycenter().y()); + trackster_barycenter_z.push_back(trackster_iterator->barycenter().z()); + trackster_barycenter_eta.push_back(trackster_iterator->barycenter().eta()); + trackster_barycenter_phi.push_back(trackster_iterator->barycenter().phi()); + trackster_EV1.push_back(trackster_iterator->eigenvalues()[0]); + trackster_EV2.push_back(trackster_iterator->eigenvalues()[1]); + trackster_EV3.push_back(trackster_iterator->eigenvalues()[2]); + trackster_eVector0_x.push_back((trackster_iterator->eigenvectors()[0]).x()); + trackster_eVector0_y.push_back((trackster_iterator->eigenvectors()[0]).y()); + trackster_eVector0_z.push_back((trackster_iterator->eigenvectors()[0]).z()); + trackster_sigmaPCA1.push_back(trackster_iterator->sigmasPCA()[0]); + trackster_sigmaPCA2.push_back(trackster_iterator->sigmasPCA()[1]); + trackster_sigmaPCA3.push_back(trackster_iterator->sigmasPCA()[2]); + + if (tracksterType_ != TracksterType::Trackster) { // is simtrackster + auto const& simclusters = *simClusters_h; + auto const& caloparticles = *caloparticles_h; + + simtrackster_timeBoundary.push_back(trackster_iterator->boundaryTime()); + + if (tracksterType_ == TracksterType::SimTracksterCP) + simtrackster_pdgID.push_back(caloparticles[trackster_iterator->seedIndex()].pdgId()); + else if (tracksterType_ == TracksterType::SimTracksterSC) + simtrackster_pdgID.push_back(simclusters[trackster_iterator->seedIndex()].pdgId()); + + using CaloObjectVariant = std::variant; + CaloObjectVariant caloObj; + if (trackster_iterator->seedID() == caloparticles_h.id()) { + caloObj = caloparticles[trackster_iterator->seedIndex()]; + } else { + caloObj = simclusters[trackster_iterator->seedIndex()]; + } + + auto const& simTrack = std::visit([](auto&& obj) { return obj.g4Tracks()[0]; }, caloObj); + auto const& caloPt = std::visit([](auto&& obj) { return obj.pt(); }, caloObj); + simtrackster_regressed_pt.push_back(caloPt); + if (simTrack.crossedBoundary()) { + simtrackster_boundaryX.push_back(simTrack.getPositionAtBoundary().x()); + simtrackster_boundaryY.push_back(simTrack.getPositionAtBoundary().y()); + simtrackster_boundaryZ.push_back(simTrack.getPositionAtBoundary().z()); + simtrackster_boundaryEta.push_back(simTrack.getPositionAtBoundary().eta()); + simtrackster_boundaryPhi.push_back(simTrack.getPositionAtBoundary().phi()); + simtrackster_boundaryPx.push_back(simTrack.getMomentumAtBoundary().x()); + simtrackster_boundaryPy.push_back(simTrack.getMomentumAtBoundary().y()); + simtrackster_boundaryPz.push_back(simTrack.getMomentumAtBoundary().z()); + } else { + simtrackster_boundaryX.push_back(-999); + simtrackster_boundaryY.push_back(-999); + simtrackster_boundaryZ.push_back(-999); + simtrackster_boundaryEta.push_back(-999); + simtrackster_boundaryPhi.push_back(-999); + simtrackster_boundaryPx.push_back(-999); + simtrackster_boundaryPy.push_back(-999); + simtrackster_boundaryPz.push_back(-999); + } + + auto const trackIdx = trackster_iterator->trackIdx(); + simtrackster_trackIdx.push_back(trackIdx); + if (trackIdx != -1) { + const auto& track = tracks[trackIdx]; + + int iSide = int(track.eta() > 0); + + const auto& fts = trajectoryStateTransform::outerFreeState((track), &detectorTools.bfield); + // to the HGCal front + const auto& tsos = detectorTools.propagator.propagate(fts, detectorTools.firstDisk_[iSide]->surface()); + if (tsos.isValid()) { + const auto& globalPos = tsos.globalPosition(); + const auto& globalMom = tsos.globalMomentum(); + simtrackster_track_boundaryX.push_back(globalPos.x()); + simtrackster_track_boundaryY.push_back(globalPos.y()); + simtrackster_track_boundaryZ.push_back(globalPos.z()); + simtrackster_track_boundaryEta.push_back(globalPos.eta()); + simtrackster_track_boundaryPhi.push_back(globalPos.phi()); + simtrackster_track_boundaryPx.push_back(globalMom.x()); + simtrackster_track_boundaryPy.push_back(globalMom.y()); + simtrackster_track_boundaryPz.push_back(globalMom.z()); + simtrackster_trackTime.push_back(track.t0()); + } else { + simtrackster_track_boundaryX.push_back(-999); + simtrackster_track_boundaryY.push_back(-999); + simtrackster_track_boundaryZ.push_back(-999); + simtrackster_track_boundaryEta.push_back(-999); + simtrackster_track_boundaryPhi.push_back(-999); + simtrackster_track_boundaryPx.push_back(-999); + simtrackster_track_boundaryPy.push_back(-999); + simtrackster_track_boundaryPz.push_back(-999); + } + } else { + simtrackster_track_boundaryX.push_back(-999); + simtrackster_track_boundaryY.push_back(-999); + simtrackster_track_boundaryZ.push_back(-999); + simtrackster_track_boundaryEta.push_back(-999); + simtrackster_track_boundaryPhi.push_back(-999); + simtrackster_track_boundaryPx.push_back(-999); + simtrackster_track_boundaryPy.push_back(-999); + simtrackster_track_boundaryPz.push_back(-999); + } + } + + std::vector id_probs; + for (size_t i = 0; i < 8; i++) + id_probs.push_back(trackster_iterator->id_probabilities(i)); + trackster_id_probabilities.push_back(id_probs); + + // Clusters + std::vector vertices_indexes; + std::vector vertices_x; + std::vector vertices_y; + std::vector vertices_z; + std::vector vertices_time; + std::vector vertices_timeErr; + std::vector vertices_energy; + std::vector vertices_correctedEnergy; + std::vector vertices_correctedEnergyUncertainty; + for (auto idx : trackster_iterator->vertices()) { + vertices_indexes.push_back(idx); + const auto& associated_cluster = clusters[idx]; + vertices_x.push_back(associated_cluster.x()); + vertices_y.push_back(associated_cluster.y()); + vertices_z.push_back(associated_cluster.z()); + vertices_energy.push_back(associated_cluster.energy()); + vertices_correctedEnergy.push_back(associated_cluster.correctedEnergy()); + vertices_correctedEnergyUncertainty.push_back(associated_cluster.correctedEnergyUncertainty()); + vertices_time.push_back(layerClustersTimes.get(idx).first); + vertices_timeErr.push_back(layerClustersTimes.get(idx).second); + } + trackster_vertices_indexes.push_back(vertices_indexes); + trackster_vertices_x.push_back(vertices_x); + trackster_vertices_y.push_back(vertices_y); + trackster_vertices_z.push_back(vertices_z); + trackster_vertices_time.push_back(vertices_time); + trackster_vertices_timeErr.push_back(vertices_timeErr); + trackster_vertices_energy.push_back(vertices_energy); + trackster_vertices_correctedEnergy.push_back(vertices_correctedEnergy); + trackster_vertices_correctedEnergyUncertainty.push_back(vertices_correctedEnergyUncertainty); + + // Multiplicity + std::vector vertices_multiplicity; + for (auto multiplicity : trackster_iterator->vertex_multiplicity()) { + vertices_multiplicity.push_back(multiplicity); + } + trackster_vertices_multiplicity.push_back(vertices_multiplicity); + } + } + +private: + TracksterType tracksterType_; + + unsigned int nTracksters; + unsigned int nClusters; + std::vector trackster_time; + std::vector trackster_timeError; + std::vector trackster_regressed_energy; + std::vector trackster_raw_energy; + std::vector trackster_raw_em_energy; + std::vector trackster_raw_pt; + std::vector trackster_raw_em_pt; + std::vector trackster_barycenter_x; + std::vector trackster_barycenter_y; + std::vector trackster_barycenter_z; + std::vector trackster_EV1; + std::vector trackster_EV2; + std::vector trackster_EV3; + std::vector trackster_eVector0_x; + std::vector trackster_eVector0_y; + std::vector trackster_eVector0_z; + std::vector trackster_sigmaPCA1; + std::vector trackster_sigmaPCA2; + std::vector trackster_sigmaPCA3; + std::vector trackster_barycenter_eta; + std::vector trackster_barycenter_phi; + + // for simtrackster + std::vector simtrackster_regressed_pt; + std::vector simtrackster_pdgID; + std::vector simtrackster_trackIdx; + std::vector simtrackster_trackTime; + std::vector simtrackster_timeBoundary; + std::vector simtrackster_boundaryX; + std::vector simtrackster_boundaryY; + std::vector simtrackster_boundaryZ; + std::vector simtrackster_boundaryEta; + std::vector simtrackster_boundaryPhi; + std::vector simtrackster_boundaryPx; + std::vector simtrackster_boundaryPy; + std::vector simtrackster_boundaryPz; + std::vector simtrackster_track_boundaryX; + std::vector simtrackster_track_boundaryY; + std::vector simtrackster_track_boundaryZ; + std::vector simtrackster_track_boundaryEta; + std::vector simtrackster_track_boundaryPhi; + std::vector simtrackster_track_boundaryPx; + std::vector simtrackster_track_boundaryPy; + std::vector simtrackster_track_boundaryPz; + + std::vector> trackster_id_probabilities; + std::vector> trackster_vertices_indexes; + std::vector> trackster_vertices_x; + std::vector> trackster_vertices_y; + std::vector> trackster_vertices_z; + std::vector> trackster_vertices_time; + std::vector> trackster_vertices_timeErr; + std::vector> trackster_vertices_energy; + std::vector> trackster_vertices_correctedEnergy; + std::vector> trackster_vertices_correctedEnergyUncertainty; + std::vector> trackster_vertices_multiplicity; +}; + +// Helper class to dump a TracksterToSimTrackster association map (dumps recoToSim and simToReco at the same time) +class TracksterToSimTracksterAssociationHelper { +public: + /** + * To be called once after tree creation. Output branches will be named prefix_recoToSim/simToReco_suffix_score/sharedE/ + * branchPrefix : for example tsCLUE3D. branchSuffix : usually one of SC or CP. + * *Do not copy/move or resize vector holding object after calling this function* + */ + void initTree(TTree* tree, std::string branchPrefix, std::string branchSuffix) { + tree->Branch((branchPrefix + "_recoToSim_" + branchSuffix).c_str(), &recoToSim); + tree->Branch((branchPrefix + "_recoToSim_" + branchSuffix + "_score").c_str(), &recoToSim_score); + tree->Branch((branchPrefix + "_recoToSim_" + branchSuffix + "_sharedE").c_str(), &recoToSim_sharedE); + tree->Branch((branchPrefix + "_simToReco_" + branchSuffix).c_str(), &simToReco); + tree->Branch((branchPrefix + "_simToReco_" + branchSuffix + "_score").c_str(), &simToReco_score); + tree->Branch((branchPrefix + "_simToReco_" + branchSuffix + "_sharedE").c_str(), &simToReco_sharedE); + } + + void clearVariables() { + recoToSim.clear(); + recoToSim_score.clear(); + recoToSim_sharedE.clear(); + simToReco.clear(); + simToReco_score.clear(); + simToReco_sharedE.clear(); + } + + void fillFromEvent(edm::Handle> tracksters_handle, + edm::Handle> simTracksters_h, + ticl::RecoToSimCollectionSimTracksters const& tsRecoSimSCMap, + ticl::SimToRecoCollectionSimTracksters const& tsSimToRecoSCMap) { + auto const& tracksters = *tracksters_handle; + auto const& simTracksters = *simTracksters_h; + + // Reco -> Sim + recoToSim.resize(tracksters.size()); + recoToSim_score.resize(tracksters.size()); + recoToSim_sharedE.resize(tracksters.size()); + for (size_t i = 0; i < tracksters.size(); ++i) { + const edm::Ref tsRef(tracksters_handle, i); + + const auto stsSC_iter = tsRecoSimSCMap.find(tsRef); + if (stsSC_iter != tsRecoSimSCMap.end()) { + const auto& stsSCassociated = stsSC_iter->val; + for (auto& sts : stsSCassociated) { + auto sts_id = (sts.first).get() - (edm::Ref(simTracksters_h, 0)).get(); + recoToSim[i].push_back(sts_id); + recoToSim_score[i].push_back(sts.second.second); + recoToSim_sharedE[i].push_back(sts.second.first); + } + } + } + + // Sim -> Reco + simToReco.resize(simTracksters.size()); + simToReco_score.resize(simTracksters.size()); + simToReco_sharedE.resize(simTracksters.size()); + for (size_t i = 0; i < simTracksters.size(); ++i) { + const edm::Ref stsSCRef(simTracksters_h, i); + + // STS-SC -> CLUE3D + const auto ts_iter = tsSimToRecoSCMap.find(stsSCRef); + if (ts_iter != tsSimToRecoSCMap.end()) { + const auto& tsAssociated = ts_iter->val; + for (auto& ts : tsAssociated) { + auto ts_idx = (ts.first).get() - (edm::Ref(tracksters_handle, 0)).get(); + simToReco[i].push_back(ts_idx); + simToReco_score[i].push_back(ts.second.second); + simToReco_sharedE[i].push_back(ts.second.first); + } + } + } + } + +private: + std::vector> recoToSim; + std::vector> recoToSim_score; + std::vector> recoToSim_sharedE; + std::vector> simToReco; + std::vector> simToReco_score; + std::vector> simToReco_sharedE; +}; + class TICLDumper : public edm::one::EDAnalyzer { public: explicit TICLDumper(const edm::ParameterSet&); @@ -69,21 +572,24 @@ class TICLDumper : public edm::one::EDAnalyzer bfieldH, - const edm::ESHandle propH); - void buildLayers(); - void analyze(const edm::Event&, const edm::EventSetup&) override; void endRun(edm::Run const& iEvent, edm::EventSetup const&) override{}; void endJob() override; // Define Tokens - const edm::EDGetTokenT> tracksters_token_; + const std::vector + tracksters_parameterSets_; ///< A parameter set for each trackster collection to dump (giving tree name, etc) + std::vector>> + tracksters_token_; ///< a token for each trackster collection to dump + std::vector + tracksters_dumperHelpers_; ///< the dumper helpers for each trackster collection to dump + std::vector tracksters_trees; ///< TTree for each trackster collection to dump + const edm::EDGetTokenT> tracksters_in_candidate_token_; const edm::EDGetTokenT> layer_clusters_token_; const edm::EDGetTokenT> ticl_candidates_token_; + const edm::EDGetTokenT> + ticl_candidates_tracksters_token_; ///< trackster collection used by TICLCandidate const edm::EDGetTokenT> tracks_token_; const edm::EDGetTokenT> tracks_mask_token_; const edm::EDGetTokenT> tracks_time_token_; @@ -101,25 +607,31 @@ class TICLDumper : public edm::one::EDAnalyzer> hgcaltracks_px_token_; const edm::EDGetTokenT> hgcaltracks_py_token_; const edm::EDGetTokenT> hgcaltracks_pz_token_; - const edm::EDGetTokenT> tracksters_merged_token_; const edm::EDGetTokenT> muons_token_; const edm::EDGetTokenT>> clustersTime_token_; const edm::EDGetTokenT> tracksterSeeds_token_; + const edm::EDGetTokenT>> superclustering_linkedResultTracksters_token; + const edm::EDGetTokenT recoSuperClusters_token; + const edm::EDGetTokenT recoSuperClusters_caloClusters_token; + const edm::EDGetTokenT> recoSuperClusters_sourceTracksters_token; edm::ESGetToken caloGeometry_token_; - const edm::EDGetTokenT> simTracksters_SC_token_; - const edm::EDGetTokenT> simTracksters_CP_token_; - const edm::EDGetTokenT> simTracksters_PU_token_; + const edm::EDGetTokenT> simTracksters_SC_token_; // needed for simticlcandidate const edm::EDGetTokenT> simTICLCandidate_token_; - const edm::EDGetTokenT tsRecoToSimSC_token_; - const edm::EDGetTokenT tsSimToRecoSC_token_; - const edm::EDGetTokenT tsRecoToSimCP_token_; - const edm::EDGetTokenT tsSimToRecoCP_token_; - const edm::EDGetTokenT MergeRecoToSimSC_token_; - const edm::EDGetTokenT MergeSimToRecoSC_token_; - const edm::EDGetTokenT MergeRecoToSimCP_token_; - const edm::EDGetTokenT MergeSimToRecoCP_token_; - const edm::EDGetTokenT MergeRecoToSimPU_token_; - const edm::EDGetTokenT MergeSimToRecoPU_token_; + + // associators + const std::vector + associations_parameterSets_; ///< A parameter set for each associator collection to dump (with treeName, etc) + std::vector> + associations_simToReco_token_; ///< The tokens for each assocation + std::vector> associations_recoToSim_token_; + std::vector + associations_dumperHelpers_; ///< the dumper helpers for each association map to dump + std::vector>> + associations_tracksterCollection_; ///< the collection of tracksters used by the associator + std::vector>> + associations_simTracksterCollection_; ///< the collection of simtracksters used by the associator + TTree* associations_tree_; + const edm::EDGetTokenT> simclusters_token_; const edm::EDGetTokenT> caloparticles_token_; @@ -128,22 +640,15 @@ class TICLDumper : public edm::one::EDAnalyzer bfield_token_; const edm::ESGetToken propagator_token_; - hgcal::RecHitTools rhtools_; edm::ESGetToken hdc_token_; - const HGCalDDDConstants* hgcons_; - std::unique_ptr firstDisk_[2]; - std::unique_ptr interfaceDisk_[2]; - edm::ESHandle bfield_; - edm::ESHandle propagator_; + std::unique_ptr detectorTools_; bool saveLCs_; - bool saveCLUE3DTracksters_; - bool saveTrackstersMerged_; - bool saveSimTrackstersSC_; - bool saveSimTrackstersCP_; + bool saveSuperclustering_; + bool saveSuperclusteringDNNScore_; + bool saveRecoSuperclusters_; bool saveTICLCandidate_; bool saveSimTICLCandidate_; bool saveTracks_; - bool saveAssociations_; // Output tree TTree* tree_; @@ -151,153 +656,24 @@ class TICLDumper : public edm::one::EDAnalyzer trackster_time; - std::vector trackster_timeError; - std::vector trackster_regressed_energy; - std::vector trackster_raw_energy; - std::vector trackster_raw_em_energy; - std::vector trackster_raw_pt; - std::vector trackster_raw_em_pt; - std::vector trackster_barycenter_x; - std::vector trackster_barycenter_y; - std::vector trackster_barycenter_z; - std::vector trackster_EV1; - std::vector trackster_EV2; - std::vector trackster_EV3; - std::vector trackster_eVector0_x; - std::vector trackster_eVector0_y; - std::vector trackster_eVector0_z; - std::vector trackster_sigmaPCA1; - std::vector trackster_sigmaPCA2; - std::vector trackster_sigmaPCA3; - std::vector trackster_barycenter_eta; - std::vector trackster_barycenter_phi; - std::vector> trackster_id_probabilities; - std::vector> trackster_vertices_indexes; - std::vector> trackster_vertices_x; - std::vector> trackster_vertices_y; - std::vector> trackster_vertices_z; - std::vector> trackster_vertices_time; - std::vector> trackster_vertices_timeErr; - std::vector> trackster_vertices_energy; - std::vector> trackster_vertices_correctedEnergy; - std::vector> trackster_vertices_correctedEnergyUncertainty; - std::vector> trackster_vertices_multiplicity; - - std::vector stsSC_trackster_time; - std::vector stsSC_trackster_timeBoundary; - std::vector stsSC_trackster_timeError; - std::vector stsSC_trackster_regressed_energy; - std::vector stsSC_trackster_regressed_pt; - std::vector stsSC_trackster_raw_energy; - std::vector stsSC_trackster_raw_em_energy; - std::vector stsSC_trackster_raw_pt; - std::vector stsSC_trackster_raw_em_pt; - std::vector stsSC_trackster_barycenter_x; - std::vector stsSC_trackster_barycenter_y; - std::vector stsSC_trackster_barycenter_z; - std::vector stsSC_trackster_barycenter_eta; - std::vector stsSC_trackster_barycenter_phi; - std::vector stsSC_trackster_EV1; - std::vector stsSC_trackster_EV2; - std::vector stsSC_trackster_EV3; - std::vector stsSC_trackster_eVector0_x; - std::vector stsSC_trackster_eVector0_y; - std::vector stsSC_trackster_eVector0_z; - std::vector stsSC_trackster_sigmaPCA1; - std::vector stsSC_trackster_sigmaPCA2; - std::vector stsSC_trackster_sigmaPCA3; - std::vector stsSC_pdgID; - std::vector stsSC_trackIdx; - std::vector stsSC_trackTime; - std::vector stsSC_boundaryX; - std::vector stsSC_boundaryY; - std::vector stsSC_boundaryZ; - std::vector stsSC_boundaryEta; - std::vector stsSC_boundaryPhi; - std::vector stsSC_boundaryPx; - std::vector stsSC_boundaryPy; - std::vector stsSC_boundaryPz; - std::vector stsSC_track_boundaryX; - std::vector stsSC_track_boundaryY; - std::vector stsSC_track_boundaryZ; - std::vector stsSC_track_boundaryEta; - std::vector stsSC_track_boundaryPhi; - std::vector stsSC_track_boundaryPx; - std::vector stsSC_track_boundaryPy; - std::vector stsSC_track_boundaryPz; - std::vector> stsSC_trackster_id_probabilities; - std::vector> stsSC_trackster_vertices_indexes; - std::vector> stsSC_trackster_vertices_x; - std::vector> stsSC_trackster_vertices_y; - std::vector> stsSC_trackster_vertices_z; - std::vector> stsSC_trackster_vertices_time; - std::vector> stsSC_trackster_vertices_timeErr; - std::vector> stsSC_trackster_vertices_energy; - std::vector> stsSC_trackster_vertices_correctedEnergy; - std::vector> stsSC_trackster_vertices_correctedEnergyUncertainty; - std::vector> stsSC_trackster_vertices_multiplicity; - std::vector stsCP_trackster_time; - std::vector stsCP_trackster_timeBoundary; - std::vector stsCP_trackster_timeError; - std::vector stsCP_trackster_regressed_energy; - std::vector stsCP_trackster_regressed_pt; - std::vector stsCP_trackster_raw_energy; - std::vector stsCP_trackster_raw_em_energy; - std::vector stsCP_trackster_raw_pt; - std::vector stsCP_trackster_raw_em_pt; - std::vector stsCP_trackster_barycenter_x; - std::vector stsCP_trackster_barycenter_y; - std::vector stsCP_trackster_barycenter_z; - std::vector stsCP_trackster_barycenter_eta; - std::vector stsCP_trackster_barycenter_phi; - std::vector stsCP_trackster_EV1; - std::vector stsCP_trackster_EV2; - std::vector stsCP_trackster_EV3; - std::vector stsCP_trackster_eVector0_x; - std::vector stsCP_trackster_eVector0_y; - std::vector stsCP_trackster_eVector0_z; - std::vector stsCP_trackster_sigmaPCA1; - std::vector stsCP_trackster_sigmaPCA2; - std::vector stsCP_trackster_sigmaPCA3; - std::vector stsCP_pdgID; - std::vector stsCP_trackIdx; - std::vector stsCP_trackTime; - std::vector stsCP_boundaryX; - std::vector stsCP_boundaryY; - std::vector stsCP_boundaryZ; - std::vector stsCP_boundaryEta; - std::vector stsCP_boundaryPhi; - std::vector stsCP_boundaryPx; - std::vector stsCP_boundaryPy; - std::vector stsCP_boundaryPz; - std::vector stsCP_track_boundaryX; - std::vector stsCP_track_boundaryY; - std::vector stsCP_track_boundaryZ; - std::vector stsCP_track_boundaryEta; - std::vector stsCP_track_boundaryPhi; - std::vector stsCP_track_boundaryPx; - std::vector stsCP_track_boundaryPy; - std::vector stsCP_track_boundaryPz; - std::vector> stsCP_trackster_id_probabilities; - std::vector> stsCP_trackster_vertices_indexes; - std::vector> stsCP_trackster_vertices_x; - std::vector> stsCP_trackster_vertices_y; - std::vector> stsCP_trackster_vertices_z; - std::vector> stsCP_trackster_vertices_time; - std::vector> stsCP_trackster_vertices_timeErr; - std::vector> stsCP_trackster_vertices_energy; - std::vector> stsCP_trackster_vertices_correctedEnergy; - std::vector> stsCP_trackster_vertices_correctedEnergyUncertainty; - std::vector> stsCP_trackster_vertices_multiplicity; + std::vector> + superclustering_linkedResultTracksters; // Map of indices from superclusteredTracksters collection back into ticlTrackstersCLUE3DEM collection + // reco::SuperCluster dump + std::vector recoSuperCluster_rawEnergy; + std::vector recoSuperCluster_energy; + std::vector recoSuperCluster_correctedEnergy; + std::vector recoSuperCluster_position_x; + std::vector recoSuperCluster_position_y; + std::vector recoSuperCluster_position_z; + std::vector recoSuperCluster_position_eta; + std::vector recoSuperCluster_position_phi; + std::vector + recoSuperCluster_seedTs; ///< Index to seed trackster (into the trackster collection used to make superclusters, given by config recoSuperClusters_sourceTracksterCollection) + std::vector> + recoSuperCluster_constituentTs; ///< Indices to all tracksters inside the supercluster (same) std::vector simTICLCandidate_raw_energy; std::vector simTICLCandidate_regressed_energy; @@ -329,77 +705,7 @@ class TICLDumper : public edm::one::EDAnalyzer> tracksters_in_candidate; std::vector track_in_candidate; - // merged tracksters - size_t nTrackstersMerged; - std::vector tracksters_merged_time; - std::vector tracksters_merged_timeError; - std::vector tracksters_merged_regressed_energy; - std::vector tracksters_merged_raw_energy; - std::vector tracksters_merged_raw_em_energy; - std::vector tracksters_merged_raw_pt; - std::vector tracksters_merged_raw_em_pt; - std::vector tracksters_merged_barycenter_x; - std::vector tracksters_merged_barycenter_y; - std::vector tracksters_merged_barycenter_z; - std::vector tracksters_merged_barycenter_eta; - std::vector tracksters_merged_barycenter_phi; - std::vector tracksters_merged_EV1; - std::vector tracksters_merged_EV2; - std::vector tracksters_merged_EV3; - std::vector tracksters_merged_eVector0_x; - std::vector tracksters_merged_eVector0_y; - std::vector tracksters_merged_eVector0_z; - std::vector tracksters_merged_sigmaPCA1; - std::vector tracksters_merged_sigmaPCA2; - std::vector tracksters_merged_sigmaPCA3; - std::vector> tracksters_merged_vertices_indexes; - std::vector> tracksters_merged_vertices_x; - std::vector> tracksters_merged_vertices_y; - std::vector> tracksters_merged_vertices_z; - std::vector> tracksters_merged_vertices_time; - std::vector> tracksters_merged_vertices_timeErr; - std::vector> tracksters_merged_vertices_energy; - std::vector> tracksters_merged_vertices_correctedEnergy; - std::vector> tracksters_merged_vertices_correctedEnergyUncertainty; - std::vector> tracksters_merged_vertices_multiplicity; - std::vector> tracksters_merged_id_probabilities; - - // associations - std::vector> trackstersCLUE3D_recoToSim_SC; - std::vector> trackstersCLUE3D_recoToSim_SC_score; - std::vector> trackstersCLUE3D_recoToSim_SC_sharedE; - std::vector> trackstersCLUE3D_simToReco_SC; - std::vector> trackstersCLUE3D_simToReco_SC_score; - std::vector> trackstersCLUE3D_simToReco_SC_sharedE; - - std::vector> trackstersCLUE3D_recoToSim_CP; - std::vector> trackstersCLUE3D_recoToSim_CP_score; - std::vector> trackstersCLUE3D_recoToSim_CP_sharedE; - std::vector> trackstersCLUE3D_simToReco_CP; - std::vector> trackstersCLUE3D_simToReco_CP_score; - std::vector> trackstersCLUE3D_simToReco_CP_sharedE; - - std::vector> MergeTracksters_recoToSim_SC; - std::vector> MergeTracksters_recoToSim_SC_score; - std::vector> MergeTracksters_recoToSim_SC_sharedE; - std::vector> MergeTracksters_simToReco_SC; - std::vector> MergeTracksters_simToReco_SC_score; - std::vector> MergeTracksters_simToReco_SC_sharedE; - - std::vector> MergeTracksters_recoToSim_CP; - std::vector> MergeTracksters_recoToSim_CP_score; - std::vector> MergeTracksters_recoToSim_CP_sharedE; - std::vector> MergeTracksters_simToReco_CP; - std::vector> MergeTracksters_simToReco_CP_score; - std::vector> MergeTracksters_simToReco_CP_sharedE; - - std::vector> MergeTracksters_recoToSim_PU; - std::vector> MergeTracksters_recoToSim_PU_score; - std::vector> MergeTracksters_recoToSim_PU_sharedE; - std::vector> MergeTracksters_simToReco_PU; - std::vector> MergeTracksters_simToReco_PU_score; - std::vector> MergeTracksters_simToReco_PU_sharedE; - + // Layer clusters std::vector cluster_seedID; std::vector cluster_energy; std::vector cluster_correctedEnergy; @@ -415,6 +721,7 @@ class TICLDumper : public edm::one::EDAnalyzer cluster_timeErr; std::vector cluster_number_of_hits; + // Tracks std::vector track_id; std::vector track_hgcal_x; std::vector track_hgcal_y; @@ -441,156 +748,33 @@ class TICLDumper : public edm::one::EDAnalyzer track_isMuon; std::vector track_isTrackerMuon; - TTree* trackster_tree_; TTree* cluster_tree_; TTree* candidate_tree_; - TTree* tracksters_merged_tree_; - TTree* associations_tree_; - TTree* simtrackstersSC_tree_; - TTree* simtrackstersCP_tree_; + TTree* superclustering_tree_; TTree* tracks_tree_; TTree* simTICLCandidate_tree; }; void TICLDumper::clearVariables() { // event info - ntracksters_ = 0; nclusters_ = 0; - trackster_time.clear(); - trackster_timeError.clear(); - trackster_regressed_energy.clear(); - trackster_raw_energy.clear(); - trackster_raw_em_energy.clear(); - trackster_raw_pt.clear(); - trackster_raw_em_pt.clear(); - trackster_barycenter_x.clear(); - trackster_barycenter_y.clear(); - trackster_barycenter_z.clear(); - trackster_EV1.clear(); - trackster_EV2.clear(); - trackster_EV3.clear(); - trackster_eVector0_x.clear(); - trackster_eVector0_y.clear(); - trackster_eVector0_z.clear(); - trackster_sigmaPCA1.clear(); - trackster_sigmaPCA2.clear(); - trackster_sigmaPCA3.clear(); - trackster_barycenter_eta.clear(); - trackster_barycenter_phi.clear(); - trackster_id_probabilities.clear(); - trackster_vertices_indexes.clear(); - trackster_vertices_x.clear(); - trackster_vertices_y.clear(); - trackster_vertices_z.clear(); - trackster_vertices_time.clear(); - trackster_vertices_timeErr.clear(); - trackster_vertices_energy.clear(); - trackster_vertices_correctedEnergy.clear(); - trackster_vertices_correctedEnergyUncertainty.clear(); - trackster_vertices_multiplicity.clear(); - - stsSC_trackster_time.clear(); - stsSC_trackster_timeBoundary.clear(); - stsSC_trackster_timeError.clear(); - stsSC_trackster_regressed_energy.clear(); - stsSC_trackster_regressed_pt.clear(); - stsSC_trackster_raw_energy.clear(); - stsSC_trackster_raw_em_energy.clear(); - stsSC_trackster_raw_pt.clear(); - stsSC_trackster_raw_em_pt.clear(); - stsSC_trackster_barycenter_x.clear(); - stsSC_trackster_barycenter_y.clear(); - stsSC_trackster_barycenter_z.clear(); - stsSC_trackster_EV1.clear(); - stsSC_trackster_EV2.clear(); - stsSC_trackster_EV3.clear(); - stsSC_trackster_eVector0_x.clear(); - stsSC_trackster_eVector0_y.clear(); - stsSC_trackster_eVector0_z.clear(); - stsSC_trackster_sigmaPCA1.clear(); - stsSC_trackster_sigmaPCA2.clear(); - stsSC_trackster_sigmaPCA3.clear(); - stsSC_trackster_barycenter_eta.clear(); - stsSC_trackster_barycenter_phi.clear(); - stsSC_pdgID.clear(); - stsSC_trackIdx.clear(); - stsSC_trackTime.clear(); - stsSC_boundaryX.clear(); - stsSC_boundaryY.clear(); - stsSC_boundaryZ.clear(); - stsSC_boundaryEta.clear(); - stsSC_boundaryPhi.clear(); - stsSC_boundaryPx.clear(); - stsSC_boundaryPy.clear(); - stsSC_boundaryPz.clear(); - stsSC_track_boundaryX.clear(); - stsSC_track_boundaryY.clear(); - stsSC_track_boundaryZ.clear(); - stsSC_track_boundaryEta.clear(); - stsSC_track_boundaryPhi.clear(); - stsSC_track_boundaryPx.clear(); - stsSC_track_boundaryPy.clear(); - stsSC_track_boundaryPz.clear(); - stsSC_trackster_id_probabilities.clear(); - stsSC_trackster_vertices_indexes.clear(); - stsSC_trackster_vertices_x.clear(); - stsSC_trackster_vertices_y.clear(); - stsSC_trackster_vertices_z.clear(); - stsSC_trackster_vertices_time.clear(); - stsSC_trackster_vertices_timeErr.clear(); - stsSC_trackster_vertices_energy.clear(); - stsSC_trackster_vertices_correctedEnergy.clear(); - stsSC_trackster_vertices_correctedEnergyUncertainty.clear(); - stsSC_trackster_vertices_multiplicity.clear(); - - stsCP_trackster_time.clear(); - stsCP_trackster_timeBoundary.clear(); - stsCP_trackster_timeError.clear(); - stsCP_trackster_regressed_energy.clear(); - stsCP_trackster_regressed_pt.clear(); - stsCP_trackster_raw_energy.clear(); - stsCP_trackster_raw_em_energy.clear(); - stsCP_trackster_raw_pt.clear(); - stsCP_trackster_raw_em_pt.clear(); - stsCP_trackster_barycenter_x.clear(); - stsCP_trackster_barycenter_y.clear(); - stsCP_trackster_barycenter_z.clear(); - stsCP_trackster_sigmaPCA1.clear(); - stsCP_trackster_sigmaPCA2.clear(); - stsCP_trackster_sigmaPCA3.clear(); - stsCP_trackster_barycenter_eta.clear(); - stsCP_trackster_barycenter_phi.clear(); - stsCP_pdgID.clear(); - stsCP_trackIdx.clear(); - stsCP_trackTime.clear(); - stsCP_boundaryX.clear(); - stsCP_boundaryY.clear(); - stsCP_boundaryZ.clear(); - stsCP_boundaryEta.clear(); - stsCP_boundaryPhi.clear(); - stsCP_boundaryPx.clear(); - stsCP_boundaryPy.clear(); - stsCP_boundaryPz.clear(); - stsCP_track_boundaryX.clear(); - stsCP_track_boundaryY.clear(); - stsCP_track_boundaryZ.clear(); - stsCP_track_boundaryEta.clear(); - stsCP_track_boundaryPhi.clear(); - stsCP_track_boundaryPx.clear(); - stsCP_track_boundaryPy.clear(); - stsCP_track_boundaryPz.clear(); - stsCP_trackster_id_probabilities.clear(); - stsCP_trackster_vertices_indexes.clear(); - stsCP_trackster_vertices_x.clear(); - stsCP_trackster_vertices_y.clear(); - stsCP_trackster_vertices_z.clear(); - stsCP_trackster_vertices_time.clear(); - stsCP_trackster_vertices_timeErr.clear(); - stsCP_trackster_vertices_energy.clear(); - stsCP_trackster_vertices_correctedEnergy.clear(); - stsCP_trackster_vertices_correctedEnergyUncertainty.clear(); - stsCP_trackster_vertices_multiplicity.clear(); + for (TracksterDumperHelper& tsDumper : tracksters_dumperHelpers_) { + tsDumper.clearVariables(); + } + + superclustering_linkedResultTracksters.clear(); + + recoSuperCluster_rawEnergy.clear(); + recoSuperCluster_energy.clear(); + recoSuperCluster_correctedEnergy.clear(); + recoSuperCluster_position_x.clear(); + recoSuperCluster_position_y.clear(); + recoSuperCluster_position_z.clear(); + recoSuperCluster_position_eta.clear(); + recoSuperCluster_position_phi.clear(); + recoSuperCluster_seedTs.clear(); + recoSuperCluster_constituentTs.clear(); simTICLCandidate_raw_energy.clear(); simTICLCandidate_regressed_energy.clear(); @@ -621,84 +805,9 @@ void TICLDumper::clearVariables() { tracksters_in_candidate.clear(); track_in_candidate.clear(); - nTrackstersMerged = 0; - tracksters_merged_time.clear(); - tracksters_merged_timeError.clear(); - tracksters_merged_regressed_energy.clear(); - tracksters_merged_raw_energy.clear(); - tracksters_merged_raw_em_energy.clear(); - tracksters_merged_raw_pt.clear(); - tracksters_merged_raw_em_pt.clear(); - tracksters_merged_barycenter_x.clear(); - tracksters_merged_barycenter_y.clear(); - tracksters_merged_barycenter_z.clear(); - tracksters_merged_barycenter_eta.clear(); - tracksters_merged_barycenter_phi.clear(); - tracksters_merged_EV1.clear(); - tracksters_merged_EV2.clear(); - tracksters_merged_EV3.clear(); - tracksters_merged_eVector0_x.clear(); - tracksters_merged_eVector0_y.clear(); - tracksters_merged_eVector0_z.clear(); - tracksters_merged_sigmaPCA1.clear(); - tracksters_merged_sigmaPCA2.clear(); - tracksters_merged_sigmaPCA3.clear(); - tracksters_merged_id_probabilities.clear(); - tracksters_merged_time.clear(); - tracksters_merged_timeError.clear(); - tracksters_merged_regressed_energy.clear(); - tracksters_merged_raw_energy.clear(); - tracksters_merged_raw_em_energy.clear(); - tracksters_merged_raw_pt.clear(); - tracksters_merged_raw_em_pt.clear(); - - tracksters_merged_vertices_indexes.clear(); - tracksters_merged_vertices_x.clear(); - tracksters_merged_vertices_y.clear(); - tracksters_merged_vertices_z.clear(); - tracksters_merged_vertices_time.clear(); - tracksters_merged_vertices_timeErr.clear(); - tracksters_merged_vertices_energy.clear(); - tracksters_merged_vertices_correctedEnergy.clear(); - tracksters_merged_vertices_correctedEnergyUncertainty.clear(); - tracksters_merged_vertices_multiplicity.clear(); - - trackstersCLUE3D_recoToSim_SC.clear(); - trackstersCLUE3D_recoToSim_SC_score.clear(); - trackstersCLUE3D_recoToSim_SC_sharedE.clear(); - trackstersCLUE3D_simToReco_SC.clear(); - trackstersCLUE3D_simToReco_SC_score.clear(); - trackstersCLUE3D_simToReco_SC_sharedE.clear(); - - trackstersCLUE3D_recoToSim_CP.clear(); - trackstersCLUE3D_recoToSim_CP_score.clear(); - trackstersCLUE3D_recoToSim_CP_sharedE.clear(); - trackstersCLUE3D_simToReco_CP.clear(); - trackstersCLUE3D_simToReco_CP_score.clear(); - trackstersCLUE3D_simToReco_CP_sharedE.clear(); - - MergeTracksters_recoToSim_SC.clear(); - MergeTracksters_recoToSim_SC_score.clear(); - MergeTracksters_recoToSim_SC_sharedE.clear(); - MergeTracksters_simToReco_SC.clear(); - MergeTracksters_simToReco_SC_score.clear(); - MergeTracksters_simToReco_SC_sharedE.clear(); - - MergeTracksters_recoToSim_CP.clear(); - MergeTracksters_recoToSim_CP_score.clear(); - MergeTracksters_recoToSim_CP_sharedE.clear(); - MergeTracksters_simToReco_CP.clear(); - MergeTracksters_simToReco_CP_score.clear(); - MergeTracksters_simToReco_CP_sharedE.clear(); - - MergeTracksters_recoToSim_PU.clear(); - MergeTracksters_recoToSim_PU_score.clear(); - MergeTracksters_recoToSim_PU_sharedE.clear(); - MergeTracksters_simToReco_PU.clear(); - MergeTracksters_simToReco_PU_score.clear(); - MergeTracksters_simToReco_PU_sharedE.clear(); - - nsimTrackstersSC = 0; + for (auto& helper : associations_dumperHelpers_) { + helper.clearVariables(); + } cluster_seedID.clear(); cluster_energy.clear(); @@ -743,11 +852,14 @@ void TICLDumper::clearVariables() { }; TICLDumper::TICLDumper(const edm::ParameterSet& ps) - : tracksters_token_(consumes>(ps.getParameter("trackstersclue3d"))), + : tracksters_parameterSets_(ps.getParameter>("tracksterCollections")), + tracksters_token_(), tracksters_in_candidate_token_( consumes>(ps.getParameter("trackstersInCand"))), layer_clusters_token_(consumes>(ps.getParameter("layerClusters"))), ticl_candidates_token_(consumes>(ps.getParameter("ticlcandidates"))), + ticl_candidates_tracksters_token_( + consumes>(ps.getParameter("ticlcandidates"))), tracks_token_(consumes>(ps.getParameter("tracks"))), tracks_time_token_(consumes>(ps.getParameter("tracksTime"))), tracks_time_quality_token_(consumes>(ps.getParameter("tracksTimeQual"))), @@ -756,40 +868,25 @@ TICLDumper::TICLDumper(const edm::ParameterSet& ps) tracks_time_mtd_token_(consumes>(ps.getParameter("tracksTimeMtd"))), tracks_time_mtd_err_token_(consumes>(ps.getParameter("tracksTimeMtdErr"))), tracks_pos_mtd_token_(consumes>(ps.getParameter("tracksPosMtd"))), - tracksters_merged_token_( - consumes>(ps.getParameter("trackstersmerged"))), muons_token_(consumes>(ps.getParameter("muons"))), clustersTime_token_( consumes>>(ps.getParameter("layer_clustersTime"))), + superclustering_linkedResultTracksters_token( + consumes>>(ps.getParameter("superclustering"))), + recoSuperClusters_token( + consumes(ps.getParameter("recoSuperClusters"))), + recoSuperClusters_caloClusters_token( + consumes(ps.getParameter("recoSuperClusters"))), + recoSuperClusters_sourceTracksters_token(consumes>( + ps.getParameter("recoSuperClusters_sourceTracksterCollection"))), caloGeometry_token_(esConsumes()), simTracksters_SC_token_( consumes>(ps.getParameter("simtrackstersSC"))), - simTracksters_CP_token_( - consumes>(ps.getParameter("simtrackstersCP"))), - simTracksters_PU_token_( - consumes>(ps.getParameter("simtrackstersPU"))), simTICLCandidate_token_( consumes>(ps.getParameter("simTICLCandidates"))), - tsRecoToSimSC_token_( - consumes(ps.getParameter("recoToSimAssociatorSC"))), - tsSimToRecoSC_token_( - consumes(ps.getParameter("simToRecoAssociatorSC"))), - tsRecoToSimCP_token_( - consumes(ps.getParameter("recoToSimAssociatorCP"))), - tsSimToRecoCP_token_( - consumes(ps.getParameter("simToRecoAssociatorCP"))), - MergeRecoToSimSC_token_(consumes( - ps.getParameter("MergerecoToSimAssociatorSC"))), - MergeSimToRecoSC_token_(consumes( - ps.getParameter("MergesimToRecoAssociatorSC"))), - MergeRecoToSimCP_token_(consumes( - ps.getParameter("MergerecoToSimAssociatorCP"))), - MergeSimToRecoCP_token_(consumes( - ps.getParameter("MergesimToRecoAssociatorCP"))), - MergeRecoToSimPU_token_(consumes( - ps.getParameter("MergerecoToSimAssociatorPU"))), - MergeSimToRecoPU_token_(consumes( - ps.getParameter("MergesimToRecoAssociatorPU"))), + associations_parameterSets_(ps.getParameter>("associators")), + // The DumperHelpers should not be moved after construction (needed by TTree branch pointers), so construct them all here + associations_dumperHelpers_(associations_parameterSets_.size()), simclusters_token_(consumes(ps.getParameter("simclusters"))), caloparticles_token_(consumes(ps.getParameter("caloparticles"))), geometry_token_(esConsumes()), @@ -799,75 +896,63 @@ TICLDumper::TICLDumper(const edm::ParameterSet& ps) propagator_token_( esConsumes(edm::ESInputTag("", propName_))), saveLCs_(ps.getParameter("saveLCs")), - saveCLUE3DTracksters_(ps.getParameter("saveCLUE3DTracksters")), - saveTrackstersMerged_(ps.getParameter("saveTrackstersMerged")), - saveSimTrackstersSC_(ps.getParameter("saveSimTrackstersSC")), - saveSimTrackstersCP_(ps.getParameter("saveSimTrackstersCP")), + saveSuperclustering_(ps.getParameter("saveSuperclustering")), + //saveSuperclusteringDNNScore_(ps.getParameter("saveSuperclusteringDNNScore")), + saveRecoSuperclusters_(ps.getParameter("saveRecoSuperclusters")), saveTICLCandidate_(ps.getParameter("saveSimTICLCandidate")), saveSimTICLCandidate_(ps.getParameter("saveSimTICLCandidate")), - saveTracks_(ps.getParameter("saveTracks")), - saveAssociations_(ps.getParameter("saveAssociations")) { + saveTracks_(ps.getParameter("saveTracks")) { std::string detectorName_ = (detector_ == "HFNose") ? "HGCalHFNoseSensitive" : "HGCalEESensitive"; hdc_token_ = esConsumes(edm::ESInputTag("", detectorName_)); + + for (edm::ParameterSet const& tracksterPset : tracksters_parameterSets_) { + tracksters_token_.push_back( + consumes>(tracksterPset.getParameter("inputTag"))); + tracksters_dumperHelpers_.emplace_back( + TracksterDumperHelper::tracksterTypeFromString(tracksterPset.getParameter("tracksterType"))); + } + + for (edm::ParameterSet const& associationPset : associations_parameterSets_) { + edm::InputTag const& inputTag = associationPset.getParameter("associatorInputTag"); + associations_recoToSim_token_.push_back(consumes( + edm::InputTag(inputTag.label(), "recoToSim", inputTag.process()))); + associations_simToReco_token_.push_back(consumes( + edm::InputTag(inputTag.label(), "simToReco", inputTag.process()))); + associations_tracksterCollection_.push_back( + consumes>(associationPset.getParameter("tracksterCollection"))); + associations_simTracksterCollection_.push_back( + consumes>(associationPset.getParameter("simTracksterCollection"))); + } }; TICLDumper::~TICLDumper() { clearVariables(); }; void TICLDumper::beginRun(edm::Run const&, edm::EventSetup const& es) { - const CaloGeometry& geom = es.getData(caloGeometry_token_); - rhtools_.setGeometry(geom); - - edm::ESHandle hdc = es.getHandle(hdc_token_); - hgcons_ = hdc.product(); - edm::ESHandle bfield_ = es.getHandle(bfield_token_); - edm::ESHandle propagator = es.getHandle(propagator_token_); - initialize(hgcons_, rhtools_, bfield_, propagator); + detectorTools_ = std::make_unique(es.getData(hdc_token_), + es.getData(caloGeometry_token_), + es.getData(bfield_token_), + es.getData(propagator_token_)); } // Define tree and branches void TICLDumper::beginJob() { edm::Service fs; - if (saveCLUE3DTracksters_) { - trackster_tree_ = fs->make("tracksters", "TICL tracksters"); - trackster_tree_->Branch("event", &ev_event_); - trackster_tree_->Branch("NClusters", &nclusters_); - trackster_tree_->Branch("NTracksters", &ntracksters_); - trackster_tree_->Branch("time", &trackster_time); - trackster_tree_->Branch("timeError", &trackster_timeError); - trackster_tree_->Branch("regressed_energy", &trackster_regressed_energy); - trackster_tree_->Branch("raw_energy", &trackster_raw_energy); - trackster_tree_->Branch("raw_em_energy", &trackster_raw_em_energy); - trackster_tree_->Branch("raw_pt", &trackster_raw_pt); - trackster_tree_->Branch("raw_em_pt", &trackster_raw_em_pt); - trackster_tree_->Branch("barycenter_x", &trackster_barycenter_x); - trackster_tree_->Branch("barycenter_y", &trackster_barycenter_y); - trackster_tree_->Branch("barycenter_z", &trackster_barycenter_z); - trackster_tree_->Branch("barycenter_eta", &trackster_barycenter_eta); - trackster_tree_->Branch("barycenter_phi", &trackster_barycenter_phi); - trackster_tree_->Branch("EV1", &trackster_EV1); - trackster_tree_->Branch("EV2", &trackster_EV2); - trackster_tree_->Branch("EV3", &trackster_EV3); - trackster_tree_->Branch("eVector0_x", &trackster_eVector0_x); - trackster_tree_->Branch("eVector0_y", &trackster_eVector0_y); - trackster_tree_->Branch("eVector0_z", &trackster_eVector0_z); - trackster_tree_->Branch("sigmaPCA1", &trackster_sigmaPCA1); - trackster_tree_->Branch("sigmaPCA2", &trackster_sigmaPCA2); - trackster_tree_->Branch("sigmaPCA3", &trackster_sigmaPCA3); - trackster_tree_->Branch("id_probabilities", &trackster_id_probabilities); - trackster_tree_->Branch("vertices_indexes", &trackster_vertices_indexes); - trackster_tree_->Branch("vertices_x", &trackster_vertices_x); - trackster_tree_->Branch("vertices_y", &trackster_vertices_y); - trackster_tree_->Branch("vertices_z", &trackster_vertices_z); - trackster_tree_->Branch("vertices_time", &trackster_vertices_time); - trackster_tree_->Branch("vertices_timeErr", &trackster_vertices_timeErr); - trackster_tree_->Branch("vertices_energy", &trackster_vertices_energy); - trackster_tree_->Branch("vertices_correctedEnergy", &trackster_vertices_correctedEnergy); - trackster_tree_->Branch("vertices_correctedEnergyUncertainty", &trackster_vertices_correctedEnergyUncertainty); - trackster_tree_->Branch("vertices_multiplicity", &trackster_vertices_multiplicity); + + // Trackster trees + for (unsigned int i = 0; i < tracksters_parameterSets_.size(); i++) { + edm::ParameterSet const& tracksterPset = tracksters_parameterSets_[i]; + TTree* tree = + fs->make(tracksterPset.getParameter("treeName").c_str(), + ("Tracksters : " + tracksterPset.getParameter("treeName") + + " (InputTag : " + tracksterPset.getParameter("inputTag").encode() + ")") + .c_str()); + tracksters_trees.push_back(tree); + tracksters_dumperHelpers_[i].initTree(tree, &eventId_); } if (saveLCs_) { cluster_tree_ = fs->make("clusters", "TICL tracksters"); + cluster_tree_->Branch("event", &eventId_); cluster_tree_->Branch("seedID", &cluster_seedID); cluster_tree_->Branch("energy", &cluster_energy); cluster_tree_->Branch("correctedEnergy", &cluster_correctedEnergy); @@ -885,6 +970,7 @@ void TICLDumper::beginJob() { } if (saveTICLCandidate_) { candidate_tree_ = fs->make("candidates", "TICL candidates"); + candidate_tree_->Branch("event", &eventId_); candidate_tree_->Branch("NCandidates", &nCandidates); candidate_tree_->Branch("candidate_charge", &candidate_charge); candidate_tree_->Branch("candidate_pdgId", &candidate_pdgId); @@ -899,205 +985,39 @@ void TICLDumper::beginJob() { candidate_tree_->Branch("track_in_candidate", &track_in_candidate); candidate_tree_->Branch("tracksters_in_candidate", &tracksters_in_candidate); } - if (saveTrackstersMerged_) { - tracksters_merged_tree_ = fs->make("trackstersMerged", "TICL tracksters merged"); - tracksters_merged_tree_->Branch("event", &ev_event_); - tracksters_merged_tree_->Branch("time", &tracksters_merged_time); - tracksters_merged_tree_->Branch("timeError", &tracksters_merged_timeError); - tracksters_merged_tree_->Branch("regressed_energy", &tracksters_merged_regressed_energy); - tracksters_merged_tree_->Branch("raw_energy", &tracksters_merged_raw_energy); - tracksters_merged_tree_->Branch("raw_em_energy", &tracksters_merged_raw_em_energy); - tracksters_merged_tree_->Branch("raw_pt", &tracksters_merged_raw_pt); - tracksters_merged_tree_->Branch("raw_em_pt", &tracksters_merged_raw_em_pt); - tracksters_merged_tree_->Branch("NTrackstersMerged", &nTrackstersMerged); - tracksters_merged_tree_->Branch("barycenter_x", &tracksters_merged_barycenter_x); - tracksters_merged_tree_->Branch("barycenter_y", &tracksters_merged_barycenter_y); - tracksters_merged_tree_->Branch("barycenter_z", &tracksters_merged_barycenter_z); - tracksters_merged_tree_->Branch("barycenter_eta", &tracksters_merged_barycenter_eta); - tracksters_merged_tree_->Branch("barycenter_phi", &tracksters_merged_barycenter_phi); - tracksters_merged_tree_->Branch("EV1", &tracksters_merged_EV1); - tracksters_merged_tree_->Branch("EV2", &tracksters_merged_EV2); - tracksters_merged_tree_->Branch("EV3", &tracksters_merged_EV3); - tracksters_merged_tree_->Branch("eVector0_x", &tracksters_merged_eVector0_x); - tracksters_merged_tree_->Branch("eVector0_y", &tracksters_merged_eVector0_y); - tracksters_merged_tree_->Branch("eVector0_z", &tracksters_merged_eVector0_z); - tracksters_merged_tree_->Branch("sigmaPCA1", &tracksters_merged_sigmaPCA1); - tracksters_merged_tree_->Branch("sigmaPCA2", &tracksters_merged_sigmaPCA2); - tracksters_merged_tree_->Branch("sigmaPCA3", &tracksters_merged_sigmaPCA3); - tracksters_merged_tree_->Branch("id_probabilities", &tracksters_merged_id_probabilities); - tracksters_merged_tree_->Branch("vertices_indexes", &tracksters_merged_vertices_indexes); - tracksters_merged_tree_->Branch("vertices_x", &tracksters_merged_vertices_x); - tracksters_merged_tree_->Branch("vertices_y", &tracksters_merged_vertices_y); - tracksters_merged_tree_->Branch("vertices_z", &tracksters_merged_vertices_z); - tracksters_merged_tree_->Branch("vertices_time", &tracksters_merged_vertices_time); - tracksters_merged_tree_->Branch("vertices_timeErr", &tracksters_merged_vertices_timeErr); - tracksters_merged_tree_->Branch("vertices_energy", &tracksters_merged_vertices_energy); - tracksters_merged_tree_->Branch("vertices_correctedEnergy", &tracksters_merged_vertices_correctedEnergy); - tracksters_merged_tree_->Branch("vertices_correctedEnergyUncertainty", - &tracksters_merged_vertices_correctedEnergyUncertainty); - tracksters_merged_tree_->Branch("vertices_multiplicity", &tracksters_merged_vertices_multiplicity); + if (saveSuperclustering_ || saveRecoSuperclusters_) { + superclustering_tree_ = fs->make("superclustering", "Superclustering in HGCAL CE-E"); + superclustering_tree_->Branch("event", &eventId_); } - if (saveAssociations_) { - associations_tree_ = fs->make("associations", "Associations"); - associations_tree_->Branch("tsCLUE3D_recoToSim_SC", &trackstersCLUE3D_recoToSim_SC); - associations_tree_->Branch("tsCLUE3D_recoToSim_SC_score", &trackstersCLUE3D_recoToSim_SC_score); - associations_tree_->Branch("tsCLUE3D_recoToSim_SC_sharedE", &trackstersCLUE3D_recoToSim_SC_sharedE); - associations_tree_->Branch("tsCLUE3D_simToReco_SC", &trackstersCLUE3D_simToReco_SC); - associations_tree_->Branch("tsCLUE3D_simToReco_SC_score", &trackstersCLUE3D_simToReco_SC_score); - associations_tree_->Branch("tsCLUE3D_simToReco_SC_sharedE", &trackstersCLUE3D_simToReco_SC_sharedE); - - associations_tree_->Branch("tsCLUE3D_recoToSim_CP", &trackstersCLUE3D_recoToSim_CP); - associations_tree_->Branch("tsCLUE3D_recoToSim_CP_score", &trackstersCLUE3D_recoToSim_CP_score); - associations_tree_->Branch("tsCLUE3D_recoToSim_CP_sharedE", &trackstersCLUE3D_recoToSim_CP_sharedE); - associations_tree_->Branch("tsCLUE3D_simToReco_CP", &trackstersCLUE3D_simToReco_CP); - associations_tree_->Branch("tsCLUE3D_simToReco_CP_score", &trackstersCLUE3D_simToReco_CP_score); - associations_tree_->Branch("tsCLUE3D_simToReco_CP_sharedE", &trackstersCLUE3D_simToReco_CP_sharedE); - - associations_tree_->Branch("Mergetstracksters_recoToSim_SC", &MergeTracksters_recoToSim_SC); - associations_tree_->Branch("Mergetstracksters_recoToSim_SC_score", &MergeTracksters_recoToSim_SC_score); - associations_tree_->Branch("Mergetstracksters_recoToSim_SC_sharedE", &MergeTracksters_recoToSim_SC_sharedE); - associations_tree_->Branch("Mergetstracksters_simToReco_SC", &MergeTracksters_simToReco_SC); - associations_tree_->Branch("Mergetstracksters_simToReco_SC_score", &MergeTracksters_simToReco_SC_score); - associations_tree_->Branch("Mergetstracksters_simToReco_SC_sharedE", &MergeTracksters_simToReco_SC_sharedE); - - associations_tree_->Branch("Mergetracksters_recoToSim_CP", &MergeTracksters_recoToSim_CP); - associations_tree_->Branch("Mergetracksters_recoToSim_CP_score", &MergeTracksters_recoToSim_CP_score); - associations_tree_->Branch("Mergetracksters_recoToSim_CP_sharedE", &MergeTracksters_recoToSim_CP_sharedE); - associations_tree_->Branch("Mergetracksters_simToReco_CP", &MergeTracksters_simToReco_CP); - associations_tree_->Branch("Mergetracksters_simToReco_CP_score", &MergeTracksters_simToReco_CP_score); - associations_tree_->Branch("Mergetracksters_simToReco_CP_sharedE", &MergeTracksters_simToReco_CP_sharedE); - - associations_tree_->Branch("Mergetracksters_recoToSim_PU", &MergeTracksters_recoToSim_PU); - associations_tree_->Branch("Mergetracksters_recoToSim_PU_score", &MergeTracksters_recoToSim_PU_score); - associations_tree_->Branch("Mergetracksters_recoToSim_PU_sharedE", &MergeTracksters_recoToSim_PU_sharedE); - associations_tree_->Branch("Mergetracksters_simToReco_PU", &MergeTracksters_simToReco_PU); - associations_tree_->Branch("Mergetracksters_simToReco_PU_score", &MergeTracksters_simToReco_PU_score); - associations_tree_->Branch("Mergetracksters_simToReco_PU_sharedE", &MergeTracksters_simToReco_PU_sharedE); + if (saveSuperclustering_) { + superclustering_tree_->Branch("linkedResultTracksters", &superclustering_linkedResultTracksters); + } + if (saveRecoSuperclusters_) { + superclustering_tree_->Branch("recoSuperCluster_rawEnergy", &recoSuperCluster_rawEnergy); + superclustering_tree_->Branch("recoSuperCluster_energy", &recoSuperCluster_energy); + superclustering_tree_->Branch("recoSuperCluster_correctedEnergy", &recoSuperCluster_correctedEnergy); + superclustering_tree_->Branch("recoSuperCluster_position_x", &recoSuperCluster_position_x); + superclustering_tree_->Branch("recoSuperCluster_position_y", &recoSuperCluster_position_y); + superclustering_tree_->Branch("recoSuperCluster_position_z", &recoSuperCluster_position_z); + superclustering_tree_->Branch("recoSuperCluster_position_eta", &recoSuperCluster_position_eta); + superclustering_tree_->Branch("recoSuperCluster_position_phi", &recoSuperCluster_position_phi); + superclustering_tree_->Branch("recoSuperCluster_seedTs", &recoSuperCluster_seedTs); + superclustering_tree_->Branch("recoSuperCluster_constituentTs", &recoSuperCluster_constituentTs); } - if (saveSimTrackstersSC_) { - simtrackstersSC_tree_ = fs->make("simtrackstersSC", "TICL simTracksters SC"); - simtrackstersSC_tree_->Branch("event", &ev_event_); - simtrackstersSC_tree_->Branch("NTracksters", &stsSC_ntracksters_); - simtrackstersSC_tree_->Branch("time", &stsSC_trackster_time); - simtrackstersSC_tree_->Branch("timeBoundary", &stsSC_trackster_timeBoundary); - simtrackstersSC_tree_->Branch("timeError", &stsSC_trackster_timeError); - simtrackstersSC_tree_->Branch("regressed_energy", &stsSC_trackster_regressed_energy); - simtrackstersSC_tree_->Branch("regressed_pt", &stsSC_trackster_regressed_pt); - simtrackstersSC_tree_->Branch("raw_energy", &stsSC_trackster_raw_energy); - simtrackstersSC_tree_->Branch("raw_em_energy", &stsSC_trackster_raw_em_energy); - simtrackstersSC_tree_->Branch("raw_pt", &stsSC_trackster_raw_pt); - simtrackstersSC_tree_->Branch("raw_em_pt", &stsSC_trackster_raw_em_pt); - simtrackstersSC_tree_->Branch("barycenter_x", &stsSC_trackster_barycenter_x); - simtrackstersSC_tree_->Branch("barycenter_y", &stsSC_trackster_barycenter_y); - simtrackstersSC_tree_->Branch("barycenter_z", &stsSC_trackster_barycenter_z); - simtrackstersSC_tree_->Branch("barycenter_eta", &stsSC_trackster_barycenter_eta); - simtrackstersSC_tree_->Branch("barycenter_phi", &stsSC_trackster_barycenter_phi); - simtrackstersSC_tree_->Branch("EV1", &stsSC_trackster_EV1); - simtrackstersSC_tree_->Branch("EV2", &stsSC_trackster_EV2); - simtrackstersSC_tree_->Branch("EV3", &stsSC_trackster_EV3); - simtrackstersSC_tree_->Branch("eVector0_x", &stsSC_trackster_eVector0_x); - simtrackstersSC_tree_->Branch("eVector0_y", &stsSC_trackster_eVector0_y); - simtrackstersSC_tree_->Branch("eVector0_z", &stsSC_trackster_eVector0_z); - simtrackstersSC_tree_->Branch("sigmaPCA1", &stsSC_trackster_sigmaPCA1); - simtrackstersSC_tree_->Branch("sigmaPCA2", &stsSC_trackster_sigmaPCA2); - simtrackstersSC_tree_->Branch("sigmaPCA3", &stsSC_trackster_sigmaPCA3); - simtrackstersSC_tree_->Branch("pdgID", &stsSC_pdgID); - simtrackstersSC_tree_->Branch("trackIdx", &stsSC_trackIdx); - simtrackstersSC_tree_->Branch("trackTime", &stsSC_trackTime); - simtrackstersSC_tree_->Branch("boundaryX", &stsSC_boundaryX); - simtrackstersSC_tree_->Branch("boundaryY", &stsSC_boundaryY); - simtrackstersSC_tree_->Branch("boundaryZ", &stsSC_boundaryZ); - simtrackstersSC_tree_->Branch("boundaryEta", &stsSC_boundaryEta); - simtrackstersSC_tree_->Branch("boundaryPhi", &stsSC_boundaryPhi); - simtrackstersSC_tree_->Branch("boundaryPx", &stsSC_boundaryPx); - simtrackstersSC_tree_->Branch("boundaryPy", &stsSC_boundaryPy); - simtrackstersSC_tree_->Branch("boundaryPz", &stsSC_boundaryPz); - simtrackstersSC_tree_->Branch("track_boundaryX", &stsSC_track_boundaryX); - simtrackstersSC_tree_->Branch("track_boundaryY", &stsSC_track_boundaryY); - simtrackstersSC_tree_->Branch("track_boundaryZ", &stsSC_track_boundaryZ); - simtrackstersSC_tree_->Branch("track_boundaryEta", &stsSC_track_boundaryEta); - simtrackstersSC_tree_->Branch("track_boundaryPhi", &stsSC_track_boundaryPhi); - simtrackstersSC_tree_->Branch("track_boundaryPx", &stsSC_track_boundaryPx); - simtrackstersSC_tree_->Branch("track_boundaryPy", &stsSC_track_boundaryPy); - simtrackstersSC_tree_->Branch("track_boundaryPz", &stsSC_track_boundaryPz); - simtrackstersSC_tree_->Branch("id_probabilities", &stsSC_trackster_id_probabilities); - simtrackstersSC_tree_->Branch("vertices_indexes", &stsSC_trackster_vertices_indexes); - simtrackstersSC_tree_->Branch("vertices_x", &stsSC_trackster_vertices_x); - simtrackstersSC_tree_->Branch("vertices_y", &stsSC_trackster_vertices_y); - simtrackstersSC_tree_->Branch("vertices_z", &stsSC_trackster_vertices_z); - simtrackstersSC_tree_->Branch("vertices_time", &stsSC_trackster_vertices_time); - simtrackstersSC_tree_->Branch("vertices_timeErr", &stsSC_trackster_vertices_timeErr); - simtrackstersSC_tree_->Branch("vertices_energy", &stsSC_trackster_vertices_energy); - simtrackstersSC_tree_->Branch("vertices_correctedEnergy", &stsSC_trackster_vertices_correctedEnergy); - simtrackstersSC_tree_->Branch("vertices_correctedEnergyUncertainty", - &stsSC_trackster_vertices_correctedEnergyUncertainty); - simtrackstersSC_tree_->Branch("vertices_multiplicity", &stsSC_trackster_vertices_multiplicity); - simtrackstersSC_tree_->Branch("NsimTrackstersSC", &nsimTrackstersSC); + if (associations_parameterSets_.size() > 0) { + associations_tree_ = fs->make("associations", "Associations"); + associations_tree_->Branch("event", &eventId_); } - if (saveSimTrackstersCP_) { - simtrackstersCP_tree_ = fs->make("simtrackstersCP", "TICL simTracksters CP"); - simtrackstersCP_tree_->Branch("event", &ev_event_); - simtrackstersCP_tree_->Branch("NTracksters", &stsCP_ntracksters_); - simtrackstersCP_tree_->Branch("time", &stsCP_trackster_time); - simtrackstersCP_tree_->Branch("timeBoundary", &stsCP_trackster_timeBoundary); - simtrackstersCP_tree_->Branch("timeError", &stsCP_trackster_timeError); - simtrackstersCP_tree_->Branch("regressed_energy", &stsCP_trackster_regressed_energy); - simtrackstersCP_tree_->Branch("regressed_pt", &stsCP_trackster_regressed_pt); - simtrackstersCP_tree_->Branch("raw_energy", &stsCP_trackster_raw_energy); - simtrackstersCP_tree_->Branch("raw_em_energy", &stsCP_trackster_raw_em_energy); - simtrackstersCP_tree_->Branch("raw_pt", &stsCP_trackster_raw_pt); - simtrackstersCP_tree_->Branch("raw_em_pt", &stsCP_trackster_raw_em_pt); - simtrackstersCP_tree_->Branch("barycenter_x", &stsCP_trackster_barycenter_x); - simtrackstersCP_tree_->Branch("barycenter_y", &stsCP_trackster_barycenter_y); - simtrackstersCP_tree_->Branch("barycenter_z", &stsCP_trackster_barycenter_z); - simtrackstersCP_tree_->Branch("barycenter_eta", &stsCP_trackster_barycenter_eta); - simtrackstersCP_tree_->Branch("barycenter_phi", &stsCP_trackster_barycenter_phi); - simtrackstersCP_tree_->Branch("pdgID", &stsCP_pdgID); - simtrackstersCP_tree_->Branch("trackIdx", &stsCP_trackIdx); - simtrackstersCP_tree_->Branch("trackTime", &stsCP_trackTime); - simtrackstersCP_tree_->Branch("boundaryX", &stsCP_boundaryX); - simtrackstersCP_tree_->Branch("boundaryY", &stsCP_boundaryY); - simtrackstersCP_tree_->Branch("boundaryZ", &stsCP_boundaryZ); - simtrackstersCP_tree_->Branch("boundaryEta", &stsCP_boundaryEta); - simtrackstersCP_tree_->Branch("boundaryPhi", &stsCP_boundaryPhi); - simtrackstersCP_tree_->Branch("boundaryPx", &stsCP_boundaryPx); - simtrackstersCP_tree_->Branch("boundaryPy", &stsCP_boundaryPy); - simtrackstersCP_tree_->Branch("boundaryPz", &stsCP_boundaryPz); - simtrackstersCP_tree_->Branch("track_boundaryX", &stsCP_track_boundaryX); - simtrackstersCP_tree_->Branch("track_boundaryY", &stsCP_track_boundaryY); - simtrackstersCP_tree_->Branch("track_boundaryZ", &stsCP_track_boundaryZ); - simtrackstersCP_tree_->Branch("track_boundaryEta", &stsCP_track_boundaryEta); - simtrackstersCP_tree_->Branch("track_boundaryPhi", &stsCP_track_boundaryPhi); - simtrackstersCP_tree_->Branch("track_boundaryPx", &stsCP_track_boundaryPx); - simtrackstersCP_tree_->Branch("track_boundaryPy", &stsCP_track_boundaryPy); - simtrackstersCP_tree_->Branch("track_boundaryPz", &stsCP_track_boundaryPz); - simtrackstersCP_tree_->Branch("EV1", &stsCP_trackster_EV1); - simtrackstersCP_tree_->Branch("EV2", &stsCP_trackster_EV2); - simtrackstersCP_tree_->Branch("EV3", &stsCP_trackster_EV3); - simtrackstersCP_tree_->Branch("eVector0_x", &stsCP_trackster_eVector0_x); - simtrackstersCP_tree_->Branch("eVector0_y", &stsCP_trackster_eVector0_y); - simtrackstersCP_tree_->Branch("eVector0_z", &stsCP_trackster_eVector0_z); - simtrackstersCP_tree_->Branch("sigmaPCA1", &stsCP_trackster_sigmaPCA1); - simtrackstersCP_tree_->Branch("sigmaPCA2", &stsCP_trackster_sigmaPCA2); - simtrackstersCP_tree_->Branch("sigmaPCA3", &stsCP_trackster_sigmaPCA3); - simtrackstersCP_tree_->Branch("id_probabilities", &stsCP_trackster_id_probabilities); - simtrackstersCP_tree_->Branch("vertices_indexes", &stsCP_trackster_vertices_indexes); - simtrackstersCP_tree_->Branch("vertices_x", &stsCP_trackster_vertices_x); - simtrackstersCP_tree_->Branch("vertices_y", &stsCP_trackster_vertices_y); - simtrackstersCP_tree_->Branch("vertices_z", &stsCP_trackster_vertices_z); - simtrackstersCP_tree_->Branch("vertices_time", &stsCP_trackster_vertices_time); - simtrackstersCP_tree_->Branch("vertices_timeErr", &stsCP_trackster_vertices_timeErr); - simtrackstersCP_tree_->Branch("vertices_energy", &stsCP_trackster_vertices_energy); - simtrackstersCP_tree_->Branch("vertices_correctedEnergy", &stsCP_trackster_vertices_correctedEnergy); - simtrackstersCP_tree_->Branch("vertices_correctedEnergyUncertainty", - &stsCP_trackster_vertices_correctedEnergyUncertainty); - simtrackstersCP_tree_->Branch("vertices_multiplicity", &stsCP_trackster_vertices_multiplicity); + for (unsigned int i = 0; i < associations_parameterSets_.size(); i++) { + associations_dumperHelpers_[i].initTree(associations_tree_, + associations_parameterSets_[i].getParameter("branchName"), + associations_parameterSets_[i].getParameter("suffix")); } if (saveTracks_) { tracks_tree_ = fs->make("tracks", "Tracks"); - tracks_tree_->Branch("event", &ev_event_); + tracks_tree_->Branch("event", &eventId_); tracks_tree_->Branch("track_id", &track_id); tracks_tree_->Branch("track_hgcal_x", &track_hgcal_x); tracks_tree_->Branch("track_hgcal_y", &track_hgcal_y); @@ -1124,6 +1044,7 @@ void TICLDumper::beginJob() { if (saveSimTICLCandidate_) { simTICLCandidate_tree = fs->make("simTICLCandidate", "Sim TICL Candidate"); + simTICLCandidate_tree->Branch("event", &eventId_); simTICLCandidate_tree->Branch("simTICLCandidate_raw_energy", &simTICLCandidate_raw_energy); simTICLCandidate_tree->Branch("simTICLCandidate_regressed_energy", &simTICLCandidate_regressed_energy); simTICLCandidate_tree->Branch("simTICLCandidate_simTracksterCPIndex", &simTICLCandidate_simTracksterCPIndex); @@ -1141,52 +1062,9 @@ void TICLDumper::beginJob() { } } -void TICLDumper::buildLayers() { - // build disks at HGCal front & EM-Had interface for track propagation - - float zVal = hgcons_->waferZ(1, true); - std::pair rMinMax = hgcons_->rangeR(zVal, true); - - float zVal_interface = rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z(); - std::pair rMinMax_interface = hgcons_->rangeR(zVal_interface, true); - - for (int iSide = 0; iSide < 2; ++iSide) { - float zSide = (iSide == 0) ? (-1. * zVal) : zVal; - firstDisk_[iSide] = - std::make_unique(Disk::build(Disk::PositionType(0, 0, zSide), - Disk::RotationType(), - SimpleDiskBounds(rMinMax.first, rMinMax.second, zSide - 0.5, zSide + 0.5)) - .get()); - - zSide = (iSide == 0) ? (-1. * zVal_interface) : zVal_interface; - interfaceDisk_[iSide] = std::make_unique( - Disk::build(Disk::PositionType(0, 0, zSide), - Disk::RotationType(), - SimpleDiskBounds(rMinMax_interface.first, rMinMax_interface.second, zSide - 0.5, zSide + 0.5)) - .get()); - } -} - -void TICLDumper::initialize(const HGCalDDDConstants* hgcons, - const hgcal::RecHitTools rhtools, - const edm::ESHandle bfieldH, - const edm::ESHandle propH) { - hgcons_ = hgcons; - rhtools_ = rhtools; - buildLayers(); - - bfield_ = bfieldH; - propagator_ = propH; -} void TICLDumper::analyze(const edm::Event& event, const edm::EventSetup& setup) { - ev_event_ += 1; + eventId_ = event.id(); clearVariables(); - auto bFieldProd = bfield_.product(); - const Propagator& prop = (*propagator_); - //get all the tracksters - edm::Handle> tracksters_handle; - event.getByToken(tracksters_token_, tracksters_handle); - const auto& tracksters = *tracksters_handle; edm::Handle> tracksters_in_candidate_handle; event.getByToken(tracksters_in_candidate_token_, tracksters_in_candidate_handle); @@ -1204,6 +1082,8 @@ void TICLDumper::analyze(const edm::Event& event, const edm::EventSetup& setup) edm::Handle> candidates_h; event.getByToken(ticl_candidates_token_, candidates_h); const auto& ticlcandidates = *candidates_h; + edm::Handle> ticlcandidates_tracksters_h = + event.getHandle(ticl_candidates_tracksters_token_); //Track edm::Handle> tracks_h; @@ -1238,452 +1118,77 @@ void TICLDumper::analyze(const edm::Event& event, const edm::EventSetup& setup) event.getByToken(tracks_pos_mtd_token_, trackPosMtd_h); const auto& trackPosMtd = *trackPosMtd_h; - //Tracksters merged - edm::Handle> tracksters_merged_h; - event.getByToken(tracksters_merged_token_, tracksters_merged_h); - const auto& trackstersmerged = *tracksters_merged_h; + // superclustering + if (saveSuperclustering_) // To support running with Mustache + superclustering_linkedResultTracksters = event.get(superclustering_linkedResultTracksters_token); // muons edm::Handle> muons_h; event.getByToken(muons_token_, muons_h); auto& muons = *muons_h; - // simTracksters from SC - edm::Handle> simTrackstersSC_h; - event.getByToken(simTracksters_SC_token_, simTrackstersSC_h); - const auto& simTrackstersSC = *simTrackstersSC_h; - - // simTracksters from CP - edm::Handle> simTrackstersCP_h; - event.getByToken(simTracksters_CP_token_, simTrackstersCP_h); - const auto& simTrackstersCP = *simTrackstersCP_h; + // recoSuperClusters + if (saveRecoSuperclusters_) { + reco::SuperClusterCollection const& recoSuperClusters = event.get(recoSuperClusters_token); + // reco::CaloClusterCollection const& recoCaloClusters = event.get(recoSuperClusters_caloClusters_token); + std::vector const& recoSuperClusters_sourceTracksters = + event.get(recoSuperClusters_sourceTracksters_token); + + // Map for fast lookup of hit to trackster index in recoSuperClusters_sourceTracksters + std::unordered_map hitToTracksterMap; + + for (unsigned ts_id = 0; ts_id < recoSuperClusters_sourceTracksters.size(); ts_id++) { + for (unsigned int lc_index : recoSuperClusters_sourceTracksters[ts_id].vertices()) { + for (auto [detId, fraction] : clusters[lc_index].hitsAndFractions()) { + bool insertionSucceeded = hitToTracksterMap.emplace(detId, ts_id).second; + assert(insertionSucceeded && "TICLDumper found tracksters sharing rechits"); + } + } + } - // simTracksters from PU - edm::Handle> simTrackstersPU_h; - event.getByToken(simTracksters_PU_token_, simTrackstersPU_h); - const auto& simTrackstersPU = *simTrackstersPU_h; + for (auto const& recoSc : recoSuperClusters) { + recoSuperCluster_rawEnergy.push_back(recoSc.rawEnergy()); + recoSuperCluster_energy.push_back(recoSc.energy()); + recoSuperCluster_correctedEnergy.push_back(recoSc.correctedEnergy()); + recoSuperCluster_position_x.push_back(recoSc.position().x()); + recoSuperCluster_position_y.push_back(recoSc.position().y()); + recoSuperCluster_position_z.push_back(recoSc.position().z()); + recoSuperCluster_position_eta.push_back(recoSc.position().eta()); + recoSuperCluster_position_phi.push_back(recoSc.position().phi()); + + // Finding the trackster that was used to create the CaloCluster, using the DetId of a hit (we assume there is no sharing of rechits between tracksters) + + // Seed trackster of the supercluster : Using the DetId of the seed rechit of the seed CaloCluster + recoSuperCluster_seedTs.push_back(hitToTracksterMap.at(recoSc.seed()->seed())); + recoSuperCluster_constituentTs.emplace_back(); + for (edm::Ptr const& caloClusterPtr : recoSc.clusters()) { + // Using the DetId of the seed rechit of the CaloCluster + recoSuperCluster_constituentTs.back().push_back(hitToTracksterMap.at(caloClusterPtr->seed())); + } + } + } edm::Handle> simTICLCandidates_h; event.getByToken(simTICLCandidate_token_, simTICLCandidates_h); const auto& simTICLCandidates = *simTICLCandidates_h; - // trackster reco to sim SC - edm::Handle tsRecoToSimSC_h; - event.getByToken(tsRecoToSimSC_token_, tsRecoToSimSC_h); - auto const& tsRecoSimSCMap = *tsRecoToSimSC_h; - - // sim simTrackster SC to reco trackster - edm::Handle tsSimToRecoSC_h; - event.getByToken(tsSimToRecoSC_token_, tsSimToRecoSC_h); - auto const& tsSimToRecoSCMap = *tsSimToRecoSC_h; - - // trackster reco to sim CP - edm::Handle tsRecoToSimCP_h; - event.getByToken(tsRecoToSimCP_token_, tsRecoToSimCP_h); - auto const& tsRecoSimCPMap = *tsRecoToSimCP_h; - - // sim simTrackster CP to reco trackster - edm::Handle tsSimToRecoCP_h; - event.getByToken(tsSimToRecoCP_token_, tsSimToRecoCP_h); - auto const& tsSimToRecoCPMap = *tsSimToRecoCP_h; - - edm::Handle mergetsRecoToSimSC_h; - event.getByToken(MergeRecoToSimSC_token_, mergetsRecoToSimSC_h); - auto const& MergetsRecoSimSCMap = *mergetsRecoToSimSC_h; - - // sim simTrackster SC to reco trackster - edm::Handle mergetsSimToRecoSC_h; - event.getByToken(MergeSimToRecoSC_token_, mergetsSimToRecoSC_h); - auto const& MergetsSimToRecoSCMap = *mergetsSimToRecoSC_h; - - // trackster reco to sim CP - edm::Handle mergetsRecoToSimCP_h; - event.getByToken(MergeRecoToSimCP_token_, mergetsRecoToSimCP_h); - auto const& MergetsRecoSimCPMap = *mergetsRecoToSimCP_h; - - // sim simTrackster CP to reco trackster - edm::Handle mergetsSimToRecoCP_h; - event.getByToken(MergeSimToRecoCP_token_, mergetsSimToRecoCP_h); - auto const& MergetsSimToRecoCPMap = *mergetsSimToRecoCP_h; - - // trackster reco to sim PU - edm::Handle mergetsRecoToSimPU_h; - event.getByToken(MergeRecoToSimPU_token_, mergetsRecoToSimPU_h); - auto const& MergetsRecoSimPUMap = *mergetsRecoToSimPU_h; - - // sim simTrackster PU to reco trackster - edm::Handle mergetsSimToRecoPU_h; - event.getByToken(MergeSimToRecoPU_token_, mergetsSimToRecoPU_h); - auto const& MergetsSimToRecoPUMap = *mergetsSimToRecoPU_h; - edm::Handle> caloparticles_h; event.getByToken(caloparticles_token_, caloparticles_h); - const auto& caloparticles = *caloparticles_h; - const auto& simclusters = event.get(simclusters_token_); + auto simclusters_h = event.getHandle(simclusters_token_); - ntracksters_ = tracksters.size(); nclusters_ = clusters.size(); - for (auto trackster_iterator = tracksters.begin(); trackster_iterator != tracksters.end(); ++trackster_iterator) { - //per-trackster analysis - trackster_time.push_back(trackster_iterator->time()); - trackster_timeError.push_back(trackster_iterator->timeError()); - trackster_regressed_energy.push_back(trackster_iterator->regressed_energy()); - trackster_raw_energy.push_back(trackster_iterator->raw_energy()); - trackster_raw_em_energy.push_back(trackster_iterator->raw_em_energy()); - trackster_raw_pt.push_back(trackster_iterator->raw_pt()); - trackster_raw_em_pt.push_back(trackster_iterator->raw_em_pt()); - trackster_barycenter_x.push_back(trackster_iterator->barycenter().x()); - trackster_barycenter_y.push_back(trackster_iterator->barycenter().y()); - trackster_barycenter_z.push_back(trackster_iterator->barycenter().z()); - trackster_barycenter_eta.push_back(trackster_iterator->barycenter().eta()); - trackster_barycenter_phi.push_back(trackster_iterator->barycenter().phi()); - trackster_EV1.push_back(trackster_iterator->eigenvalues()[0]); - trackster_EV2.push_back(trackster_iterator->eigenvalues()[1]); - trackster_EV3.push_back(trackster_iterator->eigenvalues()[2]); - trackster_eVector0_x.push_back((trackster_iterator->eigenvectors()[0]).x()); - trackster_eVector0_y.push_back((trackster_iterator->eigenvectors()[0]).y()); - trackster_eVector0_z.push_back((trackster_iterator->eigenvectors()[0]).z()); - trackster_sigmaPCA1.push_back(trackster_iterator->sigmasPCA()[0]); - trackster_sigmaPCA2.push_back(trackster_iterator->sigmasPCA()[1]); - trackster_sigmaPCA3.push_back(trackster_iterator->sigmasPCA()[2]); - std::vector id_probs; - for (size_t i = 0; i < 8; i++) - id_probs.push_back(trackster_iterator->id_probabilities(i)); - trackster_id_probabilities.push_back(id_probs); - - // Clusters - std::vector vertices_indexes; - std::vector vertices_x; - std::vector vertices_y; - std::vector vertices_z; - std::vector vertices_time; - std::vector vertices_timeErr; - std::vector vertices_energy; - std::vector vertices_correctedEnergy; - std::vector vertices_correctedEnergyUncertainty; - for (auto idx : trackster_iterator->vertices()) { - vertices_indexes.push_back(idx); - auto associated_cluster = (*layer_clusters_h)[idx]; - vertices_x.push_back(associated_cluster.x()); - vertices_y.push_back(associated_cluster.y()); - vertices_z.push_back(associated_cluster.z()); - vertices_energy.push_back(associated_cluster.energy()); - vertices_correctedEnergy.push_back(associated_cluster.correctedEnergy()); - vertices_correctedEnergyUncertainty.push_back(associated_cluster.correctedEnergyUncertainty()); - vertices_time.push_back(layerClustersTimes.get(idx).first); - vertices_timeErr.push_back(layerClustersTimes.get(idx).second); - } - trackster_vertices_indexes.push_back(vertices_indexes); - trackster_vertices_x.push_back(vertices_x); - trackster_vertices_y.push_back(vertices_y); - trackster_vertices_z.push_back(vertices_z); - trackster_vertices_time.push_back(vertices_time); - trackster_vertices_timeErr.push_back(vertices_timeErr); - trackster_vertices_energy.push_back(vertices_energy); - trackster_vertices_correctedEnergy.push_back(vertices_correctedEnergy); - trackster_vertices_correctedEnergyUncertainty.push_back(vertices_correctedEnergyUncertainty); - - // Multiplicity - std::vector vertices_multiplicity; - for (auto multiplicity : trackster_iterator->vertex_multiplicity()) { - vertices_multiplicity.push_back(multiplicity); - } - trackster_vertices_multiplicity.push_back(vertices_multiplicity); - } - - stsSC_ntracksters_ = simTrackstersSC.size(); - using CaloObjectVariant = std::variant; - for (auto trackster_iterator = simTrackstersSC.begin(); trackster_iterator != simTrackstersSC.end(); - ++trackster_iterator) { - //per-trackster analysis - stsSC_trackster_time.push_back(trackster_iterator->time()); - stsSC_trackster_timeBoundary.push_back(trackster_iterator->boundaryTime()); - stsSC_trackster_timeError.push_back(trackster_iterator->timeError()); - stsSC_trackster_regressed_energy.push_back(trackster_iterator->regressed_energy()); - stsSC_trackster_raw_energy.push_back(trackster_iterator->raw_energy()); - stsSC_trackster_raw_em_energy.push_back(trackster_iterator->raw_em_energy()); - stsSC_trackster_raw_pt.push_back(trackster_iterator->raw_pt()); - stsSC_trackster_raw_em_pt.push_back(trackster_iterator->raw_em_pt()); - stsSC_trackster_barycenter_x.push_back(trackster_iterator->barycenter().x()); - stsSC_trackster_barycenter_y.push_back(trackster_iterator->barycenter().y()); - stsSC_trackster_barycenter_z.push_back(trackster_iterator->barycenter().z()); - stsSC_trackster_barycenter_eta.push_back(trackster_iterator->barycenter().eta()); - stsSC_trackster_barycenter_phi.push_back(trackster_iterator->barycenter().phi()); - stsSC_trackster_EV1.push_back(trackster_iterator->eigenvalues()[0]); - stsSC_trackster_EV2.push_back(trackster_iterator->eigenvalues()[1]); - stsSC_trackster_EV3.push_back(trackster_iterator->eigenvalues()[2]); - stsSC_trackster_eVector0_x.push_back((trackster_iterator->eigenvectors()[0]).x()); - stsSC_trackster_eVector0_y.push_back((trackster_iterator->eigenvectors()[0]).y()); - stsSC_trackster_eVector0_z.push_back((trackster_iterator->eigenvectors()[0]).z()); - stsSC_trackster_sigmaPCA1.push_back(trackster_iterator->sigmasPCA()[0]); - stsSC_trackster_sigmaPCA2.push_back(trackster_iterator->sigmasPCA()[1]); - stsSC_trackster_sigmaPCA3.push_back(trackster_iterator->sigmasPCA()[2]); - stsSC_pdgID.push_back(simclusters[trackster_iterator->seedIndex()].pdgId()); - - CaloObjectVariant caloObj; - if (trackster_iterator->seedID() == caloparticles_h.id()) { - caloObj = caloparticles[trackster_iterator->seedIndex()]; - } else { - caloObj = simclusters[trackster_iterator->seedIndex()]; - } - - auto const& simTrack = std::visit([](auto&& obj) { return obj.g4Tracks()[0]; }, caloObj); - auto const& caloPt = std::visit([](auto&& obj) { return obj.pt(); }, caloObj); - stsSC_trackster_regressed_pt.push_back(caloPt); - if (simTrack.crossedBoundary()) { - stsSC_boundaryX.push_back(simTrack.getPositionAtBoundary().x()); - stsSC_boundaryY.push_back(simTrack.getPositionAtBoundary().y()); - stsSC_boundaryZ.push_back(simTrack.getPositionAtBoundary().z()); - stsSC_boundaryEta.push_back(simTrack.getPositionAtBoundary().eta()); - stsSC_boundaryPhi.push_back(simTrack.getPositionAtBoundary().phi()); - stsSC_boundaryPx.push_back(simTrack.getMomentumAtBoundary().x()); - stsSC_boundaryPy.push_back(simTrack.getMomentumAtBoundary().y()); - stsSC_boundaryPz.push_back(simTrack.getMomentumAtBoundary().z()); - } else { - stsSC_boundaryX.push_back(-999); - stsSC_boundaryY.push_back(-999); - stsSC_boundaryZ.push_back(-999); - stsSC_boundaryEta.push_back(-999); - stsSC_boundaryPhi.push_back(-999); - stsSC_boundaryPx.push_back(-999); - stsSC_boundaryPy.push_back(-999); - stsSC_boundaryPz.push_back(-999); - } - auto const trackIdx = trackster_iterator->trackIdx(); - stsSC_trackIdx.push_back(trackIdx); - if (trackIdx != -1) { - const auto& track = tracks[trackIdx]; - - int iSide = int(track.eta() > 0); - - const auto& fts = trajectoryStateTransform::outerFreeState((track), bFieldProd); - // to the HGCal front - const auto& tsos = prop.propagate(fts, firstDisk_[iSide]->surface()); - if (tsos.isValid()) { - const auto& globalPos = tsos.globalPosition(); - const auto& globalMom = tsos.globalMomentum(); - stsSC_track_boundaryX.push_back(globalPos.x()); - stsSC_track_boundaryY.push_back(globalPos.y()); - stsSC_track_boundaryZ.push_back(globalPos.z()); - stsSC_track_boundaryEta.push_back(globalPos.eta()); - stsSC_track_boundaryPhi.push_back(globalPos.phi()); - stsSC_track_boundaryPx.push_back(globalMom.x()); - stsSC_track_boundaryPy.push_back(globalMom.y()); - stsSC_track_boundaryPz.push_back(globalMom.z()); - stsSC_trackTime.push_back(track.t0()); - } else { - stsSC_track_boundaryX.push_back(-999); - stsSC_track_boundaryY.push_back(-999); - stsSC_track_boundaryZ.push_back(-999); - stsSC_track_boundaryEta.push_back(-999); - stsSC_track_boundaryPhi.push_back(-999); - stsSC_track_boundaryPx.push_back(-999); - stsSC_track_boundaryPy.push_back(-999); - stsSC_track_boundaryPz.push_back(-999); - } - } else { - stsSC_track_boundaryX.push_back(-999); - stsSC_track_boundaryY.push_back(-999); - stsSC_track_boundaryZ.push_back(-999); - stsSC_track_boundaryEta.push_back(-999); - stsSC_track_boundaryPhi.push_back(-999); - stsSC_track_boundaryPx.push_back(-999); - stsSC_track_boundaryPy.push_back(-999); - stsSC_track_boundaryPz.push_back(-999); - } - - std::vector id_probs; - for (size_t i = 0; i < 8; i++) - id_probs.push_back(trackster_iterator->id_probabilities(i)); - stsSC_trackster_id_probabilities.push_back(id_probs); - - // Clusters - std::vector vertices_indexes; - std::vector vertices_x; - std::vector vertices_y; - std::vector vertices_z; - std::vector vertices_time; - std::vector vertices_timeErr; - std::vector vertices_energy; - std::vector vertices_correctedEnergy; - std::vector vertices_correctedEnergyUncertainty; - for (auto idx : trackster_iterator->vertices()) { - vertices_indexes.push_back(idx); - auto associated_cluster = (*layer_clusters_h)[idx]; - vertices_x.push_back(associated_cluster.x()); - vertices_y.push_back(associated_cluster.y()); - vertices_z.push_back(associated_cluster.z()); - vertices_energy.push_back(associated_cluster.energy()); - vertices_correctedEnergy.push_back(associated_cluster.correctedEnergy()); - vertices_correctedEnergyUncertainty.push_back(associated_cluster.correctedEnergyUncertainty()); - vertices_time.push_back(layerClustersTimes.get(idx).first); - vertices_timeErr.push_back(layerClustersTimes.get(idx).second); - } - stsSC_trackster_vertices_indexes.push_back(vertices_indexes); - stsSC_trackster_vertices_x.push_back(vertices_x); - stsSC_trackster_vertices_y.push_back(vertices_y); - stsSC_trackster_vertices_z.push_back(vertices_z); - stsSC_trackster_vertices_time.push_back(vertices_time); - stsSC_trackster_vertices_timeErr.push_back(vertices_timeErr); - stsSC_trackster_vertices_energy.push_back(vertices_energy); - stsSC_trackster_vertices_correctedEnergy.push_back(vertices_correctedEnergy); - stsSC_trackster_vertices_correctedEnergyUncertainty.push_back(vertices_correctedEnergyUncertainty); - - // Multiplicity - std::vector vertices_multiplicity; - for (auto multiplicity : trackster_iterator->vertex_multiplicity()) { - vertices_multiplicity.push_back(multiplicity); - } - stsSC_trackster_vertices_multiplicity.push_back(vertices_multiplicity); - } - - stsCP_ntracksters_ = simTrackstersCP.size(); - - for (auto trackster_iterator = simTrackstersCP.begin(); trackster_iterator != simTrackstersCP.end(); - ++trackster_iterator) { - //per-trackster analysis - stsCP_trackster_time.push_back(trackster_iterator->time()); - stsCP_trackster_timeBoundary.push_back(trackster_iterator->boundaryTime()); - stsCP_trackster_timeError.push_back(trackster_iterator->timeError()); - stsCP_trackster_regressed_energy.push_back(trackster_iterator->regressed_energy()); - stsCP_trackster_raw_energy.push_back(trackster_iterator->raw_energy()); - stsCP_trackster_raw_em_energy.push_back(trackster_iterator->raw_em_energy()); - stsCP_trackster_raw_pt.push_back(trackster_iterator->raw_pt()); - stsCP_trackster_raw_em_pt.push_back(trackster_iterator->raw_em_pt()); - stsCP_trackster_barycenter_x.push_back(trackster_iterator->barycenter().x()); - stsCP_trackster_barycenter_y.push_back(trackster_iterator->barycenter().y()); - stsCP_trackster_barycenter_z.push_back(trackster_iterator->barycenter().z()); - stsCP_trackster_barycenter_eta.push_back(trackster_iterator->barycenter().eta()); - stsCP_trackster_barycenter_phi.push_back(trackster_iterator->barycenter().phi()); - stsCP_trackster_EV1.push_back(trackster_iterator->eigenvalues()[0]); - stsCP_trackster_EV2.push_back(trackster_iterator->eigenvalues()[1]); - stsCP_trackster_EV3.push_back(trackster_iterator->eigenvalues()[2]); - stsCP_trackster_eVector0_x.push_back((trackster_iterator->eigenvectors()[0]).x()); - stsCP_trackster_eVector0_y.push_back((trackster_iterator->eigenvectors()[0]).y()); - stsCP_trackster_eVector0_z.push_back((trackster_iterator->eigenvectors()[0]).z()); - stsCP_trackster_sigmaPCA1.push_back(trackster_iterator->sigmasPCA()[0]); - stsCP_trackster_sigmaPCA2.push_back(trackster_iterator->sigmasPCA()[1]); - stsCP_trackster_sigmaPCA3.push_back(trackster_iterator->sigmasPCA()[2]); - stsCP_pdgID.push_back(caloparticles[trackster_iterator->seedIndex()].pdgId()); - CaloObjectVariant caloObj; - if (trackster_iterator->seedID() == caloparticles_h.id()) { - caloObj = caloparticles[trackster_iterator->seedIndex()]; - } else { - caloObj = simclusters[trackster_iterator->seedIndex()]; - } - - auto const& simTrack = std::visit([](auto&& obj) { return obj.g4Tracks()[0]; }, caloObj); - auto const& caloPt = std::visit([](auto&& obj) { return obj.pt(); }, caloObj); - stsCP_trackster_regressed_pt.push_back(caloPt); - - if (simTrack.crossedBoundary()) { - stsCP_boundaryX.push_back(simTrack.getPositionAtBoundary().x()); - stsCP_boundaryY.push_back(simTrack.getPositionAtBoundary().y()); - stsCP_boundaryZ.push_back(simTrack.getPositionAtBoundary().z()); - stsCP_boundaryEta.push_back(simTrack.getPositionAtBoundary().eta()); - stsCP_boundaryPhi.push_back(simTrack.getPositionAtBoundary().phi()); - stsCP_boundaryPx.push_back(simTrack.getMomentumAtBoundary().x()); - stsCP_boundaryPy.push_back(simTrack.getMomentumAtBoundary().y()); - stsCP_boundaryPz.push_back(simTrack.getMomentumAtBoundary().z()); - } else { - stsCP_boundaryX.push_back(-999); - stsCP_boundaryY.push_back(-999); - stsCP_boundaryZ.push_back(-999); - stsCP_boundaryEta.push_back(-999); - stsCP_boundaryPhi.push_back(-999); - stsCP_boundaryPx.push_back(-999); - stsCP_boundaryPy.push_back(-999); - stsCP_boundaryPz.push_back(-999); - } - auto const trackIdx = trackster_iterator->trackIdx(); - stsCP_trackIdx.push_back(trackIdx); - if (trackIdx != -1) { - const auto& track = tracks[trackIdx]; - - int iSide = int(track.eta() > 0); - - const auto& fts = trajectoryStateTransform::outerFreeState((track), bFieldProd); - // to the HGCal front - const auto& tsos = prop.propagate(fts, firstDisk_[iSide]->surface()); - if (tsos.isValid()) { - const auto& globalPos = tsos.globalPosition(); - const auto& globalMom = tsos.globalMomentum(); - stsCP_track_boundaryX.push_back(globalPos.x()); - stsCP_track_boundaryY.push_back(globalPos.y()); - stsCP_track_boundaryZ.push_back(globalPos.z()); - stsCP_track_boundaryEta.push_back(globalPos.eta()); - stsCP_track_boundaryPhi.push_back(globalPos.phi()); - stsCP_track_boundaryPx.push_back(globalMom.x()); - stsCP_track_boundaryPy.push_back(globalMom.y()); - stsCP_track_boundaryPz.push_back(globalMom.z()); - stsCP_trackTime.push_back(track.t0()); - } else { - stsCP_track_boundaryX.push_back(-999); - stsCP_track_boundaryY.push_back(-999); - stsCP_track_boundaryZ.push_back(-999); - stsCP_track_boundaryEta.push_back(-999); - stsCP_track_boundaryPhi.push_back(-999); - stsCP_track_boundaryPx.push_back(-999); - stsCP_track_boundaryPy.push_back(-999); - stsCP_track_boundaryPz.push_back(-999); - } - } else { - stsCP_track_boundaryX.push_back(-999); - stsCP_track_boundaryY.push_back(-999); - stsCP_track_boundaryZ.push_back(-999); - stsCP_track_boundaryEta.push_back(-999); - stsCP_track_boundaryPhi.push_back(-999); - stsCP_track_boundaryPx.push_back(-999); - stsCP_track_boundaryPy.push_back(-999); - stsCP_track_boundaryPz.push_back(-999); - } - std::vector id_probs; - for (size_t i = 0; i < 8; i++) - id_probs.push_back(trackster_iterator->id_probabilities(i)); - stsCP_trackster_id_probabilities.push_back(id_probs); - - // Clusters - std::vector vertices_indexes; - std::vector vertices_x; - std::vector vertices_y; - std::vector vertices_z; - std::vector vertices_time; - std::vector vertices_timeErr; - std::vector vertices_energy; - std::vector vertices_correctedEnergy; - std::vector vertices_correctedEnergyUncertainty; - for (auto idx : trackster_iterator->vertices()) { - vertices_indexes.push_back(idx); - auto associated_cluster = (*layer_clusters_h)[idx]; - vertices_x.push_back(associated_cluster.x()); - vertices_y.push_back(associated_cluster.y()); - vertices_z.push_back(associated_cluster.z()); - vertices_energy.push_back(associated_cluster.energy()); - vertices_correctedEnergy.push_back(associated_cluster.correctedEnergy()); - vertices_correctedEnergyUncertainty.push_back(associated_cluster.correctedEnergyUncertainty()); - vertices_time.push_back(layerClustersTimes.get(idx).first); - vertices_timeErr.push_back(layerClustersTimes.get(idx).second); - } - stsCP_trackster_vertices_indexes.push_back(vertices_indexes); - stsCP_trackster_vertices_x.push_back(vertices_x); - stsCP_trackster_vertices_y.push_back(vertices_y); - stsCP_trackster_vertices_z.push_back(vertices_z); - stsCP_trackster_vertices_time.push_back(vertices_time); - stsCP_trackster_vertices_timeErr.push_back(vertices_timeErr); - stsCP_trackster_vertices_energy.push_back(vertices_energy); - stsCP_trackster_vertices_correctedEnergy.push_back(vertices_correctedEnergy); - stsCP_trackster_vertices_correctedEnergyUncertainty.push_back(vertices_correctedEnergyUncertainty); - - // Multiplicity - std::vector vertices_multiplicity; - for (auto multiplicity : trackster_iterator->vertex_multiplicity()) { - vertices_multiplicity.push_back(multiplicity); - } - stsCP_trackster_vertices_multiplicity.push_back(vertices_multiplicity); + // Save all the trackster collections + for (unsigned int i = 0; i < tracksters_dumperHelpers_.size(); i++) { + edm::Handle> tracksters_handle; + std::vector const& tracksters = event.get>(tracksters_token_[i]); + tracksters_dumperHelpers_[i].fillFromEvent( + tracksters, clusters, layerClustersTimes, *detectorTools_, simclusters_h, caloparticles_h, tracks); + tracksters_trees[i]->Fill(); } + const auto& simTrackstersSC_h = event.getHandle(simTracksters_SC_token_); simTICLCandidate_track_in_candidate.resize(simTICLCandidates.size(), -1); for (size_t i = 0; i < simTICLCandidates.size(); ++i) { auto const& cand = simTICLCandidates[i]; @@ -1707,9 +1212,9 @@ void TICLDumper::analyze(const edm::Event& event, const edm::EventSetup& setup) int tk_idx = trackPtr.get() - (edm::Ptr(tracks_h, 0)).get(); simTICLCandidate_track_in_candidate[i] = tk_idx; - const auto& fts = trajectoryStateTransform::outerFreeState((track), bFieldProd); + const auto& fts = trajectoryStateTransform::outerFreeState((track), &detectorTools_->bfield); // to the HGCal front - const auto& tsos = prop.propagate(fts, firstDisk_[iSide]->surface()); + const auto& tsos = detectorTools_->propagator.propagate(fts, detectorTools_->firstDisk_[iSide]->surface()); if (tsos.isValid()) { const auto& globalPos = tsos.globalPosition(); const auto& globalMom = tsos.globalMomentum(); @@ -1751,11 +1256,11 @@ void TICLDumper::analyze(const edm::Event& event, const edm::EventSetup& setup) cluster_position_eta.push_back(cluster_iterator->eta()); cluster_position_phi.push_back(cluster_iterator->phi()); auto haf = cluster_iterator->hitsAndFractions(); - auto layerId = rhtools_.getLayerWithOffset(haf[0].first); + auto layerId = detectorTools_->rhtools.getLayerWithOffset(haf[0].first); cluster_layer_id.push_back(layerId); uint32_t number_of_hits = cluster_iterator->hitsAndFractions().size(); cluster_number_of_hits.push_back(number_of_hits); - cluster_type.push_back(rhtools_.getCellType(lc_seed)); + cluster_type.push_back(detectorTools_->rhtools.getCellType(lc_seed)); cluster_timeErr.push_back(layerClustersTimes.get(c_id).second); cluster_time.push_back(layerClustersTimes.get(c_id).first); c_id += 1; @@ -1794,281 +1299,24 @@ void TICLDumper::analyze(const edm::Event& event, const edm::EventSetup& setup) track_in_candidate[i] = tk_idx; } - nTrackstersMerged = trackstersmerged.size(); - for (auto trackster_iterator = trackstersmerged.begin(); trackster_iterator != trackstersmerged.end(); - ++trackster_iterator) { - tracksters_merged_time.push_back(trackster_iterator->time()); - tracksters_merged_timeError.push_back(trackster_iterator->timeError()); - tracksters_merged_regressed_energy.push_back(trackster_iterator->regressed_energy()); - tracksters_merged_raw_energy.push_back(trackster_iterator->raw_energy()); - tracksters_merged_raw_em_energy.push_back(trackster_iterator->raw_em_energy()); - tracksters_merged_raw_pt.push_back(trackster_iterator->raw_pt()); - tracksters_merged_raw_em_pt.push_back(trackster_iterator->raw_em_pt()); - tracksters_merged_barycenter_x.push_back(trackster_iterator->barycenter().x()); - tracksters_merged_barycenter_y.push_back(trackster_iterator->barycenter().y()); - tracksters_merged_barycenter_z.push_back(trackster_iterator->barycenter().z()); - tracksters_merged_barycenter_eta.push_back(trackster_iterator->barycenter().eta()); - tracksters_merged_barycenter_phi.push_back(trackster_iterator->barycenter().phi()); - tracksters_merged_EV1.push_back(trackster_iterator->eigenvalues()[0]); - tracksters_merged_EV2.push_back(trackster_iterator->eigenvalues()[1]); - tracksters_merged_EV3.push_back(trackster_iterator->eigenvalues()[2]); - tracksters_merged_eVector0_x.push_back((trackster_iterator->eigenvectors()[0]).x()); - tracksters_merged_eVector0_y.push_back((trackster_iterator->eigenvectors()[0]).y()); - tracksters_merged_eVector0_z.push_back((trackster_iterator->eigenvectors()[0]).z()); - tracksters_merged_sigmaPCA1.push_back(trackster_iterator->sigmasPCA()[0]); - tracksters_merged_sigmaPCA2.push_back(trackster_iterator->sigmasPCA()[1]); - tracksters_merged_sigmaPCA3.push_back(trackster_iterator->sigmasPCA()[2]); - - std::vector id_probs; - for (size_t i = 0; i < 8; i++) - id_probs.push_back(trackster_iterator->id_probabilities(i)); - tracksters_merged_id_probabilities.push_back(id_probs); - - std::vector vertices_indexes; - std::vector vertices_x; - std::vector vertices_y; - std::vector vertices_z; - std::vector vertices_time; - std::vector vertices_timeErr; - std::vector vertices_energy; - std::vector vertices_correctedEnergy; - std::vector vertices_correctedEnergyUncertainty; - for (auto idx : trackster_iterator->vertices()) { - vertices_indexes.push_back(idx); - auto associated_cluster = (*layer_clusters_h)[idx]; - vertices_x.push_back(associated_cluster.x()); - vertices_y.push_back(associated_cluster.y()); - vertices_z.push_back(associated_cluster.z()); - vertices_energy.push_back(associated_cluster.energy()); - vertices_correctedEnergy.push_back(associated_cluster.correctedEnergy()); - vertices_correctedEnergyUncertainty.push_back(associated_cluster.correctedEnergyUncertainty()); - vertices_time.push_back(layerClustersTimes.get(idx).first); - vertices_timeErr.push_back(layerClustersTimes.get(idx).second); - } - tracksters_merged_vertices_indexes.push_back(vertices_indexes); - tracksters_merged_vertices_x.push_back(vertices_x); - tracksters_merged_vertices_y.push_back(vertices_y); - tracksters_merged_vertices_z.push_back(vertices_z); - tracksters_merged_vertices_time.push_back(vertices_time); - tracksters_merged_vertices_timeErr.push_back(vertices_timeErr); - tracksters_merged_vertices_energy.push_back(vertices_energy); - tracksters_merged_vertices_correctedEnergy.push_back(vertices_correctedEnergy); - tracksters_merged_vertices_correctedEnergyUncertainty.push_back(vertices_correctedEnergyUncertainty); - } - - // Tackster reco->sim associations - trackstersCLUE3D_recoToSim_SC.resize(tracksters.size()); - trackstersCLUE3D_recoToSim_SC_score.resize(tracksters.size()); - trackstersCLUE3D_recoToSim_SC_sharedE.resize(tracksters.size()); - for (size_t i = 0; i < tracksters.size(); ++i) { - const edm::Ref tsRef(tracksters_handle, i); - - // CLUE3D -> STS-SC - const auto stsSC_iter = tsRecoSimSCMap.find(tsRef); - if (stsSC_iter != tsRecoSimSCMap.end()) { - const auto& stsSCassociated = stsSC_iter->val; - for (auto& sts : stsSCassociated) { - auto sts_id = (sts.first).get() - (edm::Ref(simTrackstersSC_h, 0)).get(); - trackstersCLUE3D_recoToSim_SC[i].push_back(sts_id); - trackstersCLUE3D_recoToSim_SC_score[i].push_back(sts.second.second); - trackstersCLUE3D_recoToSim_SC_sharedE[i].push_back(sts.second.first); - } - } - } - - // SimTracksters - nsimTrackstersSC = simTrackstersSC.size(); - trackstersCLUE3D_simToReco_SC.resize(nsimTrackstersSC); - trackstersCLUE3D_simToReco_SC_score.resize(nsimTrackstersSC); - trackstersCLUE3D_simToReco_SC_sharedE.resize(nsimTrackstersSC); - for (size_t i = 0; i < nsimTrackstersSC; ++i) { - const edm::Ref stsSCRef(simTrackstersSC_h, i); - - // STS-SC -> CLUE3D - const auto ts_iter = tsSimToRecoSCMap.find(stsSCRef); - if (ts_iter != tsSimToRecoSCMap.end()) { - const auto& tsAssociated = ts_iter->val; - for (auto& ts : tsAssociated) { - auto ts_idx = (ts.first).get() - (edm::Ref(tracksters_handle, 0)).get(); - trackstersCLUE3D_simToReco_SC[i].push_back(ts_idx); - trackstersCLUE3D_simToReco_SC_score[i].push_back(ts.second.second); - trackstersCLUE3D_simToReco_SC_sharedE[i].push_back(ts.second.first); - } - } - } - - // Tackster reco->sim associations - trackstersCLUE3D_recoToSim_CP.resize(tracksters.size()); - trackstersCLUE3D_recoToSim_CP_score.resize(tracksters.size()); - trackstersCLUE3D_recoToSim_CP_sharedE.resize(tracksters.size()); - for (size_t i = 0; i < tracksters.size(); ++i) { - const edm::Ref tsRef(tracksters_handle, i); - - // CLUE3D -> STS-CP - const auto stsCP_iter = tsRecoSimCPMap.find(tsRef); - if (stsCP_iter != tsRecoSimCPMap.end()) { - const auto& stsCPassociated = stsCP_iter->val; - for (auto& sts : stsCPassociated) { - auto sts_id = (sts.first).get() - (edm::Ref(simTrackstersCP_h, 0)).get(); - trackstersCLUE3D_recoToSim_CP[i].push_back(sts_id); - trackstersCLUE3D_recoToSim_CP_score[i].push_back(sts.second.second); - trackstersCLUE3D_recoToSim_CP_sharedE[i].push_back(sts.second.first); - } - } - } - - // SimTracksters - nsimTrackstersCP = simTrackstersCP.size(); - trackstersCLUE3D_simToReco_CP.resize(nsimTrackstersCP); - trackstersCLUE3D_simToReco_CP_score.resize(nsimTrackstersCP); - trackstersCLUE3D_simToReco_CP_sharedE.resize(nsimTrackstersCP); - for (size_t i = 0; i < nsimTrackstersCP; ++i) { - const edm::Ref stsCPRef(simTrackstersCP_h, i); - - // STS-CP -> CLUE3D - const auto ts_iter = tsSimToRecoCPMap.find(stsCPRef); - if (ts_iter != tsSimToRecoCPMap.end()) { - const auto& tsAssociated = ts_iter->val; - for (auto& ts : tsAssociated) { - auto ts_idx = (ts.first).get() - (edm::Ref(tracksters_handle, 0)).get(); - trackstersCLUE3D_simToReco_CP[i].push_back(ts_idx); - trackstersCLUE3D_simToReco_CP_score[i].push_back(ts.second.second); - trackstersCLUE3D_simToReco_CP_sharedE[i].push_back(ts.second.first); - } - } - } - - // Tackster reco->sim associations - MergeTracksters_recoToSim_SC.resize(trackstersmerged.size()); - MergeTracksters_recoToSim_SC_score.resize(trackstersmerged.size()); - MergeTracksters_recoToSim_SC_sharedE.resize(trackstersmerged.size()); - for (size_t i = 0; i < trackstersmerged.size(); ++i) { - const edm::Ref tsRef(tracksters_merged_h, i); - - // CLUE3D -> STS-SC - const auto stsSC_iter = MergetsRecoSimSCMap.find(tsRef); - if (stsSC_iter != MergetsRecoSimSCMap.end()) { - const auto& stsSCassociated = stsSC_iter->val; - for (auto& sts : stsSCassociated) { - auto sts_id = (sts.first).get() - (edm::Ref(simTrackstersSC_h, 0)).get(); - MergeTracksters_recoToSim_SC[i].push_back(sts_id); - MergeTracksters_recoToSim_SC_score[i].push_back(sts.second.second); - MergeTracksters_recoToSim_SC_sharedE[i].push_back(sts.second.first); - } - } - } - - // SimTracksters - nsimTrackstersSC = simTrackstersSC.size(); - MergeTracksters_simToReco_SC.resize(nsimTrackstersSC); - MergeTracksters_simToReco_SC_score.resize(nsimTrackstersSC); - MergeTracksters_simToReco_SC_sharedE.resize(nsimTrackstersSC); - for (size_t i = 0; i < nsimTrackstersSC; ++i) { - const edm::Ref stsSCRef(simTrackstersSC_h, i); - - // STS-SC -> CLUE3D - const auto ts_iter = MergetsSimToRecoSCMap.find(stsSCRef); - if (ts_iter != MergetsSimToRecoSCMap.end()) { - const auto& tsAssociated = ts_iter->val; - for (auto& ts : tsAssociated) { - auto ts_idx = (ts.first).get() - (edm::Ref(tracksters_merged_h, 0)).get(); - MergeTracksters_simToReco_SC[i].push_back(ts_idx); - MergeTracksters_simToReco_SC_score[i].push_back(ts.second.second); - MergeTracksters_simToReco_SC_sharedE[i].push_back(ts.second.first); - } - } - } - - // Tackster reco->sim associations - MergeTracksters_recoToSim_CP.resize(trackstersmerged.size()); - MergeTracksters_recoToSim_CP_score.resize(trackstersmerged.size()); - MergeTracksters_recoToSim_CP_sharedE.resize(trackstersmerged.size()); - for (size_t i = 0; i < trackstersmerged.size(); ++i) { - const edm::Ref tsRef(tracksters_merged_h, i); - - // CLUE3D -> STS-CP - const auto stsCP_iter = MergetsRecoSimCPMap.find(tsRef); - if (stsCP_iter != MergetsRecoSimCPMap.end()) { - const auto& stsCPassociated = stsCP_iter->val; - for (auto& sts : stsCPassociated) { - auto sts_id = (sts.first).get() - (edm::Ref(simTrackstersCP_h, 0)).get(); - MergeTracksters_recoToSim_CP[i].push_back(sts_id); - MergeTracksters_recoToSim_CP_score[i].push_back(sts.second.second); - MergeTracksters_recoToSim_CP_sharedE[i].push_back(sts.second.first); - } - } - } - - // Tackster reco->sim associations - MergeTracksters_recoToSim_PU.resize(trackstersmerged.size()); - MergeTracksters_recoToSim_PU_score.resize(trackstersmerged.size()); - MergeTracksters_recoToSim_PU_sharedE.resize(trackstersmerged.size()); - for (size_t i = 0; i < trackstersmerged.size(); ++i) { - const edm::Ref tsRef(tracksters_merged_h, i); - - // CLUE3D -> STS-PU - const auto stsPU_iter = MergetsRecoSimPUMap.find(tsRef); - if (stsPU_iter != MergetsRecoSimPUMap.end()) { - const auto& stsPUassociated = stsPU_iter->val; - for (auto& sts : stsPUassociated) { - auto sts_id = (sts.first).get() - (edm::Ref(simTrackstersPU_h, 0)).get(); - MergeTracksters_recoToSim_PU[i].push_back(sts_id); - MergeTracksters_recoToSim_PU_score[i].push_back(sts.second.second); - MergeTracksters_recoToSim_PU_sharedE[i].push_back(sts.second.first); - } - } - } - - // SimTracksters - nsimTrackstersCP = simTrackstersCP.size(); - MergeTracksters_simToReco_CP.resize(nsimTrackstersCP); - MergeTracksters_simToReco_CP_score.resize(nsimTrackstersCP); - MergeTracksters_simToReco_CP_sharedE.resize(nsimTrackstersCP); - for (size_t i = 0; i < nsimTrackstersCP; ++i) { - const edm::Ref stsCPRef(simTrackstersCP_h, i); - - // STS-CP -> TrackstersMerge - const auto ts_iter = MergetsSimToRecoCPMap.find(stsCPRef); - if (ts_iter != MergetsSimToRecoCPMap.end()) { - const auto& tsAssociated = ts_iter->val; - for (auto& ts : tsAssociated) { - auto ts_idx = (ts.first).get() - (edm::Ref(tracksters_merged_h, 0)).get(); - MergeTracksters_simToReco_CP[i].push_back(ts_idx); - MergeTracksters_simToReco_CP_score[i].push_back(ts.second.second); - MergeTracksters_simToReco_CP_sharedE[i].push_back(ts.second.first); - } - } - } - - // SimTracksters - auto nsimTrackstersPU = simTrackstersPU.size(); - MergeTracksters_simToReco_PU.resize(nsimTrackstersPU); - MergeTracksters_simToReco_PU_score.resize(nsimTrackstersPU); - MergeTracksters_simToReco_PU_sharedE.resize(nsimTrackstersPU); - for (size_t i = 0; i < nsimTrackstersPU; ++i) { - const edm::Ref stsPURef(simTrackstersPU_h, i); - - // STS-PU -> Tracksters Merge - const auto ts_iter = MergetsSimToRecoPUMap.find(stsPURef); - if (ts_iter != MergetsSimToRecoPUMap.end()) { - const auto& tsAssociated = ts_iter->val; - for (auto& ts : tsAssociated) { - auto ts_idx = (ts.first).get() - (edm::Ref(tracksters_merged_h, 0)).get(); - MergeTracksters_simToReco_PU[i].push_back(ts_idx); - MergeTracksters_simToReco_PU_score[i].push_back(ts.second.second); - MergeTracksters_simToReco_PU_sharedE[i].push_back(ts.second.first); - } - } + // trackster to simTrackster associations + for (unsigned int i = 0; i < associations_dumperHelpers_.size(); i++) { + associations_dumperHelpers_[i].fillFromEvent(event.getHandle(associations_tracksterCollection_[i]), + event.getHandle(associations_simTracksterCollection_[i]), + event.get(associations_recoToSim_token_[i]), + event.get(associations_simToReco_token_[i])); } + if (associations_dumperHelpers_.size() > 0) + associations_tree_->Fill(); //Tracks for (size_t i = 0; i < tracks.size(); i++) { const auto& track = tracks[i]; reco::TrackRef trackref = reco::TrackRef(tracks_h, i); int iSide = int(track.eta() > 0); - const auto& fts = trajectoryStateTransform::outerFreeState((track), bFieldProd); + const auto& fts = trajectoryStateTransform::outerFreeState((track), &detectorTools_->bfield); // to the HGCal front - const auto& tsos = prop.propagate(fts, firstDisk_[iSide]->surface()); + const auto& tsos = detectorTools_->propagator.propagate(fts, detectorTools_->firstDisk_[iSide]->surface()); if (tsos.isValid()) { const auto& globalPos = tsos.globalPosition(); const auto& globalMom = tsos.globalMomentum(); @@ -2107,20 +1355,12 @@ void TICLDumper::analyze(const edm::Event& event, const edm::EventSetup& setup) } } - if (saveCLUE3DTracksters_) - trackster_tree_->Fill(); if (saveLCs_) cluster_tree_->Fill(); if (saveTICLCandidate_) candidate_tree_->Fill(); - if (saveTrackstersMerged_) - tracksters_merged_tree_->Fill(); - if (saveAssociations_) - associations_tree_->Fill(); - if (saveSimTrackstersSC_) - simtrackstersSC_tree_->Fill(); - if (saveSimTrackstersCP_) - simtrackstersCP_tree_->Fill(); + if (saveSuperclustering_ || saveRecoSuperclusters_) + superclustering_tree_->Fill(); if (saveTracks_) tracks_tree_->Fill(); if (saveSimTICLCandidate_) @@ -2131,8 +1371,24 @@ void TICLDumper::endJob() {} void TICLDumper::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; - desc.add("trackstersclue3d", edm::InputTag("ticlTrackstersCLUE3DHigh")); + + // Settings for dumping trackster collections + edm::ParameterSetDescription tracksterDescValidator; + tracksterDescValidator.add("treeName") + ->setComment("Name of the output tree for the trackster collection"); + tracksterDescValidator.add("inputTag")->setComment("Input tag for the trackster collection to write"); + tracksterDescValidator.ifValue( + edm::ParameterDescription( + "tracksterType", + "Trackster", + true, + edm::Comment("Type of trackster. Trackster=regular trackster (from RECO). SimTracksterCP=Simtrackster " + "from CaloParticle. SimTracksterSC=Simtrackster from SimCluster")), + edm::allowedValues("Trackster", "SimTracksterCP", "SimTracksterSC")); + desc.addVPSet("tracksterCollections", tracksterDescValidator)->setComment("Trackster collections to dump"); + desc.add("trackstersInCand", edm::InputTag("ticlTrackstersCLUE3DHigh")); + desc.add("layerClusters", edm::InputTag("hgcalMergeLayerClusters")); desc.add("layer_clustersTime", edm::InputTag("hgcalMergeLayerClusters", "timeLayerCluster")); desc.add("ticlcandidates", edm::InputTag("ticlTrackstersMerge")); @@ -2144,46 +1400,45 @@ void TICLDumper::fillDescriptions(edm::ConfigurationDescriptions& descriptions) desc.add("tracksTimeMtd", edm::InputTag("trackExtenderWithMTD:generalTracktmtd")); desc.add("tracksTimeMtdErr", edm::InputTag("trackExtenderWithMTD:generalTracksigmatmtd")); desc.add("tracksPosMtd", edm::InputTag("trackExtenderWithMTD:generalTrackmtdpos")); - desc.add("trackstersmerged", edm::InputTag("ticlTrackstersMerge")); desc.add("muons", edm::InputTag("muons1stStep")); - desc.add("simtrackstersSC", edm::InputTag("ticlSimTracksters")); - desc.add("simtrackstersCP", edm::InputTag("ticlSimTracksters", "fromCPs")); - desc.add("simtrackstersPU", edm::InputTag("ticlSimTracksters", "PU")); + desc.add("superclustering", edm::InputTag("ticlTracksterLinksSuperclusteringDNN")); + desc.add("recoSuperClusters", edm::InputTag("particleFlowSuperClusterHGCal")) + ->setComment( + "egamma supercluster collection (either from PFECALSuperClusterProducer for Mustache, or from " + "TICL->Egamma converter in case of TICL DNN superclusters)"); + desc.add("recoSuperClusters_sourceTracksterCollection", edm::InputTag("ticlTrackstersMerge")) + ->setComment( + "Trackster collection used to produce the reco::SuperCluster, used to provide a mapping back to the " + "tracksters used in superclusters"); + + desc.add("simtrackstersSC", edm::InputTag("ticlSimTracksters")) + ->setComment("SimTrackster from CaloParticle collection to use for simTICLcandidates"); desc.add("simTICLCandidates", edm::InputTag("ticlSimTracksters")); - desc.add("recoToSimAssociatorSC", - edm::InputTag("tracksterSimTracksterAssociationPRbyCLUE3D", "recoToSim")); - desc.add("simToRecoAssociatorSC", - edm::InputTag("tracksterSimTracksterAssociationPRbyCLUE3D", "simToReco")); - desc.add("recoToSimAssociatorCP", - edm::InputTag("tracksterSimTracksterAssociationLinkingbyCLUE3D", "recoToSim")); - desc.add("simToRecoAssociatorCP", - edm::InputTag("tracksterSimTracksterAssociationLinkingbyCLUE3D", "simToReco")); - desc.add("MergerecoToSimAssociatorSC", - edm::InputTag("tracksterSimTracksterAssociationPR", "recoToSim")); - desc.add("MergesimToRecoAssociatorSC", - edm::InputTag("tracksterSimTracksterAssociationPR", "simToReco")); - desc.add("MergerecoToSimAssociatorCP", - edm::InputTag("tracksterSimTracksterAssociationLinking", "recoToSim")); - desc.add("MergesimToRecoAssociatorCP", - edm::InputTag("tracksterSimTracksterAssociationLinking", "simToReco")); - desc.add("MergerecoToSimAssociatorPU", - edm::InputTag("tracksterSimTracksterAssociationLinkingPU", "recoToSim")); - desc.add("MergesimToRecoAssociatorPU", - edm::InputTag("tracksterSimTracksterAssociationLinkingPU", "simToReco")); + + // Settings for dumping trackster associators (recoToSim & simToReco) + edm::ParameterSetDescription associatorDescValidator; + associatorDescValidator.add("branchName")->setComment("Name of the output branches in the tree"); + associatorDescValidator.add("suffix")->setComment("Should be CP or SC (for the output branch name)"); + associatorDescValidator.add("associatorInputTag") + ->setComment("Input tag for the associator (do not put the instance)"); + associatorDescValidator.add("tracksterCollection") + ->setComment("Collection of tracksters used by the associator"); + associatorDescValidator.add("simTracksterCollection") + ->setComment("Collection of SimTrackster used by the associator"); + desc.addVPSet("associators", associatorDescValidator)->setComment("Tracksters to SimTracksters associators to dump"); + desc.add("simclusters", edm::InputTag("mix", "MergedCaloTruth")); desc.add("caloparticles", edm::InputTag("mix", "MergedCaloTruth")); desc.add("detector", "HGCAL"); desc.add("propagator", "PropagatorWithMaterial"); desc.add("saveLCs", true); - desc.add("saveCLUE3DTracksters", true); - desc.add("saveTrackstersMerged", true); - desc.add("saveSimTrackstersSC", true); - desc.add("saveSimTrackstersCP", true); desc.add("saveTICLCandidate", true); desc.add("saveSimTICLCandidate", true); desc.add("saveTracks", true); - desc.add("saveAssociations", true); + desc.add("saveSuperclustering", true); + desc.add("saveRecoSuperclusters", true) + ->setComment("Save superclustering Egamma collections (as reco::SuperCluster)"); descriptions.add("ticlDumper", desc); } diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingPassthrough.h b/RecoHGCal/TICL/plugins/TracksterLinkingPassthrough.h index d66f6015caa44..dd719e8973305 100644 --- a/RecoHGCal/TICL/plugins/TracksterLinkingPassthrough.h +++ b/RecoHGCal/TICL/plugins/TracksterLinkingPassthrough.h @@ -11,8 +11,10 @@ namespace ticl { class TracksterLinkingPassthrough : public TracksterLinkingAlgoBase { public: - TracksterLinkingPassthrough(const edm::ParameterSet& conf, edm::ConsumesCollector iC) - : TracksterLinkingAlgoBase(conf, iC) {} + TracksterLinkingPassthrough(const edm::ParameterSet& conf, + edm::ConsumesCollector iC, + cms::Ort::ONNXRuntime const* onnxRuntime = nullptr) + : TracksterLinkingAlgoBase(conf, iC, onnxRuntime) {} ~TracksterLinkingPassthrough() override {} @@ -20,7 +22,6 @@ namespace ticl { std::vector& resultTracksters, std::vector>& linkedResultTracksters, std::vector>& linkedTracksterIdToInputTracksterId) override; - void initialize(const HGCalDDDConstants* hgcons, const hgcal::RecHitTools rhtools, const edm::ESHandle bfieldH, diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.cc b/RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.cc index b6027ef4ebe6c..bacc9fefa31de 100644 --- a/RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.cc +++ b/RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.cc @@ -1,13 +1,19 @@ -// #include "TracksterLinkingbySuperClustering.h" #include "FWCore/ParameterSet/interface/ValidatedPluginFactoryMacros.h" #include "FWCore/ParameterSet/interface/ValidatedPluginMacros.h" #include "TracksterLinkingbyFastJet.h" +#include "TracksterLinkingbySuperClusteringDNN.h" +#include "TracksterLinkingbySuperClusteringMustache.h" #include "TracksterLinkingbySkeletons.h" #include "TracksterLinkingPassthrough.h" #include "RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.h" EDM_REGISTER_VALIDATED_PLUGINFACTORY(TracksterLinkingPluginFactory, "TracksterLinkingPluginFactory"); DEFINE_EDM_VALIDATED_PLUGIN(TracksterLinkingPluginFactory, ticl::TracksterLinkingbySkeletons, "Skeletons"); -// DEFINE_EDM_VALIDATED_PLUGIN(TracksterLinkingPluginFactory, ticl::TracksterLinkingbySuperClustering, "SuperClustering"); +DEFINE_EDM_VALIDATED_PLUGIN(TracksterLinkingPluginFactory, + ticl::TracksterLinkingbySuperClusteringDNN, + "SuperClusteringDNN"); +DEFINE_EDM_VALIDATED_PLUGIN(TracksterLinkingPluginFactory, + ticl::TracksterLinkingbySuperClusteringMustache, + "SuperClusteringMustache"); DEFINE_EDM_VALIDATED_PLUGIN(TracksterLinkingPluginFactory, ticl::TracksterLinkingbyFastJet, "FastJet"); DEFINE_EDM_VALIDATED_PLUGIN(TracksterLinkingPluginFactory, ticl::TracksterLinkingPassthrough, "Passthrough"); diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.h b/RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.h index 1b9cc19ff0122..0e4c4ee2fea7b 100644 --- a/RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.h +++ b/RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.h @@ -6,7 +6,13 @@ #include "FWCore/Framework/interface/ConsumesCollector.h" #include "RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h" -using TracksterLinkingPluginFactory = - edmplugin::PluginFactory; +namespace cms { + namespace Ort { + class ONNXRuntime; + } +} // namespace cms + +using TracksterLinkingPluginFactory = edmplugin::PluginFactory; #endif diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingbyFastJet.h b/RecoHGCal/TICL/plugins/TracksterLinkingbyFastJet.h index d2389bf15b6d9..e298aad6b6eb3 100644 --- a/RecoHGCal/TICL/plugins/TracksterLinkingbyFastJet.h +++ b/RecoHGCal/TICL/plugins/TracksterLinkingbyFastJet.h @@ -13,7 +13,9 @@ namespace ticl { class TracksterLinkingbyFastJet : public TracksterLinkingAlgoBase { public: - TracksterLinkingbyFastJet(const edm::ParameterSet& conf, edm::ConsumesCollector iC) + TracksterLinkingbyFastJet(const edm::ParameterSet& conf, + edm::ConsumesCollector iC, + cms::Ort::ONNXRuntime const* onnxRuntime = nullptr) : TracksterLinkingAlgoBase(conf, iC), radius_(conf.getParameter("radius")) { // Cluster tracksters into jets using FastJet with configurable algorithm auto algo = conf.getParameter("jet_algorithm"); diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.cc b/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.cc index 8f4df46948af9..f41a99e0a1c88 100644 --- a/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.cc +++ b/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.cc @@ -50,7 +50,9 @@ namespace { using namespace ticl; -TracksterLinkingbySkeletons::TracksterLinkingbySkeletons(const edm::ParameterSet &conf, edm::ConsumesCollector iC) +TracksterLinkingbySkeletons::TracksterLinkingbySkeletons(const edm::ParameterSet &conf, + edm::ConsumesCollector iC, + cms::Ort::ONNXRuntime const *onnxRuntime) : TracksterLinkingAlgoBase(conf, iC), timing_quality_threshold_(conf.getParameter("track_time_quality_threshold")), del_(conf.getParameter("wind")), diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.h b/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.h index d357889305bca..600d86f9ada21 100644 --- a/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.h +++ b/RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.h @@ -23,7 +23,9 @@ namespace ticl { class TracksterLinkingbySkeletons : public TracksterLinkingAlgoBase { public: - TracksterLinkingbySkeletons(const edm::ParameterSet& conf, edm::ConsumesCollector iC); + TracksterLinkingbySkeletons(const edm::ParameterSet& conf, + edm::ConsumesCollector iC, + cms::Ort::ONNXRuntime const* onnxRuntime = nullptr); ~TracksterLinkingbySkeletons() override {} diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringDNN.cc b/RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringDNN.cc new file mode 100644 index 0000000000000..49ad7f0c6b9e3 --- /dev/null +++ b/RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringDNN.cc @@ -0,0 +1,382 @@ +/* +TICL plugin for electron superclustering in HGCAL using a DNN. +DNN designed by Alessandro Tarabini. + +Inputs are CLUE3D EM tracksters. Outputs are superclusters (as vectors of IDs of trackster) +"Seed trackster" : seed of supercluster, always highest pT trackster of supercluster, normally should be an electron +"Candidate trackster" : trackster that is considered for superclustering with a seed + +Algorithm description : +1) Tracksters are ordered by decreasing pT. +2) We iterate over candidate tracksters, then over seed tracksters with higher pT than the candidate. + If the pair seed-candidate is in a compatible eta-phi window and passes some selections (seed pT, energy, etc), then we add the DNN features of the pair to a tensor for later inference. +3) We run the inference with the DNN on the pairs (in minibatches to reduce memory usage) +4) We iterate over candidate and seed pairs inference results. For each candidate, we take the seed for which the DNN score for the seed-candidate score is best. + If the score is also above a working point, then we add the candidate to the supercluster of the seed, and mask the candidate so it cannot be considered as a seed further + +The loop is first on candidate, then on seeds as it is more efficient for step 4 to find the best seed for each candidate. + +Authors : Theo Cuisset , Shamik Ghosh +Date : 11/2023 +*/ + +#include +#include +#include +#include + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/ParameterSet/interface/allowedValues.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Utilities/interface/FileInPath.h" +#include "PhysicsTools/ONNXRuntime/interface/ONNXRuntime.h" + +#include "DataFormats/HGCalReco/interface/TICLLayerTile.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" +#include "RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringDNN.h" + +using namespace ticl; + +TracksterLinkingbySuperClusteringDNN::TracksterLinkingbySuperClusteringDNN(const edm::ParameterSet& ps, + edm::ConsumesCollector iC, + cms::Ort::ONNXRuntime const* onnxRuntime) + : TracksterLinkingAlgoBase(ps, iC, onnxRuntime), + dnnInputs_(makeSuperclusteringDNNInputFromString(ps.getParameter("dnnInputsVersion"))), + inferenceBatchSize_(ps.getParameter("inferenceBatchSize")), + nnWorkingPoint_(ps.getParameter("nnWorkingPoint")), + deltaEtaWindow_(ps.getParameter("deltaEtaWindow")), + deltaPhiWindow_(ps.getParameter("deltaPhiWindow")), + seedPtThreshold_(ps.getParameter("seedPtThreshold")), + candidateEnergyThreshold_(ps.getParameter("candidateEnergyThreshold")), + explVarRatioCut_energyBoundary_(ps.getParameter("candidateEnergyThreshold")), + explVarRatioMinimum_lowEnergy_(ps.getParameter("explVarRatioMinimum_lowEnergy")), + explVarRatioMinimum_highEnergy_(ps.getParameter("explVarRatioMinimum_highEnergy")), + filterByTracksterPID_(ps.getParameter("filterByTracksterPID")), + tracksterPIDCategoriesToFilter_(ps.getParameter>("tracksterPIDCategoriesToFilter")), + PIDThreshold_(ps.getParameter("PIDThreshold")) { + assert(onnxRuntime_ && + "TracksterLinkingbySuperClusteringDNN : ONNXRuntime was not provided, the model should have been set in " + "onnxModelPath in the plugin config"); +} + +void TracksterLinkingbySuperClusteringDNN::initialize(const HGCalDDDConstants* hgcons, + const hgcal::RecHitTools rhtools, + const edm::ESHandle bfieldH, + const edm::ESHandle propH) {} + +/** + * Check if trackster passes cut on explained variance ratio. The DNN is trained only on pairs where both seed and candidate pass this cut + * Explained variance ratio is (largest PCA eigenvalue) / (sum of PCA eigenvalues) +*/ +bool TracksterLinkingbySuperClusteringDNN::checkExplainedVarianceRatioCut(ticl::Trackster const& ts) const { + float explVar_denominator = + std::accumulate(std::begin(ts.eigenvalues()), std::end(ts.eigenvalues()), 0.f, std::plus()); + if (explVar_denominator != 0.) { + float explVarRatio = ts.eigenvalues()[0] / explVar_denominator; + if (ts.raw_energy() > explVarRatioCut_energyBoundary_) + return explVarRatio >= explVarRatioMinimum_highEnergy_; + else + return explVarRatio >= explVarRatioMinimum_lowEnergy_; + } else + return false; +} + +bool TracksterLinkingbySuperClusteringDNN::trackstersPassesPIDCut(const Trackster& tst) const { + if (filterByTracksterPID_) { + float probTotal = 0.0f; + for (int cat : tracksterPIDCategoriesToFilter_) { + probTotal += tst.id_probabilities(cat); + } + return probTotal >= PIDThreshold_; + } else + return true; +} + +/** + * resultTracksters : superclusters as tracksters (ie merging of tracksters that have been superclustered together) + * outputSuperclusters : same as linkedTracksterIdToInputTracksterId. Probably should use only one of the two. + * linkedTracksterIdToInputTracksterId : maps indices from resultTracksters back into input tracksters. + * resultTracksters[i] has seed input.tracksters[linkedTracksterIdToInputTracksterId[i][0]], linked with tracksters input.tracksters[linkedTracksterIdToInputTracksterId[i][1..N]] +*/ +void TracksterLinkingbySuperClusteringDNN::linkTracksters( + const Inputs& input, + std::vector& resultTracksters, + std::vector>& outputSuperclusters, + std::vector>& linkedTracksterIdToInputTracksterId) { + // For now we use all input tracksters for superclustering. At some point there might be a filter here for EM tracksters (electromagnetic identification with DNN ?) + auto const& inputTracksters = input.tracksters; + const unsigned int tracksterCount = inputTracksters.size(); + + /* Sorting tracksters by decreasing order of pT (out-of-place sort). + inputTracksters[trackstersIndicesPt[0]], ..., inputTracksters[trackstersIndicesPt[N]] makes a list of tracksters sorted by decreasing pT + Indices into this pT sorted collection will have the suffix _pt. Thus inputTracksters[index] and inputTracksters[trackstersIndicesPt[index_pt]] are correct + */ + std::vector trackstersIndicesPt(inputTracksters.size()); + std::iota(trackstersIndicesPt.begin(), trackstersIndicesPt.end(), 0); + std::stable_sort( + trackstersIndicesPt.begin(), trackstersIndicesPt.end(), [&inputTracksters](unsigned int i1, unsigned int i2) { + return inputTracksters[i1].raw_pt() > inputTracksters[i2].raw_pt(); + }); + + /* Evaluate in minibatches since running with trackster count = 3000 leads to a short-lived ~15GB memory allocation + Also we do not know in advance how many superclustering candidate pairs there are going to be + The batch size needs to be rounded to featureCount + */ + const unsigned int miniBatchSize = + static_cast(inferenceBatchSize_) / dnnInputs_->featureCount() * dnnInputs_->featureCount(); + + std::vector> + inputTensorBatches; // DNN input features tensors, in minibatches. Outer array : minibatches, inner array : 2D (flattened) array of features (indexed by batchIndex, featureId) + // How far along in the latest tensor of inputTensorBatches are we. Set to miniBatchSize to trigger the creation of the tensor batch on first run + unsigned int candidateIndexInCurrentBatch = miniBatchSize; + // List of all (ts_seed_id; ts_cand_id) selected for DNN inference (same layout as inputTensorBatches) + // Index is in global trackster collection (not pt ordered collection) + std::vector>> tracksterIndicesUsedInDNN; + + // Use TracksterTiles to speed up search of tracksters in eta-phi window. One per endcap + std::array tracksterTilesBothEndcaps_pt; // one per endcap + for (unsigned int i_pt = 0; i_pt < trackstersIndicesPt.size(); ++i_pt) { + Trackster const& ts = inputTracksters[trackstersIndicesPt[i_pt]]; + tracksterTilesBothEndcaps_pt[ts.barycenter().eta() > 0.].fill(ts.barycenter().eta(), ts.barycenter().phi(), i_pt); + } + + // First loop on candidate tracksters (start at 1 since the highest pt trackster can only be a seed, not a candidate) + for (unsigned int ts_cand_idx_pt = 1; ts_cand_idx_pt < tracksterCount; ts_cand_idx_pt++) { + Trackster const& ts_cand = inputTracksters[trackstersIndicesPt[ts_cand_idx_pt]]; + + if (ts_cand.raw_energy() < candidateEnergyThreshold_ || + !checkExplainedVarianceRatioCut(ts_cand)) // || !trackstersPassesPIDCut(ts_cand) + continue; + + auto& tracksterTiles = tracksterTilesBothEndcaps_pt[ts_cand.barycenter().eta() > 0]; + std::array search_box = tracksterTiles.searchBoxEtaPhi(ts_cand.barycenter().Eta() - deltaEtaWindow_, + ts_cand.barycenter().Eta() + deltaEtaWindow_, + ts_cand.barycenter().Phi() - deltaPhiWindow_, + ts_cand.barycenter().Phi() + deltaPhiWindow_); + // Look for seed trackster + for (int eta_i = search_box[0]; eta_i <= search_box[1]; ++eta_i) { + for (int phi_i = search_box[2]; phi_i <= search_box[3]; ++phi_i) { + for (unsigned int ts_seed_idx_pt : + tracksterTiles[tracksterTiles.globalBin(eta_i, (phi_i % TileConstants::nPhiBins))]) { + if (ts_seed_idx_pt >= ts_cand_idx_pt) + continue; // Look only at seed tracksters with higher pT than the candidate + + Trackster const& ts_seed = inputTracksters[trackstersIndicesPt[ts_seed_idx_pt]]; + + if (ts_seed.raw_pt() < seedPtThreshold_) + break; // All further seeds will have lower pT than threshold (due to pT sorting) + + if (!checkExplainedVarianceRatioCut(ts_seed) || !trackstersPassesPIDCut(ts_seed)) + continue; + + // Check that the two tracksters are geometrically compatible for superclustering + if (std::abs(ts_seed.barycenter().Eta() - ts_cand.barycenter().Eta()) < deltaEtaWindow_ && + std::abs(deltaPhi(ts_seed.barycenter().Phi(), ts_cand.barycenter().Phi())) < deltaPhiWindow_) { + if (candidateIndexInCurrentBatch >= miniBatchSize) { + // Create new minibatch + assert(candidateIndexInCurrentBatch == miniBatchSize); + + /* Estimate how many seed-candidate pairs are remaining and don't allocate a full batch in this case. Use worst-case scenario of all pairs passing geometrical window + Also assume ts_seed_idx_pt=0 (worst case) + The last tensor of inputTensorBatches will be larger than necessary. The end of it will be uninitialized, then passed to the DNN. + We do not look at the output of the DNN on this section, so it will have no consequences. + */ + inputTensorBatches.emplace_back( + std::min(miniBatchSize, + (tracksterCount * (tracksterCount - 1) - ts_cand_idx_pt * (ts_cand_idx_pt - 1)) / 2) * + dnnInputs_->featureCount()); + + candidateIndexInCurrentBatch = 0; + tracksterIndicesUsedInDNN.emplace_back(); + } + + std::vector features = dnnInputs_->computeVector(ts_seed, ts_cand); // Compute DNN features + assert(features.size() == dnnInputs_->featureCount()); + assert((candidateIndexInCurrentBatch + 1) * dnnInputs_->featureCount() <= inputTensorBatches.back().size()); + // Copy the features into the batch (TODO : could probably avoid the copy and fill straight in the batch vector) + std::copy(features.begin(), + features.end(), + inputTensorBatches.back().begin() + candidateIndexInCurrentBatch * dnnInputs_->featureCount()); + candidateIndexInCurrentBatch++; + tracksterIndicesUsedInDNN.back().emplace_back(trackstersIndicesPt[ts_seed_idx_pt], + trackstersIndicesPt[ts_cand_idx_pt]); + } + } + } + } + } + + if (inputTensorBatches.empty()) { + LogDebug("HGCalTICLSuperclustering") + << "No superclustering candidate pairs passed preselection before DNN. There are " << tracksterCount + << " tracksters in this event."; + } + +#ifdef EDM_ML_DEBUG + if (!inputTensorBatches.empty()) { + std::ostringstream s; + // Print the first 20 seed-cndidate pairs sent for inference + for (unsigned int i = 0; + i < std::min(dnnInputs_->featureCount() * 20, static_cast(inputTensorBatches[0].size())); + i++) { + s << inputTensorBatches[0][i] << " "; + if (i != 0 && i % dnnInputs_->featureCount() == 0) + s << "],\t["; + } + LogDebug("HGCalTICLSuperclustering") << inputTensorBatches.size() + << " batches were created. First batch starts as follows : [" << s.str() + << "]"; + } +#endif + + // Run the DNN inference + std::vector> + batchOutputs; // Outer index : minibatch, inner index : inference index in minibatch, value : DNN score + for (std::vector& singleBatch : inputTensorBatches) { + // ONNXRuntime takes std::vector>& as input (non-const reference) so we have to make a new vector + std::vector> inputs_for_onnx{{std::move(singleBatch)}}; + std::vector outputs = onnxRuntime_->run( + {"input"}, inputs_for_onnx, {}, {}, inputs_for_onnx[0].size() / dnnInputs_->featureCount())[0]; + batchOutputs.push_back(std::move(outputs)); + } + + /* Build mask of tracksters already superclustered as candidates (to avoid using a trackster superclustered as candidate as a seed in further iterations). + Also mask seeds (only needed to add tracksters not in a supercluster to the output). */ + std::vector tracksterMask(tracksterCount, false); + + /* Index of the seed trackster of the previous iteration + Initialized with an id that cannot be obtained in input */ + unsigned int previousCandTrackster_idx = std::numeric_limits::max(); + unsigned int bestSeedForCurrentCandidate_idx = std::numeric_limits::max(); + float bestSeedForCurrentCandidate_dnnScore = nnWorkingPoint_; + + // Lambda to be called when there is a transition from one candidate to the next (as well as after the last iteration) + // Does the actual supercluster creation + auto onCandidateTransition = [&](unsigned ts_cand_idx) { + if (bestSeedForCurrentCandidate_idx < + std::numeric_limits::max()) { // At least one seed can be superclustered with the candidate + tracksterMask[ts_cand_idx] = true; // Mask the candidate so it is not considered as seed in later iterations + + // Look for a supercluster of the seed + std::vector>::iterator seed_supercluster_it = + std::find_if(outputSuperclusters.begin(), + outputSuperclusters.end(), + [bestSeedForCurrentCandidate_idx](std::vector const& sc) { + return sc[0] == bestSeedForCurrentCandidate_idx; + }); + + if (seed_supercluster_it == outputSuperclusters.end()) { // No supercluster exists yet for the seed. Create one. + outputSuperclusters.emplace_back(std::initializer_list{bestSeedForCurrentCandidate_idx}); + resultTracksters.emplace_back(inputTracksters[bestSeedForCurrentCandidate_idx]); + linkedTracksterIdToInputTracksterId.emplace_back( + std::initializer_list{bestSeedForCurrentCandidate_idx}); + seed_supercluster_it = outputSuperclusters.end() - 1; + tracksterMask[bestSeedForCurrentCandidate_idx] = + true; // mask the seed as well (needed to find tracksters not in any supercluster) + } + // Index of the supercluster into resultTracksters, outputSuperclusters and linkedTracksterIdToInputTracksterId collections (the indices are the same) + unsigned int indexIntoOutputTracksters = seed_supercluster_it - outputSuperclusters.begin(); + seed_supercluster_it->push_back(ts_cand_idx); + resultTracksters[indexIntoOutputTracksters].mergeTracksters(inputTracksters[ts_cand_idx]); + linkedTracksterIdToInputTracksterId[indexIntoOutputTracksters].push_back(ts_cand_idx); + + assert(outputSuperclusters.size() == resultTracksters.size() && + outputSuperclusters.size() == linkedTracksterIdToInputTracksterId.size()); + assert(seed_supercluster_it->size() == linkedTracksterIdToInputTracksterId[indexIntoOutputTracksters].size()); + + bestSeedForCurrentCandidate_idx = std::numeric_limits::max(); + bestSeedForCurrentCandidate_dnnScore = nnWorkingPoint_; + } + }; + + //Iterate over minibatches + for (unsigned int batchIndex = 0; batchIndex < batchOutputs.size(); batchIndex++) { + std::vector const& currentBatchOutputs = batchOutputs[batchIndex]; // DNN score outputs + // Iterate over seed-candidate pairs inside current minibatch + for (unsigned int indexInBatch = 0; indexInBatch < tracksterIndicesUsedInDNN[batchIndex].size(); indexInBatch++) { + assert(indexInBatch < static_cast(batchOutputs[batchIndex].size())); + + const unsigned int ts_seed_idx = tracksterIndicesUsedInDNN[batchIndex][indexInBatch].first; + const unsigned int ts_cand_idx = tracksterIndicesUsedInDNN[batchIndex][indexInBatch].second; + const float currentDnnScore = currentBatchOutputs[indexInBatch]; + + if (previousCandTrackster_idx != std::numeric_limits::max() && + ts_cand_idx != previousCandTrackster_idx) { + // There is a transition from one seed to the next (don't make a transition for the first iteration) + onCandidateTransition(previousCandTrackster_idx); + } + + if (currentDnnScore > bestSeedForCurrentCandidate_dnnScore && !tracksterMask[ts_seed_idx]) { + // Check that the DNN suggests superclustering, that this seed-candidate assoc is better than previous ones, and that the seed is not already in a supercluster as candidate + bestSeedForCurrentCandidate_idx = ts_seed_idx; + bestSeedForCurrentCandidate_dnnScore = currentDnnScore; + } + previousCandTrackster_idx = ts_cand_idx; + } + } + onCandidateTransition(previousCandTrackster_idx); + + // Adding one-trackster superclusters for all tracksters not in a supercluster already that pass the seed threshold + for (unsigned int ts_id = 0; ts_id < tracksterCount; ts_id++) { + if (!tracksterMask[ts_id] && inputTracksters[ts_id].raw_pt() >= seedPtThreshold_) { + outputSuperclusters.emplace_back(std::initializer_list{ts_id}); + resultTracksters.emplace_back(inputTracksters[ts_id]); + linkedTracksterIdToInputTracksterId.emplace_back(std::initializer_list{ts_id}); + } + } + +#ifdef EDM_ML_DEBUG + for (std::vector const& sc : outputSuperclusters) { + std::ostringstream s; + for (unsigned int trackster_id : sc) + s << trackster_id << " "; + LogDebug("HGCalTICLSuperclustering") << "Created supercluster of size " << sc.size() + << " holding tracksters (first one is seed) " << s.str(); + } +#endif +} + +void TracksterLinkingbySuperClusteringDNN::fillPSetDescription(edm::ParameterSetDescription& desc) { + TracksterLinkingAlgoBase::fillPSetDescription(desc); // adds algo_verbosity + desc.add("onnxModelPath")->setComment("Path to DNN (as ONNX model)"); + desc.ifValue(edm::ParameterDescription("dnnInputsVersion", "v2", true), + edm::allowedValues("v1", "v2")) + ->setComment( + "DNN inputs version tag. Defines which set of features is fed to the DNN. Must match with the actual DNN."); + desc.add("inferenceBatchSize", 1e5) + ->setComment( + "Size of inference batches fed to DNN. Increasing it should produce faster inference but higher memory " + "usage. " + "Has no physics impact."); + desc.add("nnWorkingPoint") + ->setComment("Working point of DNN (in [0, 1]). DNN score above WP will attempt to supercluster."); + desc.add("deltaEtaWindow", 0.1) + ->setComment( + "Size of delta eta window to consider for superclustering. Seed-candidate pairs outside this window " + "are not considered for DNN inference."); + desc.add("deltaPhiWindow", 0.5) + ->setComment( + "Size of delta phi window to consider for superclustering. Seed-candidate pairs outside this window " + "are not considered for DNN inference."); + desc.add("seedPtThreshold", 4.) + ->setComment("Minimum transverse energy of trackster to be considered as seed of a supercluster"); + desc.add("candidateEnergyThreshold", 2.) + ->setComment("Minimum energy of trackster to be considered as candidate for superclustering"); + desc.add("explVarRatioCut_energyBoundary", 50.) + ->setComment("Boundary energy between low and high energy explVarRatio cut threshold"); + desc.add("explVarRatioMinimum_lowEnergy", 0.92) + ->setComment( + "Cut on explained variance ratio of tracksters to be considered as candidate, " + "for trackster raw_energy < explVarRatioCut_energyBoundary"); + desc.add("explVarRatioMinimum_highEnergy", 0.95) + ->setComment( + "Cut on explained variance ratio of tracksters to be considered as candidate, " + "for trackster raw_energy > explVarRatioCut_energyBoundary"); + desc.add("filterByTracksterPID", true)->setComment("Filter tracksters before superclustering by PID score"); + desc.add>( + "tracksterPIDCategoriesToFilter", + {static_cast(Trackster::ParticleType::photon), static_cast(Trackster::ParticleType::electron)}) + ->setComment("List of PID particle types (ticl::Trackster::ParticleType enum) to consider for PID filtering"); + desc.add("PIDThreshold", 0.8)->setComment("PID score threshold"); +} diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringDNN.h b/RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringDNN.h new file mode 100644 index 0000000000000..25e876c5576bb --- /dev/null +++ b/RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringDNN.h @@ -0,0 +1,69 @@ +/* +TICL plugin for electron superclustering in HGCAL using a DNN. +DNN designed and trained by Alessandro Tarabini. + +Inputs are CLUE3D EM tracksters. Outputs are superclusters (as vectors of IDs of trackster) +"Seed trackster" : seed of supercluster, always highest pT trackster of supercluster, normally should be an electron +"Candidate trackster" : trackster that is considered for superclustering with a seed + +Authors : Theo Cuisset , Shamik Ghosh +Date : 11/2023 +*/ + +#ifndef RecoHGCal_TICL_TracksterLinkingSuperClustering_H +#define RecoHGCal_TICL_TracksterLinkingSuperClustering_H + +#include + +namespace cms { + namespace Ort { + class ONNXRuntime; + } +} // namespace cms + +#include "RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h" +#include "RecoHGCal/TICL/interface/SuperclusteringDNNInputs.h" + +namespace ticl { + class Trackster; + + class TracksterLinkingbySuperClusteringDNN : public TracksterLinkingAlgoBase { + public: + TracksterLinkingbySuperClusteringDNN(const edm::ParameterSet& ps, + edm::ConsumesCollector iC, + cms::Ort::ONNXRuntime const* onnxRuntime = nullptr); + /* virtual */ ~TracksterLinkingbySuperClusteringDNN() override {} + static void fillPSetDescription(edm::ParameterSetDescription& iDesc); + + void linkTracksters(const Inputs& input, + std::vector& resultTracksters, + std::vector>& linkedResultTracksters, + std::vector>& linkedTracksterIdToInputTracksterId) override; + void initialize(const HGCalDDDConstants* hgcons, + const hgcal::RecHitTools rhtools, + const edm::ESHandle bfieldH, + const edm::ESHandle propH) override; + + private: + bool checkExplainedVarianceRatioCut(ticl::Trackster const& ts) const; + bool trackstersPassesPIDCut(const Trackster& ts) const; + + std::unique_ptr dnnInputs_; // Helper class for DNN input features computation + unsigned int inferenceBatchSize_; // Size of inference batches fed to DNN + double + nnWorkingPoint_; // Working point for neural network (above this score, consider the trackster candidate for superclustering) + float deltaEtaWindow_; // Delta eta window to consider trackster seed-candidate pairs for inference + float deltaPhiWindow_; // Delta phi window + float seedPtThreshold_; // Min pT for a trackster to be considered as supercluster seed + float candidateEnergyThreshold_; // Min energy for a trackster to be superclustered as candidate + float explVarRatioCut_energyBoundary_; // Boundary energy between low and high energy explVarRatio cut threshold + float explVarRatioMinimum_lowEnergy_; // Cut on explained variance ratio of tracksters to be considered as candidate, for trackster raw_energy < explVarRatioCut_energyBoundary + float explVarRatioMinimum_highEnergy_; // Cut on explained variance ratio of tracksters to be considered as candidate, for trackster raw_energy > explVarRatioCut_energyBoundary + bool filterByTracksterPID_; + std::vector tracksterPIDCategoriesToFilter_; + float PIDThreshold_; + }; + +} // namespace ticl + +#endif diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringMustache.cc b/RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringMustache.cc new file mode 100644 index 0000000000000..72a9b891823f8 --- /dev/null +++ b/RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringMustache.cc @@ -0,0 +1,136 @@ +#include +#include +#include + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/ParameterSet/interface/allowedValues.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "DataFormats/HGCalReco/interface/TICLLayerTile.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" +#include "RecoEcal/EgammaCoreTools/interface/Mustache.h" +#include "RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringMustache.h" + +using namespace ticl; + +TracksterLinkingbySuperClusteringMustache::TracksterLinkingbySuperClusteringMustache( + const edm::ParameterSet& ps, edm::ConsumesCollector iC, cms::Ort::ONNXRuntime const* onnxRuntime) + : TracksterLinkingAlgoBase(ps, iC, onnxRuntime), + ecalMustacheSCParametersToken_(iC.esConsumes()), + ecalSCDynamicDPhiParametersToken_(iC.esConsumes()), + seedThresholdPt_(ps.getParameter("seedThresholdPt")), + candidateEnergyThreshold_(ps.getParameter("candidateEnergyThreshold")), + filterByTracksterPID_(ps.getParameter("filterByTracksterPID")), + tracksterPIDCategoriesToFilter_(ps.getParameter>("tracksterPIDCategoriesToFilter")), + PIDThreshold_(ps.getParameter("PIDThreshold")) {} + +void TracksterLinkingbySuperClusteringMustache::initialize(const HGCalDDDConstants* hgcons, + const hgcal::RecHitTools rhtools, + const edm::ESHandle bfieldH, + const edm::ESHandle propH) {} + +void TracksterLinkingbySuperClusteringMustache::setEvent(edm::Event& iEvent, edm::EventSetup const& iEventSetup) { + mustacheSCParams_ = &iEventSetup.getData(ecalMustacheSCParametersToken_); + scDynamicDPhiParams_ = &iEventSetup.getData(ecalSCDynamicDPhiParametersToken_); +} + +bool TracksterLinkingbySuperClusteringMustache::trackstersPassesPIDCut(const Trackster& tst) const { + if (filterByTracksterPID_) { + float probTotal = 0.0f; + for (int cat : tracksterPIDCategoriesToFilter_) { + probTotal += tst.id_probabilities(cat); + } + return probTotal >= PIDThreshold_; + } else + return true; +} + +/** + * resultTracksters : superclusters as tracksters (ie merging of tracksters that have been superclustered together) + * outputSuperclusters : same as linkedTracksterIdToInputTracksterId. Probably should use only one of the two. + * linkedTracksterIdToInputTracksterId : maps indices from resultTracksters back into input tracksters. + * resultTracksters[i] has seed input.tracksters[linkedTracksterIdToInputTracksterId[i][0]], linked with tracksters input.tracksters[linkedTracksterIdToInputTracksterId[i][1..N]] +*/ +void TracksterLinkingbySuperClusteringMustache::linkTracksters( + const Inputs& input, + std::vector& resultTracksters, + std::vector>& outputSuperclusters, + std::vector>& linkedTracksterIdToInputTracksterId) { + // For now we use all input tracksters for superclustering. At some point there might be a filter here for EM tracksters (electromagnetic identification with DNN ?) + auto const& inputTracksters = input.tracksters; + const unsigned int tracksterCount = inputTracksters.size(); + + /* Sorting tracksters by decreasing order of pT (out-of-place sort). + inputTracksters[trackstersIndicesPt[0]], ..., inputTracksters[trackstersIndicesPt[N]] makes a list of tracksters sorted by decreasing pT + Indices into this pT sorted collection will have the suffix _pt. Thus inputTracksters[index] and inputTracksters[trackstersIndicesPt[index_pt]] are correct + */ + std::vector trackstersIndicesPt(inputTracksters.size()); + std::iota(trackstersIndicesPt.begin(), trackstersIndicesPt.end(), 0); + std::stable_sort( + trackstersIndicesPt.begin(), trackstersIndicesPt.end(), [&inputTracksters](unsigned int i1, unsigned int i2) { + return inputTracksters[i1].raw_pt() > inputTracksters[i2].raw_pt(); + }); + + std::vector tracksterMask_pt(tracksterCount, false); // Mask for already superclustered tracksters + // We also mask tracksters that don't pass the PID cut + for (unsigned int ts_idx_pt = 0; ts_idx_pt < tracksterCount; ts_idx_pt++) { + tracksterMask_pt[ts_idx_pt] = !trackstersPassesPIDCut(inputTracksters[trackstersIndicesPt[ts_idx_pt]]); + } + + for (unsigned int ts_seed_idx_pt = 0; ts_seed_idx_pt < tracksterCount; ts_seed_idx_pt++) { + Trackster const& ts_seed = inputTracksters[trackstersIndicesPt[ts_seed_idx_pt]]; + if (ts_seed.raw_pt() <= seedThresholdPt_) + break; // Look only at seed tracksters passing threshold, take advantage of pt sorting for fast exit + if (tracksterMask_pt[ts_seed_idx_pt]) + continue; // Trackster does not pass PID cut + + outputSuperclusters.emplace_back(std::initializer_list{trackstersIndicesPt[ts_seed_idx_pt]}); + resultTracksters.emplace_back(inputTracksters[trackstersIndicesPt[ts_seed_idx_pt]]); + linkedTracksterIdToInputTracksterId.emplace_back( + std::initializer_list{trackstersIndicesPt[ts_seed_idx_pt]}); + + for (unsigned int ts_cand_idx_pt = ts_seed_idx_pt + 1; ts_cand_idx_pt < tracksterCount; ts_cand_idx_pt++) { + if (tracksterMask_pt[ts_cand_idx_pt]) + continue; // Trackster is either already superclustered or did not pass PID cut + + Trackster const& ts_cand = inputTracksters[trackstersIndicesPt[ts_cand_idx_pt]]; + + if (ts_cand.raw_energy() <= candidateEnergyThreshold_) + continue; + + const bool passes_dphi = reco::MustacheKernel::inDynamicDPhiWindow(scDynamicDPhiParams_, + ts_seed.barycenter().eta(), + ts_seed.barycenter().phi(), + ts_cand.raw_energy(), + ts_cand.barycenter().eta(), + ts_cand.barycenter().phi()); + + if (passes_dphi && reco::MustacheKernel::inMustache(mustacheSCParams_, + ts_seed.barycenter().eta(), + ts_seed.barycenter().phi(), + ts_cand.raw_energy(), + ts_cand.barycenter().eta(), + ts_cand.barycenter().phi())) { + outputSuperclusters.back().push_back(trackstersIndicesPt[ts_cand_idx_pt]); + resultTracksters.back().mergeTracksters(ts_cand); + linkedTracksterIdToInputTracksterId.back().push_back(trackstersIndicesPt[ts_cand_idx_pt]); + tracksterMask_pt[ts_cand_idx_pt] = true; + } + } + } +} + +void TracksterLinkingbySuperClusteringMustache::fillPSetDescription(edm::ParameterSetDescription& desc) { + TracksterLinkingAlgoBase::fillPSetDescription(desc); // adds algo_verbosity + desc.add("seedThresholdPt", 1.) + ->setComment("Minimum transverse energy of trackster to be considered as seed of a supercluster"); + desc.add("candidateEnergyThreshold", 0.15) + ->setComment("Minimum energy of trackster to be considered as candidate for superclustering"); + desc.add("filterByTracksterPID", true)->setComment("Filter tracksters before superclustering by PID score"); + desc.add>( + "tracksterPIDCategoriesToFilter", + {static_cast(Trackster::ParticleType::photon), static_cast(Trackster::ParticleType::electron)}) + ->setComment("List of PID particle types (ticl::Trackster::ParticleType enum) to consider for PID filtering"); + desc.add("PIDThreshold", 0.8)->setComment("PID score threshold"); +} diff --git a/RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringMustache.h b/RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringMustache.h new file mode 100644 index 0000000000000..4ceabc4de6e7f --- /dev/null +++ b/RecoHGCal/TICL/plugins/TracksterLinkingbySuperClusteringMustache.h @@ -0,0 +1,65 @@ +/* +TICL plugin for electron superclustering in HGCAL with Mustache algorithm +Authors : Theo Cuisset , Shamik Ghosh +Date : 06/2024 +*/ + +#ifndef RecoHGCal_TICL_TracksterLinkingSuperClusteringMustache_H +#define RecoHGCal_TICL_TracksterLinkingSuperClusteringMustache_H + +#include + +namespace cms { + namespace Ort { + class ONNXRuntime; + } +} // namespace cms + +#include "RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h" +#include "RecoHGCal/TICL/interface/SuperclusteringDNNInputs.h" + +#include "CondFormats/EcalObjects/interface/EcalMustacheSCParameters.h" +#include "CondFormats/DataRecord/interface/EcalMustacheSCParametersRcd.h" +#include "CondFormats/EcalObjects/interface/EcalSCDynamicDPhiParameters.h" +#include "CondFormats/DataRecord/interface/EcalSCDynamicDPhiParametersRcd.h" + +namespace ticl { + class Trackster; + + class TracksterLinkingbySuperClusteringMustache : public TracksterLinkingAlgoBase { + public: + TracksterLinkingbySuperClusteringMustache(const edm::ParameterSet& ps, + edm::ConsumesCollector iC, + cms::Ort::ONNXRuntime const* onnxRuntime = nullptr); + /* virtual */ ~TracksterLinkingbySuperClusteringMustache() override {} + static void fillPSetDescription(edm::ParameterSetDescription& iDesc); + + void linkTracksters(const Inputs& input, + std::vector& resultTracksters, + std::vector>& linkedResultTracksters, + std::vector>& linkedTracksterIdToInputTracksterId) override; + void initialize(const HGCalDDDConstants* hgcons, + const hgcal::RecHitTools rhtools, + const edm::ESHandle bfieldH, + const edm::ESHandle propH) override; + + virtual void setEvent(edm::Event& iEvent, edm::EventSetup const& iEventSetup) override; + + private: + bool trackstersPassesPIDCut(const Trackster& ts) const; + + edm::ESGetToken ecalMustacheSCParametersToken_; + edm::ESGetToken ecalSCDynamicDPhiParametersToken_; + const EcalMustacheSCParameters* mustacheSCParams_; + const EcalSCDynamicDPhiParameters* scDynamicDPhiParams_; + + float seedThresholdPt_; + float candidateEnergyThreshold_; + bool filterByTracksterPID_; + std::vector tracksterPIDCategoriesToFilter_; + float PIDThreshold_; + }; + +} // namespace ticl + +#endif diff --git a/RecoHGCal/TICL/plugins/TracksterLinksProducer.cc b/RecoHGCal/TICL/plugins/TracksterLinksProducer.cc index f929e38a525af..7cc7f58e46585 100644 --- a/RecoHGCal/TICL/plugins/TracksterLinksProducer.cc +++ b/RecoHGCal/TICL/plugins/TracksterLinksProducer.cc @@ -25,6 +25,7 @@ #include "PhysicsTools/TensorFlow/interface/TensorFlow.h" #include "PhysicsTools/TensorFlow/interface/TfGraphDefWrapper.h" #include "PhysicsTools/TensorFlow/interface/TensorFlow.h" +#include "PhysicsTools/ONNXRuntime/interface/ONNXRuntime.h" #include "RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h" #include "RecoHGCal/TICL/plugins/TracksterLinkingPluginFactory.h" @@ -45,15 +46,18 @@ #include "TrackstersPCA.h" using namespace ticl; +using cms::Ort::ONNXRuntime; -class TracksterLinksProducer : public edm::stream::EDProducer<> { +class TracksterLinksProducer : public edm::stream::EDProducer> { public: - explicit TracksterLinksProducer(const edm::ParameterSet &ps); + explicit TracksterLinksProducer(const edm::ParameterSet &ps, const ONNXRuntime *); ~TracksterLinksProducer() override{}; void produce(edm::Event &, const edm::EventSetup &) override; static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); void beginRun(edm::Run const &iEvent, edm::EventSetup const &es) override; + static std::unique_ptr initializeGlobalCache(const edm::ParameterSet &iConfig); + static void globalEndJob(const ONNXRuntime *); private: void printTrackstersDebug(const std::vector &, const char *label) const; @@ -63,6 +67,7 @@ class TracksterLinksProducer : public edm::stream::EDProducer<> { std::vector &result) const; std::unique_ptr linkingAlgo_; + std::string algoType_; std::vector>> tracksters_tokens_; const edm::EDGetTokenT> clusters_token_; @@ -94,8 +99,9 @@ class TracksterLinksProducer : public edm::stream::EDProducer<> { edm::ESGetToken hdc_token_; }; -TracksterLinksProducer::TracksterLinksProducer(const edm::ParameterSet &ps) - : clusters_token_(consumes>(ps.getParameter("layer_clusters"))), +TracksterLinksProducer::TracksterLinksProducer(const edm::ParameterSet &ps, const ONNXRuntime *onnxRuntime) + : algoType_(ps.getParameter("linkingPSet").getParameter("type")), + clusters_token_(consumes>(ps.getParameter("layer_clusters"))), clustersTime_token_( consumes>>(ps.getParameter("layer_clustersTime"))), regressionAndPid_(ps.getParameter("regressionAndPid")), @@ -129,24 +135,36 @@ TracksterLinksProducer::TracksterLinksProducer(const edm::ParameterSet &ps) // Links produces>>(); + produces>>("linkedTracksterIdToInputTracksterId"); // LayerClusters Mask produces>(); auto linkingPSet = ps.getParameter("linkingPSet"); - auto algoType = linkingPSet.getParameter("type"); - if (algoType == "Skeletons") { + if (algoType_ == "Skeletons") { std::string detectorName_ = (detector_ == "HFNose") ? "HGCalHFNoseSensitive" : "HGCalEESensitive"; hdc_token_ = esConsumes( edm::ESInputTag("", detectorName_)); } - linkingAlgo_ = TracksterLinkingPluginFactory::get()->create(algoType, linkingPSet, consumesCollector()); + linkingAlgo_ = TracksterLinkingPluginFactory::get()->create(algoType_, linkingPSet, consumesCollector(), onnxRuntime); } +std::unique_ptr TracksterLinksProducer::initializeGlobalCache(const edm::ParameterSet &iConfig) { + auto const &pluginPset = iConfig.getParameter("linkingPSet"); + if (pluginPset.exists("onnxModelPath")) + return std::make_unique(pluginPset.getParameter("onnxModelPath").fullPath()); + else + return std::unique_ptr(nullptr); +} + +void TracksterLinksProducer::globalEndJob(const ONNXRuntime *) {} + void TracksterLinksProducer::beginRun(edm::Run const &iEvent, edm::EventSetup const &es) { - edm::ESHandle hdc = es.getHandle(hdc_token_); - hgcons_ = hdc.product(); + if (algoType_ == "Skeletons") { + edm::ESHandle hdc = es.getHandle(hdc_token_); + hgcons_ = hdc.product(); + } edm::ESHandle geom = es.getHandle(geometry_token_); rhtools_.setGeometry(*geom); @@ -306,6 +324,8 @@ void TracksterLinksProducer::energyRegressionAndID(const std::vectorsetEvent(evt, es); + auto resultTracksters = std::make_unique>(); auto linkedResultTracksters = std::make_unique>>(); @@ -336,14 +356,14 @@ void TracksterLinksProducer::produce(edm::Event &evt, const edm::EventSetup &es) // Linking const typename TracksterLinkingAlgoBase::Inputs input(evt, es, layerClusters, layerClustersTimes, trackstersManager); - std::vector> linkedTracksterIdToInputTracksterId; + auto linkedTracksterIdToInputTracksterId = std::make_unique>>(); // LinkTracksters will produce a vector of vector of indices of tracksters that: // 1) are linked together if more than one // 2) are isolated if only one // Result tracksters contains the final version of the trackster collection // linkedTrackstersToInputTrackstersMap contains the mapping between the linked tracksters and the input tracksters - linkingAlgo_->linkTracksters(input, *resultTracksters, *linkedResultTracksters, linkedTracksterIdToInputTracksterId); + linkingAlgo_->linkTracksters(input, *resultTracksters, *linkedResultTracksters, *linkedTracksterIdToInputTracksterId); // Now we need to remove the tracksters that are not linked // We need to emplace_back in the resultTracksters only the tracksters that are linked @@ -357,12 +377,17 @@ void TracksterLinksProducer::produce(edm::Event &evt, const edm::EventSetup &es) if (regressionAndPid_) energyRegressionAndID(layerClusters, tfSession_, *resultTracksters); - assignPCAtoTracksters( - *resultTracksters, layerClusters, layerClustersTimes, rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z(), rhtools_, true); + assignPCAtoTracksters(*resultTracksters, + layerClusters, + layerClustersTimes, + rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z(), + rhtools_, + true); evt.put(std::move(linkedResultTracksters)); evt.put(std::move(resultMask)); evt.put(std::move(resultTracksters)); + evt.put(std::move(linkedTracksterIdToInputTracksterId), "linkedTracksterIdToInputTracksterId"); } void TracksterLinksProducer::printTrackstersDebug(const std::vector &tracksters, const char *label) const { diff --git a/RecoHGCal/TICL/python/CLUE3DHighStep_cff.py b/RecoHGCal/TICL/python/CLUE3DHighStep_cff.py index 9228830d659c1..177dc76cf7a86 100644 --- a/RecoHGCal/TICL/python/CLUE3DHighStep_cff.py +++ b/RecoHGCal/TICL/python/CLUE3DHighStep_cff.py @@ -32,7 +32,6 @@ from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5 ticl_v5.toModify(ticlTrackstersCLUE3DHigh.pluginPatternRecognitionByCLUE3D, computeLocalTime = cms.bool(True)) -ticl_v5.toModify(ticlTrackstersCLUE3DHigh.pluginPatternRecognitionByCLUE3D, doPidCut = cms.bool(False)) ticl_v5.toModify(ticlTrackstersCLUE3DHigh.pluginPatternRecognitionByCLUE3D, usePCACleaning = cms.bool(True)) ticlCLUE3DHighStepTask = cms.Task(ticlSeedingGlobal diff --git a/RecoHGCal/TICL/python/customiseForTICLv5_cff.py b/RecoHGCal/TICL/python/customiseForTICLv5_cff.py index d548d5d434b5b..f6daef64bb325 100644 --- a/RecoHGCal/TICL/python/customiseForTICLv5_cff.py +++ b/RecoHGCal/TICL/python/customiseForTICLv5_cff.py @@ -15,10 +15,11 @@ from RecoHGCal.TICL.tracksterSelectionTf_cfi import * from RecoHGCal.TICL.tracksterLinksProducer_cfi import tracksterLinksProducer as _tracksterLinksProducer +from RecoHGCal.TICL.ticlEGammaSuperClusterProducer_cfi import ticlEGammaSuperClusterProducer from RecoHGCal.TICL.ticlCandidateProducer_cfi import ticlCandidateProducer as _ticlCandidateProducer from RecoHGCal.Configuration.RecoHGCal_EventContent_cff import customiseForTICLv5EventContent from RecoHGCal.TICL.iterativeTICL_cff import ticlIterLabels, ticlIterLabelsMerge -from RecoHGCal.TICL.ticlDumper_cfi import ticlDumper +from RecoHGCal.TICL.ticlDumper_cff import ticlDumper from RecoHGCal.TICL.mergedTrackstersProducer_cfi import mergedTrackstersProducer as _mergedTrackstersProducer from SimCalorimetry.HGCalAssociatorProducers.TSToSimTSAssociation_cfi import tracksterSimTracksterAssociationLinkingbyCLUE3D as _tracksterSimTracksterAssociationLinkingbyCLUE3D from SimCalorimetry.HGCalAssociatorProducers.TSToSimTSAssociation_cfi import tracksterSimTracksterAssociationPRbyCLUE3D as _tracksterSimTracksterAssociationPRbyCLUE3D @@ -147,25 +148,75 @@ def customiseTICLv5FromReco(process, enableDumper = False): if(enableDumper): process.ticlDumper = ticlDumper.clone( - saveLCs=True, - saveCLUE3DTracksters=True, - saveTrackstersMerged=True, - saveSimTrackstersSC=True, - saveSimTrackstersCP=True, - saveTICLCandidate=True, - saveSimTICLCandidate=True, - saveTracks=True, - saveAssociations=True, - trackstersclue3d = cms.InputTag('ticlTrackstersCLUE3DHigh'), - ticlcandidates = cms.InputTag("ticlCandidate"), - trackstersmerged = cms.InputTag("ticlCandidate"), - trackstersInCand = cms.InputTag("ticlCandidate") + tracksterCollections=[ + cms.PSet( + treeName=cms.string("trackstersCLUE3DHigh"), + inputTag=cms.InputTag("ticlTrackstersCLUE3DHigh") + ), + cms.PSet( + treeName=cms.string("trackstersTiclCandidate"), + inputTag=cms.InputTag("ticlCandidate") + ), + + cms.PSet( + treeName=cms.string("simtrackstersSC"), + inputTag=cms.InputTag("ticlSimTracksters"), + tracksterType=cms.string("SimTracksterSC") + ), + cms.PSet( + treeName=cms.string("simtrackstersCP"), + inputTag=cms.InputTag("ticlSimTracksters", "fromCPs"), + tracksterType=cms.string("SimTracksterCP") + ), + ], + ticlcandidates=cms.InputTag("ticlCandidate"), + trackstersInCand=cms.InputTag("ticlCandidate"), + recoSuperClusters_sourceTracksterCollection = cms.InputTag("ticlCandidate"), + associators=[ + cms.PSet( + branchName=cms.string("tsCLUE3D"), + suffix=cms.string("SC"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationPRbyCLUE3D"), + tracksterCollection=cms.InputTag("ticlTrackstersCLUE3DHigh"), + simTracksterCollection=cms.InputTag("ticlSimTracksters") + ), + cms.PSet( + branchName=cms.string("tsCLUE3D"), + suffix=cms.string("CP"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationLinkingbyCLUE3D"), + tracksterCollection=cms.InputTag("ticlTrackstersCLUE3DHigh"), + simTracksterCollection=cms.InputTag("ticlSimTracksters", "fromCPs") + ), + + cms.PSet( + branchName=cms.string("ticlCandidate"), + suffix=cms.string("SC"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationLinking"), + tracksterCollection=cms.InputTag("ticlCandidate"), + simTracksterCollection=cms.InputTag("ticlSimTracksters", "fromCPs") + ), + cms.PSet( + branchName=cms.string("ticlCandidate"), + suffix=cms.string("CP"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationPR"), + tracksterCollection=cms.InputTag("ticlCandidate"), + simTracksterCollection=cms.InputTag("ticlSimTracksters") + ), + + cms.PSet( + branchName=cms.string("ticlCandidate"), + suffix=cms.string("PU"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationLinkingPU"), + tracksterCollection=cms.InputTag("ticlCandidate"), + simTracksterCollection=cms.InputTag("ticlSimTracksters", "PU") + ), + ] ) process.TFileService = cms.Service("TFileService", fileName=cms.string("histo.root") ) - process.FEVTDEBUGHLToutput_step = cms.EndPath(process.ticlDumper) + process.FEVTDEBUGHLToutput_step = cms.EndPath(process.ticlDumper) process.TICL_Validation = cms.Path(process.ticlSimTrackstersTask, process.hgcalAssociators) diff --git a/RecoHGCal/TICL/python/customiseTICLFromReco.py b/RecoHGCal/TICL/python/customiseTICLFromReco.py index 91d932d2914cf..fbd544532d10a 100644 --- a/RecoHGCal/TICL/python/customiseTICLFromReco.py +++ b/RecoHGCal/TICL/python/customiseTICLFromReco.py @@ -2,7 +2,7 @@ from RecoHGCal.TICL.iterativeTICL_cff import * from RecoLocalCalo.HGCalRecProducers.hgcalLayerClusters_cff import hgcalLayerClustersEE, hgcalLayerClustersHSi, hgcalLayerClustersHSci from RecoLocalCalo.HGCalRecProducers.hgcalMergeLayerClusters_cfi import hgcalMergeLayerClusters -from RecoHGCal.TICL.ticlDumper_cfi import ticlDumper +from RecoHGCal.TICL.ticlDumper_cff import ticlDumper # Validation from Validation.HGCalValidation.HGCalValidator_cfi import * from RecoLocalCalo.HGCalRecProducers.recHitMapProducer_cfi import recHitMapProducer @@ -15,6 +15,7 @@ from SimCalorimetry.HGCalAssociatorProducers.simTracksterAssociatorByEnergyScore_cfi import simTracksterAssociatorByEnergyScore as simTsAssocByEnergyScoreProducer from SimCalorimetry.HGCalAssociatorProducers.TSToSimTSAssociation_cfi import tracksterSimTracksterAssociationLinking, tracksterSimTracksterAssociationPR, tracksterSimTracksterAssociationLinkingbyCLUE3D, tracksterSimTracksterAssociationPRbyCLUE3D, tracksterSimTracksterAssociationLinkingPU, tracksterSimTracksterAssociationPRPU +from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5 def customiseTICLFromReco(process): # TensorFlow ESSource @@ -31,20 +32,11 @@ def customiseTICLFromReco(process): process.ticlLayerTileTask, process.ticlIterationsTask, process.ticlTracksterMergeTask) + ticl_v5.toModify(process.TICL, func=lambda x : x.associate(process.ticlTracksterLinksTask)) + # Validation process.TICL_ValidationProducers = cms.Task(process.recHitMapProducer, - process.lcAssocByEnergyScoreProducer, - process.layerClusterCaloParticleAssociationProducer, - process.scAssocByEnergyScoreProducer, - process.layerClusterSimClusterAssociationProducer, - process.simTsAssocByEnergyScoreProducer, - process.simTracksterHitLCAssociatorByEnergyScoreProducer, - process.tracksterSimTracksterAssociationLinking, - process.tracksterSimTracksterAssociationPR, - process.tracksterSimTracksterAssociationLinkingbyCLUE3D, - process.tracksterSimTracksterAssociationPRbyCLUE3D, - process.tracksterSimTracksterAssociationLinkingPU, - process.tracksterSimTracksterAssociationPRPU + process.hgcalAssociators ) process.TICL_Validator = cms.Task(process.hgcalValidator) @@ -68,25 +60,7 @@ def customiseTICLFromReco(process): def customiseTICLForDumper(process): - process.ticlDumper = ticlDumper.clone( - saveLCs=True, - saveCLUE3DTracksters=True, - saveTrackstersMerged=True, - saveSimTrackstersSC=True, - saveSimTrackstersCP=True, - saveTICLCandidate=True, - saveSimTICLCandidate=True, - saveTracks=True, - saveAssociations=True, - ) - - from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5 - ticl_v5.toModify(process.ticlDumper, - # trackstersclue3d = cms.InputTag('mergedTrackstersProducer'), # For future separate iterations - trackstersclue3d = cms.InputTag('ticlTrackstersCLUE3DHigh'), - ticlcandidates = cms.InputTag("ticlCandidate"), - trackstersmerged = cms.InputTag("ticlCandidate"), - trackstersInCand = cms.InputTag("ticlCandidate")) + process.ticlDumper = ticlDumper.clone() process.TFileService = cms.Service("TFileService", fileName=cms.string("histo.root") diff --git a/RecoHGCal/TICL/python/iterativeTICL_cff.py b/RecoHGCal/TICL/python/iterativeTICL_cff.py index ae13a46fc961f..7067041d41a30 100644 --- a/RecoHGCal/TICL/python/iterativeTICL_cff.py +++ b/RecoHGCal/TICL/python/iterativeTICL_cff.py @@ -17,6 +17,7 @@ from RecoHGCal.TICL.tracksterSelectionTf_cfi import * from RecoHGCal.TICL.tracksterLinksProducer_cfi import tracksterLinksProducer as _tracksterLinksProducer +from RecoHGCal.TICL.superclustering_cff import * from RecoHGCal.TICL.ticlCandidateProducer_cfi import ticlCandidateProducer as _ticlCandidateProducer from RecoHGCal.TICL.mtdSoAProducer_cfi import mtdSoAProducer as _mtdSoAProducer @@ -64,7 +65,8 @@ ''' ticlTracksterMergeTask = cms.Task(ticlTrackstersMerge) -ticlTracksterLinksTask = cms.Task(ticlTracksterLinks) +ticlTracksterLinksTask = cms.Task(ticlTracksterLinks, ticlSuperclusteringTask) + mergeTICLTask = cms.Task(ticlLayerTileTask ,ticlIterationsTask diff --git a/RecoHGCal/TICL/python/superclustering_cff.py b/RecoHGCal/TICL/python/superclustering_cff.py new file mode 100644 index 0000000000000..6f6bce1c143d2 --- /dev/null +++ b/RecoHGCal/TICL/python/superclustering_cff.py @@ -0,0 +1,46 @@ +import FWCore.ParameterSet.Config as cms + +from RecoHGCal.TICL.tracksterLinksProducer_cfi import tracksterLinksProducer as _tracksterLinksProducer +from RecoHGCal.TICL.ticlEGammaSuperClusterProducer_cfi import ticlEGammaSuperClusterProducer +from RecoEcal.EgammaClusterProducers.particleFlowSuperClusteringSequence_cff import particleFlowSuperClusterHGCal + +from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5 +from Configuration.ProcessModifiers.ticl_superclustering_dnn_cff import ticl_superclustering_dnn +from Configuration.ProcessModifiers.ticl_superclustering_mustache_pf_cff import ticl_superclustering_mustache_pf +from Configuration.ProcessModifiers.ticl_superclustering_mustache_ticl_cff import ticl_superclustering_mustache_ticl + +ticlTracksterLinksSuperclusteringDNN = _tracksterLinksProducer.clone( + linkingPSet = cms.PSet( + type=cms.string("SuperClusteringDNN"), + algo_verbosity=cms.int32(0), + onnxModelPath = cms.FileInPath("RecoHGCal/TICL/data/superclustering/supercls_v2p1.onnx"), + nnWorkingPoint=cms.double(0.3), + ), + tracksters_collections = [cms.InputTag("ticlTrackstersCLUE3DHigh")], # to be changed to ticlTrackstersCLUE3DEM once separate CLUE3D iterations are introduced +) + +ticlTracksterLinksSuperclusteringMustache = _tracksterLinksProducer.clone( + linkingPSet = cms.PSet( + type=cms.string("SuperClusteringMustache"), + algo_verbosity=cms.int32(0) + ), + tracksters_collections = [cms.InputTag("ticlTrackstersCLUE3DHigh")], # to be changed to ticlTrackstersCLUE3DEM once separate CLUE3D iterations are introduced +) + +### Superclustering : 3 options : DNN, Mustache-TICL (from tracksters), Mustache-PF (converting tracksters to PFClusters, default for ticl_v4, enable with modifier for v5) +ticlSuperclusteringTask = cms.Task() + +# DNN +_dnn_task = cms.Task(ticlTracksterLinksSuperclusteringDNN) +ticl_superclustering_dnn.toReplaceWith(ticlSuperclusteringTask, _dnn_task) +ticl_superclustering_dnn.toModify(ticlEGammaSuperClusterProducer, ticlSuperClusters=cms.InputTag("ticlTracksterLinksSuperclusteringDNN")) +ticl_superclustering_dnn.toReplaceWith(particleFlowSuperClusterHGCal, ticlEGammaSuperClusterProducer) + +# Mustache-TICL +_mustache_ticl_task = cms.Task(ticlTracksterLinksSuperclusteringMustache) +ticl_superclustering_mustache_ticl.toReplaceWith(ticlSuperclusteringTask, _mustache_ticl_task) +ticl_superclustering_mustache_ticl.toModify(ticlEGammaSuperClusterProducer, ticlSuperClusters=cms.InputTag("ticlTracksterLinksSuperclusteringMustache")) +ticl_superclustering_mustache_ticl.toReplaceWith(particleFlowSuperClusterHGCal, ticlEGammaSuperClusterProducer) + +# Mustache-PF +# (no changes to make) diff --git a/RecoHGCal/TICL/python/ticlDumper_cff.py b/RecoHGCal/TICL/python/ticlDumper_cff.py new file mode 100644 index 0000000000000..ae5a3880e1b13 --- /dev/null +++ b/RecoHGCal/TICL/python/ticlDumper_cff.py @@ -0,0 +1,192 @@ +import FWCore.ParameterSet.Config as cms +from RecoHGCal.TICL.ticlDumper_cfi import ticlDumper as ticlDumper_ + +from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5 +from Configuration.ProcessModifiers.ticl_superclustering_dnn_cff import ticl_superclustering_dnn +from Configuration.ProcessModifiers.ticl_superclustering_mustache_pf_cff import ticl_superclustering_mustache_pf +from Configuration.ProcessModifiers.ticl_superclustering_mustache_ticl_cff import ticl_superclustering_mustache_ticl + +ticlDumper = ticlDumper_.clone( + tracksterCollections = [ + cms.PSet( + treeName=cms.string("trackstersclue3d"), + inputTag=cms.InputTag("ticlTrackstersCLUE3DHigh") + ), + cms.PSet( + treeName=cms.string("trackstersmerged"), + inputTag=cms.InputTag("ticlTrackstersMerge") + ), + + cms.PSet( + treeName=cms.string("simtrackstersSC"), + inputTag=cms.InputTag("ticlSimTracksters"), + tracksterType=cms.string("SimTracksterSC") + ), + cms.PSet( + treeName=cms.string("simtrackstersCP"), + inputTag=cms.InputTag("ticlSimTracksters", "fromCPs"), + tracksterType=cms.string("SimTracksterCP") + ), + ], + + associators=[ + cms.PSet( + branchName=cms.string("tsCLUE3D"), + suffix=cms.string("SC"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationPRbyCLUE3D"), + tracksterCollection=cms.InputTag("ticlTrackstersCLUE3DHigh"), + simTracksterCollection=cms.InputTag("ticlSimTracksters") + ), + cms.PSet( + branchName=cms.string("tsCLUE3D"), + suffix=cms.string("CP"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationLinkingbyCLUE3D"), + tracksterCollection=cms.InputTag("ticlTrackstersCLUE3DHigh"), + simTracksterCollection=cms.InputTag("ticlSimTracksters", "fromCPs") + ), + + cms.PSet( + branchName=cms.string("Mergetracksters"), + suffix=cms.string("SC"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationPR"), + tracksterCollection=cms.InputTag("ticlTrackstersMerge"), + simTracksterCollection=cms.InputTag("ticlSimTracksters") + ), + cms.PSet( + branchName=cms.string("Mergetracksters"), + suffix=cms.string("CP"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationLinking"), + tracksterCollection=cms.InputTag("ticlTrackstersMerge"), + simTracksterCollection=cms.InputTag("ticlSimTracksters", "fromCPs") + ), + + cms.PSet( + branchName=cms.string("Mergetracksters"), + suffix=cms.string("PU"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationLinkingPU"), + tracksterCollection=cms.InputTag("ticlTrackstersMerge"), + simTracksterCollection=cms.InputTag("ticlSimTracksters", "PU") + ), + ] +) + +ticl_v5.toModify(ticlDumper, + tracksterCollections=[ + cms.PSet( + treeName=cms.string("trackstersCLUE3DHigh"), + inputTag=cms.InputTag("ticlTrackstersCLUE3DHigh") + ), + # For future separate CLUE3D iterations : + # cms.PSet( + # treeName=cms.string("trackstersCLUE3DEM"), + # inputTag=cms.InputTag("ticlTrackstersCLUE3DEM") + # ), + # cms.PSet( + # treeName=cms.string("trackstersCLUE3DHAD"), + # inputTag=cms.InputTag("ticlTrackstersCLUE3DHAD") + # ), + cms.PSet( + treeName=cms.string("trackstersTiclCandidate"), + inputTag=cms.InputTag("ticlCandidate") + ), + + cms.PSet( + treeName=cms.string("simtrackstersSC"), + inputTag=cms.InputTag("ticlSimTracksters"), + tracksterType=cms.string("SimTracksterSC") + ), + cms.PSet( + treeName=cms.string("simtrackstersCP"), + inputTag=cms.InputTag("ticlSimTracksters", "fromCPs"), + tracksterType=cms.string("SimTracksterCP") + ), + ], + ticlcandidates=cms.InputTag("ticlCandidate"), + trackstersInCand=cms.InputTag("ticlCandidate"), + recoSuperClusters_sourceTracksterCollection = cms.InputTag("ticlTrackstersCLUE3DHigh"), + associators=[ + cms.PSet( + branchName=cms.string("tsCLUE3D"), + suffix=cms.string("SC"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationPRbyCLUE3D"), + tracksterCollection=cms.InputTag("ticlTrackstersCLUE3DHigh"), + simTracksterCollection=cms.InputTag("ticlSimTracksters") + ), + cms.PSet( + branchName=cms.string("tsCLUE3D"), + suffix=cms.string("CP"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationLinkingbyCLUE3D"), + tracksterCollection=cms.InputTag("ticlTrackstersCLUE3DHigh"), + simTracksterCollection=cms.InputTag("ticlSimTracksters", "fromCPs") + ), + + cms.PSet( + branchName=cms.string("ticlCandidate"), + suffix=cms.string("SC"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationLinking"), + tracksterCollection=cms.InputTag("ticlCandidate"), + simTracksterCollection=cms.InputTag("ticlSimTracksters", "fromCPs") + ), + cms.PSet( + branchName=cms.string("ticlCandidate"), + suffix=cms.string("CP"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationPR"), + tracksterCollection=cms.InputTag("ticlCandidate"), + simTracksterCollection=cms.InputTag("ticlSimTracksters") + ), + + cms.PSet( + branchName=cms.string("ticlCandidate"), + suffix=cms.string("PU"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationLinkingPU"), + tracksterCollection=cms.InputTag("ticlCandidate"), + simTracksterCollection=cms.InputTag("ticlSimTracksters", "PU") + ), + ] +) + +(ticl_v5 & ticl_superclustering_mustache_pf).toModify(ticlDumper, saveSuperclustering=False, recoSuperClusters_sourceTracksterCollection=cms.InputTag("ticlCandidate")) +(ticl_v5 & ticl_superclustering_mustache_ticl).toModify(ticlDumper, func=lambda x: x.tracksterCollections.append(cms.PSet( + treeName=cms.string("trackstersSuperclusteringMustache"), + inputTag=cms.InputTag("ticlTracksterLinksSuperclusteringMustache") + )) +).toModify(ticlDumper, func=lambda x: x.associators.extend([ + cms.PSet( + branchName=cms.string("tsSuperclusters"), + suffix=cms.string("SC"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationPRSuperclustering"), + tracksterCollection=cms.InputTag("ticlTracksterLinksSuperclusteringMustache"), + simTracksterCollection=cms.InputTag("ticlSimTracksters") + ), + cms.PSet( + branchName=cms.string("tsSuperclusters"), + suffix=cms.string("CP"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationLinkingSuperclustering"), + tracksterCollection=cms.InputTag("ticlTracksterLinksSuperclusteringMustache"), + simTracksterCollection=cms.InputTag("ticlSimTracksters", "fromCPs") + ), +])).toModify(ticlDumper, + superclustering=cms.InputTag("ticlTracksterLinksSuperclusteringMustache") +) +(ticl_v5 & ticl_superclustering_dnn).toModify(ticlDumper, func=lambda x: x.tracksterCollections.append(cms.PSet( + treeName=cms.string("trackstersSuperclusteringDNN"), + inputTag=cms.InputTag("ticlTracksterLinksSuperclusteringDNN") + )) +).toModify(ticlDumper, func=lambda x: x.associators.extend([ + cms.PSet( + branchName=cms.string("tsSuperclusters"), + suffix=cms.string("SC"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationPRSuperclustering"), + tracksterCollection=cms.InputTag("ticlTracksterLinksSuperclusteringDNN"), + simTracksterCollection=cms.InputTag("ticlSimTracksters") + ), + cms.PSet( + branchName=cms.string("tsSuperclusters"), + suffix=cms.string("CP"), + associatorInputTag=cms.InputTag("tracksterSimTracksterAssociationLinkingSuperclustering"), + tracksterCollection=cms.InputTag("ticlTracksterLinksSuperclusteringDNN"), + simTracksterCollection=cms.InputTag("ticlSimTracksters", "fromCPs") + ), +])).toModify(ticlDumper, + superclustering=cms.InputTag("ticlTracksterLinksSuperclusteringDNN") +) diff --git a/RecoHGCal/TICL/src/SuperclusteringDNNInputs.cc b/RecoHGCal/TICL/src/SuperclusteringDNNInputs.cc new file mode 100644 index 0000000000000..be58724e8f5c9 --- /dev/null +++ b/RecoHGCal/TICL/src/SuperclusteringDNNInputs.cc @@ -0,0 +1,116 @@ +/** Computation of input features for superclustering DNN. Used by plugins/TracksterLinkingBySuperClustering.cc and plugins/SuperclusteringSampleDumper.cc */ +// Author: Theo Cuisset - theo.cuisset@cern.ch +// Date: 11/2023 +#include "RecoHGCal/TICL/interface/SuperclusteringDNNInputs.h" + +#include +#include +#include + +#include +#include +#include +#include + +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" + +namespace ticl { + + std::vector SuperclusteringDNNInputV1::computeVector(Trackster const& ts_base, Trackster const& ts_toCluster) { + /* We use the barycenter for most of the variables below as that is what seems to have been used by Alessandro Tarabini, + but using PCA might be better. + (It would need retraining of the DNN) + */ + return {{ + std::abs(ts_toCluster.barycenter().Eta()) - std::abs(ts_base.barycenter().Eta()), //DeltaEtaBaryc + ts_toCluster.barycenter().Phi() - ts_base.barycenter().phi(), //DeltaPhiBaryc + ts_toCluster.raw_energy(), //multi_en + ts_toCluster.barycenter().Eta(), //multi_eta + (ts_toCluster.raw_energy() * std::sin(ts_toCluster.barycenter().Theta())), //multi_pt + ts_base.barycenter().Eta(), //seedEta + ts_base.barycenter().Phi(), //seedPhi + ts_base.raw_energy(), //seedEn + (ts_base.raw_energy() * std::sin(ts_toCluster.barycenter().Theta())), //seedPt + }}; + } + + // Helper functions for angles. Adapted from ROOT (3D vectors -> 2D vectors) + template + float CosTheta2D(const Vector1& v1, const Vector2& v2) { + float arg; + float v1_r2 = v1.X() * v1.X() + v1.Y() * v1.Y(); + float v2_r2 = v2.X() * v2.X() + v2.Y() * v2.Y(); + float ptot2 = v1_r2 * v2_r2; + if (ptot2 <= 0) { + arg = 0.0; + } else { + float pdot = v1.X() * v2.X() + v1.Y() * v2.Y(); + using std::sqrt; + arg = pdot / sqrt(ptot2); + if (arg > 1.0) + arg = 1.0; + if (arg < -1.0) + arg = -1.0; + } + return arg; + } + template + inline float Angle2D(const Vector1& v1, const Vector2& v2) { + return static_cast(std::acos(CosTheta2D(v1, v2))); + } + + std::vector SuperclusteringDNNInputV2::computeVector(Trackster const& ts_base, Trackster const& ts_toCluster) { + using ROOT::Math::XYVectorF; + using ROOT::Math::XYZVectorF; + using ROOT::Math::VectorUtil::Angle; + XYZVectorF const& pca_seed_cmsFrame(ts_base.eigenvectors(0)); + XYZVectorF const& pca_cand_cmsFrame(ts_toCluster.eigenvectors(0)); + XYZVectorF xs(pca_seed_cmsFrame.Cross(XYZVectorF(0, 0, 1)).Unit()); + ROOT::Math::Rotation3D rot(xs, xs.Cross(pca_seed_cmsFrame).Unit(), pca_seed_cmsFrame); + + XYZVectorF pca_cand_seedFrame = rot(pca_cand_cmsFrame); // seed coordinates + + float explVar_denominator = std::accumulate( + std::begin(ts_toCluster.eigenvalues()), std::end(ts_toCluster.eigenvalues()), 0.f, std::plus()); + float explVarRatio = 0.; + if (explVar_denominator != 0.) { + explVarRatio = ts_toCluster.eigenvalues()[0] / explVar_denominator; + } else { + edm::LogWarning("HGCalTICLSuperclustering") + << "Sum of eigenvalues was zero for trackster. Could not compute explained variance ratio."; + } + + return {{ + std::abs(ts_toCluster.barycenter().Eta()) - std::abs(ts_base.barycenter().Eta()), //DeltaEtaBaryc + ts_toCluster.barycenter().Phi() - ts_base.barycenter().phi(), //DeltaPhiBaryc + ts_toCluster.raw_energy(), //multi_en + ts_toCluster.barycenter().Eta(), //multi_eta + (ts_toCluster.raw_energy() * std::sin(ts_toCluster.barycenter().Theta())), //multi_pt + ts_base.barycenter().Eta(), //seedEta + ts_base.barycenter().Phi(), //seedPhi + ts_base.raw_energy(), //seedEn + (ts_base.raw_energy() * std::sin(ts_toCluster.barycenter().Theta())), //seedPt + static_cast(Angle(pca_cand_cmsFrame, pca_seed_cmsFrame)), // theta : angle between seed and candidate + Angle2D(XYVectorF(pca_cand_seedFrame.x(), pca_cand_seedFrame.z()), XYVectorF(0, 1)), // theta_xz_seedFrame + Angle2D(XYVectorF(pca_cand_seedFrame.y(), pca_cand_seedFrame.z()), XYVectorF(0, 1)), // theta_yz_seedFrame + Angle2D(XYVectorF(pca_cand_cmsFrame.x(), pca_cand_cmsFrame.y()), + XYVectorF(pca_seed_cmsFrame.x(), pca_seed_cmsFrame.y())), // theta_xy_cmsFrame + Angle2D(XYVectorF(pca_cand_cmsFrame.y(), pca_cand_cmsFrame.z()), + XYVectorF(pca_seed_cmsFrame.y(), pca_seed_cmsFrame.z())), // theta_yz_cmsFrame + Angle2D(XYVectorF(pca_cand_cmsFrame.x(), pca_cand_cmsFrame.z()), + XYVectorF(pca_seed_cmsFrame.x(), pca_seed_cmsFrame.z())), // theta_xz_cmsFrame + ts_toCluster.eigenvalues()[0], // explVar + explVarRatio // explVarRatio + }}; + } + + std::unique_ptr makeSuperclusteringDNNInputFromString(std::string dnnInputVersion) { + if (dnnInputVersion == "v1") + return std::make_unique(); + else if (dnnInputVersion == "v2") + return std::make_unique(); + assert(false); + } + +} // namespace ticl \ No newline at end of file diff --git a/SimCalorimetry/HGCalAssociatorProducers/python/TSToSimTSAssociation_cfi.py b/SimCalorimetry/HGCalAssociatorProducers/python/TSToSimTSAssociation_cfi.py index a1c684519a968..40b4c8ea28c04 100644 --- a/SimCalorimetry/HGCalAssociatorProducers/python/TSToSimTSAssociation_cfi.py +++ b/SimCalorimetry/HGCalAssociatorProducers/python/TSToSimTSAssociation_cfi.py @@ -1,5 +1,8 @@ import FWCore.ParameterSet.Config as cms +from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5 +from Configuration.ProcessModifiers.ticl_superclustering_mustache_ticl_cff import ticl_superclustering_mustache_ticl + tracksterSimTracksterAssociationLinking = cms.EDProducer("TSToSimTSHitLCAssociatorEDProducer", associator = cms.InputTag('simTracksterHitLCAssociatorByEnergyScoreProducer'), label_tst = cms.InputTag("ticlTrackstersMerge"), @@ -37,6 +40,29 @@ label_cp = cms.InputTag("mix","MergedCaloTruth"), ) +tracksterSimTracksterAssociationLinkingSuperclustering = cms.EDProducer("TSToSimTSHitLCAssociatorEDProducer", + associator = cms.InputTag('simTracksterHitLCAssociatorByEnergyScoreProducer'), + label_tst = cms.InputTag("ticlTracksterLinksSuperclusteringDNN"), + label_simTst = cms.InputTag("ticlSimTracksters", "fromCPs"), + label_lcl = cms.InputTag("hgcalMergeLayerClusters"), + label_scl = cms.InputTag("mix", "MergedCaloTruth"), + label_cp = cms.InputTag("mix","MergedCaloTruth"), +) + +tracksterSimTracksterAssociationPRSuperclustering = cms.EDProducer("TSToSimTSHitLCAssociatorEDProducer", + associator = cms.InputTag('simTracksterHitLCAssociatorByEnergyScoreProducer'), + label_tst = cms.InputTag("ticlTracksterLinksSuperclusteringDNN"), + label_simTst = cms.InputTag("ticlSimTracksters"), + label_lcl = cms.InputTag("hgcalMergeLayerClusters"), + label_scl = cms.InputTag("mix", "MergedCaloTruth"), + label_cp = cms.InputTag("mix","MergedCaloTruth"), +) +(ticl_v5 & ticl_superclustering_mustache_ticl).toModify( + tracksterSimTracksterAssociationLinkingSuperclustering, label_tst = cms.InputTag("ticlTracksterLinksSuperclusteringMustache") +).toModify( + tracksterSimTracksterAssociationPRSuperclustering, label_tst = cms.InputTag("ticlTracksterLinksSuperclusteringMustache") +) + tracksterSimTracksterAssociationLinkingPU = cms.EDProducer("TSToSimTSHitLCAssociatorEDProducer", associator = cms.InputTag('simTracksterHitLCAssociatorByEnergyScoreProducer'), label_tst = cms.InputTag("ticlTrackstersMerge"), @@ -55,7 +81,6 @@ label_cp = cms.InputTag("mix","MergedCaloTruth"), ) -from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5 ''' For future separate iterations ticl_v5.toModify(tracksterSimTracksterAssociationLinkingbyCLUE3D, label_tst = cms.InputTag("mergedTrackstersProducer")) tracksterSimTracksterAssociationLinkingbyCLUE3DEM = tracksterSimTracksterAssociationLinkingbyCLUE3D.clone(label_tst = cms.InputTag("ticlTrackstersCLUE3DEM")) diff --git a/Validation/Configuration/python/hgcalSimValid_cff.py b/Validation/Configuration/python/hgcalSimValid_cff.py index 9a3ba580ad81b..a343bf2ac962b 100644 --- a/Validation/Configuration/python/hgcalSimValid_cff.py +++ b/Validation/Configuration/python/hgcalSimValid_cff.py @@ -9,7 +9,7 @@ from SimCalorimetry.HGCalAssociatorProducers.LCToSimTSAssociation_cfi import layerClusterSimTracksterAssociation as layerClusterSimTracksterAssociationProducer from SimCalorimetry.HGCalAssociatorProducers.LCToCPAssociation_cfi import layerClusterCaloParticleAssociationHFNose as layerClusterCaloParticleAssociationProducerHFNose from SimCalorimetry.HGCalAssociatorProducers.LCToSCAssociation_cfi import layerClusterSimClusterAssociationHFNose as layerClusterSimClusterAssociationProducerHFNose -from SimCalorimetry.HGCalAssociatorProducers.TSToSimTSAssociation_cfi import tracksterSimTracksterAssociationLinking, tracksterSimTracksterAssociationPR,tracksterSimTracksterAssociationLinkingbyCLUE3D, tracksterSimTracksterAssociationPRbyCLUE3D, tracksterSimTracksterAssociationLinkingPU, tracksterSimTracksterAssociationPRPU #, tracksterSimTracksterAssociationLinkingbyCLUE3DEM, tracksterSimTracksterAssociationLinkingbyCLUE3DHAD, tracksterSimTracksterAssociationPRbyCLUE3DEM, tracksterSimTracksterAssociationPRbyCLUE3DHAD +from SimCalorimetry.HGCalAssociatorProducers.TSToSimTSAssociation_cfi import tracksterSimTracksterAssociationLinking, tracksterSimTracksterAssociationPR,tracksterSimTracksterAssociationLinkingbyCLUE3D, tracksterSimTracksterAssociationPRbyCLUE3D, tracksterSimTracksterAssociationLinkingSuperclustering, tracksterSimTracksterAssociationPRSuperclustering, tracksterSimTracksterAssociationLinkingPU, tracksterSimTracksterAssociationPRPU #, tracksterSimTracksterAssociationLinkingbyCLUE3DEM, tracksterSimTracksterAssociationLinkingbyCLUE3DHAD, tracksterSimTracksterAssociationPRbyCLUE3DEM, tracksterSimTracksterAssociationPRbyCLUE3DHAD from RecoHGCal.TICL.mergedTrackstersProducer_cfi import mergedTrackstersProducer as _mergedTrackstersProducer from SimCalorimetry.HGCalAssociatorProducers.SimTauProducer_cfi import * @@ -45,6 +45,9 @@ ) from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5 +from Configuration.ProcessModifiers.ticl_superclustering_mustache_pf_cff import ticl_superclustering_mustache_pf +# Mustache-PF does not produce tracksters, therefore we cannot use the tracksterSimTracksterAssociation on superclusters +(ticl_v5 & ~ticl_superclustering_mustache_pf).toModify(hgcalAssociators, lambda x: x.add(tracksterSimTracksterAssociationLinkingSuperclustering, tracksterSimTracksterAssociationPRSuperclustering)) ''' For future separate iterations mergedTrackstersProducer = _mergedTrackstersProducer.clone() ticl_v5.toModify(hgcalAssociators, lambda x: x.add(mergedTrackstersProducer, tracksterSimTracksterAssociationLinkingbyCLUE3DEM, tracksterSimTracksterAssociationLinkingbyCLUE3DHAD, tracksterSimTracksterAssociationPRbyCLUE3DEM, tracksterSimTracksterAssociationPRbyCLUE3DHAD)) diff --git a/Validation/HGCalValidation/python/HGCalValidator_cfi.py b/Validation/HGCalValidation/python/HGCalValidator_cfi.py index 0c72904adbb52..71e06461823e6 100644 --- a/Validation/HGCalValidation/python/HGCalValidator_cfi.py +++ b/Validation/HGCalValidation/python/HGCalValidator_cfi.py @@ -106,6 +106,8 @@ phase2_hgcalV16.toModify(hgcalValidator, totallayers_to_monitor = cms.int32(47)) from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5 +from Configuration.ProcessModifiers.ticl_superclustering_dnn_cff import ticl_superclustering_dnn +from Configuration.ProcessModifiers.ticl_superclustering_mustache_ticl_cff import ticl_superclustering_mustache_ticl # labelTst_v5 = ["ticlTrackstersCLUE3DEM", "ticlTrackstersCLUE3DHAD", "ticlTracksterLinks"] # for separate CLUE3D iterations labelTst_v5 = ["ticlTrackstersCLUE3DHigh", "ticlTracksterLinks"] labelTst_v5.extend([cms.InputTag("ticlSimTracksters", "fromCPs"), cms.InputTag("ticlSimTracksters")]) @@ -119,3 +121,9 @@ ticlTrackstersMerge = cms.InputTag("ticlCandidate"), isticlv5 = cms.untracked.bool(True) ) +(ticl_v5 & ticl_superclustering_mustache_ticl).toModify( + hgcalValidator, label_tst = cms.VInputTag(labelTst_v5 + ["ticlTracksterLinksSuperclusteringMustache"]) +) +(ticl_v5 & ticl_superclustering_dnn).toModify( + hgcalValidator, label_tst = cms.VInputTag(labelTst_v5 + ["ticlTracksterLinksSuperclusteringDNN"]) +)