CLUE clustering improvements, fixing multilayer option (#1835)

Author: tvami

ComparePlots/ecal.py

@@ -121,7 +121,9 @@ def clue_cluster(d : Differ, out_dir = None) :
     """
 
     features = [
-        ('EcalClusterAnalyzer/EcalClusterAnalyzer_number_of_clusters', 'Number of CLUE clusters'),
+        ('EcalClusterAnalyzer/EcalClusterAnalyzer_number_of_clusters_first_layer', 'Number of CLUE clusters on the first layer'),
+        ('EcalClusterAnalyzer/EcalClusterAnalyzer_number_of_clusters_per_layer', 'Number of CLUE clusters per layer'),
+        ('EcalClusterAnalyzer/EcalClusterAnalyzer_number_of_clusters', 'Total number of CLUE clusters'),
         ('EcalClusterAnalyzer/EcalClusterAnalyzer_energy_percentage', 'Percentage of energy in cluster'),
         ('EcalClusterAnalyzer/EcalClusterAnalyzer_clusterless_hits_percentage', 'Percentage of hits not in a cluster'),
         ('EcalClusterAnalyzer/EcalClusterAnalyzer_sp_clue_distance', 'CLUE centroid to SP ele distance in xy-plane [mm]'),

DQM/python/dqm.py

@@ -1063,7 +1063,9 @@ def __init__(self,name='EcalClusterAnalyzer') :
         self.ecal_sp_hits_passname = ''
 
 
-        self.build1DHistogram("number_of_clusters", "Number of CLUE clusters", 5, -0.5, 4.5)
+        self.build1DHistogram("number_of_clusters_first_layer", "Number of CLUE clusters on the first layer", 5, -0.5, 4.5)
+        self.build1DHistogram("number_of_clusters_per_layer", "Number of CLUE clusters per layer", 5, -0.5, 4.5)
+        self.build1DHistogram("number_of_clusters", "Total number of CLUE clusters", 51, -0.5, 50.5)
         self.build1DHistogram("correctly_predicted_events", "Correct Cluster Count", 3, 0., 3.)
 
         # Need to mod for more than two electrons
@@ -1080,7 +1082,8 @@ def __init__(self,name='EcalClusterAnalyzer') :
         self.build2DHistogram("total_energy_vs_hits", "Total energy (edep) [MeV]", 30, 0, 150, "Hits in cluster", 20, 0, 200)
         self.build2DHistogram("total_energy_vs_purity", "Total energy (edep) [MeV]", 30, 0, 150, "Energy purity %", 21, 0, 105)
         self.build2DHistogram("sp_ele_distance_vs_purity", "SP ele distance in xy-plane [mm]", 50, 0, 250, "Energy purity %", 21, 0, 105)
-        
+        self.build2DHistogram("sp_clue_distance_vs_layer", "CLUE centroid to SP ele distance in xy-plane [mm]", 125, 0., 250., "Layer", 33, -0.5, 32.5)
+
         self.build1DHistogram("sp_clue_distance", "CLUE centroid to SP ele distance in xy-plane [mm]", 125, 0., 250.)
         self.build1DHistogram("sp_clue_x_residual", "CLUE centroid X - SP ele X [mm]", 250, -250., 250.)
         self.build1DHistogram("sp_clue_y_residual", "CLUE centroid Y - SP ele Y [mm]", 250, -250., 250.)

DQM/src/DQM/EcalClusterAnalyzer.cxx

@@ -36,9 +36,26 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) {
     nbr_of_electrons = nbr_of_electrons_;
   }
 
-  int n_ecal_clusters = ecal_clusters.size();
+  std::map<int, int> layer_cluster_count;
+  for (const auto& cluster : ecal_clusters) {
+    auto layer = cluster.getLayer();
+    layer_cluster_count[layer]++;
+  }
+
+  int total_clusters = 0;
+  for (const auto& [layer, count] : layer_cluster_count) {
+    total_clusters += count;
+  }
+
+  int n_ecal_clusters = static_cast<int>(std::round(
+      static_cast<double>(total_clusters) / layer_cluster_count.size()));
+
+  ldmx_log(info) << "Avg number of clusters per layer: " << n_ecal_clusters;
   // Fill histograms with the number of clusters
-  histograms_.fill("number_of_clusters", n_ecal_clusters);
+  histograms_.fill("number_of_clusters", total_clusters);
+  histograms_.fill("number_of_clusters_per_layer", n_ecal_clusters);
+  histograms_.fill("number_of_clusters_first_layer", layer_cluster_count[0]);
+
   // Fill simplied 3-bin histogram to check the prediction
   if (n_ecal_clusters == nbr_of_electrons) {
     // correct
@@ -150,8 +167,14 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) {
   int clustered_hits = 0;
   ldmx_log(trace) << "Loop over the clusters, N = " << n_ecal_clusters;
   for (const auto& cl : ecal_clusters) {
-    auto cluster_centroid_x = cl.getFirstLayerCentroidX();
-    auto cluster_centroid_y = cl.getFirstLayerCentroidY();
+    auto layer = cl.getLayer();
+    ldmx_log(trace) << "Cluster in layer " << layer
+                    << ", energy: " << cl.getEnergy()
+                    << ", number of hits: " << cl.getHitIDs().size();
+    auto cluster_centroid_x = cl.getCentroidX();
+    auto cluster_centroid_y = cl.getCentroidY();
+    auto cluster_rms_x = cl.getRMSX();
+    auto cluster_rms_y = cl.getRMSY();
 
     // Find the closest sp_electron_positions to the cluster centroid
     double min_distance = 9999.;
@@ -168,13 +191,17 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) {
       }
     }  // end loop on the scoring plane electron positions
     // Fill histogram with the distance to the closest scoring plane electron
-    ldmx_log(trace) << "\tCluster centroid: (" << cluster_centroid_x << ", "
-                    << cluster_centroid_y
-                    << ") mm, min distance to SP electron: " << min_distance
+    ldmx_log(trace) << "\tCluster centroid: (" << cluster_centroid_x << " +/- "
+                    << cluster_rms_x << ", " << cluster_centroid_y << " +/- "
+                    << cluster_rms_y
+                    << " mm; min distance to SP electron: " << min_distance
                     << " mm";
-    histograms_.fill("sp_clue_distance", min_distance);
-    histograms_.fill("sp_clue_x_residual", sp_clue_x_residuals);
-    histograms_.fill("sp_clue_y_residual", sp_clue_y_residuals);
+    if (layer == 0) {
+      histograms_.fill("sp_clue_distance", min_distance);
+      histograms_.fill("sp_clue_x_residual", sp_clue_x_residuals);
+      histograms_.fill("sp_clue_y_residual", sp_clue_y_residuals);
+    }
+    histograms_.fill("sp_clue_distance_vs_layer", layer, min_distance);
 
     // for each cluster
     // total number of hits coming from electron, index = electron ID

Ecal/include/Ecal/CLUE.h

@@ -29,23 +29,24 @@ class CLUE {
 
  public:
   struct Density {
-    float x_;
-    float y_;
-    float z_;
+    double x_;
+    double y_;
+    double z_;
     double total_energy_;
     int index_;
 
     // index of density this density is follower of
     // set to index of spatially closest density with higher energy; -1 if seed
     int follower_of_;
-    // separation distance to density that this is follower of
-    float delta_;
-    float z_delta_;
-
-    int cluster_id_;  // 2D cluster ID
-
+    // 2D separation distance to density that this is follower of
+    double delta_;
+    // separation in z to density that this is follower of
+    double z_delta_;
+    // 2D cluster ID
+    int cluster_id_;
+    // layer of density
     int layer_;
-
+    // hits in this density
     std::vector<const ldmx::EcalHit*> hits_;
 
     Density() {}
@@ -63,9 +64,15 @@ class CLUE {
     }
   };
 
-  float dist(double x1, double y1, double x2, double y2);
-  float floatDist(float x1, float y1, float x2, float y2);
-  float floatDist(float x1, float y1, float z1, float x2, float y2, float z2);
+  /**
+   * Euclidean distance between two points
+   * @tparam T floating point type
+   */
+  template <typename T>
+  T dist(T x1, T y1, T x2, T y2);
+  // 3D version
+  template <typename T>
+  T dist(T x1, T y1, T z1, T x2, T y2, T z2);
   std::vector<std::vector<const ldmx::EcalHit*>> createLayers(
       const std::vector<const ldmx::EcalHit*>& hits);
   float roundToDecimal(float x, int num_decimal_precision_digits);
@@ -79,7 +86,7 @@ class CLUE {
       std::vector<std::shared_ptr<Density>>& densities, bool connectingLayers,
       int layerTag = 0);
 
-  std::vector<std::shared_ptr<Density>> layerSetup();
+  std::vector<std::shared_ptr<Density>> setupForClue3D();
 
   void convertToIntermediateClusters(
       std::vector<std::vector<const ldmx::EcalHit*>>& clusters);
@@ -117,15 +124,10 @@ class CLUE {
   double deltao_;
   double dm_;
 
-  // layers in Ecal; a bit unsure if this is the correct number
-  int max_layers_{17};
+  // layers in Ecal
+  int max_layers_{32};
   int nbr_of_layers_;
-  // air between layers (a guess, but it seems to work)
-  double air_{10.};
-  // thickness of ECal layers
-  std::vector<double> layer_thickness_ = {2.,   3.5,  5.3,  5.3, 5.3, 5.3,
-                                          5.3,  5.3,  5.3,  5.3, 5.3, 10.5,
-                                          10.5, 10.5, 10.5, 10.5};
+
   std::vector<double> layer_rho_c_;
   std::vector<double> layer_delta_c_;
   // containment radius for the different layers of the ECal
@@ -144,7 +146,6 @@ class CLUE {
   std::vector<double> centroid_distances_;
   IntermediateCluster event_centroid_;
 
-  float first_layer_max_z_;
   std::vector<IntermediateCluster> first_layer_centroids_;
 
   int seed_index_{0};

Ecal/include/Ecal/Event/EcalCluster.h

@@ -57,8 +57,6 @@ class EcalCluster : public ldmx::CaloCluster {
   }
 
  private:
-  // Could add further ECal-specific info here...
-
   std::vector<unsigned int> first_layer_hit_ids_;
 
   double first_layer_centroid_x_{0};

Ecal/include/Ecal/IntermediateCluster.h

@@ -25,7 +25,8 @@ class IntermediateCluster {
   std::vector<const ldmx::EcalHit*> getHits() const { return hits_; }
   bool empty() const { return hits_.empty(); }
   void clear() { hits_.clear(); }
-  int getLayer() const { return layer_; }
+  int getLayer() const;
+  void setLayer(int layer);
 
  private:
   int layer_;

Ecal/include/Ecal/MyClusterWeight.h

@@ -10,14 +10,14 @@ namespace ecal {
 
 class MyClusterWeight {
  public:
+  // returns weighting function, where smallest, weights will be combined first
   double operator()(const IntermediateCluster& a,
-                    const IntermediateCluster&
-                        b) {  // returns weighting function, where smallest
-                              // weights will be combined first
-
-    double rmol = 10.00;    // Moliere radius_ of detector, roughly. In mm
-    double dzchar = 100.0;  // Characteristic cluster longitudinal variable TO
-                            // BE DETERMINED! in mm
+                    const IntermediateCluster& b) {
+    // Moliere radius_ of detector, roughly. In mm
+    double rmol = 10.00;
+    // Characteristic cluster longitudinal variable TO
+    // BE DETERMINED! in mm
+    double dzchar = 100.0;
 
     double a_e = a.centroid().E();
     double a_x = a.centroid().Px();

Ecal/python/ecalClusters.py

@@ -39,19 +39,13 @@ def __init__(self,name='ecalClusters') :
         # Currently only used when nbrOfLayers > 1
         self.dc = 5.
         # Minimum seed energy/maximum outlier energy
+        # Not used when nbrOfLayers > 1
         self.rhoc = 550.
         # Minimum seed separation
+        # Not used when nbrOfLayers > 1
         self.deltac = 10.
         # Minimum outlier separation
         self.deltao = 40.
         # Recluster merged clusters or not
         # No reclustering leads to more undercounting, reclustering leads to more overcounting
-        self.reclustering = False
-
-        self.build1DHistogram("nLoops", "Number of loops for clustering", 50, 0, 400) # not applicable for CLUE
-        self.build1DHistogram("nClusters", "Number of clusters", 20, 0, 20)
-        self.build1DHistogram("nHits", "Hits per cluster", 20, 0, 300)
-        self.build1DHistogram("cluster_energy", "Energy [MeV] per cluster", 100, 0, 20000)
-        self.build2DHistogram("seed_weights", "Number of seeds", 20, 0, 100, "Minimum weight", 20, 0, 10) # not applicable for CLUE
-        self.build2DHistogram("recluster", "Initial number of clusters", 20, 0, 20, "Number of clusters after reclustering", 20, 0, 20) # not applicable for existing algo
-
+        self.reclustering = False