PF hit association

DQM/include/DQM/EcalClusterAnalyzer.h

@@ -64,7 +64,13 @@ class EcalClusterAnalyzer : public framework::Analyzer {
   // Pass Name for clusters
   std::string cluster_pass_name_;
 
-  std::string ecal_sp_hits_passname_;
+  // Collection Name for Scoring Plane hits
+  std::string ecal_sp_hits_coll_name_;
+  // Pass Name for Scoring Plane hits
+  std::string ecal_sp_hits_pass_name_;
+
+  // min energy fraction from smaller contributor to consider hit "mixed"
+  double mixed_hit_cutoff_;
 };
 
 }  // namespace dqm

DQM/python/dqm.py

@@ -1060,7 +1060,9 @@ def __init__(self,name='EcalClusterAnalyzer') :
         self.cluster_coll_name = 'EcalClusters'
         self.cluster_pass_name = ''
 
-        self.ecal_sp_hits_passname = ''
+        self.ecal_sp_hits_coll_name = 'EcalScoringPlaneHits'
+        self.ecal_sp_hits_pass_name = ''
+        self.mixed_hit_cutoff = 0.05
 
 
         self.build1DHistogram("number_of_clusters_first_layer", "Number of CLUE clusters on the first layer", 5, -0.5, 4.5)
@@ -1074,20 +1076,27 @@ def __init__(self,name='EcalClusterAnalyzer') :
         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", "Number of RecHits", 20, 0., 500.)
-
+        self.build1DHistogram("unclustered_hits", "Number of hits not in a cluster", 10, 0., 200.)
+        self.build1DHistogram("unclustered_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.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("total_energy_vs_hits", "Total energy (edep) [MeV]", 30, 0., 150., "Hits in cluster", 30, 0., 300.)
+        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.)
 
