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/Configuration/python/RecoHGCal_EventContent_cff.py b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py index 3a16bea88e3b3..77b647e6ce291 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,16 +16,15 @@ #RECO content TICL_RECO = cms.PSet( outputCommands = cms.untracked.vstring( + 'keep *_ticlTrackstersTrkEM_*_*', 'keep *_ticlTrackstersEM_*_*', 'keep *_ticlTrackstersHAD_*_*', 'keep *_ticlTrackstersTrk_*_*', 'keep *_ticlTrackstersMIP_*_*', 'keep *_ticlTrackstersMerge_*_*', - 'keep *_ticlCandidateFromTracksters_*_*', 'keep *_ticlTrackstersHFNoseEM_*_*', 'keep *_ticlTrackstersHFNoseMIP_*_*', 'keep *_ticlTrackstersHFNoseMerge_*_*', - 'keep *_ticlCandidateFromTrackstersHFNose_*_*', 'keep *_pfTICL_*_*' ) ) diff --git a/RecoHGCal/TICL/plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h b/RecoHGCal/TICL/plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h new file mode 100644 index 0000000000000..5f1d9388c3fe4 --- /dev/null +++ b/RecoHGCal/TICL/plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h @@ -0,0 +1,55 @@ +// 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/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); } diff --git a/RecoHGCal/TICL/plugins/HGCGraph.cc b/RecoHGCal/TICL/plugins/HGCGraph.cc index 1599d6c98df5d..e94b87426bc95 100644 --- a/RecoHGCal/TICL/plugins/HGCGraph.cc +++ b/RecoHGCal/TICL/plugins/HGCGraph.cc @@ -19,14 +19,18 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, int deltaIPhi, float minCosTheta, float minCosPointing, + float root_doublet_max_distance_from_seed_squared, float etaLimitIncreaseWindow, - int missing_layers, + int skip_layers, int maxNumberOfLayers, float maxDeltaTime) { isOuterClusterOfDoublets_.clear(); isOuterClusterOfDoublets_.resize(layerClusters.size()); allDoublets_.clear(); theRootDoublets_.clear(); + bool checkDistanceRootDoubletVsSeed = root_doublet_max_distance_from_seed_squared < 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.) @@ -75,7 +82,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]; @@ -171,8 +178,17 @@ void HGCGraphT::makeAndConnectDoublets(const TILES &histo, minCosTheta, minCosPointing, verbosity_ > Advanced); - if (isRootDoublet) + if (isRootDoublet and checkDistanceRootDoubletVsSeed) { + 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) { theRootDoublets_.push_back(doubletId); + } } } } diff --git a/RecoHGCal/TICL/plugins/HGCGraph.h b/RecoHGCal/TICL/plugins/HGCGraph.h index 103f76b312300..f8cc626016f40 100644 --- a/RecoHGCal/TICL/plugins/HGCGraph.h +++ b/RecoHGCal/TICL/plugins/HGCGraph.h @@ -25,8 +25,9 @@ class HGCGraphT { int deltaIPhi, float minCosThetai, float maxCosPointing, + float root_doublet_max_distance_from_seed_squared, float etaLimitIncreaseWindow, - int missing_layers, + int skip_layers, int maxNumberOfLayers, float maxDeltaTime); diff --git a/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc b/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc index 195843d0015e8..27cf92d6c1265 100644 --- a/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc @@ -63,13 +63,14 @@ 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.}}; + 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,13 @@ 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)); diff --git a/RecoHGCal/TICL/plugins/PFTICLProducer.cc b/RecoHGCal/TICL/plugins/PFTICLProducer.cc index 443340944b39b..1c9f9b35035ed 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); } @@ -52,6 +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(); + 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 +87,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. diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index dcc0c039d15db..8340aa261fea8 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -26,9 +26,19 @@ 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")), etaLimitIncreaseWindow_(conf.getParameter("etaLimitIncreaseWindow")), - missing_layers_(conf.getParameter("missing_layers")), - min_clusters_per_ntuplet_(conf.getParameter("min_clusters_per_ntuplet")), + 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")), + 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")), eidOutputNameEnergy_(conf.getParameter("eid_output_name_energy")), @@ -88,17 +98,24 @@ void PatternRecognitionbyCA::makeTracksters( 1, min_cos_theta_, min_cos_pointing_, + root_doublet_max_distance_from_seed_squared_, etaLimitIncreaseWindow_, - missing_layers_, + skip_layers_, rhtools_.lastLayer(type), max_delta_time_); 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; @@ -137,29 +154,102 @@ void PatternRecognitionbyCA::makeTracksters( << input.layerClusters[outerCluster].z() << " " << tracksterId << std::endl; } } + unsigned showerMinLayerId = 99999; + 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) { - layer_cluster_usage[i]++; + 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); + } + 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_) { + int numberOfMissingLayers = 0; + unsigned int j = showerMinLayerId; + unsigned int indexInVec = 0; + for (const auto &layer : uniqueLayerIds) { + if (layer != j) { + numberOfMissingLayers++; + j++; + if (numberOfMissingLayers > max_missing_layers_in_trackster_) { + numberOfLayersInTrackster = indexInVec; + for (auto &llpair : lcIdAndLayer) { + if (llpair.second >= layer) { + effective_cluster_idx.erase(llpair.first); + } + } + break; + } + } + indexInVec++; + j++; + } + } + + 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()); + 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]); + + 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())); + tmpTracksters.push_back(tmp); + } + } + 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() < energy_em_over_total_threshold_ * t.raw_energy()); + }; + + 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_) { + selectedTrackstersIds.push_back(i); + } + } + + result.reserve(selectedTrackstersIds.size()); + + for (unsigned i = 0; i < selectedTrackstersIds.size(); ++i) { + const auto &t = tmpTracksters[selectedTrackstersIds[i]]; + for (auto const lcId : t.vertices()) { + layer_cluster_usage[lcId]++; if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Basic) - LogDebug("HGCPatternRecoByCA") << "LayerID: " << i << " count: " << (int)layer_cluster_usage[i] << std::endl; + LogDebug("HGCPatternRecoByCA") << "LayerID: " << lcId << " count: " << (int)layer_cluster_usage[lcId] + << 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); + seedToTracksterAssociation[t.seedIndex()].push_back(i); } - - 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++; + result.push_back(t); } for (auto &trackster : result) { @@ -171,7 +261,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; @@ -180,7 +269,6 @@ void PatternRecognitionbyCA::makeTracksters( } 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 f6b917cdb2899..c337dd344ec78 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.h @@ -37,8 +37,17 @@ 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_squared_; const float etaLimitIncreaseWindow_; - const int missing_layers_; + 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 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/SeedingRegionByTracks.cc b/RecoHGCal/TICL/plugins/SeedingRegionByTracks.cc index 45ecb51d13907..16d69c391f38a 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.perp2() > b.directionAtOrigin.perp2(); + }); } void SeedingRegionByTracks::buildFirstLayers() { diff --git a/RecoHGCal/TICL/plugins/TICLSeedingRegionProducer.cc b/RecoHGCal/TICL/plugins/TICLSeedingRegionProducer.cc index 152c7b7898eb4..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\") < 10"); + "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5"); desc.add("propagator", "PropagatorWithMaterial"); desc.add("algoId", 1); descriptions.add("ticlSeedingRegionProducer", desc); diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index 290eb66a605b5..d990c8291efa9 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -1,5 +1,4 @@ #include // unique_ptr - #include "FWCore/Framework/interface/stream/EDProducer.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" @@ -15,6 +14,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,13 +33,14 @@ class TrackstersMergeProducer : public edm::stream::EDProducer &, TracksterIterIndex); void printTrackstersDebug(const std::vector &, 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_; @@ -52,11 +53,18 @@ 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"))), @@ -84,11 +93,18 @@ TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps, co phi_bin_window_(ps.getParameter("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")), + 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")), - 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")), @@ -106,6 +122,7 @@ TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps, co eidSession_ = tensorflow::createSession(trackstersCache->eidGraphDef); produces>(); + produces>(); } void TrackstersMergeProducer::fillTile(TICLTracksterTiles &tracksterTile, @@ -148,15 +165,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 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; @@ -168,185 +190,345 @@ 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; 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; + const 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; + const 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]) { + const 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)); + 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(); + } + tmpCandidate.setRawEnergy(raw_energy); + 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); + + } 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()); + + tmpCandidate.setRawEnergy(raw_energy); + 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); + } + + } 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; + } + 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(); + } + tmpCandidate.setRawEnergy(raw_energy); + float momentum = std::sqrt(raw_energy * raw_energy - mpion2); + 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); + + } 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()); + tmpCandidate.setRawEnergy(raw_energy); + 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); + } + // Promote all other TRKEM tracksters as photons with their energy. + for (auto otherTracksterIdx : trackstersTRKEMwithSameSeed) { + auto tmpIndex = (otherTracksterIdx != closestTrackster) ? otherTracksterIdx : mergedIdx; + TICLCandidate photonCandidate; + const 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]; + const 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; + 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(t.raw_energy()); + float momentum = std::sqrt(t.raw_energy() * t.raw_energy() - mpion2); + 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); + } + } + } + // 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; } - tracksterHAD_idx++; } - assignPCAtoTracksters(*result, layerClusters, rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z()); - energyRegressionAndID(layerClusters, *result); - - printTrackstersDebug(*result, "TrackstersMergeProducer"); + // 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.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; + 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; + } + } - evt.put(std::move(result)); + evt.put(std::move(resultCandidates)); } void TrackstersMergeProducer::energyRegressionAndID(const std::vector &layerClusters, @@ -538,6 +720,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")); @@ -549,11 +732,18 @@ 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("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.); - 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"); diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc index c4d228b130258..38650eeab1211 100644 --- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc @@ -54,8 +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 std::string itername_; }; DEFINE_FWK_MODULE(TrackstersProducer); @@ -91,8 +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")), itername_(ps.getParameter("itername")) { if (doNose_) { layer_clusters_tiles_hfnose_token_ = @@ -116,13 +112,18 @@ 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("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); + desc.add("root_doublet_max_distance_from_seed_squared", 9999); desc.add("min_cos_pointing", -1.); - desc.add("missing_layers", 0); + 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); @@ -172,30 +173,12 @@ 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. - // 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_; - }; - - // 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) { for (auto const v : trackster.vertices()) { 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"); diff --git a/RecoHGCal/TICL/python/EMStep_cff.py b/RecoHGCal/TICL/python/EMStep_cff.py index 0c3a5b5c497f3..0aa6f9ca50fd9 100644 --- a/RecoHGCal/TICL/python/EMStep_cff.py +++ b/RecoHGCal/TICL/python/EMStep_cff.py @@ -8,10 +8,11 @@ # CLUSTER FILTERING/MASKING filteredLayerClustersEM = _filteredLayerClustersProducer.clone( - clusterFilter = "ClusterFilterByAlgoAndSize", - min_cluster_size = 2, # inclusive + clusterFilter = "ClusterFilterByAlgoAndSizeAndLayerRange", + min_cluster_size = 3, # inclusive + max_layerId = 30, # inclusive algo_number = 8, - LayerClustersInputMask = 'ticlTrackstersTrk', + LayerClustersInputMask = 'ticlTrackstersTrkEM', iteration_label = "EM" ) @@ -19,14 +20,18 @@ 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, - max_out_in_hops = 4, - missing_layers = 1, - min_clusters_per_ntuplet = 10, - min_cos_theta = 0.978, # ~12 degrees + 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, + max_missing_layers_in_trackster = 1, + min_layers_per_trackster = 10, + min_cos_theta = 0.97, # ~14 degrees min_cos_pointing = 0.9, # ~25 degrees max_delta_time = 3., itername = "EM", @@ -60,7 +65,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 883f84412b22b..72efaebc201ec 100644 --- a/RecoHGCal/TICL/python/HADStep_cff.py +++ b/RecoHGCal/TICL/python/HADStep_cff.py @@ -10,23 +10,23 @@ filteredLayerClustersHAD = _filteredLayerClustersProducer.clone( clusterFilter = "ClusterFilterByAlgoAndSize", - min_cluster_size = 2, # inclusive + min_cluster_size = 3, # 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 # pid_threshold = 0.7, - missing_layers = 1, - min_clusters_per_ntuplet = 12, + skip_layers = 1, + 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 da90839a6b90f..f7e5234b3aa9c 100644 --- a/RecoHGCal/TICL/python/MIPStep_cff.py +++ b/RecoHGCal/TICL/python/MIPStep_cff.py @@ -20,8 +20,8 @@ ticlTrackstersMIP = _trackstersProducer.clone( filtered_mask = "filteredLayerClustersMIP:MIP", seeding_regions = "ticlSeedingGlobal", - missing_layers = 3, - min_clusters_per_ntuplet = 10, + skip_layers = 3, + 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 new file mode 100644 index 0000000000000..458f489f696ec --- /dev/null +++ b/RecoHGCal/TICL/python/TrkEMStep_cff.py @@ -0,0 +1,51 @@ +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 = 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, + energy_em_over_total_threshold = 0.9, + 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 + min_cos_pointing = 0.94, # ~20 degrees + root_doublet_max_distance_from_seed_squared = 2.5e-3, # dR=0.05 + 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..739a7cafd5d2b 100644 --- a/RecoHGCal/TICL/python/TrkStep_cff.py +++ b/RecoHGCal/TICL/python/TrkStep_cff.py @@ -9,8 +9,10 @@ # CLUSTER FILTERING/MASKING filteredLayerClustersTrk = _filteredLayerClustersProducer.clone( - clusterFilter = "ClusterFilterByAlgo", + clusterFilter = "ClusterFilterByAlgoAndSize", + min_cluster_size = 3, # inclusive algo_number = 8, + LayerClustersInputMask = 'ticlTrackstersEM', iteration_label = "Trk" ) @@ -19,10 +21,11 @@ ticlTrackstersTrk = _trackstersProducer.clone( filtered_mask = "filteredLayerClustersTrk:Trk", seeding_regions = "ticlSeedingTrk", + original_mask = 'ticlTrackstersEM', filter_on_categories = [2, 4], # filter muons and charged hadrons pid_threshold = 0.0, - missing_layers = 3, - min_clusters_per_ntuplet = 10, + skip_layers = 3, + min_layers_per_trackster = 10, min_cos_theta = 0.866, # ~30 degrees min_cos_pointing = 0.798, # ~ 37 degrees max_delta_time = -1., @@ -43,3 +46,4 @@ ,ticlTrackstersTrk ,ticlMultiClustersFromTrackstersTrk) + diff --git a/RecoHGCal/TICL/python/iterativeTICL_cff.py b/RecoHGCal/TICL/python/iterativeTICL_cff.py index de31c0a2f4d3f..9af6dc68a9af6 100644 --- a/RecoHGCal/TICL/python/iterativeTICL_cff.py +++ b/RecoHGCal/TICL/python/iterativeTICL_cff.py @@ -1,11 +1,11 @@ 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 * 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 @@ -18,18 +18,14 @@ ) 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 - ,ticlMIPStepTask - ,ticlTrkStepTask + ,ticlTrkEMStepTask ,ticlEMStepTask + ,ticlTrkStepTask ,ticlHADStepTask ,ticlTracksterMergeTask ,ticlPFTask diff --git a/RecoHGCal/TICL/python/ticl_iterations.py b/RecoHGCal/TICL/python/ticl_iterations.py index b36d459990c70..6266c07807603 100644 --- a/RecoHGCal/TICL/python/ticl_iterations.py +++ b/RecoHGCal/TICL/python/ticl_iterations.py @@ -41,8 +41,8 @@ def TICL_iterations_withReco(process): process.trackstersTrk = trackstersProducer.clone( filtered_mask = "filteredLayerClustersTrk:Trk", seeding_regions = "ticlSeedingTrk", - missing_layers = 3, - min_clusters_per_ntuplet = 5, + skip_layers = 3, + min_layers_per_trackster = 5, min_cos_theta = 0.99, # ~10 degrees min_cos_pointing = 0.9 ) @@ -66,8 +66,8 @@ def TICL_iterations_withReco(process): process.trackstersMIP = trackstersProducer.clone( filtered_mask = "filteredLayerClustersMIP:MIP", seeding_regions = "ticlSeedingGlobal", - missing_layers = 3, - min_clusters_per_ntuplet = 15, + skip_layers = 3, + min_layers_per_trackster = 15, min_cos_theta = 0.99, # ~10 degrees min_cos_pointing = 0.9, out_in_dfs = False, @@ -91,8 +91,8 @@ def TICL_iterations_withReco(process): original_mask = "trackstersMIP", filtered_mask = "filteredLayerClusters:algo8", seeding_regions = "ticlSeedingGlobal", - missing_layers = 1, - min_clusters_per_ntuplet = 10, + skip_layers = 1, + min_layers_per_trackster = 10, min_cos_theta = 0.984, # ~10 degrees min_cos_pointing = 0.9 # ~26 degrees ) @@ -105,8 +105,8 @@ def TICL_iterations_withReco(process): process.trackstersHAD = trackstersProducer.clone( filtered_mask = "filteredLayerClusters:algo8", seeding_regions = "ticlSeedingGlobal", - missing_layers = 2, - min_clusters_per_ntuplet = 10, + skip_layers = 2, + min_layers_per_trackster = 10, min_cos_theta = 0.8, min_cos_pointing = 0.7 ) @@ -172,8 +172,8 @@ def TICL_iterations(process): process.trackstersMIP = trackstersProducer.clone( filtered_mask = "filteredLayerClustersMIP:MIP", seeding_regions = "ticlSeedingGlobal", - missing_layers = 3, - min_clusters_per_ntuplet = 15, + skip_layers = 3, + min_layers_per_trackster = 15, min_cos_theta = 0.99, # ~10 degrees ) @@ -193,8 +193,8 @@ def TICL_iterations(process): original_mask = "trackstersMIP", filtered_mask = "filteredLayerClusters:algo8", seeding_regions = "ticlSeedingGlobal", - missing_layers = 2, - min_clusters_per_ntuplet = 15, + skip_layers = 2, + min_layers_per_trackster = 15, min_cos_theta = 0.94, # ~20 degrees min_cos_pointing = 0.7 ) 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(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(); + 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 << 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; @@ -1979,37 +1974,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()); } } @@ -2052,6 +2050,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++; } @@ -2116,18 +2118,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 +2172,21 @@ void HGVHistoProducerAlgo::fill_multi_cluster_histos(const Histograms& histogram histograms.h_multiplicityOfLCinMCL_vs_layerclusterenergy[count]->Fill(mlp, layerClusters[lc]->energy()); } - 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()); + 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_energy[count]->Fill(multiClusters[mclId].energy()); + } } //end of loop through multiclusters