From 8b0de86aec9a2d2f913fc8343d193afbcaff0bb8 Mon Sep 17 00:00:00 2001 From: Sylvester Joosten Date: Tue, 20 Sep 2022 21:09:39 +0000 Subject: [PATCH] First working algorithm implementation, ready to finish Juggler end --- external/algorithms/CMakeLists.txt | 1 + .../algorithms/calorimetry/CMakeLists.txt | 45 ++ .../algorithms/calorimetry/ClusterRecoCoG.h | 68 +++ .../calorimetry/src/ClusterRecoCoG.cpp | 522 +++++++----------- .../core/include/algorithms/algorithm.h | 16 +- .../core/include/algorithms/property.h | 8 +- 6 files changed, 335 insertions(+), 325 deletions(-) create mode 100644 external/algorithms/calorimetry/CMakeLists.txt create mode 100644 external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h diff --git a/external/algorithms/CMakeLists.txt b/external/algorithms/CMakeLists.txt index 21217c4..bb85179 100644 --- a/external/algorithms/CMakeLists.txt +++ b/external/algorithms/CMakeLists.txt @@ -60,3 +60,4 @@ find_package(fmt REQUIRED) include(GNUInstallDirs) add_subdirectory(core) +add_subdirectory(calorimetry) diff --git a/external/algorithms/calorimetry/CMakeLists.txt b/external/algorithms/calorimetry/CMakeLists.txt new file mode 100644 index 0000000..63890cc --- /dev/null +++ b/external/algorithms/calorimetry/CMakeLists.txt @@ -0,0 +1,45 @@ +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (C) 2022 Sylvester Joosten + +################################################################################ +# Package: algorithms core utilities +################################################################################ + +set(SUBDIR "calorimetry") +set(LIBRARY "algo${SUBDIR}") +set(TARGETS ${TARGETS} ${LIBRARY} PARENT_SCOPE) + +file(GLOB SRC CONFIGURE_DEPENDS src/*.cpp) + +add_library(${LIBRARY} SHARED ${SRC}) +target_link_libraries(${LIBRARY} + PUBLIC + EDM4HEP::edm4hep + EDM4EIC::edm4eic + DD4hep::DDRec + algocore + fmt::fmt) +target_include_directories(${LIBRARY} + PUBLIC + $ + $) +set_target_properties(${LIBRARY} PROPERTIES + VERSION ${PROJECT_VERSION} + SOVERSION ${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR}) + +install(TARGETS ${LIBRARY} + EXPORT algorithmsTargets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + ARCHIVE DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT lib + INCLUDES DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}" + NAMESPACE algorithms::) + +install(DIRECTORY ${PROJECT_SOURCE_DIR}/${SUBDIR}/include/algorithms +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} COMPONENT dev) + +# TODO: Testing +#if(BUILD_TESTING) +# enable_testing() +#endif() + diff --git a/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h b/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h new file mode 100644 index 0000000..1f0b821 --- /dev/null +++ b/external/algorithms/calorimetry/include/algorithms/calorimetry/ClusterRecoCoG.h @@ -0,0 +1,68 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Sylvester Joosten, Chao, Chao Peng, Whitney Armstrong + +/* + * Reconstruct the cluster with Center of Gravity method + * Logarithmic weighting is used for mimicing energy deposit in transverse direction + * + * Author: Sylvester Joosten (ANL), Chao Peng (ANL) 09/19/2022 + */ + +#include +#include +#include + +// Data types +#include +#include +#include +#include + +namespace algorithms::calorimetry { + +using ClusterRecoCoGBase = Algorithm< + Input>, + Output>>; + +/** Clustering with center of gravity method. + * + * Reconstruct the cluster with Center of Gravity method + * Logarithmic weighting is used for mimicking energy deposit in transverse direction + * + * \ingroup reco + */ +class ClusterRecoCoG : public ClusterRecoCoGBase { +public: + using Input = ClusterRecoCoGBase::Input; + using Output = ClusterRecoCoGBase::Output; + using WeightFunc = std::function; + + // TODO: get rid of "Collection" in names + ClusterRecoCoG(std::string_view name) + : ClusterRecoCoGBase{name, + {"inputProtoClusterCollection", "mcHits"}, + {"outputClusterCollection", "outputAssociations"}} {} + + void init(); + void process(const Input&, const Output&); + +private: + edm4eic::MutableCluster reconstruct(const edm4eic::ProtoCluster&) const; + + // FIXME do we really want sampling fraction here? + Property m_sampFrac{this, "samplingFraction", 1.0}; + Property m_logWeightBase{this, "logWeightBase", 3.6}; + Property m_energyWeight{this, "energyWeight", "log"}; + Property m_moduleDimZName{this, "moduleDimZName", ""}; + // Constrain the cluster position eta to be within + // the eta of the contributing hits. This is useful to avoid edge effects + // for endcaps. + Property m_enableEtaBounds{this, "enableEtaBounds", true}; + + WeightFunc m_weightFunc; + + const GeoSvc& m_geo = GeoSvc::instance(); +}; +} // namespace algorithms::calorimetry + diff --git a/external/algorithms/calorimetry/src/ClusterRecoCoG.cpp b/external/algorithms/calorimetry/src/ClusterRecoCoG.cpp index 88dd369..3feb247 100644 --- a/external/algorithms/calorimetry/src/ClusterRecoCoG.cpp +++ b/external/algorithms/calorimetry/src/ClusterRecoCoG.cpp @@ -7,350 +7,242 @@ * * Author: Chao Peng (ANL), 09/27/2020 */ + +#include + #include -#include #include #include -#include "boost/algorithm/string/join.hpp" -#include "fmt/format.h" -#include - -#include "Gaudi/Property.h" -#include "GaudiAlg/GaudiAlgorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/PhysicalConstants.h" -#include "GaudiKernel/RndmGenerators.h" -#include "GaudiKernel/ToolHandle.h" - -#include "DDRec/CellIDPositionConverter.h" -#include "DDRec/Surface.h" -#include "DDRec/SurfaceManager.h" - -#include "JugBase/DataHandle.h" -#include "JugBase/IGeoSvc.h" +#include +#include // Event Model related classes -#include "edm4hep/MCParticle.h" -#include "edm4hep/SimCalorimeterHitCollection.h" -#include "edm4eic/ClusterCollection.h" -#include "edm4eic/MCRecoClusterParticleAssociationCollection.h" -#include "edm4eic/ProtoClusterCollection.h" #include "edm4eic/vector_utils.h" -using namespace Gaudi::Units; +namespace algorithms::calorimetry { +namespace { -namespace Jug::Reco { - -// weighting functions (with place holders for hit energy, total energy, one parameter and module -// type enum -static double constWeight(double /*E*/, double /*tE*/, double /*p*/, int /*type*/) { return 1.0; } -static double linearWeight(double E, double /*tE*/, double /*p*/, int /*type*/) { return E; } -static double logWeight(double E, double tE, double base, int /*type*/) { - return std::max(0., base + std::log(E / tE)); -} - -static const std::map> weightMethods{ - {"none", constWeight}, - {"linear", linearWeight}, - {"log", logWeight}, -}; - -/** Clustering with center of gravity method. - * - * Reconstruct the cluster with Center of Gravity method - * Logarithmic weighting is used for mimicking energy deposit in transverse direction - * - * \ingroup reco - */ -class ClusterRecoCoG : public GaudiAlgorithm { -private: - Gaudi::Property m_sampFrac{this, "samplingFraction", 1.0}; - Gaudi::Property m_logWeightBase{this, "logWeightBase", 3.6}; - Gaudi::Property m_depthCorrection{this, "depthCorrection", 0.0}; - Gaudi::Property m_energyWeight{this, "energyWeight", "log"}; - Gaudi::Property m_moduleDimZName{this, "moduleDimZName", ""}; - // Constrain the cluster position eta to be within - // the eta of the contributing hits. This is useful to avoid edge effects - // for endcaps. - Gaudi::Property m_enableEtaBounds{this, "enableEtaBounds", false}; - - DataHandle m_inputProto{"inputProtoClusterCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputClusters{"outputClusterCollection", Gaudi::DataHandle::Writer, this}; - - // Collection for MC hits when running on MC - Gaudi::Property m_mcHits{this, "mcHits", ""}; - // Optional handle to MC hits - std::unique_ptr> m_mcHits_ptr; - - // Collection for associations when running on MC - Gaudi::Property m_outputAssociations{this, "outputAssociations", ""}; - // Optional handle to MC hits - std::unique_ptr> m_outputAssociations_ptr; - - // Pointer to the geometry service - SmartIF m_geoSvc; - double m_depthCorr{0}; - std::function weightFunc; - -public: - ClusterRecoCoG(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputProtoClusterCollection", m_inputProto, ""); - declareProperty("outputClusterCollection", m_outputClusters, ""); + // weighting functions (with place holders for hit energy, total energy, one parameter + double constWeight(double /*E*/, double /*tE*/, double /*p*/) { return 1.0; } + double linearWeight(double E, double /*tE*/, double /*p*/) { return E; } + double logWeight(double E, double tE, double base) { + return std::max(0., base + std::log(E / tE)); } - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } - - // Initialize the optional MC input hit collection if requested - if (m_mcHits != "") { - m_mcHits_ptr = - std::make_unique>(m_mcHits, Gaudi::DataHandle::Reader, - this); - } - - // Initialize the optional association collection if requested - if (m_outputAssociations != "") { - m_outputAssociations_ptr = - std::make_unique>(m_outputAssociations, Gaudi::DataHandle::Writer, - this); - } - - m_geoSvc = service("GeoSvc"); - if (!m_geoSvc) { - error() << "Unable to locate Geometry Service. " - << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg; - return StatusCode::FAILURE; - } - // update depth correction if a name is provided - if (!m_moduleDimZName.value().empty()) { - m_depthCorrection = m_geoSvc->detector()->constantAsDouble(m_moduleDimZName); - } - - // select weighting method - std::string ew = m_energyWeight.value(); - // make it case-insensitive - std::transform(ew.begin(), ew.end(), ew.begin(), [](char s) { return std::tolower(s); }); - auto it = weightMethods.find(ew); - if (it == weightMethods.end()) { - error() << fmt::format("Cannot find energy weighting method {}, choose one from [{}]", m_energyWeight, - boost::algorithm::join(weightMethods | boost::adaptors::map_keys, ", ")) - << endmsg; - return StatusCode::FAILURE; - } - weightFunc = it->second; - // info() << "z_length " << depth << endmsg; - return StatusCode::SUCCESS; + const std::map weightMethods{ + {"none", constWeight}, + {"linear", linearWeight}, + {"log", logWeight}, + }; +} // namespace + +void ClusterRecoCoG::init() { + // select weighting method + std::string ew = m_energyWeight; + // make it case-insensitive + std::transform(ew.begin(), ew.end(), ew.begin(), [](char s) { return std::tolower(s); }); + if (weightMethods.count(ew)) { + std::vector keys; + std::transform(weightMethods.begin(), weightMethods.end(), std::back_inserter(keys), + [](const auto& keyvalue) { return keyvalue.first; }); + raise(fmt::format("Cannot find energy weighting method {}, choose one from {}", m_energyWeight, + keys)); } + m_weightFunc = weightMethods.at(ew); + info() << fmt::format("Energy weight method set to: {}", ew) << endmsg; +} - StatusCode execute() override { - // input collections - const auto& proto = *m_inputProto.get(); - auto& clusters = *m_outputClusters.createAndPut(); - - // Optional input MC data - const edm4hep::SimCalorimeterHitCollection* mchits = nullptr; - if (m_mcHits_ptr) { - mchits = m_mcHits_ptr->get(); - } - - // Optional output associations - edm4eic::MCRecoClusterParticleAssociationCollection* associations = nullptr; - if (m_outputAssociations_ptr) { - associations = m_outputAssociations_ptr->createAndPut(); - } - - for (const auto& pcl : proto) { - auto cl = reconstruct(pcl); - - if (msgLevel(MSG::DEBUG)) { - debug() << cl.getNhits() << " hits: " << cl.getEnergy() / GeV << " GeV, (" << cl.getPosition().x / mm << ", " - << cl.getPosition().y / mm << ", " << cl.getPosition().z / mm << ")" << endmsg; - } - clusters.push_back(cl); - - // If mcHits are available, associate cluster with MCParticle - // 1. find proto-cluster hit with largest energy deposition - // 2. find first mchit with same CellID - // 3. assign mchit's MCParticle as cluster truth - if (m_mcHits_ptr.get() != nullptr && m_outputAssociations_ptr.get() != nullptr) { - - // 1. find pclhit with largest energy deposition - auto pclhits = pcl.getHits(); - auto pclhit = std::max_element( - pclhits.begin(), - pclhits.end(), - [](const auto& pclhit1, const auto& pclhit2) { - return pclhit1.getEnergy() < pclhit2.getEnergy(); - } - ); - - // 2. find mchit with same CellID - // find_if not working, https://github.com/AIDASoft/podio/pull/273 - //auto mchit = std::find_if( - // mchits.begin(), - // mchits.end(), - // [&pclhit](const auto& mchit1) { - // return mchit1.getCellID() == pclhit->getCellID(); - // } - //); - auto mchit = mchits->begin(); - for ( ; mchit != mchits->end(); ++mchit) { - // break loop when CellID match found - if (mchit->getCellID() == pclhit->getCellID()) { - break; - } - } - if (!(mchit != mchits->end())) { - // break if no matching hit found for this CellID - warning() << "Proto-cluster has highest energy in CellID " << pclhit->getCellID() - << ", but no mc hit with that CellID was found." << endmsg; - info() << "Proto-cluster hits: " << endmsg; - for (const auto& pclhit1: pclhits) { - info() << pclhit1.getCellID() << ": " << pclhit1.getEnergy() << endmsg; - } - info() << "MC hits: " << endmsg; - for (const auto& mchit1: *mchits) { - info() << mchit1.getCellID() << ": " << mchit1.getEnergy() << endmsg; - } +void ClusterRecoCoG::process(const ClusterRecoCoG::Input& input, + const ClusterRecoCoG::Output& output) { + const auto [proto, opt_simhits] = input; + auto [clusters, opt_assoc] = output; + + for (const auto& pcl : *proto) { + auto cl = reconstruct(pcl); + + if (aboveDebugThreshold()) { + debug() << cl.getNhits() << " hits: " << cl.getEnergy() / dd4hep::GeV << " GeV, (" + << cl.getPosition().x / dd4hep::mm << ", " << cl.getPosition().y / dd4hep::mm << ", " + << cl.getPosition().z / dd4hep::mm << ")" << endmsg; + } + clusters->push_back(cl); + + // If mcHits are available, associate cluster with MCParticle + // 1. find proto-cluster hit with largest energy deposition + // 2. find first mchit with same CellID + // 3. assign mchit's MCParticle as cluster truth + if (opt_simhits && opt_assoc) { + + // 1. find pclhit with largest energy deposition + auto pclhits = pcl.getHits(); + auto pclhit = std::max_element(pclhits.begin(), pclhits.end(), + [](const auto& pclhit1, const auto& pclhit2) { + return pclhit1.getEnergy() < pclhit2.getEnergy(); + }); + + // 2. find mchit with same CellID + // find_if not working, https://github.com/AIDASoft/podio/pull/273 + // auto mchit = std::find_if( + // opt_simhits->begin(), + // opt_simhits->end(), + // [&pclhit](const auto& mchit1) { + // return mchit1.getCellID() == pclhit->getCellID(); + // } + //); + auto mchit = opt_simhits->begin(); + for (; mchit != opt_simhits->end(); ++mchit) { + // break loop when CellID match found + if (mchit->getCellID() == pclhit->getCellID()) { break; } - - // 3. find mchit's MCParticle - const auto& mcp = mchit->getContributions(0).getParticle(); - - // debug output - if (msgLevel(MSG::DEBUG)) { - debug() << "cluster has largest energy in cellID: " << pclhit->getCellID() << endmsg; - debug() << "pcl hit with highest energy " << pclhit->getEnergy() << " at index " << pclhit->getObjectID().index << endmsg; - debug() << "corresponding mc hit energy " << mchit->getEnergy() << " at index " << mchit->getObjectID().index << endmsg; - debug() << "from MCParticle index " << mcp.getObjectID().index << ", PDG " << mcp.getPDG() << ", " << edm4eic::magnitude(mcp.getMomentum()) << endmsg; + } + if (!(mchit != opt_simhits->end())) { + // error condition should not happen + // break if no matching hit found for this CellID + warning() << "Proto-cluster has highest energy in CellID " << pclhit->getCellID() + << ", but no mc hit with that CellID was found." << endmsg; + info() << "Proto-cluster hits: " << endmsg; + for (const auto& pclhit1 : pclhits) { + info() << pclhit1.getCellID() << ": " << pclhit1.getEnergy() << endmsg; } - - // set association - edm4eic::MutableMCRecoClusterParticleAssociation clusterassoc; - clusterassoc.setRecID(cl.getObjectID().index); - clusterassoc.setSimID(mcp.getObjectID().index); - clusterassoc.setWeight(1.0); - clusterassoc.setRec(cl); - //clusterassoc.setSim(mcp); - associations->push_back(clusterassoc); - } else { - if (msgLevel(MSG::DEBUG)) { - debug() << "No mcHitCollection was provided, so no truth association will be performed." << endmsg; + info() << "MC hits: " << endmsg; + for (const auto& mchit1 : *opt_simhits) { + info() << mchit1.getCellID() << ": " << mchit1.getEnergy() << endmsg; } + break; } - } - return StatusCode::SUCCESS; - } - -private: - edm4eic::MutableCluster reconstruct(const edm4eic::ProtoCluster& pcl) const { - edm4eic::MutableCluster cl; - cl.setNhits(pcl.hits_size()); - - // no hits - if (msgLevel(MSG::DEBUG)) { - debug() << "hit size = " << pcl.hits_size() << endmsg; - } - if (pcl.hits_size() == 0) { - return cl; - } - - // 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]; - if (msgLevel(MSG::DEBUG)) { - debug() << "hit energy = " << hit.getEnergy() << " hit weight: " << weight << endmsg; - } - 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; + // 3. find mchit's MCParticle + const auto& mcp = mchit->getContributions(0).getParticle(); + + // debug output + if (aboveDebugThreshold()) { + debug() << "cluster has largest energy in cellID: " << pclhit->getCellID() << endmsg; + debug() << "pcl hit with highest energy " << pclhit->getEnergy() << " at index " + << pclhit->getObjectID().index << endmsg; + debug() << "corresponding mc hit energy " << mchit->getEnergy() << " at index " + << mchit->getObjectID().index << endmsg; + debug() << "from MCParticle index " << mcp.getObjectID().index << ", PDG " << mcp.getPDG() + << ", " << edm4eic::magnitude(mcp.getMomentum()) << endmsg; } - } - 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]; - float w = weightFunc(hit.getEnergy() * weight, totalE, m_logWeightBase.value(), 0); - tw += w; - v = v + (hit.getPosition() * w); - } - if (tw == 0.) { - warning() << "zero total weights encountered, you may want to adjust your weighting parameter." << endmsg; - } - cl.setPosition(v / tw); - cl.setPositionError({}); // @TODO: Covariance matrix - - // Optionally constrain the cluster to the hit eta values - if (m_enableEtaBounds) { - const bool overflow = (edm4eic::eta(cl.getPosition()) > 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)); - if (msgLevel(MSG::DEBUG)) { - debug() << "Bound cluster position to contributing hits due to " << (overflow ? "overflow" : "underflow") - << endmsg; - } + // set association + edm4eic::MutableMCRecoClusterParticleAssociation clusterassoc; + clusterassoc.setRecID(cl.getObjectID().index); + clusterassoc.setSimID(mcp.getObjectID().index); + clusterassoc.setWeight(1.0); + clusterassoc.setRec(cl); + // clusterassoc.setSim(mcp); + opt_assoc->push_back(clusterassoc); + } else { + if (aboveDebugThreshold()) { + debug() << "No mcHitCollection was provided, so no truth association will be performed." + << endmsg; } } + } +} - // Additional convenience variables +edm4eic::MutableCluster ClusterRecoCoG::reconstruct(const edm4eic::ProtoCluster& pcl) const { + edm4eic::MutableCluster cl; + cl.setNhits(pcl.hits_size()); - // 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 + // no hits + if (aboveDebugThreshold()) { + debug() << "hit size = " << pcl.hits_size() << endmsg; + } + if (pcl.hits_size() == 0) { + return cl; + } - // 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; + // 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]; + if (aboveDebugThreshold()) { + debug() << "hit energy = " << hit.getEnergy() << " hit weight: " << weight << endmsg; + } + 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]; + float w = m_weightFunc(hit.getEnergy() * weight, totalE, m_logWeightBase.value()); + tw += w; + v = v + (hit.getPosition() * w); + } + if (tw == 0.) { + warning() << "zero total weights encountered, you may want to adjust your weighting parameter." + << endmsg; + } + cl.setPosition(v / tw); + cl.setPositionError({}); // @TODO: Covariance matrix + + // Optionally constrain the cluster to the hit eta values + if (m_enableEtaBounds) { + const bool overflow = (edm4eic::eta(cl.getPosition()) > 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)); + if (aboveDebugThreshold()) { + debug() << "Bound cluster position to contributing hits due to " + << (overflow ? "overflow" : "underflow") << endmsg; } - radius = sqrt((1. / (cl.getNhits() - 1.)) * radius); - cl.addToShapeParameters(radius); - cl.addToShapeParameters(0 /* skewness */); // skewness not yet calculated } + } - return cl; + // 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 } -}; -// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(ClusterRecoCoG) + return cl; +} -} // namespace Jug::Reco +} // namespace algorithms::calorimetry diff --git a/external/algorithms/core/include/algorithms/algorithm.h b/external/algorithms/core/include/algorithms/algorithm.h index 5e80d51..92fec5b 100644 --- a/external/algorithms/core/include/algorithms/algorithm.h +++ b/external/algorithms/core/include/algorithms/algorithm.h @@ -15,26 +15,26 @@ namespace algorithms { // T should either be the desired input type, a std::vector<> of the desired input type, // or a std::optional<> of the desired input type template struct Input : std::tuple...> { + constexpr static const size_t kSize = sizeof...(T); using value_type = std::tuple...>; using data_type = std::tuple; - constexpr static const size_t kSize = sizeof...(T); + using index_type = std::array; }; template struct Output : std::tuple...> { + constexpr static const size_t kSize = sizeof...(T); using value_type = std::tuple...>; using data_type = std::tuple; - constexpr static const size_t kSize = sizeof...(T); + using index_type = std::array; }; -template using data_name = std::array; - // TODO: C++20 Concepts version for better error handling template class Algorithm : public PropertyMixin, public LoggerMixin { public: - using Input = InputType; - using Output = OutputType; - using InputNames = data_name; - using OutputNames = data_name; + using Input = typename InputType::value_type; + using Output = typename OutputType::value_type; + using InputNames = typename InputType::index_type; + using OutputNames = typename OutputType::index_type; Algorithm(std::string_view name, const InputNames& input_names, const OutputNames& output_names) : LoggerMixin(name), m_input_names{input_names}, m_output_names{output_names} {} diff --git a/external/algorithms/core/include/algorithms/property.h b/external/algorithms/core/include/algorithms/property.h index c3a3c73..e9b3437 100644 --- a/external/algorithms/core/include/algorithms/property.h +++ b/external/algorithms/core/include/algorithms/property.h @@ -78,8 +78,8 @@ class Configurable { } Property() = delete; - Property(const Property&) = delete; - void operator=(const Property&) = delete; + Property(const Property&) = default; + Property& operator=(const Property&) = default; // Only settable by explicitly calling the ::set() member functio n // as we want the Property to mostly act as if it is constant @@ -162,4 +162,8 @@ template std::ostream& operator<<(std::ostream& os, const algorithms::Configurable::Property& p) { return os << p.value(); } + +// Make Property formateble +template +struct fmt::formatter> : fmt::formatter {}; // others needed??? TODO