+        self.build1DHistogram("sp_distance", "dR(SPhit_1, SPhit_2)", 100, -1, 202)
+        self.build1DHistogram("cluster_distance", "dR(cl_1, cl_2)", 100, -1, 202)
+        self.build1DHistogram("cluster_RMSX", "RMS(hits in cluster) X", 100, -1, 202)
+        self.build1DHistogram("cluster_RMSY", "RMS(hits in cluster) Y", 100, -1, 202)
+        self.build2DHistogram("dE_cl2_vs_cl1", "E_{cl}-E_{true}^{SP}, cluster 1 [MeV]", 100, -10000, 10000, "E_{cl}-E_{true}^{SP}, cluster 2 [MeV]", 100, -10000, 10000)
+
+        self.build2DHistogram("tag0frac_vs_SPdist", "dR(SPhit_1, SPhit_2)", 251, -1, 250,  "Fraction of mixed (purity less than {int(100*(1.-self.mixed_hit_cutoff))}%) ancestors", 200, 0, 1)
+
 ecal_dqm = [
         EcalDigiVerify(),
         EcalShowerFeatures(),

DQM/src/DQM/EcalClusterAnalyzer.cxx

@@ -16,7 +16,12 @@ void EcalClusterAnalyzer::configure(framework::config::Parameters& ps) {
   cluster_coll_name_ = ps.get<std::string>("cluster_coll_name");
   cluster_pass_name_ = ps.get<std::string>("cluster_pass_name");
 
-  ecal_sp_hits_passname_ = ps.get<std::string>("ecal_sp_hits_passname");
+  ecal_sp_hits_coll_name_ =
+      ps.getParameter<std::string>("ecal_sp_hits_coll_name");
+  ecal_sp_hits_pass_name_ =
+      ps.getParameter<std::string>("ecal_sp_hits_pass_name");
+  mixed_hit_cutoff_ = ps.getParameter<double>("mixed_hit_cutoff");
+
   return;
 }
 
@@ -74,7 +79,7 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) {
   // Determine the truth information for the recoil electron
   std::vector<std::vector<float>> sp_electron_positions;
   const auto& ecal_sp_hits{event.getCollection<ldmx::SimTrackerHit>(
-      "EcalScoringPlaneHits", ecal_sp_hits_passname_)};
+      "EcalScoringPlaneHits", ecal_sp_hits_pass_name_)};
 
   std::vector<ldmx::SimTrackerHit> sorted_sp_hits = ecal_sp_hits;
   std::sort(sorted_sp_hits.begin(), sorted_sp_hits.end(),
@@ -109,6 +114,8 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) {
                  << ", TS counted electrons: " << nbr_of_electrons
                  << ", SP electrons: " << sp_electron_positions.size();
 
+  std::map<int, float> true_energy;
+  std::map<int, float> delta_energy;
   double sp_ele_dist{9999.};
   if (nbr_of_electrons == 2 && sp_electron_positions.size() > 1) {
     // Measures sp_ele_distance between two electrons in the ECal scoring plane
@@ -120,11 +127,18 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) {
     sp_ele_dist = std::sqrt((pos1[0] - pos2[0]) * (pos1[0] - pos2[0]) +
                             (pos1[1] - pos2[1]) * (pos1[1] - pos2[1]));
 
+    histograms_.fill("sp_distance", sp_ele_dist);
+
   }  // end block about the scoring plane hits
 
   ldmx_log(trace) << "Distance between the two e- in the ECal scoring plane: "
                   << sp_ele_dist << " mm";
 
+  double tot_event_energy = 0;
+  std::vector<double> tot_origin_edep;
+  tot_origin_edep.resize(nbr_of_electrons_ + 1);
+  int n_mixed = 0;
+
   // Loop over the rechits and find the matching simhits
   ldmx_log(trace) << "Loop over the rechits and find the matching simhits";
   for (const auto& hit : ecal_rec_hits) {
@@ -133,31 +147,65 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) {
         [&hit](const auto& sim_hit) { return sim_hit.getID() == hit.getID(); });
     if (it != ecal_sim_hits.end()) {
       // if found a simhit matching this rechit
+      ldmx_log(trace) << "\tFound simhit matching rechit with ID"
+                      << hit.getID();
       int ancestor = 0;
       int prev_ancestor = 0;
       bool tagged = false;
       int tag = 0;
       std::vector<double> edep;
-      edep.resize(nbr_of_electrons + 1);
+      edep.resize(nbr_of_electrons_ + 1);
+      double eTot = 0;  // keep track of total from all counted ancestors
+      ldmx_log(trace) << "\t\tIt has " << it->getNumberOfContribs()
+                      << " contribs. ";
       for (int i = 0; i < it->getNumberOfContribs(); i++) {
         // for each contrib in this simhit
         const auto& contrib = it->getContrib(i);
         // get origin electron ID
         ancestor = contrib.origin_id_;
+        ldmx_log(trace) << "\t\t\tAncestor ID " << ancestor << " with edep "
+                        << contrib.edep_;
+        tot_event_energy += contrib.edep_;
         // store energy from this contrib at index = origin electron ID
-        if (ancestor <= nbr_of_electrons) edep[ancestor] += contrib.edep_;
+        if (ancestor <= nbr_of_electrons) {
+          edep[ancestor] += contrib.edep_;
+          tot_origin_edep[ancestor] += contrib.edep_;
+          eTot += contrib.edep_;
+        }
         if (!tagged && i != 0 && prev_ancestor != 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;
+          ldmx_log(trace) << "\t\t\t\tMixed hit! Ancestor ID changed to "
+                          << ancestor;
         }
         prev_ancestor = ancestor;
+      }  // over contribs
+      // now check if mixed really means mixed, i.e. more than small fraction
+      // from a second electron.
+      if (tagged) {
+        for (int i = 1; i < nbr_of_electrons_ + 1; i++) {
+          if (edep[i] / eTot >
+              1 - mixed_hit_cutoff_) {  // one ancestor contributes at least the
+                                        // complement to the allowed mixing
+                                        // fraction
+            tagged = false;
+            ancestor =
+                distance(edep.begin(), max_element(edep.begin(), edep.end()));
+            ldmx_log(trace)
+                << "\t\t\t\tUndid mixed hit tagging, now ancestor = "
+                << ancestor;
+            break;
+          }
+        }
       }
       if (!tagged) {
-        // if not tagged, hit was from a single electron
-        tag = prev_ancestor;
-      }
+        // if not tagged, hit was from a single electron (within acceptable
+        // purity)
+        tag = ancestor;  // prev_ancestor;
+      } else
+        n_mixed++;
       histograms_.fill("ancestors", tag);
       hit_info.insert({hit.getID(), std::make_pair(tag, edep)});
     }  // end if simhit found
