diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc index 8340aa261fea8..9235b05aabdd6 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc @@ -7,7 +7,6 @@ #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/Utilities/interface/Exception.h" #include "PatternRecognitionbyCA.h" -#include "RecoLocalCalo/HGCalRecProducers/interface/ComputeClusterTime.h" #include "TrackstersPCA.h" #include "Geometry/CaloGeometry/interface/CaloGeometry.h" @@ -117,32 +116,13 @@ void PatternRecognitionbyCA::makeTracksters( tracksterId++; std::set effective_cluster_idx; - std::pair::iterator, bool> retVal; - - std::vector times; - std::vector timeErrors; for (auto const &doublet : ntuplet) { auto innerCluster = doublets[doublet].innerClusterId(); auto outerCluster = doublets[doublet].outerClusterId(); - retVal = effective_cluster_idx.insert(innerCluster); - if (retVal.second) { - float time = input.layerClustersTime.get(innerCluster).first; - if (time > -99) { - times.push_back(time); - timeErrors.push_back(1. / pow(input.layerClustersTime.get(innerCluster).second, 2)); - } - } - - retVal = effective_cluster_idx.insert(outerCluster); - if (retVal.second) { - float time = input.layerClustersTime.get(outerCluster).first; - if (time > -99) { - times.push_back(time); - timeErrors.push_back(1. / pow(input.layerClustersTime.get(outerCluster).second, 2)); - } - } + effective_cluster_idx.insert(innerCluster); + effective_cluster_idx.insert(outerCluster); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { LogDebug("HGCPatternRecoByCA") << " New doublet " << doublet << " for trackster: " << result.size() @@ -201,16 +181,12 @@ void PatternRecognitionbyCA::makeTracksters( //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()); + tmpTracksters, input.layerClusters, input.layerClustersTime, rhtools_.getPositionLayer(rhtools_.lastLayerEE(type)).z()); // run energy regression and ID energyRegressionAndID(input.layerClusters, tmpTracksters); @@ -268,7 +244,8 @@ void PatternRecognitionbyCA::makeTracksters( tmp.swap(result); } - ticl::assignPCAtoTracksters(result, input.layerClusters, rhtools_.getPositionLayer(rhtools_.lastLayerEE(type)).z()); + ticl::assignPCAtoTracksters(result, input.layerClusters, input.layerClustersTime, rhtools_.getPositionLayer(rhtools_.lastLayerEE(type)).z()); + // run energy regression and ID energyRegressionAndID(input.layerClusters, result); diff --git a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc index 8f86de7863743..1a6268affcc4e 100644 --- a/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc @@ -47,6 +47,7 @@ class TrackstersMergeProducer : public edm::stream::EDProducer> trackstershad_token_; const edm::EDGetTokenT> seedingTrk_token_; const edm::EDGetTokenT> clusters_token_; + const edm::EDGetTokenT>> clustersTime_token_; const edm::EDGetTokenT> tracks_token_; const edm::ESGetToken geometry_token_; const bool optimiseAcrossTracksters_; @@ -87,6 +88,7 @@ TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps, co trackstershad_token_(consumes>(ps.getParameter("trackstershad"))), seedingTrk_token_(consumes>(ps.getParameter("seedingTrk"))), clusters_token_(consumes>(ps.getParameter("layer_clusters"))), + clustersTime_token_(consumes>>(ps.getParameter("layer_clustersTime"))), tracks_token_(consumes>(ps.getParameter("tracks"))), geometry_token_(esConsumes()), optimiseAcrossTracksters_(ps.getParameter("optimiseAcrossTracksters")), @@ -188,6 +190,10 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es evt.getByToken(clusters_token_, cluster_h); const auto &layerClusters = *cluster_h; + edm::Handle>> clustersTime_h; + evt.getByToken(clustersTime_token_, clustersTime_h); + const auto& layerClustersTimes = *clustersTime_h; + edm::Handle> trackstersem_h; evt.getByToken(trackstersem_token_, trackstersem_h); const auto &trackstersEM = *trackstersem_h; @@ -257,7 +263,7 @@ void TrackstersMergeProducer::produce(edm::Event &evt, const edm::EventSetup &es iterMergedTracksters.push_back(TracksterIterIndex::HAD); } - assignPCAtoTracksters(*resultTrackstersMerged, layerClusters, rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z()); + assignPCAtoTracksters(*resultTrackstersMerged, layerClusters, layerClustersTimes, rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z()); energyRegressionAndID(layerClusters, *resultTrackstersMerged); printTrackstersDebug(*resultTrackstersMerged, "TrackstersMergeProducer"); @@ -750,6 +756,7 @@ void TrackstersMergeProducer::fillDescriptions(edm::ConfigurationDescriptions &d desc.add("trackstershad", edm::InputTag("ticlTrackstersHAD")); desc.add("seedingTrk", edm::InputTag("ticlSeedingTrk")); desc.add("layer_clusters", edm::InputTag("hgcalLayerClusters")); + desc.add("layer_clustersTime", edm::InputTag("hgcalLayerClusters", "timeLayerCluster")); desc.add("tracks", edm::InputTag("generalTracks")); desc.add("optimiseAcrossTracksters", true); desc.add("eta_bin_window", 1); diff --git a/RecoHGCal/TICL/plugins/TrackstersPCA.cc b/RecoHGCal/TICL/plugins/TrackstersPCA.cc index 9a91400195da1..539e1221ccee7 100644 --- a/RecoHGCal/TICL/plugins/TrackstersPCA.cc +++ b/RecoHGCal/TICL/plugins/TrackstersPCA.cc @@ -3,12 +3,14 @@ #include "TPrincipal.h" #include +#include #include #include void ticl::assignPCAtoTracksters(std::vector &tracksters, const std::vector &layerClusters, + const edm::ValueMap> &layerClustersTime, double z_limit_em, bool energyWeight) { LogDebug("TrackstersPCA_Eigen") << "------- Eigen -------" << std::endl; @@ -40,6 +42,10 @@ void ticl::assignPCAtoTracksters(std::vector &tracksters, sigmasEigen << 0., 0., 0.; Eigen::Matrix3d covM = Eigen::Matrix3d::Zero(); + std::vector times; + std::vector timeErrors; + std::set usedLC; + for (size_t i = 0; i < N; ++i) { auto fraction = 1.f / trackster.vertex_multiplicity(i); trackster.addToRawEnergy(layerClusters[trackster.vertices(i)].energy() * fraction); @@ -52,10 +58,24 @@ void ticl::assignPCAtoTracksters(std::vector &tracksters, fillPoint(layerClusters[trackster.vertices(i)], weight); for (size_t j = 0; j < 3; ++j) barycenter[j] += point[j]; + + // Also compute timing + auto before = usedLC.size(); + usedLC.insert(trackster.vertices(i)); + if(before + 1 == usedLC.size()){ + float timeE = layerClustersTime.get(trackster.vertices(i)).second; + if (timeE > -1.) { + times.push_back(layerClustersTime.get(trackster.vertices(i)).first); + timeErrors.push_back(1. / pow(timeE, 2)); + } + } } if (energyWeight && trackster.raw_energy()) barycenter /= trackster.raw_energy(); + hgcalsimclustertime::ComputeClusterTime timeEstimator; + std::pair timeTrackster = timeEstimator.fixSizeHighestDensity(times, timeErrors); + // Compute the Covariance Matrix and the sum of the squared weights, used // to compute the correct normalization. // The barycenter has to be known. @@ -100,6 +120,7 @@ void ticl::assignPCAtoTracksters(std::vector &tracksters, // Add trackster attributes trackster.setBarycenter(ticl::Trackster::Vector(barycenter)); + trackster.setTimeAndError(timeTrackster.first, timeTrackster.second); trackster.fillPCAVariables( eigenvalues_fromEigen, eigenvectors_fromEigen, sigmas, sigmasEigen, 3, ticl::Trackster::PCAOrdering::ascending); @@ -110,6 +131,7 @@ void ticl::assignPCAtoTracksters(std::vector &tracksters, LogDebug("TrackstersPCA") << "raw_pt: " << trackster.raw_pt() << std::endl; LogDebug("TrackstersPCA") << "Means: " << barycenter[0] << ", " << barycenter[1] << ", " << barycenter[2] << std::endl; + LogDebug("TrackstersPCA") << "Time: " << trackster.time() << " +/- " << trackster.timeError() << std::endl; LogDebug("TrackstersPCA") << "EigenValues from Eigen/Tr(cov): " << eigenvalues_fromEigen[2] / covM.trace() << ", " << eigenvalues_fromEigen[1] / covM.trace() << ", " << eigenvalues_fromEigen[0] / covM.trace() << std::endl; diff --git a/RecoHGCal/TICL/plugins/TrackstersPCA.h b/RecoHGCal/TICL/plugins/TrackstersPCA.h index a95c3786413b1..819c812ca3945 100644 --- a/RecoHGCal/TICL/plugins/TrackstersPCA.h +++ b/RecoHGCal/TICL/plugins/TrackstersPCA.h @@ -3,12 +3,14 @@ #include "DataFormats/HGCalReco/interface/Trackster.h" #include "DataFormats/CaloRecHit/interface/CaloCluster.h" - +#include "DataFormats/Common/interface/ValueMap.h" +#include "RecoLocalCalo/HGCalRecProducers/interface/ComputeClusterTime.h" #include namespace ticl { void assignPCAtoTracksters(std::vector &, const std::vector &, + const edm::ValueMap> &, double, bool energyWeight = true); }