From 442d7be3afabe2c2d69ca24d7b628903eb9a6381 Mon Sep 17 00:00:00 2001 From: Dhevan Gangadharan Date: Tue, 18 Apr 2023 17:04:15 -0500 Subject: [PATCH 1/6] Move Cluster reconstruct function out of header and into .cc file. Add 3D cluster shape analysis. --- .../calorimetry/CalorimeterClusterRecoCoG.cc | 160 ++++++++++++++++++ .../calorimetry/CalorimeterClusterRecoCoG.h | 109 +----------- 2 files changed, 166 insertions(+), 103 deletions(-) diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc index 679b6dc0a1..fb812fc1d2 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc @@ -6,6 +6,7 @@ * Logarithmic weighting is used for mimicing energy deposit in transverse direction * * Author: Chao Peng (ANL), 09/27/2020 + Dhevan Gangadharan (UH): cluster profiling from Eigenvalues */ #include "CalorimeterClusterRecoCoG.h" @@ -161,3 +162,162 @@ void CalorimeterClusterRecoCoG::AlgorithmProcess() { return; } + +//------------------------------------------------------------------------ +edm4eic::Cluster* CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoCluster* pcl) const { + edm4eic::MutableCluster cl; + cl.setNhits(pcl->hits_size()); + + m_log->debug("hit size = {}", pcl->hits_size()); + + // no hits + if (pcl->hits_size() == 0) { + return nullptr; + } + + // calculate total energy, find the cell with the maximum energy deposit + float totalE = 0.; + float maxE = 0.; + // Used to optionally constrain the cluster eta to those of the contributing hits + float minHitEta = std::numeric_limits::max(); + float maxHitEta = std::numeric_limits::min(); + auto time = pcl->getHits()[0].getTime(); + auto timeError = pcl->getHits()[0].getTimeError(); + for (unsigned i = 0; i < pcl->getHits().size(); ++i) { + const auto& hit = pcl->getHits()[i]; + const auto weight = pcl->getWeights()[i]; + m_log->debug("hit energy = {} hit weight: {}", hit.getEnergy(), weight); + auto energy = hit.getEnergy() * weight; + totalE += energy; + if (energy > maxE) { + } + const float eta = edm4eic::eta(hit.getPosition()); + if (eta < minHitEta) { + minHitEta = eta; + } + if (eta > maxHitEta) { + maxHitEta = eta; + } + } + cl.setEnergy(totalE / m_sampFrac); + cl.setEnergyError(0.); + cl.setTime(time); + cl.setTimeError(timeError); + + // center of gravity with logarithmic weighting + float tw = 0.; + auto v = cl.getPosition(); + for (unsigned i = 0; i < pcl->getHits().size(); ++i) { + const auto& hit = pcl->getHits()[i]; + const auto weight = pcl->getWeights()[i]; + // _DBG_<<" -- weight = " << weight << " E=" << hit.getEnergy() << " totalE=" < maxHitEta); + const bool underflow = (edm4eic::eta(cl.getPosition()) < minHitEta); + if (overflow || underflow) { + const double newEta = overflow ? maxHitEta : minHitEta; + const double newTheta = edm4eic::etaToAngle(newEta); + const double newR = edm4eic::magnitude(cl.getPosition()); + const double newPhi = edm4eic::angleAzimuthal(cl.getPosition()); + cl.setPosition(edm4eic::sphericalToVector(newR, newTheta, newPhi)); + m_log->debug("Bound cluster position to contributing hits due to {}", (overflow ? "overflow" : "underflow")); + } + } + + // Additional convenience variables + + // best estimate on the cluster direction is the cluster position + // for simple 2D CoG clustering + cl.setIntrinsicTheta(edm4eic::anglePolar(cl.getPosition())); + cl.setIntrinsicPhi(edm4eic::angleAzimuthal(cl.getPosition())); + // TODO errors + + //_______________________________________ + // Calculate cluster profile: + // radius, + // dispersion (energy weighted radius), + // sigma_long + // sigma_short + // sigma_z + double radius = 0, dispersion = 0, lambda_1 = 0, lambda_2 = 0, lambda_3 = 0; + double w_sum = 0; + double sum_11 = 0, sum_22 = 0, sum_33 = 0; + double sum_12 = 0, sum_13 = 0, sum_23 = 0; + double sum_1 = 0, sum_2 = 0, sum_3 = 0; + + if (cl.getNhits() > 1) { + + for (const auto& hit : pcl->getHits()) { + + const auto delta = cl.getPosition() - hit.getPosition(); + radius += delta * delta; + + double w = weightFunc(hit.getEnergy(), cl.getEnergy(), m_logWeightBase, 0); + w_sum += w; + dispersion += delta * delta * w; + + double pos_1 = edm4eic::anglePolar( hit.getPosition() ); + double pos_2 = edm4eic::angleAzimuthal( hit.getPosition() ); + double pos_3 = hit.getPosition().z; + + if( m_xyClusterProfiling ) { + pos_1 = hit.getPosition().x; + pos_2 = hit.getPosition().y; + } + sum_11 += w * pos_1 * pos_1; + sum_22 += w * pos_2 * pos_2; + sum_33 += w * pos_3 * pos_3; + sum_12 += w * pos_1 * pos_2; + sum_13 += w * pos_1 * pos_3; + sum_23 += w * pos_2 * pos_3; + sum_1 += w * pos_1; + sum_2 += w * pos_2; + sum_3 += w * pos_3; + } + + if( w_sum > 0 ) { + radius = sqrt((1. / (cl.getNhits() - 1.)) * radius); + dispersion = sqrt( dispersion / w_sum ); + + // variances and covariances + double sigma_11 = sum_11 / w_sum - (sum_1/w_sum) * (sum_1/w_sum); + double sigma_22 = sum_22 / w_sum - (sum_2/w_sum) * (sum_2/w_sum); + double sigma_33 = sum_33 / w_sum - (sum_3/w_sum) * (sum_3/w_sum); + double sigma_12 = sum_12 / w_sum - (sum_1/w_sum) * (sum_2/w_sum); + double sigma_13 = sum_13 / w_sum - (sum_1/w_sum) * (sum_3/w_sum); + double sigma_23 = sum_23 / w_sum - (sum_2/w_sum) * (sum_3/w_sum); + + // covariance matrix + Eigen::MatrixXd cov3(3,3); + cov3(0,0) = sigma_11; cov3(1,1) = sigma_22; cov3(2,2) = sigma_33; + cov3(0,1) = sigma_12; cov3(0,2) = sigma_13; cov3(1,2) = sigma_23; + cov3(1,0) = cov3(0,1); cov3(2,0) = cov3(0,2); cov3(2,1) = cov3(1,2); + + // Eigenvalues correspond to cluster's 2nd moments (sigma_long, sigma_short, sigma_z) + Eigen::EigenSolver es(cov3, false); // set to true for eigenvector calculation + auto eigenValues = es.eigenvalues(); + lambda_1 = eigenValues[0].real(); // imaginary parts correspond to unphysical roots + lambda_2 = eigenValues[1].real(); + lambda_3 = eigenValues[2].real(); + } + } + + cl.addToShapeParameters( radius ); + cl.addToShapeParameters( dispersion ); + cl.addToShapeParameters( lambda_1 ); + cl.addToShapeParameters( lambda_2 ); + cl.addToShapeParameters( lambda_3 ); + + return new edm4eic::Cluster(cl); +} diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h index 6917dfa83a..ec4b007f1c 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h @@ -13,6 +13,8 @@ #include +#include + #include #include #include @@ -67,6 +69,7 @@ class CalorimeterClusterRecoCoG { // the eta of the contributing hits. This is useful to avoid edge effects // for endcaps. bool m_enableEtaBounds = false;//{this, "enableEtaBounds", false}; + bool m_xyClusterProfiling = false;//{this, "xyClusterProfiling, false}; std::shared_ptr m_geoSvc; @@ -83,108 +86,8 @@ class CalorimeterClusterRecoCoG { std::vector m_outputClusters; std::vector m_outputAssociations; - private: -edm4eic::Cluster* reconstruct(const edm4eic::ProtoCluster* pcl) const { - edm4eic::MutableCluster cl; - cl.setNhits(pcl->hits_size()); - - m_log->debug("hit size = {}", pcl->hits_size()); - - // no hits - if (pcl->hits_size() == 0) { - return nullptr; - } - - // calculate total energy, find the cell with the maximum energy deposit - float totalE = 0.; - float maxE = 0.; - // Used to optionally constrain the cluster eta to those of the contributing hits - float minHitEta = std::numeric_limits::max(); - float maxHitEta = std::numeric_limits::min(); - auto time = pcl->getHits()[0].getTime(); - auto timeError = pcl->getHits()[0].getTimeError(); - for (unsigned i = 0; i < pcl->getHits().size(); ++i) { - const auto& hit = pcl->getHits()[i]; - const auto weight = pcl->getWeights()[i]; - m_log->debug("hit energy = {} hit weight: {}", hit.getEnergy(), weight); - auto energy = hit.getEnergy() * weight; - totalE += energy; - if (energy > maxE) { - } - const float eta = edm4eic::eta(hit.getPosition()); - if (eta < minHitEta) { - minHitEta = eta; - } - if (eta > maxHitEta) { - maxHitEta = eta; - } - } - cl.setEnergy(totalE / m_sampFrac); - cl.setEnergyError(0.); - - // center of gravity with logarithmic weighting - float tw = 0.; - float timeWeighted = 0.; - float timeErrorWeighted = 0.; - auto v = cl.getPosition(); - for (unsigned i = 0; i < pcl->getHits().size(); ++i) { - const auto& hit = pcl->getHits()[i]; - const auto weight = pcl->getWeights()[i]; -// _DBG_<<" -- weight = " << weight << " E=" << hit.getEnergy() << " totalE=" < maxHitEta); - const bool underflow = (edm4eic::eta(cl.getPosition()) < minHitEta); - if (overflow || underflow) { - const double newEta = overflow ? maxHitEta : minHitEta; - const double newTheta = edm4eic::etaToAngle(newEta); - const double newR = edm4eic::magnitude(cl.getPosition()); - const double newPhi = edm4eic::angleAzimuthal(cl.getPosition()); - cl.setPosition(edm4eic::sphericalToVector(newR, newTheta, newPhi)); - m_log->debug("Bound cluster position to contributing hits due to {}", (overflow ? "overflow" : "underflow")); - } - } - - // Additional convenience variables - - // best estimate on the cluster direction is the cluster position - // for simple 2D CoG clustering - cl.setIntrinsicTheta(edm4eic::anglePolar(cl.getPosition())); - cl.setIntrinsicPhi(edm4eic::angleAzimuthal(cl.getPosition())); - // TODO errors - - // Calculate radius - // @TODO: add skewness - if (cl.getNhits() > 1) { - double radius = 0; - for (const auto& hit : pcl->getHits()) { - const auto delta = cl.getPosition() - hit.getPosition(); - radius += delta * delta; - } - radius = sqrt((1. / (cl.getNhits() - 1.)) * radius); - cl.addToShapeParameters(radius); - cl.addToShapeParameters(0 /* skewness */); // skewness not yet calculated - } -// edm4eic::Cluster c(cl); - return new edm4eic::Cluster(cl); - } - - + + edm4eic::Cluster* reconstruct(const edm4eic::ProtoCluster* pcl) const; + }; From 77063293fbd991903cbb1fb1e9b982b9e1c3b7ae Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 18 Apr 2023 22:34:04 +0000 Subject: [PATCH 2/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc | 8 ++++---- src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc index fb812fc1d2..28cae7ad15 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc @@ -245,7 +245,7 @@ edm4eic::Cluster* CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoClu //_______________________________________ // Calculate cluster profile: - // radius, + // radius, // dispersion (energy weighted radius), // sigma_long // sigma_short @@ -281,7 +281,7 @@ edm4eic::Cluster* CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoClu sum_12 += w * pos_1 * pos_2; sum_13 += w * pos_1 * pos_3; sum_23 += w * pos_2 * pos_3; - sum_1 += w * pos_1; + sum_1 += w * pos_1; sum_2 += w * pos_2; sum_3 += w * pos_3; } @@ -290,7 +290,7 @@ edm4eic::Cluster* CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoClu radius = sqrt((1. / (cl.getNhits() - 1.)) * radius); dispersion = sqrt( dispersion / w_sum ); - // variances and covariances + // variances and covariances double sigma_11 = sum_11 / w_sum - (sum_1/w_sum) * (sum_1/w_sum); double sigma_22 = sum_22 / w_sum - (sum_2/w_sum) * (sum_2/w_sum); double sigma_33 = sum_33 / w_sum - (sum_3/w_sum) * (sum_3/w_sum); @@ -311,7 +311,7 @@ edm4eic::Cluster* CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoClu lambda_2 = eigenValues[1].real(); lambda_3 = eigenValues[2].real(); } - } + } cl.addToShapeParameters( radius ); cl.addToShapeParameters( dispersion ); diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h index ec4b007f1c..6bab764d36 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h @@ -87,7 +87,7 @@ class CalorimeterClusterRecoCoG { std::vector m_outputAssociations; private: - + edm4eic::Cluster* reconstruct(const edm4eic::ProtoCluster* pcl) const; - + }; From 56dcc652bafc8f3cb87e4d9ec4bd5984103796ed Mon Sep 17 00:00:00 2001 From: Dhevan Gangadharan Date: Wed, 19 Apr 2023 14:39:13 -0500 Subject: [PATCH 3/6] Simplify eigenvalue calculation. Some other minor comments from the PR implemented. --- .../calorimetry/CalorimeterClusterRecoCoG.cc | 90 +++++++++---------- .../calorimetry/CalorimeterClusterRecoCoG.h | 2 - 2 files changed, 42 insertions(+), 50 deletions(-) diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc index 28cae7ad15..b9f2ffe1a6 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc @@ -1,24 +1,23 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Sylvester Joosten, Chao, Chao Peng, Whitney Armstrong +// Copyright (C) 2022 Sylvester Joosten, Chao, Chao Peng, Whitney Armstrong, Dhevan Gangadharan /* * Reconstruct the cluster with Center of Gravity method * Logarithmic weighting is used for mimicing energy deposit in transverse direction * * Author: Chao Peng (ANL), 09/27/2020 - Dhevan Gangadharan (UH): cluster profiling from Eigenvalues */ -#include "CalorimeterClusterRecoCoG.h" - -#include -#include +#include +#include #include #include +#include -#include -#include +#include +#include #include +#include "CalorimeterClusterRecoCoG.h" using namespace dd4hep; @@ -252,62 +251,57 @@ edm4eic::Cluster* CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoClu // sigma_z double radius = 0, dispersion = 0, lambda_1 = 0, lambda_2 = 0, lambda_3 = 0; double w_sum = 0; - double sum_11 = 0, sum_22 = 0, sum_33 = 0; - double sum_12 = 0, sum_13 = 0, sum_23 = 0; - double sum_1 = 0, sum_2 = 0, sum_3 = 0; + + Eigen::Matrix3f sum2 = Eigen::Matrix3f::Zero(); + Eigen::Vector3f sum1 = Eigen::Vector3f::Zero(); if (cl.getNhits() > 1) { for (const auto& hit : pcl->getHits()) { - const auto delta = cl.getPosition() - hit.getPosition(); - radius += delta * delta; - - double w = weightFunc(hit.getEnergy(), cl.getEnergy(), m_logWeightBase, 0); - w_sum += w; - dispersion += delta * delta * w; + float w = weightFunc(hit.getEnergy(), cl.getEnergy(), m_logWeightBase, 0); - double pos_1 = edm4eic::anglePolar( hit.getPosition() ); - double pos_2 = edm4eic::angleAzimuthal( hit.getPosition() ); - double pos_3 = hit.getPosition().z; + float pos_1 = edm4eic::anglePolar( hit.getPosition() ); + float pos_2 = edm4eic::angleAzimuthal( hit.getPosition() ); + float pos_3 = hit.getPosition().z; if( m_xyClusterProfiling ) { pos_1 = hit.getPosition().x; pos_2 = hit.getPosition().y; } - sum_11 += w * pos_1 * pos_1; - sum_22 += w * pos_2 * pos_2; - sum_33 += w * pos_3 * pos_3; - sum_12 += w * pos_1 * pos_2; - sum_13 += w * pos_1 * pos_3; - sum_23 += w * pos_2 * pos_3; - sum_1 += w * pos_1; - sum_2 += w * pos_2; - sum_3 += w * pos_3; + + const auto delta = cl.getPosition() - hit.getPosition(); + Eigen::Vector3f pos(pos_1, pos_2, pos_3); + + radius += delta * delta; + + dispersion += delta * delta * w; + + // Weighted Sum x*x, x*y, x*z, y*y, etc. + sum2 += w * pos * pos.transpose(); + + // Weighted Sum x, y, z + sum1 += w * pos; + + w_sum += w; } if( w_sum > 0 ) { radius = sqrt((1. / (cl.getNhits() - 1.)) * radius); dispersion = sqrt( dispersion / w_sum ); - - // variances and covariances - double sigma_11 = sum_11 / w_sum - (sum_1/w_sum) * (sum_1/w_sum); - double sigma_22 = sum_22 / w_sum - (sum_2/w_sum) * (sum_2/w_sum); - double sigma_33 = sum_33 / w_sum - (sum_3/w_sum) * (sum_3/w_sum); - double sigma_12 = sum_12 / w_sum - (sum_1/w_sum) * (sum_2/w_sum); - double sigma_13 = sum_13 / w_sum - (sum_1/w_sum) * (sum_3/w_sum); - double sigma_23 = sum_23 / w_sum - (sum_2/w_sum) * (sum_3/w_sum); - - // covariance matrix - Eigen::MatrixXd cov3(3,3); - cov3(0,0) = sigma_11; cov3(1,1) = sigma_22; cov3(2,2) = sigma_33; - cov3(0,1) = sigma_12; cov3(0,2) = sigma_13; cov3(1,2) = sigma_23; - cov3(1,0) = cov3(0,1); cov3(2,0) = cov3(0,2); cov3(2,1) = cov3(1,2); - - // Eigenvalues correspond to cluster's 2nd moments (sigma_long, sigma_short, sigma_z) - Eigen::EigenSolver es(cov3, false); // set to true for eigenvector calculation - auto eigenValues = es.eigenvalues(); - lambda_1 = eigenValues[0].real(); // imaginary parts correspond to unphysical roots + + // normalize matrix and vector + sum2 /= w_sum; + sum1 /= w_sum; + + // 3D covariance matrix + Eigen::Matrix3f cov3 = sum2 - sum1 * sum1.transpose(); + + // Solve for eigenvalues. Corresponds to cluster's 2nd moments (sigma_long, sigma_short, sigma_z) + Eigen::EigenSolver es(cov3, false); // set to true for eigenvector calculation + + auto eigenValues = es.eigenvalues(); // eigenvalues of symmetric real matrix are always real + lambda_1 = eigenValues[0].real(); lambda_2 = eigenValues[1].real(); lambda_3 = eigenValues[2].real(); } diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h index 6bab764d36..af5c866547 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h @@ -13,8 +13,6 @@ #include -#include - #include #include #include From c29717bffac74cc112b26a9fb92618e7112ab36d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 19 Apr 2023 19:39:22 +0000 Subject: [PATCH 4/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../calorimetry/CalorimeterClusterRecoCoG.cc | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc index b9f2ffe1a6..df1720ab59 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc @@ -269,10 +269,10 @@ edm4eic::Cluster* CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoClu pos_1 = hit.getPosition().x; pos_2 = hit.getPosition().y; } - + const auto delta = cl.getPosition() - hit.getPosition(); - Eigen::Vector3f pos(pos_1, pos_2, pos_3); - + Eigen::Vector3f pos(pos_1, pos_2, pos_3); + radius += delta * delta; dispersion += delta * delta * w; @@ -289,17 +289,17 @@ edm4eic::Cluster* CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoClu if( w_sum > 0 ) { radius = sqrt((1. / (cl.getNhits() - 1.)) * radius); dispersion = sqrt( dispersion / w_sum ); - + // normalize matrix and vector sum2 /= w_sum; sum1 /= w_sum; - + // 3D covariance matrix Eigen::Matrix3f cov3 = sum2 - sum1 * sum1.transpose(); // Solve for eigenvalues. Corresponds to cluster's 2nd moments (sigma_long, sigma_short, sigma_z) Eigen::EigenSolver es(cov3, false); // set to true for eigenvector calculation - + auto eigenValues = es.eigenvalues(); // eigenvalues of symmetric real matrix are always real lambda_1 = eigenValues[0].real(); lambda_2 = eigenValues[1].real(); From 234522c29ff9c11f4dbfdd250698178c3dbe0b1b Mon Sep 17 00:00:00 2001 From: Dhevan Gangadharan Date: Thu, 20 Apr 2023 11:01:42 -0500 Subject: [PATCH 5/6] Change theta-phi-z analysis to just theta-phi and add the 2 eigenvalues to the returend cluster shape parameters. --- .../calorimetry/CalorimeterClusterRecoCoG.cc | 76 ++++++++++--------- .../calorimetry/CalorimeterClusterRecoCoG.h | 1 - 2 files changed, 39 insertions(+), 38 deletions(-) diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc index df1720ab59..0267a5714e 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc @@ -246,14 +246,16 @@ edm4eic::Cluster* CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoClu // Calculate cluster profile: // radius, // dispersion (energy weighted radius), - // sigma_long - // sigma_short - // sigma_z - double radius = 0, dispersion = 0, lambda_1 = 0, lambda_2 = 0, lambda_3 = 0; - double w_sum = 0; + // eta-phi cluster widths (2D) + // x-y-z cluster widths (3D) + float radius = 0, dispersion = 0, w_sum = 0; - Eigen::Matrix3f sum2 = Eigen::Matrix3f::Zero(); - Eigen::Vector3f sum1 = Eigen::Vector3f::Zero(); + Eigen::Matrix2f sum2_2D = Eigen::Matrix2f::Zero(); + Eigen::Matrix3f sum2_3D = Eigen::Matrix3f::Zero(); + Eigen::Vector2f sum1_2D = Eigen::Vector2f::Zero(); + Eigen::Vector3f sum1_3D = Eigen::Vector3f::Zero(); + Eigen::Vector2cf eigenValues_2D = Eigen::Vector2cf::Zero(); + Eigen::Vector3cf eigenValues_3D = Eigen::Vector3cf::Zero(); if (cl.getNhits() > 1) { @@ -261,57 +263,57 @@ edm4eic::Cluster* CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoClu float w = weightFunc(hit.getEnergy(), cl.getEnergy(), m_logWeightBase, 0); - float pos_1 = edm4eic::anglePolar( hit.getPosition() ); - float pos_2 = edm4eic::angleAzimuthal( hit.getPosition() ); - float pos_3 = hit.getPosition().z; - - if( m_xyClusterProfiling ) { - pos_1 = hit.getPosition().x; - pos_2 = hit.getPosition().y; - } + // theta, phi + Eigen::Vector2f pos2D( edm4eic::anglePolar( hit.getPosition() ), edm4eic::angleAzimuthal( hit.getPosition() ) ); + // x, y, z + Eigen::Vector3f pos3D( hit.getPosition().x, hit.getPosition().y, hit.getPosition().z ); const auto delta = cl.getPosition() - hit.getPosition(); - Eigen::Vector3f pos(pos_1, pos_2, pos_3); - - radius += delta * delta; - - dispersion += delta * delta * w; + radius += delta * delta; + dispersion += delta * delta * w; // Weighted Sum x*x, x*y, x*z, y*y, etc. - sum2 += w * pos * pos.transpose(); + sum2_2D += w * pos2D * pos2D.transpose(); + sum2_3D += w * pos3D * pos3D.transpose(); // Weighted Sum x, y, z - sum1 += w * pos; + sum1_2D += w * pos2D; + sum1_3D += w * pos3D; w_sum += w; } if( w_sum > 0 ) { - radius = sqrt((1. / (cl.getNhits() - 1.)) * radius); + radius = sqrt((1. / (cl.getNhits() - 1.)) * radius); dispersion = sqrt( dispersion / w_sum ); - // normalize matrix and vector - sum2 /= w_sum; - sum1 /= w_sum; + // normalize matrices + sum2_2D /= w_sum; + sum2_3D /= w_sum; + sum1_2D /= w_sum; + sum1_3D /= w_sum; - // 3D covariance matrix - Eigen::Matrix3f cov3 = sum2 - sum1 * sum1.transpose(); + // 2D and 3D covariance matrices + Eigen::Matrix2f cov2 = sum2_2D - sum1_2D * sum1_2D.transpose(); + Eigen::Matrix3f cov3 = sum2_3D - sum1_3D * sum1_3D.transpose(); - // Solve for eigenvalues. Corresponds to cluster's 2nd moments (sigma_long, sigma_short, sigma_z) - Eigen::EigenSolver es(cov3, false); // set to true for eigenvector calculation + // Solve for eigenvalues. Corresponds to cluster's 2nd moments (widths) + Eigen::EigenSolver es_2D(cov2, false); // set to true for eigenvector calculation + Eigen::EigenSolver es_3D(cov3, false); // set to true for eigenvector calculation - auto eigenValues = es.eigenvalues(); // eigenvalues of symmetric real matrix are always real - lambda_1 = eigenValues[0].real(); - lambda_2 = eigenValues[1].real(); - lambda_3 = eigenValues[2].real(); + // eigenvalues of symmetric real matrix are always real + eigenValues_2D = es_2D.eigenvalues(); + eigenValues_3D = es_3D.eigenvalues(); } } cl.addToShapeParameters( radius ); cl.addToShapeParameters( dispersion ); - cl.addToShapeParameters( lambda_1 ); - cl.addToShapeParameters( lambda_2 ); - cl.addToShapeParameters( lambda_3 ); + cl.addToShapeParameters( eigenValues_2D[0].real() ); // 2D eta-phi cluster width 1 + cl.addToShapeParameters( eigenValues_2D[1].real() ); // 2D eta-phi cluster width 2 + cl.addToShapeParameters( eigenValues_3D[0].real() ); // 3D x-y-z cluster width 1 + cl.addToShapeParameters( eigenValues_3D[1].real() ); // 3D x-y-z cluster width 2 + cl.addToShapeParameters( eigenValues_3D[2].real() ); // 3D x-y-z cluster width 3 return new edm4eic::Cluster(cl); } diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h index af5c866547..903aafba1a 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h @@ -67,7 +67,6 @@ class CalorimeterClusterRecoCoG { // the eta of the contributing hits. This is useful to avoid edge effects // for endcaps. bool m_enableEtaBounds = false;//{this, "enableEtaBounds", false}; - bool m_xyClusterProfiling = false;//{this, "xyClusterProfiling, false}; std::shared_ptr m_geoSvc; From 59efc01a67683c151ebe88cea68ce2e3394a95a3 Mon Sep 17 00:00:00 2001 From: Dhevan Gangadharan Date: Thu, 20 Apr 2023 14:53:00 -0500 Subject: [PATCH 6/6] Fix misprint in comments (eta -> theta). --- src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc index 0267a5714e..166afd7e97 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc @@ -246,7 +246,7 @@ edm4eic::Cluster* CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoClu // Calculate cluster profile: // radius, // dispersion (energy weighted radius), - // eta-phi cluster widths (2D) + // theta-phi cluster widths (2D) // x-y-z cluster widths (3D) float radius = 0, dispersion = 0, w_sum = 0; @@ -309,8 +309,8 @@ edm4eic::Cluster* CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoClu cl.addToShapeParameters( radius ); cl.addToShapeParameters( dispersion ); - cl.addToShapeParameters( eigenValues_2D[0].real() ); // 2D eta-phi cluster width 1 - cl.addToShapeParameters( eigenValues_2D[1].real() ); // 2D eta-phi cluster width 2 + cl.addToShapeParameters( eigenValues_2D[0].real() ); // 2D theta-phi cluster width 1 + cl.addToShapeParameters( eigenValues_2D[1].real() ); // 2D theta-phi cluster width 2 cl.addToShapeParameters( eigenValues_3D[0].real() ); // 3D x-y-z cluster width 1 cl.addToShapeParameters( eigenValues_3D[1].real() ); // 3D x-y-z cluster width 2 cl.addToShapeParameters( eigenValues_3D[2].real() ); // 3D x-y-z cluster width 3