Implement Track-Cluster Subtraction (PFA1)

src/algorithms/calorimetry/TrackClusterSubtraction.cc

@@ -0,0 +1,50 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Derek Anderson
+
+#include <edm4eic/CalorimeterHit.h>
+#include <edm4hep/Vector3f.h>
+#include <edm4hep/utils/vector_utils.h>
+#include <fmt/core.h>
+#include <podio/ObjectID.h>
+#include <podio/RelationRange.h>
+#include <stdint.h>
+#include <cmath>
+#include <cstddef>
+#include <gsl/pointers>
+
+// algorithm definition
+#include "TrackClusterSubtraction.h"
+#include "algorithms/calorimetry/TrackClusterSubtractionConfig.h"
+
+
+
+namespace eicrecon {
+
+  // --------------------------------------------------------------------------
+  //! Initialize algorithm
+  // --------------------------------------------------------------------------
+  void TrackClusterSubtraction::init(const dd4hep::Detector* detector) {
+
+    // grab detector id
+    m_idCalo = detector -> constant<int>(m_cfg.idCalo);
+    debug("Collecting projections to detector with system id {}", m_idCalo);
+
+  }  // end 'init(dd4hep::Detector*)'
+
+
+
+  // --------------------------------------------------------------------------
+  //! Process inputs
+  // --------------------------------------------------------------------------
+  /*! TODO fill in
+   */
+  void TrackClusterSubtraction::process(
+    const TrackClusterSubtraction::Input& input,
+    const TrackClusterSubtraction::Output& output
+  ) const {
+
+    /* TODO fill in */
+
+  }  // end 'get_projections(edm4eic::CalorimeterHit&, edm4eic::TrackSegmentCollection&, VecTrkPoint&)'
+
+}  // end eicrecon namespace

src/algorithms/calorimetry/TrackClusterSubtraction.h

