From 80866b3335cfeff59a47423a31c875b81ed492b7 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Fri, 26 Jul 2024 15:18:27 -0400 Subject: [PATCH] CalorimeterClusterRecoCoG: set intrinsic{Phi,Theta} (#1511) This partially reverts effect of #1391. The clusterShape has become a way to extend EDM without going through extending it. We should sort those variable out into proper members, this is a first step. --------- Co-authored-by: Wouter Deconinck --- .../calorimetry/CalorimeterClusterRecoCoG.cc | 49 ++++++++++--------- .../CalorimeterClusterRecoCoGConfig.h | 2 + src/detectors/BEMC/BEMC.cc | 1 + src/detectors/FHCAL/FHCAL.cc | 4 ++ src/detectors/ZDC/ZDC.cc | 5 ++ .../CalorimeterClusterRecoCoG_factory.h | 1 + .../calorimetry_CalorimeterClusterRecoCoG.cc | 38 ++++++-------- 7 files changed, 53 insertions(+), 47 deletions(-) diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc index 60716fc61e..151483ac67 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc @@ -228,12 +228,6 @@ std::optional CalorimeterClusterRecoCoG::reconstruct(co // Additional convenience variables - // best estimate on the cluster direction is the cluster position - // for simple 2D CoG clustering - cl.setIntrinsicTheta(edm4hep::utils::anglePolar(cl.getPosition())); - cl.setIntrinsicPhi(edm4hep::utils::angleAzimuthal(cl.getPosition())); - // TODO errors - //_______________________________________ // Calculate cluster profile: // radius, @@ -248,12 +242,11 @@ std::optional CalorimeterClusterRecoCoG::reconstruct(co Eigen::Vector3f sum1_3D = Eigen::Vector3f::Zero(); Eigen::Vector2cf eigenValues_2D = Eigen::Vector2cf::Zero(); Eigen::Vector3cf eigenValues_3D = Eigen::Vector3cf::Zero(); - //the axis is the direction of the eigenvalue corresponding to the largest eigenvalue. - double axis_x=0, axis_y=0, axis_z=0; - if (cl.getNhits() > 1) { + // the axis is the direction of the eigenvalue corresponding to the largest eigenvalue. + edm4hep::Vector3f axis; + if (cl.getNhits() > 1) { for (const auto& hit : pcl.getHits()) { - float w = weightFunc(hit.getEnergy(), totalE, logWeightBase, 0); // theta, phi @@ -306,19 +299,15 @@ std::optional CalorimeterClusterRecoCoG::reconstruct(co return std::real(a) < std::real(b); } ); - auto axis = eigenvectors.col(std::distance( + auto axis_eigen = eigenvectors.col(std::distance( eigenValues_3D.begin(), max_eigenvalue_it )); - axis_x=axis(0,0).real(); - axis_y=axis(1,0).real(); - axis_z=axis(2,0).real(); - double norm=sqrt(axis_x*axis_x+axis_y*axis_y+axis_z*axis_z); - if (norm!=0){ - axis_x/=norm; - axis_y/=norm; - axis_z/=norm; - } + axis = { + axis_eigen(0,0).real(), + axis_eigen(1,0).real(), + axis_eigen(2,0).real(), + }; } } @@ -329,10 +318,22 @@ std::optional CalorimeterClusterRecoCoG::reconstruct(co 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 - //last 3 shape parameters are the components of the axis direction - cl.addToShapeParameters( axis_x ); - cl.addToShapeParameters( axis_y ); - cl.addToShapeParameters( axis_z ); + + + double dot_product = cl.getPosition() * axis; + if (dot_product < 0) { + axis = -1 * axis; + } + + if (m_cfg.longitudinalShowerInfoAvailable) { + cl.setIntrinsicTheta(edm4hep::utils::anglePolar(axis)); + cl.setIntrinsicPhi(edm4hep::utils::angleAzimuthal(axis)); + // TODO intrinsicDirectionError + } else { + cl.setIntrinsicTheta(NAN); + cl.setIntrinsicPhi(NAN); + } + return std::move(cl); } diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h b/src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h index f6ea18e2e2..1b0624462e 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h @@ -23,6 +23,8 @@ namespace eicrecon { std::vector logWeightBaseCoeffs{}; double logWeightBase_Eref = 50 * dd4hep::MeV; + bool longitudinalShowerInfoAvailable = false; + // Constrain the cluster position eta to be within // the eta of the contributing hits. This is useful to avoid edge effects // for endcaps. diff --git a/src/detectors/BEMC/BEMC.cc b/src/detectors/BEMC/BEMC.cc index 10dbd0bee2..2d652a5ed6 100644 --- a/src/detectors/BEMC/BEMC.cc +++ b/src/detectors/BEMC/BEMC.cc @@ -97,6 +97,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 6.2, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false }, app // TODO: Remove me once fixed diff --git a/src/detectors/FHCAL/FHCAL.cc b/src/detectors/FHCAL/FHCAL.cc index 2af4015ba5..0d78240fce 100644 --- a/src/detectors/FHCAL/FHCAL.cc +++ b/src/detectors/FHCAL/FHCAL.cc @@ -117,6 +117,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 0.0257, .logWeightBase = 3.6, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = true }, app // TODO: Remove me once fixed @@ -134,6 +135,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 0.0257, .logWeightBase = 6.2, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false, }, app // TODO: Remove me once fixed @@ -230,6 +232,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 4.5, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false }, app // TODO: Remove me once fixed @@ -247,6 +250,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 4.5, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false, }, app // TODO: Remove me once fixed diff --git a/src/detectors/ZDC/ZDC.cc b/src/detectors/ZDC/ZDC.cc index fb88ee0305..3a3ae384ea 100644 --- a/src/detectors/ZDC/ZDC.cc +++ b/src/detectors/ZDC/ZDC.cc @@ -85,6 +85,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 3.6, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false }, app // TODO: Remove me once fixed @@ -102,6 +103,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 6.2, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false, }, app // TODO: Remove me once fixed @@ -193,6 +195,7 @@ extern "C" { .sampFrac = 0.0203, .logWeightBaseCoeffs={5.0,0.65,0.31}, .logWeightBase_Eref=50*dd4hep::GeV, + .longitudinalShowerInfoAvailable = true, }, app // TODO: Remove me once fixed )); @@ -228,6 +231,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 3.6, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false }, app // TODO: Remove me once fixed @@ -244,6 +248,7 @@ extern "C" { .energyWeight = "log", .sampFrac = 0.0203, .logWeightBase = 6.2, + .longitudinalShowerInfoAvailable = true, .enableEtaBounds = false, }, app // TODO: Remove me once fixed diff --git a/src/factories/calorimetry/CalorimeterClusterRecoCoG_factory.h b/src/factories/calorimetry/CalorimeterClusterRecoCoG_factory.h index 7d4b275e38..5dc7a6a14d 100644 --- a/src/factories/calorimetry/CalorimeterClusterRecoCoG_factory.h +++ b/src/factories/calorimetry/CalorimeterClusterRecoCoG_factory.h @@ -29,6 +29,7 @@ class CalorimeterClusterRecoCoG_factory : public JOmniFactory m_logWeightBase {this, "logWeightBase", config().logWeightBase}; ParameterRef> m_logWeightBaseCoeffs {this, "logWeightBaseCoeffs", config().logWeightBaseCoeffs}; ParameterRef m_logWeightBase_Eref {this, "logWeightBase_Eref", config().logWeightBase_Eref}; + ParameterRef m_longitudinalShowerInfoAvailable {this, "longitudinalShowerInfoAvailable", config().longitudinalShowerInfoAvailable}; ParameterRef m_enableEtaBounds {this, "enableEtaBounds", config().enableEtaBounds}; Service m_algorithmsInit {this}; diff --git a/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc b/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc index 32f11942d3..86358ee4b5 100644 --- a/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc +++ b/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc @@ -4,6 +4,8 @@ #include // for MeV, mm, keV, ns #include // for AssertionHandler, operator""_catch_sr, StringRef, REQUIRE, operator<, operator==, operator>, TEST_CASE +#include +#include #include // for CalorimeterHitCollection, MutableCalorimeterHit, CalorimeterHitMutableCollectionIterator #include #include @@ -11,7 +13,6 @@ #include #include // for Vector3f #include -#include #include // for level_enum #include // for logger #include // for default_logger @@ -28,6 +29,8 @@ using eicrecon::CalorimeterClusterRecoCoGConfig; using edm4eic::CalorimeterHit; TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) { + const float EPSILON = 1e-5; + CalorimeterClusterRecoCoG algo("CalorimeterClusterRecoCoG"); std::shared_ptr logger = spdlog::default_logger()->clone("CalorimeterClusterRecoCoG"); @@ -35,9 +38,10 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) CalorimeterClusterRecoCoGConfig cfg; cfg.energyWeight = "log"; - cfg.sampFrac = 0.0203, - cfg.logWeightBaseCoeffs={5.0,0.65,0.31}, - cfg.logWeightBase_Eref=50*dd4hep::GeV, + cfg.sampFrac = 0.0203; + cfg.logWeightBaseCoeffs = {5.0,0.65,0.31}; + cfg.logWeightBase_Eref = 50*dd4hep::GeV; + cfg.longitudinalShowerInfoAvailable = true; algo.applyConfig(cfg); @@ -51,7 +55,7 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) //create a protocluster with 3 hits auto pclust = pclust_coll.create(); - edm4hep::Vector3f position({0,0,0*dd4hep::mm}); + edm4hep::Vector3f position({0, 0, 1 *dd4hep::mm}); auto hit1 = hits_coll.create(); hit1.setCellID(0); @@ -63,8 +67,9 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) hit1.setDimension({0,0,0}); hit1.setLocal(position); pclust.addToHits(hit1); + pclust.addToWeights(1); - position={0,0, 1*dd4hep::mm}; + position={-1 * dd4hep::mm, 0, 2 * dd4hep::mm}; auto hit2 = hits_coll.create(); hit2.setCellID(0); hit2.setEnergy(0.1*dd4hep::GeV); @@ -75,19 +80,7 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) hit2.setDimension({0,0,0}); hit2.setLocal(position); pclust.addToHits(hit2); - - position={0,0, 2*dd4hep::mm}; - auto hit3 = hits_coll.create();//0, 0.1*dd4hep::GeV, 0,0,0,position, {0,0,0}, 0,0, position); - hit3.setCellID(0); - hit3.setEnergy(0.1*dd4hep::GeV); - hit3.setEnergyError(0); - hit3.setTime(0); - hit3.setTimeError(0); - hit3.setPosition(position); - hit3.setDimension({0,0,0}); - hit3.setLocal(position); - pclust.addToHits(hit3); - pclust.addToWeights(1);pclust.addToWeights(1);pclust.addToWeights(1); + pclust.addToWeights(1); // Constructing input and output as per the algorithm's expected signature auto input = std::make_tuple(&pclust_coll, &simhits); @@ -97,10 +90,9 @@ TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) for (auto clust : *clust_coll){ - // require that this cluster's axis is 0,0,1 - REQUIRE(clust.getShapeParameters()[7] == 0); - REQUIRE(clust.getShapeParameters()[8] == 0); - REQUIRE(abs(clust.getShapeParameters()[9]) == 1); + REQUIRE_THAT(clust.getIntrinsicTheta(), Catch::Matchers::WithinAbs(M_PI / 4, EPSILON)); + // std::abs() checks if we land on -M_PI + REQUIRE_THAT(std::abs(clust.getIntrinsicPhi()), Catch::Matchers::WithinAbs(M_PI, EPSILON)); }