@@ -166,6 +214,22 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) {
   // Loop over the clusters
   int clustered_hits = 0;
   ldmx_log(trace) << "Loop over the clusters, N = " << n_ecal_clusters;
+  histograms_.fill("tag0frac_vs_SPdist", sp_ele_dist,
+                   (float)n_mixed / ecal_rec_hits.size());
+  ldmx_log(debug) << "Got " << n_mixed << " mixed hits, a fraction of "
+                  << (float)n_mixed / ecal_rec_hits.size();
+
+  if (ecal_clusters.size() >= 2) {
+    float dR = std::sqrt(std::pow((ecal_clusters[0].getCentroidX() -
+                                   ecal_clusters[1].getCentroidX()),
+                                  2) +
+                         std::pow((ecal_clusters[0].getCentroidY() -
+                                   ecal_clusters[1].getCentroidY()),
+                                  2));
+    histograms_.fill("cluster_distance", dR);
+    ldmx_log(trace) << "Gt cluster distance (0,1) = " << dR;
+  }
+
   for (const auto& cl : ecal_clusters) {
     auto layer = cl.getLayer();
     ldmx_log(trace) << "Cluster in layer " << layer
@@ -212,7 +276,7 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) {
     energy_from_electron.resize(nbr_of_electrons + 1);
     double energy_sum = 0.;
     double n_sum = 0.;
-
+    ldmx_log(trace) << "Looping over hits in the cluster";
     const auto& hit_ids = cl.getHitIDs();
     for (const auto& id : hit_ids) {
       // for each hit in cluster, find previously stored info
@@ -250,6 +314,11 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) {
       // get largest energy contribution
       double max_energy_contribution = *max_element(
           energy_from_electron.begin(), energy_from_electron.end());
+      std::string to_log;
+      for (auto nb : energy_from_electron)
+        to_log.append(std::to_string(nb) + " ");
+      ldmx_log(debug) << "Energies vector is " << to_log;
+
       // energy purity = largest contribution / all energy
       histograms_.fill("energy_percentage",
                        100. * (max_energy_contribution / energy_sum));
@@ -272,12 +341,25 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) {
                                   n_hits_from_electron.end());
       histograms_.fill("same_ancestor", 100. * (n_max / n_sum));
     }
+    // find the main contributor
+    auto elt = distance(
+        energy_from_electron.begin(),
+        max_element(energy_from_electron.begin(), energy_from_electron.end()));
+    ldmx_log(debug) << "Found that the maximum contributing trackID is " << elt;
+    delta_energy[elt] = cl.getEnergy() - true_energy[elt];
+    //    delta_energy[2]=e[2]-true_energy[2];
+    histograms_.fill("cluster_RMSX", cl.getRMSX());
   }  // end loop on the clusters
+  std::string more_log;
+  for (auto nb : tot_origin_edep) more_log.append(std::to_string(nb) + " ");
+  ldmx_log(debug) << "Edep per ancestor in event is " << more_log;
+  ldmx_log(debug) << "Total energy deposited in event: " << tot_event_energy;
 
-  histograms_.fill("clusterless_hits", (ecal_rec_hits.size() - clustered_hits));
+  histograms_.fill("dE_cl2_vs_cl1", delta_energy[1], delta_energy[2]);
+  histograms_.fill("unclustered_hits", (ecal_rec_hits.size() - clustered_hits));
   histograms_.fill("total_rechits_in_event", ecal_rec_hits.size());
   histograms_.fill(
-      "clusterless_hits_percentage",
+      "unclustered_hits_percentage",
       100. * (ecal_rec_hits.size() - clustered_hits) / ecal_rec_hits.size());
 }
 

Ecal/include/Ecal/CLUE.h

@@ -79,6 +79,9 @@ class CLUE {
   std::vector<std::shared_ptr<Density>> setup(
       const std::vector<const ldmx::EcalHit*>& hits);
 
+  // get distance between clusters in the first layers, proxy for electron sep.
+  void electronSeparation(std::vector<ldmx::EcalHit> hits);
+
   // connectingLayers marks if we're currently doing 3D clustering (i.e.
   // connecting seeds between layers) otherwise, layerTag tells us which layer
   // number we're working on