Updated CLUE clustering to pull RoC values from CSV file, if available.

Ecal/include/Ecal/CLUE.h

@@ -8,8 +8,10 @@
 #define ECAL_CLUE_H_
 
 #include <math.h>
+#include <stdlib.h>
 
 #include <algorithm>
+#include <fstream>
 #include <iomanip>
 #include <iostream>
 #include <limits>
@@ -88,16 +90,16 @@ class CLUE {
   // number we're working on
   std::vector<std::vector<const ldmx::EcalHit*>> clustering(
       std::vector<std::shared_ptr<Density>>& densities, bool connectingLayers,
-      int layerTag = 0);
+      int layerTag = 0, std::string roc_file_name = "");
 
   std::vector<std::shared_ptr<Density>> setupForClue3D();
 
   void convertToIntermediateClusters(
       std::vector<std::vector<const ldmx::EcalHit*>>& clusters);
 
   void cluster(const std::vector<ldmx::EcalHit>& hits, double dc, double rc,
-               double deltac, double deltao, int nbrOfLayers,
-               bool reclustering);
+               double deltac, double deltao, int nbrOfLayers, bool reclustering,
+               std::string roc_file_name);
 
   std::vector<double> getCentroidDistances() const {
     return centroid_distances_;

Ecal/include/Ecal/EcalClusterProducer.h

@@ -74,6 +74,7 @@ class EcalClusterProducer : public framework::Producer {
   bool clue_;
   int nbr_of_layers_;
   bool reclustering_;
+  std::string roc_file_name_;
 
   /** The name of the cluster algorithm used. */
   TString algo_name_;

Ecal/python/ecal_clusters.py

@@ -7,7 +7,7 @@
 """
 
 from LDMX.Framework import Processor, processor
-
+from .make_path import make_roc_path
 
 @processor("ecal::EcalClusterProducer", "Ecal")
 class EcalClusterProducer(Processor):
@@ -45,10 +45,12 @@ class EcalClusterProducer(Processor):
         not used when nbr_of_layers > 1
     deltao: float
         minimum outlier separation
-    recluster: bool
+    reclustering: bool
         recluster merged clusters or not
         No reclustering leads to more undercounting
         reclustering leads to more overcounting
+    roc_file: str
+        Used in CLUE to set the roc values, if available
     """
 
     rec_hit_coll_name: str = "EcalRecHits"
@@ -65,3 +67,4 @@ class EcalClusterProducer(Processor):
     deltac: float = 10.0
     deltao: float = 40.0
     reclustering: bool = False
+    roc_file: str = make_roc_path("RoC_v14_8gev")

Ecal/src/Ecal/CLUE.cxx

@@ -1,4 +1,3 @@
-
 #include "Ecal/CLUE.h"
 
 #include <cmath>
@@ -254,13 +253,52 @@ std::vector<std::shared_ptr<CLUE::Density>> CLUE::setup(
 // number we're working on
 std::vector<std::vector<const ldmx::EcalHit*>> CLUE::clustering(
     std::vector<std::shared_ptr<CLUE::Density>>& densities,
-    bool connecting_layers, int layer_index) {
+    bool connecting_layers, int layer_index, std::string roc_file_name) {
   ldmx_log(trace) << "--- CLUSTERING ---";
   ldmx_log(trace) << "Number of densities: " << densities.size()
                   << "; connecting_layers: " << connecting_layers
                   << "; layer_index: " << layer_index;
+
   if (!connecting_layers && nbr_of_layers_ > 1) {
     // if layerwise clustering override rhoc_ and deltac_ with per-layer values
+
+    // If possible, overwrite hard-coded roc values with values imported from
+    // the CSV file These should only need to be read once
+    static std::vector<double> radius_from_file;
+    if (radius_from_file.empty() && roc_file_name != "") {
+      ldmx_log(info) << "Attempting to use RoC values from CSV file.";
+      // File reading algorithm adapted from EcalVetoProcessor.cxx
+      if (!std::ifstream(roc_file_name).good()) {
+        EXCEPTION_RAISE("CLUE", "The specified RoC file '" + roc_file_name +
+                                    "' does not exist!");
+      } else {
+        std::ifstream rocfile(roc_file_name);
+        std::string line, value;
+
+        // Throw away the first (header) line in the file
+        std::getline(rocfile, line);
+        // In EcalVetoProcessor, the RoC values are a 2D array for many angle
+        // ranges
+        // As near as I (CJ) can tell, the RoC values in CLUE are only for the
+        // first angle range, 0<theta<10. So we only read the first line of the
+        // CSV file's values
+
+        std::getline(rocfile, line);
+        std::stringstream ss(line);
+        int values_read = 0;
+        while (std::getline(ss, value, ',')) {
+          values_read++;
+          if (values_read < 5)
+            continue;  // First few entries in each line of RoC file are not RoC
+                       // values
+          float f_value = (value != "") ? std::stof(value) : -1.0;
+          radius_from_file.push_back(f_value);
+        }
+      }
+    }
+
+    if (!radius_from_file.empty()) radius_ = radius_from_file;
+
     rhoc_ = layer_rho_c_[layer_index];
     ldmx_log(trace) << "Setting rho_c on layer " << layer_index << " to "
                     << rhoc_;
@@ -278,7 +316,6 @@ std::vector<std::vector<const ldmx::EcalHit*>> CLUE::clustering(
     // deltac_ = 100.;
     // rhoc_ = 1000.;
   }
-
   bool energy_overload = false;
   double max_energy = 10000.;
   clustering_loops_ = 0;
@@ -581,7 +618,7 @@ void CLUE::convertToIntermediateClusters(
 
 void CLUE::cluster(const std::vector<ldmx::EcalHit>& unsorted_hits, double dc,
                    double rc, double delta_c, double delta_o, int nbr_of_layers,
-                   bool reclustering) {
+                   bool reclustering, std::string roc_file_name) {
   ldmx_log(info) << "Starting CLUE clustering with parameters:" << "dc " << dc
                  << ", rc " << rc << ", delta_c " << delta_c << ", delta_o "
                  << delta_o << ", nbr_of_layers " << nbr_of_layers
@@ -636,7 +673,7 @@ void CLUE::cluster(const std::vector<ldmx::EcalHit>& unsorted_hits, double dc,
     for (int i = 0; i < layers.size(); i++) {
       ldmx_log(trace) << "--- LAYER " << i + 1 << " ---";
       auto densities = setup(layers[i]);
-      auto clusters = clustering(densities, false, i);
+      auto clusters = clustering(densities, false, i, roc_file_name);
       convertToIntermediateClusters(clusters);
       // clustering without 3D
     }
@@ -648,7 +685,7 @@ void CLUE::cluster(const std::vector<ldmx::EcalHit>& unsorted_hits, double dc,
   } else {
     ldmx_log(debug) << "Only one layer, doing 2D clustering";
     auto densities = setup(hits);
-    auto clusters = clustering(densities, false);
+    auto clusters = clustering(densities, false, 0, roc_file_name);
     convertToIntermediateClusters(clusters);
   }
 }

Ecal/src/Ecal/EcalClusterProducer.cxx

@@ -31,6 +31,7 @@ void EcalClusterProducer::configure(framework::config::Parameters& ps) {
   clue_ = ps.get<bool>("clue");
   nbr_of_layers_ = ps.get<int>("nbr_of_layers");
   reclustering_ = ps.get<bool>("reclustering");
+  roc_file_name_ = ps.get<std::string>("roc_file");
 }
 
 void EcalClusterProducer::produce(framework::Event& event) {
@@ -46,7 +47,7 @@ void EcalClusterProducer::produce(framework::Event& event) {
     ldmx_log(info) << "Using CLUE clustering algorithm";
     CLUE cf;
     cf.cluster(ecal_hits, dc_, rhoc_, deltac_, deltao_, nbr_of_layers_,
-               reclustering_);
+               reclustering_, roc_file_name_);
     ldmx_log(debug) << "CLUE algorithm finished clustering";
     std::vector<IntermediateCluster> interm_cluster = cf.getClusters();
     std::vector<IntermediateCluster> f_interm_cluster =