Skip to content

Commit

Permalink
Merge pull request #2 from amartelli/timeForMerged_11_2_X_forLeo
Browse files Browse the repository at this point in the history 
Time for merged tracksters
  • Loading branch information
@lecriste
lecriste authored Nov 10, 2020
2 parents 28193e0 + df5983d commit 6f18708
Show file tree
Hide file tree
Showing 4 changed files with 38 additions and 30 deletions.
33 changes: 5 additions & 28 deletions RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -117,32 +116,13 @@ void PatternRecognitionbyCA<TILES>::makeTracksters(
tracksterId++;

std::set<unsigned int> effective_cluster_idx;
std::pair<std::set<unsigned int>::iterator, bool> retVal;

std::vector<float> times;
std::vector<float> 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<TILES>::algo_verbosity_ > PatternRecognitionAlgoBaseT<TILES>::Advanced) {
LogDebug("HGCPatternRecoByCA") << " New doublet " << doublet << " for trackster: " << result.size()
Expand Down Expand Up @@ -201,16 +181,12 @@ void PatternRecognitionbyCA<TILES>::makeTracksters(
//if a seeding region does not lead to any trackster
tmp.setSeed(input.regions[0].collectionID, seedIndices[tracksterId]);

std::pair<float, float> 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);
Expand Down Expand Up @@ -268,7 +244,8 @@ void PatternRecognitionbyCA<TILES>::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);

Expand Down
9 changes: 8 additions & 1 deletion RecoHGCal/TICL/plugins/TrackstersMergeProducer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ class TrackstersMergeProducer : public edm::stream::EDProducer<edm::GlobalCache<
const edm::EDGetTokenT<std::vector<Trackster>> trackstershad_token_;
const edm::EDGetTokenT<std::vector<TICLSeedingRegion>> seedingTrk_token_;
const edm::EDGetTokenT<std::vector<reco::CaloCluster>> clusters_token_;
const edm::EDGetTokenT<edm::ValueMap<std::pair<float, float>>> clustersTime_token_;
const edm::EDGetTokenT<std::vector<reco::Track>> tracks_token_;
const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometry_token_;
const bool optimiseAcrossTracksters_;
Expand Down Expand Up @@ -87,6 +88,7 @@ TrackstersMergeProducer::TrackstersMergeProducer(const edm::ParameterSet &ps, co
trackstershad_token_(consumes<std::vector<Trackster>>(ps.getParameter<edm::InputTag>("trackstershad"))),
seedingTrk_token_(consumes<std::vector<TICLSeedingRegion>>(ps.getParameter<edm::InputTag>("seedingTrk"))),
clusters_token_(consumes<std::vector<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("layer_clusters"))),
clustersTime_token_(consumes<edm::ValueMap<std::pair<float, float>>>(ps.getParameter<edm::InputTag>("layer_clustersTime"))),
tracks_token_(consumes<std::vector<reco::Track>>(ps.getParameter<edm::InputTag>("tracks"))),
geometry_token_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
optimiseAcrossTracksters_(ps.getParameter<bool>("optimiseAcrossTracksters")),
Expand Down Expand Up @@ -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<edm::ValueMap<std::pair<float, float>>> clustersTime_h;
evt.getByToken(clustersTime_token_, clustersTime_h);
const auto& layerClustersTimes = *clustersTime_h;

edm::Handle<std::vector<Trackster>> trackstersem_h;
evt.getByToken(trackstersem_token_, trackstersem_h);
const auto &trackstersEM = *trackstersem_h;
Expand Down Expand Up @@ -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");
Expand Down Expand Up @@ -750,6 +756,7 @@ void TrackstersMergeProducer::fillDescriptions(edm::ConfigurationDescriptions &d
desc.add<edm::InputTag>("trackstershad", edm::InputTag("ticlTrackstersHAD"));
desc.add<edm::InputTag>("seedingTrk", edm::InputTag("ticlSeedingTrk"));
desc.add<edm::InputTag>("layer_clusters", edm::InputTag("hgcalLayerClusters"));
desc.add<edm::InputTag>("layer_clustersTime", edm::InputTag("hgcalLayerClusters", "timeLayerCluster"));
desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
desc.add<bool>("optimiseAcrossTracksters", true);
desc.add<int>("eta_bin_window", 1);
Expand Down
22 changes: 22 additions & 0 deletions RecoHGCal/TICL/plugins/TrackstersPCA.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,14 @@
#include "TPrincipal.h"

#include <iostream>
#include <set>

#include <Eigen/Core>
#include <Eigen/Dense>

void ticl::assignPCAtoTracksters(std::vector<Trackster> &tracksters,
const std::vector<reco::CaloCluster> &layerClusters,
const edm::ValueMap<std::pair<float, float>> &layerClustersTime,
double z_limit_em,
bool energyWeight) {
LogDebug("TrackstersPCA_Eigen") << "------- Eigen -------" << std::endl;
Expand Down Expand Up @@ -40,6 +42,10 @@ void ticl::assignPCAtoTracksters(std::vector<Trackster> &tracksters,
sigmasEigen << 0., 0., 0.;
Eigen::Matrix3d covM = Eigen::Matrix3d::Zero();

std::vector<float> times;
std::vector<float> timeErrors;
std::set<uint32_t> 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);
Expand All @@ -52,10 +58,24 @@ void ticl::assignPCAtoTracksters(std::vector<Trackster> &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<float, float> 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.
Expand Down Expand Up @@ -100,6 +120,7 @@ void ticl::assignPCAtoTracksters(std::vector<Trackster> &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);

Expand All @@ -110,6 +131,7 @@ void ticl::assignPCAtoTracksters(std::vector<Trackster> &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;
Expand Down
4 changes: 3 additions & 1 deletion RecoHGCal/TICL/plugins/TrackstersPCA.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <vector>

namespace ticl {
void assignPCAtoTracksters(std::vector<Trackster> &,
const std::vector<reco::CaloCluster> &,
const edm::ValueMap<std::pair<float, float>> &,
double,
bool energyWeight = true);
}
Expand Down

0 comments on commit 6f18708

Please sign in to comment.