Add CLUE clustering (#1699)

* manually copy over CLUE clustering files from Ella's developments

  git switch ella-dev-clustering
  git rebase trunk
  git switch 1411-add-clue-clustering
  git checkout ella-dev-clustering -- Ecal/

Ella Viirola originally authored these files while doing summer research
into ECal clustering at Lund University

Co-authored-by: Lysarina

* store the track ID of the 'origin' particle which may be more precise than the incident

copied over from work done originally by Ella Viirola during Summer 2024
at Lund.

Co-authored-by: Lysarina

* bump ClassDef version after introducing new member

* use variable that was determined earlier

* move CLUE into logging, transition member variables to snake_case

* dont override, that interacts poorly with ROOT's ClassDef macro

* purge debug_ flag, drop debug messages to trace

* Apply clang-format

* parameter naming in python to snake_case

* add example config to test-run clustering

needed to patch python parameter naming and now it works

* merge Working clusters and prefer pointers to save on copying

* Apply clang-format

* remove unused variables

* rename WorkingCluster to IntermediateCluster for clarity

* move cluster analyzer into DQM

* move max origin track ID from hardcode to config parameter

add a bit more documentation to the SD implementation

* add ecal clustering to ecal_pn, inclusive, it_pileup samples

* clang-format

* fix typo in ecal_pn config

* update declaration macro from Factory updates

---------

Co-authored-by: github-actions[bot] <github-actions[bot]@users.noreply.github.com>

Author: tomeichlersmith

.github/validation_samples/ecal_pn/config.py

@@ -32,6 +32,7 @@
 import LDMX.Ecal.ecal_hardcoded_conditions
 import LDMX.Ecal.digi as ecal_digi
 import LDMX.Ecal.vetos as ecal_vetos
+import LDMX.Ecal.ecalClusters as ecal_cluster
 
 # Load the HCAL modules
 import LDMX.Hcal.HcalGeometry
@@ -83,6 +84,7 @@
 p.sequence.extend([
         ecal_digi.EcalDigiProducer(),
         ecal_digi.EcalRecProducer(), 
+        ecal_cluster.EcalClusterProducer(),
         ecal_veto,
         hcal_digi_reco,
         hcal_veto,
@@ -91,6 +93,7 @@
         trigScintTrack,
         count, TriggerProcessor('trigger', 8000.),
         dqm.PhotoNuclearDQM(),
+        dqm.EcalClusterAnalyzer()
         ])
 
 p.sequence.extend(dqm.all_dqm)

.github/validation_samples/inclusive/config.py

@@ -32,6 +32,7 @@
 import LDMX.Ecal.ecal_hardcoded_conditions
 import LDMX.Ecal.digi as ecal_digi
 import LDMX.Ecal.vetos as ecal_vetos
+import LDMX.Ecal.ecalClusters as ecal_cluster
 
 # Load the HCAL modules
 import LDMX.Hcal.HcalGeometry
@@ -83,6 +84,7 @@
 p.sequence.extend([
         ecal_digi.EcalDigiProducer(),
         ecal_digi.EcalRecProducer(), 
+        ecal_cluster.EcalClusterProducer(),
         ecalVeto,
         hcal_digi_reco,
         hcal_veto,
@@ -91,6 +93,7 @@
         trigScintTrack,
         count, TriggerProcessor('trigger', 8000.),
         dqm.PhotoNuclearDQM(),
+        dqm.EcalClusterAnalyzer()
         ])
 
 p.sequence.extend(dqm.all_dqm)

.github/validation_samples/it_pileup/config.py

@@ -48,6 +48,7 @@
 
 from LDMX.Ecal import digi as eDigi
 from LDMX.Ecal import vetos
+import LDMX.Ecal.ecalClusters as ecal_cluster
 from LDMX.Hcal import digi as hDigi
 
 # this is hardwired into the code to be appended to the sim hits collections
@@ -194,14 +195,15 @@
 p.sequence.extend(full_tracking_sequence.dqm_sequence)
 
 p.sequence.extend([
-    ecalDigi, ecalReco, ecalVeto,
+    ecalDigi, ecalReco, ecalVeto, ecal_cluster.EcalClusterProducer(),
     hcal_digi_reco, 
     hcal_veto,
     *ts_digis,
     *ts_clusters,
     trigScintTrack,
     count, TriggerProcessor('trigger', 8000.),
     dqm.PhotoNuclearDQM(),
+    dqm.EcalClusterAnalyzer()
 ])
 
 p.sequence.extend(dqm_with_overlay)

DQM/include/DQM/EcalClusterAnalyzer.h

@@ -0,0 +1,52 @@
+/**
+ * @file EcalClusterAnalyzer.h
+ * @brief Analysis of cluster performance
+ * @author Ella Viirola, Lund University
+ */
+
+#ifndef DQM_ECALCLUSTERANALYZER_H
+#define DQM_ECALCLUSTERANALYZER_H
+
+// LDMX Framework
+#include "Framework/Configure/Parameters.h"
+#include "Framework/EventProcessor.h"
+
+namespace dqm {
+
+/**
+ * @class EcalClusterAnalyzer
+ * @brief
+ */
+class EcalClusterAnalyzer : public framework::Analyzer {
+ public:
+  EcalClusterAnalyzer(const std::string& name, framework::Process& process)
+      : Analyzer(name, process) {}
+  ~EcalClusterAnalyzer() override = default;
+  void configure(framework::config::Parameters& ps) override;
+  void analyze(const framework::Event& event) override;
+
+ private:
+  int nbr_of_electrons_;
+
+  // Collection Name for SimHits
+  std::string ecal_sim_hit_coll_;
+
+  // Pass Name for SimHits
+  std::string ecal_sim_hit_pass_;
+
+  // Collection Name for RecHits
+  std::string rec_hit_coll_name_;
+
+  // Pass Name for RecHits
+  std::string rec_hit_pass_name_;
+
+  // Collection Name for clusters
+  std::string cluster_coll_name_;
+
+  // Pass Name for clusters
+  std::string cluster_pass_name_;
+};
+
+}  // namespace dqm
+
+#endif

DQM/python/dqm.py

@@ -804,6 +804,38 @@ def __init__(self, name='SampleValidation') :
         self.build1DHistogram("pdgid_hardbremdaughters", "ID of hard brem daughters", 20, 0, 20)
         self.build1DHistogram("startZ_hardbremdaughters", "Start z position of hard brem daughters  [mm]", 200, -1000, 1000)
 
+class EcalClusterAnalyzer(ldmxcfg.Analyzer) :
+    """Analyze clustering"""
+
+    def __init__(self,name='EcalClusterAnalyzer') :
+        super().__init__(name, "dqm::EcalClusterAnalyzer", 'DQM')
+
+        self.nbr_of_electrons = 2
+
+        self.ecal_sim_hit_coll = "EcalSimHits"
+        self.ecal_sim_hit_pass = "" #use whatever pass is available
+
+        # Pass name for ecal digis and rec hits
+        self.rec_hit_coll_name = 'EcalRecHits'
+        self.rec_hit_pass_name = ''
+
+        self.cluster_coll_name = 'ecalClusters'
+        self.cluster_pass_name = ''
+        
+        # Need to mod for more than two electrons
+        self.build1DHistogram("ancestors", "Ancestors of particles", 4, 0, 3)
+
+        self.build1DHistogram("same_ancestor", "Percentage of hits in cluster coming from the electron that produced most hits", 21, 0, 105)
+        self.build1DHistogram("energy_percentage", "Percentage of energy in cluster coming from the electron that produced most of energy", 21, 0, 105)
+        self.build1DHistogram("mixed_hit_energy", "Percentage of total energy coming from hits with energy contributions from more than one electron", 21, 0, 105)
+        self.build1DHistogram("clusterless_hits", "Number of hits not in a cluster", 10, 0, 200)
+        self.build1DHistogram("clusterless_hits_percentage", "Percentage of hits not in a cluster", 21, 0, 105)
+        self.build1DHistogram("total_rechits_in_event", "Rechits per event", 20, 0, 500)
+        self.build1DHistogram("correctly_predicted_events", "Correctly predicted events", 3, 0, 3)
+
+        self.build2DHistogram("total_energy_vs_hits", "Total energy (edep)", 30, 0, 150, "Hits in cluster", 20, 0, 200)
+        self.build2DHistogram("total_energy_vs_purity", "Total energy (edep)", 30, 0, 150, "Energy purity %", 21, 0, 105)
+        self.build2DHistogram("distance_energy_purity", "Distance in xy-plane", 20, 0, 220, "Energy purity %", 21, 0, 105)
 
 ecal_dqm = [
         EcalDigiVerify(),
@@ -855,4 +887,4 @@ def __init__(self, name='SampleValidation') :
         ]
 
 
-all_dqm = ecal_dqm + hcal_dqm + trigScint_dqm + trigger_dqm
+all_dqm = ecal_dqm + hcal_dqm + trigScint_dqm + trigger_dqm

DQM/src/DQM/EcalClusterAnalyzer.cxx

@@ -0,0 +1,185 @@
+/**
+ * @file EcalClusterAnalyzer.cxx
+ * @brief Analysis of cluster performance
+ * @author Ella Viirola, Lund University
+ */
+
+#include "DQM/EcalClusterAnalyzer.h"
+
+#include <algorithm>
+#include <fstream>
+#include <iostream>
+
+#include "DetDescr/SimSpecialID.h"
+#include "Ecal/Event/EcalCluster.h"
+#include "Ecal/Event/EcalHit.h"
+#include "SimCore/Event/SimCalorimeterHit.h"
+#include "SimCore/Event/SimTrackerHit.h"
+
+namespace dqm {
+
+void EcalClusterAnalyzer::configure(framework::config::Parameters& ps) {
+  nbr_of_electrons_ = ps.getParameter<int>("nbr_of_electrons");
+
+  ecal_sim_hit_coll_ = ps.getParameter<std::string>("ecal_sim_hit_coll");
+  ecal_sim_hit_pass_ = ps.getParameter<std::string>("ecal_sim_hit_pass");
+
+  rec_hit_coll_name_ = ps.getParameter<std::string>("rec_hit_coll_name");
+  rec_hit_pass_name_ = ps.getParameter<std::string>("rec_hit_pass_name");
+
+  cluster_coll_name_ = ps.getParameter<std::string>("cluster_coll_name");
+  cluster_pass_name_ = ps.getParameter<std::string>("cluster_pass_name");
+  return;
+}
+
+void EcalClusterAnalyzer::analyze(const framework::Event& event) {
+  const auto& ecal_rec_hits{event.getCollection<ldmx::EcalHit>(
+      rec_hit_coll_name_, rec_hit_pass_name_)};
+  const auto& ecal_sim_hits{event.getCollection<ldmx::SimCalorimeterHit>(
+      ecal_sim_hit_coll_, ecal_sim_hit_pass_)};
+  const auto& ecal_clusters{event.getCollection<ldmx::EcalCluster>(
+      cluster_coll_name_, cluster_pass_name_)};
+
+  if (ecal_clusters.size() == nbr_of_electrons_)
+    histograms_.fill("correctly_predicted_events", 1);  // correct
+  else if (ecal_clusters.size() < nbr_of_electrons_)
+    histograms_.fill("correctly_predicted_events", 0);  // undercounting
+  else if (ecal_clusters.size() > nbr_of_electrons_)
+    histograms_.fill("correctly_predicted_events", 2);  // overcounting
+
+  std::unordered_map<int, std::pair<int, std::vector<double>>> hitInfo;
+  hitInfo.reserve(ecal_rec_hits.size());
+
+  double dist;
+  if (nbr_of_electrons_ == 2) {
+    // Measures distance between two electrons in the ECal scoring plane
+    // TODO: generalize for n electrons
+    std::vector<float> pos1;
+    std::vector<float> pos2;
+    bool p1 = false;
+    bool p2 = false;
+
+    const auto& ecal_sp_hits{
+        event.getCollection<ldmx::SimTrackerHit>("EcalScoringPlaneHits")};
+    for (const ldmx::SimTrackerHit& spHit : ecal_sp_hits) {
+      if (spHit.getTrackID() == 1) {
+        pos1 = spHit.getPosition();
+        p1 = true;
+      } else if (spHit.getTrackID() == 2) {
+        pos2 = spHit.getPosition();
+        p2 = true;
+      }
+    }
+    if (p1 && p2)
+      dist = std::sqrt(std::pow((pos1[0] - pos2[0]), 2) +
+                       std::pow((pos1[1] - pos2[1]), 2));
+  }
+
+  for (const auto& hit : ecal_rec_hits) {
+    auto it = std::find_if(
+        ecal_sim_hits.begin(), ecal_sim_hits.end(),
+        [&hit](const auto& simHit) { return simHit.getID() == hit.getID(); });
+    if (it != ecal_sim_hits.end()) {
+      // if found a simhit matching this rechit
+      int ancestor = 0;
+      int prevAncestor = 0;
+      bool tagged = false;
+      int tag = 0;
+      std::vector<double> edep;
+      edep.resize(nbr_of_electrons_ + 1);
+      for (int i = 0; i < it->getNumberOfContribs(); i++) {
+        // for each contrib in this simhit
+        const auto& c = it->getContrib(i);
+        // get origin electron ID
+        ancestor = c.originID;
+        // store energy from this contrib at index = origin electron ID
+        if (ancestor <= nbr_of_electrons_) edep[ancestor] += c.edep;
+        if (!tagged && i != 0 && prevAncestor != ancestor) {
+          // if origin electron ID does not match previous origin electron ID
+          // this hit has contributions from several electrons, ie mixed case
+          tag = 0;
+          tagged = true;
+        }
+        prevAncestor = ancestor;
+      }
+      if (!tagged) {
+        // if not tagged, hit was from a single electron
+        tag = prevAncestor;
+      }
+      histograms_.fill("ancestors", tag);
+      hitInfo.insert({hit.getID(), std::make_pair(tag, edep)});
+    }
+  }
+
+  int clusteredHits = 0;
+
+  for (const auto& cl : ecal_clusters) {
+    // for each cluster
+    // total number of hits coming from electron, index = electron ID
+    std::vector<double> n;
+    n.resize(nbr_of_electrons_ + 1);
+    // total number of energy coming from electron, index = electron ID
+    std::vector<double> e;
+    e.resize(nbr_of_electrons_ + 1);
+    double eSum = 0.;
+    double nSum = 0.;
+
+    const auto& hitIDs = cl.getHitIDs();
+    for (const auto& id : hitIDs) {
+      // for each hit in cluster, find previously stored info
+      auto it = hitInfo.find(id);
+      if (it != hitInfo.end()) {
+        auto t = it->second;
+        auto eId = t.first;        // origin electron ID (or 0 for mixed)
+        auto energies = t.second;  // energy vector
+        n[eId]++;  // increment number of hits coming from this electron
+        nSum++;
+
+        double hitESum = 0.;
+        for (int i = 1; i < nbr_of_electrons_ + 1; i++) {
+          // loop through energy vector
+          if (energies[i] > 0.) {
+            hitESum += energies[i];
+            // add energy from electron i in this hit to total energy from
+            // electron i in cluster
+            e[i] += energies[i];
+          }
+        }
+        // if mixed hit, add the total energy of this hit to mixed hit energy
+        // counter
+        if (eId == 0) e[0] += hitESum;
+        eSum += hitESum;
+
+        clusteredHits++;
+      }
+    }
+
+    if (eSum > 0) {
+      // get largest energy contribution
+      double eMax = *max_element(e.begin(), e.end());
+      // energy purity = largest contribution / all energy
+      histograms_.fill("energy_percentage", 100. * (eMax / eSum));
+      if (e[0] > 0.) histograms_.fill("mixed_hit_energy", 100. * (e[0] / eSum));
+
+      histograms_.fill("total_energy_vs_hits", eSum, cl.getHitIDs().size());
+      histograms_.fill("total_energy_vs_purity", eSum, 100. * (eMax / eSum));
+
+      if (nbr_of_electrons_ == 2)
+        histograms_.fill("distance_energy_purity", dist, 100. * (eMax / eSum));
+    }
+    if (nSum > 0) {
+      double nMax = *max_element(n.begin(), n.end());
+      histograms_.fill("same_ancestor", 100. * (nMax / nSum));
+    }
+  }
+
+  histograms_.fill("clusterless_hits", (ecal_rec_hits.size() - clusteredHits));
+  histograms_.fill("total_rechits_in_event", ecal_rec_hits.size());
+  histograms_.fill(
+      "clusterless_hits_percentage",
+      100. * (ecal_rec_hits.size() - clusteredHits) / ecal_rec_hits.size());
+}
+
+}  // namespace dqm
+
+DECLARE_ANALYZER(dqm::EcalClusterAnalyzer)

Ecal/exampleConfigs/cluster.py

@@ -0,0 +1,13 @@
+from LDMX.Framework import ldmxcfg
+p = ldmxcfg.Process('cluster')
+import sys
+p.inputFiles = sys.argv[1:]
+p.outputFiles = [ 'clusters.root' ]
+p.histogramFile = 'h_clusters.root'
+
+from LDMX.Ecal.ecalClusters import *
+p.sequence = [
+    EcalClusterProducer(),
+    EcalClusterAnalyzer()
+]
+