From 89b12fbc3caccda60ddc7d2275f3ddf9816377ae Mon Sep 17 00:00:00 2001 From: Felice Date: Tue, 29 Sep 2020 18:23:47 +0200 Subject: [PATCH 01/66] Created a ClusterFilter by Algo, Size and Layer Range --- .../ClusterFilterByAlgoAndSizeAndLayerRange.h | 56 +++++++++++++++++++ RecoHGCal/TICL/plugins/filters.cc | 4 ++ 2 files changed, 60 insertions(+) create mode 100644 RecoHGCal/TICL/plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h diff --git a/RecoHGCal/TICL/plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h b/RecoHGCal/TICL/plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h new file mode 100644 index 0000000000000..f1c901ade9404 --- /dev/null +++ b/RecoHGCal/TICL/plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h @@ -0,0 +1,56 @@ +// Authors: Marco Rovere - marco.rovere@cern.ch, Felice Pantaleo - felice.pantaleo@cern.ch +// Date: 09/2020 + +#ifndef RecoHGCal_TICL_ClusterFilterByAlgoAndSizeAndLayerRange_H__ +#define RecoHGCal_TICL_ClusterFilterByAlgoAndSizeAndLayerRange_H__ + +#include "DataFormats/CaloRecHit/interface/CaloCluster.h" +#include "ClusterFilterBase.h" + +#include +#include + +// Filter clusters that belong to a specific algorithm +namespace ticl { + class ClusterFilterByAlgoAndSizeAndLayerRange final : public ClusterFilterBase { + public: + ClusterFilterByAlgoAndSizeAndLayerRange(const edm::ParameterSet& ps) + : ClusterFilterBase(ps), + algo_number_(ps.getParameter("algo_number")), + min_cluster_size_(ps.getParameter("min_cluster_size")), + max_cluster_size_(ps.getParameter("max_cluster_size")), + min_layerId_(ps.getParameter("min_layerId")), + max_layerId_(ps.getParameter("max_layerId")) {} + ~ClusterFilterByAlgoAndSizeAndLayerRange() override{}; + + void filter(const std::vector& layerClusters, + const HgcalClusterFilterMask& availableLayerClusters, + std::vector& layerClustersMask, + hgcal::RecHitTools& rhtools) const override { + auto filteredLayerClusters = std::make_unique(); + for (auto const& cl : availableLayerClusters) { + auto const& layerCluster = layerClusters[cl.first]; + auto const& haf = layerCluster.hitsAndFractions(); + auto layerId = rhtools.getLayerWithOffset(haf[0].first); + + if (layerCluster.algo() == algo_number_ and layerId <= max_layerId_ and + layerId >= min_layerId_ and haf.size()<= max_cluster_size_ and + (haf.size() >= min_cluster_size_ or !(rhtools.isSilicon(haf[0].first)) ) ) { + filteredLayerClusters->emplace_back(cl); + } else { + layerClustersMask[cl.first] = 0.; + } + } + } + + private: + int algo_number_; + unsigned int min_cluster_size_; + unsigned int max_cluster_size_; + unsigned int min_layerId_; + unsigned int max_layerId_; + + }; +} // namespace ticl + +#endif diff --git a/RecoHGCal/TICL/plugins/filters.cc b/RecoHGCal/TICL/plugins/filters.cc index edbd834f4df90..54a89147dfb1d 100644 --- a/RecoHGCal/TICL/plugins/filters.cc +++ b/RecoHGCal/TICL/plugins/filters.cc @@ -7,9 +7,13 @@ #include "ClusterFilterByAlgo.h" #include "ClusterFilterByAlgoAndSize.h" #include "ClusterFilterBySize.h" +#include "ClusterFilterByAlgoAndSizeAndLayerRange.h" using namespace ticl; DEFINE_EDM_PLUGIN(ClusterFilterFactory, ClusterFilterByAlgo, "ClusterFilterByAlgo"); DEFINE_EDM_PLUGIN(ClusterFilterFactory, ClusterFilterByAlgoAndSize, "ClusterFilterByAlgoAndSize"); DEFINE_EDM_PLUGIN(ClusterFilterFactory, ClusterFilterBySize, "ClusterFilterBySize"); +DEFINE_EDM_PLUGIN(ClusterFilterFactory, + ClusterFilterByAlgoAndSizeAndLayerRange, + "ClusterFilterByAlgoAndSizeAndLayerRange"); From f39f3a981c20781d1c82f1b843c3053f61491179 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 30 Sep 2020 14:39:54 +0200 Subject: [PATCH 02/66] include an Electromagnetic TrackSeeded iteration as first iteration --- RecoHGCal/TICL/python/EMStep_cff.py | 4 +- RecoHGCal/TICL/python/HADStep_cff.py | 4 +- RecoHGCal/TICL/python/TrkEMStep_cff.py | 47 ++++++++++++++++++++++ RecoHGCal/TICL/python/TrkStep_cff.py | 2 + RecoHGCal/TICL/python/iterativeTICL_cff.py | 4 +- 5 files changed, 56 insertions(+), 5 deletions(-) create mode 100644 RecoHGCal/TICL/python/TrkEMStep_cff.py diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 0c3a5b5c497f3..7034a05bd6b86 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -11,7 +11,7 @@ clusterFilter = "ClusterFilterByAlgoAndSize", min_cluster_size = 2, # inclusive algo_number = 8, - LayerClustersInputMask = 'ticlTrackstersTrk', + LayerClustersInputMask = 'ticlTrackstersTrkEM', iteration_label = "EM" ) @@ -19,7 +19,7 @@ ticlTrackstersEM = _trackstersProducer.clone( filtered_mask = "filteredLayerClustersEM:EM", - original_mask = 'ticlTrackstersTrk', + original_mask = 'ticlTrackstersTrkEM', seeding_regions = "ticlSeedingGlobal", filter_on_categories = [0, 1], pid_threshold = 0.8, diff --git a/RecoHGCal/TICL/python/HADStep_cff.py b/RecoHGCal/TICL/python/HADStep_cff.py index 883f84412b22b..7bea230bd0c78 100644 --- a/RecoHGCal/TICL/python/HADStep_cff.py +++ b/RecoHGCal/TICL/python/HADStep_cff.py @@ -13,14 +13,14 @@ min_cluster_size = 2, # inclusive algo_number = 8, iteration_label = "HAD", - LayerClustersInputMask = "ticlTrackstersEM" + LayerClustersInputMask = "ticlTrackstersTrk" ) # CA - PATTERN RECOGNITION ticlTrackstersHAD = _trackstersProducer.clone( filtered_mask = "filteredLayerClustersHAD:HAD", - original_mask = 'ticlTrackstersEM', + original_mask = 'ticlTrackstersTrk', seeding_regions = "ticlSeedingGlobal", # For the moment we mask everything w/o requirements since we are last # filter_on_categories = [5], # filter neutral hadrons diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py new file mode 100644 index 0000000000000..1369fdc7dc118 --- /dev/null +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -0,0 +1,47 @@ +import FWCore.ParameterSet.Config as cms + +from RecoHGCal.TICL.TICLSeedingRegions_cff import ticlSeedingTrk +from RecoHGCal.TICL.trackstersProducer_cfi import trackstersProducer as _trackstersProducer +from RecoHGCal.TICL.filteredLayerClustersProducer_cfi import filteredLayerClustersProducer as _filteredLayerClustersProducer +from RecoHGCal.TICL.multiClustersFromTrackstersProducer_cfi import multiClustersFromTrackstersProducer as _multiClustersFromTrackstersProducer + +# CLUSTER FILTERING/MASKING + +filteredLayerClustersTrkEM = _filteredLayerClustersProducer.clone( + clusterFilter = "ClusterFilterByAlgoAndSizeAndLayerRange", + min_cluster_size = 2, # inclusive + min_layerId = 3, # inclusive + max_layerId = 30, # inclusive + algo_number = 8, + iteration_label = "TrkEM" +) + +# CA - PATTERN RECOGNITION + +ticlTrackstersTrkEM = _trackstersProducer.clone( + filtered_mask = cms.InputTag("filteredLayerClustersTrkEM", "TrkEM"), + seeding_regions = "ticlSeedingTrk", + filter_on_categories = [0, 1], + pid_threshold = 0.5, + max_out_in_hops = 4, + missing_layers = 1, + min_clusters_per_ntuplet = 10, + min_cos_theta = 0.978, # ~12 degrees + min_cos_pointing = 0.9, # ~25 degrees + max_delta_time = 3., + itername = "TrkEM", + algo_verbosity = 0, +) + + +# MULTICLUSTERS + +ticlMultiClustersFromTrackstersTrkEM = _multiClustersFromTrackstersProducer.clone( + Tracksters = "ticlTrackstersTrkEM" +) + +ticlTrkEMStepTask = cms.Task(ticlSeedingTrk + ,filteredLayerClustersTrkEM + ,ticlTrackstersTrkEM + ,ticlMultiClustersFromTrackstersTrkEM) + diff --git a/RecoHGCal/TICL/python/TrkStep_cff.py b/RecoHGCal/TICL/python/TrkStep_cff.py index ce33a5eab880f..67868d22dd341 100644 --- a/RecoHGCal/TICL/python/TrkStep_cff.py +++ b/RecoHGCal/TICL/python/TrkStep_cff.py @@ -11,6 +11,7 @@ filteredLayerClustersTrk = _filteredLayerClustersProducer.clone( clusterFilter = "ClusterFilterByAlgo", algo_number = 8, + LayerClustersInputMask = 'ticlTrackstersEM', iteration_label = "Trk" ) @@ -19,6 +20,7 @@ ticlTrackstersTrk = _trackstersProducer.clone( filtered_mask = "filteredLayerClustersTrk:Trk", seeding_regions = "ticlSeedingTrk", + original_mask = 'ticlTrackstersTrkEM', filter_on_categories = [2, 4], # filter muons and charged hadrons pid_threshold = 0.0, missing_layers = 3, diff --git a/RecoHGCal/TICL/python/iterativeTICL_cff.py b/RecoHGCal/TICL/python/iterativeTICL_cff.py index de31c0a2f4d3f..60f3788ed701b 100644 --- a/RecoHGCal/TICL/python/iterativeTICL_cff.py +++ b/RecoHGCal/TICL/python/iterativeTICL_cff.py @@ -1,6 +1,7 @@ import FWCore.ParameterSet.Config as cms from RecoHGCal.TICL.MIPStep_cff import * +from RecoHGCal.TICL.TrkEMStep_cff import * from RecoHGCal.TICL.TrkStep_cff import * from RecoHGCal.TICL.EMStep_cff import * from RecoHGCal.TICL.HADStep_cff import * @@ -28,8 +29,9 @@ iterTICLTask = cms.Task(ticlLayerTileTask ,ticlMIPStepTask - ,ticlTrkStepTask + ,ticlTrkEMStepTask ,ticlEMStepTask + ,ticlTrkStepTask ,ticlHADStepTask ,ticlTracksterMergeTask ,ticlPFTask From 2e437da2a6e468e60e9c9b341cb6da5a5155e421 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 30 Sep 2020 14:40:32 +0200 Subject: [PATCH 03/66] Add the range of layers to the filter parameters --- RecoHGCal/TICL/plugins/FilteredLayerClustersProducer.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/RecoHGCal/TICL/plugins/FilteredLayerClustersProducer.cc b/RecoHGCal/TICL/plugins/FilteredLayerClustersProducer.cc index bd26080cc959a..e3f0204cf47ba 100644 --- a/RecoHGCal/TICL/plugins/FilteredLayerClustersProducer.cc +++ b/RecoHGCal/TICL/plugins/FilteredLayerClustersProducer.cc @@ -65,6 +65,8 @@ void FilteredLayerClustersProducer::fillDescriptions(edm::ConfigurationDescripti desc.add("algo_number", 9); desc.add("min_cluster_size", 0); desc.add("max_cluster_size", 9999); + desc.add("min_layerId", 0); + desc.add("max_layerId", 9999); descriptions.add("filteredLayerClustersProducer", desc); } From eceff2738321508de8285c24eb39c687ccbbedb3 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 30 Sep 2020 14:41:01 +0200 Subject: [PATCH 04/66] Inject tracksters from the track seeded electromagnetic interation in the merged collection --- RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index 290eb66a605b5..2cb916a41a6dd 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -538,6 +538,7 @@ void TrackstersMergeProducer::printTrackstersDebug(const std::vector void TrackstersMergeProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { edm::ParameterSetDescription desc; + desc.add("tracksterstrkem", edm::InputTag("ticlTrackstersTrkEM")); desc.add("trackstersem", edm::InputTag("ticlTrackstersEM")); desc.add("tracksterstrk", edm::InputTag("ticlTrackstersTrk")); desc.add("trackstershad", edm::InputTag("ticlTrackstersHAD")); From d53fa4c7db359e5926694424adf4d33070a30c75 Mon Sep 17 00:00:00 2001 From: Felice Date: Fri, 2 Oct 2020 09:30:29 +0200 Subject: [PATCH 05/66] create empty multiclusters --- .../MultiClustersFromTrackstersProducer.cc | 22 ++++++++++--------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc b/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc index 195843d0015e8..938481fd887e0 100644 --- a/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc @@ -65,11 +65,12 @@ void MultiClustersFromTrackstersProducer::produce(edm::Event& evt, const edm::Ev std::for_each(std::begin(tracksters), std::end(tracksters), [&](auto const& trackster) { // Do not create a multicluster if the trackster has no layer clusters. // This could happen when a seed leads to no trackster and a dummy one is produced. + + std::array baricenter{{0., 0., 0.}}; + double total_weight = 0.; + reco::HGCalMultiCluster temp; + int counter = 0; if (!trackster.vertices().empty()) { - std::array baricenter{{0., 0., 0.}}; - double total_weight = 0.; - reco::HGCalMultiCluster temp; - int counter = 0; std::for_each(std::begin(trackster.vertices()), std::end(trackster.vertices()), [&](unsigned int idx) { temp.push_back(clusterPtrs[idx]); auto fraction = 1.f / trackster.vertex_multiplicity(counter++); @@ -86,13 +87,14 @@ void MultiClustersFromTrackstersProducer::produce(edm::Event& evt, const edm::Ev std::begin(baricenter), std::end(baricenter), std::begin(baricenter), [&total_weight](double val) -> double { return val / total_weight; }); - temp.setEnergy(total_weight); - temp.setCorrectedEnergy(total_weight); - temp.setPosition(math::XYZPoint(baricenter[0], baricenter[1], baricenter[2])); - temp.setAlgoId(reco::CaloCluster::hgcal_em); - temp.setTime(trackster.time(), trackster.timeError()); - multiclusters->push_back(temp); } + temp.setEnergy(total_weight); + temp.setCorrectedEnergy(total_weight); + temp.setPosition(math::XYZPoint(baricenter[0], baricenter[1], baricenter[2])); + temp.setAlgoId(reco::CaloCluster::hgcal_em); + temp.setTime(trackster.time(), trackster.timeError()); + multiclusters->push_back(temp); + }); evt.put(std::move(multiclusters)); From a394f3e177965e40bac43616202f988f927f0eda Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Mon, 5 Oct 2020 15:04:40 +0200 Subject: [PATCH 06/66] Add TrkEM tracksters into the Merge collection --- .../TICL/plugins/TrackstersMergeProducer.cc | 25 ++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index 2cb916a41a6dd..be103455d401b 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -40,6 +40,7 @@ class TrackstersMergeProducer : public edm::stream::EDProducer &, const char *label) const; void dumpTrackster(const Trackster &) const; + const edm::EDGetTokenT> tracksterstrkem_token_; const edm::EDGetTokenT> trackstersem_token_; const edm::EDGetTokenT> tracksterstrk_token_; const edm::EDGetTokenT> trackstershad_token_; @@ -72,7 +73,8 @@ class TrackstersMergeProducer : public edm::stream::EDProducer>(ps.getParameter("trackstersem"))), + : tracksterstrkem_token_(consumes>(ps.getParameter("tracksterstrkem"))), + trackstersem_token_(consumes>(ps.getParameter("trackstersem"))), tracksterstrk_token_(consumes>(ps.getParameter("tracksterstrk"))), trackstershad_token_(consumes>(ps.getParameter("trackstershad"))), seedingTrk_token_(consumes>(ps.getParameter("seedingTrk"))), @@ -153,6 +155,7 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es TICLTracksterTiles tracksterTile; std::vector usedTrackstersEM; + std::vector usedTrackstersTRKEM; std::vector usedTrackstersTRK; std::vector usedTrackstersHAD; std::vector usedSeeds; @@ -170,6 +173,11 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es const auto &trackstersEM = *trackstersem_h; usedTrackstersEM.resize(trackstersEM.size(), false); + edm::Handle> tracksterstrkem_h; + evt.getByToken(tracksterstrkem_token_, tracksterstrkem_h); + const auto &trackstersTRKEM = *tracksterstrkem_h; + usedTrackstersTRKEM.resize(trackstersTRKEM.size(), false); + edm::Handle> tracksterstrk_h; evt.getByToken(tracksterstrk_token_, tracksterstrk_h); const auto &trackstersTRK = *tracksterstrk_h; @@ -341,6 +349,21 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es tracksterHAD_idx++; } + auto tracksterTRKEM_idx = 0; + for (auto const &t : trackstersTRKEM) { + if (debug_) { + LogDebug("TrackstersMergeProducer") << " Considering trackster " << tracksterTRKEM_idx + << " as used: " << usedTrackstersTRKEM[tracksterTRKEM_idx] << std::endl; + } + if (!usedTrackstersTRKEM[tracksterTRKEM_idx]) { + LogDebug("TrackstersMergeProducer") + << " Creating a photon from EM Trackster with track energy " << t.raw_energy() << " and direction " + << t.eigenvectors(0).eta() << ", " << t.eigenvectors(0).phi() << std::endl; + result->push_back(t); + } + tracksterTRKEM_idx++; + } + assignPCAtoTracksters(*result, layerClusters, rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z()); energyRegressionAndID(layerClusters, *result); From 108265729386446d8a70ea1ef74ab9641b957202 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Mon, 5 Oct 2020 15:05:57 +0200 Subject: [PATCH 07/66] Select trackster-PID also on EM/Total trackster energy ratio --- RecoHGCal/TICL/plugins/TrackstersProducer.cc | 7 +++++-- RecoHGCal/TICL/python/EMStep_cff.py | 3 ++- RecoHGCal/TICL/python/TrkEMStep_cff.py | 1 + 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc index c4d228b130258..6c42d263e575a 100644 --- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc @@ -56,6 +56,7 @@ class TrackstersProducer : public edm::stream::EDProducer> seeding_regions_token_; const std::vector filter_on_categories_; const double pid_threshold_; + const double energy_em_over_total_threshold_; const std::string itername_; }; DEFINE_FWK_MODULE(TrackstersProducer); @@ -93,6 +94,7 @@ TrackstersProducer::TrackstersProducer(const edm::ParameterSet& ps, const Tracks consumes>(ps.getParameter("seeding_regions"))), filter_on_categories_(ps.getParameter>("filter_on_categories")), pid_threshold_(ps.getParameter("pid_threshold")), + energy_em_over_total_threshold_(ps.getParameter("energy_em_over_total_threshold")), itername_(ps.getParameter("itername")) { if (doNose_) { layer_clusters_tiles_hfnose_token_ = @@ -117,6 +119,7 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri desc.add("seeding_regions", edm::InputTag("ticlSeedingRegionProducer")); desc.add>("filter_on_categories", {0}); desc.add("pid_threshold", 0.); // make default such that no filtering is applied + desc.add("energy_em_over_total_threshold", -1.); // make default such that no filtering is applied desc.add("algo_verbosity", 0); desc.add("min_cos_theta", 0.915); desc.add("min_cos_pointing", -1.); @@ -179,7 +182,7 @@ void TrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) { std::copy( std::begin(original_layerclusters_mask), std::end(original_layerclusters_mask), std::back_inserter(*output_mask)); - // Filter results based on PID criteria. + // Filter results based on PID criteria or EM/Total energy ratio. // We want to **keep** tracksters whose cumulative // probability summed up over the selected categories // is greater than the chosen threshold. Therefore @@ -190,7 +193,7 @@ void TrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) { for (auto index : filter_on_categories_) { cumulative_prob += t.id_probabilities(index); } - return cumulative_prob <= pid_threshold_; + return (cumulative_prob <= pid_threshold_) && (t.raw_em_energy()/t.raw_energy() < energy_em_over_total_threshold_); }; // Actually filter results and shrink size to fit diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 7034a05bd6b86..49546e6207ff0 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -22,7 +22,8 @@ original_mask = 'ticlTrackstersTrkEM', seeding_regions = "ticlSeedingGlobal", filter_on_categories = [0, 1], - pid_threshold = 0.8, + pid_threshold = 0.5, + energy_em_over_total_threshold = 0.85, max_out_in_hops = 4, missing_layers = 1, min_clusters_per_ntuplet = 10, diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index 1369fdc7dc118..05aaeac6298e0 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -23,6 +23,7 @@ seeding_regions = "ticlSeedingTrk", filter_on_categories = [0, 1], pid_threshold = 0.5, + energy_em_over_total_threshold = 0.85, max_out_in_hops = 4, missing_layers = 1, min_clusters_per_ntuplet = 10, From 0101b451732cd27b5482ade59216b57b825e8342 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Mon, 5 Oct 2020 16:24:31 +0200 Subject: [PATCH 08/66] Fix merge conflict. --- RecoHGCal/TICL/python/TrkStep_cff.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/RecoHGCal/TICL/python/TrkStep_cff.py b/RecoHGCal/TICL/python/TrkStep_cff.py index 67868d22dd341..3b5d527405d50 100644 --- a/RecoHGCal/TICL/python/TrkStep_cff.py +++ b/RecoHGCal/TICL/python/TrkStep_cff.py @@ -25,6 +25,8 @@ pid_threshold = 0.0, missing_layers = 3, min_clusters_per_ntuplet = 10, +# min_cos_theta = 0.978, # same as em +# min_cos_pointing = 0.9, # same as em min_cos_theta = 0.866, # ~30 degrees min_cos_pointing = 0.798, # ~ 37 degrees max_delta_time = -1., @@ -45,3 +47,4 @@ ,ticlTrackstersTrk ,ticlMultiClustersFromTrackstersTrk) + From b039481d3d22b13f9d3ce28d8b08966209e77e2d Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 7 Oct 2020 16:00:40 +0200 Subject: [PATCH 09/66] Fix input collection for TRK iteration to be the output of EM iteration --- RecoHGCal/TICL/python/TrkStep_cff.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/python/TrkStep_cff.py b/RecoHGCal/TICL/python/TrkStep_cff.py index 3b5d527405d50..cd9dcd9b7d0ff 100644 --- a/RecoHGCal/TICL/python/TrkStep_cff.py +++ b/RecoHGCal/TICL/python/TrkStep_cff.py @@ -20,7 +20,7 @@ ticlTrackstersTrk = _trackstersProducer.clone( filtered_mask = "filteredLayerClustersTrk:Trk", seeding_regions = "ticlSeedingTrk", - original_mask = 'ticlTrackstersTrkEM', + original_mask = 'ticlTrackstersEM', filter_on_categories = [2, 4], # filter muons and charged hadrons pid_threshold = 0.0, missing_layers = 3, From 23b03976cd9dbc81dc18fd8a7ba934f06a19bc52 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 7 Oct 2020 16:22:16 +0200 Subject: [PATCH 10/66] Sorting seeding regions by track momentum --- RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc b/RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc index 45ecb51d13907..55b9cf82da219 100644 --- a/RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc +++ b/RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc @@ -59,6 +59,10 @@ void SeedingRegionByTracks::makeRegions(const edm::Event &ev, result.emplace_back(tsos.globalPosition(), tsos.globalMomentum(), iSide, i, trkId); } } + // sorting seeding region by descending momentum + std::sort(result.begin(), result.end(), [](const TICLSeedingRegion& a, const TICLSeedingRegion& b) { + return a.directionAtOrigin.mag2() > b.directionAtOrigin.mag2(); + }); } void SeedingRegionByTracks::buildFirstLayers() { From d69212b42535077ddd2e341506648005b54896f7 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 7 Oct 2020 16:31:06 +0200 Subject: [PATCH 11/66] Use only clusters with minimum size of 3 --- RecoHGCal/TICL/python/EMStep_cff.py | 2 +- RecoHGCal/TICL/python/HADStep_cff.py | 2 +- RecoHGCal/TICL/python/TrkEMStep_cff.py | 2 +- RecoHGCal/TICL/python/TrkStep_cff.py | 3 ++- 4 files changed, 5 insertions(+), 4 deletions(-) diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 49546e6207ff0..3c1fa84a67bae 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -9,7 +9,7 @@ filteredLayerClustersEM = _filteredLayerClustersProducer.clone( clusterFilter = "ClusterFilterByAlgoAndSize", - min_cluster_size = 2, # inclusive + min_cluster_size = 3, # inclusive algo_number = 8, LayerClustersInputMask = 'ticlTrackstersTrkEM', iteration_label = "EM" diff --git a/RecoHGCal/TICL/python/HADStep_cff.py b/RecoHGCal/TICL/python/HADStep_cff.py index 7bea230bd0c78..d631f7196d137 100644 --- a/RecoHGCal/TICL/python/HADStep_cff.py +++ b/RecoHGCal/TICL/python/HADStep_cff.py @@ -10,7 +10,7 @@ filteredLayerClustersHAD = _filteredLayerClustersProducer.clone( clusterFilter = "ClusterFilterByAlgoAndSize", - min_cluster_size = 2, # inclusive + min_cluster_size = 3, # inclusive algo_number = 8, iteration_label = "HAD", LayerClustersInputMask = "ticlTrackstersTrk" diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index 05aaeac6298e0..05f71337f305c 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -9,7 +9,7 @@ filteredLayerClustersTrkEM = _filteredLayerClustersProducer.clone( clusterFilter = "ClusterFilterByAlgoAndSizeAndLayerRange", - min_cluster_size = 2, # inclusive + min_cluster_size = 3, # inclusive min_layerId = 3, # inclusive max_layerId = 30, # inclusive algo_number = 8, diff --git a/RecoHGCal/TICL/python/TrkStep_cff.py b/RecoHGCal/TICL/python/TrkStep_cff.py index cd9dcd9b7d0ff..e35a692ff194b 100644 --- a/RecoHGCal/TICL/python/TrkStep_cff.py +++ b/RecoHGCal/TICL/python/TrkStep_cff.py @@ -9,7 +9,8 @@ # CLUSTER FILTERING/MASKING filteredLayerClustersTrk = _filteredLayerClustersProducer.clone( - clusterFilter = "ClusterFilterByAlgo", + clusterFilter = "ClusterFilterByAlgoAndSize", + min_cluster_size = 3, # inclusive algo_number = 8, LayerClustersInputMask = 'ticlTrackstersEM', iteration_label = "Trk" From cb843b5c69e5e5f8951f1fe6d2c21d037e0ca60b Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 7 Oct 2020 16:31:56 +0200 Subject: [PATCH 12/66] Filter only with a maxLayerId in TRKEM iteration. --- RecoHGCal/TICL/python/TrkEMStep_cff.py | 1 - 1 file changed, 1 deletion(-) diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index 05f71337f305c..9ada6832ebd3a 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -10,7 +10,6 @@ filteredLayerClustersTrkEM = _filteredLayerClustersProducer.clone( clusterFilter = "ClusterFilterByAlgoAndSizeAndLayerRange", min_cluster_size = 3, # inclusive - min_layerId = 3, # inclusive max_layerId = 30, # inclusive algo_number = 8, iteration_label = "TrkEM" From 8790b01102e0f063c461f61da6af4b3082b9fac0 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 7 Oct 2020 17:15:07 +0200 Subject: [PATCH 13/66] Remove MIP iteration --- RecoHGCal/TICL/python/iterativeTICL_cff.py | 1 - 1 file changed, 1 deletion(-) diff --git a/RecoHGCal/TICL/python/iterativeTICL_cff.py b/RecoHGCal/TICL/python/iterativeTICL_cff.py index 60f3788ed701b..b6a9b8048629a 100644 --- a/RecoHGCal/TICL/python/iterativeTICL_cff.py +++ b/RecoHGCal/TICL/python/iterativeTICL_cff.py @@ -28,7 +28,6 @@ ticlPFTask = cms.Task(ticlCandidateFromTracksters, pfTICL) iterTICLTask = cms.Task(ticlLayerTileTask - ,ticlMIPStepTask ,ticlTrkEMStepTask ,ticlEMStepTask ,ticlTrkStepTask From 41636a6382de551baaf7de2b77a8f6620d0bdd8f Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 7 Oct 2020 18:41:59 +0200 Subject: [PATCH 14/66] clang format --- .../plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h | 7 +++---- .../TICL/plugins/MultiClustersFromTrackstersProducer.cc | 3 +-- RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc | 8 ++++---- RecoHGCal/TICL/plugins/TrackstersProducer.cc | 5 +++-- 4 files changed, 11 insertions(+), 12 deletions(-) diff --git a/RecoHGCal/TICL/plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h b/RecoHGCal/TICL/plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h index f1c901ade9404..5f1d9388c3fe4 100644 --- a/RecoHGCal/TICL/plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h +++ b/RecoHGCal/TICL/plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h @@ -33,9 +33,9 @@ namespace ticl { auto const& haf = layerCluster.hitsAndFractions(); auto layerId = rhtools.getLayerWithOffset(haf[0].first); - if (layerCluster.algo() == algo_number_ and layerId <= max_layerId_ and - layerId >= min_layerId_ and haf.size()<= max_cluster_size_ and - (haf.size() >= min_cluster_size_ or !(rhtools.isSilicon(haf[0].first)) ) ) { + if (layerCluster.algo() == algo_number_ and layerId <= max_layerId_ and layerId >= min_layerId_ and + haf.size() <= max_cluster_size_ and + (haf.size() >= min_cluster_size_ or !(rhtools.isSilicon(haf[0].first)))) { filteredLayerClusters->emplace_back(cl); } else { layerClustersMask[cl.first] = 0.; @@ -49,7 +49,6 @@ namespace ticl { unsigned int max_cluster_size_; unsigned int min_layerId_; unsigned int max_layerId_; - }; } // namespace ticl diff --git a/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc b/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc index 938481fd887e0..c438a92be6ef9 100644 --- a/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc @@ -65,7 +65,7 @@ void MultiClustersFromTrackstersProducer::produce(edm::Event& evt, const edm::Ev std::for_each(std::begin(tracksters), std::end(tracksters), [&](auto const& trackster) { // Do not create a multicluster if the trackster has no layer clusters. // This could happen when a seed leads to no trackster and a dummy one is produced. - + std::array baricenter{{0., 0., 0.}}; double total_weight = 0.; reco::HGCalMultiCluster temp; @@ -94,7 +94,6 @@ void MultiClustersFromTrackstersProducer::produce(edm::Event& evt, const edm::Ev temp.setAlgoId(reco::CaloCluster::hgcal_em); temp.setTime(trackster.time(), trackster.timeError()); multiclusters->push_back(temp); - }); evt.put(std::move(multiclusters)); diff --git a/RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc b/RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc index 55b9cf82da219..a907a69c750a6 100644 --- a/RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc +++ b/RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc @@ -59,10 +59,10 @@ void SeedingRegionByTracks::makeRegions(const edm::Event &ev, result.emplace_back(tsos.globalPosition(), tsos.globalMomentum(), iSide, i, trkId); } } - // sorting seeding region by descending momentum - std::sort(result.begin(), result.end(), [](const TICLSeedingRegion& a, const TICLSeedingRegion& b) { - return a.directionAtOrigin.mag2() > b.directionAtOrigin.mag2(); - }); + // sorting seeding region by descending momentum + std::sort(result.begin(), result.end(), [](const TICLSeedingRegion &a, const TICLSeedingRegion &b) { + return a.directionAtOrigin.mag2() > b.directionAtOrigin.mag2(); + }); } void SeedingRegionByTracks::buildFirstLayers() { diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc index 6c42d263e575a..aa47c332d187e 100644 --- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc @@ -118,7 +118,7 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri desc.add("layer_clusters_hfnose_tiles", edm::InputTag("ticlLayerTileHFNose")); desc.add("seeding_regions", edm::InputTag("ticlSeedingRegionProducer")); desc.add>("filter_on_categories", {0}); - desc.add("pid_threshold", 0.); // make default such that no filtering is applied + desc.add("pid_threshold", 0.); // make default such that no filtering is applied desc.add("energy_em_over_total_threshold", -1.); // make default such that no filtering is applied desc.add("algo_verbosity", 0); desc.add("min_cos_theta", 0.915); @@ -193,7 +193,8 @@ void TrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) { for (auto index : filter_on_categories_) { cumulative_prob += t.id_probabilities(index); } - return (cumulative_prob <= pid_threshold_) && (t.raw_em_energy()/t.raw_energy() < energy_em_over_total_threshold_); + return (cumulative_prob <= pid_threshold_) && + (t.raw_em_energy() / t.raw_energy() < energy_em_over_total_threshold_); }; // Actually filter results and shrink size to fit From 9eb885d6a11f175fbf27e90cd5a463ea42bd06c6 Mon Sep 17 00:00:00 2001 From: Felice Date: Thu, 8 Oct 2020 15:30:57 +0200 Subject: [PATCH 15/66] Limit Global EM iteration to layer 30 --- RecoHGCal/TICL/python/EMStep_cff.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 3c1fa84a67bae..47550f184bce3 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -8,8 +8,9 @@ # CLUSTER FILTERING/MASKING filteredLayerClustersEM = _filteredLayerClustersProducer.clone( - clusterFilter = "ClusterFilterByAlgoAndSize", + clusterFilter = "ClusterFilterByAlgoAndSizeAndLayerRange", min_cluster_size = 3, # inclusive + max_layerId = 30, # inclusive algo_number = 8, LayerClustersInputMask = 'ticlTrackstersTrkEM', iteration_label = "EM" From dbbaaae6f74b174842e42bd6561c26cb0000d5ca Mon Sep 17 00:00:00 2001 From: Felice Date: Thu, 8 Oct 2020 15:31:57 +0200 Subject: [PATCH 16/66] Limit track seeded iterations to tracks with at most 2 missing outer hits in the Tracker --- RecoHGCal/TICL/plugins/TICLSeedingRegionProducer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/TICLSeedingRegionProducer.cc b/RecoHGCal/TICL/plugins/TICLSeedingRegionProducer.cc index 152c7b7898eb4..e42abbd1c43c7 100644 --- a/RecoHGCal/TICL/plugins/TICLSeedingRegionProducer.cc +++ b/RecoHGCal/TICL/plugins/TICLSeedingRegionProducer.cc @@ -59,7 +59,7 @@ void TICLSeedingRegionProducer::fillDescriptions(edm::ConfigurationDescriptions& desc.add("tracks", edm::InputTag("generalTracks")); desc.add("cutTk", "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && " - "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 10"); + "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 3"); desc.add("propagator", "PropagatorWithMaterial"); desc.add("algoId", 1); descriptions.add("ticlSeedingRegionProducer", desc); From b576c9eb867efe27e1effe73fda4aff6ee00b647 Mon Sep 17 00:00:00 2001 From: Felice Date: Thu, 8 Oct 2020 18:43:44 +0200 Subject: [PATCH 17/66] Add TrkEM products to EventContent --- RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py index 3a16bea88e3b3..5b0ecc4d4b9c9 100644 --- a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py +++ b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py @@ -7,6 +7,7 @@ 'keep *_ticlMultiClustersFromTrackstersEM_*_*', 'keep *_ticlMultiClustersFromTrackstersHAD_*_*', 'keep *_ticlMultiClustersFromTrackstersTrk_*_*', + 'keep *_ticlMultiClustersFromTrackstersTrkEM_*_*', 'keep *_ticlMultiClustersFromTrackstersMIP_*_*', 'keep *_ticlMultiClustersFromTrackstersMerge_*_*', ) @@ -15,6 +16,7 @@ #RECO content TICL_RECO = cms.PSet( outputCommands = cms.untracked.vstring( + 'keep *_ticlTrackstersTrkEM_*_*', 'keep *_ticlTrackstersEM_*_*', 'keep *_ticlTrackstersHAD_*_*', 'keep *_ticlTrackstersTrk_*_*', From 72725e3fa46ecf67c9971e80336600b4976ed026 Mon Sep 17 00:00:00 2001 From: Felice Date: Thu, 8 Oct 2020 18:45:12 +0200 Subject: [PATCH 18/66] Remove missing layers cut for EM based iterations --- RecoHGCal/TICL/python/EMStep_cff.py | 2 +- RecoHGCal/TICL/python/TrkEMStep_cff.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 47550f184bce3..391a2280b6d14 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -26,7 +26,7 @@ pid_threshold = 0.5, energy_em_over_total_threshold = 0.85, max_out_in_hops = 4, - missing_layers = 1, + missing_layers = 0, min_clusters_per_ntuplet = 10, min_cos_theta = 0.978, # ~12 degrees min_cos_pointing = 0.9, # ~25 degrees diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index 9ada6832ebd3a..e7cda493daf87 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -24,7 +24,7 @@ pid_threshold = 0.5, energy_em_over_total_threshold = 0.85, max_out_in_hops = 4, - missing_layers = 1, + missing_layers = 0, min_clusters_per_ntuplet = 10, min_cos_theta = 0.978, # ~12 degrees min_cos_pointing = 0.9, # ~25 degrees From 918c3a4cb69852c2e668ad87873b1f2b1837e7bf Mon Sep 17 00:00:00 2001 From: Felice Date: Thu, 8 Oct 2020 19:44:56 +0200 Subject: [PATCH 19/66] Introducing a shower_start_max_layer cut --- .../TICL/plugins/PatternRecognitionbyCA.cc | 41 +++++++++++-------- .../TICL/plugins/PatternRecognitionbyCA.h | 2 + RecoHGCal/TICL/plugins/TrackstersProducer.cc | 4 ++ RecoHGCal/TICL/python/EMStep_cff.py | 1 + RecoHGCal/TICL/python/TrkEMStep_cff.py | 1 + 5 files changed, 32 insertions(+), 17 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index dcc0c039d15db..c740ab342e285 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -28,6 +28,7 @@ PatternRecognitionbyCA::PatternRecognitionbyCA(const edm::ParameterSet &c min_cos_pointing_(conf.getParameter("min_cos_pointing")), etaLimitIncreaseWindow_(conf.getParameter("etaLimitIncreaseWindow")), missing_layers_(conf.getParameter("missing_layers")), + shower_start_max_layer_(conf.getParameter("shower_start_max_layer")), min_clusters_per_ntuplet_(conf.getParameter("min_clusters_per_ntuplet")), max_delta_time_(conf.getParameter("max_delta_time")), eidInputName_(conf.getParameter("eid_input_name")), @@ -137,29 +138,35 @@ void PatternRecognitionbyCA::makeTracksters( << input.layerClusters[outerCluster].z() << " " << tracksterId << std::endl; } } + unsigned showerMinLayerId = 99999; for (auto const i : effective_cluster_idx) { layer_cluster_usage[i]++; + auto const& haf = input.layerClusters[i].hitsAndFractions(); + showerMinLayerId = std::min(rhtools_.getLayerWithOffset(haf[0].first),showerMinLayerId); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) LogDebug("HGCPatternRecoByCA") << "LayerID: " << i << " count: " << (int)layer_cluster_usage[i] << std::endl; } - // Put back indices, in the form of a Trackster, into the results vector - Trackster tmp; - tmp.vertices().reserve(effective_cluster_idx.size()); - tmp.vertex_multiplicity().resize(effective_cluster_idx.size(), 0); - //regions and seedIndices can have different size - //if a seeding region does not lead to any trackster - tmp.setSeed(input.regions[0].collectionID, seedIndices[tracksterId]); - if (isRegionalIter) { - seedToTracksterAssociation[tmp.seedIndex()].push_back(tracksterId); - } - std::pair timeTrackster(-99., -1.); - hgcalsimclustertime::ComputeClusterTime timeEstimator; - timeTrackster = timeEstimator.fixSizeHighestDensity(times, timeErrors); - tmp.setTimeAndError(timeTrackster.first, timeTrackster.second); - std::copy(std::begin(effective_cluster_idx), std::end(effective_cluster_idx), std::back_inserter(tmp.vertices())); - result.push_back(tmp); - tracksterId++; + if(showerMinLayerId <= shower_start_max_layer_){ + // Put back indices, in the form of a Trackster, into the results vector + Trackster tmp; + tmp.vertices().reserve(effective_cluster_idx.size()); + tmp.vertex_multiplicity().resize(effective_cluster_idx.size(), 0); + //regions and seedIndices can have different size + //if a seeding region does not lead to any trackster + tmp.setSeed(input.regions[0].collectionID, seedIndices[tracksterId]); + if (isRegionalIter) { + seedToTracksterAssociation[tmp.seedIndex()].push_back(tracksterId); + } + + std::pair timeTrackster(-99., -1.); + hgcalsimclustertime::ComputeClusterTime timeEstimator; + timeTrackster = timeEstimator.fixSizeHighestDensity(times, timeErrors); + tmp.setTimeAndError(timeTrackster.first, timeTrackster.second); + std::copy(std::begin(effective_cluster_idx), std::end(effective_cluster_idx), std::back_inserter(tmp.vertices())); + result.push_back(tmp); + tracksterId++; + } } for (auto &trackster : result) { diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h index f6b917cdb2899..0ed31dac06f44 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h @@ -39,6 +39,8 @@ namespace ticl { const float min_cos_pointing_; const float etaLimitIncreaseWindow_; const int missing_layers_; + const unsigned int shower_start_max_layer_; + const int min_clusters_per_ntuplet_; const float max_delta_time_; const std::string eidInputName_; diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc index aa47c332d187e..16462d76612e7 100644 --- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc @@ -57,6 +57,8 @@ class TrackstersProducer : public edm::stream::EDProducer filter_on_categories_; const double pid_threshold_; const double energy_em_over_total_threshold_; + const int shower_start_max_layer_; + const std::string itername_; }; DEFINE_FWK_MODULE(TrackstersProducer); @@ -95,6 +97,7 @@ TrackstersProducer::TrackstersProducer(const edm::ParameterSet& ps, const Tracks filter_on_categories_(ps.getParameter>("filter_on_categories")), pid_threshold_(ps.getParameter("pid_threshold")), energy_em_over_total_threshold_(ps.getParameter("energy_em_over_total_threshold")), + shower_start_max_layer_(ps.getParameter("shower_start_max_layer")), itername_(ps.getParameter("itername")) { if (doNose_) { layer_clusters_tiles_hfnose_token_ = @@ -120,6 +123,7 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri desc.add>("filter_on_categories", {0}); desc.add("pid_threshold", 0.); // make default such that no filtering is applied desc.add("energy_em_over_total_threshold", -1.); // make default such that no filtering is applied + desc.add("shower_start_max_layer", 9999); // make default such that no filtering is applied desc.add("algo_verbosity", 0); desc.add("min_cos_theta", 0.915); desc.add("min_cos_pointing", -1.); diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 391a2280b6d14..5b19a704bd877 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -25,6 +25,7 @@ filter_on_categories = [0, 1], pid_threshold = 0.5, energy_em_over_total_threshold = 0.85, + shower_start_max_layer = 3, #inclusive max_out_in_hops = 4, missing_layers = 0, min_clusters_per_ntuplet = 10, diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index e7cda493daf87..fcdfef7063199 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -23,6 +23,7 @@ filter_on_categories = [0, 1], pid_threshold = 0.5, energy_em_over_total_threshold = 0.85, + shower_start_max_layer = 3, #inclusive max_out_in_hops = 4, missing_layers = 0, min_clusters_per_ntuplet = 10, From 152492522adfee961b4679acc0fd22c4a5f24332 Mon Sep 17 00:00:00 2001 From: Felice Date: Fri, 9 Oct 2020 10:33:27 +0200 Subject: [PATCH 20/66] raising shower_start_max_layer from 3 to 7 in EM iterations --- RecoHGCal/TICL/python/EMStep_cff.py | 2 +- RecoHGCal/TICL/python/TrkEMStep_cff.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 5b19a704bd877..e8fe29101959f 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -25,7 +25,7 @@ filter_on_categories = [0, 1], pid_threshold = 0.5, energy_em_over_total_threshold = 0.85, - shower_start_max_layer = 3, #inclusive + shower_start_max_layer = 7, #inclusive max_out_in_hops = 4, missing_layers = 0, min_clusters_per_ntuplet = 10, diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index fcdfef7063199..ec73608400394 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -23,7 +23,7 @@ filter_on_categories = [0, 1], pid_threshold = 0.5, energy_em_over_total_threshold = 0.85, - shower_start_max_layer = 3, #inclusive + shower_start_max_layer = 7, #inclusive max_out_in_hops = 4, missing_layers = 0, min_clusters_per_ntuplet = 10, From cdb0e6a98e4c43665eaa29282cce2aed725636e9 Mon Sep 17 00:00:00 2001 From: Felice Date: Fri, 9 Oct 2020 11:25:26 +0200 Subject: [PATCH 21/66] Count usage of layer clusters, only for those tracksters passing the selection --- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index c740ab342e285..653edd5710858 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -140,14 +140,16 @@ void PatternRecognitionbyCA::makeTracksters( } unsigned showerMinLayerId = 99999; for (auto const i : effective_cluster_idx) { - layer_cluster_usage[i]++; auto const& haf = input.layerClusters[i].hitsAndFractions(); showerMinLayerId = std::min(rhtools_.getLayerWithOffset(haf[0].first),showerMinLayerId); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) - LogDebug("HGCPatternRecoByCA") << "LayerID: " << i << " count: " << (int)layer_cluster_usage[i] << std::endl; } if(showerMinLayerId <= shower_start_max_layer_){ + for (auto const i : effective_cluster_idx) { + layer_cluster_usage[i]++; + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) + LogDebug("HGCPatternRecoByCA") << "LayerID: " << i << " count: " << (int)layer_cluster_usage[i] << std::endl; + } // Put back indices, in the form of a Trackster, into the results vector Trackster tmp; tmp.vertices().reserve(effective_cluster_idx.size()); From 936c10041aff5f4844e659fbcc7a071b2463fcf2 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 7 Oct 2020 18:41:41 +0200 Subject: [PATCH 22/66] Protect from empty multiclusters --- .../plugins/PFClusterFromHGCalMultiCluster.cc | 3 + .../src/HGVHistoProducerAlgo.cc | 377 +++++++++--------- 2 files changed, 194 insertions(+), 186 deletions(-) diff --git a/RecoParticleFlow/PFClusterProducer/plugins/PFClusterFromHGCalMultiCluster.cc b/RecoParticleFlow/PFClusterProducer/plugins/PFClusterFromHGCalMultiCluster.cc index 33d6362ce04f4..dd4352de9719e 100644 --- a/RecoParticleFlow/PFClusterProducer/plugins/PFClusterFromHGCalMultiCluster.cc +++ b/RecoParticleFlow/PFClusterProducer/plugins/PFClusterFromHGCalMultiCluster.cc @@ -43,6 +43,9 @@ void PFClusterFromHGCalMultiCluster::buildClusters(const edm::Handle CPEnergyInMCL; + int maxCPId_byNumberOfHits = -1; + unsigned int maxCPNumberOfHitsInMCL = 0; + int maxCPId_byEnergy = -1; + float maxEnergySharedMCLandCP = 0.f; + float energyFractionOfMCLinCP = 0.f; + float energyFractionOfCPinMCL = 0.f; + + //In case of matched rechit-simhit, so matched + //caloparticle-layercluster-multicluster, he counts and saves the number of + //rechits related to the maximum energy CaloParticle out of all + //CaloParticles related to that layer cluster and multicluster. + + std::unordered_map occurrencesCPinMCL; + unsigned int numberOfNoiseHitsInMCL = 0; + unsigned int numberOfHaloHitsInMCL = 0; + unsigned int numberOfHitsInMCL = 0; + + //number of hits related to that cluster. + unsigned int numberOfHitsInLC = hits_and_fractions.size(); + numberOfHitsInMCL += numberOfHitsInLC; + std::unordered_map CPEnergyInLC; + + //hitsToCaloParticleId is a vector of ints, one for each rechit of the + //layer cluster under study. If negative, there is no simhit from any CaloParticle related. + //If positive, at least one CaloParticle has been found with matched simhit. + //In more detail: + // 1. hitsToCaloParticleId[hitId] = -3 + // TN: These represent Halo Cells(N) that have not been + // assigned to any CaloParticle (hence the T). + // 2. hitsToCaloParticleId[hitId] = -2 + // FN: There represent Halo Cells(N) that have been assigned + // to a CaloParticle (hence the F, since those should have not been marked as halo) + // 3. hitsToCaloParticleId[hitId] = -1 + // FP: These represent Real Cells(P) that have not been + // assigned to any CaloParticle (hence the F, since these are fakes) + // 4. hitsToCaloParticleId[hitId] >= 0 + // TP There represent Real Cells(P) that have been assigned + // to a CaloParticle (hence the T) + + std::vector hitsToCaloParticleId(numberOfHitsInLC); + //det id of the first hit just to make the lcLayerId variable + //which maps the layers in -z: 0->51 and in +z: 52->103 + const auto firstHitDetId = hits_and_fractions[0].first; + int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId) + + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1; - std::unordered_map CPEnergyInMCL; - int maxCPId_byNumberOfHits = -1; - unsigned int maxCPNumberOfHitsInMCL = 0; - int maxCPId_byEnergy = -1; - float maxEnergySharedMCLandCP = 0.f; - float energyFractionOfMCLinCP = 0.f; - float energyFractionOfCPinMCL = 0.f; - - //In case of matched rechit-simhit, so matched - //caloparticle-layercluster-multicluster, he counts and saves the number of - //rechits related to the maximum energy CaloParticle out of all - //CaloParticles related to that layer cluster and multicluster. - - std::unordered_map occurrencesCPinMCL; - unsigned int numberOfNoiseHitsInMCL = 0; - unsigned int numberOfHaloHitsInMCL = 0; - unsigned int numberOfHitsInMCL = 0; - - //number of hits related to that cluster. - unsigned int numberOfHitsInLC = hits_and_fractions.size(); - numberOfHitsInMCL += numberOfHitsInLC; - std::unordered_map CPEnergyInLC; - - //hitsToCaloParticleId is a vector of ints, one for each rechit of the - //layer cluster under study. If negative, there is no simhit from any CaloParticle related. - //If positive, at least one CaloParticle has been found with matched simhit. - //In more detail: - // 1. hitsToCaloParticleId[hitId] = -3 - // TN: These represent Halo Cells(N) that have not been - // assigned to any CaloParticle (hence the T). - // 2. hitsToCaloParticleId[hitId] = -2 - // FN: There represent Halo Cells(N) that have been assigned - // to a CaloParticle (hence the F, since those should have not been marked as halo) - // 3. hitsToCaloParticleId[hitId] = -1 - // FP: These represent Real Cells(P) that have not been - // assigned to any CaloParticle (hence the F, since these are fakes) - // 4. hitsToCaloParticleId[hitId] >= 0 - // TP There represent Real Cells(P) that have been assigned - // to a CaloParticle (hence the T) - - std::vector hitsToCaloParticleId(numberOfHitsInLC); - //det id of the first hit just to make the lcLayerId variable - //which maps the layers in -z: 0->51 and in +z: 52->103 - const auto firstHitDetId = hits_and_fractions[0].first; - int lcLayerId = - recHitTools_->getLayerWithOffset(firstHitDetId) + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1; - - //Loop though the hits of the layer cluster under study - for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) { - DetId rh_detid = hits_and_fractions[hitId].first; - auto rhFraction = hits_and_fractions[hitId].second; + //Loop though the hits of the layer cluster under study + for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) { + DetId rh_detid = hits_and_fractions[hitId].first; + auto rhFraction = hits_and_fractions[hitId].second; - //Since the hit is belonging to the layer cluster, it must also be in the rechits map. - std::unordered_map::const_iterator itcheck = hitMap.find(rh_detid); - const HGCRecHit* hit = itcheck->second; + //Since the hit is belonging to the layer cluster, it must also be in the rechits map. + std::unordered_map::const_iterator itcheck = hitMap.find(rh_detid); + const HGCRecHit* hit = itcheck->second; - //Make a map that will connect a detid (that belongs to a rechit of the layer cluster under study, - //no need to save others) with: - //1. the layer clusters that have rechits in that detid - //2. the fraction of the rechit of each layer cluster that contributes to that detid. - //So, something like: - //detid: (layer cluster 1, hit fraction) , (layer cluster 2, hit fraction), (layer cluster 3, hit fraction) ... - //here comparing with the calo particle map above the - auto hit_find_in_LC = detIdToMultiClusterId_Map.find(rh_detid); - if (hit_find_in_LC == detIdToMultiClusterId_Map.end()) { - detIdToMultiClusterId_Map[rh_detid] = std::vector(); - } - detIdToMultiClusterId_Map[rh_detid].emplace_back( - HGVHistoProducerAlgo::detIdInfoInMultiCluster{mclId, mclId, rhFraction}); + //Make a map that will connect a detid (that belongs to a rechit of the layer cluster under study, + //no need to save others) with: + //1. the layer clusters that have rechits in that detid + //2. the fraction of the rechit of each layer cluster that contributes to that detid. + //So, something like: + //detid: (layer cluster 1, hit fraction) , (layer cluster 2, hit fraction), (layer cluster 3, hit fraction) ... + //here comparing with the calo particle map above the + auto hit_find_in_LC = detIdToMultiClusterId_Map.find(rh_detid); + if (hit_find_in_LC == detIdToMultiClusterId_Map.end()) { + detIdToMultiClusterId_Map[rh_detid] = std::vector(); + } + detIdToMultiClusterId_Map[rh_detid].emplace_back( + HGVHistoProducerAlgo::detIdInfoInMultiCluster{mclId, mclId, rhFraction}); - //Check whether the rechit of the layer cluster under study has a sim hit in the same cell. - auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid); + //Check whether the rechit of the layer cluster under study has a sim hit in the same cell. + auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid); - // if the fraction is zero or the hit does not belong to any calo - // particle, set the caloparticleId for the hit to -1 this will - // contribute to the number of noise hits + // if the fraction is zero or the hit does not belong to any calo + // particle, set the caloparticleId for the hit to -1 this will + // contribute to the number of noise hits - // MR Remove the case in which the fraction is 0, since this could be a - // real hit that has been marked as halo. - if (rhFraction == 0.) { - hitsToCaloParticleId[hitId] = -2; - numberOfHaloHitsInMCL++; - } - if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) { - hitsToCaloParticleId[hitId] -= 1; - } else { - auto maxCPEnergyInLC = 0.f; - auto maxCPId = -1; - for (auto& h : hit_find_in_CP->second) { - auto shared_fraction = std::min(rhFraction, h.fraction); - //We are in the case where there are calo particles with simhits connected via detid with the rechit under study - //So, from all layers clusters, find the rechits that are connected with a calo particle and save/calculate the - //energy of that calo particle as the sum over all rechits of the rechits energy weighted - //by the caloparticle's fraction related to that rechit. - CPEnergyInMCL[h.clusterId] += shared_fraction * hit->energy(); - //Same but for layer clusters for the cell association per layer - CPEnergyInLC[h.clusterId] += shared_fraction * hit->energy(); - //Here cPOnLayer[caloparticle][layer] describe above is set. - //Here for multi clusters with matched rechit the CP fraction times hit energy is added and saved . - cPOnLayer[h.clusterId][lcLayerId].layerClusterIdToEnergyAndScore[mclId].first += - shared_fraction * hit->energy(); - cPOnLayer[h.clusterId][lcLayerId].layerClusterIdToEnergyAndScore[mclId].second = FLT_MAX; - //cpsInMultiCluster[multicluster][CPids] - //Connects a multi cluster with all related caloparticles. - cpsInMultiCluster[mclId].emplace_back(std::make_pair(h.clusterId, FLT_MAX)); - //From all CaloParticles related to a layer cluster, he saves id and energy of the calo particle - //that after simhit-rechit matching in layer has the maximum energy. - if (shared_fraction > maxCPEnergyInLC) { - //energy is used only here. cpid is saved for multiclusters - maxCPEnergyInLC = CPEnergyInLC[h.clusterId]; - maxCPId = h.clusterId; + // MR Remove the case in which the fraction is 0, since this could be a + // real hit that has been marked as halo. + if (rhFraction == 0.) { + hitsToCaloParticleId[hitId] = -2; + numberOfHaloHitsInMCL++; + } + if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) { + hitsToCaloParticleId[hitId] -= 1; + } else { + auto maxCPEnergyInLC = 0.f; + auto maxCPId = -1; + for (auto& h : hit_find_in_CP->second) { + auto shared_fraction = std::min(rhFraction, h.fraction); + //We are in the case where there are calo particles with simhits connected via detid with the rechit under study + //So, from all layers clusters, find the rechits that are connected with a calo particle and save/calculate the + //energy of that calo particle as the sum over all rechits of the rechits energy weighted + //by the caloparticle's fraction related to that rechit. + CPEnergyInMCL[h.clusterId] += shared_fraction * hit->energy(); + //Same but for layer clusters for the cell association per layer + CPEnergyInLC[h.clusterId] += shared_fraction * hit->energy(); + //Here cPOnLayer[caloparticle][layer] describe above is set. + //Here for multi clusters with matched rechit the CP fraction times hit energy is added and saved . + cPOnLayer[h.clusterId][lcLayerId].layerClusterIdToEnergyAndScore[mclId].first += + shared_fraction * hit->energy(); + cPOnLayer[h.clusterId][lcLayerId].layerClusterIdToEnergyAndScore[mclId].second = FLT_MAX; + //cpsInMultiCluster[multicluster][CPids] + //Connects a multi cluster with all related caloparticles. + cpsInMultiCluster[mclId].emplace_back(std::make_pair(h.clusterId, FLT_MAX)); + //From all CaloParticles related to a layer cluster, he saves id and energy of the calo particle + //that after simhit-rechit matching in layer has the maximum energy. + if (shared_fraction > maxCPEnergyInLC) { + //energy is used only here. cpid is saved for multiclusters + maxCPEnergyInLC = CPEnergyInLC[h.clusterId]; + maxCPId = h.clusterId; + } } + //Keep in mind here maxCPId could be zero. So, below ask for negative not including zero to count noise. + hitsToCaloParticleId[hitId] = maxCPId; } - //Keep in mind here maxCPId could be zero. So, below ask for negative not including zero to count noise. - hitsToCaloParticleId[hitId] = maxCPId; - } - } //end of loop through rechits of the layer cluster. + } //end of loop through rechits of the layer cluster. - //Loop through all rechits to count how many of them are noise and how many are matched. - //In case of matched rechit-simhit, he counts and saves the number of rechits related to the maximum energy CaloParticle. - for (auto& c : hitsToCaloParticleId) { - if (c < 0) { - numberOfNoiseHitsInMCL++; - } else { - occurrencesCPinMCL[c]++; + //Loop through all rechits to count how many of them are noise and how many are matched. + //In case of matched rechit-simhit, he counts and saves the number of rechits related to the maximum energy CaloParticle. + for (auto& c : hitsToCaloParticleId) { + if (c < 0) { + numberOfNoiseHitsInMCL++; + } else { + occurrencesCPinMCL[c]++; + } } - } - //Below from all maximum energy CaloParticles, he saves the one with the largest amount - //of related rechits. - for (auto& c : occurrencesCPinMCL) { - if (c.second > maxCPNumberOfHitsInMCL) { - maxCPId_byNumberOfHits = c.first; - maxCPNumberOfHitsInMCL = c.second; + //Below from all maximum energy CaloParticles, he saves the one with the largest amount + //of related rechits. + for (auto& c : occurrencesCPinMCL) { + if (c.second > maxCPNumberOfHitsInMCL) { + maxCPId_byNumberOfHits = c.first; + maxCPNumberOfHitsInMCL = c.second; + } } - } - //Find the CaloParticle that has the maximum energy shared with the multicluster under study. - for (auto& c : CPEnergyInMCL) { - if (c.second > maxEnergySharedMCLandCP) { - maxCPId_byEnergy = c.first; - maxEnergySharedMCLandCP = c.second; - } - } - //The energy of the CaloParticle that found to have the maximum energy shared with the multicluster under study. - float totalCPEnergyFromLayerCP = 0.f; - if (maxCPId_byEnergy >= 0) { - //Loop through all layers - for (unsigned int j = 0; j < layers * 2; ++j) { - totalCPEnergyFromLayerCP = totalCPEnergyFromLayerCP + cPOnLayer[maxCPId_byEnergy][j].energy; + //Find the CaloParticle that has the maximum energy shared with the multicluster under study. + for (auto& c : CPEnergyInMCL) { + if (c.second > maxEnergySharedMCLandCP) { + maxCPId_byEnergy = c.first; + maxEnergySharedMCLandCP = c.second; + } } - energyFractionOfCPinMCL = maxEnergySharedMCLandCP / totalCPEnergyFromLayerCP; - if (multiClusters[mclId].energy() > 0.f) { - energyFractionOfMCLinCP = maxEnergySharedMCLandCP / multiClusters[mclId].energy(); + //The energy of the CaloParticle that found to have the maximum energy shared with the multicluster under study. + float totalCPEnergyFromLayerCP = 0.f; + if (maxCPId_byEnergy >= 0) { + //Loop through all layers + for (unsigned int j = 0; j < layers * 2; ++j) { + totalCPEnergyFromLayerCP = totalCPEnergyFromLayerCP + cPOnLayer[maxCPId_byEnergy][j].energy; + } + energyFractionOfCPinMCL = maxEnergySharedMCLandCP / totalCPEnergyFromLayerCP; + if (multiClusters[mclId].energy() > 0.f) { + energyFractionOfMCLinCP = maxEnergySharedMCLandCP / multiClusters[mclId].energy(); + } } - } - - LogDebug("HGCalValidator") << std::setw(12) << "multiCluster" - << "\t" //LogDebug("HGCalValidator") - << std::setw(10) << "mulcl energy" - << "\t" << std::setw(5) << "nhits" - << "\t" << std::setw(12) << "noise hits" - << "\t" << std::setw(22) << "maxCPId_byNumberOfHits" - << "\t" << std::setw(8) << "nhitsCP" - << "\t" << std::setw(16) << "maxCPId_byEnergy" - << "\t" << std::setw(23) << "maxEnergySharedMCLandCP" - << "\t" << std::setw(22) << "totalCPEnergyFromAllLayerCP" - << "\t" << std::setw(22) << "energyFractionOfMCLinCP" - << "\t" << std::setw(25) << "energyFractionOfCPinMCL" - << "\t" << std::endl; - LogDebug("HGCalValidator") << std::setw(12) << mclId << "\t" //LogDebug("HGCalValidator") - << std::setw(10) << multiClusters[mclId].energy() << "\t" << std::setw(5) - << numberOfHitsInMCL << "\t" << std::setw(12) << numberOfNoiseHitsInMCL << "\t" - << std::setw(22) << maxCPId_byNumberOfHits << "\t" << std::setw(8) - << maxCPNumberOfHitsInMCL << "\t" << std::setw(16) << maxCPId_byEnergy << "\t" - << std::setw(23) << maxEnergySharedMCLandCP << "\t" << std::setw(22) - << totalCPEnergyFromLayerCP << "\t" << std::setw(22) << energyFractionOfMCLinCP << "\t" - << std::setw(25) << energyFractionOfCPinMCL << std::endl; - - } //end of loop through multi clusters + LogDebug("HGCalValidator") << std::setw(12) << "multiCluster" + << "\t" //LogDebug("HGCalValidator") + << std::setw(10) << "mulcl energy" + << "\t" << std::setw(5) << "nhits" + << "\t" << std::setw(12) << "noise hits" + << "\t" << std::setw(22) << "maxCPId_byNumberOfHits" + << "\t" << std::setw(8) << "nhitsCP" + << "\t" << std::setw(16) << "maxCPId_byEnergy" + << "\t" << std::setw(23) << "maxEnergySharedMCLandCP" + << "\t" << std::setw(22) << "totalCPEnergyFromAllLayerCP" + << "\t" << std::setw(22) << "energyFractionOfMCLinCP" + << "\t" << std::setw(25) << "energyFractionOfCPinMCL" + << "\t" << std::endl; + LogDebug("HGCalValidator") << std::setw(12) << mclId << "\t" //LogDebug("HGCalValidator") + << std::setw(10) << multiClusters[mclId].energy() << "\t" << std::setw(5) + << numberOfHitsInMCL << "\t" << std::setw(12) << numberOfNoiseHitsInMCL << "\t" + << std::setw(22) << maxCPId_byNumberOfHits << "\t" << std::setw(8) + << maxCPNumberOfHitsInMCL << "\t" << std::setw(16) << maxCPId_byEnergy << "\t" + << std::setw(23) << maxEnergySharedMCLandCP << "\t" << std::setw(22) + << totalCPEnergyFromLayerCP << "\t" << std::setw(22) << energyFractionOfMCLinCP << "\t" + << std::setw(25) << energyFractionOfCPinMCL << std::endl; + + } //end of loop through multi clusters + } //Loop through multiclusters for (unsigned int mclId = 0; mclId < nMultiClusters; ++mclId) { const auto& hits_and_fractions = multiClusters[mclId].hitsAndFractions(); @@ -1792,7 +1792,6 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist score->second, sharedeneCPallLayers / multiClusters[mclId].energy()); } } - auto assocFakeMerge = std::count_if(std::begin(cpsInMultiCluster[mclId]), std::end(cpsInMultiCluster[mclId]), [](const auto& obj) { return obj.second < ScoreCutMCLtoCPFakeMerge_; }); @@ -2116,18 +2115,20 @@ void HGVHistoProducerAlgo::fill_multi_cluster_histos(const Histograms& histogram //Since we want to also check for non contiguous multiclusters bool contimulti = false; //Observe that we start from 1 and go up to size - 1 element. - for (unsigned int i = 1; i < multicluster_layers_vec.size() - 1; ++i) { - if ((multicluster_layers_vec[i - 1] + 1 == multicluster_layers_vec[i]) && - (multicluster_layers_vec[i + 1] - 1 == multicluster_layers_vec[i])) { - //So, this is a multicluster with 3 contiguous layers per event - if (multiclusterInZplus) { - tncontmclpz++; - } - if (multiclusterInZminus) { - tncontmclmz++; + if (multicluster_layers_vec.size() >= 3) { + for (unsigned int i = 1; i < multicluster_layers_vec.size() - 1; ++i) { + if ((multicluster_layers_vec[i - 1] + 1 == multicluster_layers_vec[i]) && + (multicluster_layers_vec[i + 1] - 1 == multicluster_layers_vec[i])) { + //So, this is a multicluster with 3 contiguous layers per event + if (multiclusterInZplus) { + tncontmclpz++; + } + if (multiclusterInZminus) { + tncontmclmz++; + } + contimulti = true; + break; } - contimulti = true; - break; } } //Count non contiguous multiclusters @@ -2168,16 +2169,20 @@ void HGVHistoProducerAlgo::fill_multi_cluster_histos(const Histograms& histogram histograms.h_multiplicityOfLCinMCL_vs_layerclusterenergy[count]->Fill(mlp, layerClusters[lc]->energy()); } + if (!multicluster_layers.empty()) { + histograms.h_multicluster_x[count]->Fill(multiClusters[mclId].x()); + histograms.h_multicluster_y[count]->Fill(multiClusters[mclId].y()); + histograms.h_multicluster_z[count]->Fill(multiClusters[mclId].z()); + histograms.h_multicluster_eta[count]->Fill(multiClusters[mclId].eta()); + histograms.h_multicluster_phi[count]->Fill(multiClusters[mclId].phi()); + + histograms.h_multicluster_firstlayer[count]->Fill((float)*multicluster_layers.begin()); + histograms.h_multicluster_lastlayer[count]->Fill((float)*multicluster_layers.rbegin()); + histograms.h_multicluster_layersnum[count]->Fill((float)multicluster_layers.size()); + } histograms.h_multicluster_pt[count]->Fill(multiClusters[mclId].pt()); - histograms.h_multicluster_eta[count]->Fill(multiClusters[mclId].eta()); - histograms.h_multicluster_phi[count]->Fill(multiClusters[mclId].phi()); + histograms.h_multicluster_energy[count]->Fill(multiClusters[mclId].energy()); - histograms.h_multicluster_x[count]->Fill(multiClusters[mclId].x()); - histograms.h_multicluster_y[count]->Fill(multiClusters[mclId].y()); - histograms.h_multicluster_z[count]->Fill(multiClusters[mclId].z()); - histograms.h_multicluster_firstlayer[count]->Fill((float)*multicluster_layers.begin()); - histograms.h_multicluster_lastlayer[count]->Fill((float)*multicluster_layers.rbegin()); - histograms.h_multicluster_layersnum[count]->Fill((float)multicluster_layers.size()); } //end of loop through multiclusters From 56d168acdb70f5d42de39a4cd09cfafbe5772754 Mon Sep 17 00:00:00 2001 From: Felice Date: Mon, 12 Oct 2020 11:20:49 +0200 Subject: [PATCH 23/66] Increasing missing outer hits from 3 to 5 --- RecoHGCal/TICL/plugins/TICLSeedingRegionProducer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/TICLSeedingRegionProducer.cc b/RecoHGCal/TICL/plugins/TICLSeedingRegionProducer.cc index e42abbd1c43c7..b886b69aef8e0 100644 --- a/RecoHGCal/TICL/plugins/TICLSeedingRegionProducer.cc +++ b/RecoHGCal/TICL/plugins/TICLSeedingRegionProducer.cc @@ -59,7 +59,7 @@ void TICLSeedingRegionProducer::fillDescriptions(edm::ConfigurationDescriptions& desc.add("tracks", edm::InputTag("generalTracks")); desc.add("cutTk", "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && " - "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 3"); + "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5"); desc.add("propagator", "PropagatorWithMaterial"); desc.add("algoId", 1); descriptions.add("ticlSeedingRegionProducer", desc); From 9f62bea1253786bf9cbb6a0798097c048d902d3f Mon Sep 17 00:00:00 2001 From: Felice Date: Mon, 12 Oct 2020 19:11:20 +0200 Subject: [PATCH 24/66] increase missing layer from 0 to 1 and loosen min_cos_theta to 0.961 in EM iterations --- RecoHGCal/TICL/python/EMStep_cff.py | 4 ++-- RecoHGCal/TICL/python/TrkEMStep_cff.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index e8fe29101959f..f3b9495e19177 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -27,9 +27,9 @@ energy_em_over_total_threshold = 0.85, shower_start_max_layer = 7, #inclusive max_out_in_hops = 4, - missing_layers = 0, + missing_layers = 1, min_clusters_per_ntuplet = 10, - min_cos_theta = 0.978, # ~12 degrees + min_cos_theta = 0.961, # ~16 degrees min_cos_pointing = 0.9, # ~25 degrees max_delta_time = 3., itername = "EM", diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index ec73608400394..b9be4b751ad63 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -25,9 +25,9 @@ energy_em_over_total_threshold = 0.85, shower_start_max_layer = 7, #inclusive max_out_in_hops = 4, - missing_layers = 0, + missing_layers = 1, min_clusters_per_ntuplet = 10, - min_cos_theta = 0.978, # ~12 degrees + min_cos_theta = 0.961, # ~16 degrees min_cos_pointing = 0.9, # ~25 degrees max_delta_time = 3., itername = "TrkEM", From d6af9c3ebe08d37953d5a77b41da06a0c02898db Mon Sep 17 00:00:00 2001 From: Felice Date: Mon, 12 Oct 2020 19:14:20 +0200 Subject: [PATCH 25/66] tighten min_cos_theta from 0.961 to 0.978 in EM iterations --- RecoHGCal/TICL/python/EMStep_cff.py | 2 +- RecoHGCal/TICL/python/TrkEMStep_cff.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index f3b9495e19177..24e90fe3d468f 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -29,7 +29,7 @@ max_out_in_hops = 4, missing_layers = 1, min_clusters_per_ntuplet = 10, - min_cos_theta = 0.961, # ~16 degrees + min_cos_theta = 0.978, # ~12 degrees min_cos_pointing = 0.9, # ~25 degrees max_delta_time = 3., itername = "EM", diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index b9be4b751ad63..be891f413713e 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -27,7 +27,7 @@ max_out_in_hops = 4, missing_layers = 1, min_clusters_per_ntuplet = 10, - min_cos_theta = 0.961, # ~16 degrees + min_cos_theta = 0.978, # ~12 degrees min_cos_pointing = 0.9, # ~25 degrees max_delta_time = 3., itername = "TrkEM", From 0b56bd390e7817e5d5babf4787676171644f1477 Mon Sep 17 00:00:00 2001 From: Felice Date: Tue, 13 Oct 2020 11:53:50 +0200 Subject: [PATCH 26/66] Rename parameter missing_layers to skip_layers. --- RecoHGCal/TICL/plugins/HGCGraph.cc | 4 ++-- RecoHGCal/TICL/plugins/HGCGraph.h | 2 +- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc | 4 ++-- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h | 2 +- RecoHGCal/TICL/plugins/TrackstersProducer.cc | 2 +- RecoHGCal/TICL/python/EMStep_cff.py | 2 +- RecoHGCal/TICL/python/HADStep_cff.py | 2 +- RecoHGCal/TICL/python/MIPStep_cff.py | 2 +- RecoHGCal/TICL/python/TrkEMStep_cff.py | 2 +- RecoHGCal/TICL/python/TrkStep_cff.py | 2 +- RecoHGCal/TICL/python/ticl_iterations.py | 12 ++++++------ 11 files changed, 18 insertions(+), 18 deletions(-) diff --git a/RecoHGCal/TICL/plugins/HGCGraph.cc b/RecoHGCal/TICL/plugins/HGCGraph.cc index 1599d6c98df5d..1be81fb9a6f1b 100644 --- a/RecoHGCal/TICL/plugins/HGCGraph.cc +++ b/RecoHGCal/TICL/plugins/HGCGraph.cc @@ -20,7 +20,7 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, float minCosTheta, float minCosPointing, float etaLimitIncreaseWindow, - int missing_layers, + int skip_layers, int maxNumberOfLayers, float maxDeltaTime) { isOuterClusterOfDoublets_.clear(); @@ -75,7 +75,7 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, } for (int il = 0; il < maxNumberOfLayers - 1; ++il) { - for (int outer_layer = 0; outer_layer < std::min(1 + missing_layers, maxNumberOfLayers - 1 - il); ++outer_layer) { + for (int outer_layer = 0; outer_layer < std::min(1 + skip_layers, maxNumberOfLayers - 1 - il); ++outer_layer) { int currentInnerLayerId = il + maxNumberOfLayers * zSide; int currentOuterLayerId = currentInnerLayerId + 1 + outer_layer; auto const &outerLayerHisto = histo[currentOuterLayerId]; diff --git a/RecoHGCal/TICL/plugins/HGCGraph.h b/RecoHGCal/TICL/plugins/HGCGraph.h index 103f76b312300..e9ee6da110821 100644 --- a/RecoHGCal/TICL/plugins/HGCGraph.h +++ b/RecoHGCal/TICL/plugins/HGCGraph.h @@ -26,7 +26,7 @@ class HGCGraphT { float minCosThetai, float maxCosPointing, float etaLimitIncreaseWindow, - int missing_layers, + int skip_layers, int maxNumberOfLayers, float maxDeltaTime); diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index 653edd5710858..60a6e4c8f25ca 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -27,7 +27,7 @@ PatternRecognitionbyCA::PatternRecognitionbyCA(const edm::ParameterSet &c min_cos_theta_(conf.getParameter("min_cos_theta")), min_cos_pointing_(conf.getParameter("min_cos_pointing")), etaLimitIncreaseWindow_(conf.getParameter("etaLimitIncreaseWindow")), - missing_layers_(conf.getParameter("missing_layers")), + skip_layers_(conf.getParameter("skip_layers")), shower_start_max_layer_(conf.getParameter("shower_start_max_layer")), min_clusters_per_ntuplet_(conf.getParameter("min_clusters_per_ntuplet")), max_delta_time_(conf.getParameter("max_delta_time")), @@ -90,7 +90,7 @@ void PatternRecognitionbyCA::makeTracksters( min_cos_theta_, min_cos_pointing_, etaLimitIncreaseWindow_, - missing_layers_, + skip_layers_, rhtools_.lastLayer(type), max_delta_time_); diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h index 0ed31dac06f44..6cf977675e8dc 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h @@ -38,7 +38,7 @@ namespace ticl { const float min_cos_theta_; const float min_cos_pointing_; const float etaLimitIncreaseWindow_; - const int missing_layers_; + const int skip_layers_; const unsigned int shower_start_max_layer_; const int min_clusters_per_ntuplet_; diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc index 16462d76612e7..eb00738963849 100644 --- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc @@ -127,7 +127,7 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri desc.add("algo_verbosity", 0); desc.add("min_cos_theta", 0.915); desc.add("min_cos_pointing", -1.); - desc.add("missing_layers", 0); + desc.add("skip_layers", 0); desc.add("etaLimitIncreaseWindow", 2.1); desc.add("min_clusters_per_ntuplet", 10); desc.add("max_delta_time", 3.); //nsigma diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 24e90fe3d468f..90ad96a0f2a23 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -27,7 +27,7 @@ energy_em_over_total_threshold = 0.85, shower_start_max_layer = 7, #inclusive max_out_in_hops = 4, - missing_layers = 1, + skip_layers = 1, min_clusters_per_ntuplet = 10, min_cos_theta = 0.978, # ~12 degrees min_cos_pointing = 0.9, # ~25 degrees diff --git a/RecoHGCal/TICL/python/HADStep_cff.py b/RecoHGCal/TICL/python/HADStep_cff.py index d631f7196d137..e1eed6eb81917 100644 --- a/RecoHGCal/TICL/python/HADStep_cff.py +++ b/RecoHGCal/TICL/python/HADStep_cff.py @@ -25,7 +25,7 @@ # For the moment we mask everything w/o requirements since we are last # filter_on_categories = [5], # filter neutral hadrons # pid_threshold = 0.7, - missing_layers = 1, + skip_layers = 1, min_clusters_per_ntuplet = 12, min_cos_theta = 0.866, # ~30 degrees min_cos_pointing = 0.819, # ~35 degrees diff --git a/RecoHGCal/TICL/python/MIPStep_cff.py b/RecoHGCal/TICL/python/MIPStep_cff.py index da90839a6b90f..75146b1783ffe 100644 --- a/RecoHGCal/TICL/python/MIPStep_cff.py +++ b/RecoHGCal/TICL/python/MIPStep_cff.py @@ -20,7 +20,7 @@ ticlTrackstersMIP = _trackstersProducer.clone( filtered_mask = "filteredLayerClustersMIP:MIP", seeding_regions = "ticlSeedingGlobal", - missing_layers = 3, + skip_layers = 3, min_clusters_per_ntuplet = 10, min_cos_theta = 0.99, # ~10 degrees min_cos_pointing = 0.5, diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index be891f413713e..53e04dbff080b 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -25,7 +25,7 @@ energy_em_over_total_threshold = 0.85, shower_start_max_layer = 7, #inclusive max_out_in_hops = 4, - missing_layers = 1, + skip_layers = 1, min_clusters_per_ntuplet = 10, min_cos_theta = 0.978, # ~12 degrees min_cos_pointing = 0.9, # ~25 degrees diff --git a/RecoHGCal/TICL/python/TrkStep_cff.py b/RecoHGCal/TICL/python/TrkStep_cff.py index e35a692ff194b..af995cdb995e1 100644 --- a/RecoHGCal/TICL/python/TrkStep_cff.py +++ b/RecoHGCal/TICL/python/TrkStep_cff.py @@ -24,7 +24,7 @@ original_mask = 'ticlTrackstersEM', filter_on_categories = [2, 4], # filter muons and charged hadrons pid_threshold = 0.0, - missing_layers = 3, + skip_layers = 3, min_clusters_per_ntuplet = 10, # min_cos_theta = 0.978, # same as em # min_cos_pointing = 0.9, # same as em diff --git a/RecoHGCal/TICL/python/ticl_iterations.py b/RecoHGCal/TICL/python/ticl_iterations.py index b36d459990c70..d6b3ea87d0ef9 100644 --- a/RecoHGCal/TICL/python/ticl_iterations.py +++ b/RecoHGCal/TICL/python/ticl_iterations.py @@ -41,7 +41,7 @@ def TICL_iterations_withReco(process): process.trackstersTrk = trackstersProducer.clone( filtered_mask = "filteredLayerClustersTrk:Trk", seeding_regions = "ticlSeedingTrk", - missing_layers = 3, + skip_layers = 3, min_clusters_per_ntuplet = 5, min_cos_theta = 0.99, # ~10 degrees min_cos_pointing = 0.9 @@ -66,7 +66,7 @@ def TICL_iterations_withReco(process): process.trackstersMIP = trackstersProducer.clone( filtered_mask = "filteredLayerClustersMIP:MIP", seeding_regions = "ticlSeedingGlobal", - missing_layers = 3, + skip_layers = 3, min_clusters_per_ntuplet = 15, min_cos_theta = 0.99, # ~10 degrees min_cos_pointing = 0.9, @@ -91,7 +91,7 @@ def TICL_iterations_withReco(process): original_mask = "trackstersMIP", filtered_mask = "filteredLayerClusters:algo8", seeding_regions = "ticlSeedingGlobal", - missing_layers = 1, + skip_layers = 1, min_clusters_per_ntuplet = 10, min_cos_theta = 0.984, # ~10 degrees min_cos_pointing = 0.9 # ~26 degrees @@ -105,7 +105,7 @@ def TICL_iterations_withReco(process): process.trackstersHAD = trackstersProducer.clone( filtered_mask = "filteredLayerClusters:algo8", seeding_regions = "ticlSeedingGlobal", - missing_layers = 2, + skip_layers = 2, min_clusters_per_ntuplet = 10, min_cos_theta = 0.8, min_cos_pointing = 0.7 @@ -172,7 +172,7 @@ def TICL_iterations(process): process.trackstersMIP = trackstersProducer.clone( filtered_mask = "filteredLayerClustersMIP:MIP", seeding_regions = "ticlSeedingGlobal", - missing_layers = 3, + skip_layers = 3, min_clusters_per_ntuplet = 15, min_cos_theta = 0.99, # ~10 degrees ) @@ -193,7 +193,7 @@ def TICL_iterations(process): original_mask = "trackstersMIP", filtered_mask = "filteredLayerClusters:algo8", seeding_regions = "ticlSeedingGlobal", - missing_layers = 2, + skip_layers = 2, min_clusters_per_ntuplet = 15, min_cos_theta = 0.94, # ~20 degrees min_cos_pointing = 0.7 From 7d80a780f48ba1f257be31a95ddfdded68940160 Mon Sep 17 00:00:00 2001 From: Felice Date: Tue, 13 Oct 2020 13:33:05 +0200 Subject: [PATCH 27/66] Retuning of PR cuts in EM iterations --- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc | 1 + RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h | 1 + RecoHGCal/TICL/plugins/TrackstersProducer.cc | 1 + RecoHGCal/TICL/python/EMStep_cff.py | 7 ++++--- RecoHGCal/TICL/python/TrkEMStep_cff.py | 6 +++--- 5 files changed, 10 insertions(+), 6 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index 60a6e4c8f25ca..a1fe363ce6dd0 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -28,6 +28,7 @@ PatternRecognitionbyCA::PatternRecognitionbyCA(const edm::ParameterSet &c min_cos_pointing_(conf.getParameter("min_cos_pointing")), etaLimitIncreaseWindow_(conf.getParameter("etaLimitIncreaseWindow")), skip_layers_(conf.getParameter("skip_layers")), + max_missing_layers_in_trackster_(conf.getParameter("max_missing_layers_in_trackster")), shower_start_max_layer_(conf.getParameter("shower_start_max_layer")), min_clusters_per_ntuplet_(conf.getParameter("min_clusters_per_ntuplet")), max_delta_time_(conf.getParameter("max_delta_time")), diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h index 6cf977675e8dc..1147c6bfe197e 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h @@ -39,6 +39,7 @@ namespace ticl { const float min_cos_pointing_; const float etaLimitIncreaseWindow_; const int skip_layers_; + const int max_missing_layers_in_trackster_; const unsigned int shower_start_max_layer_; const int min_clusters_per_ntuplet_; diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc index eb00738963849..aa63e248bc6f1 100644 --- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc @@ -128,6 +128,7 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri desc.add("min_cos_theta", 0.915); desc.add("min_cos_pointing", -1.); desc.add("skip_layers", 0); + desc.add("max_missing_layers_in_trackster", 9999); desc.add("etaLimitIncreaseWindow", 2.1); desc.add("min_clusters_per_ntuplet", 10); desc.add("max_delta_time", 3.); //nsigma diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 90ad96a0f2a23..1ff7184506bc3 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -26,11 +26,12 @@ pid_threshold = 0.5, energy_em_over_total_threshold = 0.85, shower_start_max_layer = 7, #inclusive - max_out_in_hops = 4, + max_out_in_hops = 1, skip_layers = 1, + max_missing_layers_in_trackster = 1, min_clusters_per_ntuplet = 10, - min_cos_theta = 0.978, # ~12 degrees - min_cos_pointing = 0.9, # ~25 degrees + min_cos_theta = 0.97, # ~14 degrees + min_cos_pointing = 0.94, # ~20 degrees max_delta_time = 3., itername = "EM", algo_verbosity = 0, diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index 53e04dbff080b..ae41eaffb6719 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -24,11 +24,11 @@ pid_threshold = 0.5, energy_em_over_total_threshold = 0.85, shower_start_max_layer = 7, #inclusive - max_out_in_hops = 4, + max_out_in_hops = 1, skip_layers = 1, min_clusters_per_ntuplet = 10, - min_cos_theta = 0.978, # ~12 degrees - min_cos_pointing = 0.9, # ~25 degrees + min_cos_theta = 0.97, # ~14 degrees + min_cos_pointing = 0.94, # ~20 degrees max_delta_time = 3., itername = "TrkEM", algo_verbosity = 0, From 42413ccac794bbb6e2a181ec8f96e56bd343be23 Mon Sep 17 00:00:00 2001 From: Felice Date: Tue, 13 Oct 2020 17:19:18 +0200 Subject: [PATCH 28/66] Introducing two cuts: a max_missing_layers_in_trackster and a minimum number of different layers in a trackster --- .../TICL/plugins/PatternRecognitionbyCA.cc | 47 +++++++++++++++++-- .../TICL/plugins/PatternRecognitionbyCA.h | 3 +- RecoHGCal/TICL/plugins/TrackstersProducer.cc | 4 +- RecoHGCal/TICL/python/EMStep_cff.py | 4 +- RecoHGCal/TICL/python/HADStep_cff.py | 2 +- RecoHGCal/TICL/python/MIPStep_cff.py | 4 +- RecoHGCal/TICL/python/TrkEMStep_cff.py | 2 +- RecoHGCal/TICL/python/TrkStep_cff.py | 2 +- RecoHGCal/TICL/python/ticl_iterations.py | 12 ++--- 9 files changed, 60 insertions(+), 20 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index a1fe363ce6dd0..a7b85c8e72f68 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -30,7 +30,8 @@ PatternRecognitionbyCA::PatternRecognitionbyCA(const edm::ParameterSet &c skip_layers_(conf.getParameter("skip_layers")), max_missing_layers_in_trackster_(conf.getParameter("max_missing_layers_in_trackster")), shower_start_max_layer_(conf.getParameter("shower_start_max_layer")), - min_clusters_per_ntuplet_(conf.getParameter("min_clusters_per_ntuplet")), + min_layers_per_trackster_(conf.getParameter("min_layers_per_trackster")), + min_clusters_per_ntuplet_(min_layers_per_trackster_), max_delta_time_(conf.getParameter("max_delta_time")), eidInputName_(conf.getParameter("eid_input_name")), eidOutputNameEnergy_(conf.getParameter("eid_output_name_energy")), @@ -46,6 +47,9 @@ PatternRecognitionbyCA::PatternRecognitionbyCA(const edm::ParameterSet &c << "PatternRecognitionbyCA received an empty graph definition from the global cache"; } eidSession_ = tensorflow::createSession(trackstersCache->eidGraphDef); + if (max_missing_layers_in_trackster_ < 100) { + check_missing_layers_ = true; + } } template @@ -140,12 +144,47 @@ void PatternRecognitionbyCA::makeTracksters( } } unsigned showerMinLayerId = 99999; + std::vector layerIds; + std::vector uniqueLayerIds; + uniqueLayerIds.reserve(effective_cluster_idx.size()); + std::vector> lcIdAndLayer; + lcIdAndLayer.reserve(effective_cluster_idx.size()); for (auto const i : effective_cluster_idx) { - auto const& haf = input.layerClusters[i].hitsAndFractions(); - showerMinLayerId = std::min(rhtools_.getLayerWithOffset(haf[0].first),showerMinLayerId); + auto const &haf = input.layerClusters[i].hitsAndFractions(); + auto layerId = rhtools_.getLayerWithOffset(haf[0].first); + showerMinLayerId = std::min(layerId, showerMinLayerId); + uniqueLayerIds.push_back(layerId); + lcIdAndLayer.emplace_back(i, layerId); + } + + if (check_missing_layers_) { + std::sort(uniqueLayerIds.begin(), uniqueLayerIds.end()); + uniqueLayerIds.erase(std::unique(uniqueLayerIds.begin(), uniqueLayerIds.end()), uniqueLayerIds.end()); + int numberOfMissingLayers = 0; + unsigned int j = showerMinLayerId; + unsigned int indexInVec = 0; + for (auto &layer : uniqueLayerIds) { + if (layer != j) { + numberOfMissingLayers++; + j++; + if (numberOfMissingLayers == max_missing_layers_in_trackster_) { + uniqueLayerIds.erase(uniqueLayerIds.begin() + indexInVec, uniqueLayerIds.end()); + for (auto &llpair : lcIdAndLayer) { + if (llpair.second >= layer) { + effective_cluster_idx.erase(llpair.first); + } + } + break; + } + } + indexInVec++; + j++; + } } - if(showerMinLayerId <= shower_start_max_layer_){ + bool selected = + (uniqueLayerIds.size() >= min_layers_per_trackster_) and (showerMinLayerId <= shower_start_max_layer_); + if (selected) { for (auto const i : effective_cluster_idx) { layer_cluster_usage[i]++; if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h index 1147c6bfe197e..d67526585bd74 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h @@ -40,8 +40,9 @@ namespace ticl { const float etaLimitIncreaseWindow_; const int skip_layers_; const int max_missing_layers_in_trackster_; + bool check_missing_layers_ = false; const unsigned int shower_start_max_layer_; - + const unsigned int min_layers_per_trackster_; const int min_clusters_per_ntuplet_; const float max_delta_time_; const std::string eidInputName_; diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc index aa63e248bc6f1..5aa97b4a9b0f6 100644 --- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc @@ -123,14 +123,14 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri desc.add>("filter_on_categories", {0}); desc.add("pid_threshold", 0.); // make default such that no filtering is applied desc.add("energy_em_over_total_threshold", -1.); // make default such that no filtering is applied - desc.add("shower_start_max_layer", 9999); // make default such that no filtering is applied + desc.add("shower_start_max_layer", 9999); // make default such that no filtering is applied desc.add("algo_verbosity", 0); desc.add("min_cos_theta", 0.915); desc.add("min_cos_pointing", -1.); desc.add("skip_layers", 0); desc.add("max_missing_layers_in_trackster", 9999); desc.add("etaLimitIncreaseWindow", 2.1); - desc.add("min_clusters_per_ntuplet", 10); + desc.add("min_layers_per_trackster", 10); desc.add("max_delta_time", 3.); //nsigma desc.add("out_in_dfs", true); desc.add("max_out_in_hops", 10); diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 1ff7184506bc3..e2b5b48e814b0 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -29,7 +29,7 @@ max_out_in_hops = 1, skip_layers = 1, max_missing_layers_in_trackster = 1, - min_clusters_per_ntuplet = 10, + min_layers_per_trackster = 10, min_cos_theta = 0.97, # ~14 degrees min_cos_pointing = 0.94, # ~20 degrees max_delta_time = 3., @@ -64,7 +64,7 @@ filtered_mask = "filteredLayerClustersHFNoseEM:EMn", seeding_regions = "ticlSeedingGlobalHFNose", time_layerclusters = "hgcalLayerClustersHFNose:timeLayerCluster", - min_clusters_per_ntuplet = 6 + min_layers_per_trackster = 6 ) ticlHFNoseEMStepTask = cms.Task(ticlSeedingGlobalHFNose diff --git a/RecoHGCal/TICL/python/HADStep_cff.py b/RecoHGCal/TICL/python/HADStep_cff.py index e1eed6eb81917..72efaebc201ec 100644 --- a/RecoHGCal/TICL/python/HADStep_cff.py +++ b/RecoHGCal/TICL/python/HADStep_cff.py @@ -26,7 +26,7 @@ # filter_on_categories = [5], # filter neutral hadrons # pid_threshold = 0.7, skip_layers = 1, - min_clusters_per_ntuplet = 12, + min_layers_per_trackster = 12, min_cos_theta = 0.866, # ~30 degrees min_cos_pointing = 0.819, # ~35 degrees max_delta_time = -1, diff --git a/RecoHGCal/TICL/python/MIPStep_cff.py b/RecoHGCal/TICL/python/MIPStep_cff.py index 75146b1783ffe..f7e5234b3aa9c 100644 --- a/RecoHGCal/TICL/python/MIPStep_cff.py +++ b/RecoHGCal/TICL/python/MIPStep_cff.py @@ -21,7 +21,7 @@ filtered_mask = "filteredLayerClustersMIP:MIP", seeding_regions = "ticlSeedingGlobal", skip_layers = 3, - min_clusters_per_ntuplet = 10, + min_layers_per_trackster = 10, min_cos_theta = 0.99, # ~10 degrees min_cos_pointing = 0.5, out_in_dfs = False, @@ -55,7 +55,7 @@ filtered_mask = "filteredLayerClustersHFNoseMIP:MIPn", seeding_regions = "ticlSeedingGlobalHFNose", time_layerclusters = "hgcalLayerClustersHFNose:timeLayerCluster", - min_clusters_per_ntuplet = 6 + min_layers_per_trackster = 6 ) ticlHFNoseMIPStepTask = cms.Task(ticlSeedingGlobalHFNose diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index ae41eaffb6719..6208c4e018d9f 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -26,7 +26,7 @@ shower_start_max_layer = 7, #inclusive max_out_in_hops = 1, skip_layers = 1, - min_clusters_per_ntuplet = 10, + min_layers_per_trackster = 10, min_cos_theta = 0.97, # ~14 degrees min_cos_pointing = 0.94, # ~20 degrees max_delta_time = 3., diff --git a/RecoHGCal/TICL/python/TrkStep_cff.py b/RecoHGCal/TICL/python/TrkStep_cff.py index af995cdb995e1..697441571b52e 100644 --- a/RecoHGCal/TICL/python/TrkStep_cff.py +++ b/RecoHGCal/TICL/python/TrkStep_cff.py @@ -25,7 +25,7 @@ filter_on_categories = [2, 4], # filter muons and charged hadrons pid_threshold = 0.0, skip_layers = 3, - min_clusters_per_ntuplet = 10, + min_layers_per_trackster = 10, # min_cos_theta = 0.978, # same as em # min_cos_pointing = 0.9, # same as em min_cos_theta = 0.866, # ~30 degrees diff --git a/RecoHGCal/TICL/python/ticl_iterations.py b/RecoHGCal/TICL/python/ticl_iterations.py index d6b3ea87d0ef9..6266c07807603 100644 --- a/RecoHGCal/TICL/python/ticl_iterations.py +++ b/RecoHGCal/TICL/python/ticl_iterations.py @@ -42,7 +42,7 @@ def TICL_iterations_withReco(process): filtered_mask = "filteredLayerClustersTrk:Trk", seeding_regions = "ticlSeedingTrk", skip_layers = 3, - min_clusters_per_ntuplet = 5, + min_layers_per_trackster = 5, min_cos_theta = 0.99, # ~10 degrees min_cos_pointing = 0.9 ) @@ -67,7 +67,7 @@ def TICL_iterations_withReco(process): filtered_mask = "filteredLayerClustersMIP:MIP", seeding_regions = "ticlSeedingGlobal", skip_layers = 3, - min_clusters_per_ntuplet = 15, + min_layers_per_trackster = 15, min_cos_theta = 0.99, # ~10 degrees min_cos_pointing = 0.9, out_in_dfs = False, @@ -92,7 +92,7 @@ def TICL_iterations_withReco(process): filtered_mask = "filteredLayerClusters:algo8", seeding_regions = "ticlSeedingGlobal", skip_layers = 1, - min_clusters_per_ntuplet = 10, + min_layers_per_trackster = 10, min_cos_theta = 0.984, # ~10 degrees min_cos_pointing = 0.9 # ~26 degrees ) @@ -106,7 +106,7 @@ def TICL_iterations_withReco(process): filtered_mask = "filteredLayerClusters:algo8", seeding_regions = "ticlSeedingGlobal", skip_layers = 2, - min_clusters_per_ntuplet = 10, + min_layers_per_trackster = 10, min_cos_theta = 0.8, min_cos_pointing = 0.7 ) @@ -173,7 +173,7 @@ def TICL_iterations(process): filtered_mask = "filteredLayerClustersMIP:MIP", seeding_regions = "ticlSeedingGlobal", skip_layers = 3, - min_clusters_per_ntuplet = 15, + min_layers_per_trackster = 15, min_cos_theta = 0.99, # ~10 degrees ) @@ -194,7 +194,7 @@ def TICL_iterations(process): filtered_mask = "filteredLayerClusters:algo8", seeding_regions = "ticlSeedingGlobal", skip_layers = 2, - min_clusters_per_ntuplet = 15, + min_layers_per_trackster = 15, min_cos_theta = 0.94, # ~20 degrees min_cos_pointing = 0.7 ) From b93c8adeec32d4ac57c311692cfc3bfe8bcba00a Mon Sep 17 00:00:00 2001 From: Felice Date: Tue, 13 Oct 2020 17:49:17 +0200 Subject: [PATCH 29/66] Refactoring --- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index a7b85c8e72f68..4935f96323679 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -156,10 +156,11 @@ void PatternRecognitionbyCA::makeTracksters( uniqueLayerIds.push_back(layerId); lcIdAndLayer.emplace_back(i, layerId); } - + std::sort(uniqueLayerIds.begin(), uniqueLayerIds.end()); + uniqueLayerIds.erase(std::unique(uniqueLayerIds.begin(), uniqueLayerIds.end()), uniqueLayerIds.end()); + unsigned int numberOfLayersInTrackster = uniqueLayerIds.size(); if (check_missing_layers_) { - std::sort(uniqueLayerIds.begin(), uniqueLayerIds.end()); - uniqueLayerIds.erase(std::unique(uniqueLayerIds.begin(), uniqueLayerIds.end()), uniqueLayerIds.end()); + int numberOfMissingLayers = 0; unsigned int j = showerMinLayerId; unsigned int indexInVec = 0; @@ -168,7 +169,7 @@ void PatternRecognitionbyCA::makeTracksters( numberOfMissingLayers++; j++; if (numberOfMissingLayers == max_missing_layers_in_trackster_) { - uniqueLayerIds.erase(uniqueLayerIds.begin() + indexInVec, uniqueLayerIds.end()); + numberOfLayersInTrackster = indexInVec; for (auto &llpair : lcIdAndLayer) { if (llpair.second >= layer) { effective_cluster_idx.erase(llpair.first); @@ -183,7 +184,7 @@ void PatternRecognitionbyCA::makeTracksters( } bool selected = - (uniqueLayerIds.size() >= min_layers_per_trackster_) and (showerMinLayerId <= shower_start_max_layer_); + (numberOfLayersInTrackster >= min_layers_per_trackster_) and (showerMinLayerId <= shower_start_max_layer_); if (selected) { for (auto const i : effective_cluster_idx) { layer_cluster_usage[i]++; From 12518390ff75594eda5eb55faef77f134816489d Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Tue, 13 Oct 2020 19:06:44 +0200 Subject: [PATCH 30/66] Update cuts --- RecoHGCal/TICL/python/EMStep_cff.py | 8 ++++---- RecoHGCal/TICL/python/TrkEMStep_cff.py | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index e2b5b48e814b0..0346c7d762639 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -24,14 +24,14 @@ seeding_regions = "ticlSeedingGlobal", filter_on_categories = [0, 1], pid_threshold = 0.5, - energy_em_over_total_threshold = 0.85, - shower_start_max_layer = 7, #inclusive + energy_em_over_total_threshold = 0.9, + shower_start_max_layer = 5, #inclusive max_out_in_hops = 1, - skip_layers = 1, + skip_layers = 2, max_missing_layers_in_trackster = 1, min_layers_per_trackster = 10, min_cos_theta = 0.97, # ~14 degrees - min_cos_pointing = 0.94, # ~20 degrees + min_cos_pointing = 0.9, # ~25 degrees max_delta_time = 3., itername = "EM", algo_verbosity = 0, diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index 6208c4e018d9f..67b5d91d73ae9 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -22,10 +22,10 @@ seeding_regions = "ticlSeedingTrk", filter_on_categories = [0, 1], pid_threshold = 0.5, - energy_em_over_total_threshold = 0.85, - shower_start_max_layer = 7, #inclusive + energy_em_over_total_threshold = 0.9, + shower_start_max_layer = 5, #inclusive max_out_in_hops = 1, - skip_layers = 1, + skip_layers = 2, min_layers_per_trackster = 10, min_cos_theta = 0.97, # ~14 degrees min_cos_pointing = 0.94, # ~20 degrees From 4c54e8b01fdb329da7daf09b0981ce67c77c478e Mon Sep 17 00:00:00 2001 From: Felice Date: Tue, 13 Oct 2020 19:23:42 +0200 Subject: [PATCH 31/66] Make max_missing_layers_in_trackster cut inclusive --- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index 4935f96323679..088c8f9ac8884 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -168,7 +168,7 @@ void PatternRecognitionbyCA::makeTracksters( if (layer != j) { numberOfMissingLayers++; j++; - if (numberOfMissingLayers == max_missing_layers_in_trackster_) { + if (numberOfMissingLayers > max_missing_layers_in_trackster_) { numberOfLayersInTrackster = indexInVec; for (auto &llpair : lcIdAndLayer) { if (llpair.second >= layer) { From 97b951b8a6068a87b7295d8954192c04a5b902a8 Mon Sep 17 00:00:00 2001 From: Felice Date: Tue, 13 Oct 2020 19:24:18 +0200 Subject: [PATCH 32/66] code format --- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index 088c8f9ac8884..3f714e33ee80e 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -160,7 +160,6 @@ void PatternRecognitionbyCA::makeTracksters( uniqueLayerIds.erase(std::unique(uniqueLayerIds.begin(), uniqueLayerIds.end()), uniqueLayerIds.end()); unsigned int numberOfLayersInTrackster = uniqueLayerIds.size(); if (check_missing_layers_) { - int numberOfMissingLayers = 0; unsigned int j = showerMinLayerId; unsigned int indexInVec = 0; From 01a0f3e50f81c229ce2014e000d03b0d3295e55a Mon Sep 17 00:00:00 2001 From: Felice Date: Tue, 13 Oct 2020 21:10:55 +0200 Subject: [PATCH 33/66] Move trackster energy calculation after the selections are applied --- .../TICL/plugins/PatternRecognitionbyCA.cc | 68 +++++++++++++++---- .../TICL/plugins/PatternRecognitionbyCA.h | 3 + RecoHGCal/TICL/plugins/TrackstersProducer.cc | 26 ------- 3 files changed, 58 insertions(+), 39 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index 3f714e33ee80e..2c3c0869ac694 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -31,6 +31,9 @@ PatternRecognitionbyCA::PatternRecognitionbyCA(const edm::ParameterSet &c max_missing_layers_in_trackster_(conf.getParameter("max_missing_layers_in_trackster")), shower_start_max_layer_(conf.getParameter("shower_start_max_layer")), min_layers_per_trackster_(conf.getParameter("min_layers_per_trackster")), + filter_on_categories_(conf.getParameter>("filter_on_categories")), + pid_threshold_(conf.getParameter("pid_threshold")), + energy_em_over_total_threshold_(conf.getParameter("energy_em_over_total_threshold")), min_clusters_per_ntuplet_(min_layers_per_trackster_), max_delta_time_(conf.getParameter("max_delta_time")), eidInputName_(conf.getParameter("eid_input_name")), @@ -104,6 +107,10 @@ void PatternRecognitionbyCA::makeTracksters( const auto &doublets = theGraph_->getAllDoublets(); int tracksterId = 0; + // container for holding tracksters before selection + std::vector tmpTracksters; + tmpTracksters.reserve(foundNtuplets.size()); + for (auto const &ntuplet : foundNtuplets) { std::set effective_cluster_idx; std::pair::iterator, bool> retVal; @@ -181,35 +188,71 @@ void PatternRecognitionbyCA::makeTracksters( j++; } } + bool selected = (numberOfLayersInTrackster >= min_layers_per_trackster_) and (showerMinLayerId <= shower_start_max_layer_); if (selected) { - for (auto const i : effective_cluster_idx) { - layer_cluster_usage[i]++; - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) - LogDebug("HGCPatternRecoByCA") << "LayerID: " << i << " count: " << (int)layer_cluster_usage[i] << std::endl; - } // Put back indices, in the form of a Trackster, into the results vector Trackster tmp; tmp.vertices().reserve(effective_cluster_idx.size()); - tmp.vertex_multiplicity().resize(effective_cluster_idx.size(), 0); + tmp.vertex_multiplicity().resize(effective_cluster_idx.size(), 1); //regions and seedIndices can have different size //if a seeding region does not lead to any trackster - tmp.setSeed(input.regions[0].collectionID, seedIndices[tracksterId]); - if (isRegionalIter) { - seedToTracksterAssociation[tmp.seedIndex()].push_back(tracksterId); - } std::pair timeTrackster(-99., -1.); hgcalsimclustertime::ComputeClusterTime timeEstimator; timeTrackster = timeEstimator.fixSizeHighestDensity(times, timeErrors); tmp.setTimeAndError(timeTrackster.first, timeTrackster.second); std::copy(std::begin(effective_cluster_idx), std::end(effective_cluster_idx), std::back_inserter(tmp.vertices())); - result.push_back(tmp); + tmpTracksters.push_back(tmp); tracksterId++; } } + ticl::assignPCAtoTracksters(tmpTracksters, input.layerClusters, rhtools_.getPositionLayer(rhtools_.lastLayerEE(type)).z()); + + // run energy regression and ID + energyRegressionAndID(input.layerClusters, tmpTracksters); + // Filter results based on PID criteria or EM/Total energy ratio. + // We want to **keep** tracksters whose cumulative + // probability summed up over the selected categories + // is greater than the chosen threshold. Therefore + // the filtering function should **discard** all + // tracksters **below** the threshold. + auto filter_on_pids = [&](Trackster& t) -> bool { + auto cumulative_prob = 0.; + for (auto index : filter_on_categories_) { + cumulative_prob += t.id_probabilities(index); + } + return (cumulative_prob <= pid_threshold_) && + (t.raw_em_energy() / t.raw_energy() < energy_em_over_total_threshold_); + }; + + std::vector selectedTrackstersIds; + for(unsigned i = 0; i< tmpTracksters.size(); ++i) + { + if(!filter_on_pids(tmpTracksters[i])){ + selectedTrackstersIds.push_back(i); + } + } + + result.reserve(selectedTrackstersIds.size()); + + for(unsigned i = 0; i::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) + LogDebug("HGCPatternRecoByCA") << "LayerID: " << lcId << " count: " << (int)layer_cluster_usage[lcId] << std::endl; + } + t.setSeed(input.regions[0].collectionID, seedIndices[i]); + if (isRegionalIter) { + seedToTracksterAssociation[t.seedIndex()].push_back(i); + } + result.push_back(t); + } + for (auto &trackster : result) { assert(trackster.vertices().size() <= trackster.vertex_multiplicity().size()); @@ -220,7 +263,6 @@ void PatternRecognitionbyCA::makeTracksters( << " count: " << (int)trackster.vertex_multiplicity(i) << std::endl; } } - // Now decide if the tracksters from the track-based iterations have to be merged if (oneTracksterPerTrackSeed_) { std::vector tmp; @@ -228,8 +270,8 @@ void PatternRecognitionbyCA::makeTracksters( tmp.swap(result); } - ticl::assignPCAtoTracksters(result, input.layerClusters, rhtools_.getPositionLayer(rhtools_.lastLayerEE(type)).z()); + ticl::assignPCAtoTracksters(result, input.layerClusters, rhtools_.getPositionLayer(rhtools_.lastLayerEE(type)).z()); // run energy regression and ID energyRegressionAndID(input.layerClusters, result); diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h index d67526585bd74..1e42d65339ccf 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h @@ -43,6 +43,9 @@ namespace ticl { bool check_missing_layers_ = false; const unsigned int shower_start_max_layer_; const unsigned int min_layers_per_trackster_; + const std::vector filter_on_categories_; + const double pid_threshold_; + const double energy_em_over_total_threshold_; const int min_clusters_per_ntuplet_; const float max_delta_time_; const std::string eidInputName_; diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc index 5aa97b4a9b0f6..8c260b9b35a45 100644 --- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc @@ -54,11 +54,6 @@ class TrackstersProducer : public edm::stream::EDProducer layer_clusters_tiles_token_; edm::EDGetTokenT layer_clusters_tiles_hfnose_token_; const edm::EDGetTokenT> seeding_regions_token_; - const std::vector filter_on_categories_; - const double pid_threshold_; - const double energy_em_over_total_threshold_; - const int shower_start_max_layer_; - const std::string itername_; }; DEFINE_FWK_MODULE(TrackstersProducer); @@ -94,10 +89,6 @@ TrackstersProducer::TrackstersProducer(const edm::ParameterSet& ps, const Tracks consumes>>(ps.getParameter("time_layerclusters"))), seeding_regions_token_( consumes>(ps.getParameter("seeding_regions"))), - filter_on_categories_(ps.getParameter>("filter_on_categories")), - pid_threshold_(ps.getParameter("pid_threshold")), - energy_em_over_total_threshold_(ps.getParameter("energy_em_over_total_threshold")), - shower_start_max_layer_(ps.getParameter("shower_start_max_layer")), itername_(ps.getParameter("itername")) { if (doNose_) { layer_clusters_tiles_hfnose_token_ = @@ -180,30 +171,13 @@ void TrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) { myAlgo_->makeTracksters(input, *result, seedToTrackstersAssociation); } - // Now update the global mask and put it into the event output_mask->reserve(original_layerclusters_mask.size()); // Copy over the previous state std::copy( std::begin(original_layerclusters_mask), std::end(original_layerclusters_mask), std::back_inserter(*output_mask)); - // Filter results based on PID criteria or EM/Total energy ratio. - // We want to **keep** tracksters whose cumulative - // probability summed up over the selected categories - // is greater than the chosen threshold. Therefore - // the filtering function should **discard** all - // tracksters **below** the threshold. - auto filter_on_pids = [&](Trackster& t) -> bool { - auto cumulative_prob = 0.; - for (auto index : filter_on_categories_) { - cumulative_prob += t.id_probabilities(index); - } - return (cumulative_prob <= pid_threshold_) && - (t.raw_em_energy() / t.raw_energy() < energy_em_over_total_threshold_); - }; - // Actually filter results and shrink size to fit - result->erase(std::remove_if(result->begin(), result->end(), filter_on_pids), result->end()); // Mask the used elements, accordingly for (auto const& trackster : *result) { From 86ecf556e823ab6314ff603fb694a33e96d02589 Mon Sep 17 00:00:00 2001 From: Felice Date: Tue, 13 Oct 2020 21:21:08 +0200 Subject: [PATCH 34/66] Add a new max_longitudinal_sigmaPCA cut for EM iterations, with an initial value of 10 --- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc | 3 ++- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h | 1 + RecoHGCal/TICL/plugins/TrackstersProducer.cc | 1 + RecoHGCal/TICL/python/EMStep_cff.py | 1 + RecoHGCal/TICL/python/TrkEMStep_cff.py | 1 + 5 files changed, 6 insertions(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index 2c3c0869ac694..a62507130e206 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -34,6 +34,7 @@ PatternRecognitionbyCA::PatternRecognitionbyCA(const edm::ParameterSet &c filter_on_categories_(conf.getParameter>("filter_on_categories")), pid_threshold_(conf.getParameter("pid_threshold")), energy_em_over_total_threshold_(conf.getParameter("energy_em_over_total_threshold")), + max_longitudinal_sigmaPCA_(conf.getParameter("max_longitudinal_sigmaPCA")), min_clusters_per_ntuplet_(min_layers_per_trackster_), max_delta_time_(conf.getParameter("max_delta_time")), eidInputName_(conf.getParameter("eid_input_name")), @@ -231,7 +232,7 @@ void PatternRecognitionbyCA::makeTracksters( std::vector selectedTrackstersIds; for(unsigned i = 0; i< tmpTracksters.size(); ++i) { - if(!filter_on_pids(tmpTracksters[i])){ + if(!filter_on_pids(tmpTracksters[i]) and tmpTracksters[i].sigmasPCA()[0] < max_longitudinal_sigmaPCA_){ selectedTrackstersIds.push_back(i); } } diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h index 1e42d65339ccf..5a0229d5aa1e4 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h @@ -46,6 +46,7 @@ namespace ticl { const std::vector filter_on_categories_; const double pid_threshold_; const double energy_em_over_total_threshold_; + const double max_longitudinal_sigmaPCA_; const int min_clusters_per_ntuplet_; const float max_delta_time_; const std::string eidInputName_; diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc index 8c260b9b35a45..3729b9ec80cad 100644 --- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc @@ -114,6 +114,7 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri desc.add>("filter_on_categories", {0}); desc.add("pid_threshold", 0.); // make default such that no filtering is applied desc.add("energy_em_over_total_threshold", -1.); // make default such that no filtering is applied + desc.add("max_longitudinal_sigmaPCA", 9999); desc.add("shower_start_max_layer", 9999); // make default such that no filtering is applied desc.add("algo_verbosity", 0); desc.add("min_cos_theta", 0.915); diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 0346c7d762639..0aa6f9ca50fd9 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -25,6 +25,7 @@ filter_on_categories = [0, 1], pid_threshold = 0.5, energy_em_over_total_threshold = 0.9, + max_longitudinal_sigmaPCA = 10, shower_start_max_layer = 5, #inclusive max_out_in_hops = 1, skip_layers = 2, diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index 67b5d91d73ae9..a2a38163512ee 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -23,6 +23,7 @@ filter_on_categories = [0, 1], pid_threshold = 0.5, energy_em_over_total_threshold = 0.9, + max_longitudinal_sigmaPCA = 10, shower_start_max_layer = 5, #inclusive max_out_in_hops = 1, skip_layers = 2, From f497df561dfea97a64ec6d1c79464795d9934004 Mon Sep 17 00:00:00 2001 From: Felice Date: Tue, 13 Oct 2020 21:50:40 +0200 Subject: [PATCH 35/66] Add a new cut: root_doublet_max_distance_from_seed for TrkEM iteration, initial value 0.05 --- RecoHGCal/TICL/plugins/HGCGraph.cc | 28 +++++++++++++++---- RecoHGCal/TICL/plugins/HGCGraph.h | 1 + .../TICL/plugins/PatternRecognitionbyCA.cc | 2 ++ .../TICL/plugins/PatternRecognitionbyCA.h | 1 + RecoHGCal/TICL/plugins/TrackstersProducer.cc | 1 + RecoHGCal/TICL/python/TrkEMStep_cff.py | 1 + 6 files changed, 28 insertions(+), 6 deletions(-) diff --git a/RecoHGCal/TICL/plugins/HGCGraph.cc b/RecoHGCal/TICL/plugins/HGCGraph.cc index 1be81fb9a6f1b..892cfb8a277b1 100644 --- a/RecoHGCal/TICL/plugins/HGCGraph.cc +++ b/RecoHGCal/TICL/plugins/HGCGraph.cc @@ -19,6 +19,7 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, int deltaIPhi, float minCosTheta, float minCosPointing, + float root_doublet_max_distance_from_seed, float etaLimitIncreaseWindow, int skip_layers, int maxNumberOfLayers, @@ -27,6 +28,9 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, isOuterClusterOfDoublets_.resize(layerClusters.size()); allDoublets_.clear(); theRootDoublets_.clear(); + bool checkDistanceRootDoubletVsSeed = root_doublet_max_distance_from_seed < 9999; + float origin_eta; + float origin_phi; for (const auto &r : regions) { bool isGlobal = (r.index == -1); auto zSide = r.zSide; @@ -37,19 +41,22 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, startPhiBin = 0; endEtaBin = nEtaBins; endPhiBin = nPhiBins; + origin_eta = 0; + origin_phi = 0; } else { auto firstLayerOnZSide = maxNumberOfLayers * zSide; const auto &firstLayerHisto = histo[firstLayerOnZSide]; - - int entryEtaBin = firstLayerHisto.etaBin(r.origin.eta()); - int entryPhiBin = firstLayerHisto.phiBin(r.origin.phi()); + origin_eta = r.origin.eta(); + origin_phi = r.origin.phi(); + int entryEtaBin = firstLayerHisto.etaBin(origin_eta); + int entryPhiBin = firstLayerHisto.phiBin(origin_phi); // For track-seeded iterations, if the impact point is below a certain // eta-threshold, i.e., it has higher eta, make the initial search // window bigger in both eta and phi by one bin, to contain better low // energy showers. auto etaWindow = deltaIEta; auto phiWindow = deltaIPhi; - if (std::abs(r.origin.eta()) > etaLimitIncreaseWindow) { + if (std::abs(origin_eta) > etaLimitIncreaseWindow) { etaWindow++; phiWindow++; LogDebug("HGCGraph") << "Limit of Eta for increase: " << etaLimitIncreaseWindow @@ -60,9 +67,9 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, startPhiBin = entryPhiBin - phiWindow; endPhiBin = entryPhiBin + phiWindow + 1; if (verbosity_ > Guru) { - LogDebug("HGCGraph") << " Entrance eta, phi: " << r.origin.eta() << ", " << r.origin.phi() + LogDebug("HGCGraph") << " Entrance eta, phi: " << origin_eta << ", " << origin_phi << " entryEtaBin: " << entryEtaBin << " entryPhiBin: " << entryPhiBin - << " globalBin: " << firstLayerHisto.globalBin(r.origin.eta(), r.origin.phi()) + << " globalBin: " << firstLayerHisto.globalBin(origin_eta, origin_phi) << " on layer: " << firstLayerOnZSide << " startEtaBin: " << startEtaBin << " endEtaBin: " << endEtaBin << " startPhiBin: " << startPhiBin << " endPhiBin: " << endPhiBin << " phiBin(0): " << firstLayerHisto.phiBin(0.) @@ -171,8 +178,17 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, minCosTheta, minCosPointing, verbosity_ > Advanced); + if(isRootDoublet and checkDistanceRootDoubletVsSeed) + { + auto lcEta = layerClusters[innerClusterId].eta(); + auto lcPhi = layerClusters[innerClusterId].phi(); + if(std::hypot(lcEta-origin_eta, lcPhi-origin_phi) > root_doublet_max_distance_from_seed) + isRootDoublet = false; + } if (isRootDoublet) + { theRootDoublets_.push_back(doubletId); + } } } } diff --git a/RecoHGCal/TICL/plugins/HGCGraph.h b/RecoHGCal/TICL/plugins/HGCGraph.h index e9ee6da110821..836b024179c04 100644 --- a/RecoHGCal/TICL/plugins/HGCGraph.h +++ b/RecoHGCal/TICL/plugins/HGCGraph.h @@ -25,6 +25,7 @@ class HGCGraphT { int deltaIPhi, float minCosThetai, float maxCosPointing, + float root_doublet_max_distance_from_seed, float etaLimitIncreaseWindow, int skip_layers, int maxNumberOfLayers, diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index a62507130e206..b1c32f617b32a 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -26,6 +26,7 @@ PatternRecognitionbyCA::PatternRecognitionbyCA(const edm::ParameterSet &c max_out_in_hops_(conf.getParameter("max_out_in_hops")), min_cos_theta_(conf.getParameter("min_cos_theta")), min_cos_pointing_(conf.getParameter("min_cos_pointing")), + root_doublet_max_distance_from_seed_(conf.getParameter("root_doublet_max_distance_from_seed")), etaLimitIncreaseWindow_(conf.getParameter("etaLimitIncreaseWindow")), skip_layers_(conf.getParameter("skip_layers")), max_missing_layers_in_trackster_(conf.getParameter("max_missing_layers_in_trackster")), @@ -98,6 +99,7 @@ void PatternRecognitionbyCA::makeTracksters( 1, min_cos_theta_, min_cos_pointing_, + root_doublet_max_distance_from_seed_, etaLimitIncreaseWindow_, skip_layers_, rhtools_.lastLayer(type), diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h index 5a0229d5aa1e4..c96bf265e3f99 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h @@ -37,6 +37,7 @@ namespace ticl { const unsigned int max_out_in_hops_; const float min_cos_theta_; const float min_cos_pointing_; + const float root_doublet_max_distance_from_seed_; const float etaLimitIncreaseWindow_; const int skip_layers_; const int max_missing_layers_in_trackster_; diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc index 3729b9ec80cad..c9a15d6dc5bf1 100644 --- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc @@ -118,6 +118,7 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri desc.add("shower_start_max_layer", 9999); // make default such that no filtering is applied desc.add("algo_verbosity", 0); desc.add("min_cos_theta", 0.915); + desc.add("root_doublet_max_distance_from_seed", 9999); desc.add("min_cos_pointing", -1.); desc.add("skip_layers", 0); desc.add("max_missing_layers_in_trackster", 9999); diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index a2a38163512ee..55195ad60fd66 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -30,6 +30,7 @@ min_layers_per_trackster = 10, min_cos_theta = 0.97, # ~14 degrees min_cos_pointing = 0.94, # ~20 degrees + root_doublet_max_distance_from_seed = 0.05, max_delta_time = 3., itername = "TrkEM", algo_verbosity = 0, From 6d284044dfdf6967a70dbfbaf38aa8142eb575ba Mon Sep 17 00:00:00 2001 From: Felice Date: Tue, 13 Oct 2020 21:52:25 +0200 Subject: [PATCH 36/66] clang format --- RecoHGCal/TICL/plugins/HGCGraph.cc | 8 +++---- .../TICL/plugins/PatternRecognitionbyCA.cc | 23 ++++++++----------- RecoHGCal/TICL/plugins/TrackstersProducer.cc | 4 +--- 3 files changed, 14 insertions(+), 21 deletions(-) diff --git a/RecoHGCal/TICL/plugins/HGCGraph.cc b/RecoHGCal/TICL/plugins/HGCGraph.cc index 892cfb8a277b1..a69996468399a 100644 --- a/RecoHGCal/TICL/plugins/HGCGraph.cc +++ b/RecoHGCal/TICL/plugins/HGCGraph.cc @@ -178,15 +178,13 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, minCosTheta, minCosPointing, verbosity_ > Advanced); - if(isRootDoublet and checkDistanceRootDoubletVsSeed) - { + if (isRootDoublet and checkDistanceRootDoubletVsSeed) { auto lcEta = layerClusters[innerClusterId].eta(); auto lcPhi = layerClusters[innerClusterId].phi(); - if(std::hypot(lcEta-origin_eta, lcPhi-origin_phi) > root_doublet_max_distance_from_seed) + if (std::hypot(lcEta - origin_eta, lcPhi - origin_phi) > root_doublet_max_distance_from_seed) isRootDoublet = false; } - if (isRootDoublet) - { + if (isRootDoublet) { theRootDoublets_.push_back(doubletId); } } diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index b1c32f617b32a..f0bc70c1cad4e 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -191,7 +191,6 @@ void PatternRecognitionbyCA::makeTracksters( j++; } } - bool selected = (numberOfLayersInTrackster >= min_layers_per_trackster_) and (showerMinLayerId <= shower_start_max_layer_); @@ -212,17 +211,18 @@ void PatternRecognitionbyCA::makeTracksters( tracksterId++; } } - ticl::assignPCAtoTracksters(tmpTracksters, input.layerClusters, rhtools_.getPositionLayer(rhtools_.lastLayerEE(type)).z()); + ticl::assignPCAtoTracksters( + tmpTracksters, input.layerClusters, rhtools_.getPositionLayer(rhtools_.lastLayerEE(type)).z()); // run energy regression and ID energyRegressionAndID(input.layerClusters, tmpTracksters); - // Filter results based on PID criteria or EM/Total energy ratio. + // Filter results based on PID criteria or EM/Total energy ratio. // We want to **keep** tracksters whose cumulative // probability summed up over the selected categories // is greater than the chosen threshold. Therefore // the filtering function should **discard** all // tracksters **below** the threshold. - auto filter_on_pids = [&](Trackster& t) -> bool { + auto filter_on_pids = [&](Trackster &t) -> bool { auto cumulative_prob = 0.; for (auto index : filter_on_categories_) { cumulative_prob += t.id_probabilities(index); @@ -232,22 +232,21 @@ void PatternRecognitionbyCA::makeTracksters( }; std::vector selectedTrackstersIds; - for(unsigned i = 0; i< tmpTracksters.size(); ++i) - { - if(!filter_on_pids(tmpTracksters[i]) and tmpTracksters[i].sigmasPCA()[0] < max_longitudinal_sigmaPCA_){ + for (unsigned i = 0; i < tmpTracksters.size(); ++i) { + if (!filter_on_pids(tmpTracksters[i]) and tmpTracksters[i].sigmasPCA()[0] < max_longitudinal_sigmaPCA_) { selectedTrackstersIds.push_back(i); } } result.reserve(selectedTrackstersIds.size()); - for(unsigned i = 0; i::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) - LogDebug("HGCPatternRecoByCA") << "LayerID: " << lcId << " count: " << (int)layer_cluster_usage[lcId] << std::endl; + LogDebug("HGCPatternRecoByCA") << "LayerID: " << lcId << " count: " << (int)layer_cluster_usage[lcId] + << std::endl; } t.setSeed(input.regions[0].collectionID, seedIndices[i]); if (isRegionalIter) { @@ -256,7 +255,6 @@ void PatternRecognitionbyCA::makeTracksters( result.push_back(t); } - for (auto &trackster : result) { assert(trackster.vertices().size() <= trackster.vertex_multiplicity().size()); for (size_t i = 0; i < trackster.vertices().size(); ++i) { @@ -273,7 +271,6 @@ void PatternRecognitionbyCA::makeTracksters( tmp.swap(result); } - ticl::assignPCAtoTracksters(result, input.layerClusters, rhtools_.getPositionLayer(rhtools_.lastLayerEE(type)).z()); // run energy regression and ID energyRegressionAndID(input.layerClusters, result); diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc index c9a15d6dc5bf1..7202385174a30 100644 --- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc @@ -115,7 +115,7 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri desc.add("pid_threshold", 0.); // make default such that no filtering is applied desc.add("energy_em_over_total_threshold", -1.); // make default such that no filtering is applied desc.add("max_longitudinal_sigmaPCA", 9999); - desc.add("shower_start_max_layer", 9999); // make default such that no filtering is applied + desc.add("shower_start_max_layer", 9999); // make default such that no filtering is applied desc.add("algo_verbosity", 0); desc.add("min_cos_theta", 0.915); desc.add("root_doublet_max_distance_from_seed", 9999); @@ -179,8 +179,6 @@ void TrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) { std::copy( std::begin(original_layerclusters_mask), std::end(original_layerclusters_mask), std::back_inserter(*output_mask)); - - // Mask the used elements, accordingly for (auto const& trackster : *result) { for (auto const v : trackster.vertices()) { From 9638aa733eedab90832824cc25bd9a6e1cb35e99 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 14 Oct 2020 11:59:23 +0200 Subject: [PATCH 37/66] fix trackster to seed association --- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index f0bc70c1cad4e..ba2a9fad94716 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -108,13 +108,15 @@ void PatternRecognitionbyCA::makeTracksters( theGraph_->findNtuplets(foundNtuplets, seedIndices, min_clusters_per_ntuplet_, out_in_dfs_, max_out_in_hops_); //#ifdef FP_DEBUG const auto &doublets = theGraph_->getAllDoublets(); - int tracksterId = 0; + int tracksterId = -1; // container for holding tracksters before selection std::vector tmpTracksters; tmpTracksters.reserve(foundNtuplets.size()); for (auto const &ntuplet : foundNtuplets) { + tracksterId++; + std::set effective_cluster_idx; std::pair::iterator, bool> retVal; @@ -201,6 +203,7 @@ void PatternRecognitionbyCA::makeTracksters( tmp.vertex_multiplicity().resize(effective_cluster_idx.size(), 1); //regions and seedIndices can have different size //if a seeding region does not lead to any trackster + t.setSeed(input.regions[0].collectionID, seedIndices[i]); std::pair timeTrackster(-99., -1.); hgcalsimclustertime::ComputeClusterTime timeEstimator; @@ -208,7 +211,6 @@ void PatternRecognitionbyCA::makeTracksters( tmp.setTimeAndError(timeTrackster.first, timeTrackster.second); std::copy(std::begin(effective_cluster_idx), std::end(effective_cluster_idx), std::back_inserter(tmp.vertices())); tmpTracksters.push_back(tmp); - tracksterId++; } } ticl::assignPCAtoTracksters( @@ -248,7 +250,6 @@ void PatternRecognitionbyCA::makeTracksters( LogDebug("HGCPatternRecoByCA") << "LayerID: " << lcId << " count: " << (int)layer_cluster_usage[lcId] << std::endl; } - t.setSeed(input.regions[0].collectionID, seedIndices[i]); if (isRegionalIter) { seedToTracksterAssociation[t.seedIndex()].push_back(i); } From 4317f3ff75c6495714c2c8e46720e94ee7d421da Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 14 Oct 2020 12:04:07 +0200 Subject: [PATCH 38/66] Setting max_missing_layers_in_trackster to 2 from 0 in EM iterations --- RecoHGCal/TICL/python/EMStep_cff.py | 1 + RecoHGCal/TICL/python/TrkEMStep_cff.py | 1 + 2 files changed, 2 insertions(+) diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 0aa6f9ca50fd9..01cf65ad17f08 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -28,6 +28,7 @@ max_longitudinal_sigmaPCA = 10, shower_start_max_layer = 5, #inclusive max_out_in_hops = 1, + max_missing_layers_in_trackster = 2, skip_layers = 2, max_missing_layers_in_trackster = 1, min_layers_per_trackster = 10, diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index 55195ad60fd66..9f9c4c5f94fa1 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -26,6 +26,7 @@ max_longitudinal_sigmaPCA = 10, shower_start_max_layer = 5, #inclusive max_out_in_hops = 1, + max_missing_layers_in_trackster = 2, skip_layers = 2, min_layers_per_trackster = 10, min_cos_theta = 0.97, # ~14 degrees From 5bdb7707d121b7248a8d845549ada3cc969d5ea1 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Wed, 14 Oct 2020 12:29:33 +0200 Subject: [PATCH 39/66] Fix compilation error in C++ --- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index ba2a9fad94716..40abc16853eb9 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -203,7 +203,7 @@ void PatternRecognitionbyCA::makeTracksters( tmp.vertex_multiplicity().resize(effective_cluster_idx.size(), 1); //regions and seedIndices can have different size //if a seeding region does not lead to any trackster - t.setSeed(input.regions[0].collectionID, seedIndices[i]); + tmp.setSeed(input.regions[0].collectionID, seedIndices[tracksterId]); std::pair timeTrackster(-99., -1.); hgcalsimclustertime::ComputeClusterTime timeEstimator; From 57b342cb6ebfc6551a0cc43c9d39f63a21905011 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Wed, 14 Oct 2020 12:29:58 +0200 Subject: [PATCH 40/66] Remove duplicate python parameter --- RecoHGCal/TICL/python/EMStep_cff.py | 1 - 1 file changed, 1 deletion(-) diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 01cf65ad17f08..0aa6f9ca50fd9 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -28,7 +28,6 @@ max_longitudinal_sigmaPCA = 10, shower_start_max_layer = 5, #inclusive max_out_in_hops = 1, - max_missing_layers_in_trackster = 2, skip_layers = 2, max_missing_layers_in_trackster = 1, min_layers_per_trackster = 10, From efbafb6d96825d7cb2414435881002c963bb36b1 Mon Sep 17 00:00:00 2001 From: Felice Date: Fri, 16 Oct 2020 13:56:00 +0200 Subject: [PATCH 41/66] Sorting seeding regions by pT rather than p --- RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc b/RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc index a907a69c750a6..16d69c391f38a 100644 --- a/RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc +++ b/RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc @@ -61,7 +61,7 @@ void SeedingRegionByTracks::makeRegions(const edm::Event &ev, } // sorting seeding region by descending momentum std::sort(result.begin(), result.end(), [](const TICLSeedingRegion &a, const TICLSeedingRegion &b) { - return a.directionAtOrigin.mag2() > b.directionAtOrigin.mag2(); + return a.directionAtOrigin.perp2() > b.directionAtOrigin.perp2(); }); } From 8d6e738b6d3397f1747d90e5313c180bc431282f Mon Sep 17 00:00:00 2001 From: Felice Date: Fri, 16 Oct 2020 17:36:36 +0200 Subject: [PATCH 42/66] Added TICL PF interpretation. --- .../HGCalReco/interface/TICLLayerTile.h | 8 + .../TICL/plugins/TrackstersMergeProducer.cc | 486 ++++++++++++------ 2 files changed, 323 insertions(+), 171 deletions(-) diff --git a/DataFormats/HGCalReco/interface/TICLLayerTile.h b/DataFormats/HGCalReco/interface/TICLLayerTile.h index 83f07a3b18598..5d9cc8fb51a1c 100644 --- a/DataFormats/HGCalReco/interface/TICLLayerTile.h +++ b/DataFormats/HGCalReco/interface/TICLLayerTile.h @@ -33,6 +33,14 @@ class TICLLayerTileT { return phiBin; } + std::array searchBoxEtaPhi(float etaMin, float etaMax, float phiMin, float phiMax) const { + int etaBinMin = etaBin(etaMin); + int etaBinMax = etaBin(etaMax); + int phiBinMin = phiBin(phiMin); + int phiBinMax = phiBin(phiMax); + return std::array({{etaBinMin, etaBinMax, phiBinMin, phiBinMax}}); + } + int globalBin(int etaBin, int phiBin) const { return phiBin + etaBin * T::nPhiBins; } int globalBin(double eta, double phi) const { return phiBin(phi) + etaBin(eta) * T::nPhiBins; } diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index be103455d401b..4c7ec55367b8f 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -1,5 +1,5 @@ #include // unique_ptr - +#include // cosh #include "FWCore/Framework/interface/stream/EDProducer.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" @@ -15,6 +15,7 @@ #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" #include "RecoHGCal/TICL/plugins/GlobalCache.h" #include "PhysicsTools/TensorFlow/interface/TensorFlow.h" +#include "DataFormats/HGCalReco/interface/TICLCandidate.h" #include "TrackstersPCA.h" @@ -33,7 +34,7 @@ class TrackstersMergeProducer : public edm::stream::EDProducer &, TracksterIterIndex); @@ -53,11 +54,15 @@ class TrackstersMergeProducer : public edm::stream::EDProducer("phi_bin_window")), pt_sigma_high_(ps.getParameter("pt_sigma_high")), pt_sigma_low_(ps.getParameter("pt_sigma_low")), + halo_max_distance2_(ps.getParameter("halo_max_distance2")), + mergeTRK_max_dR_(ps.getParameter("mergeTRK_max_dR")), cosangle_align_(ps.getParameter("cosangle_align")), e_over_h_threshold_(ps.getParameter("e_over_h_threshold")), pt_neutral_threshold_(ps.getParameter("pt_neutral_threshold")), - resol_calo_offset_(ps.getParameter("resol_calo_offset")), - resol_calo_scale_(ps.getParameter("resol_calo_scale")), + resol_calo_offset_had_(ps.getParameter("resol_calo_offset_had")), + resol_calo_scale_had_(ps.getParameter("resol_calo_scale_had")), + resol_calo_offset_em_(ps.getParameter("resol_calo_offset_em")), + resol_calo_scale_em_(ps.getParameter("resol_calo_scale_em")), debug_(ps.getParameter("debug")), eidInputName_(ps.getParameter("eid_input_name")), eidOutputNameEnergy_(ps.getParameter("eid_output_name_energy")), @@ -108,6 +117,7 @@ TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps, co eidSession_ = tensorflow::createSession(trackstersCache->eidGraphDef); produces>(); + produces>(); } void TrackstersMergeProducer::fillTile(TICLTracksterTiles &tracksterTile, @@ -150,16 +160,20 @@ void TrackstersMergeProducer::dumpTrackster(const Trackster &t) const { void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es) { edm::ESHandle geom = es.getHandle(geometry_token_); rhtools_.setGeometry(*geom); - auto result = std::make_unique>(); - auto mergedTrackstersTRK = std::make_unique>(); + auto resultTrackstersMerged = std::make_unique>(); + auto resultCandidates = std::make_unique>(); TICLTracksterTiles tracksterTile; - std::vector usedTrackstersEM; - std::vector usedTrackstersTRKEM; - std::vector usedTrackstersTRK; - std::vector usedTrackstersHAD; + std::vector usedTrackstersMerged; + std::vector indexInMergedCollTRKEM; + std::vector indexInMergedCollEM; + std::vector indexInMergedCollTRK; + std::vector indexInMergedCollHAD; std::vector usedSeeds; + // associating seed to the index of the trackster in the merged collection and the iteration that found it + std::map>> seedToTracksterAssociator; + std::vector iterMergedTracksters; edm::Handle> track_h; evt.getByToken(tracks_token_, track_h); const auto &tracks = *track_h; @@ -171,205 +185,331 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es edm::Handle> trackstersem_h; evt.getByToken(trackstersem_token_, trackstersem_h); const auto &trackstersEM = *trackstersem_h; - usedTrackstersEM.resize(trackstersEM.size(), false); edm::Handle> tracksterstrkem_h; evt.getByToken(tracksterstrkem_token_, tracksterstrkem_h); const auto &trackstersTRKEM = *tracksterstrkem_h; - usedTrackstersTRKEM.resize(trackstersTRKEM.size(), false); edm::Handle> tracksterstrk_h; evt.getByToken(tracksterstrk_token_, tracksterstrk_h); const auto &trackstersTRK = *tracksterstrk_h; - usedTrackstersTRK.resize(trackstersTRK.size(), false); edm::Handle> trackstershad_h; evt.getByToken(trackstershad_token_, trackstershad_h); const auto &trackstersHAD = *trackstershad_h; - usedTrackstersHAD.resize(trackstersHAD.size(), false); edm::Handle> seedingTrk_h; evt.getByToken(seedingTrk_token_, seedingTrk_h); const auto &seedingTrk = *seedingTrk_h; - usedSeeds.resize(seedingTrk.size(), false); + usedSeeds.resize(tracks.size(), false); + fillTile(tracksterTile, trackstersTRKEM, TracksterIterIndex::TRKEM); fillTile(tracksterTile, trackstersEM, TracksterIterIndex::EM); fillTile(tracksterTile, trackstersTRK, TracksterIterIndex::TRK); fillTile(tracksterTile, trackstersHAD, TracksterIterIndex::HAD); - auto seedId = 0; - for (auto const &s : seedingTrk) { - tracksterTile.fill(TracksterIterIndex::SEED, s.origin.eta(), s.origin.phi(), seedId++); - - if (debug_) { - LogDebug("TrackstersMergeProducer") - << "Seed index: " << seedId << " internal index: " << s.index << " origin: " << s.origin - << " mom: " << s.directionAtOrigin << " pt: " << std::sqrt(s.directionAtOrigin.perp2()) - << " zSide: " << s.zSide << " collectionID: " << s.collectionID << " track pt " << tracks[s.index].pt() - << std::endl; - } - } + auto totalNumberOfTracksters = + trackstersTRKEM.size() + trackstersTRK.size() + trackstersEM.size() + trackstersHAD.size(); + resultTrackstersMerged->reserve(totalNumberOfTracksters); + iterMergedTracksters.reserve(totalNumberOfTracksters); + usedTrackstersMerged.resize(totalNumberOfTracksters, false); + indexInMergedCollTRKEM.reserve(trackstersTRKEM.size()); + indexInMergedCollEM.reserve(trackstersEM.size()); + indexInMergedCollTRK.reserve(trackstersTRK.size()); + indexInMergedCollHAD.reserve(trackstersHAD.size()); if (debug_) { - printTrackstersDebug(trackstersTRK, "tracksterTRK"); + printTrackstersDebug(trackstersTRKEM, "tracksterTRKEM"); printTrackstersDebug(trackstersEM, "tracksterEM"); + printTrackstersDebug(trackstersTRK, "tracksterTRK"); printTrackstersDebug(trackstersHAD, "tracksterHAD"); } - int tracksterTRK_idx = 0; - int tracksterHAD_idx = 0; - if (optimiseAcrossTracksters_) { - for (auto const &t : trackstersTRK) { - if (debug_) { - int entryEtaBin = tracksterTile[TracksterIterIndex::TRK].etaBin(t.barycenter().eta()); - int entryPhiBin = tracksterTile[TracksterIterIndex::TRK].phiBin(t.barycenter().phi()); - int bin = tracksterTile[TracksterIterIndex::TRK].globalBin(t.barycenter().eta(), t.barycenter().phi()); - LogDebug("TrackstersMergeProducer") - << "TrackstersMergeProducer Tracking obj: " << t.barycenter() << " in bin " << bin << " etaBin " - << entryEtaBin << " phiBin " << entryPhiBin << std::endl; - dumpTrackster(t); - } - auto const &track = tracks[t.seedIndex()]; - auto trk_pt = (float)track.pt(); - auto diff_pt = t.raw_pt() - trk_pt; - auto pt_err = trk_pt * resol_calo_scale_ + resol_calo_offset_; - auto w_cal = 1. / (pt_err * pt_err); - auto w_trk = 1. / (track.ptError() * track.ptError()); - auto diff_pt_sigmas = diff_pt / pt_err; - auto e_over_h = (t.raw_em_pt() / ((t.raw_pt() - t.raw_em_pt()) != 0. ? (t.raw_pt() - t.raw_em_pt()) : 1.)); - LogDebug("TrackstersMergeProducer") - << "trackster_pt: " << t.raw_pt() << std::endl - << "track pt (inner): " << track.pt() << std::endl - << "track eta (inner): " << track.eta() << std::endl - << "track _phi (inner): " << track.phi() << std::endl - << "track pt (outer): " << std::sqrt(track.outerMomentum().perp2()) << std::endl - << "track eta (outer): " << track.outerMomentum().eta() << std::endl - << "track _phi (outer): " << track.outerMomentum().phi() << std::endl - << "pt_err_track: " << track.ptError() << std::endl - << "diff_pt: " << diff_pt << std::endl - << "pt_err: " << pt_err << std::endl - << "diff_pt_sigmas: " << diff_pt_sigmas << std::endl - << "w_cal: " << w_cal << std::endl - << "w_trk: " << w_trk << std::endl - << "average_pt: " << (t.raw_pt() * w_cal + trk_pt * w_trk) / (w_cal + w_trk) << std::endl - << "e_over_h: " << e_over_h << std::endl; - - // If the energy is unbalanced and higher in Calo ==> balance it by - // emitting gammas/neutrals - if (diff_pt_sigmas > pt_sigma_high_) { - if (e_over_h > e_over_h_threshold_) { - auto gamma_pt = std::min(diff_pt, t.raw_em_pt()); - // Create gamma with calo direction - LogDebug("TrackstersMergeProducer") - << " Creating a photon from TRK Trackster with energy " << gamma_pt << " and direction " - << t.eigenvectors(0).eta() << ", " << t.eigenvectors(0).phi() << std::endl; - diff_pt -= gamma_pt; - if (diff_pt > pt_neutral_threshold_) { - // Create also a neutral on top, with calo direction and diff_pt as energy - LogDebug("TrackstersMergeProducer") - << " Adding also a neutral hadron from TRK Trackster with energy " << diff_pt << " and direction " - << t.eigenvectors(0).eta() << ", " << t.eigenvectors(0).phi() << std::endl; - } - } else { - // Create neutral with calo direction - LogDebug("TrackstersMergeProducer") - << " Creating a neutral hadron from TRK Trackster with energy " << diff_pt << " and direction " - << t.eigenvectors(0).eta() << ", " << t.eigenvectors(0).phi() << std::endl; - } - } - // If the energy is in the correct ball-park (this also covers the - // previous case after the neutral emission), use weighted averages - if (diff_pt_sigmas > -pt_sigma_low_) { - // Create either an electron or a charged hadron, using the weighted - // average information between the track and the cluster for the - // energy, the direction of the track at the vertex. The PID is simply - // passed over to the final trackster, while the energy is changed. - auto average_pt = (w_cal * t.raw_pt() + trk_pt * w_trk) / (w_cal + w_trk); - LogDebug("TrackstersMergeProducer") - << " Creating electron/charged hadron from TRK Trackster with weighted p_t " << average_pt - << " and direction " << track.eta() << ", " << track.phi() << std::endl; + for (auto const &t : trackstersTRKEM) { + indexInMergedCollTRKEM.push_back(resultTrackstersMerged->size()); + seedToTracksterAssociator[t.seedIndex()].emplace_back(resultTrackstersMerged->size(), TracksterIterIndex::TRKEM); + resultTrackstersMerged->push_back(t); + iterMergedTracksters.push_back(TracksterIterIndex::TRKEM); + } - } - // If the energy of the calo is too low, just use track-only information - else { - // Create either an electron or a charged hadron, using the track - // information only. - LogDebug("TrackstersMergeProducer") - << " Creating electron/charged hadron from TRK Trackster with track p_t " << trk_pt << " and direction " - << track.eta() << ", " << track.phi() << std::endl; - } - result->push_back(t); - usedTrackstersTRK[tracksterTRK_idx] = true; - tracksterTRK_idx++; - } + for (auto const &t : trackstersEM) { + indexInMergedCollEM.push_back(resultTrackstersMerged->size()); + resultTrackstersMerged->push_back(t); + iterMergedTracksters.push_back(TracksterIterIndex::EM); } - tracksterTRK_idx = 0; for (auto const &t : trackstersTRK) { - if (debug_) { - LogDebug("TrackstersMergeProducer") << " Considering trackster " << tracksterTRK_idx - << " as used: " << usedTrackstersTRK[tracksterTRK_idx] << std::endl; - } - if (!usedTrackstersTRK[tracksterTRK_idx]) { - LogDebug("TrackstersMergeProducer") - << " Creating a charge hadron from TRK Trackster with track energy " << t.raw_energy() << " and direction " - << t.eigenvectors(0).eta() << ", " << t.eigenvectors(0).phi() << std::endl; - result->push_back(t); - } - tracksterTRK_idx++; + indexInMergedCollTRK.push_back(resultTrackstersMerged->size()); + seedToTracksterAssociator[t.seedIndex()].emplace_back(resultTrackstersMerged->size(), TracksterIterIndex::TRK); + resultTrackstersMerged->push_back(t); + iterMergedTracksters.push_back(TracksterIterIndex::TRK); } - auto tracksterEM_idx = 0; - for (auto const &t : trackstersEM) { - if (debug_) { - LogDebug("TrackstersMergeProducer") << " Considering trackster " << tracksterEM_idx - << " as used: " << usedTrackstersEM[tracksterEM_idx] << std::endl; - } - if (!usedTrackstersEM[tracksterEM_idx]) { - LogDebug("TrackstersMergeProducer") - << " Creating a photon from EM Trackster with track energy " << t.raw_energy() << " and direction " - << t.eigenvectors(0).eta() << ", " << t.eigenvectors(0).phi() << std::endl; - result->push_back(t); - } - tracksterEM_idx++; + for (auto const &t : trackstersHAD) { + indexInMergedCollHAD.push_back(resultTrackstersMerged->size()); + resultTrackstersMerged->push_back(t); + iterMergedTracksters.push_back(TracksterIterIndex::HAD); } - tracksterHAD_idx = 0; - for (auto const &t : trackstersHAD) { - if (debug_) { - LogDebug("TrackstersMergeProducer") << " Considering trackster " << tracksterHAD_idx - << " as used: " << usedTrackstersHAD[tracksterHAD_idx] << std::endl; + assignPCAtoTracksters(*resultTrackstersMerged, layerClusters, rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z()); + energyRegressionAndID(layerClusters, *resultTrackstersMerged); + + printTrackstersDebug(*resultTrackstersMerged, "TrackstersMergeProducer"); + + auto trackstersMergedHandle = evt.put(std::move(resultTrackstersMerged)); + + // TICL Candidate creation + // We start from neutrals first + + // Photons + for (unsigned i = 0; i < trackstersEM.size(); ++i) { + auto mergedIdx = indexInMergedCollEM[i]; + usedTrackstersMerged[mergedIdx] = true; + auto &t = trackstersEM[i]; //trackster + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); + tmpCandidate.setCharge(0); + tmpCandidate.setPdgId(22); + tmpCandidate.setRawEnergy(t.raw_energy()); + math::XYZTLorentzVector p4(t.raw_energy() * t.barycenter().unit().x(), + t.raw_energy() * t.barycenter().unit().y(), + t.raw_energy() * t.barycenter().unit().z(), + t.raw_energy()); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + } + + // Neutral Hadrons + constexpr float mpion2 = 0.13957f * 0.13957f; + for (unsigned i = 0; i < trackstersHAD.size(); ++i) { + auto mergedIdx = indexInMergedCollHAD[i]; + usedTrackstersMerged[mergedIdx] = true; + auto &t = trackstersHAD[i]; //trackster + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); + tmpCandidate.setCharge(0); + tmpCandidate.setPdgId(130); + tmpCandidate.setRawEnergy(t.raw_energy()); + float momentum = std::sqrt(t.raw_energy() * t.raw_energy() - mpion2); + math::XYZTLorentzVector p4(momentum * t.barycenter().unit().x(), + momentum * t.barycenter().unit().y(), + momentum * t.barycenter().unit().z(), + t.raw_energy()); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + } + + // Charged Particles + for (unsigned i = 0; i < trackstersTRKEM.size(); ++i) { + auto mergedIdx = indexInMergedCollTRKEM[i]; + if (!usedTrackstersMerged[mergedIdx]) { + auto &t = trackstersTRKEM[i]; //trackster + auto trackIdx = t.seedIndex(); + auto const &track = tracks[trackIdx]; + if (!usedSeeds[trackIdx] and t.raw_energy() > 0) { + usedSeeds[trackIdx] = true; + usedTrackstersMerged[mergedIdx] = true; + + std::vector trackstersTRKwithSameSeed; + std::vector trackstersTRKEMwithSameSeed; + + for (const auto &tracksterIterationPair : seedToTracksterAssociator[trackIdx]) { + if (tracksterIterationPair.first != mergedIdx and !usedTrackstersMerged[tracksterIterationPair.first] and + trackstersMergedHandle->at(tracksterIterationPair.first).raw_energy() > 0.) { + if (tracksterIterationPair.second == TracksterIterIndex::TRKEM) { + trackstersTRKEMwithSameSeed.push_back(tracksterIterationPair.first); + } else if (tracksterIterationPair.second == TracksterIterIndex::TRK) { + trackstersTRKwithSameSeed.push_back(tracksterIterationPair.first); + } + } + } + + float tracksterTotalRawPt = t.raw_pt(); + std::vector haloTrackstersTRKIdx; + bool foundCompatibleTRK = false; + + for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { + usedTrackstersMerged[otherTracksterIdx] = true; + tracksterTotalRawPt += trackstersMergedHandle->at(otherTracksterIdx).raw_pt(); + + // Check the X,Y,Z barycenter and merge if they are very close (halo) + if ((t.barycenter() - trackstersMergedHandle->at(otherTracksterIdx).barycenter()).mag2() < + halo_max_distance2_) { + haloTrackstersTRKIdx.push_back(otherTracksterIdx); + + } else { + foundCompatibleTRK = true; + } + } + + //check if there is 1-to-1 relationship + if (trackstersTRKEMwithSameSeed.empty()) { + if (foundCompatibleTRK) { + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); + tmpCandidate.setPdgId(211 * track.charge()); + for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); + } + float p = tracksterTotalRawPt * cosh(t.barycenter().eta()); + float energy = std::sqrt(p * p + mpion2); + tmpCandidate.setRawEnergy(energy); + math::XYZTLorentzVector p4(p * track.momentum().unit().x(), + p * track.momentum().unit().y(), + p * track.momentum().unit().z(), + energy); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + + } else { + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); + for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); + } + tmpCandidate.setPdgId(11 * track.charge()); + float p = tracksterTotalRawPt * cosh(t.barycenter().eta()); + tmpCandidate.setRawEnergy(p); + math::XYZTLorentzVector p4( + p * track.momentum().unit().x(), p * track.momentum().unit().y(), p * track.momentum().unit().z(), p); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + } + + } else { + // if 1-to-many find closest trackster in momentum + int closestTrackster = mergedIdx; + float minPtDiff = std::abs(t.raw_pt() - track.pt()); + for (auto otherTracksterIdx : trackstersTRKEMwithSameSeed) { + auto thisPt = tracksterTotalRawPt + trackstersMergedHandle->at(otherTracksterIdx).raw_pt() - t.raw_pt(); + closestTrackster = std::abs(thisPt - track.pt()) < minPtDiff ? otherTracksterIdx : closestTrackster; + usedTrackstersMerged[otherTracksterIdx] = true; + } + tracksterTotalRawPt += trackstersMergedHandle->at(closestTrackster).raw_pt() - t.raw_pt(); + + if (foundCompatibleTRK) { + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, closestTrackster)); + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); + tmpCandidate.setPdgId(211 * track.charge()); + for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); + } + float p = tracksterTotalRawPt * cosh(t.barycenter().eta()); + float energy = std::sqrt(p * p + mpion2); + tmpCandidate.setRawEnergy(energy); + math::XYZTLorentzVector p4(p * track.momentum().unit().x(), + p * track.momentum().unit().y(), + p * track.momentum().unit().z(), + energy); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + + } else { + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, closestTrackster)); + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); + for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); + } + tmpCandidate.setPdgId(11 * track.charge()); + float p = tracksterTotalRawPt * cosh(t.barycenter().eta()); + tmpCandidate.setRawEnergy(p); + math::XYZTLorentzVector p4( + p * track.momentum().unit().x(), p * track.momentum().unit().y(), p * track.momentum().unit().z(), p); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + } + // Promote all other TRKEM tracksters as photons with their energy. + for (auto otherTracksterIdx : trackstersTRKEMwithSameSeed) { + auto tmpIndex = (otherTracksterIdx != closestTrackster) ? otherTracksterIdx : mergedIdx; + TICLCandidate photonCandidate; + auto &otherTrackster = trackstersMergedHandle->at(tmpIndex); + auto gammaEnergy = otherTrackster.raw_energy(); + photonCandidate.setCharge(0); + photonCandidate.setPdgId(22); + photonCandidate.setRawEnergy(gammaEnergy); + math::XYZTLorentzVector gammaP4(gammaEnergy * otherTrackster.barycenter().unit().x(), + gammaEnergy * otherTrackster.barycenter().unit().y(), + gammaEnergy * otherTrackster.barycenter().unit().z(), + gammaEnergy); + photonCandidate.setP4(gammaP4); + photonCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, tmpIndex)); + resultCandidates->push_back(photonCandidate); + } + } + } } - if (!usedTrackstersHAD[tracksterHAD_idx]) { - LogDebug("TrackstersMergeProducer") - << " Creating a neutral hadron from HAD Trackster with track energy " << t.raw_energy() << " and direction " - << t.eigenvectors(0).eta() << ", " << t.eigenvectors(0).phi() << std::endl; - result->push_back(t); + } //end of loop over trackstersTRKEM + + for (unsigned i = 0; i < trackstersTRK.size(); ++i) { + auto mergedIdx = indexInMergedCollTRK[i]; + auto &t = trackstersTRK[i]; //trackster + + if (!usedTrackstersMerged[mergedIdx] and t.raw_energy() > 0) { + auto trackIdx = t.seedIndex(); + auto const &track = tracks[trackIdx]; + if (!usedSeeds[trackIdx]) { + usedSeeds[trackIdx] = true; + usedTrackstersMerged[mergedIdx] = true; + float p = t.raw_pt() * cosh(t.barycenter().eta()); + float energy = std::sqrt(p * p + mpion2); + TICLCandidate tmpCandidate; + tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); + tmpCandidate.setPdgId(211 * track.charge()); + tmpCandidate.setRawEnergy(energy); + math::XYZTLorentzVector p4( + p * track.momentum().unit().x(), p * track.momentum().unit().y(), p * track.momentum().unit().z(), energy); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + } } - tracksterHAD_idx++; } - - auto tracksterTRKEM_idx = 0; - for (auto const &t : trackstersTRKEM) { - if (debug_) { - LogDebug("TrackstersMergeProducer") << " Considering trackster " << tracksterTRKEM_idx - << " as used: " << usedTrackstersTRKEM[tracksterTRKEM_idx] << std::endl; + // For all seeds that have 0-energy tracksters whose track is not marked as used, create a charged hadron with the track information. + for (auto const &s : seedingTrk) { + if (usedSeeds[s.index] == false) { + auto const &track = tracks[s.index]; + // emit a charged hadron + TICLCandidate tmpCandidate; + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, s.index)); + tmpCandidate.setPdgId(211 * track.charge()); + float energy = std::sqrt(track.pt() * track.pt() + mpion2); + tmpCandidate.setRawEnergy(energy); + math::XYZTLorentzVector p4(track.momentum().x(), track.momentum().y(), track.momentum().z(), energy); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + usedSeeds[s.index] = true; } - if (!usedTrackstersTRKEM[tracksterTRKEM_idx]) { - LogDebug("TrackstersMergeProducer") - << " Creating a photon from EM Trackster with track energy " << t.raw_energy() << " and direction " - << t.eigenvectors(0).eta() << ", " << t.eigenvectors(0).phi() << std::endl; - result->push_back(t); + } + // for all general tracks (high purity, pt > 1), check if they have been used: if not, promote them as charged hadrons + for (unsigned i = 0; i < tracks.size(); ++i) { + auto const &track = tracks[i]; + if (track.pt() > 1.f and track.quality(reco::TrackBase::highPurity) and track.missingOuterHits() < 5 and + std::abs(track.outerEta()) > 1.48 and std::abs(track.outerEta()) < 3.0 and usedSeeds[i] == false) { + // emit a charged hadron + TICLCandidate tmpCandidate; + tmpCandidate.setCharge(track.charge()); + tmpCandidate.setTrackPtr(edm::Ptr(track_h, i)); + tmpCandidate.setPdgId(211 * track.charge()); + float energy = std::sqrt(track.pt() * track.pt() + mpion2); + tmpCandidate.setRawEnergy(energy); + math::XYZTLorentzVector p4(track.momentum().x(), track.momentum().y(), track.momentum().z(), energy); + tmpCandidate.setP4(p4); + resultCandidates->push_back(tmpCandidate); + usedSeeds[i] = true; } - tracksterTRKEM_idx++; } - assignPCAtoTracksters(*result, layerClusters, rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z()); - energyRegressionAndID(layerClusters, *result); - - printTrackstersDebug(*result, "TrackstersMergeProducer"); - - evt.put(std::move(result)); + evt.put(std::move(resultCandidates)); } void TrackstersMergeProducer::energyRegressionAndID(const std::vector &layerClusters, @@ -573,11 +713,15 @@ void TrackstersMergeProducer::fillDescriptions(edm::ConfigurationDescriptions &d desc.add("phi_bin_window", 1); desc.add("pt_sigma_high", 2.); desc.add("pt_sigma_low", 2.); + desc.add("halo_max_distance2", 4.); + desc.add("mergeTRK_max_dR", 0.03); desc.add("cosangle_align", 0.9945); desc.add("e_over_h_threshold", 1.); desc.add("pt_neutral_threshold", 2.); - desc.add("resol_calo_offset", 1.5); - desc.add("resol_calo_scale", 0.15); + desc.add("resol_calo_offset_had", 1.5); + desc.add("resol_calo_scale_had", 0.15); + desc.add("resol_calo_offset_em", 1.5); + desc.add("resol_calo_scale_em", 0.15); desc.add("debug", true); desc.add("eid_graph_path", "RecoHGCal/TICL/data/tf_models/energy_id_v0.pb"); desc.add("eid_input_name", "input"); From 20adab600042e56410f5e8989e3f5dd6edb0e94a Mon Sep 17 00:00:00 2001 From: Felice Date: Mon, 19 Oct 2020 17:46:58 +0200 Subject: [PATCH 43/66] Feed PFTICL the correct InputTag --- RecoHGCal/TICL/plugins/PFTICLProducer.cc | 2 +- RecoHGCal/TICL/python/iterativeTICL_cff.py | 8 ++------ 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PFTICLProducer.cc b/RecoHGCal/TICL/plugins/PFTICLProducer.cc index 443340944b39b..74d69a350b7b4 100644 --- a/RecoHGCal/TICL/plugins/PFTICLProducer.cc +++ b/RecoHGCal/TICL/plugins/PFTICLProducer.cc @@ -36,7 +36,7 @@ PFTICLProducer::PFTICLProducer(const edm::ParameterSet& conf) void PFTICLProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; - desc.add("ticlCandidateSrc", edm::InputTag("ticlCandidateFromTracksters")); + desc.add("ticlCandidateSrc", edm::InputTag("ticlTrackstersMerge")); descriptions.add("pfTICLProducer", desc); } diff --git a/RecoHGCal/TICL/python/iterativeTICL_cff.py b/RecoHGCal/TICL/python/iterativeTICL_cff.py index b6a9b8048629a..c6085a16e0370 100644 --- a/RecoHGCal/TICL/python/iterativeTICL_cff.py +++ b/RecoHGCal/TICL/python/iterativeTICL_cff.py @@ -19,13 +19,9 @@ ) ticlTracksterMergeTask = cms.Task(ticlTrackstersMerge, ticlMultiClustersFromTrackstersMerge) -ticlCandidateFromTracksters = _ticlCandidateFromTrackstersProducer.clone( - tracksterCollections = ["ticlTrackstersMerge"], - # A possible alternative for momentum computation: - # momentumPlugin = dict(plugin="TracksterP4FromTrackAndPCA") - ) + pfTICL = _pfTICLProducer.clone() -ticlPFTask = cms.Task(ticlCandidateFromTracksters, pfTICL) +ticlPFTask = cms.Task(pfTICL) iterTICLTask = cms.Task(ticlLayerTileTask ,ticlTrkEMStepTask From 6c32c60c232982470d20b398ee1b203bd5d6f37a Mon Sep 17 00:00:00 2001 From: Felice Date: Mon, 19 Oct 2020 17:48:49 +0200 Subject: [PATCH 44/66] Adapt event content --- RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py index 5b0ecc4d4b9c9..77b647e6ce291 100644 --- a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py +++ b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py @@ -22,11 +22,9 @@ 'keep *_ticlTrackstersTrk_*_*', 'keep *_ticlTrackstersMIP_*_*', 'keep *_ticlTrackstersMerge_*_*', - 'keep *_ticlCandidateFromTracksters_*_*', 'keep *_ticlTrackstersHFNoseEM_*_*', 'keep *_ticlTrackstersHFNoseMIP_*_*', 'keep *_ticlTrackstersHFNoseMerge_*_*', - 'keep *_ticlCandidateFromTrackstersHFNose_*_*', 'keep *_pfTICL_*_*' ) ) From d3a76b32036f58fe01c89e32f0d1d983669625d4 Mon Sep 17 00:00:00 2001 From: Felice Date: Fri, 23 Oct 2020 14:30:40 +0200 Subject: [PATCH 45/66] Address review comments --- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc | 13 ++++--------- RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc | 8 ++++---- .../HGCalValidation/src/HGVHistoProducerAlgo.cc | 4 ++-- 3 files changed, 10 insertions(+), 15 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index 40abc16853eb9..10ba4b0cf1e87 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -30,6 +30,7 @@ PatternRecognitionbyCA::PatternRecognitionbyCA(const edm::ParameterSet &c etaLimitIncreaseWindow_(conf.getParameter("etaLimitIncreaseWindow")), skip_layers_(conf.getParameter("skip_layers")), max_missing_layers_in_trackster_(conf.getParameter("max_missing_layers_in_trackster")), + check_missing_layers_(max_missing_layers_in_trackster_ < 100), shower_start_max_layer_(conf.getParameter("shower_start_max_layer")), min_layers_per_trackster_(conf.getParameter("min_layers_per_trackster")), filter_on_categories_(conf.getParameter>("filter_on_categories")), @@ -52,9 +53,6 @@ PatternRecognitionbyCA::PatternRecognitionbyCA(const edm::ParameterSet &c << "PatternRecognitionbyCA received an empty graph definition from the global cache"; } eidSession_ = tensorflow::createSession(trackstersCache->eidGraphDef); - if (max_missing_layers_in_trackster_ < 100) { - check_missing_layers_ = true; - } } template @@ -156,7 +154,6 @@ void PatternRecognitionbyCA::makeTracksters( } } unsigned showerMinLayerId = 99999; - std::vector layerIds; std::vector uniqueLayerIds; uniqueLayerIds.reserve(effective_cluster_idx.size()); std::vector> lcIdAndLayer; @@ -175,7 +172,7 @@ void PatternRecognitionbyCA::makeTracksters( int numberOfMissingLayers = 0; unsigned int j = showerMinLayerId; unsigned int indexInVec = 0; - for (auto &layer : uniqueLayerIds) { + for (const auto &layer : uniqueLayerIds) { if (layer != j) { numberOfMissingLayers++; j++; @@ -194,9 +191,7 @@ void PatternRecognitionbyCA::makeTracksters( } } - bool selected = - (numberOfLayersInTrackster >= min_layers_per_trackster_) and (showerMinLayerId <= shower_start_max_layer_); - if (selected) { + if ((numberOfLayersInTrackster >= min_layers_per_trackster_) and (showerMinLayerId <= shower_start_max_layer_)) { // Put back indices, in the form of a Trackster, into the results vector Trackster tmp; tmp.vertices().reserve(effective_cluster_idx.size()); @@ -243,7 +238,7 @@ void PatternRecognitionbyCA::makeTracksters( result.reserve(selectedTrackstersIds.size()); for (unsigned i = 0; i < selectedTrackstersIds.size(); ++i) { - auto &t = tmpTracksters[selectedTrackstersIds[i]]; + const auto &t = tmpTracksters[selectedTrackstersIds[i]]; for (auto const lcId : t.vertices()) { layer_cluster_usage[lcId]++; if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index 4c7ec55367b8f..4bcfbe4dd58a0 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -265,7 +265,7 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es for (unsigned i = 0; i < trackstersEM.size(); ++i) { auto mergedIdx = indexInMergedCollEM[i]; usedTrackstersMerged[mergedIdx] = true; - auto &t = trackstersEM[i]; //trackster + const auto &t = trackstersEM[i]; //trackster TICLCandidate tmpCandidate; tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); tmpCandidate.setCharge(0); @@ -284,7 +284,7 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es for (unsigned i = 0; i < trackstersHAD.size(); ++i) { auto mergedIdx = indexInMergedCollHAD[i]; usedTrackstersMerged[mergedIdx] = true; - auto &t = trackstersHAD[i]; //trackster + const auto &t = trackstersHAD[i]; //trackster TICLCandidate tmpCandidate; tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); tmpCandidate.setCharge(0); @@ -303,7 +303,7 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es for (unsigned i = 0; i < trackstersTRKEM.size(); ++i) { auto mergedIdx = indexInMergedCollTRKEM[i]; if (!usedTrackstersMerged[mergedIdx]) { - auto &t = trackstersTRKEM[i]; //trackster + const auto &t = trackstersTRKEM[i]; //trackster auto trackIdx = t.seedIndex(); auto const &track = tracks[trackIdx]; if (!usedSeeds[trackIdx] and t.raw_energy() > 0) { @@ -450,7 +450,7 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es for (unsigned i = 0; i < trackstersTRK.size(); ++i) { auto mergedIdx = indexInMergedCollTRK[i]; - auto &t = trackstersTRK[i]; //trackster + const auto &t = trackstersTRK[i]; //trackster if (!usedTrackstersMerged[mergedIdx] and t.raw_energy() > 0) { auto trackIdx = t.seedIndex(); diff --git a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc index bfe0dc34fafc5..5adc6768be73d 100644 --- a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc +++ b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc @@ -1616,7 +1616,7 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist cPOnLayer[h.clusterId][lcLayerId].layerClusterIdToEnergyAndScore[mclId].second = FLT_MAX; //cpsInMultiCluster[multicluster][CPids] //Connects a multi cluster with all related caloparticles. - cpsInMultiCluster[mclId].emplace_back(std::make_pair(h.clusterId, FLT_MAX)); + cpsInMultiCluster[mclId].emplace_back(h.clusterId, FLT_MAX); //From all CaloParticles related to a layer cluster, he saves id and energy of the calo particle //that after simhit-rechit matching in layer has the maximum energy. if (shared_fraction > maxCPEnergyInLC) { @@ -1633,7 +1633,7 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist //Loop through all rechits to count how many of them are noise and how many are matched. //In case of matched rechit-simhit, he counts and saves the number of rechits related to the maximum energy CaloParticle. - for (auto& c : hitsToCaloParticleId) { + for (auto c : hitsToCaloParticleId) { if (c < 0) { numberOfNoiseHitsInMCL++; } else { From f2dda06578678baaf6373d64552f0dbab9f74399 Mon Sep 17 00:00:00 2001 From: Felice Date: Fri, 23 Oct 2020 14:32:36 +0200 Subject: [PATCH 46/66] remove unused import --- RecoHGCal/TICL/python/iterativeTICL_cff.py | 1 - 1 file changed, 1 deletion(-) diff --git a/RecoHGCal/TICL/python/iterativeTICL_cff.py b/RecoHGCal/TICL/python/iterativeTICL_cff.py index c6085a16e0370..9af6dc68a9af6 100644 --- a/RecoHGCal/TICL/python/iterativeTICL_cff.py +++ b/RecoHGCal/TICL/python/iterativeTICL_cff.py @@ -6,7 +6,6 @@ from RecoHGCal.TICL.EMStep_cff import * from RecoHGCal.TICL.HADStep_cff import * from RecoHGCal.TICL.ticlLayerTileProducer_cfi import ticlLayerTileProducer -from RecoHGCal.TICL.ticlCandidateFromTrackstersProducer_cfi import ticlCandidateFromTrackstersProducer as _ticlCandidateFromTrackstersProducer from RecoHGCal.TICL.pfTICLProducer_cfi import pfTICLProducer as _pfTICLProducer from RecoHGCal.TICL.trackstersMergeProducer_cfi import trackstersMergeProducer as _trackstersMergeProducer from RecoHGCal.TICL.multiClustersFromTrackstersProducer_cfi import multiClustersFromTrackstersProducer as _multiClustersFromTrackstersProducer From 882ff98a8fe81cb66f4ef0c0be2fb090564c439c Mon Sep 17 00:00:00 2001 From: Felice Date: Tue, 27 Oct 2020 17:19:30 +0100 Subject: [PATCH 47/66] Fix Validation for empty tracksters --- .../src/HGVHistoProducerAlgo.cc | 248 +++++++++--------- 1 file changed, 128 insertions(+), 120 deletions(-) diff --git a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc index 5adc6768be73d..cfbcbce357200 100644 --- a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc +++ b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc @@ -1697,106 +1697,106 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist //Loop through multiclusters for (unsigned int mclId = 0; mclId < nMultiClusters; ++mclId) { const auto& hits_and_fractions = multiClusters[mclId].hitsAndFractions(); + if (!hits_and_fractions.empty()) { + // find the unique caloparticles id contributing to the multi clusters + //cpsInMultiCluster[multicluster][CPids] + std::sort(cpsInMultiCluster[mclId].begin(), cpsInMultiCluster[mclId].end()); + auto last = std::unique(cpsInMultiCluster[mclId].begin(), cpsInMultiCluster[mclId].end()); + cpsInMultiCluster[mclId].erase(last, cpsInMultiCluster[mclId].end()); + + if (multiClusters[mclId].energy() == 0. && !cpsInMultiCluster[mclId].empty()) { + //Loop through all CaloParticles contributing to multicluster mclId. + for (auto& cpPair : cpsInMultiCluster[mclId]) { + //In case of a multi cluster with zero energy but related CaloParticles the score is set to 1. + cpPair.second = 1.; + // LogDebug("HGCalValidator") << "multiCluster Id: \t" << mclId + // << "\t CP id: \t" << cpPair.first + // << "\t score \t" << cpPair.second + // << "\n"; + LogDebug("HGCalValidator") << "multiCluster Id: \t" << mclId << "\t CP id: \t" << cpPair.first + << "\t score \t" << cpPair.second << std::endl; + histograms.h_score_multicl2caloparticle[count]->Fill(cpPair.second); + } + continue; + } + + // Compute the correct normalization + float invMultiClusterEnergyWeight = 0.f; + for (auto const& haf : multiClusters[mclId].hitsAndFractions()) { + invMultiClusterEnergyWeight += + (haf.second * hitMap.at(haf.first)->energy()) * (haf.second * hitMap.at(haf.first)->energy()); + } + invMultiClusterEnergyWeight = 1.f / invMultiClusterEnergyWeight; + + unsigned int numberOfHitsInLC = hits_and_fractions.size(); + for (unsigned int i = 0; i < numberOfHitsInLC; ++i) { + DetId rh_detid = hits_and_fractions[i].first; + float rhFraction = hits_and_fractions[i].second; + bool hitWithNoCP = false; + + auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid); + if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) + hitWithNoCP = true; + auto itcheck = hitMap.find(rh_detid); + const HGCRecHit* hit = itcheck->second; + float hitEnergyWeight = hit->energy() * hit->energy(); + + for (auto& cpPair : cpsInMultiCluster[mclId]) { + float cpFraction = 0.f; + if (!hitWithNoCP) { + auto findHitIt = std::find(detIdToCaloParticleId_Map[rh_detid].begin(), + detIdToCaloParticleId_Map[rh_detid].end(), + HGVHistoProducerAlgo::detIdInfoInCluster{cpPair.first, 0.f}); + if (findHitIt != detIdToCaloParticleId_Map[rh_detid].end()) { + cpFraction = findHitIt->fraction; + } + } + if (cpPair.second == FLT_MAX) { + cpPair.second = 0.f; + } + cpPair.second += + (rhFraction - cpFraction) * (rhFraction - cpFraction) * hitEnergyWeight * invMultiClusterEnergyWeight; + } + } //end of loop through rechits of layer cluster - // find the unique caloparticles id contributing to the multi clusters - //cpsInMultiCluster[multicluster][CPids] - std::sort(cpsInMultiCluster[mclId].begin(), cpsInMultiCluster[mclId].end()); - auto last = std::unique(cpsInMultiCluster[mclId].begin(), cpsInMultiCluster[mclId].end()); - cpsInMultiCluster[mclId].erase(last, cpsInMultiCluster[mclId].end()); + //In case of a multi cluster with some energy but none related CaloParticles print some info. + if (cpsInMultiCluster[mclId].empty()) + LogDebug("HGCalValidator") << "multiCluster Id: \t" << mclId << "\tCP id:\t-1 " + << "\t score \t-1" + << "\n"; - if (multiClusters[mclId].energy() == 0. && !cpsInMultiCluster[mclId].empty()) { - //Loop through all CaloParticles contributing to multicluster mclId. + auto score = std::min_element(std::begin(cpsInMultiCluster[mclId]), + std::end(cpsInMultiCluster[mclId]), + [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; }); for (auto& cpPair : cpsInMultiCluster[mclId]) { - //In case of a multi cluster with zero energy but related CaloParticles the score is set to 1. - cpPair.second = 1.; // LogDebug("HGCalValidator") << "multiCluster Id: \t" << mclId // << "\t CP id: \t" << cpPair.first // << "\t score \t" << cpPair.second // << "\n"; LogDebug("HGCalValidator") << "multiCluster Id: \t" << mclId << "\t CP id: \t" << cpPair.first << "\t score \t" << cpPair.second << std::endl; - histograms.h_score_multicl2caloparticle[count]->Fill(cpPair.second); - } - continue; - } - - // Compute the correct normalization - float invMultiClusterEnergyWeight = 0.f; - for (auto const& haf : multiClusters[mclId].hitsAndFractions()) { - invMultiClusterEnergyWeight += - (haf.second * hitMap.at(haf.first)->energy()) * (haf.second * hitMap.at(haf.first)->energy()); - } - invMultiClusterEnergyWeight = 1.f / invMultiClusterEnergyWeight; - - unsigned int numberOfHitsInLC = hits_and_fractions.size(); - for (unsigned int i = 0; i < numberOfHitsInLC; ++i) { - DetId rh_detid = hits_and_fractions[i].first; - float rhFraction = hits_and_fractions[i].second; - bool hitWithNoCP = false; - - auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid); - if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) - hitWithNoCP = true; - auto itcheck = hitMap.find(rh_detid); - const HGCRecHit* hit = itcheck->second; - float hitEnergyWeight = hit->energy() * hit->energy(); - - for (auto& cpPair : cpsInMultiCluster[mclId]) { - float cpFraction = 0.f; - if (!hitWithNoCP) { - auto findHitIt = std::find(detIdToCaloParticleId_Map[rh_detid].begin(), - detIdToCaloParticleId_Map[rh_detid].end(), - HGVHistoProducerAlgo::detIdInfoInCluster{cpPair.first, 0.f}); - if (findHitIt != detIdToCaloParticleId_Map[rh_detid].end()) { - cpFraction = findHitIt->fraction; - } + if (cpPair.first == score->first) { + histograms.h_score_multicl2caloparticle[count]->Fill(score->second); } - if (cpPair.second == FLT_MAX) { - cpPair.second = 0.f; + float sharedeneCPallLayers = 0.; + //Loop through all layers + for (unsigned int j = 0; j < layers * 2; ++j) { + auto const& cp_linked = cPOnLayer[cpPair.first][j].layerClusterIdToEnergyAndScore[mclId]; + sharedeneCPallLayers += cp_linked.first; + } //end of loop through layers + LogDebug("HGCalValidator") << "sharedeneCPallLayers " << sharedeneCPallLayers << std::endl; + if (cpPair.first == score->first) { + histograms.h_sharedenergy_multicl2caloparticle[count]->Fill(sharedeneCPallLayers / + multiClusters[mclId].energy()); + histograms.h_energy_vs_score_multicl2caloparticle[count]->Fill( + score->second, sharedeneCPallLayers / multiClusters[mclId].energy()); } - cpPair.second += - (rhFraction - cpFraction) * (rhFraction - cpFraction) * hitEnergyWeight * invMultiClusterEnergyWeight; - } - } //end of loop through rechits of layer cluster - - //In case of a multi cluster with some energy but none related CaloParticles print some info. - if (cpsInMultiCluster[mclId].empty()) - LogDebug("HGCalValidator") << "multiCluster Id: \t" << mclId << "\tCP id:\t-1 " - << "\t score \t-1" - << "\n"; - - auto score = std::min_element(std::begin(cpsInMultiCluster[mclId]), - std::end(cpsInMultiCluster[mclId]), - [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; }); - for (auto& cpPair : cpsInMultiCluster[mclId]) { - // LogDebug("HGCalValidator") << "multiCluster Id: \t" << mclId - // << "\t CP id: \t" << cpPair.first - // << "\t score \t" << cpPair.second - // << "\n"; - LogDebug("HGCalValidator") << "multiCluster Id: \t" << mclId << "\t CP id: \t" << cpPair.first << "\t score \t" - << cpPair.second << std::endl; - if (cpPair.first == score->first) { - histograms.h_score_multicl2caloparticle[count]->Fill(score->second); - } - float sharedeneCPallLayers = 0.; - //Loop through all layers - for (unsigned int j = 0; j < layers * 2; ++j) { - auto const& cp_linked = cPOnLayer[cpPair.first][j].layerClusterIdToEnergyAndScore[mclId]; - sharedeneCPallLayers += cp_linked.first; - } //end of loop through layers - LogDebug("HGCalValidator") << "sharedeneCPallLayers " << sharedeneCPallLayers << std::endl; - if (cpPair.first == score->first) { - histograms.h_sharedenergy_multicl2caloparticle[count]->Fill(sharedeneCPallLayers / - multiClusters[mclId].energy()); - histograms.h_energy_vs_score_multicl2caloparticle[count]->Fill( - score->second, sharedeneCPallLayers / multiClusters[mclId].energy()); } + auto assocFakeMerge = std::count_if(std::begin(cpsInMultiCluster[mclId]), + std::end(cpsInMultiCluster[mclId]), + [](const auto& obj) { return obj.second < ScoreCutMCLtoCPFakeMerge_; }); + tracksters_fakemerge[mclId] = assocFakeMerge; } - auto assocFakeMerge = std::count_if(std::begin(cpsInMultiCluster[mclId]), - std::end(cpsInMultiCluster[mclId]), - [](const auto& obj) { return obj.second < ScoreCutMCLtoCPFakeMerge_; }); - tracksters_fakemerge[mclId] = assocFakeMerge; - } //end of loop through multiclusters std::unordered_map> score3d; @@ -1978,37 +1978,40 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist // reco-level, namely fake-rate an merge-rate. In this loop we should *not* // restrict only to the selected caloParaticles. for (unsigned int mclId = 0; mclId < nMultiClusters; ++mclId) { - auto assocFakeMerge = tracksters_fakemerge[mclId]; - auto assocDuplicate = tracksters_duplicate[mclId]; - if (assocDuplicate) { - histograms.h_numDup_multicl_eta[count]->Fill(multiClusters[mclId].eta()); - histograms.h_numDup_multicl_phi[count]->Fill(multiClusters[mclId].phi()); - } - if (assocFakeMerge > 0) { - histograms.h_num_multicl_eta[count]->Fill(multiClusters[mclId].eta()); - histograms.h_num_multicl_phi[count]->Fill(multiClusters[mclId].phi()); - auto best = std::min_element(std::begin(cpsInMultiCluster[mclId]), - std::end(cpsInMultiCluster[mclId]), - [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; }); - - //This is the shared energy taking the best caloparticle in each layer - float sharedeneCPallLayers = 0.; - //Loop through all layers - for (unsigned int j = 0; j < layers * 2; ++j) { - auto const& best_cp_linked = cPOnLayer[best->first][j].layerClusterIdToEnergyAndScore[mclId]; - sharedeneCPallLayers += best_cp_linked.first; - } //end of loop through layers - histograms.h_sharedenergy_multicl2caloparticle_vs_eta[count]->Fill( - multiClusters[mclId].eta(), sharedeneCPallLayers / multiClusters[mclId].energy()); - histograms.h_sharedenergy_multicl2caloparticle_vs_phi[count]->Fill( - multiClusters[mclId].phi(), sharedeneCPallLayers / multiClusters[mclId].energy()); - } - if (assocFakeMerge >= 2) { - histograms.h_numMerge_multicl_eta[count]->Fill(multiClusters[mclId].eta()); - histograms.h_numMerge_multicl_phi[count]->Fill(multiClusters[mclId].phi()); + const auto& hits_and_fractions = multiClusters[mclId].hitsAndFractions(); + if (!hits_and_fractions.empty()) { + auto assocFakeMerge = tracksters_fakemerge[mclId]; + auto assocDuplicate = tracksters_duplicate[mclId]; + if (assocDuplicate) { + histograms.h_numDup_multicl_eta[count]->Fill(multiClusters[mclId].eta()); + histograms.h_numDup_multicl_phi[count]->Fill(multiClusters[mclId].phi()); + } + if (assocFakeMerge > 0) { + histograms.h_num_multicl_eta[count]->Fill(multiClusters[mclId].eta()); + histograms.h_num_multicl_phi[count]->Fill(multiClusters[mclId].phi()); + auto best = std::min_element(std::begin(cpsInMultiCluster[mclId]), + std::end(cpsInMultiCluster[mclId]), + [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; }); + + //This is the shared energy taking the best caloparticle in each layer + float sharedeneCPallLayers = 0.; + //Loop through all layers + for (unsigned int j = 0; j < layers * 2; ++j) { + auto const& best_cp_linked = cPOnLayer[best->first][j].layerClusterIdToEnergyAndScore[mclId]; + sharedeneCPallLayers += best_cp_linked.first; + } //end of loop through layers + histograms.h_sharedenergy_multicl2caloparticle_vs_eta[count]->Fill( + multiClusters[mclId].eta(), sharedeneCPallLayers / multiClusters[mclId].energy()); + histograms.h_sharedenergy_multicl2caloparticle_vs_phi[count]->Fill( + multiClusters[mclId].phi(), sharedeneCPallLayers / multiClusters[mclId].energy()); + } + if (assocFakeMerge >= 2) { + histograms.h_numMerge_multicl_eta[count]->Fill(multiClusters[mclId].eta()); + histograms.h_numMerge_multicl_phi[count]->Fill(multiClusters[mclId].phi()); + } + histograms.h_denom_multicl_eta[count]->Fill(multiClusters[mclId].eta()); + histograms.h_denom_multicl_phi[count]->Fill(multiClusters[mclId].phi()); } - histograms.h_denom_multicl_eta[count]->Fill(multiClusters[mclId].eta()); - histograms.h_denom_multicl_phi[count]->Fill(multiClusters[mclId].phi()); } } @@ -2051,6 +2054,10 @@ void HGVHistoProducerAlgo::fill_multi_cluster_histos(const Histograms& histogram for (unsigned int mclId = 0; mclId < nMultiClusters; ++mclId) { const auto layerClusters = multiClusters[mclId].clusters(); auto nLayerClusters = layerClusters.size(); + + if (nLayerClusters == 0) + continue; + if (multiClusters[mclId].z() < 0.) { tnmclmz++; } @@ -2179,10 +2186,11 @@ void HGVHistoProducerAlgo::fill_multi_cluster_histos(const Histograms& histogram histograms.h_multicluster_firstlayer[count]->Fill((float)*multicluster_layers.begin()); histograms.h_multicluster_lastlayer[count]->Fill((float)*multicluster_layers.rbegin()); histograms.h_multicluster_layersnum[count]->Fill((float)multicluster_layers.size()); - } - histograms.h_multicluster_pt[count]->Fill(multiClusters[mclId].pt()); - histograms.h_multicluster_energy[count]->Fill(multiClusters[mclId].energy()); + histograms.h_multicluster_pt[count]->Fill(multiClusters[mclId].pt()); + + histograms.h_multicluster_energy[count]->Fill(multiClusters[mclId].energy()); + } } //end of loop through multiclusters From 202b86505f08f7e941e3225f20d87378133691e0 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 28 Oct 2020 13:21:17 +0100 Subject: [PATCH 48/66] modify comment to reflect code --- RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc b/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc index c438a92be6ef9..27cf92d6c1265 100644 --- a/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc @@ -63,7 +63,7 @@ void MultiClustersFromTrackstersProducer::produce(edm::Event& evt, const edm::Ev } std::for_each(std::begin(tracksters), std::end(tracksters), [&](auto const& trackster) { - // Do not create a multicluster if the trackster has no layer clusters. + // Create an empty multicluster if the trackster has no layer clusters. // This could happen when a seed leads to no trackster and a dummy one is produced. std::array baricenter{{0., 0., 0.}}; From 5289d77206d07df5c33c40a60038d4ee04d70272 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 28 Oct 2020 13:42:24 +0100 Subject: [PATCH 49/66] replace parameter root_doublet_max_distance_from_seed with its square --- RecoHGCal/TICL/plugins/HGCGraph.cc | 12 +++++++----- RecoHGCal/TICL/plugins/HGCGraph.h | 2 +- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc | 4 ++-- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h | 2 +- RecoHGCal/TICL/plugins/TrackstersProducer.cc | 2 +- RecoHGCal/TICL/python/TrkEMStep_cff.py | 2 +- 6 files changed, 13 insertions(+), 11 deletions(-) diff --git a/RecoHGCal/TICL/plugins/HGCGraph.cc b/RecoHGCal/TICL/plugins/HGCGraph.cc index a69996468399a..e94b87426bc95 100644 --- a/RecoHGCal/TICL/plugins/HGCGraph.cc +++ b/RecoHGCal/TICL/plugins/HGCGraph.cc @@ -19,7 +19,7 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, int deltaIPhi, float minCosTheta, float minCosPointing, - float root_doublet_max_distance_from_seed, + float root_doublet_max_distance_from_seed_squared, float etaLimitIncreaseWindow, int skip_layers, int maxNumberOfLayers, @@ -28,7 +28,7 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, isOuterClusterOfDoublets_.resize(layerClusters.size()); allDoublets_.clear(); theRootDoublets_.clear(); - bool checkDistanceRootDoubletVsSeed = root_doublet_max_distance_from_seed < 9999; + bool checkDistanceRootDoubletVsSeed = root_doublet_max_distance_from_seed_squared < 9999; float origin_eta; float origin_phi; for (const auto &r : regions) { @@ -179,9 +179,11 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, minCosPointing, verbosity_ > Advanced); if (isRootDoublet and checkDistanceRootDoubletVsSeed) { - auto lcEta = layerClusters[innerClusterId].eta(); - auto lcPhi = layerClusters[innerClusterId].phi(); - if (std::hypot(lcEta - origin_eta, lcPhi - origin_phi) > root_doublet_max_distance_from_seed) + auto dEtaSquared = (layerClusters[innerClusterId].eta() - origin_eta); + dEtaSquared *= dEtaSquared; + auto dPhiSquared = (layerClusters[innerClusterId].phi() - origin_phi); + dPhiSquared *= dPhiSquared; + if (dEtaSquared + dPhiSquared > root_doublet_max_distance_from_seed_squared) isRootDoublet = false; } if (isRootDoublet) { diff --git a/RecoHGCal/TICL/plugins/HGCGraph.h b/RecoHGCal/TICL/plugins/HGCGraph.h index 836b024179c04..f8cc626016f40 100644 --- a/RecoHGCal/TICL/plugins/HGCGraph.h +++ b/RecoHGCal/TICL/plugins/HGCGraph.h @@ -25,7 +25,7 @@ class HGCGraphT { int deltaIPhi, float minCosThetai, float maxCosPointing, - float root_doublet_max_distance_from_seed, + float root_doublet_max_distance_from_seed_squared, float etaLimitIncreaseWindow, int skip_layers, int maxNumberOfLayers, diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index 10ba4b0cf1e87..777abf2c7ffe8 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -26,7 +26,7 @@ PatternRecognitionbyCA::PatternRecognitionbyCA(const edm::ParameterSet &c max_out_in_hops_(conf.getParameter("max_out_in_hops")), min_cos_theta_(conf.getParameter("min_cos_theta")), min_cos_pointing_(conf.getParameter("min_cos_pointing")), - root_doublet_max_distance_from_seed_(conf.getParameter("root_doublet_max_distance_from_seed")), + root_doublet_max_distance_from_seed_squared_(conf.getParameter("root_doublet_max_distance_from_seed_squared")), etaLimitIncreaseWindow_(conf.getParameter("etaLimitIncreaseWindow")), skip_layers_(conf.getParameter("skip_layers")), max_missing_layers_in_trackster_(conf.getParameter("max_missing_layers_in_trackster")), @@ -97,7 +97,7 @@ void PatternRecognitionbyCA::makeTracksters( 1, min_cos_theta_, min_cos_pointing_, - root_doublet_max_distance_from_seed_, + root_doublet_max_distance_from_seed_squared_, etaLimitIncreaseWindow_, skip_layers_, rhtools_.lastLayer(type), diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h index c96bf265e3f99..c337dd344ec78 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h @@ -37,7 +37,7 @@ namespace ticl { const unsigned int max_out_in_hops_; const float min_cos_theta_; const float min_cos_pointing_; - const float root_doublet_max_distance_from_seed_; + const float root_doublet_max_distance_from_seed_squared_; const float etaLimitIncreaseWindow_; const int skip_layers_; const int max_missing_layers_in_trackster_; diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc index 7202385174a30..38650eeab1211 100644 --- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc @@ -118,7 +118,7 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri desc.add("shower_start_max_layer", 9999); // make default such that no filtering is applied desc.add("algo_verbosity", 0); desc.add("min_cos_theta", 0.915); - desc.add("root_doublet_max_distance_from_seed", 9999); + desc.add("root_doublet_max_distance_from_seed_squared", 9999); desc.add("min_cos_pointing", -1.); desc.add("skip_layers", 0); desc.add("max_missing_layers_in_trackster", 9999); diff --git a/RecoHGCal/TICL/python/TrkEMStep_cff.py b/RecoHGCal/TICL/python/TrkEMStep_cff.py index 9f9c4c5f94fa1..458f489f696ec 100644 --- a/RecoHGCal/TICL/python/TrkEMStep_cff.py +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -31,7 +31,7 @@ min_layers_per_trackster = 10, min_cos_theta = 0.97, # ~14 degrees min_cos_pointing = 0.94, # ~20 degrees - root_doublet_max_distance_from_seed = 0.05, + root_doublet_max_distance_from_seed_squared = 2.5e-3, # dR=0.05 max_delta_time = 3., itername = "TrkEM", algo_verbosity = 0, From ec4a82b7ff71f429f53a6d2a435b304e86bc9c8a Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 28 Oct 2020 13:43:12 +0100 Subject: [PATCH 50/66] Optimization --- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index 777abf2c7ffe8..eb1554c446970 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -225,7 +225,7 @@ void PatternRecognitionbyCA::makeTracksters( cumulative_prob += t.id_probabilities(index); } return (cumulative_prob <= pid_threshold_) && - (t.raw_em_energy() / t.raw_energy() < energy_em_over_total_threshold_); + (t.raw_em_energy() < energy_em_over_total_threshold_*t.raw_energy()); }; std::vector selectedTrackstersIds; From d771885862064e3fda716e24ad692695d9b90065 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 28 Oct 2020 13:46:14 +0100 Subject: [PATCH 51/66] remove unused parameter --- RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc | 3 --- 1 file changed, 3 deletions(-) diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index 4bcfbe4dd58a0..9248b7cee365d 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -55,7 +55,6 @@ class TrackstersMergeProducer : public edm::stream::EDProducer("pt_sigma_high")), pt_sigma_low_(ps.getParameter("pt_sigma_low")), halo_max_distance2_(ps.getParameter("halo_max_distance2")), - mergeTRK_max_dR_(ps.getParameter("mergeTRK_max_dR")), cosangle_align_(ps.getParameter("cosangle_align")), e_over_h_threshold_(ps.getParameter("e_over_h_threshold")), pt_neutral_threshold_(ps.getParameter("pt_neutral_threshold")), @@ -714,7 +712,6 @@ void TrackstersMergeProducer::fillDescriptions(edm::ConfigurationDescriptions &d desc.add("pt_sigma_high", 2.); desc.add("pt_sigma_low", 2.); desc.add("halo_max_distance2", 4.); - desc.add("mergeTRK_max_dR", 0.03); desc.add("cosangle_align", 0.9945); desc.add("e_over_h_threshold", 1.); desc.add("pt_neutral_threshold", 2.); From 69128ffcfe107405b6c168c55dc3e7a10ffe7c69 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 28 Oct 2020 13:49:00 +0100 Subject: [PATCH 52/66] remove commented out code --- RecoHGCal/TICL/python/TrkStep_cff.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/RecoHGCal/TICL/python/TrkStep_cff.py b/RecoHGCal/TICL/python/TrkStep_cff.py index 697441571b52e..739a7cafd5d2b 100644 --- a/RecoHGCal/TICL/python/TrkStep_cff.py +++ b/RecoHGCal/TICL/python/TrkStep_cff.py @@ -26,8 +26,6 @@ pid_threshold = 0.0, skip_layers = 3, min_layers_per_trackster = 10, -# min_cos_theta = 0.978, # same as em -# min_cos_pointing = 0.9, # same as em min_cos_theta = 0.866, # ~30 degrees min_cos_pointing = 0.798, # ~ 37 degrees max_delta_time = -1., From 0d4b6f002dd652ff7da42d5b00cfb86965ec54c0 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 28 Oct 2020 14:28:48 +0100 Subject: [PATCH 53/66] optimization --- RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index 9248b7cee365d..4f9338ae43cd5 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -1,5 +1,4 @@ #include // unique_ptr -#include // cosh #include "FWCore/Framework/interface/stream/EDProducer.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" @@ -351,7 +350,8 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); } - float p = tracksterTotalRawPt * cosh(t.barycenter().eta()); + const auto& barycenter = t.barycenter(); + float p = tracksterTotalRawPt * std::sqrt(1+(barycenter.z()*barycenter.z())/(barycenter.x()*barycenter.x()+(barycenter.y()*barycenter.y()))); float energy = std::sqrt(p * p + mpion2); tmpCandidate.setRawEnergy(energy); math::XYZTLorentzVector p4(p * track.momentum().unit().x(), @@ -370,7 +370,8 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); } tmpCandidate.setPdgId(11 * track.charge()); - float p = tracksterTotalRawPt * cosh(t.barycenter().eta()); + const auto& barycenter = t.barycenter(); + float p = tracksterTotalRawPt * std::sqrt(1+(barycenter.z()*barycenter.z())/(barycenter.x()*barycenter.x()+(barycenter.y()*barycenter.y()))); tmpCandidate.setRawEnergy(p); math::XYZTLorentzVector p4( p * track.momentum().unit().x(), p * track.momentum().unit().y(), p * track.momentum().unit().z(), p); @@ -398,7 +399,8 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); } - float p = tracksterTotalRawPt * cosh(t.barycenter().eta()); + const auto& barycenter = t.barycenter(); + float p = tracksterTotalRawPt * std::sqrt(1+(barycenter.z()*barycenter.z())/(barycenter.x()*barycenter.x()+(barycenter.y()*barycenter.y()))); float energy = std::sqrt(p * p + mpion2); tmpCandidate.setRawEnergy(energy); math::XYZTLorentzVector p4(p * track.momentum().unit().x(), @@ -456,7 +458,8 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es if (!usedSeeds[trackIdx]) { usedSeeds[trackIdx] = true; usedTrackstersMerged[mergedIdx] = true; - float p = t.raw_pt() * cosh(t.barycenter().eta()); + const auto& barycenter = t.barycenter(); + float p = t.raw_pt() * std::sqrt(1+(barycenter.z()*barycenter.z())/(barycenter.x()*barycenter.x()+(barycenter.y()*barycenter.y()))); float energy = std::sqrt(p * p + mpion2); TICLCandidate tmpCandidate; tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); From a61f591c426bab2c5a8708139ff3a6d7c5eb9be6 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 28 Oct 2020 14:30:37 +0100 Subject: [PATCH 54/66] code format --- .../TICL/plugins/TrackstersMergeProducer.cc | 23 ++++++++++++------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index 4f9338ae43cd5..b44e5be141a7a 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -350,8 +350,10 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); } - const auto& barycenter = t.barycenter(); - float p = tracksterTotalRawPt * std::sqrt(1+(barycenter.z()*barycenter.z())/(barycenter.x()*barycenter.x()+(barycenter.y()*barycenter.y()))); + const auto &barycenter = t.barycenter(); + float p = tracksterTotalRawPt * + std::sqrt(1 + (barycenter.z() * barycenter.z()) / + (barycenter.x() * barycenter.x() + (barycenter.y() * barycenter.y()))); float energy = std::sqrt(p * p + mpion2); tmpCandidate.setRawEnergy(energy); math::XYZTLorentzVector p4(p * track.momentum().unit().x(), @@ -370,8 +372,10 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); } tmpCandidate.setPdgId(11 * track.charge()); - const auto& barycenter = t.barycenter(); - float p = tracksterTotalRawPt * std::sqrt(1+(barycenter.z()*barycenter.z())/(barycenter.x()*barycenter.x()+(barycenter.y()*barycenter.y()))); + const auto &barycenter = t.barycenter(); + float p = tracksterTotalRawPt * + std::sqrt(1 + (barycenter.z() * barycenter.z()) / + (barycenter.x() * barycenter.x() + (barycenter.y() * barycenter.y()))); tmpCandidate.setRawEnergy(p); math::XYZTLorentzVector p4( p * track.momentum().unit().x(), p * track.momentum().unit().y(), p * track.momentum().unit().z(), p); @@ -399,8 +403,10 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); } - const auto& barycenter = t.barycenter(); - float p = tracksterTotalRawPt * std::sqrt(1+(barycenter.z()*barycenter.z())/(barycenter.x()*barycenter.x()+(barycenter.y()*barycenter.y()))); + const auto &barycenter = t.barycenter(); + float p = tracksterTotalRawPt * + std::sqrt(1 + (barycenter.z() * barycenter.z()) / + (barycenter.x() * barycenter.x() + (barycenter.y() * barycenter.y()))); float energy = std::sqrt(p * p + mpion2); tmpCandidate.setRawEnergy(energy); math::XYZTLorentzVector p4(p * track.momentum().unit().x(), @@ -458,8 +464,9 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es if (!usedSeeds[trackIdx]) { usedSeeds[trackIdx] = true; usedTrackstersMerged[mergedIdx] = true; - const auto& barycenter = t.barycenter(); - float p = t.raw_pt() * std::sqrt(1+(barycenter.z()*barycenter.z())/(barycenter.x()*barycenter.x()+(barycenter.y()*barycenter.y()))); + const auto &barycenter = t.barycenter(); + float p = t.raw_pt() * std::sqrt(1 + (barycenter.z() * barycenter.z()) / + (barycenter.x() * barycenter.x() + (barycenter.y() * barycenter.y()))); float energy = std::sqrt(p * p + mpion2); TICLCandidate tmpCandidate; tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); From b4a8cecc228f9dc235bc9f60be5718454c515f65 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 28 Oct 2020 14:30:44 +0100 Subject: [PATCH 55/66] code format --- RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index eb1554c446970..8340aa261fea8 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -26,7 +26,8 @@ PatternRecognitionbyCA::PatternRecognitionbyCA(const edm::ParameterSet &c max_out_in_hops_(conf.getParameter("max_out_in_hops")), min_cos_theta_(conf.getParameter("min_cos_theta")), min_cos_pointing_(conf.getParameter("min_cos_pointing")), - root_doublet_max_distance_from_seed_squared_(conf.getParameter("root_doublet_max_distance_from_seed_squared")), + root_doublet_max_distance_from_seed_squared_( + conf.getParameter("root_doublet_max_distance_from_seed_squared")), etaLimitIncreaseWindow_(conf.getParameter("etaLimitIncreaseWindow")), skip_layers_(conf.getParameter("skip_layers")), max_missing_layers_in_trackster_(conf.getParameter("max_missing_layers_in_trackster")), @@ -225,7 +226,7 @@ void PatternRecognitionbyCA::makeTracksters( cumulative_prob += t.id_probabilities(index); } return (cumulative_prob <= pid_threshold_) && - (t.raw_em_energy() < energy_em_over_total_threshold_*t.raw_energy()); + (t.raw_em_energy() < energy_em_over_total_threshold_ * t.raw_energy()); }; std::vector selectedTrackstersIds; From 214708151c71e45c269a02c12da3184e1376e5e7 Mon Sep 17 00:00:00 2001 From: Felice Date: Wed, 28 Oct 2020 14:30:49 +0100 Subject: [PATCH 56/66] code format --- Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc | 4 ---- 1 file changed, 4 deletions(-) diff --git a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc index cfbcbce357200..9e45bb02da062 100644 --- a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc +++ b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc @@ -1709,10 +1709,6 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist for (auto& cpPair : cpsInMultiCluster[mclId]) { //In case of a multi cluster with zero energy but related CaloParticles the score is set to 1. cpPair.second = 1.; - // LogDebug("HGCalValidator") << "multiCluster Id: \t" << mclId - // << "\t CP id: \t" << cpPair.first - // << "\t score \t" << cpPair.second - // << "\n"; LogDebug("HGCalValidator") << "multiCluster Id: \t" << mclId << "\t CP id: \t" << cpPair.first << "\t score \t" << cpPair.second << std::endl; histograms.h_score_multicl2caloparticle[count]->Fill(cpPair.second); From 1c730e156f85688fedea57d04a0b718a4d7e0bf4 Mon Sep 17 00:00:00 2001 From: Felice Date: Fri, 30 Oct 2020 15:19:26 +0100 Subject: [PATCH 57/66] optimization --- RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index b44e5be141a7a..55362a44c65e3 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -425,7 +425,10 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); } tmpCandidate.setPdgId(11 * track.charge()); - float p = tracksterTotalRawPt * cosh(t.barycenter().eta()); + const auto &barycenter = t.barycenter(); + float p = tracksterTotalRawPt * + std::sqrt(1 + (barycenter.z() * barycenter.z()) / + (barycenter.x() * barycenter.x() + (barycenter.y() * barycenter.y()))); tmpCandidate.setRawEnergy(p); math::XYZTLorentzVector p4( p * track.momentum().unit().x(), p * track.momentum().unit().y(), p * track.momentum().unit().z(), p); From 520d3ad986d1de303a87938a0a218df01cf38510 Mon Sep 17 00:00:00 2001 From: Felice Date: Fri, 30 Oct 2020 17:57:01 +0100 Subject: [PATCH 58/66] Simplify energy assignment --- .../TICL/plugins/TrackstersMergeProducer.cc | 80 +++++++++---------- 1 file changed, 40 insertions(+), 40 deletions(-) diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index 55362a44c65e3..b2a2165566789 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -344,41 +344,41 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es if (foundCompatibleTRK) { TICLCandidate tmpCandidate; tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); + double raw_energy = t.raw_energy(); + tmpCandidate.setCharge(track.charge()); tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); tmpCandidate.setPdgId(211 * track.charge()); for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); + raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); + } - const auto &barycenter = t.barycenter(); - float p = tracksterTotalRawPt * - std::sqrt(1 + (barycenter.z() * barycenter.z()) / - (barycenter.x() * barycenter.x() + (barycenter.y() * barycenter.y()))); - float energy = std::sqrt(p * p + mpion2); - tmpCandidate.setRawEnergy(energy); - math::XYZTLorentzVector p4(p * track.momentum().unit().x(), - p * track.momentum().unit().y(), - p * track.momentum().unit().z(), - energy); + tmpCandidate.setRawEnergy(raw_energy); + math::XYZTLorentzVector p4(raw_energy * t.barycenter().unit().x(), + raw_energy * t.barycenter().unit().y(), + raw_energy * t.barycenter().unit().z(), + raw_energy); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); } else { TICLCandidate tmpCandidate; tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); + double raw_energy = t.raw_energy(); tmpCandidate.setCharge(track.charge()); tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); + raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); } tmpCandidate.setPdgId(11 * track.charge()); - const auto &barycenter = t.barycenter(); - float p = tracksterTotalRawPt * - std::sqrt(1 + (barycenter.z() * barycenter.z()) / - (barycenter.x() * barycenter.x() + (barycenter.y() * barycenter.y()))); - tmpCandidate.setRawEnergy(p); - math::XYZTLorentzVector p4( - p * track.momentum().unit().x(), p * track.momentum().unit().y(), p * track.momentum().unit().z(), p); + + tmpCandidate.setRawEnergy(raw_energy); + math::XYZTLorentzVector p4(raw_energy * t.barycenter().unit().x(), + raw_energy * t.barycenter().unit().y(), + raw_energy * t.barycenter().unit().z(), + raw_energy); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); } @@ -390,48 +390,49 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es for (auto otherTracksterIdx : trackstersTRKEMwithSameSeed) { auto thisPt = tracksterTotalRawPt + trackstersMergedHandle->at(otherTracksterIdx).raw_pt() - t.raw_pt(); closestTrackster = std::abs(thisPt - track.pt()) < minPtDiff ? otherTracksterIdx : closestTrackster; - usedTrackstersMerged[otherTracksterIdx] = true; } tracksterTotalRawPt += trackstersMergedHandle->at(closestTrackster).raw_pt() - t.raw_pt(); + usedTrackstersMerged[closestTrackster] = true; if (foundCompatibleTRK) { TICLCandidate tmpCandidate; tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, closestTrackster)); + double raw_energy = trackstersMergedHandle->at(closestTrackster).raw_energy(); + tmpCandidate.setCharge(track.charge()); tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); tmpCandidate.setPdgId(211 * track.charge()); for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); + raw_energy +=trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); } - const auto &barycenter = t.barycenter(); - float p = tracksterTotalRawPt * - std::sqrt(1 + (barycenter.z() * barycenter.z()) / - (barycenter.x() * barycenter.x() + (barycenter.y() * barycenter.y()))); - float energy = std::sqrt(p * p + mpion2); - tmpCandidate.setRawEnergy(energy); - math::XYZTLorentzVector p4(p * track.momentum().unit().x(), - p * track.momentum().unit().y(), - p * track.momentum().unit().z(), - energy); + tmpCandidate.setRawEnergy(raw_energy); + float momentum = std::sqrt(raw_energy * raw_energy - mpion2); + const auto &barycenter = trackstersMergedHandle->at(closestTrackster).barycenter(); + math::XYZTLorentzVector p4(momentum * barycenter.unit().x(), + momentum * barycenter.unit().y(), + momentum * barycenter.unit().z(), + t.raw_energy()); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); } else { TICLCandidate tmpCandidate; tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, closestTrackster)); + double raw_energy = trackstersMergedHandle->at(closestTrackster).raw_energy(); + tmpCandidate.setCharge(track.charge()); tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); + raw_energy +=trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); } tmpCandidate.setPdgId(11 * track.charge()); - const auto &barycenter = t.barycenter(); - float p = tracksterTotalRawPt * - std::sqrt(1 + (barycenter.z() * barycenter.z()) / - (barycenter.x() * barycenter.x() + (barycenter.y() * barycenter.y()))); - tmpCandidate.setRawEnergy(p); - math::XYZTLorentzVector p4( - p * track.momentum().unit().x(), p * track.momentum().unit().y(), p * track.momentum().unit().z(), p); + tmpCandidate.setRawEnergy(raw_energy); + math::XYZTLorentzVector p4(raw_energy * t.barycenter().unit().x(), + raw_energy * t.barycenter().unit().y(), + raw_energy * t.barycenter().unit().z(), + raw_energy); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); } @@ -468,17 +469,16 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es usedSeeds[trackIdx] = true; usedTrackstersMerged[mergedIdx] = true; const auto &barycenter = t.barycenter(); - float p = t.raw_pt() * std::sqrt(1 + (barycenter.z() * barycenter.z()) / - (barycenter.x() * barycenter.x() + (barycenter.y() * barycenter.y()))); - float energy = std::sqrt(p * p + mpion2); TICLCandidate tmpCandidate; tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); tmpCandidate.setCharge(track.charge()); tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); tmpCandidate.setPdgId(211 * track.charge()); - tmpCandidate.setRawEnergy(energy); + tmpCandidate.setRawEnergy(t.raw_energy()); + float momentum = std::sqrt(t.raw_energy() * t.raw_energy() - mpion2); + math::XYZTLorentzVector p4( - p * track.momentum().unit().x(), p * track.momentum().unit().y(), p * track.momentum().unit().z(), energy); + momentum * barycenter.unit().x(), momentum * barycenter.unit().y(), momentum * barycenter.unit().z(), t.raw_energy()); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); } From 5dcb7272f5ffd70df03e3bd0a269868f1d30d3a8 Mon Sep 17 00:00:00 2001 From: Felice Date: Fri, 30 Oct 2020 17:57:35 +0100 Subject: [PATCH 59/66] Provide EcalEnergy and HcalEnergy to JetMET for calibration --- RecoHGCal/TICL/plugins/PFTICLProducer.cc | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/RecoHGCal/TICL/plugins/PFTICLProducer.cc b/RecoHGCal/TICL/plugins/PFTICLProducer.cc index 74d69a350b7b4..19a7a8538a17a 100644 --- a/RecoHGCal/TICL/plugins/PFTICLProducer.cc +++ b/RecoHGCal/TICL/plugins/PFTICLProducer.cc @@ -52,6 +52,17 @@ void PFTICLProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSet const auto abs_pdg_id = std::abs(ticl_cand.pdgId()); const auto charge = ticl_cand.charge(); const auto& four_mom = ticl_cand.p4(); + // const auto raw_energy = 0.; + double ecal_energy = 0.; + + for(const auto t: ticl_cand.tracksters()) + { + double ecal_energy_fraction = t->raw_em_pt()/t->raw_pt(); + ecal_energy += t->raw_energy() * ecal_energy_fraction; + } + double hcal_energy = ticl_cand.rawEnergy() - ecal_energy; + // fix for floating point rounding could go slightly below 0 + hcal_energy = hcal_energy < 0? 0:hcal_energy; reco::PFCandidate::ParticleType part_type; switch (abs_pdg_id) { @@ -78,6 +89,8 @@ void PFTICLProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSet candidates->emplace_back(charge, four_mom, part_type); auto& candidate = candidates->back(); + candidate.setEcalEnergy(ecal_energy, ecal_energy); + candidate.setHcalEnergy(hcal_energy, hcal_energy); if (candidate.charge()) { // otherwise PFCandidate throws // Construct edm::Ref from edm::Ptr. As of now, assumes type to be reco::Track. To be extended (either via // dynamic type checking or configuration) if additional track types are needed. From c432855656901b16c4aa548f2efc9e9f2fdbded0 Mon Sep 17 00:00:00 2001 From: Felice Date: Fri, 30 Oct 2020 17:58:51 +0100 Subject: [PATCH 60/66] code format --- RecoHGCal/TICL/plugins/PFTICLProducer.cc | 8 ++--- .../TICL/plugins/TrackstersMergeProducer.cc | 35 ++++++++++--------- 2 files changed, 21 insertions(+), 22 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PFTICLProducer.cc b/RecoHGCal/TICL/plugins/PFTICLProducer.cc index 19a7a8538a17a..09d244908330d 100644 --- a/RecoHGCal/TICL/plugins/PFTICLProducer.cc +++ b/RecoHGCal/TICL/plugins/PFTICLProducer.cc @@ -52,17 +52,15 @@ void PFTICLProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSet const auto abs_pdg_id = std::abs(ticl_cand.pdgId()); const auto charge = ticl_cand.charge(); const auto& four_mom = ticl_cand.p4(); - // const auto raw_energy = 0.; double ecal_energy = 0.; - for(const auto t: ticl_cand.tracksters()) - { - double ecal_energy_fraction = t->raw_em_pt()/t->raw_pt(); + for (const auto t : ticl_cand.tracksters()) { + double ecal_energy_fraction = t->raw_em_pt() / t->raw_pt(); ecal_energy += t->raw_energy() * ecal_energy_fraction; } double hcal_energy = ticl_cand.rawEnergy() - ecal_energy; // fix for floating point rounding could go slightly below 0 - hcal_energy = hcal_energy < 0? 0:hcal_energy; + hcal_energy = hcal_energy < 0 ? 0 : hcal_energy; reco::PFCandidate::ParticleType part_type; switch (abs_pdg_id) { diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index b2a2165566789..d214365c53719 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -352,13 +352,12 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); - } tmpCandidate.setRawEnergy(raw_energy); math::XYZTLorentzVector p4(raw_energy * t.barycenter().unit().x(), - raw_energy * t.barycenter().unit().y(), - raw_energy * t.barycenter().unit().z(), - raw_energy); + raw_energy * t.barycenter().unit().y(), + raw_energy * t.barycenter().unit().z(), + raw_energy); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); @@ -376,9 +375,9 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es tmpCandidate.setRawEnergy(raw_energy); math::XYZTLorentzVector p4(raw_energy * t.barycenter().unit().x(), - raw_energy * t.barycenter().unit().y(), - raw_energy * t.barycenter().unit().z(), - raw_energy); + raw_energy * t.barycenter().unit().y(), + raw_energy * t.barycenter().unit().z(), + raw_energy); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); } @@ -404,15 +403,15 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es tmpCandidate.setPdgId(211 * track.charge()); for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); - raw_energy +=trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); + raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); } tmpCandidate.setRawEnergy(raw_energy); float momentum = std::sqrt(raw_energy * raw_energy - mpion2); const auto &barycenter = trackstersMergedHandle->at(closestTrackster).barycenter(); math::XYZTLorentzVector p4(momentum * barycenter.unit().x(), - momentum * barycenter.unit().y(), - momentum * barycenter.unit().z(), - t.raw_energy()); + momentum * barycenter.unit().y(), + momentum * barycenter.unit().z(), + t.raw_energy()); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); @@ -425,14 +424,14 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es tmpCandidate.setTrackPtr(edm::Ptr(track_h, trackIdx)); for (auto otherTracksterIdx : trackstersTRKwithSameSeed) { tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, otherTracksterIdx)); - raw_energy +=trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); + raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); } tmpCandidate.setPdgId(11 * track.charge()); tmpCandidate.setRawEnergy(raw_energy); math::XYZTLorentzVector p4(raw_energy * t.barycenter().unit().x(), - raw_energy * t.barycenter().unit().y(), - raw_energy * t.barycenter().unit().z(), - raw_energy); + raw_energy * t.barycenter().unit().y(), + raw_energy * t.barycenter().unit().z(), + raw_energy); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); } @@ -477,8 +476,10 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es tmpCandidate.setRawEnergy(t.raw_energy()); float momentum = std::sqrt(t.raw_energy() * t.raw_energy() - mpion2); - math::XYZTLorentzVector p4( - momentum * barycenter.unit().x(), momentum * barycenter.unit().y(), momentum * barycenter.unit().z(), t.raw_energy()); + math::XYZTLorentzVector p4(momentum * barycenter.unit().x(), + momentum * barycenter.unit().y(), + momentum * barycenter.unit().z(), + t.raw_energy()); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); } From 8a48807a91ee6939084368f948bc32d16ae543d3 Mon Sep 17 00:00:00 2001 From: Felice Date: Sun, 1 Nov 2020 11:22:25 +0100 Subject: [PATCH 61/66] assign track momentum direction at point of closest approach instead of trackster's barycenter in case of charged particles --- .../TICL/plugins/TrackstersMergeProducer.cc | 35 +++++++++---------- 1 file changed, 16 insertions(+), 19 deletions(-) diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index d214365c53719..d2143b17d618c 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -354,9 +354,9 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es raw_energy += trackstersMergedHandle->at(otherTracksterIdx).raw_energy(); } tmpCandidate.setRawEnergy(raw_energy); - math::XYZTLorentzVector p4(raw_energy * t.barycenter().unit().x(), - raw_energy * t.barycenter().unit().y(), - raw_energy * t.barycenter().unit().z(), + math::XYZTLorentzVector p4(raw_energy * track.momentum().unit().x(), + raw_energy * track.momentum().unit().y(), + raw_energy * track.momentum().unit().z(), raw_energy); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); @@ -374,9 +374,9 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es tmpCandidate.setPdgId(11 * track.charge()); tmpCandidate.setRawEnergy(raw_energy); - math::XYZTLorentzVector p4(raw_energy * t.barycenter().unit().x(), - raw_energy * t.barycenter().unit().y(), - raw_energy * t.barycenter().unit().z(), + math::XYZTLorentzVector p4(raw_energy * track.momentum().unit().x(), + raw_energy * track.momentum().unit().y(), + raw_energy * track.momentum().unit().z(), raw_energy); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); @@ -407,11 +407,10 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es } tmpCandidate.setRawEnergy(raw_energy); float momentum = std::sqrt(raw_energy * raw_energy - mpion2); - const auto &barycenter = trackstersMergedHandle->at(closestTrackster).barycenter(); - math::XYZTLorentzVector p4(momentum * barycenter.unit().x(), - momentum * barycenter.unit().y(), - momentum * barycenter.unit().z(), - t.raw_energy()); + math::XYZTLorentzVector p4(momentum * track.momentum().unit().x(), + momentum * track.momentum().unit().y(), + momentum * track.momentum().unit().z(), + raw_energy); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); @@ -428,9 +427,9 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es } tmpCandidate.setPdgId(11 * track.charge()); tmpCandidate.setRawEnergy(raw_energy); - math::XYZTLorentzVector p4(raw_energy * t.barycenter().unit().x(), - raw_energy * t.barycenter().unit().y(), - raw_energy * t.barycenter().unit().z(), + math::XYZTLorentzVector p4(raw_energy * track.momentum().unit().x(), + raw_energy * track.momentum().unit().y(), + raw_energy * track.momentum().unit().z(), raw_energy); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); @@ -467,7 +466,6 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es if (!usedSeeds[trackIdx]) { usedSeeds[trackIdx] = true; usedTrackstersMerged[mergedIdx] = true; - const auto &barycenter = t.barycenter(); TICLCandidate tmpCandidate; tmpCandidate.addTrackster(edm::Ptr(trackstersMergedHandle, mergedIdx)); tmpCandidate.setCharge(track.charge()); @@ -475,10 +473,9 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es tmpCandidate.setPdgId(211 * track.charge()); tmpCandidate.setRawEnergy(t.raw_energy()); float momentum = std::sqrt(t.raw_energy() * t.raw_energy() - mpion2); - - math::XYZTLorentzVector p4(momentum * barycenter.unit().x(), - momentum * barycenter.unit().y(), - momentum * barycenter.unit().z(), + math::XYZTLorentzVector p4(momentum * track.momentum().unit().x(), + momentum * track.momentum().unit().y(), + momentum * track.momentum().unit().z(), t.raw_energy()); tmpCandidate.setP4(p4); resultCandidates->push_back(tmpCandidate); From abc84da0735ccd850b313111396a8ec7361f4423 Mon Sep 17 00:00:00 2001 From: Felice Date: Sun, 1 Nov 2020 12:25:01 +0100 Subject: [PATCH 62/66] avoid copy --- RecoHGCal/TICL/plugins/PFTICLProducer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/PFTICLProducer.cc b/RecoHGCal/TICL/plugins/PFTICLProducer.cc index 09d244908330d..1c9f9b35035ed 100644 --- a/RecoHGCal/TICL/plugins/PFTICLProducer.cc +++ b/RecoHGCal/TICL/plugins/PFTICLProducer.cc @@ -54,7 +54,7 @@ void PFTICLProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSet const auto& four_mom = ticl_cand.p4(); double ecal_energy = 0.; - for (const auto t : ticl_cand.tracksters()) { + for (const auto& t : ticl_cand.tracksters()) { double ecal_energy_fraction = t->raw_em_pt() / t->raw_pt(); ecal_energy += t->raw_energy() * ecal_energy_fraction; } From 2443c12ee4c5c0ce20db8522cbec55c648501b8a Mon Sep 17 00:00:00 2001 From: Felice Date: Mon, 2 Nov 2020 17:58:17 +0100 Subject: [PATCH 63/66] Add missing const. --- RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index d2143b17d618c..ca06c277bb05b 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -438,7 +438,7 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es for (auto otherTracksterIdx : trackstersTRKEMwithSameSeed) { auto tmpIndex = (otherTracksterIdx != closestTrackster) ? otherTracksterIdx : mergedIdx; TICLCandidate photonCandidate; - auto &otherTrackster = trackstersMergedHandle->at(tmpIndex); + const auto &otherTrackster = trackstersMergedHandle->at(tmpIndex); auto gammaEnergy = otherTrackster.raw_energy(); photonCandidate.setCharge(0); photonCandidate.setPdgId(22); From a4029c0bcd1a23882b89d608bf6e25fac6bba871 Mon Sep 17 00:00:00 2001 From: Felice Date: Mon, 2 Nov 2020 19:08:23 +0100 Subject: [PATCH 64/66] Make track cuts configurable --- .../TICL/plugins/TrackstersMergeProducer.cc | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index ca06c277bb05b..34fa47485ee74 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -54,6 +54,10 @@ class TrackstersMergeProducer : public edm::stream::EDProducer("pt_sigma_high")), pt_sigma_low_(ps.getParameter("pt_sigma_low")), halo_max_distance2_(ps.getParameter("halo_max_distance2")), + track_min_pt_(ps.getParameter("track_min_pt")), + track_min_eta_(ps.getParameter("track_min_eta")), + track_max_eta_(ps.getParameter("track_max_eta")), + track_max_missing_outerhits_(ps.getParameter("track_max_missing_outerhits")), cosangle_align_(ps.getParameter("cosangle_align")), e_over_h_threshold_(ps.getParameter("e_over_h_threshold")), pt_neutral_threshold_(ps.getParameter("pt_neutral_threshold")), @@ -499,11 +507,12 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es usedSeeds[s.index] = true; } } + // for all general tracks (high purity, pt > 1), check if they have been used: if not, promote them as charged hadrons for (unsigned i = 0; i < tracks.size(); ++i) { auto const &track = tracks[i]; - if (track.pt() > 1.f and track.quality(reco::TrackBase::highPurity) and track.missingOuterHits() < 5 and - std::abs(track.outerEta()) > 1.48 and std::abs(track.outerEta()) < 3.0 and usedSeeds[i] == false) { + if (track.pt() > track_min_pt_ and track.quality(reco::TrackBase::highPurity) and track_max_missing_outerhits_ < 5 and + std::abs(track.outerEta()) > track_min_eta_ and std::abs(track.outerEta()) < track_max_eta_ and usedSeeds[i] == false) { // emit a charged hadron TICLCandidate tmpCandidate; tmpCandidate.setCharge(track.charge()); @@ -723,6 +732,10 @@ void TrackstersMergeProducer::fillDescriptions(edm::ConfigurationDescriptions &d desc.add("pt_sigma_high", 2.); desc.add("pt_sigma_low", 2.); desc.add("halo_max_distance2", 4.); + desc.add("track_min_pt", 1.); + desc.add("track_min_eta", 1.48); + desc.add("track_max_eta", 3.); + desc.add("track_max_missing_outerhits", 5); desc.add("cosangle_align", 0.9945); desc.add("e_over_h_threshold", 1.); desc.add("pt_neutral_threshold", 2.); From bbbacf1c98af9f3e974b1291f1ef521921c8fb7d Mon Sep 17 00:00:00 2001 From: Felice Date: Mon, 2 Nov 2020 19:11:40 +0100 Subject: [PATCH 65/66] code format --- RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index 34fa47485ee74..b6ea44679a759 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -511,8 +511,9 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es // for all general tracks (high purity, pt > 1), check if they have been used: if not, promote them as charged hadrons for (unsigned i = 0; i < tracks.size(); ++i) { auto const &track = tracks[i]; - if (track.pt() > track_min_pt_ and track.quality(reco::TrackBase::highPurity) and track_max_missing_outerhits_ < 5 and - std::abs(track.outerEta()) > track_min_eta_ and std::abs(track.outerEta()) < track_max_eta_ and usedSeeds[i] == false) { + if (track.pt() > track_min_pt_ and track.quality(reco::TrackBase::highPurity) and + track_max_missing_outerhits_ < 5 and std::abs(track.outerEta()) > track_min_eta_ and + std::abs(track.outerEta()) < track_max_eta_ and usedSeeds[i] == false) { // emit a charged hadron TICLCandidate tmpCandidate; tmpCandidate.setCharge(track.charge()); From 92c59aac9e13ac34526f972536a03747009b0680 Mon Sep 17 00:00:00 2001 From: Felice Date: Tue, 3 Nov 2020 11:30:00 +0100 Subject: [PATCH 66/66] fix configurable parameter --- RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index b6ea44679a759..d990c8291efa9 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -512,7 +512,7 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es for (unsigned i = 0; i < tracks.size(); ++i) { auto const &track = tracks[i]; if (track.pt() > track_min_pt_ and track.quality(reco::TrackBase::highPurity) and - track_max_missing_outerhits_ < 5 and std::abs(track.outerEta()) > track_min_eta_ and + track.missingOuterHits() < track_max_missing_outerhits_ and std::abs(track.outerEta()) > track_min_eta_ and std::abs(track.outerEta()) < track_max_eta_ and usedSeeds[i] == false) { // emit a charged hadron TICLCandidate tmpCandidate;