From 5bbda0b8990f115a434aaae643e3588917e0a79c Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Tue, 8 Feb 2022 22:03:15 +0100 Subject: [PATCH 01/25] Use absolute metric in CLUE3D --- .../plugins/PatternRecognitionbyCLUE3D.cc | 56 +++++++++++++++---- .../TICL/plugins/PatternRecognitionbyCLUE3D.h | 6 ++ 2 files changed, 51 insertions(+), 11 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 3ebf8b0517aca..4a20739139745 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -93,8 +93,8 @@ template void PatternRecognitionbyCLUE3D::dumpClusters(const std::vector> &layerIdx2layerandSoa, const int eventNumber) const { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "[evt, layer, x, y, eta, phi, cells, energy, radius, rho, delta, " - "isSeed, clusterIdx, layerClusterOriginalIdx"; + edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "[evt, layer, x, y, z, r_over_absz, eta, phi, cells, energy, radius, rho, delta, " + "isSeed, clusterIdx, layerClusterOriginalIdx, SOAidx"; } for (unsigned int layer = 0; layer < clusters_.size(); layer++) { @@ -107,10 +107,11 @@ void PatternRecognitionbyCLUE3D::dumpClusters(const std::vector::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "ClusterInfo: " << eventNumber << ", " << layer << ", " << v << ", " << thisLayer.y[num] << ", " + << thisLayer.z[num] << ", " << thisLayer.r_over_absz[num] << ", " << thisLayer.eta[num] << ", " << thisLayer.phi[num] << ", " << thisLayer.cells[num] << ", " << thisLayer.energy[num] << ", " << thisLayer.radius[num] << ", " << thisLayer.rho[num] << ", " << thisLayer.delta[num] << ", " << thisLayer.isSeed[num] << ", " << thisLayer.clusterIndex[num] << ", " - << thisLayer.layerClusterOriginalIdx[num]; + << thisLayer.layerClusterOriginalIdx[num] << ", " << num; } ++num; } @@ -154,7 +155,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( for (auto const &lc : input.layerClusters) { if (input.mask[layerIdx] == 0.) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Skipping masked clustrer: " << layerIdx; + edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Skipping masked cluster: " << layerIdx; } layerIdx2layerandSoa.emplace_back(-1, -1); layerIdx++; @@ -189,6 +190,8 @@ void PatternRecognitionbyCLUE3D::makeTracksters( float radius_y = sqrt((sum_sqr_y - (sum_y * sum_y) * invClsize) * invClsize); clusters_[layer].x.emplace_back(lc.x()); clusters_[layer].y.emplace_back(lc.y()); + clusters_[layer].z.emplace_back(lc.z()); + clusters_[layer].r_over_absz.emplace_back(sqrt(lc.x()*lc.x()+lc.y()*lc.y())/std::abs(lc.z())); clusters_[layer].radius.emplace_back(radius_x + radius_y); clusters_[layer].eta.emplace_back(lc.eta()); clusters_[layer].phi.emplace_back(lc.phi()); @@ -436,17 +439,23 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( auto isReachable = [&](float x1, float x2, float y1, float y2, float delta_sqr) -> bool { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecogntionbyCLUE3D") + << " is Reachable: " << ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)) << " vs " << delta_sqr << "[" << ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) < delta_sqr) << "]\n"; } return (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) < delta_sqr; }; + auto distance_debug = [&](float x1, float x2, float y1, float y2) -> float { + return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)); + }; + for (unsigned int i = 0; i < numberOfClusters; i++) { // We need to partition the two sides of the HGCAL detector auto lastLayerPerSide = static_cast(rhtools_.lastLayer(false)); int minLayer = 0; int maxLayer = 2 * lastLayerPerSide - 1; + const float thisClusterAbsZNormalization = 1./(clustersOnLayer.z[i]*clustersOnLayer.z[i]); if (layerId < lastLayerPerSide) { minLayer = std::max(layerId - densitySiblingLayers_, minLayer); maxLayer = std::min(layerId + densitySiblingLayers_, lastLayerPerSide - 1); @@ -523,10 +532,22 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( clustersLayer.eta[layerandSoa.second], clustersLayer.phi[layerandSoa.second]); } + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + auto dist = distance_debug(clustersOnLayer.r_over_absz[i], + clustersLayer.r_over_absz[layerandSoa.second], + clustersOnLayer.r_over_absz[i]*std::abs(clustersOnLayer.phi[i]), + clustersLayer.r_over_absz[layerandSoa.second]*std::abs(clustersLayer.phi[layerandSoa.second])); + edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Distance[cm]: " << (dist*clustersOnLayer.z[i]); + } + if (isReachable(clustersOnLayer.r_over_absz[i], + clustersLayer.r_over_absz[layerandSoa.second], + clustersOnLayer.r_over_absz[i]*std::abs(clustersOnLayer.phi[i]), + clustersLayer.r_over_absz[layerandSoa.second]*std::abs(clustersLayer.phi[layerandSoa.second]), thisClusterAbsZNormalization*2.6f*2.6f)) { + /* if (reco::deltaR2(clustersOnLayer.eta[i], clustersOnLayer.phi[i], clustersLayer.eta[layerandSoa.second], - clustersLayer.phi[layerandSoa.second]) < densityEtaPhiDistanceSqr_) { + clustersLayer.phi[layerandSoa.second]) < densityEtaPhiDistanceSqr_) {*/ auto energyToAdd = (clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx ? 1.f : 0.5f) * clustersLayer.energy[layerandSoa.second]; clustersOnLayer.rho[i] += energyToAdd; @@ -552,6 +573,10 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( auto &clustersOnLayer = clusters_[layerId]; unsigned int numberOfClusters = clustersOnLayer.x.size(); + auto distance_debug = [&](float x1, float x2, float y1, float y2) -> float { + return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)); + }; + for (unsigned int i = 0; i < numberOfClusters; i++) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecogntionbyCLUE3D") @@ -595,10 +620,19 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( if ((layerandSoa.first == -1) && (layerandSoa.second == -1)) continue; auto const &clustersOnOtherLayer = clusters_[layerandSoa.first]; - float dist = reco::deltaR2(clustersOnLayer.eta[i], + auto dist = distance_debug(clustersOnLayer.r_over_absz[i], + clustersOnOtherLayer.r_over_absz[layerandSoa.second], + clustersOnLayer.r_over_absz[i]*std::abs(clustersOnLayer.phi[i]), + clustersOnOtherLayer.r_over_absz[layerandSoa.second]*std::abs(clustersOnOtherLayer.phi[layerandSoa.second])); + // Rescale the distance in cm on the current layer + dist *=clustersOnLayer.z[i]; + // TODO(rovere): add also z to the distance metric? + /* + float dist = reco::deltaR2(clustersOnLayer.eta[i], clustersOnLayer.phi[i], clustersOnOtherLayer.eta[layerandSoa.second], clustersOnOtherLayer.phi[layerandSoa.second]); + */ bool foundHigher = (clustersOnOtherLayer.rho[layerandSoa.second] > clustersOnLayer.rho[i]) || (clustersOnOtherLayer.rho[layerandSoa.second] == clustersOnLayer.rho[i] && clustersOnOtherLayer.layerClusterOriginalIdx[layerandSoa.second] > @@ -608,7 +642,7 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( << "Searching nearestHigher on " << currentLayer << " with rho: " << clustersOnOtherLayer.rho[layerandSoa.second] << " on layerIdxInSOA: " << layerandSoa.first << ", " << layerandSoa.second - << " with distance: " << sqrt(dist) << " foundHigher: " << foundHigher; + << " with distance: " << dist << " foundHigher: " << foundHigher; } if (foundHigher && dist <= i_delta) { // update i_delta @@ -624,11 +658,11 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( bool foundNearestHigherInEtaPhiCylinder = (i_delta != maxDelta); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecogntionbyCLUE3D") - << "i_delta: " << sqrt(i_delta) << " passed: " << foundNearestHigherInEtaPhiCylinder << " " + << "i_delta: " << i_delta << " passed: " << foundNearestHigherInEtaPhiCylinder << " " << i_nearestHigher.first << " " << i_nearestHigher.second; } if (foundNearestHigherInEtaPhiCylinder) { - clustersOnLayer.delta[i] = sqrt(i_delta); + clustersOnLayer.delta[i] = i_delta; clustersOnLayer.nearestHigher[i] = i_nearestHigher; } else { // otherwise delta is guaranteed to be larger outlierDeltaFactor_*delta_c @@ -653,8 +687,8 @@ int PatternRecognitionbyCLUE3D::findAndAssignTracksters( // initialize clusterIndex clustersOnLayer.clusterIndex[i] = -1; bool isSeed = - (clustersOnLayer.delta[i] > criticalEtaPhiDistance_) && (clustersOnLayer.rho[i] >= criticalDensity_); - bool isOutlier = (clustersOnLayer.delta[i] > outlierMultiplier_ * criticalEtaPhiDistance_) && + (clustersOnLayer.delta[i] > 2.6) && (clustersOnLayer.rho[i] >= criticalDensity_); + bool isOutlier = (clustersOnLayer.delta[i] > outlierMultiplier_ * 2.6) && (clustersOnLayer.rho[i] < criticalDensity_); if (isSeed) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h index 3802273fd872b..da0efefd8380f 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h @@ -28,6 +28,8 @@ namespace ticl { struct ClustersOnLayer { std::vector x; std::vector y; + std::vector z; + std::vector r_over_absz; std::vector radius; std::vector eta; std::vector phi; @@ -46,6 +48,8 @@ namespace ticl { void clear() { x.clear(); y.clear(); + z.clear(); + r_over_absz.clear(); radius.clear(); eta.clear(); phi.clear(); @@ -63,6 +67,8 @@ namespace ticl { void shrink_to_fit() { x.shrink_to_fit(); y.shrink_to_fit(); + z.shrink_to_fit(); + r_over_absz.shrink_to_fit(); radius.shrink_to_fit(); eta.shrink_to_fit(); phi.shrink_to_fit(); From 213a135b70a8ee2c8ff5ac8f01fa4a824599cdb0 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Mon, 14 Feb 2022 11:29:29 +0100 Subject: [PATCH 02/25] Transverse and longitudinal critical distances In the computation of the nearest cluster with higher local density, the distance, in cm, along the Z direction is also taken into account. This should follow more closely the real 'energy-flow' within real showers, avoid the ping-pong case of connections going back and forth while slowly approaching the shower axis. The longitudinal separation, express in number of layers, is also taken into account while promoting a layercluster as a seed. The requirement is in 'or' with the traditional one on the transverse distance. A layercluster now, in order to be promoted as a seed, needs to be separated at least threshold_xy distance on the transverse plane from its nearest-higher OR separated at least threshold_layers from it. This will avoid excessive fragmentation of showers along the Z direction due to 'energy fluctuations' in the overall 'energy flow'. A new variable is now considered to promote a layer cluster as a seed. It represent the ratio between the energy of the layer cluster and the local density associated to it. A reasonable value for this variable should be set to 1/(densitySiblingLayers+1), assuming a somewhat uniform local density distribution. Minor improvements in the printouts. --- .../plugins/PatternRecognitionbyCLUE3D.cc | 238 ++++++++++-------- .../TICL/plugins/PatternRecognitionbyCLUE3D.h | 9 +- 2 files changed, 140 insertions(+), 107 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 4a20739139745..acfa961d2b799 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -22,10 +22,15 @@ PatternRecognitionbyCLUE3D::PatternRecognitionbyCLUE3D(const edm::Paramet : PatternRecognitionAlgoBaseT(conf, iC), caloGeomToken_(iC.esConsumes()), criticalDensity_(conf.getParameter("criticalDensity")), + criticalSelfDensity_(conf.getParameter("criticalSelfDensity")), densitySiblingLayers_(conf.getParameter("densitySiblingLayers")), densityEtaPhiDistanceSqr_(conf.getParameter("densityEtaPhiDistanceSqr")), + densityXYDistanceSqr_(conf.getParameter("densityXYDistanceSqr")), densityOnSameLayer_(conf.getParameter("densityOnSameLayer")), + useAbsoluteProjectiveScale_(conf.getParameter("useAbsoluteProjectiveScale")), criticalEtaPhiDistance_(conf.getParameter("criticalEtaPhiDistance")), + criticalXYDistance_(conf.getParameter("criticalXYDistance")), + criticalZDistanceLyr_(conf.getParameter("criticalZDistanceLyr")), outlierMultiplier_(conf.getParameter("outlierMultiplier")), minNumLayerCluster_(conf.getParameter("minNumLayerCluster")), eidInputName_(conf.getParameter("eid_input_name")), @@ -64,7 +69,7 @@ void PatternRecognitionbyCLUE3D::dumpTracksters(const std::vector::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "[evt, tracksterId, cells, prob_photon, prob_ele, prob_chad, prob_nhad, layer_i, x_i, y_i, eta_i, phi_i, " - "energy_i, radius_i, rho_i, delta_i, isSeed_i"; + "energy_i, radius_i, rho_i, delta_tr, delta_lyr, isSeed_i"; } int num = 0; @@ -75,14 +80,25 @@ void PatternRecognitionbyCLUE3D::dumpTracksters(const std::vector::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecogntionbyCLUE3D") - << "TracksterInfo: " << eventNumber << sep << num << sep << t.vertices().size() << sep - << t.id_probability(ticl::Trackster::ParticleType::photon) << sep - << t.id_probability(ticl::Trackster::ParticleType::electron) << sep - << t.id_probability(ticl::Trackster::ParticleType::charged_hadron) << sep - << t.id_probability(ticl::Trackster::ParticleType::neutral_hadron) << sep << lyrIdx << sep - << thisLayer.x[soaIdx] << sep << thisLayer.y[soaIdx] << sep << thisLayer.eta[soaIdx] << sep - << thisLayer.phi[soaIdx] << sep << thisLayer.energy[soaIdx] << sep << thisLayer.radius[soaIdx] << sep - << thisLayer.rho[soaIdx] << sep << thisLayer.delta[soaIdx] << sep << thisLayer.isSeed[soaIdx] << '\n'; + << "TracksterInfo: " + << std::setw(4) << eventNumber << sep + << std::setw(4) << num << sep + << std::setw(4) << t.vertices().size() << sep + << std::setw(8) << t.id_probability(ticl::Trackster::ParticleType::photon) << sep + << std::setw(8) << t.id_probability(ticl::Trackster::ParticleType::electron) << sep + << std::setw(8) << t.id_probability(ticl::Trackster::ParticleType::charged_hadron) << sep + << std::setw(8) << t.id_probability(ticl::Trackster::ParticleType::neutral_hadron) << sep + << std::setw(4) << lyrIdx << sep + << std::setw(10) << thisLayer.x[soaIdx] << sep + << std::setw(10) << thisLayer.y[soaIdx] << sep + << std::setw(10) << thisLayer.eta[soaIdx] << sep + << std::setw(10) << thisLayer.phi[soaIdx] << sep + << std::setw(10) << thisLayer.energy[soaIdx] << sep + << std::setw(10) << thisLayer.radius[soaIdx] << sep + << std::setw(10) << thisLayer.rho[soaIdx] << sep + << std::setw(10) << thisLayer.delta[soaIdx].first << sep + << std::setw(10) << thisLayer.delta[soaIdx].second << sep + << std::setw(4) << thisLayer.isSeed[soaIdx]; } } num++; @@ -93,25 +109,37 @@ template void PatternRecognitionbyCLUE3D::dumpClusters(const std::vector> &layerIdx2layerandSoa, const int eventNumber) const { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "[evt, layer, x, y, z, r_over_absz, eta, phi, cells, energy, radius, rho, delta, " - "isSeed, clusterIdx, layerClusterOriginalIdx, SOAidx"; + edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "[evt, layer, isSeed, x, y, z, r_over_absz, eta, phi, cells, energy, energy/rho, rho, delta_tr, delta_lyr, " + " nearestHighLayer, nearestHighSoaIdx, radius, clusterIdx, layerClusterOriginalIdx, SOAidx"; } for (unsigned int layer = 0; layer < clusters_.size(); layer++) { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "On Layer " << layer; - } auto const &thisLayer = clusters_[layer]; int num = 0; for (auto v : thisLayer.x) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecogntionbyCLUE3D") - << "ClusterInfo: " << eventNumber << ", " << layer << ", " << v << ", " << thisLayer.y[num] << ", " - << thisLayer.z[num] << ", " << thisLayer.r_over_absz[num] << ", " - << thisLayer.eta[num] << ", " << thisLayer.phi[num] << ", " << thisLayer.cells[num] << ", " - << thisLayer.energy[num] << ", " << thisLayer.radius[num] << ", " << thisLayer.rho[num] << ", " - << thisLayer.delta[num] << ", " << thisLayer.isSeed[num] << ", " << thisLayer.clusterIndex[num] << ", " - << thisLayer.layerClusterOriginalIdx[num] << ", " << num; + << "ClusterInfo: " << std::setw(8) << eventNumber << ", " + << std::setw(4) << layer << ", " + << std::setw(4) << thisLayer.isSeed[num] << ", " + << std::setw(8) << v << ", " + << std::setw(8) << thisLayer.y[num] << ", " + << std::setw(8) << thisLayer.z[num] << ", " + << std::setw(8) << thisLayer.r_over_absz[num] << ", " + << std::setw(8) << thisLayer.eta[num] << ", " + << std::setw(8) << thisLayer.phi[num] << ", " + << std::setw(4) << thisLayer.cells[num] << ", " + << std::setw(10) << thisLayer.energy[num] << ", " + << std::setw(10) << (thisLayer.energy[num]/thisLayer.rho[num]) << ", " + << std::setw(10) << thisLayer.rho[num] << ", " + << std::setw(10) << thisLayer.delta[num].first << ", " + << std::setw(10) << thisLayer.delta[num].second << ", " + << std::setw(10) << thisLayer.nearestHigher[num].first << ", " + << std::setw(10) << thisLayer.nearestHigher[num].second << ", " + << std::setw(10) << thisLayer.radius[num] << ", " + << std::setw(10) << thisLayer.clusterIndex[num] << ", " + << std::setw(10) << thisLayer.layerClusterOriginalIdx[num] << ", " + << std::setw(4) << num; } ++num; } @@ -202,7 +230,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( clusters_[layer].layerClusterOriginalIdx.emplace_back(layerIdx++); clusters_[layer].nearestHigher.emplace_back(-1, -1); clusters_[layer].rho.emplace_back(0.f); - clusters_[layer].delta.emplace_back(std::numeric_limits::max()); + clusters_[layer].delta.emplace_back(std::make_pair(std::numeric_limits::max(), std::numeric_limits::max())); } for (unsigned int layer = 0; layer < clusters_.size(); layer++) { clusters_[layer].followers.resize(clusters_[layer].x.size()); @@ -436,16 +464,10 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( auto &clustersOnLayer = clusters_[layerId]; unsigned int numberOfClusters = clustersOnLayer.x.size(); - auto isReachable = [&](float x1, float x2, float y1, float y2, float delta_sqr) -> bool { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") - << " is Reachable: " - << ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)) << " vs " << delta_sqr << "[" - << ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) < delta_sqr) << "]\n"; - } - return (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) < delta_sqr; + auto isReachable = [](float r0, float r1, float phi0, float phi1, float delta_sqr) -> bool { + auto delta_phi = reco::deltaPhi(phi0, phi1); + return (r0 - r1) * (r0 - r1) + r1 * r1 * delta_phi * delta_phi < delta_sqr; }; - auto distance_debug = [&](float x1, float x2, float y1, float y2) -> float { return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)); }; @@ -455,7 +477,6 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( auto lastLayerPerSide = static_cast(rhtools_.lastLayer(false)); int minLayer = 0; int maxLayer = 2 * lastLayerPerSide - 1; - const float thisClusterAbsZNormalization = 1./(clustersOnLayer.z[i]*clustersOnLayer.z[i]); if (layerId < lastLayerPerSide) { minLayer = std::max(layerId - densitySiblingLayers_, minLayer); maxLayer = std::min(layerId + densitySiblingLayers_, lastLayerPerSide - 1); @@ -513,47 +534,38 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "OtherEta: " << clustersLayer.eta[layerandSoa.second]; edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "OtherPhi: " << clustersLayer.phi[layerandSoa.second]; } - // extend by 26 mm, roughly 2 cells, more wrt sum of radii - float delta = clustersOnLayer.radius[i] + clustersLayer.radius[layerandSoa.second] + 2.6f; - if (densityOnSameLayer_ && onSameLayer) { - if (isReachable(clustersOnLayer.x[i], - clustersLayer.x[layerandSoa.second], - clustersOnLayer.y[i], - clustersLayer.y[layerandSoa.second], - delta * delta)) { - clustersOnLayer.rho[i] += (clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx ? 1.f : 0.2f) * - clustersLayer.energy[layerandSoa.second]; - } + bool reachable = false; + if (useAbsoluteProjectiveScale_) { + reachable = isReachable(clustersOnLayer.r_over_absz[i]*clustersOnLayer.z[i], + clustersLayer.r_over_absz[layerandSoa.second]*clustersOnLayer.z[i], + clustersOnLayer.phi[i], + clustersLayer.phi[layerandSoa.second], + densityXYDistanceSqr_); } else { - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Distance: " - << reco::deltaR2(clustersOnLayer.eta[i], - clustersOnLayer.phi[i], - clustersLayer.eta[layerandSoa.second], - clustersLayer.phi[layerandSoa.second]); - } - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - auto dist = distance_debug(clustersOnLayer.r_over_absz[i], - clustersLayer.r_over_absz[layerandSoa.second], - clustersOnLayer.r_over_absz[i]*std::abs(clustersOnLayer.phi[i]), - clustersLayer.r_over_absz[layerandSoa.second]*std::abs(clustersLayer.phi[layerandSoa.second])); - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Distance[cm]: " << (dist*clustersOnLayer.z[i]); - } - if (isReachable(clustersOnLayer.r_over_absz[i], - clustersLayer.r_over_absz[layerandSoa.second], - clustersOnLayer.r_over_absz[i]*std::abs(clustersOnLayer.phi[i]), - clustersLayer.r_over_absz[layerandSoa.second]*std::abs(clustersLayer.phi[layerandSoa.second]), thisClusterAbsZNormalization*2.6f*2.6f)) { - /* - if (reco::deltaR2(clustersOnLayer.eta[i], - clustersOnLayer.phi[i], - clustersLayer.eta[layerandSoa.second], - clustersLayer.phi[layerandSoa.second]) < densityEtaPhiDistanceSqr_) {*/ - auto energyToAdd = (clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx ? 1.f : 0.5f) * - clustersLayer.energy[layerandSoa.second]; - clustersOnLayer.rho[i] += energyToAdd; - edm::LogVerbatim("PatternRecogntionbyCLUE3D") - << "Adding " << energyToAdd << " partial " << clustersOnLayer.rho[i]; - } + reachable = (reco::deltaR2(clustersOnLayer.eta[i], + clustersOnLayer.phi[i], + clustersLayer.eta[layerandSoa.second], + clustersLayer.phi[layerandSoa.second]) < densityEtaPhiDistanceSqr_); + } + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Distance[eta,phi]: " + << reco::deltaR2(clustersOnLayer.eta[i], + clustersOnLayer.phi[i], + clustersLayer.eta[layerandSoa.second], + clustersLayer.phi[layerandSoa.second]); + auto dist = distance_debug(clustersOnLayer.r_over_absz[i], + clustersLayer.r_over_absz[layerandSoa.second], + clustersOnLayer.r_over_absz[i]*std::abs(clustersOnLayer.phi[i]), + clustersLayer.r_over_absz[layerandSoa.second]*std::abs(clustersLayer.phi[layerandSoa.second])); + edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Distance[cm]: " << (dist*clustersOnLayer.z[i]); + } + if (reachable) { + float factor_same_layer_different_cluster = (onSameLayer && !densityOnSameLayer_) ? 0.f : 0.5f; + auto energyToAdd = (clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx ? 1.f : 0.5f*factor_same_layer_different_cluster) * + clustersLayer.energy[layerandSoa.second]; + clustersOnLayer.rho[i] += energyToAdd; + edm::LogVerbatim("PatternRecogntionbyCLUE3D") + << "Adding " << energyToAdd << " partial " << clustersOnLayer.rho[i]; } } // end of loop on possible compatible clusters } // end of loop over phi-bin region @@ -573,8 +585,9 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( auto &clustersOnLayer = clusters_[layerId]; unsigned int numberOfClusters = clustersOnLayer.x.size(); - auto distance_debug = [&](float x1, float x2, float y1, float y2) -> float { - return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)); + auto distanceSqr = [](float r0, float r1, float phi0, float phi1) -> float { + auto delta_phi = reco::deltaPhi(phi0, phi1); + return (r0 - r1) * (r0 - r1) + r1 * r1 * delta_phi * delta_phi; }; for (unsigned int i = 0; i < numberOfClusters; i++) { @@ -595,9 +608,10 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( minLayer = std::max(layerId - densitySiblingLayers_, lastLayerPerSide + 1); maxLayer = std::min(layerId + densitySiblingLayers_, maxLayer); } - float maxDelta = std::numeric_limits::max(); + constexpr float maxDelta = std::numeric_limits::max(); float i_delta = maxDelta; std::pair i_nearestHigher(-1, -1); + std::pair nearest_distances(maxDelta, std::numeric_limits::max()); for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) { const auto &tileOnLayer = tiles[currentLayer]; int etaWindow = 3; @@ -620,33 +634,38 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( if ((layerandSoa.first == -1) && (layerandSoa.second == -1)) continue; auto const &clustersOnOtherLayer = clusters_[layerandSoa.first]; - auto dist = distance_debug(clustersOnLayer.r_over_absz[i], - clustersOnOtherLayer.r_over_absz[layerandSoa.second], - clustersOnLayer.r_over_absz[i]*std::abs(clustersOnLayer.phi[i]), - clustersOnOtherLayer.r_over_absz[layerandSoa.second]*std::abs(clustersOnOtherLayer.phi[layerandSoa.second])); - // Rescale the distance in cm on the current layer - dist *=clustersOnLayer.z[i]; - // TODO(rovere): add also z to the distance metric? - /* - float dist = reco::deltaR2(clustersOnLayer.eta[i], - clustersOnLayer.phi[i], - clustersOnOtherLayer.eta[layerandSoa.second], - clustersOnOtherLayer.phi[layerandSoa.second]); - */ + auto dist = maxDelta; + auto dist_transverse = maxDelta; + int dist_layers = std::abs(layerandSoa.first - layerId); + if (useAbsoluteProjectiveScale_) { + dist_transverse = distanceSqr(clustersOnLayer.r_over_absz[i]*clustersOnLayer.z[i], + clustersOnOtherLayer.r_over_absz[layerandSoa.second]*clustersOnLayer.z[i], + clustersOnLayer.phi[i], + clustersOnOtherLayer.phi[layerandSoa.second]); + // Add Z-scale to the final distance + dist = dist_transverse + (clustersOnLayer.z[i] - clustersOnOtherLayer.z[layerandSoa.second])*(clustersOnLayer.z[i] - clustersOnOtherLayer.z[layerandSoa.second]); + } else { + dist = reco::deltaR2(clustersOnLayer.eta[i], + clustersOnLayer.phi[i], + clustersOnOtherLayer.eta[layerandSoa.second], + clustersOnOtherLayer.phi[layerandSoa.second]); + dist_transverse = dist; + } bool foundHigher = (clustersOnOtherLayer.rho[layerandSoa.second] > clustersOnLayer.rho[i]) || - (clustersOnOtherLayer.rho[layerandSoa.second] == clustersOnLayer.rho[i] && - clustersOnOtherLayer.layerClusterOriginalIdx[layerandSoa.second] > - clustersOnLayer.layerClusterOriginalIdx[i]); + (clustersOnOtherLayer.rho[layerandSoa.second] == clustersOnLayer.rho[i] && + clustersOnOtherLayer.layerClusterOriginalIdx[layerandSoa.second] > + clustersOnLayer.layerClusterOriginalIdx[i]); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecogntionbyCLUE3D") - << "Searching nearestHigher on " << currentLayer - << " with rho: " << clustersOnOtherLayer.rho[layerandSoa.second] - << " on layerIdxInSOA: " << layerandSoa.first << ", " << layerandSoa.second - << " with distance: " << dist << " foundHigher: " << foundHigher; + << "Searching nearestHigher on " << currentLayer + << " with rho: " << clustersOnOtherLayer.rho[layerandSoa.second] + << " on layerIdxInSOA: " << layerandSoa.first << ", " << layerandSoa.second + << " with distance: " << sqrt(dist) << " foundHigher: " << foundHigher; } if (foundHigher && dist <= i_delta) { // update i_delta i_delta = dist; + nearest_distances = std::make_pair(sqrt(dist_transverse), dist_layers); // update i_nearestHigher i_nearestHigher = layerandSoa; } @@ -655,19 +674,20 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( } // End of loop on eta bins } // End of loop on layers - bool foundNearestHigherInEtaPhiCylinder = (i_delta != maxDelta); + bool foundNearestInFiducialVolume = (i_delta != maxDelta); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecogntionbyCLUE3D") - << "i_delta: " << i_delta << " passed: " << foundNearestHigherInEtaPhiCylinder << " " - << i_nearestHigher.first << " " << i_nearestHigher.second; + << "i_delta: " << i_delta << " passed: " << foundNearestInFiducialVolume << " " + << i_nearestHigher.first << " " << i_nearestHigher.second + << " distances: " << nearest_distances.first << ", " << nearest_distances.second; } - if (foundNearestHigherInEtaPhiCylinder) { - clustersOnLayer.delta[i] = i_delta; + if (foundNearestInFiducialVolume) { + clustersOnLayer.delta[i] = nearest_distances; clustersOnLayer.nearestHigher[i] = i_nearestHigher; } else { // otherwise delta is guaranteed to be larger outlierDeltaFactor_*delta_c // we can safely maximize delta to be maxDelta - clustersOnLayer.delta[i] = maxDelta; + clustersOnLayer.delta[i] = std::make_pair(maxDelta, std::numeric_limits::max()); clustersOnLayer.nearestHigher[i] = {-1, -1}; } } // End of loop on clusters @@ -679,6 +699,7 @@ int PatternRecognitionbyCLUE3D::findAndAssignTracksters( unsigned int nTracksters = 0; std::vector> localStack; + auto critical_transverse_distance = useAbsoluteProjectiveScale_ ? criticalXYDistance_ : criticalEtaPhiDistance_; // find cluster seeds and outlier for (unsigned int layer = 0; layer < 2 * rhtools_.lastLayer(); layer++) { auto &clustersOnLayer = clusters_[layer]; @@ -687,8 +708,10 @@ int PatternRecognitionbyCLUE3D::findAndAssignTracksters( // initialize clusterIndex clustersOnLayer.clusterIndex[i] = -1; bool isSeed = - (clustersOnLayer.delta[i] > 2.6) && (clustersOnLayer.rho[i] >= criticalDensity_); - bool isOutlier = (clustersOnLayer.delta[i] > outlierMultiplier_ * 2.6) && + (clustersOnLayer.delta[i].first > critical_transverse_distance || clustersOnLayer.delta[i].second > criticalZDistanceLyr_) + && (clustersOnLayer.rho[i] >= criticalDensity_) + && (clustersOnLayer.energy[i]/clustersOnLayer.rho[i] > criticalSelfDensity_); + bool isOutlier = (clustersOnLayer.delta[i].first > outlierMultiplier_ * critical_transverse_distance) && (clustersOnLayer.rho[i] < criticalDensity_); if (isSeed) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { @@ -711,7 +734,7 @@ int PatternRecognitionbyCLUE3D::findAndAssignTracksters( if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Found Outlier on Layer " << layer << " SOAidx: " << i << " with rho: " << clustersOnLayer.rho[i] - << " and delta: " << clustersOnLayer.delta[i]; + << " and delta: " << clustersOnLayer.delta[i].first << ", " << clustersOnLayer.delta[i].second; } } } @@ -738,10 +761,15 @@ template void PatternRecognitionbyCLUE3D::fillPSetDescription(edm::ParameterSetDescription &iDesc) { iDesc.add("algo_verbosity", 0); iDesc.add("criticalDensity", 4)->setComment("in GeV"); - iDesc.add("densitySiblingLayers", 3); + iDesc.add("criticalSelfDensity", 0.15 /* roughly 1/(densitySiblingLayers+1) */)->setComment("Minimum ratio of self_energy/local_density to become a seed."); + iDesc.add("densitySiblingLayers", 5)->setComment("inclusive, layers to consider while computing local density"); iDesc.add("densityEtaPhiDistanceSqr", 0.0008); + iDesc.add("densityXYDistanceSqr", 16 /*6.76*/)->setComment("in cm, 2.6*2.6, distance on the transverse plane to consider for local density"); iDesc.add("densityOnSameLayer", false); - iDesc.add("criticalEtaPhiDistance", 0.035); + iDesc.add("useAbsoluteProjectiveScale", true)->setComment("Express all cuts in terms of r/z*z_0{,phi} projective variables"); + iDesc.add("criticalEtaPhiDistance", 0.035)->setComment("Minimal distance in eta,phi space from nearestHigher to become a seed"); + iDesc.add("criticalXYDistance", 4.0)->setComment("Minimal distance in cm on the XY plane from nearestHigher to become a seed"); + iDesc.add("criticalZDistanceLyr", 3)->setComment("Minimal distance in layers along the Z axis from nearestHigher to become a seed"); iDesc.add("outlierMultiplier", 2); iDesc.add("minNumLayerCluster", 2)->setComment("Not Inclusive"); iDesc.add("eid_input_name", "input"); diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h index da0efefd8380f..436a3116e7cf3 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h @@ -38,7 +38,7 @@ namespace ticl { std::vector energy; std::vector rho; - std::vector delta; + std::vector> delta; std::vector> nearestHigher; std::vector clusterIndex; std::vector layerClusterOriginalIdx; @@ -103,10 +103,15 @@ namespace ticl { edm::ESGetToken caloGeomToken_; const double criticalDensity_; + const double criticalSelfDensity_; const int densitySiblingLayers_; const double densityEtaPhiDistanceSqr_; - const double densityOnSameLayer_; + const double densityXYDistanceSqr_; + const bool densityOnSameLayer_; + const bool useAbsoluteProjectiveScale_; const double criticalEtaPhiDistance_; + const double criticalXYDistance_; + const int criticalZDistanceLyr_; const double outlierMultiplier_; const int minNumLayerCluster_; const std::vector filter_on_categories_; From 8e0677df5ff6891d91e8eb263a86f658caa80f2e Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Tue, 15 Feb 2022 12:19:14 +0100 Subject: [PATCH 03/25] Add ntuplizer ad-hoc label. Extend criticalZDistanceLyr to 5 --- RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index acfa961d2b799..c2017b14ecf72 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -79,8 +79,7 @@ void PatternRecognitionbyCLUE3D::dumpTracksters(const std::vector::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") - << "TracksterInfo: " + edm::LogVerbatim("PatternRecogntionbyCLUE3D_NTP") << std::setw(4) << eventNumber << sep << std::setw(4) << num << sep << std::setw(4) << t.vertices().size() << sep @@ -769,7 +768,7 @@ void PatternRecognitionbyCLUE3D::fillPSetDescription(edm::ParameterSetDes iDesc.add("useAbsoluteProjectiveScale", true)->setComment("Express all cuts in terms of r/z*z_0{,phi} projective variables"); iDesc.add("criticalEtaPhiDistance", 0.035)->setComment("Minimal distance in eta,phi space from nearestHigher to become a seed"); iDesc.add("criticalXYDistance", 4.0)->setComment("Minimal distance in cm on the XY plane from nearestHigher to become a seed"); - iDesc.add("criticalZDistanceLyr", 3)->setComment("Minimal distance in layers along the Z axis from nearestHigher to become a seed"); + iDesc.add("criticalZDistanceLyr", 5)->setComment("Minimal distance in layers along the Z axis from nearestHigher to become a seed"); iDesc.add("outlierMultiplier", 2); iDesc.add("minNumLayerCluster", 2)->setComment("Not Inclusive"); iDesc.add("eid_input_name", "input"); From 669adf6a1023d64c919e1e7e05a8487b63cb5d57 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Tue, 22 Feb 2022 11:01:45 +0100 Subject: [PATCH 04/25] Make kernel factor configurable --- RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc | 8 +++++--- RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h | 1 + 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index c2017b14ecf72..ad8842d23cc1f 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -26,6 +26,7 @@ PatternRecognitionbyCLUE3D::PatternRecognitionbyCLUE3D(const edm::Paramet densitySiblingLayers_(conf.getParameter("densitySiblingLayers")), densityEtaPhiDistanceSqr_(conf.getParameter("densityEtaPhiDistanceSqr")), densityXYDistanceSqr_(conf.getParameter("densityXYDistanceSqr")), + kernelDensityFactor_(conf.getParameter("kernelDensityFactor")), densityOnSameLayer_(conf.getParameter("densityOnSameLayer")), useAbsoluteProjectiveScale_(conf.getParameter("useAbsoluteProjectiveScale")), criticalEtaPhiDistance_(conf.getParameter("criticalEtaPhiDistance")), @@ -559,8 +560,8 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Distance[cm]: " << (dist*clustersOnLayer.z[i]); } if (reachable) { - float factor_same_layer_different_cluster = (onSameLayer && !densityOnSameLayer_) ? 0.f : 0.5f; - auto energyToAdd = (clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx ? 1.f : 0.5f*factor_same_layer_different_cluster) * + float factor_same_layer_different_cluster = (onSameLayer && !densityOnSameLayer_) ? 0.f : 1.f; + auto energyToAdd = (clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx ? 1.f : kernelDensityFactor_*factor_same_layer_different_cluster) * clustersLayer.energy[layerandSoa.second]; clustersOnLayer.rho[i] += energyToAdd; edm::LogVerbatim("PatternRecogntionbyCLUE3D") @@ -761,9 +762,10 @@ void PatternRecognitionbyCLUE3D::fillPSetDescription(edm::ParameterSetDes iDesc.add("algo_verbosity", 0); iDesc.add("criticalDensity", 4)->setComment("in GeV"); iDesc.add("criticalSelfDensity", 0.15 /* roughly 1/(densitySiblingLayers+1) */)->setComment("Minimum ratio of self_energy/local_density to become a seed."); - iDesc.add("densitySiblingLayers", 5)->setComment("inclusive, layers to consider while computing local density"); + iDesc.add("densitySiblingLayers", 5)->setComment("inclusive, layers to consider while computing local density and searching for nearestHigher higher"); iDesc.add("densityEtaPhiDistanceSqr", 0.0008); iDesc.add("densityXYDistanceSqr", 16 /*6.76*/)->setComment("in cm, 2.6*2.6, distance on the transverse plane to consider for local density"); + iDesc.add("kernelDensityFactor", 0.2)->setComment("Kernel factor to be applied to other LC while computing the local density"); iDesc.add("densityOnSameLayer", false); iDesc.add("useAbsoluteProjectiveScale", true)->setComment("Express all cuts in terms of r/z*z_0{,phi} projective variables"); iDesc.add("criticalEtaPhiDistance", 0.035)->setComment("Minimal distance in eta,phi space from nearestHigher to become a seed"); diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h index 436a3116e7cf3..9adb531298eaf 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h @@ -107,6 +107,7 @@ namespace ticl { const int densitySiblingLayers_; const double densityEtaPhiDistanceSqr_; const double densityXYDistanceSqr_; + const double kernelDensityFactor_; const bool densityOnSameLayer_; const bool useAbsoluteProjectiveScale_; const double criticalEtaPhiDistance_; From da782a615c650f218051e97c33277d560bc0479d Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Tue, 22 Feb 2022 11:05:48 +0100 Subject: [PATCH 05/25] Nearest on same layer configurable Add the possibility to avoid searching for the nearest LayerCluster with higher local density on the same layer. This will guide the pattern recognition more in the longitudinal energy-flow direction, rather than in the transverse plane. If anything related on a plane should be gathered together, we assume CLUE will have done that already. If, instead, CLUE did an excellent job separating sibling showers, we do not want to give up on that while running CLUE3D by linking together LC on the same layer. --- RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc | 4 ++++ RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h | 1 + 2 files changed, 5 insertions(+) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index ad8842d23cc1f..646508a62e4d6 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -28,6 +28,7 @@ PatternRecognitionbyCLUE3D::PatternRecognitionbyCLUE3D(const edm::Paramet densityXYDistanceSqr_(conf.getParameter("densityXYDistanceSqr")), kernelDensityFactor_(conf.getParameter("kernelDensityFactor")), densityOnSameLayer_(conf.getParameter("densityOnSameLayer")), + nearestHigherOnSameLayer_(conf.getParameter("nearestHigherOnSameLayer")), useAbsoluteProjectiveScale_(conf.getParameter("useAbsoluteProjectiveScale")), criticalEtaPhiDistance_(conf.getParameter("criticalEtaPhiDistance")), criticalXYDistance_(conf.getParameter("criticalXYDistance")), @@ -613,6 +614,8 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( std::pair i_nearestHigher(-1, -1); std::pair nearest_distances(maxDelta, std::numeric_limits::max()); for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) { + if (! nearestHigherOnSameLayer_ && (layerId == currentLayer)) + continue; const auto &tileOnLayer = tiles[currentLayer]; int etaWindow = 3; int phiWindow = 3; @@ -767,6 +770,7 @@ void PatternRecognitionbyCLUE3D::fillPSetDescription(edm::ParameterSetDes iDesc.add("densityXYDistanceSqr", 16 /*6.76*/)->setComment("in cm, 2.6*2.6, distance on the transverse plane to consider for local density"); iDesc.add("kernelDensityFactor", 0.2)->setComment("Kernel factor to be applied to other LC while computing the local density"); iDesc.add("densityOnSameLayer", false); + iDesc.add("nearestHigherOnSameLayer", false)->setComment("Allow the nearestHigher to be located on the same layer"); iDesc.add("useAbsoluteProjectiveScale", true)->setComment("Express all cuts in terms of r/z*z_0{,phi} projective variables"); iDesc.add("criticalEtaPhiDistance", 0.035)->setComment("Minimal distance in eta,phi space from nearestHigher to become a seed"); iDesc.add("criticalXYDistance", 4.0)->setComment("Minimal distance in cm on the XY plane from nearestHigher to become a seed"); diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h index 9adb531298eaf..de374d4e4c7a7 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h @@ -109,6 +109,7 @@ namespace ticl { const double densityXYDistanceSqr_; const double kernelDensityFactor_; const bool densityOnSameLayer_; + const bool nearestHigherOnSameLayer_; const bool useAbsoluteProjectiveScale_; const double criticalEtaPhiDistance_; const double criticalXYDistance_; From b1b84ef7c7b6c202bf18b3e7de56a0ba9f54073e Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Tue, 22 Feb 2022 11:14:34 +0100 Subject: [PATCH 06/25] Possibility to compute linear local density A boolean flag will control if the local energy density should be rescaled taking into account the effective extension along the Z direction of the area used to compute it. An additional vector is added to the SOA structure to save the position of each layer in HGCAL. This makes the computation of the linear density easier. --- .../plugins/PatternRecognitionbyCLUE3D.cc | 32 +++++++++++++++++-- .../TICL/plugins/PatternRecognitionbyCLUE3D.h | 5 +++ 2 files changed, 34 insertions(+), 3 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 646508a62e4d6..64a4c7d8b1e0c 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -30,6 +30,7 @@ PatternRecognitionbyCLUE3D::PatternRecognitionbyCLUE3D(const edm::Paramet densityOnSameLayer_(conf.getParameter("densityOnSameLayer")), nearestHigherOnSameLayer_(conf.getParameter("nearestHigherOnSameLayer")), useAbsoluteProjectiveScale_(conf.getParameter("useAbsoluteProjectiveScale")), + rescaleDensityByZ_(conf.getParameter("rescaleDensityByZ")), criticalEtaPhiDistance_(conf.getParameter("criticalEtaPhiDistance")), criticalXYDistance_(conf.getParameter("criticalXYDistance")), criticalZDistanceLyr_(conf.getParameter("criticalZDistanceLyr")), @@ -97,6 +98,7 @@ void PatternRecognitionbyCLUE3D::dumpTracksters(const std::vector::dumpClusters(const std::vector::makeTracksters( const int eventNumber = input.ev.eventAuxiliary().event(); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "New Event"; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "New Event"; } edm::EventSetup const &es = input.es; const CaloGeometry &geom = es.getData(caloGeomToken_); rhtools_.setGeometry(geom); + // Assume identical Z-positioning between positive and negative sides. + // Also, layers inside the HGCAL geometry start from 1. + for (unsigned int i = 0; i < rhtools_.lastLayer(); ++i) { + layersPosZ_.push_back(rhtools_.getPositionLayer(i+1).z()); + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Layer " << i << " located at Z: " << layersPosZ_.back(); + } + } + clusters_.clear(); clusters_.resize(2 * rhtools_.lastLayer(false)); std::vector> layerIdx2layerandSoa; //used everywhere also to propagate cluster masking @@ -231,6 +243,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( clusters_[layer].layerClusterOriginalIdx.emplace_back(layerIdx++); clusters_[layer].nearestHigher.emplace_back(-1, -1); clusters_[layer].rho.emplace_back(0.f); + clusters_[layer].z_extension.emplace_back(0.f); clusters_[layer].delta.emplace_back(std::make_pair(std::numeric_limits::max(), std::numeric_limits::max())); } for (unsigned int layer = 0; layer < clusters_.size(); layer++) { @@ -485,6 +498,8 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( minLayer = std::max(layerId - densitySiblingLayers_, lastLayerPerSide); maxLayer = std::min(layerId + densitySiblingLayers_, maxLayer); } + float deltaLayersZ = std::abs(layersPosZ_[maxLayer % lastLayerPerSide] - layersPosZ_[minLayer % lastLayerPerSide]); + for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "RefLayer: " << layerId << " SoaIDX: " << i; @@ -565,13 +580,23 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( auto energyToAdd = (clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx ? 1.f : kernelDensityFactor_*factor_same_layer_different_cluster) * clustersLayer.energy[layerandSoa.second]; clustersOnLayer.rho[i] += energyToAdd; - edm::LogVerbatim("PatternRecogntionbyCLUE3D") - << "Adding " << energyToAdd << " partial " << clustersOnLayer.rho[i]; + clustersOnLayer.z_extension[i] = deltaLayersZ; + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + edm::LogVerbatim("PatternRecognitionbyCLUE3D") + << "Adding " << energyToAdd << " partial " << clustersOnLayer.rho[i]; + } } } // end of loop on possible compatible clusters } // end of loop over phi-bin region } // end of loop over eta-bin region } // end of loop on the sibling layers + if (rescaleDensityByZ_) { + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Rescaling original density: " << clustersOnLayer.rho[i] + << " by Z: " << deltaLayersZ << " to final density/cm: " << clustersOnLayer.rho[i]/deltaLayersZ; + } + clustersOnLayer.rho[i] /= deltaLayersZ; + } } // end of loop over clusters on this layer if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecogntionbyCLUE3D") << std::endl; @@ -772,6 +797,7 @@ void PatternRecognitionbyCLUE3D::fillPSetDescription(edm::ParameterSetDes iDesc.add("densityOnSameLayer", false); iDesc.add("nearestHigherOnSameLayer", false)->setComment("Allow the nearestHigher to be located on the same layer"); iDesc.add("useAbsoluteProjectiveScale", true)->setComment("Express all cuts in terms of r/z*z_0{,phi} projective variables"); + iDesc.add("rescaleDensityByZ", false)->setComment("Rescale local density by the extension of the Z 'volume' explored. The transvere dimension is, at present, fixed and factored out."); iDesc.add("criticalEtaPhiDistance", 0.035)->setComment("Minimal distance in eta,phi space from nearestHigher to become a seed"); iDesc.add("criticalXYDistance", 4.0)->setComment("Minimal distance in cm on the XY plane from nearestHigher to become a seed"); iDesc.add("criticalZDistanceLyr", 5)->setComment("Minimal distance in layers along the Z axis from nearestHigher to become a seed"); diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h index de374d4e4c7a7..073400ad5a83e 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h @@ -37,6 +37,7 @@ namespace ticl { std::vector energy; std::vector rho; + std::vector z_extension; std::vector> delta; std::vector> nearestHigher; @@ -56,6 +57,7 @@ namespace ticl { cells.clear(); energy.clear(); rho.clear(); + z_extension.clear(); delta.clear(); nearestHigher.clear(); clusterIndex.clear(); @@ -75,6 +77,7 @@ namespace ticl { cells.shrink_to_fit(); energy.shrink_to_fit(); rho.shrink_to_fit(); + z_extension.shrink_to_fit(); delta.shrink_to_fit(); nearestHigher.shrink_to_fit(); clusterIndex.shrink_to_fit(); @@ -100,6 +103,7 @@ namespace ticl { void dumpTiles(const TILES&) const; std::vector clusters_; + std::vector layersPosZ_; edm::ESGetToken caloGeomToken_; const double criticalDensity_; @@ -111,6 +115,7 @@ namespace ticl { const bool densityOnSameLayer_; const bool nearestHigherOnSameLayer_; const bool useAbsoluteProjectiveScale_; + const bool rescaleDensityByZ_;; const double criticalEtaPhiDistance_; const double criticalXYDistance_; const int criticalZDistanceLyr_; From 8db27fa41821ba9fe91c792382d8d6de36cba67d Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Tue, 22 Feb 2022 11:17:10 +0100 Subject: [PATCH 07/25] Fix typo in LogVerbatim --- .../plugins/PatternRecognitionbyCLUE3D.cc | 82 +++++++++---------- 1 file changed, 41 insertions(+), 41 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 64a4c7d8b1e0c..7d582a2bf869e 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -56,7 +56,7 @@ void PatternRecognitionbyCLUE3D::dumpTiles(const TILES &tiles) const { int iphi = ((phi % nPhiBin + nPhiBin) % nPhiBin); if (!tiles[layer][offset + iphi].empty()) { if (this->algo_verbosity_ > this->Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Layer: " << layer << " ieta: " << ieta << " phi: " << phi + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Layer: " << layer << " ieta: " << ieta << " phi: " << phi << " " << tiles[layer][offset + iphi].size(); } } @@ -70,9 +70,9 @@ void PatternRecognitionbyCLUE3D::dumpTracksters(const std::vector &tracksters) const { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "[evt, tracksterId, cells, prob_photon, prob_ele, prob_chad, prob_nhad, layer_i, x_i, y_i, eta_i, phi_i, " - "energy_i, radius_i, rho_i, delta_tr, delta_lyr, isSeed_i"; + "energy_i, radius_i, rho_i, z_extension, delta_tr, delta_lyr, isSeed_i"; } int num = 0; @@ -82,7 +82,7 @@ void PatternRecognitionbyCLUE3D::dumpTracksters(const std::vector::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D_NTP") + edm::LogVerbatim("PatternRecognitionbyCLUE3D_NTP") << std::setw(4) << eventNumber << sep << std::setw(4) << num << sep << std::setw(4) << t.vertices().size() << sep @@ -112,7 +112,7 @@ template void PatternRecognitionbyCLUE3D::dumpClusters(const std::vector> &layerIdx2layerandSoa, const int eventNumber) const { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "[evt, layer, isSeed, x, y, z, r_over_absz, eta, phi, cells, energy, energy/rho, rho, delta_tr, delta_lyr, " + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "[evt, layer, isSeed, x, y, z, r_over_absz, eta, phi, cells, energy, energy/rho, rho, z_extension, delta_tr, delta_lyr, " " nearestHighLayer, nearestHighSoaIdx, radius, clusterIdx, layerClusterOriginalIdx, SOAidx"; } @@ -121,7 +121,7 @@ void PatternRecognitionbyCLUE3D::dumpClusters(const std::vector::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "ClusterInfo: " << std::setw(8) << eventNumber << ", " << std::setw(4) << layer << ", " << std::setw(4) << thisLayer.isSeed[num] << ", " @@ -154,7 +154,7 @@ void PatternRecognitionbyCLUE3D::dumpClusters(const std::vector::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "lcIdx: " << lcIdx << " on Layer: " << layerandSoa.first << " SOA: " << layerandSoa.second; } } @@ -196,7 +196,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( for (auto const &lc : input.layerClusters) { if (input.mask[layerIdx] == 0.) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Skipping masked cluster: " << layerIdx; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Skipping masked cluster: " << layerIdx; } layerIdx2layerandSoa.emplace_back(-1, -1); layerIdx++; @@ -262,7 +262,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( auto nTracksters = findAndAssignTracksters(input.tiles, layerIdx2layerandSoa); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Reconstructed " << nTracksters << " tracksters" << std::endl; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Reconstructed " << nTracksters << " tracksters" << std::endl; dumpClusters(layerIdx2layerandSoa, eventNumber); } @@ -272,15 +272,15 @@ void PatternRecognitionbyCLUE3D::makeTracksters( for (unsigned int layer = 0; layer < clusters_.size(); ++layer) { const auto &thisLayer = clusters_[layer]; if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Examining Layer: " << layer; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Examining Layer: " << layer; } for (unsigned int lc = 0; lc < thisLayer.x.size(); ++lc) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Trackster " << thisLayer.clusterIndex[lc]; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Trackster " << thisLayer.clusterIndex[lc]; } if (thisLayer.clusterIndex[lc] >= 0) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << " adding lcIdx: " << thisLayer.layerClusterOriginalIdx[lc]; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << " adding lcIdx: " << thisLayer.layerClusterOriginalIdx[lc]; } result[thisLayer.clusterIndex[lc]].vertices().push_back(thisLayer.layerClusterOriginalIdx[lc]); result[thisLayer.clusterIndex[lc]].vertex_multiplicity().push_back(1); @@ -311,10 +311,10 @@ void PatternRecognitionbyCLUE3D::makeTracksters( energyRegressionAndID(input.layerClusters, input.tfSession, result); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { for (auto const &t : result) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Barycenter: " << t.barycenter(); - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "LCs: " << t.vertices().size(); - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Energy: " << t.raw_energy(); - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Regressed: " << t.regressed_energy(); + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Barycenter: " << t.barycenter(); + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "LCs: " << t.vertices().size(); + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Energy: " << t.raw_energy(); + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Regressed: " << t.regressed_energy(); } } @@ -326,7 +326,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( // Reset internal clusters_ structure of array for next event reset(); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << std::endl; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << std::endl; } } @@ -502,13 +502,13 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "RefLayer: " << layerId << " SoaIDX: " << i; - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "NextLayer: " << currentLayer; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "RefLayer: " << layerId << " SoaIDX: " << i; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "NextLayer: " << currentLayer; } const auto &tileOnLayer = tiles[currentLayer]; bool onSameLayer = (currentLayer == layerId); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "onSameLayer: " << onSameLayer; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "onSameLayer: " << onSameLayer; } const int etaWindow = 2; const int phiWindow = 2; @@ -517,21 +517,21 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( int phiBinMin = tileOnLayer.phiBin(clustersOnLayer.phi[i]) - phiWindow; int phiBinMax = tileOnLayer.phiBin(clustersOnLayer.phi[i]) + phiWindow; if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "eta: " << clustersOnLayer.eta[i]; - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "phi: " << clustersOnLayer.phi[i]; - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "etaBinMin: " << etaBinMin << ", etaBinMax: " << etaBinMax; - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "phiBinMin: " << phiBinMin << ", phiBinMax: " << phiBinMax; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "eta: " << clustersOnLayer.eta[i]; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "phi: " << clustersOnLayer.phi[i]; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "etaBinMin: " << etaBinMin << ", etaBinMax: " << etaBinMax; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "phiBinMin: " << phiBinMin << ", phiBinMax: " << phiBinMax; } for (int ieta = etaBinMin; ieta <= etaBinMax; ++ieta) { auto offset = ieta * nPhiBin; if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "offset: " << offset; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "offset: " << offset; } for (int iphi_it = phiBinMin; iphi_it <= phiBinMax; ++iphi_it) { int iphi = ((iphi_it % nPhiBin + nPhiBin) % nPhiBin); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "iphi: " << iphi; - edm::LogVerbatim("PatternRecogntionbyCLUE3D") + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "iphi: " << iphi; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Entries in tileBin: " << tileOnLayer[offset + iphi].size(); } for (auto otherClusterIdx : tileOnLayer[offset + iphi]) { @@ -539,16 +539,16 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( // Skip masked layer clusters if ((layerandSoa.first == -1) && (layerandSoa.second == -1)) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Skipping masked layerIdx " << otherClusterIdx; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Skipping masked layerIdx " << otherClusterIdx; } continue; } auto const &clustersLayer = clusters_[layerandSoa.first]; if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "OtherLayer: " << layerandSoa.first << " SoaIDX: " << layerandSoa.second; - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "OtherEta: " << clustersLayer.eta[layerandSoa.second]; - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "OtherPhi: " << clustersLayer.phi[layerandSoa.second]; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "OtherEta: " << clustersLayer.eta[layerandSoa.second]; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "OtherPhi: " << clustersLayer.phi[layerandSoa.second]; } bool reachable = false; if (useAbsoluteProjectiveScale_) { @@ -564,7 +564,7 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( clustersLayer.phi[layerandSoa.second]) < densityEtaPhiDistanceSqr_); } if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Distance[eta,phi]: " + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[eta,phi]: " << reco::deltaR2(clustersOnLayer.eta[i], clustersOnLayer.phi[i], clustersLayer.eta[layerandSoa.second], @@ -573,7 +573,7 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( clustersLayer.r_over_absz[layerandSoa.second], clustersOnLayer.r_over_absz[i]*std::abs(clustersOnLayer.phi[i]), clustersLayer.r_over_absz[layerandSoa.second]*std::abs(clustersLayer.phi[layerandSoa.second])); - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "Distance[cm]: " << (dist*clustersOnLayer.z[i]); + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[cm]: " << (dist*clustersOnLayer.z[i]); } if (reachable) { float factor_same_layer_different_cluster = (onSameLayer && !densityOnSameLayer_) ? 0.f : 1.f; @@ -599,7 +599,7 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( } } // end of loop over clusters on this layer if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") << std::endl; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << std::endl; } } @@ -618,7 +618,7 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( for (unsigned int i = 0; i < numberOfClusters; i++) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Starting searching nearestHigher on " << layerId << " with rho: " << clustersOnLayer.rho[i] << " at eta, phi: " << tiles[layerId].etaBin(clustersOnLayer.eta[i]) << ", " << tiles[layerId].etaBin(clustersOnLayer.phi[i]); @@ -653,7 +653,7 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( for (int iphi_it = phiBinMin; iphi_it <= phiBinMax; ++iphi_it) { int iphi = ((iphi_it % nPhiBin + nPhiBin) % nPhiBin); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Searching nearestHigher on " << currentLayer << " eta, phi: " << ieta << ", " << iphi_it; } for (auto otherClusterIdx : tileOnLayer[offset + iphi]) { @@ -684,7 +684,7 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( clustersOnOtherLayer.layerClusterOriginalIdx[layerandSoa.second] > clustersOnLayer.layerClusterOriginalIdx[i]); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Searching nearestHigher on " << currentLayer << " with rho: " << clustersOnOtherLayer.rho[layerandSoa.second] << " on layerIdxInSOA: " << layerandSoa.first << ", " << layerandSoa.second @@ -704,7 +704,7 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( bool foundNearestInFiducialVolume = (i_delta != maxDelta); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "i_delta: " << i_delta << " passed: " << foundNearestInFiducialVolume << " " << i_nearestHigher.first << " " << i_nearestHigher.second << " distances: " << nearest_distances.first << ", " << nearest_distances.second; @@ -743,7 +743,7 @@ int PatternRecognitionbyCLUE3D::findAndAssignTracksters( (clustersOnLayer.rho[i] < criticalDensity_); if (isSeed) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Found seed on Layer " << layer << " SOAidx: " << i << " assigned ClusterIdx: " << nTracksters; } clustersOnLayer.clusterIndex[i] = nTracksters++; @@ -752,7 +752,7 @@ int PatternRecognitionbyCLUE3D::findAndAssignTracksters( } else if (!isOutlier) { auto [lyrIdx, soaIdx] = clustersOnLayer.nearestHigher[i]; if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Found follower on Layer " << layer << " SOAidx: " << i << " attached to cluster on layer: " << lyrIdx << " SOAidx: " << soaIdx; } @@ -760,7 +760,7 @@ int PatternRecognitionbyCLUE3D::findAndAssignTracksters( clusters_[lyrIdx].followers[soaIdx].emplace_back(layer, i); } else { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecogntionbyCLUE3D") + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Found Outlier on Layer " << layer << " SOAidx: " << i << " with rho: " << clustersOnLayer.rho[i] << " and delta: " << clustersOnLayer.delta[i].first << ", " << clustersOnLayer.delta[i].second; } From d004c62e8f08b77c81fdc75d2597ef79a1c6b928 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Tue, 22 Feb 2022 11:17:49 +0100 Subject: [PATCH 08/25] Add energy filter to the SimTrackster producer --- RecoHGCal/TICL/plugins/SimTrackstersProducer.cc | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc b/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc index 2abe1a0f7be40..eb14025345ac5 100644 --- a/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc @@ -57,6 +57,7 @@ class SimTrackstersProducer : public edm::stream::EDProducer<> { const edm::ESGetToken geom_token_; hgcal::RecHitTools rhtools_; const double fractionCut_; + const double energyCut_; }; DEFINE_FWK_MODULE(SimTrackstersProducer); @@ -73,7 +74,8 @@ SimTrackstersProducer::SimTrackstersProducer(const edm::ParameterSet& ps) associatorMapCaloParticleToReco_token_( consumes(ps.getParameter("layerClusterCaloParticleAssociator"))), geom_token_(esConsumes()), - fractionCut_(ps.getParameter("fractionCut")) { + fractionCut_(ps.getParameter("fractionCut")), + energyCut_(ps.getParameter("energyCut")) { produces(); produces>(); produces("fromCPs"); @@ -95,6 +97,7 @@ void SimTrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& des desc.add("layerClusterCaloParticleAssociator", edm::InputTag("layerClusterCaloParticleAssociationProducer")); desc.add("fractionCut", 0.); + desc.add("energyCut", 2.); descriptions.addWithDefaultLabel(desc); } @@ -136,6 +139,9 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) if (cp.g4Tracks()[0].crossedBoundary()) { regr_energy = cp.g4Tracks()[0].getMomentumAtBoundary().energy(); + if (regr_energy < energyCut_) + continue; + addTrackster(cpIndex, lcVec, inputClusterMask, @@ -156,6 +162,9 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) auto const& sc = *(scRef); auto const scIndex = &sc - &simclusters[0]; + if (sc.g4Tracks()[0].getMomentumAtBoundary().energy() < energyCut_) + continue; + addTrackster(scIndex, lcVec, inputClusterMask, From 8506c26da4c4199c0411a591a080966c57edd336 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Tue, 22 Feb 2022 11:18:44 +0100 Subject: [PATCH 09/25] WIP single CLUE iteration --- RecoHGCal/TICL/python/CLUE3DHighStep_cff.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/python/CLUE3DHighStep_cff.py b/RecoHGCal/TICL/python/CLUE3DHighStep_cff.py index 8b72a481af200..e0666aab0ed32 100644 --- a/RecoHGCal/TICL/python/CLUE3DHighStep_cff.py +++ b/RecoHGCal/TICL/python/CLUE3DHighStep_cff.py @@ -22,7 +22,10 @@ itername = "CLUE3DHigh", patternRecognitionBy = "CLUE3D", pluginPatternRecognitionByCLUE3D = dict ( - criticalEtaPhiDistance = 0.025 + criticalDensity = 0.6, + criticalEtaPhiDistance = 0.025, + kernelDensityFactor = 0., + algo_verbosity = 0 ) ) From d8aee61ebed8d2b5f88e4a7e2dfa0366411f7398 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Tue, 1 Mar 2022 08:44:39 +0100 Subject: [PATCH 10/25] Add clusterDimesionXY boolean to CLUE3D The compatibility between layerClusters in order to consider them into the computation of the local-density has been historically driven by a fixed value, either in the eta-phi space or in the r/z*z_ref space. This has the drawback of completely ignoring the physical extension of the LayerCluster for which the local-density is computed. The clusterDimensionXY will, instead, take that into account. If True, the projected r/z coordinate at the plane of the Layerclusters will be compared to an estimate of its physical extension and, if compatible, will be taken into account, otherwise it will be disregarded. This will automatically configure the CLUE3D local-density computation to the properties of each layerClusters and will, de-facto, remove one configuration parameter from the algorithm itself. The estimate of the cluster dimension is rather crude but in the assumption of a "circular" layer cluster it does a good job. For single hit clusters, a specific treatment will be added in the future. --- RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc | 9 +++++++++ RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h | 1 + 2 files changed, 10 insertions(+) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 7d582a2bf869e..8beaeb546ab82 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -30,6 +30,7 @@ PatternRecognitionbyCLUE3D::PatternRecognitionbyCLUE3D(const edm::Paramet densityOnSameLayer_(conf.getParameter("densityOnSameLayer")), nearestHigherOnSameLayer_(conf.getParameter("nearestHigherOnSameLayer")), useAbsoluteProjectiveScale_(conf.getParameter("useAbsoluteProjectiveScale")), + useClusterDimensionXY_(conf.getParameter("useClusterDimensionXY")), rescaleDensityByZ_(conf.getParameter("rescaleDensityByZ")), criticalEtaPhiDistance_(conf.getParameter("criticalEtaPhiDistance")), criticalXYDistance_(conf.getParameter("criticalXYDistance")), @@ -552,11 +553,18 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( } bool reachable = false; if (useAbsoluteProjectiveScale_) { + if (useClusterDimensionXY_) { + auto deltaR_sqr = (clustersOnLayer.r_over_absz[i]*clustersOnLayer.z[i] - clustersLayer.r_over_absz[layerandSoa.second]*clustersOnLayer.z[i]) * + (clustersOnLayer.r_over_absz[i]*clustersOnLayer.z[i] - clustersLayer.r_over_absz[layerandSoa.second]*clustersOnLayer.z[i]); + reachable = deltaR_sqr < 4.*clustersOnLayer.radius[i]*clustersOnLayer.radius[i]; + } + else { reachable = isReachable(clustersOnLayer.r_over_absz[i]*clustersOnLayer.z[i], clustersLayer.r_over_absz[layerandSoa.second]*clustersOnLayer.z[i], clustersOnLayer.phi[i], clustersLayer.phi[layerandSoa.second], densityXYDistanceSqr_); + } } else { reachable = (reco::deltaR2(clustersOnLayer.eta[i], clustersOnLayer.phi[i], @@ -797,6 +805,7 @@ void PatternRecognitionbyCLUE3D::fillPSetDescription(edm::ParameterSetDes iDesc.add("densityOnSameLayer", false); iDesc.add("nearestHigherOnSameLayer", false)->setComment("Allow the nearestHigher to be located on the same layer"); iDesc.add("useAbsoluteProjectiveScale", true)->setComment("Express all cuts in terms of r/z*z_0{,phi} projective variables"); + iDesc.add("useClusterDimensionXY", true)->setComment("Boolean. If true use the estimated cluster radius to determine the cluster compatibility while computing the local density"); iDesc.add("rescaleDensityByZ", false)->setComment("Rescale local density by the extension of the Z 'volume' explored. The transvere dimension is, at present, fixed and factored out."); iDesc.add("criticalEtaPhiDistance", 0.035)->setComment("Minimal distance in eta,phi space from nearestHigher to become a seed"); iDesc.add("criticalXYDistance", 4.0)->setComment("Minimal distance in cm on the XY plane from nearestHigher to become a seed"); diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h index 073400ad5a83e..76f8c51290ea2 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h @@ -115,6 +115,7 @@ namespace ticl { const bool densityOnSameLayer_; const bool nearestHigherOnSameLayer_; const bool useAbsoluteProjectiveScale_; + const bool useClusterDimensionXY_; const bool rescaleDensityByZ_;; const double criticalEtaPhiDistance_; const double criticalXYDistance_; From 998f7a23b02b1ed6fda305a73be1f3e10d0fbdfd Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Tue, 1 Mar 2022 10:46:11 +0100 Subject: [PATCH 11/25] Single Cell treatment in CLUE3D LayerClusters made of a single cell have been specially handled in order to assign a meaningful physical dimension to them. For the Silicon case, the size of the cell as returned by the geometry is used. For the Scintillator case, the extension in eta-phi is retrieved from the geometry and converted to a proper dimension in cm at the correct r and z. For simplicity the Scintillators tiles are assumed to be perfect square, therefore the phi span is used also in the eta direction. --- .../plugins/PatternRecognitionbyCLUE3D.cc | 21 +++++++++++++++++++ RecoHGCal/TICL/python/CLUE3DHighStep_cff.py | 2 +- 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 8beaeb546ab82..8988a01f8cd3c 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -230,6 +230,27 @@ void PatternRecognitionbyCLUE3D::makeTracksters( // below, too. float radius_x = sqrt((sum_sqr_x - (sum_x * sum_x) * invClsize) * invClsize); float radius_y = sqrt((sum_sqr_y - (sum_y * sum_y) * invClsize) * invClsize); + + // The case of single cell layer clusters has to be handled differently. + + if (invClsize == 1.) { + auto detId = lc.hitsAndFractions()[0].first; + // Silicon case + if (rhtools_.isSilicon(detId)) { + radius_x = radius_y = rhtools_.getRadiusToSide(detId); + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Single cell cluster in silicon: " << radius_x << ", " << radius_y; + } + } else { + auto const &point = rhtools_.getPosition(detId); + auto const &eta_phi_window = rhtools_.getScintDEtaDPhi(detId); + radius_x = radius_y = point.perp() * eta_phi_window.second; + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Single cell cluster in scintillator: " << radius_x << ", " << radius_y; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Single cell cluster eta-phi span: " << eta_phi_window.first << ", " << eta_phi_window.second; + } + } + } clusters_[layer].x.emplace_back(lc.x()); clusters_[layer].y.emplace_back(lc.y()); clusters_[layer].z.emplace_back(lc.z()); diff --git a/RecoHGCal/TICL/python/CLUE3DHighStep_cff.py b/RecoHGCal/TICL/python/CLUE3DHighStep_cff.py index e0666aab0ed32..a8a8a0371970f 100644 --- a/RecoHGCal/TICL/python/CLUE3DHighStep_cff.py +++ b/RecoHGCal/TICL/python/CLUE3DHighStep_cff.py @@ -24,7 +24,7 @@ pluginPatternRecognitionByCLUE3D = dict ( criticalDensity = 0.6, criticalEtaPhiDistance = 0.025, - kernelDensityFactor = 0., + kernelDensityFactor = 0.2, algo_verbosity = 0 ) From 8cad68711deb59e024de211df55b893c46c0515a Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Tue, 8 Mar 2022 12:37:14 +0100 Subject: [PATCH 12/25] Code Format --- .../plugins/PatternRecognitionbyCLUE3D.cc | 237 +++++++++--------- .../TICL/plugins/PatternRecognitionbyCLUE3D.h | 3 +- 2 files changed, 125 insertions(+), 115 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 8988a01f8cd3c..9d4dd7e03add2 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -58,7 +58,7 @@ void PatternRecognitionbyCLUE3D::dumpTiles(const TILES &tiles) const { if (!tiles[layer][offset + iphi].empty()) { if (this->algo_verbosity_ > this->Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Layer: " << layer << " ieta: " << ieta << " phi: " << phi - << " " << tiles[layer][offset + iphi].size(); + << " " << tiles[layer][offset + iphi].size(); } } } @@ -84,25 +84,17 @@ void PatternRecognitionbyCLUE3D::dumpTracksters(const std::vector::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D_NTP") - << std::setw(4) << eventNumber << sep - << std::setw(4) << num << sep - << std::setw(4) << t.vertices().size() << sep - << std::setw(8) << t.id_probability(ticl::Trackster::ParticleType::photon) << sep - << std::setw(8) << t.id_probability(ticl::Trackster::ParticleType::electron) << sep - << std::setw(8) << t.id_probability(ticl::Trackster::ParticleType::charged_hadron) << sep - << std::setw(8) << t.id_probability(ticl::Trackster::ParticleType::neutral_hadron) << sep - << std::setw(4) << lyrIdx << sep - << std::setw(10) << thisLayer.x[soaIdx] << sep - << std::setw(10) << thisLayer.y[soaIdx] << sep - << std::setw(10) << thisLayer.eta[soaIdx] << sep - << std::setw(10) << thisLayer.phi[soaIdx] << sep - << std::setw(10) << thisLayer.energy[soaIdx] << sep - << std::setw(10) << thisLayer.radius[soaIdx] << sep - << std::setw(10) << thisLayer.rho[soaIdx] << sep - << std::setw(10) << thisLayer.z_extension[soaIdx] << sep - << std::setw(10) << thisLayer.delta[soaIdx].first << sep - << std::setw(10) << thisLayer.delta[soaIdx].second << sep - << std::setw(4) << thisLayer.isSeed[soaIdx]; + << std::setw(4) << eventNumber << sep << std::setw(4) << num << sep << std::setw(4) << t.vertices().size() + << sep << std::setw(8) << t.id_probability(ticl::Trackster::ParticleType::photon) << sep << std::setw(8) + << t.id_probability(ticl::Trackster::ParticleType::electron) << sep << std::setw(8) + << t.id_probability(ticl::Trackster::ParticleType::charged_hadron) << sep << std::setw(8) + << t.id_probability(ticl::Trackster::ParticleType::neutral_hadron) << sep << std::setw(4) << lyrIdx << sep + << std::setw(10) << thisLayer.x[soaIdx] << sep << std::setw(10) << thisLayer.y[soaIdx] << sep + << std::setw(10) << thisLayer.eta[soaIdx] << sep << std::setw(10) << thisLayer.phi[soaIdx] << sep + << std::setw(10) << thisLayer.energy[soaIdx] << sep << std::setw(10) << thisLayer.radius[soaIdx] << sep + << std::setw(10) << thisLayer.rho[soaIdx] << sep << std::setw(10) << thisLayer.z_extension[soaIdx] << sep + << std::setw(10) << thisLayer.delta[soaIdx].first << sep << std::setw(10) << thisLayer.delta[soaIdx].second + << sep << std::setw(4) << thisLayer.isSeed[soaIdx]; } } num++; @@ -113,8 +105,10 @@ template void PatternRecognitionbyCLUE3D::dumpClusters(const std::vector> &layerIdx2layerandSoa, const int eventNumber) const { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "[evt, layer, isSeed, x, y, z, r_over_absz, eta, phi, cells, energy, energy/rho, rho, z_extension, delta_tr, delta_lyr, " - " nearestHighLayer, nearestHighSoaIdx, radius, clusterIdx, layerClusterOriginalIdx, SOAidx"; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") + << "[evt, layer, isSeed, x, y, z, r_over_absz, eta, phi, cells, energy, energy/rho, rho, z_extension, " + "delta_tr, delta_lyr, " + " nearestHighLayer, nearestHighSoaIdx, radius, clusterIdx, layerClusterOriginalIdx, SOAidx"; } for (unsigned int layer = 0; layer < clusters_.size(); layer++) { @@ -123,28 +117,17 @@ void PatternRecognitionbyCLUE3D::dumpClusters(const std::vector::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "ClusterInfo: " << std::setw(8) << eventNumber << ", " - << std::setw(4) << layer << ", " - << std::setw(4) << thisLayer.isSeed[num] << ", " - << std::setw(8) << v << ", " - << std::setw(8) << thisLayer.y[num] << ", " - << std::setw(8) << thisLayer.z[num] << ", " - << std::setw(8) << thisLayer.r_over_absz[num] << ", " - << std::setw(8) << thisLayer.eta[num] << ", " - << std::setw(8) << thisLayer.phi[num] << ", " - << std::setw(4) << thisLayer.cells[num] << ", " - << std::setw(10) << thisLayer.energy[num] << ", " - << std::setw(10) << (thisLayer.energy[num]/thisLayer.rho[num]) << ", " - << std::setw(10) << thisLayer.rho[num] << ", " - << std::setw(10) << thisLayer.z_extension[num] << ", " - << std::setw(10) << thisLayer.delta[num].first << ", " - << std::setw(10) << thisLayer.delta[num].second << ", " - << std::setw(10) << thisLayer.nearestHigher[num].first << ", " - << std::setw(10) << thisLayer.nearestHigher[num].second << ", " - << std::setw(10) << thisLayer.radius[num] << ", " - << std::setw(10) << thisLayer.clusterIndex[num] << ", " - << std::setw(10) << thisLayer.layerClusterOriginalIdx[num] << ", " - << std::setw(4) << num; + << "ClusterInfo: " << std::setw(8) << eventNumber << ", " << std::setw(4) << layer << ", " << std::setw(4) + << thisLayer.isSeed[num] << ", " << std::setw(8) << v << ", " << std::setw(8) << thisLayer.y[num] << ", " + << std::setw(8) << thisLayer.z[num] << ", " << std::setw(8) << thisLayer.r_over_absz[num] << ", " + << std::setw(8) << thisLayer.eta[num] << ", " << std::setw(8) << thisLayer.phi[num] << ", " << std::setw(4) + << thisLayer.cells[num] << ", " << std::setw(10) << thisLayer.energy[num] << ", " << std::setw(10) + << (thisLayer.energy[num] / thisLayer.rho[num]) << ", " << std::setw(10) << thisLayer.rho[num] << ", " + << std::setw(10) << thisLayer.z_extension[num] << ", " << std::setw(10) << thisLayer.delta[num].first + << ", " << std::setw(10) << thisLayer.delta[num].second << ", " << std::setw(10) + << thisLayer.nearestHigher[num].first << ", " << std::setw(10) << thisLayer.nearestHigher[num].second + << ", " << std::setw(10) << thisLayer.radius[num] << ", " << std::setw(10) << thisLayer.clusterIndex[num] + << ", " << std::setw(10) << thisLayer.layerClusterOriginalIdx[num] << ", " << std::setw(4) << num; } ++num; } @@ -182,10 +165,10 @@ void PatternRecognitionbyCLUE3D::makeTracksters( // Assume identical Z-positioning between positive and negative sides. // Also, layers inside the HGCAL geometry start from 1. for (unsigned int i = 0; i < rhtools_.lastLayer(); ++i) { - layersPosZ_.push_back(rhtools_.getPositionLayer(i+1).z()); - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Layer " << i << " located at Z: " << layersPosZ_.back(); - } + layersPosZ_.push_back(rhtools_.getPositionLayer(i + 1).z()); + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Layer " << i << " located at Z: " << layersPosZ_.back(); + } } clusters_.clear(); @@ -239,22 +222,25 @@ void PatternRecognitionbyCLUE3D::makeTracksters( if (rhtools_.isSilicon(detId)) { radius_x = radius_y = rhtools_.getRadiusToSide(detId); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Single cell cluster in silicon: " << radius_x << ", " << radius_y; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") + << "Single cell cluster in silicon: " << radius_x << ", " << radius_y; } } else { auto const &point = rhtools_.getPosition(detId); auto const &eta_phi_window = rhtools_.getScintDEtaDPhi(detId); radius_x = radius_y = point.perp() * eta_phi_window.second; if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Single cell cluster in scintillator: " << radius_x << ", " << radius_y; - edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Single cell cluster eta-phi span: " << eta_phi_window.first << ", " << eta_phi_window.second; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") + << "Single cell cluster in scintillator: " << radius_x << ", " << radius_y; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") + << "Single cell cluster eta-phi span: " << eta_phi_window.first << ", " << eta_phi_window.second; } } } clusters_[layer].x.emplace_back(lc.x()); clusters_[layer].y.emplace_back(lc.y()); clusters_[layer].z.emplace_back(lc.z()); - clusters_[layer].r_over_absz.emplace_back(sqrt(lc.x()*lc.x()+lc.y()*lc.y())/std::abs(lc.z())); + clusters_[layer].r_over_absz.emplace_back(sqrt(lc.x() * lc.x() + lc.y() * lc.y()) / std::abs(lc.z())); clusters_[layer].radius.emplace_back(radius_x + radius_y); clusters_[layer].eta.emplace_back(lc.eta()); clusters_[layer].phi.emplace_back(lc.phi()); @@ -266,7 +252,8 @@ void PatternRecognitionbyCLUE3D::makeTracksters( clusters_[layer].nearestHigher.emplace_back(-1, -1); clusters_[layer].rho.emplace_back(0.f); clusters_[layer].z_extension.emplace_back(0.f); - clusters_[layer].delta.emplace_back(std::make_pair(std::numeric_limits::max(), std::numeric_limits::max())); + clusters_[layer].delta.emplace_back( + std::make_pair(std::numeric_limits::max(), std::numeric_limits::max())); } for (unsigned int layer = 0; layer < clusters_.size(); layer++) { clusters_[layer].followers.resize(clusters_[layer].x.size()); @@ -575,44 +562,48 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( bool reachable = false; if (useAbsoluteProjectiveScale_) { if (useClusterDimensionXY_) { - auto deltaR_sqr = (clustersOnLayer.r_over_absz[i]*clustersOnLayer.z[i] - clustersLayer.r_over_absz[layerandSoa.second]*clustersOnLayer.z[i]) * - (clustersOnLayer.r_over_absz[i]*clustersOnLayer.z[i] - clustersLayer.r_over_absz[layerandSoa.second]*clustersOnLayer.z[i]); - reachable = deltaR_sqr < 4.*clustersOnLayer.radius[i]*clustersOnLayer.radius[i]; - } - else { - reachable = isReachable(clustersOnLayer.r_over_absz[i]*clustersOnLayer.z[i], - clustersLayer.r_over_absz[layerandSoa.second]*clustersOnLayer.z[i], - clustersOnLayer.phi[i], - clustersLayer.phi[layerandSoa.second], - densityXYDistanceSqr_); + auto deltaR_sqr = (clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i] - + clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i]) * + (clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i] - + clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i]); + reachable = deltaR_sqr < 4. * clustersOnLayer.radius[i] * clustersOnLayer.radius[i]; + } else { + reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i], + clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i], + clustersOnLayer.phi[i], + clustersLayer.phi[layerandSoa.second], + densityXYDistanceSqr_); } } else { reachable = (reco::deltaR2(clustersOnLayer.eta[i], - clustersOnLayer.phi[i], - clustersLayer.eta[layerandSoa.second], - clustersLayer.phi[layerandSoa.second]) < densityEtaPhiDistanceSqr_); + clustersOnLayer.phi[i], + clustersLayer.eta[layerandSoa.second], + clustersLayer.phi[layerandSoa.second]) < densityEtaPhiDistanceSqr_); } if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[eta,phi]: " - << reco::deltaR2(clustersOnLayer.eta[i], - clustersOnLayer.phi[i], - clustersLayer.eta[layerandSoa.second], - clustersLayer.phi[layerandSoa.second]); - auto dist = distance_debug(clustersOnLayer.r_over_absz[i], + << reco::deltaR2(clustersOnLayer.eta[i], + clustersOnLayer.phi[i], + clustersLayer.eta[layerandSoa.second], + clustersLayer.phi[layerandSoa.second]); + auto dist = distance_debug( + clustersOnLayer.r_over_absz[i], clustersLayer.r_over_absz[layerandSoa.second], - clustersOnLayer.r_over_absz[i]*std::abs(clustersOnLayer.phi[i]), - clustersLayer.r_over_absz[layerandSoa.second]*std::abs(clustersLayer.phi[layerandSoa.second])); - edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[cm]: " << (dist*clustersOnLayer.z[i]); + clustersOnLayer.r_over_absz[i] * std::abs(clustersOnLayer.phi[i]), + clustersLayer.r_over_absz[layerandSoa.second] * std::abs(clustersLayer.phi[layerandSoa.second])); + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[cm]: " << (dist * clustersOnLayer.z[i]); } if (reachable) { float factor_same_layer_different_cluster = (onSameLayer && !densityOnSameLayer_) ? 0.f : 1.f; - auto energyToAdd = (clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx ? 1.f : kernelDensityFactor_*factor_same_layer_different_cluster) * - clustersLayer.energy[layerandSoa.second]; + auto energyToAdd = (clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx + ? 1.f + : kernelDensityFactor_ * factor_same_layer_different_cluster) * + clustersLayer.energy[layerandSoa.second]; clustersOnLayer.rho[i] += energyToAdd; clustersOnLayer.z_extension[i] = deltaLayersZ; if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "Adding " << energyToAdd << " partial " << clustersOnLayer.rho[i]; + << "Adding " << energyToAdd << " partial " << clustersOnLayer.rho[i]; } } } // end of loop on possible compatible clusters @@ -621,12 +612,13 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( } // end of loop on the sibling layers if (rescaleDensityByZ_) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Rescaling original density: " << clustersOnLayer.rho[i] - << " by Z: " << deltaLayersZ << " to final density/cm: " << clustersOnLayer.rho[i]/deltaLayersZ; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") + << "Rescaling original density: " << clustersOnLayer.rho[i] << " by Z: " << deltaLayersZ + << " to final density/cm: " << clustersOnLayer.rho[i] / deltaLayersZ; } clustersOnLayer.rho[i] /= deltaLayersZ; } - } // end of loop over clusters on this layer + } // end of loop over clusters on this layer if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << std::endl; } @@ -668,7 +660,7 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( std::pair i_nearestHigher(-1, -1); std::pair nearest_distances(maxDelta, std::numeric_limits::max()); for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) { - if (! nearestHigherOnSameLayer_ && (layerId == currentLayer)) + if (!nearestHigherOnSameLayer_ && (layerId == currentLayer)) continue; const auto &tileOnLayer = tiles[currentLayer]; int etaWindow = 3; @@ -695,29 +687,30 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( auto dist_transverse = maxDelta; int dist_layers = std::abs(layerandSoa.first - layerId); if (useAbsoluteProjectiveScale_) { - dist_transverse = distanceSqr(clustersOnLayer.r_over_absz[i]*clustersOnLayer.z[i], - clustersOnOtherLayer.r_over_absz[layerandSoa.second]*clustersOnLayer.z[i], - clustersOnLayer.phi[i], - clustersOnOtherLayer.phi[layerandSoa.second]); + dist_transverse = distanceSqr(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i], + clustersOnOtherLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i], + clustersOnLayer.phi[i], + clustersOnOtherLayer.phi[layerandSoa.second]); // Add Z-scale to the final distance - dist = dist_transverse + (clustersOnLayer.z[i] - clustersOnOtherLayer.z[layerandSoa.second])*(clustersOnLayer.z[i] - clustersOnOtherLayer.z[layerandSoa.second]); + dist = dist_transverse + (clustersOnLayer.z[i] - clustersOnOtherLayer.z[layerandSoa.second]) * + (clustersOnLayer.z[i] - clustersOnOtherLayer.z[layerandSoa.second]); } else { dist = reco::deltaR2(clustersOnLayer.eta[i], - clustersOnLayer.phi[i], - clustersOnOtherLayer.eta[layerandSoa.second], - clustersOnOtherLayer.phi[layerandSoa.second]); + clustersOnLayer.phi[i], + clustersOnOtherLayer.eta[layerandSoa.second], + clustersOnOtherLayer.phi[layerandSoa.second]); dist_transverse = dist; } bool foundHigher = (clustersOnOtherLayer.rho[layerandSoa.second] > clustersOnLayer.rho[i]) || - (clustersOnOtherLayer.rho[layerandSoa.second] == clustersOnLayer.rho[i] && - clustersOnOtherLayer.layerClusterOriginalIdx[layerandSoa.second] > - clustersOnLayer.layerClusterOriginalIdx[i]); + (clustersOnOtherLayer.rho[layerandSoa.second] == clustersOnLayer.rho[i] && + clustersOnOtherLayer.layerClusterOriginalIdx[layerandSoa.second] > + clustersOnLayer.layerClusterOriginalIdx[i]); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "Searching nearestHigher on " << currentLayer - << " with rho: " << clustersOnOtherLayer.rho[layerandSoa.second] - << " on layerIdxInSOA: " << layerandSoa.first << ", " << layerandSoa.second - << " with distance: " << sqrt(dist) << " foundHigher: " << foundHigher; + << "Searching nearestHigher on " << currentLayer + << " with rho: " << clustersOnOtherLayer.rho[layerandSoa.second] + << " on layerIdxInSOA: " << layerandSoa.first << ", " << layerandSoa.second + << " with distance: " << sqrt(dist) << " foundHigher: " << foundHigher; } if (foundHigher && dist <= i_delta) { // update i_delta @@ -734,9 +727,9 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( bool foundNearestInFiducialVolume = (i_delta != maxDelta); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "i_delta: " << i_delta << " passed: " << foundNearestInFiducialVolume << " " - << i_nearestHigher.first << " " << i_nearestHigher.second - << " distances: " << nearest_distances.first << ", " << nearest_distances.second; + << "i_delta: " << i_delta << " passed: " << foundNearestInFiducialVolume << " " << i_nearestHigher.first + << " " << i_nearestHigher.second << " distances: " << nearest_distances.first << ", " + << nearest_distances.second; } if (foundNearestInFiducialVolume) { clustersOnLayer.delta[i] = nearest_distances; @@ -764,10 +757,10 @@ int PatternRecognitionbyCLUE3D::findAndAssignTracksters( for (unsigned int i = 0; i < numberOfClusters; i++) { // initialize clusterIndex clustersOnLayer.clusterIndex[i] = -1; - bool isSeed = - (clustersOnLayer.delta[i].first > critical_transverse_distance || clustersOnLayer.delta[i].second > criticalZDistanceLyr_) - && (clustersOnLayer.rho[i] >= criticalDensity_) - && (clustersOnLayer.energy[i]/clustersOnLayer.rho[i] > criticalSelfDensity_); + bool isSeed = (clustersOnLayer.delta[i].first > critical_transverse_distance || + clustersOnLayer.delta[i].second > criticalZDistanceLyr_) && + (clustersOnLayer.rho[i] >= criticalDensity_) && + (clustersOnLayer.energy[i] / clustersOnLayer.rho[i] > criticalSelfDensity_); bool isOutlier = (clustersOnLayer.delta[i].first > outlierMultiplier_ * critical_transverse_distance) && (clustersOnLayer.rho[i] < criticalDensity_); if (isSeed) { @@ -818,19 +811,35 @@ template void PatternRecognitionbyCLUE3D::fillPSetDescription(edm::ParameterSetDescription &iDesc) { iDesc.add("algo_verbosity", 0); iDesc.add("criticalDensity", 4)->setComment("in GeV"); - iDesc.add("criticalSelfDensity", 0.15 /* roughly 1/(densitySiblingLayers+1) */)->setComment("Minimum ratio of self_energy/local_density to become a seed."); - iDesc.add("densitySiblingLayers", 5)->setComment("inclusive, layers to consider while computing local density and searching for nearestHigher higher"); + iDesc.add("criticalSelfDensity", 0.15 /* roughly 1/(densitySiblingLayers+1) */) + ->setComment("Minimum ratio of self_energy/local_density to become a seed."); + iDesc.add("densitySiblingLayers", 5) + ->setComment( + "inclusive, layers to consider while computing local density and searching for nearestHigher higher"); iDesc.add("densityEtaPhiDistanceSqr", 0.0008); - iDesc.add("densityXYDistanceSqr", 16 /*6.76*/)->setComment("in cm, 2.6*2.6, distance on the transverse plane to consider for local density"); - iDesc.add("kernelDensityFactor", 0.2)->setComment("Kernel factor to be applied to other LC while computing the local density"); + iDesc.add("densityXYDistanceSqr", 16 /*6.76*/) + ->setComment("in cm, 2.6*2.6, distance on the transverse plane to consider for local density"); + iDesc.add("kernelDensityFactor", 0.2) + ->setComment("Kernel factor to be applied to other LC while computing the local density"); iDesc.add("densityOnSameLayer", false); - iDesc.add("nearestHigherOnSameLayer", false)->setComment("Allow the nearestHigher to be located on the same layer"); - iDesc.add("useAbsoluteProjectiveScale", true)->setComment("Express all cuts in terms of r/z*z_0{,phi} projective variables"); - iDesc.add("useClusterDimensionXY", true)->setComment("Boolean. If true use the estimated cluster radius to determine the cluster compatibility while computing the local density"); - iDesc.add("rescaleDensityByZ", false)->setComment("Rescale local density by the extension of the Z 'volume' explored. The transvere dimension is, at present, fixed and factored out."); - iDesc.add("criticalEtaPhiDistance", 0.035)->setComment("Minimal distance in eta,phi space from nearestHigher to become a seed"); - iDesc.add("criticalXYDistance", 4.0)->setComment("Minimal distance in cm on the XY plane from nearestHigher to become a seed"); - iDesc.add("criticalZDistanceLyr", 5)->setComment("Minimal distance in layers along the Z axis from nearestHigher to become a seed"); + iDesc.add("nearestHigherOnSameLayer", false) + ->setComment("Allow the nearestHigher to be located on the same layer"); + iDesc.add("useAbsoluteProjectiveScale", true) + ->setComment("Express all cuts in terms of r/z*z_0{,phi} projective variables"); + iDesc.add("useClusterDimensionXY", true) + ->setComment( + "Boolean. If true use the estimated cluster radius to determine the cluster compatibility while computing " + "the local density"); + iDesc.add("rescaleDensityByZ", false) + ->setComment( + "Rescale local density by the extension of the Z 'volume' explored. The transvere dimension is, at present, " + "fixed and factored out."); + iDesc.add("criticalEtaPhiDistance", 0.035) + ->setComment("Minimal distance in eta,phi space from nearestHigher to become a seed"); + iDesc.add("criticalXYDistance", 4.0) + ->setComment("Minimal distance in cm on the XY plane from nearestHigher to become a seed"); + iDesc.add("criticalZDistanceLyr", 5) + ->setComment("Minimal distance in layers along the Z axis from nearestHigher to become a seed"); iDesc.add("outlierMultiplier", 2); iDesc.add("minNumLayerCluster", 2)->setComment("Not Inclusive"); iDesc.add("eid_input_name", "input"); diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h index 76f8c51290ea2..87d8ec56f5ce4 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h @@ -116,7 +116,8 @@ namespace ticl { const bool nearestHigherOnSameLayer_; const bool useAbsoluteProjectiveScale_; const bool useClusterDimensionXY_; - const bool rescaleDensityByZ_;; + const bool rescaleDensityByZ_; + ; const double criticalEtaPhiDistance_; const double criticalXYDistance_; const int criticalZDistanceLyr_; From d8965e9eed8f100cbe0f76e05a404cc471a8c329 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Tue, 8 Mar 2022 14:06:13 +0100 Subject: [PATCH 13/25] Revert "Add energy filter to SimTrackster producer" This reverts commit d004c62e8f08b77c81fdc75d2597ef79a1c6b928. For the time being we do not integrate this development in release. Most likely a more flexible mechanism will be added to selectively filter SimTrackster. --- RecoHGCal/TICL/plugins/SimTrackstersProducer.cc | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc b/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc index eb14025345ac5..2abe1a0f7be40 100644 --- a/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/SimTrackstersProducer.cc @@ -57,7 +57,6 @@ class SimTrackstersProducer : public edm::stream::EDProducer<> { const edm::ESGetToken geom_token_; hgcal::RecHitTools rhtools_; const double fractionCut_; - const double energyCut_; }; DEFINE_FWK_MODULE(SimTrackstersProducer); @@ -74,8 +73,7 @@ SimTrackstersProducer::SimTrackstersProducer(const edm::ParameterSet& ps) associatorMapCaloParticleToReco_token_( consumes(ps.getParameter("layerClusterCaloParticleAssociator"))), geom_token_(esConsumes()), - fractionCut_(ps.getParameter("fractionCut")), - energyCut_(ps.getParameter("energyCut")) { + fractionCut_(ps.getParameter("fractionCut")) { produces(); produces>(); produces("fromCPs"); @@ -97,7 +95,6 @@ void SimTrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& des desc.add("layerClusterCaloParticleAssociator", edm::InputTag("layerClusterCaloParticleAssociationProducer")); desc.add("fractionCut", 0.); - desc.add("energyCut", 2.); descriptions.addWithDefaultLabel(desc); } @@ -139,9 +136,6 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) if (cp.g4Tracks()[0].crossedBoundary()) { regr_energy = cp.g4Tracks()[0].getMomentumAtBoundary().energy(); - if (regr_energy < energyCut_) - continue; - addTrackster(cpIndex, lcVec, inputClusterMask, @@ -162,9 +156,6 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) auto const& sc = *(scRef); auto const scIndex = &sc - &simclusters[0]; - if (sc.g4Tracks()[0].getMomentumAtBoundary().energy() < energyCut_) - continue; - addTrackster(scIndex, lcVec, inputClusterMask, From d1b2d6c6fd8d10ca3619e10f40026ab7de85d8ac Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Sat, 19 Mar 2022 11:10:00 +0100 Subject: [PATCH 14/25] Better debug formatting --- .../plugins/PatternRecognitionbyCLUE3D.cc | 76 +++++++++++++------ .../TICL/plugins/PatternRecognitionbyCLUE3D.h | 2 +- 2 files changed, 55 insertions(+), 23 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 9d4dd7e03add2..01fa25aa21e26 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -102,13 +102,14 @@ void PatternRecognitionbyCLUE3D::dumpTracksters(const std::vector -void PatternRecognitionbyCLUE3D::dumpClusters(const std::vector> &layerIdx2layerandSoa, +void PatternRecognitionbyCLUE3D::dumpClusters(const TILES &tiles, + const std::vector> &layerIdx2layerandSoa, const int eventNumber) const { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "[evt, layer, isSeed, x, y, z, r_over_absz, eta, phi, cells, energy, energy/rho, rho, z_extension, " - "delta_tr, delta_lyr, " - " nearestHighLayer, nearestHighSoaIdx, radius, clusterIdx, layerClusterOriginalIdx, SOAidx"; + << "[evt, lyr, Seed, x, y, z, r/|z|, eta, phi, etab, phib, cells, enrgy, e/rho, rho, z_ext, " + "dlt_tr,\t dlt_lyr, " + " nestHL, nestHSoaIdx, radius, clIdx, lClOrigIdx, SOAidx"; } for (unsigned int layer = 0; layer < clusters_.size(); layer++) { @@ -117,17 +118,42 @@ void PatternRecognitionbyCLUE3D::dumpClusters(const std::vector::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "ClusterInfo: " << std::setw(8) << eventNumber << ", " << std::setw(4) << layer << ", " << std::setw(4) - << thisLayer.isSeed[num] << ", " << std::setw(8) << v << ", " << std::setw(8) << thisLayer.y[num] << ", " - << std::setw(8) << thisLayer.z[num] << ", " << std::setw(8) << thisLayer.r_over_absz[num] << ", " - << std::setw(8) << thisLayer.eta[num] << ", " << std::setw(8) << thisLayer.phi[num] << ", " << std::setw(4) - << thisLayer.cells[num] << ", " << std::setw(10) << thisLayer.energy[num] << ", " << std::setw(10) - << (thisLayer.energy[num] / thisLayer.rho[num]) << ", " << std::setw(10) << thisLayer.rho[num] << ", " - << std::setw(10) << thisLayer.z_extension[num] << ", " << std::setw(10) << thisLayer.delta[num].first - << ", " << std::setw(10) << thisLayer.delta[num].second << ", " << std::setw(10) - << thisLayer.nearestHigher[num].first << ", " << std::setw(10) << thisLayer.nearestHigher[num].second - << ", " << std::setw(10) << thisLayer.radius[num] << ", " << std::setw(10) << thisLayer.clusterIndex[num] - << ", " << std::setw(10) << thisLayer.layerClusterOriginalIdx[num] << ", " << std::setw(4) << num; + .format("{:>4d},", eventNumber) + .format("{:>4d},", layer) + .format("{:>5d}, ", thisLayer.isSeed[num]) + .format("{:05.3f}, ", v) + .format("{:05.3f}, ", thisLayer.y[num]) + .format("{:05.3f}, ", thisLayer.z[num]) + .format("{:05.3f}, ", thisLayer.r_over_absz[num]) + .format("{:05.3f}, ", thisLayer.eta[num]) + .format("{:05.3f}, ", thisLayer.phi[num]) + .format("{:>5d}, ", tiles[layer].etaBin(thisLayer.eta[num])) + .format("{:>5d}, ", tiles[layer].phiBin(thisLayer.phi[num])) + .format("{:>4d}, ", thisLayer.cells[num]) + .format("{:05.3f}, ", thisLayer.energy[num]) + .format("{:05.3f}, ", thisLayer.energy[num]/thisLayer.rho[num]) + .format("{:05.3f}, ", thisLayer.rho[num]) + .format("{:05.3f}, ", thisLayer.z_extension[num]) + .format("{:05.3g},", thisLayer.delta[num].first) + .format("\t{:>12d}, ", thisLayer.delta[num].second) + .format("{:>4d}, ", thisLayer.nearestHigher[num].first) + .format("{:>8d}, ", thisLayer.nearestHigher[num].second) + .format("{:05.3f}, ", thisLayer.radius[num]) + .format("{:>4d},", thisLayer.clusterIndex[num]) + .format("{:>4d},", thisLayer.layerClusterOriginalIdx[num]) + .format("{:>4d}, ", num) + .format("ClusterInfo"); +// << "ClusterInfo: " << std::setw(8) << eventNumber << ", " << std::setw(4) << layer << ", " << std::setw(4) +// << thisLayer.isSeed[num] << ", " << std::setw(8) << v << ", " << std::setw(8) << thisLayer.y[num] << ", " +// << std::setw(8) << thisLayer.z[num] << ", " << std::setw(8) << thisLayer.r_over_absz[num] << ", " +// << std::setw(8) << thisLayer.eta[num] << ", " << std::setw(8) << thisLayer.phi[num] << ", " << std::setw(4) +// << thisLayer.cells[num] << ", " << std::setw(10) << thisLayer.energy[num] << ", " << std::setw(10) +// << (thisLayer.energy[num] / thisLayer.rho[num]) << ", " << std::setw(10) << thisLayer.rho[num] << ", " +// << std::setw(10) << thisLayer.z_extension[num] << ", " << std::setw(10) << thisLayer.delta[num].first +// << ", " << std::setw(10) << thisLayer.delta[num].second << ", " << std::setw(10) +// << thisLayer.nearestHigher[num].first << ", " << std::setw(10) << thisLayer.nearestHigher[num].second +// << ", " << std::setw(10) << thisLayer.radius[num] << ", " << std::setw(10) << thisLayer.clusterIndex[num] +// << ", " << std::setw(10) << thisLayer.layerClusterOriginalIdx[num] << ", " << std::setw(4) << num; } ++num; } @@ -213,6 +239,10 @@ void PatternRecognitionbyCLUE3D::makeTracksters( // below, too. float radius_x = sqrt((sum_sqr_x - (sum_x * sum_x) * invClsize) * invClsize); float radius_y = sqrt((sum_sqr_y - (sum_y * sum_y) * invClsize) * invClsize); + if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { + edm::LogVerbatim("PatternRecognitionbyCLUE3D") + .format("Cluster rx: {:5.3f}, ry: {:5.3f}, r: {:5.3f}, cells: {:>4d}", radius_x, radius_y, (radius_x+radius_y), lc.hitsAndFractions().size()); + } // The case of single cell layer clusters has to be handled differently. @@ -223,7 +253,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( radius_x = radius_y = rhtools_.getRadiusToSide(detId); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "Single cell cluster in silicon: " << radius_x << ", " << radius_y; + .format("Single cell cluster in silicon. rx: {:5.3f}, ry: {:5.3f}", radius_x, radius_y); } } else { auto const &point = rhtools_.getPosition(detId); @@ -231,9 +261,8 @@ void PatternRecognitionbyCLUE3D::makeTracksters( radius_x = radius_y = point.perp() * eta_phi_window.second; if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "Single cell cluster in scintillator: " << radius_x << ", " << radius_y; - edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "Single cell cluster eta-phi span: " << eta_phi_window.first << ", " << eta_phi_window.second; + .format("Single cell cluster in scintillator. rx: {:5.3f}, ry: {:5.3f}", radius_x, radius_y) + .format("eta-span: {5.3f}, phi-span: {5.3f}", eta_phi_window.first, eta_phi_window.second); } } } @@ -272,7 +301,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( auto nTracksters = findAndAssignTracksters(input.tiles, layerIdx2layerandSoa); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Reconstructed " << nTracksters << " tracksters" << std::endl; - dumpClusters(layerIdx2layerandSoa, eventNumber); + dumpClusters(input.tiles, layerIdx2layerandSoa, eventNumber); } // Build Trackster @@ -592,6 +621,8 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( clustersOnLayer.r_over_absz[i] * std::abs(clustersOnLayer.phi[i]), clustersLayer.r_over_absz[layerandSoa.second] * std::abs(clustersLayer.phi[layerandSoa.second])); edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[cm]: " << (dist * clustersOnLayer.z[i]); + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Energy Other: " << clustersLayer.energy[layerandSoa.second]; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Cluster radius: " << clustersOnLayer.radius[i]; } if (reachable) { float factor_same_layer_different_cluster = (onSameLayer && !densityOnSameLayer_) ? 0.f : 1.f; @@ -642,7 +673,7 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Starting searching nearestHigher on " << layerId << " with rho: " << clustersOnLayer.rho[i] << " at eta, phi: " << tiles[layerId].etaBin(clustersOnLayer.eta[i]) << ", " - << tiles[layerId].etaBin(clustersOnLayer.phi[i]); + << tiles[layerId].phiBin(clustersOnLayer.phi[i]); } // We need to partition the two sides of the HGCAL detector auto lastLayerPerSide = static_cast(rhtools_.lastLayer(false)); @@ -675,7 +706,8 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( int iphi = ((iphi_it % nPhiBin + nPhiBin) % nPhiBin); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "Searching nearestHigher on " << currentLayer << " eta, phi: " << ieta << ", " << iphi_it; + << "Searching nearestHigher on " << currentLayer << " eta, phi: " << ieta << ", " << iphi_it + << " " << iphi << " " << offset << " " << (offset+iphi); } for (auto otherClusterIdx : tileOnLayer[offset + iphi]) { auto const &layerandSoa = layerIdx2layerandSoa[otherClusterIdx]; diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h index 87d8ec56f5ce4..d84266063d0b4 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h @@ -96,7 +96,7 @@ namespace ticl { void calculateLocalDensity(const TILES&, const int layerId, const std::vector>&); void calculateDistanceToHigher(const TILES&, const int layerId, const std::vector>&); int findAndAssignTracksters(const TILES&, const std::vector>&); - void dumpClusters(const std::vector>& layerIdx2layerandSoa, const int) const; + void dumpClusters(const TILES &tiles, const std::vector>& layerIdx2layerandSoa, const int) const; void dumpTracksters(const std::vector>& layerIdx2layerandSoa, const int, const std::vector&) const; From 90e1069ed2e2f90444be219a52f8fcbf24c31ab6 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Sat, 19 Mar 2022 11:10:54 +0100 Subject: [PATCH 15/25] Fix r_z projective distance calculation --- RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 01fa25aa21e26..a6112a43cdfcc 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -591,11 +591,11 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( bool reachable = false; if (useAbsoluteProjectiveScale_) { if (useClusterDimensionXY_) { - auto deltaR_sqr = (clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i] - - clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i]) * - (clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i] - - clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i]); - reachable = deltaR_sqr < 4. * clustersOnLayer.radius[i] * clustersOnLayer.radius[i]; + reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i], + clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i], + clustersOnLayer.phi[i], + clustersLayer.phi[layerandSoa.second], + clustersOnLayer.radius[i] * clustersOnLayer.radius[i]); } else { reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i], clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i], From 9c9b364b2dd34579c68de497090521913b64f263 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Sat, 19 Mar 2022 11:11:44 +0100 Subject: [PATCH 16/25] Shrink eta-phi search window to find nearest-higher --- RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index a6112a43cdfcc..d26c489bbd664 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -694,8 +694,8 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( if (!nearestHigherOnSameLayer_ && (layerId == currentLayer)) continue; const auto &tileOnLayer = tiles[currentLayer]; - int etaWindow = 3; - int phiWindow = 3; + int etaWindow = 1; + int phiWindow = 1; int etaBinMin = std::max(tileOnLayer.etaBin(clustersOnLayer.eta[i]) - etaWindow, 0); int etaBinMax = std::min(tileOnLayer.etaBin(clustersOnLayer.eta[i]) + etaWindow, nEtaBin); int phiBinMin = tileOnLayer.phiBin(clustersOnLayer.phi[i]) - phiWindow; From 58db226368b76f148d66e90f281d1b75257952a3 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Sat, 19 Mar 2022 11:12:42 +0100 Subject: [PATCH 17/25] Remove Z from distance to higher in projective coordinates --- RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index d26c489bbd664..8adbd47659b11 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -723,9 +723,15 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( clustersOnOtherLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i], clustersOnLayer.phi[i], clustersOnOtherLayer.phi[layerandSoa.second]); + // TODO: does that really make sense? + // The real logic should be: + // Among all the LC that falls within my "reachability" distance + // in the XY plane, take the closest in Z with higher local + // density. // Add Z-scale to the final distance - dist = dist_transverse + (clustersOnLayer.z[i] - clustersOnOtherLayer.z[layerandSoa.second]) * - (clustersOnLayer.z[i] - clustersOnOtherLayer.z[layerandSoa.second]); + dist = dist_transverse;// + (clustersOnLayer.z[i] - +// clustersOnOtherLayer.z[layerandSoa.second]) * +// (clustersOnLayer.z[i] - clustersOnOtherLayer.z[layerandSoa.second]); } else { dist = reco::deltaR2(clustersOnLayer.eta[i], clustersOnLayer.phi[i], From f13732e43ca5219c5fad3ae8982253f6aa95f6b2 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Sat, 19 Mar 2022 11:13:18 +0100 Subject: [PATCH 18/25] Disable cluster size compatibility feature --- RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 8adbd47659b11..d17bc209fad3b 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -864,7 +864,7 @@ void PatternRecognitionbyCLUE3D::fillPSetDescription(edm::ParameterSetDes ->setComment("Allow the nearestHigher to be located on the same layer"); iDesc.add("useAbsoluteProjectiveScale", true) ->setComment("Express all cuts in terms of r/z*z_0{,phi} projective variables"); - iDesc.add("useClusterDimensionXY", true) + iDesc.add("useClusterDimensionXY", false) ->setComment( "Boolean. If true use the estimated cluster radius to determine the cluster compatibility while computing " "the local density"); From d47a1380fc9349607cf00aa98a8a10d5780e829c Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Sat, 19 Mar 2022 11:13:39 +0100 Subject: [PATCH 19/25] Further parameter tuning --- RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index d17bc209fad3b..0dbdebd93d9f4 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -851,11 +851,11 @@ void PatternRecognitionbyCLUE3D::fillPSetDescription(edm::ParameterSetDes iDesc.add("criticalDensity", 4)->setComment("in GeV"); iDesc.add("criticalSelfDensity", 0.15 /* roughly 1/(densitySiblingLayers+1) */) ->setComment("Minimum ratio of self_energy/local_density to become a seed."); - iDesc.add("densitySiblingLayers", 5) + iDesc.add("densitySiblingLayers", 3) ->setComment( "inclusive, layers to consider while computing local density and searching for nearestHigher higher"); iDesc.add("densityEtaPhiDistanceSqr", 0.0008); - iDesc.add("densityXYDistanceSqr", 16 /*6.76*/) + iDesc.add("densityXYDistanceSqr", 3.24 /*6.76*/) ->setComment("in cm, 2.6*2.6, distance on the transverse plane to consider for local density"); iDesc.add("kernelDensityFactor", 0.2) ->setComment("Kernel factor to be applied to other LC while computing the local density"); @@ -872,9 +872,9 @@ void PatternRecognitionbyCLUE3D::fillPSetDescription(edm::ParameterSetDes ->setComment( "Rescale local density by the extension of the Z 'volume' explored. The transvere dimension is, at present, " "fixed and factored out."); - iDesc.add("criticalEtaPhiDistance", 0.035) + iDesc.add("criticalEtaPhiDistance", 0.025) ->setComment("Minimal distance in eta,phi space from nearestHigher to become a seed"); - iDesc.add("criticalXYDistance", 4.0) + iDesc.add("criticalXYDistance", 1.8) ->setComment("Minimal distance in cm on the XY plane from nearestHigher to become a seed"); iDesc.add("criticalZDistanceLyr", 5) ->setComment("Minimal distance in layers along the Z axis from nearestHigher to become a seed"); From f01dfb580db29ce146ea72717687dc601d1023ae Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Sat, 19 Mar 2022 23:41:05 +0100 Subject: [PATCH 20/25] Treat Scintillators properly --- .../plugins/PatternRecognitionbyCLUE3D.cc | 24 ++++++++++++++----- .../TICL/plugins/PatternRecognitionbyCLUE3D.h | 3 +++ 2 files changed, 21 insertions(+), 6 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 0dbdebd93d9f4..f7cbee2a62e5f 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -216,6 +216,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( int layer = rhtools_.getLayerWithOffset(firstHitDetId) - 1 + rhtools_.lastLayer(false) * ((rhtools_.zside(firstHitDetId) + 1) >> 1); assert(layer >= 0); + auto detId = lc.hitsAndFractions()[0].first; layerIdx2layerandSoa.emplace_back(layer, clusters_[layer].x.size()); float sum_x = 0.; @@ -247,7 +248,6 @@ void PatternRecognitionbyCLUE3D::makeTracksters( // The case of single cell layer clusters has to be handled differently. if (invClsize == 1.) { - auto detId = lc.hitsAndFractions()[0].first; // Silicon case if (rhtools_.isSilicon(detId)) { radius_x = radius_y = rhtools_.getRadiusToSide(detId); @@ -274,6 +274,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( clusters_[layer].eta.emplace_back(lc.eta()); clusters_[layer].phi.emplace_back(lc.phi()); clusters_[layer].cells.push_back(lc.hitsAndFractions().size()); + clusters_[layer].isSilicon.push_back(rhtools_.isSilicon(detId)); clusters_[layer].energy.emplace_back(lc.energy()); clusters_[layer].isSeed.push_back(false); clusters_[layer].clusterIndex.emplace_back(-1); @@ -597,11 +598,22 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( clustersLayer.phi[layerandSoa.second], clustersOnLayer.radius[i] * clustersOnLayer.radius[i]); } else { - reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i], - clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i], - clustersOnLayer.phi[i], - clustersLayer.phi[layerandSoa.second], - densityXYDistanceSqr_); + // Still differentiate between silicon and Scintillator. + // Silicon has yet to be studied further. + if (clustersOnLayer.isSilicon[i]) { + reachable = isReachable(clustersOnLayer.r_over_absz[i] * + clustersOnLayer.z[i], + clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i], + clustersOnLayer.phi[i], + clustersLayer.phi[layerandSoa.second], + densityXYDistanceSqr_); + } else { + reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i], + clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i], + clustersOnLayer.phi[i], + clustersLayer.phi[layerandSoa.second], + clustersOnLayer.radius[i] * clustersOnLayer.radius[i]); + } } } else { reachable = (reco::deltaR2(clustersOnLayer.eta[i], diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h index d84266063d0b4..fc2327fb7ae64 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h @@ -34,6 +34,7 @@ namespace ticl { std::vector eta; std::vector phi; std::vector cells; + std::vector isSilicon; std::vector energy; std::vector rho; @@ -55,6 +56,7 @@ namespace ticl { eta.clear(); phi.clear(); cells.clear(); + isSilicon.clear(); energy.clear(); rho.clear(); z_extension.clear(); @@ -75,6 +77,7 @@ namespace ticl { eta.shrink_to_fit(); phi.shrink_to_fit(); cells.shrink_to_fit(); + isSilicon.shrink_to_fit(); energy.shrink_to_fit(); rho.shrink_to_fit(); z_extension.shrink_to_fit(); From 03c61f187ce8aef914a32369a9984d5ad18eaa9a Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Sun, 20 Mar 2022 11:52:13 +0100 Subject: [PATCH 21/25] Remove formtting from LogVerbatim --- .../plugins/PatternRecognitionbyCLUE3D.cc | 64 +++++++------------ 1 file changed, 24 insertions(+), 40 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index f7cbee2a62e5f..e8456272af3d6 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -118,42 +118,18 @@ void PatternRecognitionbyCLUE3D::dumpClusters(const TILES &tiles, for (auto v : thisLayer.x) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - .format("{:>4d},", eventNumber) - .format("{:>4d},", layer) - .format("{:>5d}, ", thisLayer.isSeed[num]) - .format("{:05.3f}, ", v) - .format("{:05.3f}, ", thisLayer.y[num]) - .format("{:05.3f}, ", thisLayer.z[num]) - .format("{:05.3f}, ", thisLayer.r_over_absz[num]) - .format("{:05.3f}, ", thisLayer.eta[num]) - .format("{:05.3f}, ", thisLayer.phi[num]) - .format("{:>5d}, ", tiles[layer].etaBin(thisLayer.eta[num])) - .format("{:>5d}, ", tiles[layer].phiBin(thisLayer.phi[num])) - .format("{:>4d}, ", thisLayer.cells[num]) - .format("{:05.3f}, ", thisLayer.energy[num]) - .format("{:05.3f}, ", thisLayer.energy[num]/thisLayer.rho[num]) - .format("{:05.3f}, ", thisLayer.rho[num]) - .format("{:05.3f}, ", thisLayer.z_extension[num]) - .format("{:05.3g},", thisLayer.delta[num].first) - .format("\t{:>12d}, ", thisLayer.delta[num].second) - .format("{:>4d}, ", thisLayer.nearestHigher[num].first) - .format("{:>8d}, ", thisLayer.nearestHigher[num].second) - .format("{:05.3f}, ", thisLayer.radius[num]) - .format("{:>4d},", thisLayer.clusterIndex[num]) - .format("{:>4d},", thisLayer.layerClusterOriginalIdx[num]) - .format("{:>4d}, ", num) - .format("ClusterInfo"); -// << "ClusterInfo: " << std::setw(8) << eventNumber << ", " << std::setw(4) << layer << ", " << std::setw(4) -// << thisLayer.isSeed[num] << ", " << std::setw(8) << v << ", " << std::setw(8) << thisLayer.y[num] << ", " -// << std::setw(8) << thisLayer.z[num] << ", " << std::setw(8) << thisLayer.r_over_absz[num] << ", " -// << std::setw(8) << thisLayer.eta[num] << ", " << std::setw(8) << thisLayer.phi[num] << ", " << std::setw(4) -// << thisLayer.cells[num] << ", " << std::setw(10) << thisLayer.energy[num] << ", " << std::setw(10) -// << (thisLayer.energy[num] / thisLayer.rho[num]) << ", " << std::setw(10) << thisLayer.rho[num] << ", " -// << std::setw(10) << thisLayer.z_extension[num] << ", " << std::setw(10) << thisLayer.delta[num].first -// << ", " << std::setw(10) << thisLayer.delta[num].second << ", " << std::setw(10) -// << thisLayer.nearestHigher[num].first << ", " << std::setw(10) << thisLayer.nearestHigher[num].second -// << ", " << std::setw(10) << thisLayer.radius[num] << ", " << std::setw(10) << thisLayer.clusterIndex[num] -// << ", " << std::setw(10) << thisLayer.layerClusterOriginalIdx[num] << ", " << std::setw(4) << num; + << std::setw(4) << eventNumber << ", " << std::setw(4) << layer << ", " << std::setw(5) + << thisLayer.isSeed[num] << ", " << std::setw(5) << v << ", " << std::setw(5) << thisLayer.y[num] << ", " + << std::setw(5) << thisLayer.z[num] << ", " << std::setw(5) << thisLayer.r_over_absz[num] << ", " + << std::setw(5) << thisLayer.eta[num] << ", " << std::setw(5) << thisLayer.phi[num] << ", " + << std::setw(5) << tiles[layer].etaBin(thisLayer.eta[num]) << ", " << std::setw(5) << tiles[layer].phiBin(thisLayer.phi[num]) << ", " + << std::setw(4) << thisLayer.cells[num] << ", " << std::setw(5) << thisLayer.energy[num] << ", " << std::setw(5) + << (thisLayer.energy[num] / thisLayer.rho[num]) << ", " << std::setw(5) << thisLayer.rho[num] << ", " + << std::setw(5) << thisLayer.z_extension[num] << ", " << std::setw(5) << thisLayer.delta[num].first + << ", " << std::setw(10) << thisLayer.delta[num].second << ", " << std::setw(4) + << thisLayer.nearestHigher[num].first << ", " << std::setw(10) << thisLayer.nearestHigher[num].second + << ", " << std::setw(5) << thisLayer.radius[num] << ", " << std::setw(5) << thisLayer.clusterIndex[num] + << ", " << std::setw(4) << thisLayer.layerClusterOriginalIdx[num] << ", " << std::setw(4) << num << ", ClusterInfo"; } ++num; } @@ -242,7 +218,10 @@ void PatternRecognitionbyCLUE3D::makeTracksters( float radius_y = sqrt((sum_sqr_y - (sum_y * sum_y) * invClsize) * invClsize); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - .format("Cluster rx: {:5.3f}, ry: {:5.3f}, r: {:5.3f}, cells: {:>4d}", radius_x, radius_y, (radius_x+radius_y), lc.hitsAndFractions().size()); + << "cluster rx: " << std::setw(5) << radius_x + << ", ry: " << std::setw(5) << radius_y + << ", r: " << std::setw(5) << (radius_x+radius_y) + << ", cells: " << std::setw(4) << lc.hitsAndFractions().size(); } // The case of single cell layer clusters has to be handled differently. @@ -253,7 +232,9 @@ void PatternRecognitionbyCLUE3D::makeTracksters( radius_x = radius_y = rhtools_.getRadiusToSide(detId); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - .format("Single cell cluster in silicon. rx: {:5.3f}, ry: {:5.3f}", radius_x, radius_y); + << "Single cell cluster in silicon, rx: " + << std::setw(5) << radius_x << ", ry: " + << std::setw(5) << radius_y; } } else { auto const &point = rhtools_.getPosition(detId); @@ -261,8 +242,11 @@ void PatternRecognitionbyCLUE3D::makeTracksters( radius_x = radius_y = point.perp() * eta_phi_window.second; if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - .format("Single cell cluster in scintillator. rx: {:5.3f}, ry: {:5.3f}", radius_x, radius_y) - .format("eta-span: {5.3f}, phi-span: {5.3f}", eta_phi_window.first, eta_phi_window.second); + << "Single cell cluster in scintillator. rx: " + << std::setw(5) << radius_x << ", ry: " + << std::setw(5) << radius_y << ", eta-span: " + << std::setw(5) << eta_phi_window.first << ", phi-span: " + << std::setw(5) << eta_phi_window.second; } } } From a04e3c005ddd0a10e3a99cbbc23abc49f9f0353e Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Sun, 20 Mar 2022 15:31:14 +0100 Subject: [PATCH 22/25] Improve debug formatting --- .../plugins/PatternRecognitionbyCLUE3D.cc | 44 ++++++++++++------- 1 file changed, 27 insertions(+), 17 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index e8456272af3d6..59c87c6861af5 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -107,8 +107,8 @@ void PatternRecognitionbyCLUE3D::dumpClusters(const TILES &tiles, const int eventNumber) const { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "[evt, lyr, Seed, x, y, z, r/|z|, eta, phi, etab, phib, cells, enrgy, e/rho, rho, z_ext, " - "dlt_tr,\t dlt_lyr, " + << "[evt, lyr, Seed, x, y, z, r/|z|, eta, phi, etab, phib, cells, enrgy, e/rho, rho, z_ext, " + " dlt_tr, dlt_lyr, " " nestHL, nestHSoaIdx, radius, clIdx, lClOrigIdx, SOAidx"; } @@ -118,18 +118,31 @@ void PatternRecognitionbyCLUE3D::dumpClusters(const TILES &tiles, for (auto v : thisLayer.x) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << std::setw(4) << eventNumber << ", " << std::setw(4) << layer << ", " << std::setw(5) - << thisLayer.isSeed[num] << ", " << std::setw(5) << v << ", " << std::setw(5) << thisLayer.y[num] << ", " - << std::setw(5) << thisLayer.z[num] << ", " << std::setw(5) << thisLayer.r_over_absz[num] << ", " - << std::setw(5) << thisLayer.eta[num] << ", " << std::setw(5) << thisLayer.phi[num] << ", " - << std::setw(5) << tiles[layer].etaBin(thisLayer.eta[num]) << ", " << std::setw(5) << tiles[layer].phiBin(thisLayer.phi[num]) << ", " - << std::setw(4) << thisLayer.cells[num] << ", " << std::setw(5) << thisLayer.energy[num] << ", " << std::setw(5) - << (thisLayer.energy[num] / thisLayer.rho[num]) << ", " << std::setw(5) << thisLayer.rho[num] << ", " - << std::setw(5) << thisLayer.z_extension[num] << ", " << std::setw(5) << thisLayer.delta[num].first - << ", " << std::setw(10) << thisLayer.delta[num].second << ", " << std::setw(4) - << thisLayer.nearestHigher[num].first << ", " << std::setw(10) << thisLayer.nearestHigher[num].second - << ", " << std::setw(5) << thisLayer.radius[num] << ", " << std::setw(5) << thisLayer.clusterIndex[num] - << ", " << std::setw(4) << thisLayer.layerClusterOriginalIdx[num] << ", " << std::setw(4) << num << ", ClusterInfo"; + << std::setw(4) << eventNumber << ", " + << std::setw(3) << layer << ", " + << std::setw(4) << thisLayer.isSeed[num] << ", " + << std::setprecision(3) << std::fixed << v << ", " + << thisLayer.y[num] << ", " + << thisLayer.z[num] << ", " + << thisLayer.r_over_absz[num] << ", " + << thisLayer.eta[num] << ", " + << thisLayer.phi[num] << ", " + << std::setw(5) << tiles[layer].etaBin(thisLayer.eta[num]) << ", " + << std::setw(5) << tiles[layer].phiBin(thisLayer.phi[num]) << ", " + << std::setw(4) << thisLayer.cells[num] << ", " + << std::setprecision(3) << thisLayer.energy[num] << ", " + << (thisLayer.energy[num] / thisLayer.rho[num]) << ", " + << thisLayer.rho[num] << ", " + << thisLayer.z_extension[num] << ", " + << std::scientific << thisLayer.delta[num].first << ", " + << std::setw(10) << thisLayer.delta[num].second << ", " + << std::setw(5) << thisLayer.nearestHigher[num].first << ", " + << std::setw(10) << thisLayer.nearestHigher[num].second << ", " + << std::defaultfloat << std::setprecision(3) << thisLayer.radius[num] << ", " + << std::setw(5) << thisLayer.clusterIndex[num] << ", " + << std::setw(4) << thisLayer.layerClusterOriginalIdx[num] << ", " + << std::setw(4) << num + << ", ClusterInfo"; } ++num; } @@ -646,9 +659,6 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( clustersOnLayer.rho[i] /= deltaLayersZ; } } // end of loop over clusters on this layer - if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecognitionbyCLUE3D") << std::endl; - } } template From ea85d25e34affc34a62fa83552a958008904436f Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Mon, 21 Mar 2022 10:43:14 +0100 Subject: [PATCH 23/25] Use a different seeding logic for scintillator --- RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 59c87c6861af5..24e25049006c9 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -805,6 +805,12 @@ int PatternRecognitionbyCLUE3D::findAndAssignTracksters( clustersOnLayer.delta[i].second > criticalZDistanceLyr_) && (clustersOnLayer.rho[i] >= criticalDensity_) && (clustersOnLayer.energy[i] / clustersOnLayer.rho[i] > criticalSelfDensity_); + if (!clustersOnLayer.isSilicon[i]) { + isSeed = (clustersOnLayer.delta[i].first > clustersOnLayer.radius[i] || + clustersOnLayer.delta[i].second > criticalZDistanceLyr_) && + (clustersOnLayer.rho[i] >= criticalDensity_) && + (clustersOnLayer.energy[i] / clustersOnLayer.rho[i] > criticalSelfDensity_); + } bool isOutlier = (clustersOnLayer.delta[i].first > outlierMultiplier_ * critical_transverse_distance) && (clustersOnLayer.rho[i] < criticalDensity_); if (isSeed) { From 5e91c01edc597eae38f13ccb40573f5f6206345c Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Mon, 21 Mar 2022 14:46:47 +0100 Subject: [PATCH 24/25] Code format --- .../plugins/PatternRecognitionbyCLUE3D.cc | 105 ++++++++---------- .../TICL/plugins/PatternRecognitionbyCLUE3D.h | 4 +- 2 files changed, 47 insertions(+), 62 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 24e25049006c9..d99f5b4bcfbf3 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -106,10 +106,10 @@ void PatternRecognitionbyCLUE3D::dumpClusters(const TILES &tiles, const std::vector> &layerIdx2layerandSoa, const int eventNumber) const { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "[evt, lyr, Seed, x, y, z, r/|z|, eta, phi, etab, phib, cells, enrgy, e/rho, rho, z_ext, " - " dlt_tr, dlt_lyr, " - " nestHL, nestHSoaIdx, radius, clIdx, lClOrigIdx, SOAidx"; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "[evt, lyr, Seed, x, y, z, r/|z|, eta, phi, " + "etab, phib, cells, enrgy, e/rho, rho, z_ext, " + " dlt_tr, dlt_lyr, " + " nestHL, nestHSoaIdx, radius, clIdx, lClOrigIdx, SOAidx"; } for (unsigned int layer = 0; layer < clusters_.size(); layer++) { @@ -118,31 +118,19 @@ void PatternRecognitionbyCLUE3D::dumpClusters(const TILES &tiles, for (auto v : thisLayer.x) { if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << std::setw(4) << eventNumber << ", " - << std::setw(3) << layer << ", " - << std::setw(4) << thisLayer.isSeed[num] << ", " - << std::setprecision(3) << std::fixed << v << ", " - << thisLayer.y[num] << ", " - << thisLayer.z[num] << ", " - << thisLayer.r_over_absz[num] << ", " - << thisLayer.eta[num] << ", " - << thisLayer.phi[num] << ", " - << std::setw(5) << tiles[layer].etaBin(thisLayer.eta[num]) << ", " - << std::setw(5) << tiles[layer].phiBin(thisLayer.phi[num]) << ", " - << std::setw(4) << thisLayer.cells[num] << ", " - << std::setprecision(3) << thisLayer.energy[num] << ", " - << (thisLayer.energy[num] / thisLayer.rho[num]) << ", " - << thisLayer.rho[num] << ", " - << thisLayer.z_extension[num] << ", " - << std::scientific << thisLayer.delta[num].first << ", " - << std::setw(10) << thisLayer.delta[num].second << ", " - << std::setw(5) << thisLayer.nearestHigher[num].first << ", " - << std::setw(10) << thisLayer.nearestHigher[num].second << ", " - << std::defaultfloat << std::setprecision(3) << thisLayer.radius[num] << ", " - << std::setw(5) << thisLayer.clusterIndex[num] << ", " - << std::setw(4) << thisLayer.layerClusterOriginalIdx[num] << ", " - << std::setw(4) << num - << ", ClusterInfo"; + << std::setw(4) << eventNumber << ", " << std::setw(3) << layer << ", " << std::setw(4) + << thisLayer.isSeed[num] << ", " << std::setprecision(3) << std::fixed << v << ", " << thisLayer.y[num] + << ", " << thisLayer.z[num] << ", " << thisLayer.r_over_absz[num] << ", " << thisLayer.eta[num] << ", " + << thisLayer.phi[num] << ", " << std::setw(5) << tiles[layer].etaBin(thisLayer.eta[num]) << ", " + << std::setw(5) << tiles[layer].phiBin(thisLayer.phi[num]) << ", " << std::setw(4) << thisLayer.cells[num] + << ", " << std::setprecision(3) << thisLayer.energy[num] << ", " + << (thisLayer.energy[num] / thisLayer.rho[num]) << ", " << thisLayer.rho[num] << ", " + << thisLayer.z_extension[num] << ", " << std::scientific << thisLayer.delta[num].first << ", " + << std::setw(10) << thisLayer.delta[num].second << ", " << std::setw(5) + << thisLayer.nearestHigher[num].first << ", " << std::setw(10) << thisLayer.nearestHigher[num].second + << ", " << std::defaultfloat << std::setprecision(3) << thisLayer.radius[num] << ", " << std::setw(5) + << thisLayer.clusterIndex[num] << ", " << std::setw(4) << thisLayer.layerClusterOriginalIdx[num] << ", " + << std::setw(4) << num << ", ClusterInfo"; } ++num; } @@ -231,10 +219,9 @@ void PatternRecognitionbyCLUE3D::makeTracksters( float radius_y = sqrt((sum_sqr_y - (sum_y * sum_y) * invClsize) * invClsize); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "cluster rx: " << std::setw(5) << radius_x - << ", ry: " << std::setw(5) << radius_y - << ", r: " << std::setw(5) << (radius_x+radius_y) - << ", cells: " << std::setw(4) << lc.hitsAndFractions().size(); + << "cluster rx: " << std::setw(5) << radius_x << ", ry: " << std::setw(5) << radius_y + << ", r: " << std::setw(5) << (radius_x + radius_y) << ", cells: " << std::setw(4) + << lc.hitsAndFractions().size(); } // The case of single cell layer clusters has to be handled differently. @@ -244,10 +231,8 @@ void PatternRecognitionbyCLUE3D::makeTracksters( if (rhtools_.isSilicon(detId)) { radius_x = radius_y = rhtools_.getRadiusToSide(detId); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { - edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "Single cell cluster in silicon, rx: " - << std::setw(5) << radius_x << ", ry: " - << std::setw(5) << radius_y; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Single cell cluster in silicon, rx: " << std::setw(5) + << radius_x << ", ry: " << std::setw(5) << radius_y; } } else { auto const &point = rhtools_.getPosition(detId); @@ -255,11 +240,9 @@ void PatternRecognitionbyCLUE3D::makeTracksters( radius_x = radius_y = point.perp() * eta_phi_window.second; if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "Single cell cluster in scintillator. rx: " - << std::setw(5) << radius_x << ", ry: " - << std::setw(5) << radius_y << ", eta-span: " - << std::setw(5) << eta_phi_window.first << ", phi-span: " - << std::setw(5) << eta_phi_window.second; + << "Single cell cluster in scintillator. rx: " << std::setw(5) << radius_x << ", ry: " << std::setw(5) + << radius_y << ", eta-span: " << std::setw(5) << eta_phi_window.first << ", phi-span: " << std::setw(5) + << eta_phi_window.second; } } } @@ -598,18 +581,17 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( // Still differentiate between silicon and Scintillator. // Silicon has yet to be studied further. if (clustersOnLayer.isSilicon[i]) { - reachable = isReachable(clustersOnLayer.r_over_absz[i] * - clustersOnLayer.z[i], - clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i], - clustersOnLayer.phi[i], - clustersLayer.phi[layerandSoa.second], - densityXYDistanceSqr_); + reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i], + clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i], + clustersOnLayer.phi[i], + clustersLayer.phi[layerandSoa.second], + densityXYDistanceSqr_); } else { reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i], - clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i], - clustersOnLayer.phi[i], - clustersLayer.phi[layerandSoa.second], - clustersOnLayer.radius[i] * clustersOnLayer.radius[i]); + clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i], + clustersOnLayer.phi[i], + clustersLayer.phi[layerandSoa.second], + clustersOnLayer.radius[i] * clustersOnLayer.radius[i]); } } } else { @@ -630,7 +612,8 @@ void PatternRecognitionbyCLUE3D::calculateLocalDensity( clustersOnLayer.r_over_absz[i] * std::abs(clustersOnLayer.phi[i]), clustersLayer.r_over_absz[layerandSoa.second] * std::abs(clustersLayer.phi[layerandSoa.second])); edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[cm]: " << (dist * clustersOnLayer.z[i]); - edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Energy Other: " << clustersLayer.energy[layerandSoa.second]; + edm::LogVerbatim("PatternRecognitionbyCLUE3D") + << "Energy Other: " << clustersLayer.energy[layerandSoa.second]; edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Cluster radius: " << clustersOnLayer.radius[i]; } if (reachable) { @@ -712,8 +695,8 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( int iphi = ((iphi_it % nPhiBin + nPhiBin) % nPhiBin); if (PatternRecognitionAlgoBaseT::algo_verbosity_ > PatternRecognitionAlgoBaseT::Advanced) { edm::LogVerbatim("PatternRecognitionbyCLUE3D") - << "Searching nearestHigher on " << currentLayer << " eta, phi: " << ieta << ", " << iphi_it - << " " << iphi << " " << offset << " " << (offset+iphi); + << "Searching nearestHigher on " << currentLayer << " eta, phi: " << ieta << ", " << iphi_it << " " + << iphi << " " << offset << " " << (offset + iphi); } for (auto otherClusterIdx : tileOnLayer[offset + iphi]) { auto const &layerandSoa = layerIdx2layerandSoa[otherClusterIdx]; @@ -735,9 +718,9 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( // in the XY plane, take the closest in Z with higher local // density. // Add Z-scale to the final distance - dist = dist_transverse;// + (clustersOnLayer.z[i] - -// clustersOnOtherLayer.z[layerandSoa.second]) * -// (clustersOnLayer.z[i] - clustersOnOtherLayer.z[layerandSoa.second]); + dist = dist_transverse; // + (clustersOnLayer.z[i] - + // clustersOnOtherLayer.z[layerandSoa.second]) * + // (clustersOnLayer.z[i] - clustersOnOtherLayer.z[layerandSoa.second]); } else { dist = reco::deltaR2(clustersOnLayer.eta[i], clustersOnLayer.phi[i], @@ -807,9 +790,9 @@ int PatternRecognitionbyCLUE3D::findAndAssignTracksters( (clustersOnLayer.energy[i] / clustersOnLayer.rho[i] > criticalSelfDensity_); if (!clustersOnLayer.isSilicon[i]) { isSeed = (clustersOnLayer.delta[i].first > clustersOnLayer.radius[i] || - clustersOnLayer.delta[i].second > criticalZDistanceLyr_) && - (clustersOnLayer.rho[i] >= criticalDensity_) && - (clustersOnLayer.energy[i] / clustersOnLayer.rho[i] > criticalSelfDensity_); + clustersOnLayer.delta[i].second > criticalZDistanceLyr_) && + (clustersOnLayer.rho[i] >= criticalDensity_) && + (clustersOnLayer.energy[i] / clustersOnLayer.rho[i] > criticalSelfDensity_); } bool isOutlier = (clustersOnLayer.delta[i].first > outlierMultiplier_ * critical_transverse_distance) && (clustersOnLayer.rho[i] < criticalDensity_); diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h index fc2327fb7ae64..d0a0d379e2fb8 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h @@ -99,7 +99,9 @@ namespace ticl { void calculateLocalDensity(const TILES&, const int layerId, const std::vector>&); void calculateDistanceToHigher(const TILES&, const int layerId, const std::vector>&); int findAndAssignTracksters(const TILES&, const std::vector>&); - void dumpClusters(const TILES &tiles, const std::vector>& layerIdx2layerandSoa, const int) const; + void dumpClusters(const TILES& tiles, + const std::vector>& layerIdx2layerandSoa, + const int) const; void dumpTracksters(const std::vector>& layerIdx2layerandSoa, const int, const std::vector&) const; From b81d94073005772a293746face14f3a52f4ba957 Mon Sep 17 00:00:00 2001 From: Marco Rovere Date: Mon, 21 Mar 2022 14:54:20 +0100 Subject: [PATCH 25/25] Cleanup commented code --- RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc | 9 +-------- RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h | 1 - 2 files changed, 1 insertion(+), 9 deletions(-) diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index d99f5b4bcfbf3..16b0ab927f483 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -712,15 +712,8 @@ void PatternRecognitionbyCLUE3D::calculateDistanceToHigher( clustersOnOtherLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i], clustersOnLayer.phi[i], clustersOnOtherLayer.phi[layerandSoa.second]); - // TODO: does that really make sense? - // The real logic should be: - // Among all the LC that falls within my "reachability" distance - // in the XY plane, take the closest in Z with higher local - // density. // Add Z-scale to the final distance - dist = dist_transverse; // + (clustersOnLayer.z[i] - - // clustersOnOtherLayer.z[layerandSoa.second]) * - // (clustersOnLayer.z[i] - clustersOnOtherLayer.z[layerandSoa.second]); + dist = dist_transverse; } else { dist = reco::deltaR2(clustersOnLayer.eta[i], clustersOnLayer.phi[i], diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h index d0a0d379e2fb8..8e77ff55e605e 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h @@ -122,7 +122,6 @@ namespace ticl { const bool useAbsoluteProjectiveScale_; const bool useClusterDimensionXY_; const bool rescaleDensityByZ_; - ; const double criticalEtaPhiDistance_; const double criticalXYDistance_; const int criticalZDistanceLyr_;