added densities to ECalCluster for debugging

LDMX-Software · Oct 16, 2024 · 005dba8 · 005dba8
1 parent b6f2c5f
commit 005dba8
Show file tree

Hide file tree

Showing 3 changed files with 78 additions and 2 deletions.
diff --git a/Ecal/include/Ecal/CLUE.h b/Ecal/include/Ecal/CLUE.h
@@ -4,7 +4,7 @@
  * @author Ella Viirola, Lund University
  */
 
-#ifndef ECAL_CLUE_H_
+#ifndef ECAL_CLUE_H_ 
 #define ECAL_CLUE_H_
 
 #include <math.h>
@@ -20,6 +20,7 @@
 
 #include "Ecal/Event/EcalHit.h"
 #include "Ecal/WorkingEcalCluster.h"
+#include "SimCore/Event/SimCalorimeterHit.h"
 
 namespace ecal {
 
@@ -44,6 +45,7 @@ class CLUE {
     int layer;
 
     std::vector<ldmx::EcalHit> hits;
+    std::vector<unsigned int> hit_origins;
 
     Density() {}
 
@@ -58,6 +60,45 @@ class CLUE {
       z = 0.;
       hits = {};
     }
+
+    void findHitOrigins(const std::vector<ldmx::SimCalorimeterHit>& ecalSimHits) {
+      std::vector<unsigned int> vecIDs;
+      for (const auto& hit : this->hits) {
+        auto id = hit.getID();
+        int tag = 0;
+        auto it = std::find_if(ecalSimHits.begin(), ecalSimHits.end(),
+          [&](const auto& simHit) { return simHit.getID() == id; });
+
+        if (it != ecalSimHits.end()) {
+          int ancestor = 0;
+          int prevAncestor = 0;
+          bool tagged = false;
+          tag = 0;
+          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;
+            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;
+          }
+        }
+        else {tag = -1;}
+        vecIDs.push_back(tag);
+      }
+      hit_origins = vecIDs; 
+    }
   };
 
   float dist(double x1, double y1, double x2, double y2) {
@@ -627,6 +668,13 @@ class CLUE {
       convertToWorkingClusters(clusters);
     } else {
       auto densities = setup(hits);
+      densities_ = densities;
+      //densities->findHitOrigins();
+      //for (auto& density: densities) {
+      //    densities_.push_back(density->totalEnergy);
+      //    auto count = std::count_if(density->hit_origins.begin(), density->hit_origins.end(),[&](auto const& val){ return val >= 10; });
+      //    if (count >= 1) {densities_secondaries_.push_back(density->totalEnergy);}
+      //}
       auto clusters = clustering(densities, false);
       convertToWorkingClusters(clusters);
     }
@@ -642,6 +690,10 @@ class CLUE {
 
   std::vector<WorkingEcalCluster> getClusters() const { return finalClusters_; }
 
+  std::vector<std::shared_ptr<Density>> getDensities() const { return densities_; }
+  std::vector<double> getDensitiesSecond() const { return densities_secondaries_; }
+
+
   // First layer centroids are available for potential future combination with
   // TS
   std::vector<WorkingEcalCluster> getFirstLayerCentroids() const {
@@ -691,6 +743,8 @@ class CLUE {
 
   int seedIndex_{0};
   std::vector<std::vector<std::shared_ptr<Density>>> seeds_;
+  std::vector<std::shared_ptr<Density>> densities_;
+  std::vector<double> densities_secondaries_;
 
   int initialClusterNbr_{-1};
   std::vector<WorkingEcalCluster> finalClusters_;

diff --git a/Ecal/include/Ecal/Event/EcalCluster.h b/Ecal/include/Ecal/Event/EcalCluster.h
@@ -44,6 +44,9 @@ class EcalCluster : public ldmx::CaloCluster {
 
   void findHitOrigins(const std::vector<ldmx::SimCalorimeterHit>& ecalSimHits);
 
+  void setDensities(const std::vector<double>& densities) {densities_ = densities;}
+  void setDensitiesSecond(const std::vector<double>& densities) {densities_secondaries_ = densities;}
+
   bool operator<(const EcalCluster& rhs) const {
     return this->getEnergy() < rhs.getEnergy();
   }
@@ -68,6 +71,8 @@ class EcalCluster : public ldmx::CaloCluster {
   std::vector<unsigned int> firstLayerHitIDs_;
   std::vector<unsigned int> hitOriginIDs_;
 
+  std::vector<double> densities_;
+  std::vector<double> densities_secondaries_;
 
   double firstLayerCentroidX_{0};
   double firstLayerCentroidY_{0};

diff --git a/Ecal/src/Ecal/EcalClusterProducer.cxx b/Ecal/src/Ecal/EcalClusterProducer.cxx
@@ -57,9 +57,20 @@ void EcalClusterProducer::produce(framework::Event& event) {
     CLUE cf;
 
     cf.cluster(ecalHits, dc_, rhoc_, deltac_, deltao_, nbrOfLayers_,
-               reclustering_, debug_);
+               reclustering_, debug_); 
+
     std::vector<WorkingEcalCluster> wcVec = cf.getClusters();
     std::vector<WorkingEcalCluster> fWcVec = cf.getFirstLayerCentroids();
+    std::vector<std::shared_ptr<CLUE::Density>> densities = cf.getDensities();
+    std::vector<double> densities_energy;
+    std::vector<double> densities_secondaries_energy;
+    for (auto& density: densities) {
+       density->findHitOrigins(ecalSimHits);
+       densities_energy.push_back(density->totalEnergy);
+       std::cout << density->totalEnergy << std::endl;
+       auto count = std::count_if(density->hit_origins.begin(), density->hit_origins.end(),[&](auto const& val){ return val >= 10; });
+       if (count >= 1) {densities_secondaries_energy.push_back(density->totalEnergy);}
+    }
 
     auto nLoops = cf.getNLoops();
     histograms_.fill("nLoops", nLoops);
@@ -83,6 +94,12 @@ void EcalClusterProducer::produce(framework::Event& event) {
       cluster.addFirstLayerHits(fWcVec[aWC].getHits());
       cluster.findHitOrigins(ecalSimHits);
 
+      if (aWC == 0) {
+        cluster.setDensities(densities_energy);
+        cluster.setDensitiesSecond(densities_secondaries_energy);
+      } 
+
+
       histograms_.fill("nHits", wcVec[aWC].getHits().size());
       histograms_.fill("cluster_energy", wcVec[aWC].centroid().E());