diff --git a/.github/validation_samples/ecal_pn/config.py b/.github/validation_samples/ecal_pn/config.py index 5c29801a8..51c3ecb1a 100644 --- a/.github/validation_samples/ecal_pn/config.py +++ b/.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) diff --git a/.github/validation_samples/inclusive/config.py b/.github/validation_samples/inclusive/config.py index fafdcb4b3..7e8b9e8b6 100644 --- a/.github/validation_samples/inclusive/config.py +++ b/.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) diff --git a/.github/validation_samples/it_pileup/config.py b/.github/validation_samples/it_pileup/config.py index 7ce271dbf..12c761499 100644 --- a/.github/validation_samples/it_pileup/config.py +++ b/.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,7 +195,7 @@ 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, @@ -202,6 +203,7 @@ trigScintTrack, count, TriggerProcessor('trigger', 8000.), dqm.PhotoNuclearDQM(), + dqm.EcalClusterAnalyzer() ]) p.sequence.extend(dqm_with_overlay) diff --git a/DQM/include/DQM/EcalClusterAnalyzer.h b/DQM/include/DQM/EcalClusterAnalyzer.h new file mode 100644 index 000000000..a6647afa6 --- /dev/null +++ b/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 diff --git a/DQM/python/dqm.py b/DQM/python/dqm.py index a8da9bc02..f267abf46 100644 --- a/DQM/python/dqm.py +++ b/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 \ No newline at end of file +all_dqm = ecal_dqm + hcal_dqm + trigScint_dqm + trigger_dqm diff --git a/DQM/src/DQM/EcalClusterAnalyzer.cxx b/DQM/src/DQM/EcalClusterAnalyzer.cxx new file mode 100644 index 000000000..9813b381b --- /dev/null +++ b/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 +#include +#include + +#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("nbr_of_electrons"); + + ecal_sim_hit_coll_ = ps.getParameter("ecal_sim_hit_coll"); + ecal_sim_hit_pass_ = ps.getParameter("ecal_sim_hit_pass"); + + rec_hit_coll_name_ = ps.getParameter("rec_hit_coll_name"); + rec_hit_pass_name_ = ps.getParameter("rec_hit_pass_name"); + + cluster_coll_name_ = ps.getParameter("cluster_coll_name"); + cluster_pass_name_ = ps.getParameter("cluster_pass_name"); + return; +} + +void EcalClusterAnalyzer::analyze(const framework::Event& event) { + const auto& ecal_rec_hits{event.getCollection( + rec_hit_coll_name_, rec_hit_pass_name_)}; + const auto& ecal_sim_hits{event.getCollection( + ecal_sim_hit_coll_, ecal_sim_hit_pass_)}; + const auto& ecal_clusters{event.getCollection( + 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>> 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 pos1; + std::vector pos2; + bool p1 = false; + bool p2 = false; + + const auto& ecal_sp_hits{ + event.getCollection("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 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 n; + n.resize(nbr_of_electrons_ + 1); + // total number of energy coming from electron, index = electron ID + std::vector 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) diff --git a/Ecal/exampleConfigs/cluster.py b/Ecal/exampleConfigs/cluster.py new file mode 100644 index 000000000..04c840c3d --- /dev/null +++ b/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() +] + diff --git a/Ecal/include/Ecal/CLUE.h b/Ecal/include/Ecal/CLUE.h new file mode 100644 index 000000000..5700d8aa5 --- /dev/null +++ b/Ecal/include/Ecal/CLUE.h @@ -0,0 +1,159 @@ +/** + * @file CLUE.h + * @brief A version of CLUE (CMS) for clustering in ECal + * @author Ella Viirola, Lund University + */ + +#ifndef ECAL_CLUE_H_ +#define ECAL_CLUE_H_ + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Ecal/Event/EcalHit.h" +#include "Ecal/IntermediateCluster.h" +#include "Framework/Logger.h" + +namespace ecal { + +class CLUE { + enableLogging("CLUE"); + + public: + struct Density { + float x; + float y; + float 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 + + int layer; + + std::vector hits; + + Density() {} + + Density(float xx, float yy) : x(xx), y(yy) { + total_energy = 0.; + index = -1; + follower_of = -1; + delta = std::numeric_limits::max(); + z_delta = std::numeric_limits::max(); + cluster_id = -1; + layer = -1; + z = 0.; + hits = {}; + } + }; + + 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); + std::vector> createLayers( + const std::vector& hits); + float roundToDecimal(float x, int num_decimal_precision_digits); + std::vector> setup( + const std::vector& 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 + std::vector> clustering( + std::vector>& densities, bool connectingLayers, + int layerTag = 0); + + std::vector> layerSetup(); + + void convertToIntermediateClusters( + std::vector>& clusters); + + void cluster(const std::vector& hits, double dc, double rc, + double deltac, double deltao, int nbrOfLayers, + bool reclustering); + + std::vector getCentroidDistances() const { + return centroid_distances_; + } + + int getNLoops() const { return clustering_loops_; } + + int getInitialClusterNbr() const { return initial_cluster_nbr_; } + + std::vector getClusters() const { + return final_clusters_; + } + + // First layer centroids are available for potential future combination with + // TS + std::vector getFirstLayerCentroids() const { + return first_layer_centroids_; + } + + private: + int clustering_loops_; + + bool reclustering_; + + double dc_; + double rhoc_; + double deltac_; + double deltao_; + double dm_; + + // layers in Ecal; a bit unsure if this is the correct number + int max_layers_{17}; + int nbr_of_layers_; + // air between layers (a guess, but it seems to work) + double air_{10.}; + // thickness of ECal layers + std::vector 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 layer_rho_c_; + std::vector layer_delta_c; + // containment radius for the different layers of the ECal + std::vector radius{ + 5.723387467629167, 5.190678018534044, 5.927290663506518, + 6.182560329200212, 7.907549398117859, 8.606100542857211, + 10.93381822596916, 12.043201938160239, 14.784548371508041, + 16.102403056546482, 18.986402399412817, 20.224453740305716, + 23.048820910305643, 24.11202594672678, 26.765135236851666, + 27.78700483852502, 30.291794353801293, 31.409870873194464, + 33.91006482486666, 35.173073672355926, 38.172422630271, + 40.880288341493205, 44.696485719120005, 49.23802839743545, + 53.789910813378675, 60.87843355562641, 66.32931132415688, + 75.78117972604727, 86.04697356716805, 96.90360704034346}; + + std::vector centroid_distances_; + IntermediateCluster event_centroid_; + + float first_layer_max_z_; + std::vector first_layer_centroids_; + + int seed_index_{0}; + std::vector>> seeds_; + + int initial_cluster_nbr_{-1}; + std::vector final_clusters_; + std::vector> layer_centroid_separations_; +}; +} // namespace ecal + +#endif diff --git a/Ecal/include/Ecal/EcalClusterProducer.h b/Ecal/include/Ecal/EcalClusterProducer.h index b98f49853..1843f53be 100644 --- a/Ecal/include/Ecal/EcalClusterProducer.h +++ b/Ecal/include/Ecal/EcalClusterProducer.h @@ -10,8 +10,6 @@ //----------// // ROOT // //----------// -#include "TCanvas.h" -#include "TFile.h" #include "TH1F.h" #include "TH2F.h" @@ -21,12 +19,13 @@ #include "DetDescr/DetectorID.h" #include "DetDescr/EcalGeometry.h" #include "DetDescr/EcalID.h" +#include "Ecal/CLUE.h" #include "Ecal/Event/ClusterAlgoResult.h" #include "Ecal/Event/EcalCluster.h" #include "Ecal/Event/EcalHit.h" +#include "Ecal/IntermediateCluster.h" #include "Ecal/MyClusterWeight.h" #include "Ecal/TemplatedClusterFinder.h" -#include "Ecal/WorkingCluster.h" #include "Framework/Configure/Parameters.h" #include "Framework/EventProcessor.h" @@ -45,8 +44,7 @@ namespace ecal { class EcalClusterProducer : public framework::Producer { public: EcalClusterProducer(const std::string& name, framework::Process& process); - - virtual ~EcalClusterProducer(); + ~EcalClusterProducer() override = default; /** * Configure the processor using the given user specified parameters. @@ -55,17 +53,28 @@ class EcalClusterProducer : public framework::Producer { */ void configure(framework::config::Parameters& parameters) override; - virtual void produce(framework::Event& event) override; + void produce(framework::Event& event) override; private: - double seedThreshold_{0}; + double seed_threshold_{0}; double cutoff_{0}; - std::string digisPassName_; - std::string algoCollName_; - std::string clusterCollName_; + + double dc_{0}; + double rhoc_{0}; + double deltac_{0}; + double deltao_{0}; + + std::string rec_hit_coll_name_; + std::string rec_hit_pass_name_; + std::string algo_coll_name_; + std::string cluster_coll_name_; + + bool CLUE_; + int nbr_of_layers_; + bool reclustering_; /** The name of the cluster algorithm used. */ - TString algoName_; + TString algo_name_; }; } // namespace ecal diff --git a/Ecal/include/Ecal/Event/EcalCluster.h b/Ecal/include/Ecal/Event/EcalCluster.h index 94eb4b8c5..6c6ade345 100644 --- a/Ecal/include/Ecal/Event/EcalCluster.h +++ b/Ecal/include/Ecal/Event/EcalCluster.h @@ -34,15 +34,37 @@ class EcalCluster : public ldmx::CaloCluster { * @param hit The digi hit's entry number in the events digi * collection. */ - void addHits(const std::vector hitsVec); + void addHits(const std::vector& hits); + + void addFirstLayerHits(const std::vector& hits); bool operator<(const EcalCluster& rhs) const { return this->getEnergy() < rhs.getEnergy(); } + void setFirstLayerCentroidXYZ(double x, double y, double z) { + first_layer_centroid_x_ = x; + first_layer_centroid_y_ = y; + first_layer_centroid_z_ = z; + } + + double getFirstLayerCentroidX() const { return first_layer_centroid_x_; } + double getFirstLayerCentroidY() const { return first_layer_centroid_y_; } + double getFirstLayerCentroidZ() const { return first_layer_centroid_z_; } + + std::vector getFirstLayerHitIDs() const { + return first_layer_hit_IDs_; + } + private: // Could add further ECal-specific info here... + std::vector first_layer_hit_IDs_; + + double first_layer_centroid_x_{0}; + double first_layer_centroid_y_{0}; + double first_layer_centroid_z_{0}; + ClassDef(EcalCluster, 1); }; } // namespace ldmx diff --git a/Ecal/include/Ecal/IntermediateCluster.h b/Ecal/include/Ecal/IntermediateCluster.h new file mode 100644 index 000000000..403eab029 --- /dev/null +++ b/Ecal/include/Ecal/IntermediateCluster.h @@ -0,0 +1,37 @@ +/* + IntermediateCluster -- In-memory tool for working on clusters + */ +#ifndef ECAL_WORKINGCLUSTER_H_ +#define ECAL_WORKINGCLUSTER_H_ + +#include +#include + +#include "DetDescr/EcalGeometry.h" +#include "Ecal/Event/EcalHit.h" +#include "TLorentzVector.h" + +namespace ecal { + +class IntermediateCluster { + public: + IntermediateCluster(const ldmx::EcalHit& eh, int layer = -1); + IntermediateCluster() = default; + ~IntermediateCluster() = default; + void add(const ldmx::EcalHit& eh); + void add(const ldmx::EcalHit* eh); + void add(const IntermediateCluster& wc); + const ROOT::Math::XYZTVector& centroid() const { return centroid_; } + std::vector getHits() const { return hits_; } + bool empty() const { return hits_.empty(); } + void clear() { hits_.clear(); } + int getLayer() const { return layer_; } + + private: + int layer_; + std::vector hits_; + ROOT::Math::XYZTVector centroid_; +}; +} // namespace ecal + +#endif diff --git a/Ecal/include/Ecal/MyClusterWeight.h b/Ecal/include/Ecal/MyClusterWeight.h index 263c4861b..a32998b37 100644 --- a/Ecal/include/Ecal/MyClusterWeight.h +++ b/Ecal/include/Ecal/MyClusterWeight.h @@ -4,16 +4,16 @@ #include -#include "Ecal/WorkingCluster.h" +#include "Ecal/IntermediateCluster.h" namespace ecal { class MyClusterWeight { public: - double operator()( - const WorkingCluster& a, - const WorkingCluster& b) { // 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 diff --git a/Ecal/include/Ecal/TemplatedClusterFinder.h b/Ecal/include/Ecal/TemplatedClusterFinder.h index a8bf3f08e..0a2138754 100644 --- a/Ecal/include/Ecal/TemplatedClusterFinder.h +++ b/Ecal/include/Ecal/TemplatedClusterFinder.h @@ -7,9 +7,10 @@ #include +#include #include -#include "Ecal/WorkingCluster.h" +#include "Ecal/IntermediateCluster.h" #include "TH2F.h" namespace ecal { @@ -18,19 +19,21 @@ template class TemplatedClusterFinder { public: - void add(const ldmx::EcalHit* eh, const ldmx::EcalGeometry& hex) { - clusters_.push_back(WorkingCluster(eh, hex)); + void add(const ldmx::EcalHit& eh) { + clusters_.push_back(IntermediateCluster(eh)); } - static bool compClusters(const WorkingCluster& a, const WorkingCluster& b) { + static bool compClusters(const IntermediateCluster& a, + const IntermediateCluster& b) { return a.centroid().E() > b.centroid().E(); } void cluster(double seed_threshold, double cutoff) { int ncluster = clusters_.size(); double minwgt = cutoff; - + // Sort after highest energy std::sort(clusters_.begin(), clusters_.end(), compClusters); + loops_ = 0; do { bool any = false; size_t mi(0), mj(0); @@ -51,7 +54,8 @@ class TemplatedClusterFinder { for (size_t j = i + 1; j < clusters_.size(); j++) { if (clusters_[j].empty() || - (!iseed && clusters_[j].centroid().E() < seed_threshold)) + (!iseed && clusters_[j].centroid().E() < + seed_threshold)) // this never happens continue; double wgt = wgt_(clusters_[i], clusters_[j]); if (!any || wgt < minwgt) { @@ -77,25 +81,37 @@ class TemplatedClusterFinder { // decrement cluster count ncluster--; } - + loops_++; } while (minwgt < cutoff && ncluster > 1); finalwgt_ = minwgt; + for (const auto& cl : clusters_) { + if (!cl.empty() && cl.centroid().E() >= seed_threshold) + finalClusters_.push_back(cl); + else if (cl.centroid().E() < seed_threshold) + break; // clusters are sorted, so should be safe to break + } } double getYMax() const { return finalwgt_; } int getNSeeds() const { return nseeds_; } + int getNLoops() const { return loops_; } + std::map getWeights() const { return transitionWeights_; } - std::vector getClusters() const { return clusters_; } + std::vector getClusters() const { + return finalClusters_; + } private: WeightClass wgt_; double finalwgt_; int nseeds_; + int loops_; std::map transitionWeights_; - std::vector clusters_; + std::vector clusters_; + std::vector finalClusters_; }; } // namespace ecal diff --git a/Ecal/include/Ecal/WorkingCluster.h b/Ecal/include/Ecal/WorkingCluster.h deleted file mode 100644 index 92c2fed9b..000000000 --- a/Ecal/include/Ecal/WorkingCluster.h +++ /dev/null @@ -1,40 +0,0 @@ -/* - WorkingCluster -- In-memory tool for working on clusters - */ -#ifndef ECAL_WORKINGCLUSTER_H_ -#define ECAL_WORKINGCLUSTER_H_ - -#include -#include - -#include "DetDescr/EcalGeometry.h" -#include "Ecal/Event/EcalHit.h" -#include "TLorentzVector.h" - -namespace ecal { - -class WorkingCluster { - public: - WorkingCluster(const ldmx::EcalHit* eh, const ldmx::EcalGeometry& geom); - - ~WorkingCluster(){}; - - void add(const ldmx::EcalHit* eh, const ldmx::EcalGeometry& geom); - - void add(const WorkingCluster& wc); - - const TLorentzVector& centroid() const { return centroid_; } - - std::vector getHits() const { return hits_; } - - bool empty() const { return hits_.empty(); } - - void clear() { hits_.clear(); } - - private: - std::vector hits_; - TLorentzVector centroid_; -}; -} // namespace ecal - -#endif diff --git a/Ecal/python/ecalClusters.py b/Ecal/python/ecalClusters.py index e1b5bbb97..1e03d3e31 100644 --- a/Ecal/python/ecalClusters.py +++ b/Ecal/python/ecalClusters.py @@ -12,20 +12,46 @@ class EcalClusterProducer(ldmxcfg.Producer) : """Configure the clustering""" def __init__(self,name='ecalClusters') : - super().__init__(name,"ecal::EcalClusterProducer") + super().__init__(name,"ecal::EcalClusterProducer", 'Ecal') + # Pass name rec hits + self.rec_hit_coll_name = 'EcalRecHits' + self.rec_hit_pass_name = '' - self.cutoff = 10. - self.seedThreshold = 100.0 #MeV + # Name of the cluster collection to make + self.cluster_coll_name = "ecalClusters" - # Pass name for ecal digis - self.digisPassName = "recon" + # --- EXISTING ALGORITHM --- + self.cutoff = 10. + self.seed_threshold = 100.0 #MeV # Name of the algo to save to the root file - self.algoName = "MyClusterAlgo" - - # Name of the cluster collection to make - self.clusterCollName = "ecalClusters" - + self.algo_name = "MyClusterAlgo" # Name of the cluster algo collection to make - self.algoCollName = "ClusterAlgoResult" + self.algo_coll_name = "ClusterAlgoResult" + + # --- CLUE ALGORITHM --- + # Enable CLUE algorithm + self.CLUE = True + # Nbr of layers to perform CLUE on + # = 1 collapses all hits into same z-dimension, gives best results atm + self.nbr_of_layers = 1 + # Cutoff distance in calculation of local density + # Currently only used when nbrOfLayers > 1 + self.dc = 5. + # Minimum seed energy/maximum outlier energy + self.rhoc = 550. + # Minimum seed separation + 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 diff --git a/Ecal/src/Ecal/CLUE.cxx b/Ecal/src/Ecal/CLUE.cxx new file mode 100644 index 000000000..69af52a40 --- /dev/null +++ b/Ecal/src/Ecal/CLUE.cxx @@ -0,0 +1,564 @@ + +#include "Ecal/CLUE.h" + +namespace ecal { + +float CLUE::dist(double x1, double y1, double x2, double y2) { + return pow(pow(x1 - x2, 2) + pow(y1 - y2, 2), 0.5); +} + +float CLUE::floatDist(float x1, float y1, float x2, float y2) { + return powf(powf(x1 - x2, 2) + powf(y1 - y2, 2), 0.5); +} + +float CLUE::floatDist(float x1, float y1, float z1, float x2, float y2, + float z2) { + return powf(powf(x1 - x2, 2) + powf(y1 - y2, 2) + powf(z1 - z2, 2), 0.5); +} + +/* Old code, idea was to do electron reclustering based on first layer + centroids' distances to each other I.e. if electrons are close together => + likely merged => recluster Did not quite work and I don't remember the idea + anymore but leaving the code here for inspo */ + +// void electronSeparation(std::vector hits) { +// std::vector layerThickness = +// { 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 }; double air = 10.; std::sort(hits.begin(), hits.end(), [](const +// ldmx::EcalHit& a, const ldmx::EcalHit& b) { +// return a.getZPos() < b.getZPos(); +// }); +// std::vector firstLayers; +// std::vector firstLayerClusters; +// int layerTag = 0; +// double layerZ = hits[0].getZPos(); +// for (const auto& hit : hits) { +// if (hit.getZPos() > layerZ + layerThickness[layerTag] + air) { +// layerTag++; +// // if (layerTag > limit) break; +// break; +// } +// firstLayers.push_back(hit); +// firstLayerClusters.push_back(IntermediateEcalCluster(hit, layerTag)); + +// } +// bool merge = false; +// do { +// merge = false; +// for (int i = 0; i < firstLayerClusters.size(); i++) { +// if (firstLayerClusters[i].empty()) continue; +// // if (firstLayerClusters[i].centroid().E() >= seedThreshold_) { +// for (int j = i + 1; j < firstLayerClusters.size(); j++) { +// if (firstLayerClusters[j].empty()) continue; +// if (dist(firstLayerClusters[i].centroid().Px(), +// firstLayerClusters[i].centroid().Py(), +// firstLayerClusters[j].centroid().Px(), +// firstLayerClusters[j].centroid().Py()) < 8.) { +// firstLayerClusters[i].add(firstLayerClusters[j]); +// firstLayerClusters[j].clear(); +// merge = true; +// } +// } +// // } else break; +// } +// } while (merge); +// ldmx_log(trace) << "--- ELECTRON SEPARATION ---"; +// for (int i = 0; i < firstLayerClusters.size(); i++) { +// if (firstLayerClusters[i].empty()) continue; +// ldmx_log(trace) << " Cluster " << i << " x: " +// << firstLayerClusters[i].centroid().Px() << " y: " +// << firstLayerClusters[i].centroid().Py(); +// for (int j = i + 1; +// j < firstLayerClusters.size(); j++) { +// if (firstLayerClusters[j].empty()) continue; +// auto d = dist(firstLayerClusters[i].centroid().Px(), +// firstLayerClusters[i].centroid().Py(), +// firstLayerClusters[j].centroid().Px(), +// firstLayerClusters[j].centroid().Py()); +// ldmx_log(trace) << "Dist to cluster " << j << ": " << d; +// } +// } +// } + +std::vector> CLUE::createLayers( + const std::vector& hits) { + ldmx_log(trace) << "--- LAYER CREATION ---"; + std::vector> layers; + + int layerTag = 0; + int trueLayer = 0; + double layerZ = hits[0]->getZPos(); + double trueLayerZ = layerZ; + double maxZ = hits[hits.size() - 1]->getZPos(); + layers.push_back({}); + double layerSeparation = (maxZ - layerZ) / nbr_of_layers_; + ldmx_log(trace) << " Layer separation: " << layerSeparation; + ldmx_log(trace) << " Creating layer 0"; + + double highestEnergy = 0.; + int rhocFactor = 2.; + for (const auto& hit : hits) { + // If z of hit is in new layer, both calculated and real (we don't want to + // split in the middle of actual ecal layer) + if (layerTag != nbr_of_layers_ && + hit->getZPos() > (layerZ + layerSeparation) && + hit->getZPos() > trueLayerZ + layer_thickness_[trueLayer] + air_) { + layerZ = hit->getZPos(); + layers.push_back({}); + // Set seed threshold for layer to highest energy of layer / factor + // TODO: decide division factor + layer_rho_c_.push_back(highestEnergy / rhocFactor); + layerTag++; + ldmx_log(trace) << " Highest energy: " << highestEnergy << "\n" + << " Creating layer " << layerTag; + highestEnergy = 0.; + } + if (hit->getZPos() > trueLayerZ + layer_thickness_[trueLayer] + air_) { + // keep track of true layers + if (trueLayer == 0) first_layer_max_z_ = hit->getZPos(); + trueLayer++; + trueLayerZ = hit->getZPos(); + if (nbr_of_layers_ < 2) return layers; + } + layers[layerTag].push_back(hit); + if (hit->getEnergy() > highestEnergy) highestEnergy = hit->getEnergy(); + } + layer_rho_c_.push_back(highestEnergy / rhocFactor); + return layers; +} + +float CLUE::roundToDecimal(float x, int num_decimal_precision_digits) { + float power_of_10 = std::pow(10, num_decimal_precision_digits); + return std::round(x * power_of_10) / power_of_10; +} + +std::vector> CLUE::setup( + const std::vector& hits) { + std::vector> densities; + std::map, std::shared_ptr> densityMap; + event_centroid_ = IntermediateCluster(); + ldmx_log(trace) << "--- SETUP ---"; + ldmx_log(trace) << "Building densities"; + for (const auto& hit : hits) { + // collapse z dimension + float x = roundToDecimal(hit->getXPos(), 4); + float y = roundToDecimal(hit->getYPos(), 4); + ldmx_log(trace) << " New hit { x: " << x << " y: " << y << "}"; + std::pair coords; + if (dc_ != 0 && nbr_of_layers_ > 1) { + // if more than one layer, divide hits into densities with side dc_ + double i = std::ceil(std::abs(x) / dc_); + double j = std::ceil(std::abs(y) / dc_); + if (x < 0) { + i = -i; + x = (i + 0.5) * dc_; + } else + x = (i - 0.5) * dc_; + if (y < 0) { + j = -j; + y = (j + 0.5) * dc_; + } else + y = (1 - 0.5) * dc_; // set x,y to middle of box + coords = {i, j}; + ldmx_log(trace) << " Index " << i << ", " << j << "; x: " << x + << " y: " << y; + } else { + // if just one layer, have all densities with the same x,y be in same + // density + coords = {x, y}; + } + if (densityMap.find(coords) == densityMap.end()) { + densityMap.emplace(coords, std::make_shared(x, y)); + ldmx_log(trace) << " New density created"; + } else { + ldmx_log(trace) << " Found density with x: " << densityMap[coords]->x + << " y: " << densityMap[coords]->y; + } + densityMap[coords]->hits.push_back(hit); + densityMap[coords]->total_energy += hit->getEnergy(); + densityMap[coords]->z += hit->getZPos(); + + event_centroid_.add(hit); + } + + densities.reserve(densityMap.size()); + for (const auto& entry : densityMap) { + densities.push_back(std::move(entry.second)); + } + // sort according to energy + std::sort(densities.begin(), densities.end(), + [](const std::shared_ptr& a, + const std::shared_ptr& b) { + return a->total_energy > b->total_energy; + }); + + ldmx_log(trace) << "Decide parents"; + + // decide delta and followerOf + for (int i = 0; i < densities.size(); i++) { + densities[i]->index = i; + densities[i]->z = + densities[i]->z / densities[i]->hits.size(); // avg z position + ldmx_log(trace) << " Index: " << i << "; x: " << densities[i]->x + << "; y: " << densities[i]->y + << "; Energy: " << densities[i]->total_energy; + // loop through all higher energy densities + for (int j = 0; j < i; j++) { + float d = floatDist(densities[i]->x, densities[i]->y, densities[j]->x, + densities[j]->y); + // condition energyJ > energyI but this should be baked in as we sorted + // according to energy + if (d < dm_ && d < densities[i]->delta) { + densities[i]->delta = d; + densities[i]->follower_of = j; + ldmx_log(trace) << " New parent, index " << j + << "; delta: " << std::setprecision(20) << d; + } + } + } + return densities; +} + +// 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 +std::vector> CLUE::clustering( + std::vector>& densities, + bool connectingLayers, int layerTag) { + ldmx_log(trace) << "--- CLUSTERING ---"; + + if (!connectingLayers && nbr_of_layers_ > 1) { + // if doing layerwise clustering + rhoc_ = layer_rho_c_[layerTag]; + ldmx_log(trace) << "Setting rhoc on layer " << layerTag << " to " << rhoc_; + if (layerTag * 2 - 1 < radius.size()) { + deltac_ = radius[layerTag * 2 - 1]; + ldmx_log(trace) << "Setting deltac on layer " << layerTag << " to " + << deltac_; + } + } else if (connectingLayers) { + // if currently doing 3D clustering + // These values need to be modded; very random + deltao_ = 200.; + deltac_ = 100.; + rhoc_ = 1000.; + } + + bool energyOverload = false; + double maxEnergy = 10000.; + clustering_loops_ = 0; + double deltacMod = deltac_; + double centroidRadius = 10.; + + // stores seeds of this layer + std::vector>& layerSeeds = seeds_[layerTag]; + + // stores hits in cluster + std::vector> clusters; + // keeps track of which densities have merged; only used if reclustering + std::vector mergedDensities; // index = cluster id + mergedDensities.resize(densities.size()); + // keeps track of cluster energies + std::vector clusterEnergies; + do { + // while no cluster has merged + if (energyOverload) { + // makes delta c smaller if clusters have merged + deltacMod = deltacMod / 1.1; + ldmx_log(trace) << "Energy overload, new deltacmod: " << deltacMod; + energyOverload = false; + } + + clustering_loops_++; + ldmx_log(trace) << "Clustering loop " << clustering_loops_; + + // cluster index + int k = 0; + + layerSeeds.clear(); + layerSeeds.reserve(densities.size()); + clusters.clear(); + clusters.reserve(densities.size()); + clusterEnergies.clear(); + clusterEnergies.reserve(densities.size()); + + std::stack clusterStack; + // stores followers of densities at corr index + std::vector> followers; + followers.resize(densities.size()); + + // Mark as seed, follower or outlier + for (int j = 0; j < densities.size(); j++) { + // funky line to generalize this function for both 2D and 3D case + int i = densities[j]->index; + ldmx_log(trace) << " Index: " << i << "; x: " << densities[i]->x + << "; y: " << densities[i]->y + << "; Energy: " << densities[i]->total_energy; + ldmx_log(trace) << " Parent ID: " << densities[i]->follower_of + << "; Delta: " << densities[i]->delta; + + bool isSeed; + if (deltacMod != deltac_ && mergedDensities[densities[i]->cluster_id] && + floatDist(densities[i]->x, densities[i]->y, + event_centroid_.centroid().x(), + event_centroid_.centroid().y()) < centroidRadius) { + // if energy has been overloaded and this density belongs to cluster + // that was overloaded and this density is close enough to event + // centroid use modded delta c + isSeed = densities[i]->total_energy > rhoc_ && + densities[i]->delta > deltacMod; + } else + isSeed = + densities[i]->total_energy > rhoc_ && densities[i]->delta > deltac_; + + if (isSeed) { + ldmx_log(trace) << " Distance to centroid: " + << floatDist(densities[i]->x, densities[i]->y, + event_centroid_.centroid().Px(), + event_centroid_.centroid().Py()); + } + + bool isOutlier = + densities[i]->total_energy < rhoc_ && densities[i]->delta > deltao_; + + densities[i]->cluster_id = -1; + if (isSeed) { + ldmx_log(trace) << " SEED, cluster id " << k; + densities[i]->cluster_id = k; + k++; + clusterStack.push(i); + clusters.push_back(densities[i]->hits); + clusterEnergies.push_back(densities[i]->total_energy); + layerSeeds.push_back(densities[i]); + } else if (!isOutlier) { + ldmx_log(trace) << " Follower"; + int& parentIndex = densities[i]->follower_of; + if (parentIndex != -1) followers[parentIndex].push_back(i); + } else { + ldmx_log(trace) << " Outlier"; + } + } + + mergedDensities.clear(); + mergedDensities.resize(densities.size()); + + // Go through all seeds and add followers, then follower's followers, etc. + while (clusterStack.size() > 0) { + auto& d = densities[clusterStack.top()]; + clusterStack.pop(); + auto& cid = d->cluster_id; + // for indices of followers of dp + for (const auto& j : followers[d->index]) { + auto& f = densities[j]; + // set clusterindex of follower to clusterindex of d + f->cluster_id = cid; + clusterEnergies[cid] += f->total_energy; + if (reclustering_ && clusterEnergies[cid] > maxEnergy && + deltacMod > 0.5 && clustering_loops_ < 100) { + // if reclustering is on and cluster energy is too high and + // deltacmod is not too low and we haven't tried for too long + mergedDensities[cid] = true; + if (!energyOverload && clustering_loops_ == 99) + ldmx_log(warn) << "Merging clusters, max cluster loops hit"; + energyOverload = true; + if (clustering_loops_ != 1) + goto endwhile; // don't break on first loop to save initial + // cluster number + } + clusters[cid].insert(std::end(clusters[cid]), std::begin(f->hits), + std::end(f->hits)); + // add follower to stack, so its followers can also get correct + // clusterindex + clusterStack.push(j); + } + } + // for first clusteringloop, we want to save number of clusters before + // reclustering + if (clustering_loops_ == 1 && energyOverload) + initial_cluster_nbr_ = clusters.size(); + endwhile:; + } while (energyOverload); + // if we have more than one layer and we are not currently doing CLUE3D + if (!connectingLayers && nbr_of_layers_ > 1) { + // Overwrite seed densities' properties with cluster properties + // Might be cleaner to just create new densities for cluster seeds + for (auto& seed : layerSeeds) { + seed->delta = std::numeric_limits::max(); + seed->hits = clusters[seed->cluster_id]; + seed->total_energy = clusterEnergies[seed->cluster_id]; + seed->index = seed_index_; + seed_index_++; + } + // Sort seeds in layer based on energy + std::sort(layerSeeds.begin(), layerSeeds.end(), + [](const std::shared_ptr& a, + const std::shared_ptr& b) { + return a->total_energy > b->total_energy; + }); + } + return clusters; +} + +std::vector> CLUE::layerSetup() { + ldmx_log(trace) << "--- LAYER SETUP ---"; + std::vector> densities; + layer_rho_c_.clear(); + for (int layer = 0; layer < nbr_of_layers_; layer++) { + ldmx_log(trace) << " LAYER " << layer; + auto& currentLayer = seeds_[layer]; + double highestEnergy = 0.; + for (int i = 0; i < currentLayer.size(); i++) { + // for each seed in layer + currentLayer[i]->layer = layer; + if (currentLayer[i]->total_energy > highestEnergy) + highestEnergy = currentLayer[i]->total_energy; + ldmx_log(trace) << " Density with index " << currentLayer[i]->index + << ", energy: " << currentLayer[i]->total_energy; + int depth = 1; + // decide delta and followerof from seeds in previous and next layer + // do { + // depth++; + if (layer - depth >= 0) { + // look at previous layer + auto& previousLayer = seeds_[layer - depth]; + for (int j = 0; j < previousLayer.size(); j++) { + // for each seed in previous layer + auto d = floatDist(currentLayer[i]->x, currentLayer[i]->y, + previousLayer[j]->x, previousLayer[j]->y); + auto dz = std::abs(currentLayer[i]->z - previousLayer[i]->z); + ldmx_log(trace) << " Delta to index " << previousLayer[j]->index + << ": " << std::setprecision(20) << d; + ldmx_log(trace) << " DeltaZ to index " << previousLayer[j]->index + << ": " << std::setprecision(20) << dz << std::endl; + if (previousLayer[j]->total_energy > currentLayer[i]->total_energy && + d < currentLayer[i]->delta && dz < currentLayer[i]->z_delta) { + ldmx_log(trace) + << " New parent: index " << previousLayer[j]->index + << " on layer " << layer - depth << "; energy " + << previousLayer[j]->total_energy; + ldmx_log(trace) + << " New delta: " << std::setprecision(20) << d; + ldmx_log(trace) + << " New deltaZ: " << std::setprecision(20) << dz; + currentLayer[i]->delta = d; + currentLayer[i]->z_delta = dz; + currentLayer[i]->follower_of = previousLayer[j]->index; + } else if (previousLayer[j]->total_energy < + currentLayer[i]->total_energy) { + break; + } + } + } + if (layer + depth < nbr_of_layers_) { + auto& nextLayer = seeds_[layer + depth]; + for (int j = 0; j < nextLayer.size(); j++) { + auto d = floatDist(currentLayer[i]->x, currentLayer[i]->y, + nextLayer[j]->x, nextLayer[j]->y); + auto dz = std::abs(currentLayer[i]->z - nextLayer[i]->z); + ldmx_log(trace) << " Delta to index " << nextLayer[j]->index + << ": " << std::setprecision(20) << d; + ldmx_log(trace) << " DeltaZ to index " << nextLayer[j]->index + << ": " << std::setprecision(20) << dz; + if (nextLayer[j]->total_energy > currentLayer[i]->total_energy && + d < currentLayer[i]->delta && dz < currentLayer[i]->z_delta) { + ldmx_log(trace) << " New parent: index " << nextLayer[j]->index + << " on layer " << layer + depth << "; energy " + << nextLayer[j]->total_energy; + ldmx_log(trace) + << " New delta: " << std::setprecision(20) << d; + ldmx_log(trace) + << " New deltaZ: " << std::setprecision(20) << dz; + currentLayer[i]->delta = d; + currentLayer[i]->z_delta = dz; + currentLayer[i]->follower_of = nextLayer[j]->index; + } else if (nextLayer[j]->total_energy < + currentLayer[i]->total_energy) { + break; + } + } + } + // } while (currentLayer[i]->layerFollowerOf == -1 && (layer - depth >= + // 0 || layer + depth < nbrOfLayers_)); + densities.push_back(currentLayer[i]); + } + layer_rho_c_.push_back(highestEnergy / 2); + } + return densities; +} + +void CLUE::convertToIntermediateClusters( + std::vector>& clusters) { + // Convert to workingecalclusters to ensure compatibility with + // EcalClusterProducer + for (const auto& vec : clusters) { + auto c = IntermediateCluster(); + auto fc = IntermediateCluster(); + for (const auto& hit : vec) { + c.add(hit); + // if hit is in first layer, add to first layer cluster + if (hit->getZPos() < first_layer_max_z_) fc.add(hit); + } + final_clusters_.push_back(c); + first_layer_centroids_.push_back(fc); + const auto& d = + dist(c.centroid().Px(), c.centroid().Py(), + event_centroid_.centroid().Px(), event_centroid_.centroid().Py()); + centroid_distances_.push_back(d); + } +} + +void CLUE::cluster(const std::vector& unsorted_hits, double dc, + double rc, double deltac, double deltao, int nbrOfLayers, + bool reclustering) { + // cutoff distance for local density + dc_ = dc; + // min density to promote as seed/max density to demote as outlier + rhoc_ = rc; + // min separation distance for seeds + deltac_ = deltac; + // min separation distance for outliers + deltao_ = deltao; + dm_ = std::max(deltac, deltao); + + reclustering_ = reclustering; // Recluster merged clusters or not + nbr_of_layers_ = nbrOfLayers; + + if (nbr_of_layers_ < 1) + nbr_of_layers_ = max_layers_; // anything below 1 => include all layers + else if (nbr_of_layers_ > max_layers_) + nbr_of_layers_ = max_layers_; + + // first copy *addresses* so we are only ever passing around pointers + std::vector hits; + hits.reserve(unsorted_hits.size()); + for (const auto& eh : unsorted_hits) { + hits.push_back(&eh); + } + // sort hits by Z position + std::sort(hits.begin(), hits.end(), + [](const ldmx::EcalHit* a, const ldmx::EcalHit* b) { + return a->getZPos() < b->getZPos(); + }); + + seeds_.resize(nbrOfLayers); + const auto layers = createLayers(hits); + if (nbr_of_layers_ > 1) { + for (int i = 0; i < layers.size(); i++) { + ldmx_log(trace) << "--- LAYER " << i << " ---"; + auto densities = setup(layers[i]); + auto clusters = clustering(densities, false, i); + // convertToIntermediateClusters(clusters); // uncomment for layer + // clustering without 3D + } + // Below for CLUE3D, comment for just layer clustering + auto densities = layerSetup(); + auto clusters = clustering(densities, true); + convertToIntermediateClusters(clusters); + } else { + auto densities = setup(hits); + auto clusters = clustering(densities, false); + convertToIntermediateClusters(clusters); + } +} + +} // namespace ecal diff --git a/Ecal/src/Ecal/EcalClusterProducer.cxx b/Ecal/src/Ecal/EcalClusterProducer.cxx index 55a633154..e26216356 100644 --- a/Ecal/src/Ecal/EcalClusterProducer.cxx +++ b/Ecal/src/Ecal/EcalClusterProducer.cxx @@ -6,80 +6,132 @@ #include "Ecal/EcalClusterProducer.h" +#include + namespace ecal { EcalClusterProducer::EcalClusterProducer(const std::string& name, framework::Process& process) : Producer(name, process) {} -EcalClusterProducer::~EcalClusterProducer() {} - void EcalClusterProducer::configure(framework::config::Parameters& parameters) { cutoff_ = parameters.getParameter("cutoff"); - seedThreshold_ = parameters.getParameter("seedThreshold"); - digisPassName_ = parameters.getParameter("digisPassName"); - algoCollName_ = parameters.getParameter("algoCollName"); - algoName_ = parameters.getParameter("algoName"); - clusterCollName_ = parameters.getParameter("clusterCollName"); + seed_threshold_ = parameters.getParameter("seed_threshold"); + + dc_ = parameters.getParameter("dc"); + rhoc_ = parameters.getParameter("rhoc"); + deltac_ = parameters.getParameter("deltac"); + deltao_ = parameters.getParameter("deltao"); + + rec_hit_coll_name_ = + parameters.getParameter("rec_hit_coll_name"); + rec_hit_pass_name_ = + parameters.getParameter("rec_hit_pass_name"); + algo_coll_name_ = parameters.getParameter("algo_coll_name"); + algo_name_ = parameters.getParameter("algo_name"); + cluster_coll_name_ = + parameters.getParameter("cluster_coll_name"); + + CLUE_ = parameters.getParameter("CLUE"); + nbr_of_layers_ = parameters.getParameter("nbr_of_layers"); + reclustering_ = parameters.getParameter("reclustering"); } void EcalClusterProducer::produce(framework::Event& event) { - // Get the Ecal Geometry - const auto& geometry = getCondition( - ldmx::EcalGeometry::CONDITIONS_OBJECT_NAME); - - TemplatedClusterFinder cf; - - std::vector ecalHits = - event.getCollection("ecalDigis", digisPassName_); - int nEcalDigis = ecalHits.size(); - - // Don't do anything if there are no ECal digis! - if (!(nEcalDigis > 0)) { + const auto& ecal_hits{event.getCollection(rec_hit_coll_name_, + rec_hit_pass_name_)}; + if (ecal_hits.size() == 0) { + // don't do anything if there are no ECal hits return; } - for (ldmx::EcalHit& hit : ecalHits) { - // Skip zero energy digis. - if (hit.getEnergy() == 0) { - continue; + if (CLUE_) { + CLUE cf; + cf.cluster(ecal_hits, dc_, rhoc_, deltac_, deltao_, nbr_of_layers_, + reclustering_); + std::vector wcVec = cf.getClusters(); + std::vector fWcVec = cf.getFirstLayerCentroids(); + + auto nLoops = cf.getNLoops(); + histograms_.fill("nLoops", nLoops); + histograms_.fill("nClusters", wcVec.size()); + if (reclustering_) + histograms_.fill("recluster", cf.getInitialClusterNbr(), wcVec.size()); + + std::vector ecalClusters; + for (int aWC = 0; aWC < wcVec.size(); aWC++) { + ldmx::EcalCluster cluster; + + cluster.setEnergy(wcVec[aWC].centroid().E()); + cluster.setCentroidXYZ(wcVec[aWC].centroid().Px(), + wcVec[aWC].centroid().Py(), + wcVec[aWC].centroid().Pz()); + cluster.setFirstLayerCentroidXYZ(fWcVec[aWC].centroid().Px(), + fWcVec[aWC].centroid().Py(), + fWcVec[aWC].centroid().Pz()); + cluster.setNHits(wcVec[aWC].getHits().size()); + cluster.addHits(wcVec[aWC].getHits()); + cluster.addFirstLayerHits(fWcVec[aWC].getHits()); + + histograms_.fill("nHits", wcVec[aWC].getHits().size()); + histograms_.fill("cluster_energy", wcVec[aWC].centroid().E()); + + ecalClusters.push_back(cluster); } - cf.add(&hit, geometry); - } + event.add(cluster_coll_name_, ecalClusters); + } else { + TemplatedClusterFinder cf; - cf.cluster(seedThreshold_, cutoff_); - std::vector wcVec = cf.getClusters(); + for (const ldmx::EcalHit& hit : ecal_hits) { + // Skip zero energy digis. + if (hit.getEnergy() == 0) { + continue; + } + cf.add(hit); + } - std::map cWeights = cf.getWeights(); + cf.cluster(seed_threshold_, cutoff_); + std::vector wcVec = cf.getClusters(); - ldmx::ClusterAlgoResult algoResult; - algoResult.set(algoName_, 3, cWeights.rbegin()->first); - algoResult.setAlgoVar(0, cutoff_); - algoResult.setAlgoVar(1, seedThreshold_); - algoResult.setAlgoVar(2, cf.getNSeeds()); + auto nLoops = cf.getNLoops(); + histograms_.fill("nLoops", nLoops); + histograms_.fill("nClusters", wcVec.size()); - std::map::iterator it = cWeights.begin(); - for (it = cWeights.begin(); it != cWeights.end(); it++) { - algoResult.setWeight(it->first, it->second / 100); - } + std::map cWeights = cf.getWeights(); - std::vector ecalClusters; - for (int aWC = 0; aWC < wcVec.size(); aWC++) { - ldmx::EcalCluster cluster; + ldmx::ClusterAlgoResult algoResult; + algoResult.set(algo_name_, 3, cWeights.rbegin()->first); + algoResult.setAlgoVar(0, cutoff_); + algoResult.setAlgoVar(1, seed_threshold_); + algoResult.setAlgoVar(2, cf.getNSeeds()); - cluster.setEnergy(wcVec[aWC].centroid().E()); - cluster.setCentroidXYZ(wcVec[aWC].centroid().Px(), - wcVec[aWC].centroid().Py(), - wcVec[aWC].centroid().Pz()); - cluster.setNHits(wcVec[aWC].getHits().size()); - cluster.addHits(wcVec[aWC].getHits()); + std::map::iterator it = cWeights.begin(); + for (it = cWeights.begin(); it != cWeights.end(); it++) { + algoResult.setWeight(it->first, it->second / 100); + histograms_.fill("seed_weights", it->first, it->second); + } - ecalClusters.push_back(cluster); - } + std::vector ecalClusters; + for (int aWC = 0; aWC < wcVec.size(); aWC++) { + ldmx::EcalCluster cluster; + + cluster.setEnergy(wcVec[aWC].centroid().E()); + cluster.setCentroidXYZ(wcVec[aWC].centroid().Px(), + wcVec[aWC].centroid().Py(), + wcVec[aWC].centroid().Pz()); + cluster.setNHits(wcVec[aWC].getHits().size()); + cluster.addHits(wcVec[aWC].getHits()); + + histograms_.fill("nHits", wcVec[aWC].getHits().size()); + histograms_.fill("cluster_energy", wcVec[aWC].centroid().E()); - event.add(clusterCollName_, ecalClusters); - event.add(algoCollName_, algoResult); + ecalClusters.push_back(cluster); + } + + event.add(cluster_coll_name_, ecalClusters); + event.add(algo_coll_name_, algoResult); + } } } // namespace ecal diff --git a/Ecal/src/Ecal/Event/EcalCluster.cxx b/Ecal/src/Ecal/Event/EcalCluster.cxx index e18d0c41a..d1ef04874 100644 --- a/Ecal/src/Ecal/Event/EcalCluster.cxx +++ b/Ecal/src/Ecal/Event/EcalCluster.cxx @@ -1,17 +1,28 @@ #include "Ecal/Event/EcalCluster.h" -ClassImp(ldmx::EcalCluster) +ClassImp(ldmx::EcalCluster); - namespace ldmx { - EcalCluster::EcalCluster() {} +namespace ldmx { - EcalCluster::~EcalCluster() { Clear(); } +EcalCluster::EcalCluster() {} - void EcalCluster::addHits(const std::vector hitsVec) { - std::vector vecIDs; - for (int iHit = 0; iHit < hitsVec.size(); iHit++) { - vecIDs.push_back(hitsVec[iHit]->getID()); - } - setIDs(vecIDs); +EcalCluster::~EcalCluster() { Clear(); } + +void EcalCluster::addHits(const std::vector& hits) { + std::vector ids; + ids.reserve(hits.size()); + for (const auto& h : hits) { + ids.push_back(h->getID()); + } + setIDs(ids); +} + +void EcalCluster::addFirstLayerHits(const std::vector& hits) { + first_layer_hit_IDs_.clear(); + first_layer_hit_IDs_.reserve(hits.size()); + for (const auto& h : hits) { + first_layer_hit_IDs_.push_back(h->getID()); } +} + } // namespace ldmx diff --git a/Ecal/src/Ecal/IntermediateCluster.cxx b/Ecal/src/Ecal/IntermediateCluster.cxx new file mode 100644 index 000000000..45c2d874b --- /dev/null +++ b/Ecal/src/Ecal/IntermediateCluster.cxx @@ -0,0 +1,47 @@ + +#include "Ecal/IntermediateCluster.h" + +#include + +namespace ecal { + +IntermediateCluster::IntermediateCluster(const ldmx::EcalHit& eh, int layer) + : layer_{layer}, hits_{}, centroid_{} { + add(eh); +} + +void IntermediateCluster::add(const ldmx::EcalHit& eh) { + hits_.push_back(&eh); + + double hitE = eh.getEnergy(); + double hitX = eh.getXPos(); + double hitY = eh.getYPos(); + double hitZ = eh.getZPos(); + + double newE = hitE + centroid_.E(); + centroid_.SetXYZT((centroid_.x() * centroid_.E() + hitE * hitX) / newE, + (centroid_.y() * centroid_.E() + hitE * hitY) / newE, + (centroid_.z() * centroid_.E() + hitE * hitZ) / newE, newE); +} + +void IntermediateCluster::add(const ldmx::EcalHit* eh) { + if (eh != nullptr) add(*eh); +} + +void IntermediateCluster::add(const IntermediateCluster& wc) { + double newE = wc.centroid().E() + centroid_.E(); + centroid_.SetXYZT( + (centroid_.x() * centroid_.E() + wc.centroid().x() * wc.centroid().E()) / + newE, + (centroid_.y() * centroid_.E() + wc.centroid().y() * wc.centroid().E()) / + newE, + (centroid_.z() * centroid_.E() + wc.centroid().z() * wc.centroid().E()) / + newE, + newE); + + for (const auto eh : wc.getHits()) { + hits_.push_back(eh); + } +} + +} // namespace ecal diff --git a/Ecal/src/Ecal/WorkingCluster.cxx b/Ecal/src/Ecal/WorkingCluster.cxx deleted file mode 100644 index 9fdc3eb08..000000000 --- a/Ecal/src/Ecal/WorkingCluster.cxx +++ /dev/null @@ -1,55 +0,0 @@ -/* - WorkingCluster -- In-memory tool for working on clusters - */ - -#include "Ecal/WorkingCluster.h" - -#include - -namespace ecal { - -WorkingCluster::WorkingCluster(const ldmx::EcalHit* eh, - const ldmx::EcalGeometry& hex) { - add(eh, hex); -} - -void WorkingCluster::add(const ldmx::EcalHit* eh, - const ldmx::EcalGeometry& hex) { - double hitE = eh->getEnergy(); - - /// The ID number is implicitly converted to EcalID - auto [hitX, hitY, hitZ] = hex.getPosition(eh->getID()); - - double newE = hitE + centroid_.E(); - double newCentroidX = (centroid_.Px() * centroid_.E() + hitE * hitX) / newE; - double newCentroidY = (centroid_.Py() * centroid_.E() + hitE * hitY) / newE; - double newCentroidZ = (centroid_.Pz() * centroid_.E() + hitE * hitZ) / newE; - - centroid_.SetPxPyPzE(newCentroidX, newCentroidY, newCentroidZ, newE); - - hits_.push_back(eh); -} - -void WorkingCluster::add(const WorkingCluster& wc) { - double clusterE = wc.centroid().E(); - double centroidX = wc.centroid().Px(); - double centroidY = wc.centroid().Py(); - double centroidZ = wc.centroid().Pz(); - - double newE = clusterE + centroid_.E(); - double newCentroidX = - (centroid_.Px() * centroid_.E() + clusterE * centroidX) / newE; - double newCentroidY = - (centroid_.Py() * centroid_.E() + clusterE * centroidY) / newE; - double newCentroidZ = - (centroid_.Pz() * centroid_.E() + clusterE * centroidZ) / newE; - - centroid_.SetPxPyPzE(newCentroidX, newCentroidY, newCentroidZ, newE); - - std::vector clusterHits = wc.getHits(); - - for (size_t i = 0; i < clusterHits.size(); i++) { - hits_.push_back(clusterHits[i]); - } -} -} // namespace ecal diff --git a/Recon/include/Recon/PFEcalClusterProducer.h b/Recon/include/Recon/PFEcalClusterProducer.h index 9eedaffd0..29a719a5d 100644 --- a/Recon/include/Recon/PFEcalClusterProducer.h +++ b/Recon/include/Recon/PFEcalClusterProducer.h @@ -11,9 +11,6 @@ #include "DetDescr/EcalGeometry.h" #include "Ecal/Event/EcalCluster.h" #include "Ecal/Event/EcalHit.h" -#include "Ecal/MyClusterWeight.h" -#include "Ecal/TemplatedClusterFinder.h" -#include "Ecal/WorkingCluster.h" #include "Framework/Configure/Parameters.h" // Needed to import parameters from configuration file #include "Framework/Event.h" #include "Framework/EventProcessor.h" //Needed to declare processor diff --git a/SimCore/include/SimCore/Event/SimCalorimeterHit.h b/SimCore/include/SimCore/Event/SimCalorimeterHit.h index 444c5184a..80417d7bf 100644 --- a/SimCore/include/SimCore/Event/SimCalorimeterHit.h +++ b/SimCore/include/SimCore/Event/SimCalorimeterHit.h @@ -66,6 +66,9 @@ class SimCalorimeterHit { /// Time this contributor made the hit (global Geant4 time) float time{0}; + + /// Saves ID of electron incident particle came from + int originID{-1}; }; /** @@ -241,7 +244,7 @@ class SimCalorimeterHit { * @param time The time of the hit [ns]. */ void addContrib(int incidentID, int trackID, int pdgCode, float edep, - float time); + float time, int originID = -1); /** * Get a hit contribution by index. @@ -357,6 +360,11 @@ class SimCalorimeterHit { */ std::vector timeContribs_; + /** + * The list of origin IDs contributing to the hit. + */ + std::vector originContribs_; + /** * The number of hit contributions. */ @@ -404,7 +412,7 @@ class SimCalorimeterHit { /** * ROOT class definition. */ - ClassDef(SimCalorimeterHit, 5) + ClassDef(SimCalorimeterHit, 6) }; } // namespace ldmx diff --git a/SimCore/include/SimCore/SDs/EcalSD.h b/SimCore/include/SimCore/SDs/EcalSD.h index 12ae1ed3a..486b32e6a 100644 --- a/SimCore/include/SimCore/SDs/EcalSD.h +++ b/SimCore/include/SimCore/SDs/EcalSD.h @@ -83,6 +83,8 @@ class EcalSD : public SensitiveDetector { bool enableHitContribs_; /// compress hit contribs bool compressHitContribs_; + /// maximum track ID to be considered an "origin" + int max_origin_track_id_; }; } // namespace simcore diff --git a/SimCore/python/sensitive_detectors.py b/SimCore/python/sensitive_detectors.py index 08f549e79..17dcfd0a3 100644 --- a/SimCore/python/sensitive_detectors.py +++ b/SimCore/python/sensitive_detectors.py @@ -111,6 +111,7 @@ def __init__(self) : super().__init__('ecal_sd', 'simcore::EcalSD','SimCore_SDs') self.enableHitContribs = True self.compressHitContribs = True + self.max_origin_track_id = 6 class TrigScintSD(simcfg.SensitiveDetector) : """Trigger Scintillaotr Sensitive Detector diff --git a/SimCore/src/SimCore/Event/SimCalorimeterHit.cxx b/SimCore/src/SimCore/Event/SimCalorimeterHit.cxx index 0ba29b012..b066d404a 100644 --- a/SimCore/src/SimCore/Event/SimCalorimeterHit.cxx +++ b/SimCore/src/SimCore/Event/SimCalorimeterHit.cxx @@ -35,12 +35,13 @@ ClassImp(ldmx::SimCalorimeterHit) } void SimCalorimeterHit::addContrib(int incidentID, int trackID, int pdgCode, - float edep, float time) { + float edep, float time, int originID) { incidentIDContribs_.push_back(incidentID); trackIDContribs_.push_back(trackID); pdgCodeContribs_.push_back(pdgCode); edepContribs_.push_back(edep); timeContribs_.push_back(time); + originContribs_.push_back(originID); edep_ += edep; if (time < time_ || time_ == 0) { time_ = time; @@ -55,6 +56,7 @@ ClassImp(ldmx::SimCalorimeterHit) contrib.edep = edepContribs_.at(i); contrib.time = timeContribs_.at(i); contrib.pdgCode = pdgCodeContribs_.at(i); + contrib.originID = originContribs_.at(i); return contrib; } diff --git a/SimCore/src/SimCore/SDs/EcalSD.cxx b/SimCore/src/SimCore/SDs/EcalSD.cxx index b0244650e..d9cba4ade 100644 --- a/SimCore/src/SimCore/SDs/EcalSD.cxx +++ b/SimCore/src/SimCore/SDs/EcalSD.cxx @@ -21,6 +21,7 @@ EcalSD::EcalSD(const std::string& name, simcore::ConditionsInterface& ci, : SensitiveDetector(name, ci, p) { enableHitContribs_ = p.getParameter("enableHitContribs"); compressHitContribs_ = p.getParameter("compressHitContribs"); + max_origin_track_id_ = p.getParameter("max_origin_track_id"); } G4bool EcalSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) { @@ -112,8 +113,19 @@ G4bool EcalSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) { if (compressHitContribs_ and contrib_i != -1) { hit.updateContrib(contrib_i, edep, time); } else { - hit.addContrib(getTrackMap().findIncident(track_id), track_id, pdg, edep, - time); + auto map{getTrackMap()}; + auto incident{map.findIncident(track_id)}; + // default "origin" is just the same as incident + // "origin" checks if a hit "originates" from one of the earliest + // track IDs (i.e. probably one of the primaries) + int origin{incident}; + for (int i{1}; i < max_origin_track_id_; ++i) { + if (map.isDescendant(track_id, i, 100)) { + origin = i; + break; + } + } + hit.addContrib(incident, track_id, pdg, edep, time, origin); } } else { // no hit contribs and hit already exists