diff --git a/src/algorithms/CMakeLists.txt b/src/algorithms/CMakeLists.txt index 9adcb8fb49..ebebfd4699 100644 --- a/src/algorithms/CMakeLists.txt +++ b/src/algorithms/CMakeLists.txt @@ -3,6 +3,7 @@ cmake_minimum_required(VERSION 3.16) add_subdirectory(interfaces) add_subdirectory(calorimetry) add_subdirectory(tracking) +add_subdirectory(particle) add_subdirectory(pid) add_subdirectory(pid_lut) add_subdirectory(digi) diff --git a/src/algorithms/particle/CMakeLists.txt b/src/algorithms/particle/CMakeLists.txt new file mode 100644 index 0000000000..9075ac6dac --- /dev/null +++ b/src/algorithms/particle/CMakeLists.txt @@ -0,0 +1,18 @@ +cmake_minimum_required(VERSION 3.16) + +set(PLUGIN_NAME "algorithms_particle") + +# Function creates ${PLUGIN_NAME}_plugin and ${PLUGIN_NAME}_library targets +# Setting default includes, libraries and installation paths +plugin_add(${PLUGIN_NAME} WITH_SHARED_LIBRARY WITHOUT_PLUGIN) + +# The macro grabs sources as *.cc *.cpp *.c and headers as *.h *.hh *.hpp Then +# correctly sets sources for ${_name}_plugin and ${_name}_library targets Adds +# headers to the correct installation directory +plugin_glob_all(${PLUGIN_NAME}) + +# Find dependencies +plugin_add_algorithms(${PLUGIN_NAME}) +plugin_add_dd4hep(${PLUGIN_NAME}) +plugin_add_event_model(${PLUGIN_NAME}) +plugin_add_eigen3(${PLUGIN_NAME}) diff --git a/src/algorithms/particle/PFTools.h b/src/algorithms/particle/PFTools.h new file mode 100644 index 0000000000..f46cd7d47e --- /dev/null +++ b/src/algorithms/particle/PFTools.h @@ -0,0 +1,63 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2025 Derek Anderson + +#pragma once + +#include +#include +#include +#include +#include +#include + + + +namespace eicrecon { + + // -------------------------------------------------------------------------- + //! Particle Flow tools namespace + // -------------------------------------------------------------------------- + /*! This namespace collects a variety of useful types and methods + * used throughout particle flow algorithms. + */ + namespace PFTools { + + // ------------------------------------------------------------------------ + //! Comparator struct for clusters + // ------------------------------------------------------------------------ + /*! Organizes protoclusters by their ObjectID's in decreasing collection + * ID first, and second by decreasing index second. + * + * TODO should also order by energy... + */ + struct CompareClust { + + bool operator() (const edm4eic::Cluster& lhs, const edm4eic::Cluster& rhs) const { + if (lhs.getObjectID().collectionID == rhs.getObjectID().collectionID) { + return (lhs.getObjectID().index < rhs.getObjectID().index); + } else { + return (lhs.getObjectID().collectionID < rhs.getObjectID().collectionID); + } + } + + }; // end CompareCluster + + + + // ------------------------------------------------------------------------ + //! Convenience types + // ------------------------------------------------------------------------ + typedef std::vector> MatrixF; + typedef std::vector VecMatrixF; + typedef std::vector VecTrk; + typedef std::vector VecProj; + typedef std::vector VecSeg; + typedef std::vector VecClust; + typedef std::set SetClust; + typedef std::map MapToVecTrk; + typedef std::map MapToVecSeg; + typedef std::map MapToVecProj; + typedef std::map MapToVecClust; + + } // end PFTools namespace +} // end eicrecon namespace diff --git a/src/algorithms/particle/TrackClusterSubtractor.cc b/src/algorithms/particle/TrackClusterSubtractor.cc new file mode 100644 index 0000000000..fc0e9f7853 --- /dev/null +++ b/src/algorithms/particle/TrackClusterSubtractor.cc @@ -0,0 +1,224 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Derek Anderson + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "TrackClusterSubtractor.h" +#include "algorithms/particle/TrackClusterSubtractorConfig.h" + + + +namespace eicrecon { + + // -------------------------------------------------------------------------- + //! Initialize algorithm + // -------------------------------------------------------------------------- + void TrackClusterSubtractor::init() { + + //... nothing to do ...// + + } // end 'init(dd4hep::Detector*)' + + + + // -------------------------------------------------------------------------- + //! Process inputs + // -------------------------------------------------------------------------- + /*! Subtract energy of matched tracks via the following algorithm. + * 1. Build a map of each cluster onto a list of matched + * track projections. + * 2. For each cluster, subtract the sum of momenta of + * all matched tracks scaled by the specified fraction + * from the cluster's energy. + * 3. If subtracted energy is greater than 0, copy cluster + * into output with new subtracted energy. + */ + void TrackClusterSubtractor::process( + const TrackClusterSubtractor::Input& input, + const TrackClusterSubtractor::Output& output + ) const { + + // grab inputs/outputs + const auto [in_matches, in_projections] = input; + auto [out_sub_clusters, out_remain_clusters, out_sub_matches] = output; + + // exit if no matched tracks in collection + if (in_matches->size() == 0) { + debug("No matched tracks in collection."); + return; + } + + // ------------------------------------------------------------------------ + // 1. Build map of clusters onto projections + // ------------------------------------------------------------------------ + PFTools::MapToVecSeg mapClustToProj; + for (const auto& match : *in_matches) { + for (const auto& project : *in_projections) { + + // pick out corresponding projection from track + if (match.getTrack() != project.getTrack()) { + continue; + } else { + mapClustToProj[ match.getCluster() ].push_back( project ); + } + + } // end projection loop + } // end track-cluster match loop + debug("Built map of clusters-onto-tracks, size = {}", mapClustToProj.size()); + + // ------------------------------------------------------------------------ + // 2. Subtract energy for tracks + // ------------------------------------------------------------------------ + for (const auto& [cluster, projects] : mapClustToProj) { + + // do subtraction + const double eTrkSum = sum_track_energy(projects); + const double eToSub = m_cfg.fracEnergyToSub * eTrkSum; + const double eSub = cluster.getEnergy() - eToSub; + trace("Subtracted {} GeV from cluster with {} GeV", eToSub, cluster.getEnergy()); + + // check if consistent with zero, + // set eSub accordingly + const bool isZero = is_zero(eSub); + const double eSubToUse = isZero ? 0. : eSub; + + // calculate energy fractions + const double remainFrac = eSubToUse / cluster.getEnergy(); + const double subtractFrac = 1. - remainFrac; + + // scale subtracted cluster energy + edm4eic::MutableCluster sub_clust = cluster.clone(); + sub_clust.setEnergy( subtractFrac * cluster.getEnergy() ); + out_sub_clusters->push_back(sub_clust); + trace( + "Created subtracted cluster with {} GeV (originally {} GeV)", + sub_clust.getEnergy(), + cluster.getEnergy() + ); + + // create track cluster matches + for (const auto& project : projects) { + edm4eic::MutableTrackClusterMatch match = out_sub_matches->create(); + match.setCluster( sub_clust ); + match.setTrack( project.getTrack() ); + match.setWeight( 1.0 ); // FIXME placeholder + trace( + "Matched subtracted cluster {} to track {}", + sub_clust.getObjectID().index, + project.getTrack().getObjectID().index + ); + } + + // if NOT consistent with zero, write + // out remnant cluster + if (!isZero) { + edm4eic::MutableCluster remain_clust = cluster.clone(); + remain_clust.setEnergy( remainFrac * cluster.getEnergy() ); + out_remain_clusters->push_back(remain_clust); + trace( + "Created remnant cluster with {} GeV", + remain_clust.getEnergy() + ); + } + + } // end cluster-to-projections loop + debug("Finished subtraction, {} remnant clusters", out_remain_clusters->size()); + + } // end 'get_projections(edm4eic::CalorimeterHit&, edm4eic::TrackSegmentCollection&, VecTrkPoint&)' + + + + // -------------------------------------------------------------------------- + //! Sum energy of tracks + // -------------------------------------------------------------------------- + /*! Sums energy of tracks projected to the surface in the + * calorimeter specified by `surfaceToUse`. Uses PDG of + * track to select mass for energy; if not available, + * uses mass set by `defaultMassPdg`. + */ + double TrackClusterSubtractor::sum_track_energy(const PFTools::VecSeg& projects) const { + + double eSum = 0.; + for (const auto& project : projects) { + + // measure momentum at specified surface + double momentum = 0.; + for (const auto& point : project.getPoints()) { + if (point.surface != m_cfg.surfaceToUse) { + continue; + } else { + momentum = edm4hep::utils::magnitude( point.momentum ); + break; + } + } + + // get mass based on track pdg + double mass = m_parSvc.particle(m_cfg.defaultMassPdg).mass; + if (project.getTrack().getPdg() != 0) { + mass = m_parSvc.particle(project.getTrack().getPdg()).mass; + } + + // increment sum + eSum += std::sqrt((momentum*momentum) + (mass*mass)); + } + + // output debugging and exit + trace("Sum of track energy = {} GeV", eSum); + return eSum; + + } // end 'sum_track_energy(PFTools::VecSeg&)' + + + + // -------------------------------------------------------------------------- + //! Is difference consistent with zero? + // -------------------------------------------------------------------------- + /*! Checks if provided difference is consistent with zero, + * either checking if difference is within an epsilon + * (if `doNSigmaCut` is false), or if difference is within + * `nSigmaMax` of zero (if `doNSigmaCut` is true) based on + * the provided tracker and calorimeter resolutions. + */ + bool TrackClusterSubtractor::is_zero(const double difference) const { + + // if < 0, automatically return true + if (difference < 0) { + return true; + } + + // calculate nSigma + const double totReso = std::sqrt((m_cfg.trkReso*m_cfg.trkReso) + (m_cfg.calReso*m_cfg.calReso)); + const double nSigma = difference / totReso; + + // do appropriate comparison + bool isZero = false; + if (m_cfg.doNSigmaCut) { + isZero = (nSigma < m_cfg.nSigmaMax); + trace( + "Difference of {} GeV consistent with zero: nSigma = {} < {}", + difference, + nSigma, + m_cfg.nSigmaMax + ); + } else { + isZero = std::abs(difference) < std::numeric_limits::epsilon(); + trace( + "Difference of {} GeV consistent with zero within an epsilon", + difference + ); + } + return isZero; + + } // end 'is_zero(double)' + +} // end eicrecon namespace diff --git a/src/algorithms/particle/TrackClusterSubtractor.h b/src/algorithms/particle/TrackClusterSubtractor.h new file mode 100644 index 0000000000..210c470df0 --- /dev/null +++ b/src/algorithms/particle/TrackClusterSubtractor.h @@ -0,0 +1,86 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2025 Derek Anderson + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "PFTools.h" +#include "TrackClusterSubtractorConfig.h" +#include "algorithms/interfaces/ParticleSvc.h" +#include "algorithms/interfaces/WithPodConfig.h" + + +namespace eicrecon { + + // -------------------------------------------------------------------------- + //! Algorithm input/output + // -------------------------------------------------------------------------- + using TrackClusterSubtractorAlgorithm = algorithms::Algorithm< + algorithms::Input< + edm4eic::TrackClusterMatchCollection, + edm4eic::TrackSegmentCollection + >, + algorithms::Output< + edm4eic::ClusterCollection, + edm4eic::ClusterCollection, + edm4eic::TrackClusterMatchCollection + > + >; + + + + // -------------------------------------------------------------------------- + //! Track-Cluster Subtraction + // -------------------------------------------------------------------------- + /*! An algorithm which takes a collection of clusters and their matched + * tracks, subtracts the sum of all tracks pointing to the cluster, + * and outputs the remnant cluster and their matched tracks. + */ + class TrackClusterSubtractor + : public TrackClusterSubtractorAlgorithm + , public WithPodConfig + { + + public: + + // ctor + TrackClusterSubtractor(std::string_view name) : + TrackClusterSubtractorAlgorithm { + name, + {"inputTrackClusterMatches", "inputTrackProjections"}, + {"outputSubtractedClusterCollection", + "outputRemnantClusterCollection", + "outputTrackSubtractedClusterMatches"}, + "Subtracts energy of tracks pointing to clusters." + } {} + + // public methods + void init(); + void process(const Input&, const Output&) const final; + + private: + + // private methods + double sum_track_energy(const PFTools::VecSeg& projects) const; + bool is_zero(const double difference) const; + + // services + const algorithms::ParticleSvc& m_parSvc = algorithms::ParticleSvc::instance(); + + }; // end TrackClusterSubtractor + +} // end eicrecon namespace diff --git a/src/algorithms/particle/TrackClusterSubtractorConfig.h b/src/algorithms/particle/TrackClusterSubtractorConfig.h new file mode 100644 index 0000000000..96d2092534 --- /dev/null +++ b/src/algorithms/particle/TrackClusterSubtractorConfig.h @@ -0,0 +1,26 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2025 Derek Anderson + +#pragma once + +#include + +namespace eicrecon { + + struct TrackClusterSubtractorConfig { + + // general parameters + double fracEnergyToSub = 1.0; ///< fraction of track energy to subtract + int32_t defaultMassPdg = 211; ///< default mass to use for track energy + uint64_t surfaceToUse = 1; ///< index of surface to use for measuring momentum + + // parameters for resolution-based + // comparison + bool doNSigmaCut = false; ///< turn on/off checking against resolutions + uint32_t nSigmaMax = 1; ///< max no. of sigma to be consistent w/ zero + double trkReso = 1.0; ///< tracking momentum resolution to use + double calReso = 1.0; ///< calorimeter energy resolution to use + + }; // end TrackClusterSubtractorConfig + +} // end eicrecon namespace diff --git a/src/factories/CMakeLists.txt b/src/factories/CMakeLists.txt index c4934b6e50..13e01c1a22 100644 --- a/src/factories/CMakeLists.txt +++ b/src/factories/CMakeLists.txt @@ -1,5 +1,6 @@ add_subdirectory(calorimetry) add_subdirectory(digi) add_subdirectory(fardetectors) +add_subdirectory(particle) add_subdirectory(tracking) add_subdirectory(meta) diff --git a/src/factories/particle/CMakeLists.txt b/src/factories/particle/CMakeLists.txt new file mode 100644 index 0000000000..591839f6df --- /dev/null +++ b/src/factories/particle/CMakeLists.txt @@ -0,0 +1,2 @@ +set(PLUGIN_NAME "factories_particle") +plugin_headers_only(${PLUGIN_NAME}) diff --git a/src/factories/particle/TrackClusterSubtractor_factory.h b/src/factories/particle/TrackClusterSubtractor_factory.h new file mode 100644 index 0000000000..1369718efe --- /dev/null +++ b/src/factories/particle/TrackClusterSubtractor_factory.h @@ -0,0 +1,75 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Derek Anderson + +#pragma once + +// c++ utilities +#include +// dd4hep utilities +#include +// eicrecon components +#include "extensions/jana/JOmniFactory.h" +#include "services/geometry/dd4hep/DD4hep_service.h" +#include "services/algorithms_init/AlgorithmsInit_service.h" +#include "algorithms/particle/TrackClusterSubtractor.h" + +namespace eicrecon { + + class TrackClusterSubtractor_factory + : public JOmniFactory + { + + public: + + using AlgoT = eicrecon::TrackClusterSubtractor; + + private: + + // algorithm to run + std::unique_ptr m_algo; + + // input collections + PodioInput m_track_cluster_match_input {this}; + PodioInput m_track_projections_input {this}; + + // output collections + PodioOutput m_subtract_clusters_output {this}; + PodioOutput m_remnant_clusters_output {this}; + PodioOutput m_track_sub_cluster_match_output {this}; + + // parameter bindings + ParameterRef m_fracEnergyToSub {this, "fracEnergyToSub", config().fracEnergyToSub}; + ParameterRef m_defaultMassPdg {this, "defaultMassPdg", config().defaultMassPdg}; + ParameterRef m_surfaceToUse {this, "surfaceToUse", config().surfaceToUse}; + ParameterRef m_doNSigmaCut {this, "doNSigmaCut", config().doNSigmaCut}; + ParameterRef m_nSigmaMax {this, "nSigmaMax", config().nSigmaMax}; + ParameterRef m_trkReso {this, "trkReso", config().trkReso}; + ParameterRef m_calReso {this, "calReso", config().calReso}; + + // services + Service m_algoInitSvc {this}; + + public: + + void Configure() { + m_algo = std::make_unique(GetPrefix()); + m_algo->applyConfig( config() ); + m_algo->init(); + } + + void ChangeRun(int64_t run_number) { + //... nothing to do ...// + } + + void Process(int64_t run_number, uint64_t event_number) { + m_algo->process( + {m_track_cluster_match_input(), m_track_projections_input()}, + {m_subtract_clusters_output().get(), + m_remnant_clusters_output().get(), + m_track_sub_cluster_match_output().get()} + ); + } + + }; // end TrackClusterSubtractor_factory + +} // end eicrecon namespace diff --git a/src/global/particle/CMakeLists.txt b/src/global/particle/CMakeLists.txt new file mode 100644 index 0000000000..e5c5509550 --- /dev/null +++ b/src/global/particle/CMakeLists.txt @@ -0,0 +1,23 @@ +# Automatically set plugin name the same as the directory name Don't forget +# string(REPLACE " " "_" PLUGIN_NAME ${PLUGIN_NAME}) if this dir has spaces in +# its name +get_filename_component(PLUGIN_NAME ${CMAKE_CURRENT_LIST_DIR} NAME) + +# Function creates ${PLUGIN_NAME}_plugin and ${PLUGIN_NAME}_library targets +# Setting default includes, libraries and installation paths +plugin_add(${PLUGIN_NAME}) + +# The macro grabs sources as *.cc *.cpp *.c and headers as *.h *.hh *.hpp Then +# correctly sets sources for ${_name}_plugin and ${_name}_library targets Adds +# headers to the correct installation directory +plugin_glob_all(${PLUGIN_NAME}) + +# Find dependencies +plugin_add_event_model(${PLUGIN_NAME}) +plugin_add_dd4hep(${PLUGIN_NAME}) + +# Add include directories (works same as target_include_directories) +# plugin_include_directories(${PLUGIN_NAME} SYSTEM PUBLIC ...) + +# Add libraries (works same as target_include_directories) +plugin_link_libraries(${PLUGIN_NAME} algorithms_particle_library) diff --git a/src/global/particle/particle.cc b/src/global/particle/particle.cc new file mode 100644 index 0000000000..32100ca16e --- /dev/null +++ b/src/global/particle/particle.cc @@ -0,0 +1,145 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Derek Anderson + +#include +#include +#include +#include + +#include "extensions/jana/JOmniFactoryGeneratorT.h" +#include "factories/particle/TrackClusterMergeSplitter_factory.h" +#include "factories/particle/TrackClusterSubtractor_factory.h" + +extern "C" { + + void InitPlugin(JApplication *app) { + + using namespace eicrecon; + + InitJANAPlugin(app); + + // ==================================================================== + // PFAlpha: baseline PF implementation + // ==================================================================== + + // -------------------------------------------------------------------- + // PFA (0) connection: split/merge clusters accordingly + // -------------------------------------------------------------------- + + /* TODO move here when ready */ + + // -------------------------------------------------------------------- + // PFA (1a) arbitration: apply track correction to clusters + // -------------------------------------------------------------------- + + // backward ----------------------------------------------------------- + + app->Add( + new JOmniFactoryGeneratorT( + "EcalEndcapNSubtractedClusters", + {"EcalEndcapNTrackSplitMergeClusterMatches", + "CalorimeterTrackProjections"}, + {"EcalEndcapNSubtractedClusters", + "EcalEndcapNRemnantClusters", + "EcalEndcapNTrackSubtractedClusterMatches"}, + { + .fracEnergyToSub = 1.0, + .defaultMassPdg = 211, + .surfaceToUse = 1, + }, + app // TODO: remove me once fixed + ) + ); + + app->Add( + new JOmniFactoryGeneratorT( + "HcalEndcapNSubtractedClusters", + {"HcalEndcapNTrackSplitMergeClusterMatches", + "CalorimeterTrackProjections"}, + {"HcalEndcapNSubtractedClusters", + "HcalEndcapNRemnantClusters", + "HcalEndcapNTrackSubtractedClusterMatches"}, + { + .fracEnergyToSub = 1.0, + .defaultMassPdg = 211, + .surfaceToUse = 1 + }, + app // TODO: remove me once fixed + ) + ); + + // central ------------------------------------------------------------ + + app->Add( + new JOmniFactoryGeneratorT( + "HcalBarrelSubtractedClusters", + {"HcalBarrelTrackSplitMergeClusterMatches", + "CalorimeterTrackProjections"}, + {"HcalBarrelSubtractedClusters", + "HcalBarrelRemnantClusters", + "HcalBarrelTrackSubtractedClusterMatches"}, + { + .fracEnergyToSub = 1.0, + .defaultMassPdg = 211, + .surfaceToUse = 1 + }, + app // TODO: remove me once fixed + ) + ); + + // forward ------------------------------------------------------------ + + app->Add( + new JOmniFactoryGeneratorT( + "EcalEndcapPSubtractedClusters", + {"EcalEndcapPTrackSplitMergeClusterMatches", + "CalorimeterTrackProjections"}, + {"EcalEndcapPSubtractedClusters", + "EcalEndcapPRemnantClusters", + "EcalEndcapPTrackSubtractedClusterMatches"}, + { + .fracEnergyToSub = 1.0, + .defaultMassPdg = 211, + .surfaceToUse = 1 + }, + app // TODO: remove me once fixed + ) + ); + + app->Add( + new JOmniFactoryGeneratorT( + "LFHCALSubtractedClusters", + {"LFHCALTrackSplitMergeClusterMatches", + "CalorimeterTrackProjections"}, + {"LFHCALSubtractedClusters", + "LFHCALRemnantClusters", + "LFHCALTrackSubtractedClusterMatches"}, + { + .fracEnergyToSub = 1.0, + .defaultMassPdg = 211, + .surfaceToUse = 1 + }, + app // TODO: remove me once fixed + ) + ); + + // -------------------------------------------------------------------- + // PFA (1b) arbitration: form charged candidates + // -------------------------------------------------------------------- + + /* TODO add here */ + + // -------------------------------------------------------------------- + // PFA (2) arbitration: combine remnants, form neutral candidates + // -------------------------------------------------------------------- + + /* TODO add here */ + + // -------------------------------------------------------------------- + // PFA (3) regression: convert candidates to reco particles + // -------------------------------------------------------------------- + + /* TODO add here */ + + } +} diff --git a/src/utilities/eicrecon/eicrecon.cc b/src/utilities/eicrecon/eicrecon.cc index 798b30c640..7f732b5573 100644 --- a/src/utilities/eicrecon/eicrecon.cc +++ b/src/utilities/eicrecon/eicrecon.cc @@ -24,6 +24,7 @@ std::vector EICRECON_DEFAULT_PLUGINS = { "beam", "reco", "tracking", + "particle", "pid", "EEMC", "BEMC",