@@ -0,0 +1,119 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Derek Anderson
+
+#pragma once
+
+#include <DD4hep/Detector.h>
+#include <algorithms/algorithm.h>
+#include <edm4eic/ClusterCollection.h>
+#include <edm4eic/TrackPoint.h>
+#include <edm4eic/TrackSegmentCollection.h>
+#include <edm4hep/Vector3f.h>
+#include <podio/ObjectID.h>
+#include <algorithm>
+#include <map>
+#include <set>
+#include <string>
+#include <string_view>
+#include <vector>
+
+// for algorithm configuration
+#include "TrackClusterSubtractionConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+
+
+namespace eicrecon {
+
+  // --------------------------------------------------------------------------
+  //! Comparator struct for clusters
+  // --------------------------------------------------------------------------
+  /*! Organizes protoclusters by their ObjectID's in decreasing collection
+   *  ID first, and second by decreasing index second.
+   *
+   *  TODO order by energy first,then by objectID
+   */
+  struct CompareCluster {
+
+    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
+  // --------------------------------------------------------------------------
+  using VecProj = std::vector<edm4eic::TrackPoint>;
+  using VecClust = std::vector<edm4eic::Cluster>;
+  using SetClust = std::set<edm4eic::Cluster, CompareProto>;
+  using MapToVecProj = std::map<edm4eic::Cluster, VecProj, CompareProto>;
+  using MapToVecClust = std::map<edm4eic::Cluster, VecClust, CompareProto>;
+
+
+
+  // --------------------------------------------------------------------------
+  //! Algorithm input/output
+  // --------------------------------------------------------------------------
+  using TrackClusterSubtractionAlgorithm = algorithms::Algorithm<
+    algorithms::Input<
+      edm4eic::ClusterCollection,
+      edm4eic::ClusterCollection,
+      edm4eic::TrackSegmentCollection
+    >,
+    algorithms::Output<
+      edm4eic::ClusterCollection,
+      edm4eic::ClusterCollection
+    >
+  >;
+
+
+
+  // --------------------------------------------------------------------------
+  //! Track-Cluster Subtraction
+  // --------------------------------------------------------------------------
+  /*! An algorithm which takes collections of EM and HCal clusters, matches
+   *  track projections, subtracts the sum of the energy the projections
+   *  to the clusters, and returns the remnants of the subtracted clusters.
+   */
+  class TrackClusterSubtraction :
+    public TrackClusterSubtractionAlgorithm,
+    public WithPodConfig<TrackClusterSubtractionConfig>
+  {
+
+    public:
+
+      // ctor
+      TrackClusterSubtraction(std::string_view name) :
+        TrackClusterSubtractionAlgorithm {
+          name,
+          {
+            "InputEMCalClusterCollection",
+            "InputHCalClusterCollection",
+            "InputTrackProjections"
+          },
+          {"OutputClusterCollection"},
+          "Subtracts energy of tracks pointing to clusters."
+        } {}
+
+      // public methods
+      void init(const dd4hep::Detector* detector);
+      void process (const Input&, const Output&) const final;
+
+    private:
+
+      // private methods
+      /* TODO fill in */
+
+      // calorimeter id
+      int m_idCalo {0};
+
+  };  // end TrackClusterSubtraction
+
+}  // end eicrecon namespace

src/algorithms/calorimetry/TrackClusterSubtractionConfig.h

@@ -0,0 +1,21 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Derek Anderson
+
+#pragma once
+
+#include <string>
+
+namespace eicrecon {
+
+  struct TrackClusterSubtractionConfig {
+
+    std::string idCalo = "HcalBarrel_ID";  // id of calorimeter to match projections to
+
+    /* TODO parameters will go here */
+
+    // scale for hit-track distance
+    double transverseEnergyProfileScale = 1.0;
+
+  };  // end TrackClusterSubtractionConfig
+
+}  // end eicrecon namespace

src/factories/calorimetry/TrackClusterSubtraction_factory.h

@@ -0,0 +1,66 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Derek Anderson
+
+#pragma once
+
+// c++ utilities
+#include <string>
+// dd4hep utilities
+#include <DD4hep/Detector.h>
+// eicrecon components
+#include "extensions/jana/JOmniFactory.h"
+#include "services/geometry/dd4hep/DD4hep_service.h"
+#include "services/algorithms_init/AlgorithmsInit_service.h"
+#include "algorithms/calorimetry/TrackClusterSubtraction.h"
+
+namespace eicrecon {
+
+  class TrackClusterSubtraction_factory : public JOmniFactory<TrackClusterSubtraction_factory, TrackClusterSubtractionConfig> {
+
+    public:
+
+      using AlgoT = eicrecon::TrackClusterSubtraction;
+
+    private:
+
+      // algorithm to run
+      std::unique_ptr<AlgoT> m_algo;
+
+      // input collections
+      PodioInput<edm4eic::Cluster> m_emclusters_input {this};
+      PodioInput<edm4eic::Cluster> m_hclusters_input {this};
+      PodioInput<edm4eic::TrackSegment> m_track_projections_input {this};
+
+      // output collections
+      PodioOutput<edm4eic::Cluster> m_emclusters_output {this};
+      PodioOutput<edm4eic::Cluster> m_hclusters_output {this};
+
+      // parameter bindings
+      /* TODO fill in */
+
+      // services
+      Service<DD4hep_service> m_geoSvc {this};
+      Service<AlgorithmsInit_service> m_algoInitSvc {this};
+
+    public:
+
+      void Configure() {
+        m_algo = std::make_unique<AlgoT>(GetPrefix());
+        m_algo->applyConfig( config() );
+        m_algo->init(m_geoSvc().detector());
+      }
+
+      void ChangeRun(int64_t run_number) {
+        /* nothing to do here */
+      }
+
+      void Process(int64_t run_number, uint64_t event_number) {
+        m_algo->process(
+          {m_emclusters_input(), m_hclusters_input(), m_track_projections_input()},
+          {m_emclusters_output().get(), m_hclusters_output().get()}
+        );
+      }
+
+  };  // end TrackClusterSubtraction_factory
+
+}  // end eicrecon namespace