From 335aef81376178c9f81975102305e09c373085b1 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 13:42:26 +0200 Subject: [PATCH 01/47] Remove comment about unused hit association --- Recon/include/Recon/Event/CaloCluster.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Recon/include/Recon/Event/CaloCluster.h b/Recon/include/Recon/Event/CaloCluster.h index 9e7144d933..426656e1fa 100644 --- a/Recon/include/Recon/Event/CaloCluster.h +++ b/Recon/include/Recon/Event/CaloCluster.h @@ -124,7 +124,7 @@ class CaloCluster { double getEDYDZ() const { return errDYDZ_; } - // get hit rawIDs (unused) + // get hit rawIDs const std::vector& getHitIDs() const { return hitIDs_; } // ability to store limited hit info From d3446c4b44b24e8e3eb80c68394e258790fbac5d Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 13:45:01 +0200 Subject: [PATCH 02/47] Add hit, cluster, track indices to pf cand --- Recon/include/Recon/Event/PFCandidate.h | 45 +++++++++++++++++++++++++ Recon/src/Recon/ParticleFlow.cxx | 4 ++- 2 files changed, 48 insertions(+), 1 deletion(-) diff --git a/Recon/include/Recon/Event/PFCandidate.h b/Recon/include/Recon/Event/PFCandidate.h index 49a24fe93a..96a5f67285 100644 --- a/Recon/include/Recon/Event/PFCandidate.h +++ b/Recon/include/Recon/Event/PFCandidate.h @@ -9,6 +9,9 @@ // ROOT #include "TObject.h" //For ClassDef +// ldmx-sw objects +//#include "Ecal/Event/EcalHit.h" +//#include "Hcal/Event/HcalHit.h" namespace ldmx { @@ -49,6 +52,11 @@ class PFCandidate { posHcalZ_ = z; } + // associate component indices to the pf candidate + void setTrackIndex(int x) { track_idx_ = x; } + void setEcalIndex(int x) { ecal_idx_ = x; } + void setHcalIndex(int x) { hcal_idx_ = x; } + void setTrackPxPyPz(float x, float y, float z) { trackPx_ = x; trackPy_ = y; @@ -103,6 +111,19 @@ class PFCandidate { void setTruthEnergy(double x) { truthEnergy_ = x; } void setTruthPdgId(int x) { truthPdgId_ = x; } + /** + * Take in the ecal hits that make up the candidate. + * @param hit The digi hit's entry number in the events digi + * collection. + */ + //void setEcalHits(const std::vector hits) { ecal_hits_ = hits; } + + /** + * Take in the hcal hits that make up the candidate. + * @param hit The digi hit's entry number in the events digi + * collection. + */ + // void setHcalHits(const std::vector hits) { hcal_hits_ = hits; } /* Getters */ @@ -119,6 +140,10 @@ class PFCandidate { std::vector getHcalPositionXYZ() const { return {posHcalX_, posHcalY_, posHcalZ_}; } + // associate component indices to the pf candidate + int getTrackIndex() { return track_idx_; } + int getEcalIndex() { return ecal_idx_ ; } + int getHcalIndex() { return hcal_idx_ ; } std::vector getTrackPxPyPz() const { return {trackPx_, trackPy_, trackPz_}; @@ -160,6 +185,21 @@ class PFCandidate { double getTruthEnergy() { return truthEnergy_; } int getTruthPdgId() { return truthPdgId_; } + /** + * Take in the ecal hits that make up the candidate. + * @param hit The digi hit's entry number in the events digi + * collection. + */ + //std::vector getEcalHits() { return ecal_hits_; } + + /** + * Take in the hcal hits that make up the candidate. + * @param hit The digi hit's entry number in the events digi + * collection. + */ + // std::vector getHcalHits() { return hcal_hits_; } + + private: /* Particle ID enum */ int pid_{0}; @@ -223,6 +263,11 @@ class PFCandidate { double truthEnergy_{0}; int truthPdgId_{0}; + /* Indices of the components making up the PFlow object */ + int track_idx_{-1}; + int ecal_idx_{-1}; + int hcal_idx_{-1}; + /* The ROOT class definition. */ ClassDef(PFCandidate, 1); }; diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx index 1736f3714d..8a77a4f9f1 100644 --- a/Recon/src/Recon/ParticleFlow.cxx +++ b/Recon/src/Recon/ParticleFlow.cxx @@ -312,7 +312,7 @@ void ParticleFlow::produce(framework::Event& event) { for (int i = 0; i < tracks.size(); i++) { ldmx::PFCandidate cand; fillCandTrack(cand, tracks[i]); // append track info to candidate - + cand.setTrackIndex(i); if (!tkIsEMLinked[i]) { // chargedUnmatch.push_back(cand); } else { // if track is linked with ECal cluster @@ -334,6 +334,7 @@ void ParticleFlow::produce(framework::Event& event) { ldmx::PFCandidate cand; fillCandEMCalo(cand, ecalClusters[i]); + cand.setEcalIndex(i); if (EMIsHadLinked[tkEMPairs[i]]) { fillCandHadCalo(cand, hcalClusters[EMHadPairs[i]]); // emMatch.push_back(cand); @@ -347,6 +348,7 @@ void ParticleFlow::produce(framework::Event& event) { if (HadIsEMLinked[i]) continue; ldmx::PFCandidate cand; fillCandHadCalo(cand, hcalClusters[i]); + cand.setHcalIndex(i); // hadOnly.push_back(cand); pfCands.push_back(cand); } From b31e0e39ee61b58a0fad31a5ddab69ded5feebb5 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 15:16:03 +0200 Subject: [PATCH 03/47] Add first attempts of pileup hit finding --- Recon/include/Recon/PileupFinder.h | 55 +++++++++++++++++ Recon/src/Recon/PileupFinder.cxx | 97 ++++++++++++++++++++++++++++++ 2 files changed, 152 insertions(+) create mode 100644 Recon/include/Recon/PileupFinder.h create mode 100644 Recon/src/Recon/PileupFinder.cxx diff --git a/Recon/include/Recon/PileupFinder.h b/Recon/include/Recon/PileupFinder.h new file mode 100644 index 0000000000..657f22bd37 --- /dev/null +++ b/Recon/include/Recon/PileupFinder.h @@ -0,0 +1,55 @@ +/** + * @file PileupFinder.h + * @brief Simple PFlow algorithm + * @author Christian Herwig, Fermilab + */ + +#ifndef PARTICLEFLOW_H +#define PARTICLEFLOW_H + +// LDMX Framework +#include "Ecal/Event/EcalCluster.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 +#include "Hcal/Event/HcalCluster.h" +#include "Recon/Event/CaloCluster.h" +#include "Recon/Event/PFCandidate.h" +#include "SimCore/Event/SimParticle.h" +#include "SimCore/Event/SimTrackerHit.h" +#include "TGraph.h" + +namespace recon { + +/** + * @class PileupFinder + * @brief + */ +class PileupFinder : public framework::Producer { + public: + PileupFinder(const std::string& name, framework::Process& process) + : framework::Producer(name, process) {} + + virtual void configure(framework::config::Parameters& ps); + + virtual void produce(framework::Event& event); + + virtual void onProcessEnd(); + + private: + + // name of collections for PF input object to be passed + std::string rec_hit_coll_name_; + std::string rec_hit_pass_name_; + std::string pf_cand_coll_name_; + std::string pf_cand_pass_name_; + std::string cluster_coll_name_; + std::string cluster_pass_name_; + // name of collection for pileup-free output hit coll + std::string output_rec_hit_coll_name_; + // configuration + +}; +} // namespace recon + +#endif /* PILEUPFINDER_H */ diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx new file mode 100644 index 0000000000..554c162f64 --- /dev/null +++ b/Recon/src/Recon/PileupFinder.cxx @@ -0,0 +1,97 @@ +#include "Recon/PileupFinder.h" + +#include + +namespace recon { + +void PileupFinder::configure(framework::config::Parameters& ps) { + // I/O + rec_hit_coll_name_ = ps.getParameter("rec_hit_coll_name"); + rec_hit_pass_name_ = ps.getParameter("rec_hit_pass_name"); + pf_cand_coll_name_ = ps.getParameter("pf_cand_coll_name"); + pf_cand_pass_name_ = ps.getParameter("pf_cand_pass_name"); + cluster_coll_name_ = ps.getParameter("cluster_coll_name"); + cluster_pass_name_ = ps.getParameter("cluster_pass_name"); + output_rec_hit_coll_name_ = ps.getParameter("output_rec_hit_coll_name"); + // Algorithm configuration + min_mom_ = ps.getParameter("min_momentum"); +} + +// get pileup candidates from PFlow and make a cleaned-up hit collection +void PileupFinder::produce(framework::Event& event) { + if (!event.exists(rec_hit_coll_name_,rec_hit_pass_name_)) { //ecal rechits + ldmx_log(error) << "Unable to find (one) collection named " << rec_hit_coll_name << + "_" << rec_hit_pass_name; + return; + } + if (!event.exists(pf_cand_coll_name_, pf_cand_pass_name_ )) { + ldmx_log(error) << "Unable to find (one) collection named " << pf_cand_coll_name << + "_" << pf_cand_pass_name; + return; + } + if (!event.exists(cluster_coll_name_, cluster_pass_name_ )) { + ldmx_log(error) << "Unable to find (one) collection named " << cluster_coll_name << + "_" << cluster_pass_name; + return; + } + + const auto& ecal_hits{event.getCollection(rec_hit_coll_name_, + rec_hit_pass_name_)}; + + const auto& pf_cands{event.getCollection(pf_cand_coll_name_, + pf_cand_pass_name_)}; + + const auto& clusters{event.getCollection(cluster_coll_name_, + cluster_pass_name_)}; + // get PID 3 and 7 -- the ones with track and ecal matching + // get the high-momentum track ones from there -- pileup candidate! + // get the list of hits associated with pileup candidates + // if a rechit is not on that list, add to output collection. + std::vector output_hits; + + for (const auto& pf_cand : pfCandidates){ + if(pf_cand.getPID()==3 || pf_cand.getPID()==7){ + + double mom_vec = pf_cand.getTrackPxPyPz(); + float mom = mom_vec[0]*mom_vec[0]+mom_vec[1]*mom_vec[1]+mom_vec[2]*mom_vec[2]; + mom = sqrt(mom); + + if (mom < min_mom_) + continue; + ldmx_log(trace) << "Got pileup candidate with PID = " << pf_cand.getPID() << " and momentum = " << mom << " MeV." + + // now! use the hit-candidate association to get the associated ecal hits. + + + // TODO: + // write a header for this file --> DONE + // clean up all the methods below here --> DONE + // compile, possibly incrememnt particleflow candidate class imp nb + // try out hit association (print list?) + int pf_cl_idx = pf_cand.getEcalIndex(); + auto cl = clusters[pf_cl_idx]; + auto hitIds = cl.getHitIDs(); + // make a collection without pileup hits + for (auto hit : ecal_hits ) { + + int *foundIndex = std::find(std::begin(hitIDs), std::end(hitIDs), hit.id()); + // When the element is not found, std::find returns the end of the range + if (foundIndex == std::end(hitIDs)) { //hit not found in the pileup cluster + output_hits.emplace_back(hit); //keep it + } + }// if trk/ecal matched + }// over PF objects + + } + +void PileupFinder::onProcessEnd() { + ldmx_log(debug) << "Process ends!"; + delete eCorr_; + delete hCorr_; + + return; +} + +} // namespace recon + +DECLARE_PRODUCER(recon::PileupFinder); From 2477dbcbc0426cc28814e9ab27f67c3937339152 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 15:26:38 +0200 Subject: [PATCH 04/47] Fix compilation problems --- Recon/include/Recon/PileupFinder.h | 1 + Recon/src/Recon/PileupFinder.cxx | 30 +++++++++++++++--------------- 2 files changed, 16 insertions(+), 15 deletions(-) diff --git a/Recon/include/Recon/PileupFinder.h b/Recon/include/Recon/PileupFinder.h index 657f22bd37..120b1dada4 100644 --- a/Recon/include/Recon/PileupFinder.h +++ b/Recon/include/Recon/PileupFinder.h @@ -49,6 +49,7 @@ class PileupFinder : public framework::Producer { std::string output_rec_hit_coll_name_; // configuration + float min_mom_{0.}; //MeV }; } // namespace recon diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index 554c162f64..5ab4f16f54 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -20,18 +20,18 @@ void PileupFinder::configure(framework::config::Parameters& ps) { // get pileup candidates from PFlow and make a cleaned-up hit collection void PileupFinder::produce(framework::Event& event) { if (!event.exists(rec_hit_coll_name_,rec_hit_pass_name_)) { //ecal rechits - ldmx_log(error) << "Unable to find (one) collection named " << rec_hit_coll_name << - "_" << rec_hit_pass_name; + ldmx_log(error) << "Unable to find (one) collection named " << rec_hit_coll_name_ << + "_" << rec_hit_pass_name_; return; } if (!event.exists(pf_cand_coll_name_, pf_cand_pass_name_ )) { - ldmx_log(error) << "Unable to find (one) collection named " << pf_cand_coll_name << - "_" << pf_cand_pass_name; + ldmx_log(error) << "Unable to find (one) collection named " << pf_cand_coll_name_ << + "_" << pf_cand_pass_name_; return; } if (!event.exists(cluster_coll_name_, cluster_pass_name_ )) { - ldmx_log(error) << "Unable to find (one) collection named " << cluster_coll_name << - "_" << cluster_pass_name; + ldmx_log(error) << "Unable to find (one) collection named " << cluster_coll_name_ << + "_" << cluster_pass_name_; return; } @@ -49,7 +49,7 @@ void PileupFinder::produce(framework::Event& event) { // if a rechit is not on that list, add to output collection. std::vector output_hits; - for (const auto& pf_cand : pfCandidates){ + for (const auto& pf_cand : pf_cands){ if(pf_cand.getPID()==3 || pf_cand.getPID()==7){ double mom_vec = pf_cand.getTrackPxPyPz(); @@ -70,20 +70,20 @@ void PileupFinder::produce(framework::Event& event) { // try out hit association (print list?) int pf_cl_idx = pf_cand.getEcalIndex(); auto cl = clusters[pf_cl_idx]; - auto hitIds = cl.getHitIDs(); + auto hitIDs = cl.getHitIDs(); // make a collection without pileup hits for (auto hit : ecal_hits ) { - - int *foundIndex = std::find(std::begin(hitIDs), std::end(hitIDs), hit.id()); + + int *foundIndex = std::find(std::begin(hitIDs), std::end(hitIDs), hit.getID()); // When the element is not found, std::find returns the end of the range if (foundIndex == std::end(hitIDs)) { //hit not found in the pileup cluster output_hits.emplace_back(hit); //keep it - } - }// if trk/ecal matched + } + }// over hits + }// if trk/ecal matched }// over PF objects - - } - + +} void PileupFinder::onProcessEnd() { ldmx_log(debug) << "Process ends!"; delete eCorr_; From 1d739da485dcd114f2ed6020b4eb71b6f54499d9 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 15:29:05 +0200 Subject: [PATCH 05/47] Fix more compilation problems --- Recon/src/Recon/PileupFinder.cxx | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index 5ab4f16f54..c4468e9435 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -52,23 +52,23 @@ void PileupFinder::produce(framework::Event& event) { for (const auto& pf_cand : pf_cands){ if(pf_cand.getPID()==3 || pf_cand.getPID()==7){ - double mom_vec = pf_cand.getTrackPxPyPz(); + std::vector mom_vec = pf_cand.getTrackPxPyPz(); float mom = mom_vec[0]*mom_vec[0]+mom_vec[1]*mom_vec[1]+mom_vec[2]*mom_vec[2]; mom = sqrt(mom); if (mom < min_mom_) continue; - ldmx_log(trace) << "Got pileup candidate with PID = " << pf_cand.getPID() << " and momentum = " << mom << " MeV." - - // now! use the hit-candidate association to get the associated ecal hits. - - - // TODO: - // write a header for this file --> DONE - // clean up all the methods below here --> DONE - // compile, possibly incrememnt particleflow candidate class imp nb - // try out hit association (print list?) - int pf_cl_idx = pf_cand.getEcalIndex(); + ldmx_log(trace) << "Got pileup candidate with PID = " << pf_cand.getPID() << " and momentum = " << mom << " MeV." ; + + // now! use the hit-candidate association to get the associated ecal hits. + + + // TODO: + // write a header for this file --> DONE + // clean up all the methods below here --> DONE + // compile, possibly incrememnt particleflow candidate class imp nb + // try out hit association (print list?) + int pf_cl_idx = pf_cand.getEcalIndex(); auto cl = clusters[pf_cl_idx]; auto hitIDs = cl.getHitIDs(); // make a collection without pileup hits @@ -86,8 +86,6 @@ void PileupFinder::produce(framework::Event& event) { } void PileupFinder::onProcessEnd() { ldmx_log(debug) << "Process ends!"; - delete eCorr_; - delete hCorr_; return; } From efdae1a3d6f40f6c55b9af11a9899b1742b1e637 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 15:30:53 +0200 Subject: [PATCH 06/47] Fix more compilation problems --- Recon/src/Recon/PileupFinder.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index c4468e9435..0e93ee7eb5 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -74,7 +74,7 @@ void PileupFinder::produce(framework::Event& event) { // make a collection without pileup hits for (auto hit : ecal_hits ) { - int *foundIndex = std::find(std::begin(hitIDs), std::end(hitIDs), hit.getID()); + auto foundIndex = std::find(std::begin(hitIDs), std::end(hitIDs), hit.getID()); // When the element is not found, std::find returns the end of the range if (foundIndex == std::end(hitIDs)) { //hit not found in the pileup cluster output_hits.emplace_back(hit); //keep it From 1148c14a785139193c1224da3b4a88095026dc13 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 16:14:39 +0200 Subject: [PATCH 07/47] Make idx getters const --- Recon/include/Recon/Event/PFCandidate.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Recon/include/Recon/Event/PFCandidate.h b/Recon/include/Recon/Event/PFCandidate.h index 96a5f67285..4d9a8d9da3 100644 --- a/Recon/include/Recon/Event/PFCandidate.h +++ b/Recon/include/Recon/Event/PFCandidate.h @@ -141,9 +141,9 @@ class PFCandidate { return {posHcalX_, posHcalY_, posHcalZ_}; } // associate component indices to the pf candidate - int getTrackIndex() { return track_idx_; } - int getEcalIndex() { return ecal_idx_ ; } - int getHcalIndex() { return hcal_idx_ ; } + int getTrackIndex() const { return track_idx_; } + int getEcalIndex() const { return ecal_idx_ ; } + int getHcalIndex() const { return hcal_idx_ ; } std::vector getTrackPxPyPz() const { return {trackPx_, trackPy_, trackPz_}; From beb87aab2d32d4a132bea026e0923eb62b587dc3 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 16:25:47 +0200 Subject: [PATCH 08/47] Add python config --- Recon/python/pileupFinder.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 Recon/python/pileupFinder.py diff --git a/Recon/python/pileupFinder.py b/Recon/python/pileupFinder.py new file mode 100644 index 0000000000..db4e8e9318 --- /dev/null +++ b/Recon/python/pileupFinder.py @@ -0,0 +1,23 @@ +"""Configuration for pfReco producers + +Sets all parameters to reasonable defaults. + +Examples +-------- + from LDMX.Recon.pfReco import pileupFinder + p.sequence.append( pileupFinder ) +""" + +from LDMX.Framework import ldmxcfg + +class pileupFinder(ldmxcfg.Producer) : + """Configuration for pileup finding from particle flow objects""" + def __init__(self, name='PileupFinder') : + super().__init__(name, 'recon::PileupFinder','Recon') + self.rec_hit_coll_name = 'EcalRecHits' + self.rec_hit_pass_name = '' + self.cluster_coll_name = 'PFEcalClusters' + self.pf_cand_coll_name = 'PFCandidates' + self.pf_cand_pass_name = '' + self.output_rec_hit_coll_name = 'EcalRecHitsNoPileup' + self.min_mom = 4000. From 1d30f7262dcf81715a2081f52218d3fa6b29a3b4 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 19:43:25 +0200 Subject: [PATCH 09/47] Add missing parameter in python config --- Recon/python/pileupFinder.py | 1 + 1 file changed, 1 insertion(+) diff --git a/Recon/python/pileupFinder.py b/Recon/python/pileupFinder.py index db4e8e9318..b125dc3315 100644 --- a/Recon/python/pileupFinder.py +++ b/Recon/python/pileupFinder.py @@ -17,6 +17,7 @@ def __init__(self, name='PileupFinder') : self.rec_hit_coll_name = 'EcalRecHits' self.rec_hit_pass_name = '' self.cluster_coll_name = 'PFEcalClusters' + self.cluster_pass_name = '' self.pf_cand_coll_name = 'PFCandidates' self.pf_cand_pass_name = '' self.output_rec_hit_coll_name = 'EcalRecHitsNoPileup' From 682036be7e2869026c8ffb517a7b4b9c680ab985 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 19:49:58 +0200 Subject: [PATCH 10/47] Another missing parameter in python config --- Recon/python/pileupFinder.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Recon/python/pileupFinder.py b/Recon/python/pileupFinder.py index b125dc3315..e9d28dc873 100644 --- a/Recon/python/pileupFinder.py +++ b/Recon/python/pileupFinder.py @@ -21,4 +21,4 @@ def __init__(self, name='PileupFinder') : self.pf_cand_coll_name = 'PFCandidates' self.pf_cand_pass_name = '' self.output_rec_hit_coll_name = 'EcalRecHitsNoPileup' - self.min_mom = 4000. + self.min_momentum = 4000. From 780721bbe5cb152fd48727f71d63ea3e12ebb917 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 19:57:13 +0200 Subject: [PATCH 11/47] Fix python double vs float nuisance --- Recon/include/Recon/PileupFinder.h | 2 +- Recon/src/Recon/PileupFinder.cxx | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Recon/include/Recon/PileupFinder.h b/Recon/include/Recon/PileupFinder.h index 120b1dada4..745df7fdc0 100644 --- a/Recon/include/Recon/PileupFinder.h +++ b/Recon/include/Recon/PileupFinder.h @@ -49,7 +49,7 @@ class PileupFinder : public framework::Producer { std::string output_rec_hit_coll_name_; // configuration - float min_mom_{0.}; //MeV + double min_mom_{0.}; //MeV }; } // namespace recon diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index 0e93ee7eb5..7984f8f53b 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -14,7 +14,7 @@ void PileupFinder::configure(framework::config::Parameters& ps) { cluster_pass_name_ = ps.getParameter("cluster_pass_name"); output_rec_hit_coll_name_ = ps.getParameter("output_rec_hit_coll_name"); // Algorithm configuration - min_mom_ = ps.getParameter("min_momentum"); + min_mom_ = ps.getParameter("min_momentum"); } // get pileup candidates from PFlow and make a cleaned-up hit collection From 7361e40d97d3268c9cb248d9b807fceadd0d9d88 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 20:57:40 +0200 Subject: [PATCH 12/47] Add output collection to the bus --- Recon/src/Recon/PileupFinder.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index 7984f8f53b..b086d0a3c2 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -51,7 +51,7 @@ void PileupFinder::produce(framework::Event& event) { for (const auto& pf_cand : pf_cands){ if(pf_cand.getPID()==3 || pf_cand.getPID()==7){ - + //we have both ecal cluster and track std::vector mom_vec = pf_cand.getTrackPxPyPz(); float mom = mom_vec[0]*mom_vec[0]+mom_vec[1]*mom_vec[1]+mom_vec[2]*mom_vec[2]; mom = sqrt(mom); @@ -83,6 +83,8 @@ void PileupFinder::produce(framework::Event& event) { }// if trk/ecal matched }// over PF objects + event.add(output_rec_hit_coll_name_, output_hits); + } void PileupFinder::onProcessEnd() { ldmx_log(debug) << "Process ends!"; From 3e6edf240ce643f4c3e529e3337e1954866ddab7 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Wed, 9 Jul 2025 11:39:42 +0200 Subject: [PATCH 13/47] Increment class def nb --- Recon/include/Recon/Event/PFCandidate.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Recon/include/Recon/Event/PFCandidate.h b/Recon/include/Recon/Event/PFCandidate.h index 4d9a8d9da3..17206e56fb 100644 --- a/Recon/include/Recon/Event/PFCandidate.h +++ b/Recon/include/Recon/Event/PFCandidate.h @@ -269,7 +269,7 @@ class PFCandidate { int hcal_idx_{-1}; /* The ROOT class definition. */ - ClassDef(PFCandidate, 1); + ClassDef(PFCandidate, 2); }; } // namespace ldmx From e982cc219d437ce7441ac37eca9fbaf3e0bca8a0 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Wed, 9 Jul 2025 11:40:23 +0200 Subject: [PATCH 14/47] Check that index was set before doing lookup --- Recon/src/Recon/PileupFinder.cxx | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index b086d0a3c2..d93802b38a 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -66,9 +66,12 @@ void PileupFinder::produce(framework::Event& event) { // TODO: // write a header for this file --> DONE // clean up all the methods below here --> DONE - // compile, possibly incrememnt particleflow candidate class imp nb + // compile, possibly increment particleflow candidate class imp nb --> DONE // try out hit association (print list?) int pf_cl_idx = pf_cand.getEcalIndex(); + ldmx_log(trace) << "Got Ecal cluster with index " << pf_cl_idx << " while cluster array length is " << clusters.size(); + if (pf_cl_idx < 0 ) // was never set + continue; auto cl = clusters[pf_cl_idx]; auto hitIDs = cl.getHitIDs(); // make a collection without pileup hits @@ -78,6 +81,7 @@ void PileupFinder::produce(framework::Event& event) { // When the element is not found, std::find returns the end of the range if (foundIndex == std::end(hitIDs)) { //hit not found in the pileup cluster output_hits.emplace_back(hit); //keep it + ldmx_log(trace) << "Got no-pileup hit! "; } }// over hits }// if trk/ecal matched From d5f93c6366736ee28ade3a3ca5e498ad03cf52a5 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Wed, 9 Jul 2025 11:56:52 +0200 Subject: [PATCH 15/47] Add more cluster index setting --- Recon/src/Recon/ParticleFlow.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx index 8a77a4f9f1..dc1dd40cfb 100644 --- a/Recon/src/Recon/ParticleFlow.cxx +++ b/Recon/src/Recon/ParticleFlow.cxx @@ -275,6 +275,7 @@ void ParticleFlow::produce(framework::Event& event) { EMIsTkLinked[em_idx] = true; tkIsEMLinked[i] = true; tkEMPairs[i] = em_idx; + cand.setEcalIndex(em_idx); break; } } From 21097d5c518d5b8f0165ff2ab0a1952dd5b23137 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Wed, 9 Jul 2025 12:01:39 +0200 Subject: [PATCH 16/47] Improve cluster index setting --- Recon/src/Recon/ParticleFlow.cxx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx index dc1dd40cfb..89e901286a 100644 --- a/Recon/src/Recon/ParticleFlow.cxx +++ b/Recon/src/Recon/ParticleFlow.cxx @@ -275,7 +275,6 @@ void ParticleFlow::produce(framework::Event& event) { EMIsTkLinked[em_idx] = true; tkIsEMLinked[i] = true; tkEMPairs[i] = em_idx; - cand.setEcalIndex(em_idx); break; } } @@ -318,6 +317,7 @@ void ParticleFlow::produce(framework::Event& event) { // chargedUnmatch.push_back(cand); } else { // if track is linked with ECal cluster fillCandEMCalo(cand, ecalClusters[tkEMPairs[i]]); + cand.setEcalIndex(tkEMPairs[i]); if (EMIsHadLinked[tkEMPairs[i]]) { // if ECal is linked with HCal // cluster fillCandHadCalo(cand, hcalClusters[EMHadPairs[tkEMPairs[i]]]); @@ -332,7 +332,6 @@ void ParticleFlow::produce(framework::Event& event) { for (int i = 0; i < ecalClusters.size(); i++) { // already linked with ECal in the previous step if (EMIsTkLinked[i]) continue; - ldmx::PFCandidate cand; fillCandEMCalo(cand, ecalClusters[i]); cand.setEcalIndex(i); From 26a8f32e683112bc5805deeaf1d59fbb31c7a5e0 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Thu, 17 Jul 2025 13:41:27 +0200 Subject: [PATCH 17/47] Cleanup naming conventions and verbosity --- .../include/TrigScint/TruthHitProducer.h | 16 +++--- TrigScint/python/truthHits.py | 5 +- TrigScint/src/TrigScint/TruthHitProducer.cxx | 54 ++++++++----------- 3 files changed, 31 insertions(+), 44 deletions(-) diff --git a/TrigScint/include/TrigScint/TruthHitProducer.h b/TrigScint/include/TrigScint/TruthHitProducer.h index 7a9aa28261..50d2035ca1 100644 --- a/TrigScint/include/TrigScint/TruthHitProducer.h +++ b/TrigScint/include/TrigScint/TruthHitProducer.h @@ -52,24 +52,22 @@ class TruthHitProducer : public framework::Producer { */ void produce(framework::Event &event) override; - /// Class to set the verbosity level. - // TODO: Make use of the global verbose parameter. - bool verbose_{false}; + private: /// Name of the input collection containing the sim hits - std::string inputCollection_; + std::string input_collection_; /// Name of the pass that the input collection is on (empty string means take /// any pass) - std::string inputPassName_; + std::string input_pass_name_; + + /// Name of the pass that the input simparticles collection is on + std::string sim_particles_pass_name_; /// Name of the output collection that will be used to store the /// selected sim hits - std::string outputCollection_; + std::string output_collection_; - private: - std::string sim_particles_passname_; - std::string input_collection_events_passname_; }; // TruthHitProducer diff --git a/TrigScint/python/truthHits.py b/TrigScint/python/truthHits.py index e9d9023493..933f9dd355 100644 --- a/TrigScint/python/truthHits.py +++ b/TrigScint/python/truthHits.py @@ -20,10 +20,7 @@ def __init__(self,name) : self.input_collection="TriggerPadUpSimHits" self.input_pass_name="" #take any pass self.output_collection="truthBeamElectronsUp" - self.verbose = False - - self.sim_particles_passname = "" - self.input_collection_events_passname = "" + self.sim_particles_pass_name = "" def up() : """Get the beam electron truth hits for the trigger pad upstream of target""" diff --git a/TrigScint/src/TrigScint/TruthHitProducer.cxx b/TrigScint/src/TrigScint/TruthHitProducer.cxx index 5120a701f4..0689aac192 100644 --- a/TrigScint/src/TrigScint/TruthHitProducer.cxx +++ b/TrigScint/src/TrigScint/TruthHitProducer.cxx @@ -8,38 +8,31 @@ TruthHitProducer::TruthHitProducer(const std::string &name, : Producer(name, process) {} void TruthHitProducer::configure(framework::config::Parameters ¶meters) { - inputCollection_ = parameters.getParameter("input_collection"); - inputPassName_ = parameters.getParameter("input_pass_name"); - outputCollection_ = parameters.getParameter("output_collection"); - sim_particles_passname_ = - parameters.getParameter("sim_particles_passname"); - input_collection_events_passname_ = - parameters.getParameter("input_collection_events_passname"); + input_collection_ = parameters.getParameter("input_collection"); + input_pass_name_ = parameters.getParameter("input_pass_name"); + output_collection_ = parameters.getParameter("output_collection"); + sim_particles_pass_name_ = + parameters.getParameter("sim_particles_pass_name"); - verbose_ = parameters.getParameter("verbose"); - - if (verbose_) { ldmx_log(info) << "In TruthHitProducer: configure done!"; ldmx_log(info) << "Got parameters: " << "\nInput collection: " - << inputCollection_ - << "\nInput pass name: " << inputPassName_ - << "\nOutput collection: " << outputCollection_ - << "\nVerbose: " << verbose_; - } + << input_collection_ + << "\nInput pass name: " << input_pass_name_ + << "\nOutput collection: " << output_collection_; } void TruthHitProducer::produce(framework::Event &event) { // Check if the collection exists. If not, don't bother processing the event. - if (!event.exists(inputCollection_, input_collection_events_passname_)) { - ldmx_log(error) << "No input collection called " << inputCollection_ + if (!event.exists(input_collection_, input_pass_name_)) { + ldmx_log(error) << "No input collection called " << input_collection_ << " found; skipping!"; return; } // looper over sim hits and aggregate energy depositions for each detID const auto simHits{event.getCollection( - inputCollection_, inputPassName_)}; + input_collection_, input_pass_name_)}; auto particleMap{event.getMap( - "SimParticles", sim_particles_passname_)}; + "SimParticles", sim_particles_pass_name_)}; std::vector truthBeamElectrons; @@ -49,15 +42,14 @@ void TruthHitProducer::produce(framework::Event &event) { // check if hit is from beam electron and, if so, add to output collection for (int i = 0; i < simHit.getNumberOfContribs(); i++) { auto contrib = simHit.getContrib(i); - if (verbose_) { - ldmx_log(debug) << "contrib " << i << " trackID: " << contrib.trackID - << " pdgID: " << contrib.pdgCode - << " edep: " << contrib.edep; - ldmx_log(debug) << "\t particle id: " - << particleMap[contrib.trackID].getPdgID() - << " particle status: " - << particleMap[contrib.trackID].getGenStatus(); - } + ldmx_log(trace) << "contrib " << i << " trackID: " << contrib.trackID + << " pdgID: " << contrib.pdgCode + << " edep: " << contrib.edep; + ldmx_log(trace) << "\t particle id: " + << particleMap[contrib.trackID].getPdgID() + << " particle status: " + << particleMap[contrib.trackID].getGenStatus(); + // if the trackID is in the map if (particleMap.find(contrib.trackID) != particleMap.end()) { // beam electron (PDGID = 11, genStatus == 1) @@ -67,9 +59,9 @@ void TruthHitProducer::produce(framework::Event &event) { } } if (keep) truthBeamElectrons.push_back(simHit); - } - } - event.add(outputCollection_, truthBeamElectrons); + }//over simhit contribs + }//over simhits + event.add(output_collection_, truthBeamElectrons); } } // namespace trigscint From dcc6075d9896da56c32bd23e76bbf872d8da0326 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Wed, 17 Sep 2025 09:43:33 +0200 Subject: [PATCH 18/47] Add analysis histograms --- DQM/include/DQM/EcalClusterAnalyzer.h | 9 ++- DQM/python/dqm.py | 20 +++-- DQM/src/DQM/EcalClusterAnalyzer.cxx | 111 ++++++++++++++++++++------ 3 files changed, 111 insertions(+), 29 deletions(-) diff --git a/DQM/include/DQM/EcalClusterAnalyzer.h b/DQM/include/DQM/EcalClusterAnalyzer.h index a7fab81738..e22948c23c 100644 --- a/DQM/include/DQM/EcalClusterAnalyzer.h +++ b/DQM/include/DQM/EcalClusterAnalyzer.h @@ -46,7 +46,14 @@ 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 diff --git a/DQM/python/dqm.py b/DQM/python/dqm.py index f51920247b..160901aa65 100644 --- a/DQM/python/dqm.py +++ b/DQM/python/dqm.py @@ -876,23 +876,33 @@ 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 # 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("mixed_hit_energy", f"Percentage of total energy coming from hits with at least {100*self.mixed_hit_cutoff}% energy contributions from additional electron(s)", 21, 0, 105) + 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.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_hits", "Total energy (edep)", 30, 0, 150, "Hits in cluster", 30, 0, 300) 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) + 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(), diff --git a/DQM/src/DQM/EcalClusterAnalyzer.cxx b/DQM/src/DQM/EcalClusterAnalyzer.cxx index 5c73911ae9..d303a00fc8 100644 --- a/DQM/src/DQM/EcalClusterAnalyzer.cxx +++ b/DQM/src/DQM/EcalClusterAnalyzer.cxx @@ -30,8 +30,12 @@ void EcalClusterAnalyzer::configure(framework::config::Parameters& ps) { cluster_coll_name_ = ps.getParameter("cluster_coll_name"); cluster_pass_name_ = ps.getParameter("cluster_pass_name"); - ecal_sp_hits_passname_ = - ps.getParameter("ecal_sp_hits_passname"); + ecal_sp_hits_coll_name_ = + ps.getParameter("ecal_sp_hits_coll_name"); + ecal_sp_hits_pass_name_ = + ps.getParameter("ecal_sp_hits_pass_name"); + mixed_hit_cutoff_ = ps.getParameter("mixed_hit_cutoff"); + return; } @@ -53,8 +57,11 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { std::unordered_map>> hitInfo; hitInfo.reserve(ecal_rec_hits.size()); - double dist; - if (nbr_of_electrons_ == 2) { + double dist = -1; + std::map true_energy; + std::map delta_energy; + + // if (nbr_of_electrons_ == 2) { // Measures distance between two electrons in the ECal scoring plane // TODO: generalize for n electrons std::vector pos1; @@ -63,59 +70,102 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { bool p2 = false; const auto& ecal_sp_hits{event.getCollection( - "EcalScoringPlaneHits", ecal_sp_hits_passname_)}; + ecal_sp_hits_coll_name_, ecal_sp_hits_pass_name_)}; for (const ldmx::SimTrackerHit& spHit : ecal_sp_hits) { - if (spHit.getTrackID() == 1) { + //here could require pdgID = 11, but that would mean always only storing electron energies, along with their incoming separation + true_energy[spHit.getTrackID()] = spHit.getEnergy(); + ldmx_log(trace) << "For TrackID " << spHit.getTrackID() << ", got SP energy " << spHit.getEnergy(); + // TODO: thistrackID assumption breaks with overlay. use origin id? use parents and PDGID + if (spHit.getTrackID() == 1) { pos1 = spHit.getPosition(); p1 = true; + } else if (spHit.getTrackID() == 2) { pos2 = spHit.getPosition(); p2 = true; } } - if (p1 && p2) + if (p1 && p2) { dist = std::sqrt(std::pow((pos1[0] - pos2[0]), 2) + std::pow((pos1[1] - pos2[1]), 2)); - } - + histograms_.fill("SP_distance", dist); + ldmx_log(debug) << "Found two electrons! Distance = " << dist; + } + //} + double totEventEnergy=0; + std::vector totOriginEdep; + totOriginEdep.resize(nbr_of_electrons_ + 1); + int nMixed = 0; 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 + ldmx_log(trace) << "\tFound simhit matching rechit with ID" << hit.getID(); int ancestor = 0; int prevAncestor = 0; bool tagged = false; int tag = 0; std::vector edep; 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& c = it->getContrib(i); // get origin electron ID ancestor = c.originID; + ldmx_log(trace) << "\t\t\tAncestor ID " << ancestor << " with edep " << c.edep; + totEventEnergy+=c.edep; // store energy from this contrib at index = origin electron ID - if (ancestor <= nbr_of_electrons_) edep[ancestor] += c.edep; + if (ancestor <= nbr_of_electrons_) { + edep[ancestor] += c.edep; + totOriginEdep[ancestor] += c.edep; + eTot+= 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; + ldmx_log(trace) << "\t\t\t\tMixed hit! Ancestor ID changed to " << ancestor; } prevAncestor = 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 = prevAncestor; + // if not tagged, hit was from a single electron (within acceptable purity) + tag = ancestor; //prevAncestor; } + else nMixed++; histograms_.fill("ancestors", tag); hitInfo.insert({hit.getID(), std::make_pair(tag, edep)}); } } int clusteredHits = 0; + histograms_.fill("tag0frac_vs_SPdist",dist, (float)nMixed/ecal_rec_hits.size()); + ldmx_log(debug) << "Got " << nMixed << " mixed hits, a fraction of " << (float)nMixed/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; + } + ldmx_log(trace) << "Looping over ecal clusters"; for (const auto& cl : ecal_clusters) { // for each cluster // total number of hits coming from electron, index = electron ID @@ -126,7 +176,7 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { e.resize(nbr_of_electrons_ + 1); double eSum = 0.; double nSum = 0.; - + ldmx_log(trace) << "Looping over hits in the cluster"; const auto& hitIDs = cl.getHitIDs(); for (const auto& id : hitIDs) { // for each hit in cluster, find previously stored info @@ -155,11 +205,15 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { clusteredHits++; } - } - + }//over hits + + // get largest energy contribution + double eMax = *max_element(e.begin(), e.end()); + std::string toLog; + for (auto nb : e) toLog.append(std::to_string(nb)+" "); + ldmx_log(debug) << "Energies vector is " << toLog ; + 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)); @@ -174,12 +228,23 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { double nMax = *max_element(n.begin(), n.end()); histograms_.fill("same_ancestor", 100. * (nMax / nSum)); } - } - - histograms_.fill("clusterless_hits", (ecal_rec_hits.size() - clusteredHits)); + // find the main contributor + auto elt = distance(e.begin(), max_element(e.begin(), e.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()); + histograms_.fill("cluster_RMSY", cl.getRMSY()); + }// over clusters + std::string moreLog; + for (auto nb : totOriginEdep) moreLog.append(std::to_string(nb)+" "); + ldmx_log(debug) << "Edep per ancestor in event is " << moreLog ; + ldmx_log(debug) << "Total energy deposited in event: " << totEventEnergy ; + + histograms_.fill("dE_cl2_vs_cl1", delta_energy[1], delta_energy[2]); + histograms_.fill("unclustered_hits", (ecal_rec_hits.size() - clusteredHits)); histograms_.fill("total_rechits_in_event", ecal_rec_hits.size()); - histograms_.fill( - "clusterless_hits_percentage", + histograms_.fill("unclustered_hits_percentage", 100. * (ecal_rec_hits.size() - clusteredHits) / ecal_rec_hits.size()); } From 7117e56467447f7bff85184b2e6d610c4d17734f Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Wed, 17 Sep 2025 09:44:58 +0200 Subject: [PATCH 19/47] Add pu finding improvements --- Ecal/include/Ecal/CLUE.h | 5 +- Ecal/src/Ecal/CLUE.cxx | 179 ++++++++++--------- Recon/src/Recon/PileupFinder.cxx | 32 ++-- TrigScint/src/TrigScint/TruthHitProducer.cxx | 6 + 4 files changed, 121 insertions(+), 101 deletions(-) diff --git a/Ecal/include/Ecal/CLUE.h b/Ecal/include/Ecal/CLUE.h index 5700d8aa5e..b64387d7e2 100644 --- a/Ecal/include/Ecal/CLUE.h +++ b/Ecal/include/Ecal/CLUE.h @@ -62,7 +62,7 @@ class CLUE { hits = {}; } }; - + /* helper functions to calculate distances */ 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); @@ -72,6 +72,9 @@ class CLUE { std::vector> setup( const std::vector& hits); + // get distance between clusters in the first layers, proxy for electron sep. + void electronSeparation(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 diff --git a/Ecal/src/Ecal/CLUE.cxx b/Ecal/src/Ecal/CLUE.cxx index bd4e00c19c..55de87740d 100644 --- a/Ecal/src/Ecal/CLUE.cxx +++ b/Ecal/src/Ecal/CLUE.cxx @@ -21,64 +21,68 @@ float CLUE::floatDist(float x1, float y1, float z1, float x2, float y2, 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; -// } -// } -// } + void CLUE::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.; + // sort hits in z + 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(IntermediateCluster(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) { @@ -93,6 +97,7 @@ std::vector> CLUE::createLayers( layers.push_back({}); double layerSeparation = (maxZ - layerZ) / nbr_of_layers_; ldmx_log(trace) << " Layer separation: " << layerSeparation; + ldmx_log(trace) << " Number of layers: " << nbr_of_layers_; ldmx_log(trace) << " Creating layer 0"; double highestEnergy = 0.; @@ -143,7 +148,8 @@ std::vector> CLUE::setup( // collapse z dimension float x = roundToDecimal(hit->getXPos(), 4); float y = roundToDecimal(hit->getYPos(), 4); - ldmx_log(trace) << " New hit { x: " << x << " y: " << y << "}"; + float z = roundToDecimal(hit->getZPos(), 4); + ldmx_log(trace) << " New hit { x: " << x << " y: " << y << "}" << " (and z: " << z << ")"; std::pair coords; if (dc_ != 0 && nbr_of_layers_ > 1) { // if more than one layer, divide hits into densities with side dc_ @@ -158,7 +164,7 @@ std::vector> CLUE::setup( j = -j; y = (j + 0.5) * dc_; } else // TODO i think this 1 should be a j - y = (1 - 0.5) * dc_; // set x,y to middle of box + y = (j - 0.5) * dc_; // set x,y to middle of box coords = {i, j}; ldmx_log(trace) << " Index " << i << ", " << j << "; x: " << x << " y: " << y; @@ -169,9 +175,9 @@ std::vector> CLUE::setup( } if (densityMap.find(coords) == densityMap.end()) { densityMap.emplace(coords, std::make_shared(x, y)); - ldmx_log(trace) << " New density created"; + ldmx_log(trace) << " * New density created"; } else { - ldmx_log(trace) << " Found density with x: " << densityMap[coords]->x + ldmx_log(trace) << " --> Found density with x: " << densityMap[coords]->x << " y: " << densityMap[coords]->y; } densityMap[coords]->hits.push_back(hit); @@ -212,7 +218,7 @@ std::vector> CLUE::setup( densities[i]->delta = d; densities[i]->follower_of = j; ldmx_log(trace) << " New parent, index " << j - << "; delta: " << std::setprecision(20) << d; + << "; delta: " << std::setprecision(10) << d; } } } @@ -312,7 +318,7 @@ std::vector> CLUE::clustering( densities[i]->total_energy > rhoc_ && densities[i]->delta > deltac_; if (isSeed) { - ldmx_log(trace) << " Distance to centroid: " + ldmx_log(trace) << " Distance to event centroid: " << floatDist(densities[i]->x, densities[i]->y, event_centroid_.centroid().Px(), event_centroid_.centroid().Py()); @@ -334,6 +340,7 @@ std::vector> CLUE::clustering( ldmx_log(trace) << " Follower"; int& parentIndex = densities[i]->follower_of; if (parentIndex != -1) followers[parentIndex].push_back(i); + else ldmx_log(error) << " Somehow found a follower with parent index -1: id = " << i; } else { ldmx_log(trace) << " Outlier"; } @@ -413,33 +420,37 @@ std::vector> CLUE::layerSetup() { 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; + << ", energy: " << currentLayer[i]->total_energy + << " position (x,y)= {" << currentLayer[i]->x << "," + << currentLayer[i]->y << ")"; int depth = 1; // decide delta and followerof from seeds in previous and next layer // do { // depth++; if (layer - depth >= 0) { // look at previous layer + ldmx_log(trace) << " In previous layer... "; auto& previousLayer = seeds_[layer - depth]; + ldmx_log(trace) << " Got " << previousLayer.size() << " seeds"; 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); + auto dz = std::abs(currentLayer[i]->z - previousLayer[j]->z); ldmx_log(trace) << " Delta to index " << previousLayer[j]->index - << ": " << std::setprecision(20) << d; + << ": " << std::setprecision(10) << d; ldmx_log(trace) << " DeltaZ to index " << previousLayer[j]->index - << ": " << std::setprecision(20) << dz << std::endl; + << ": " << std::setprecision(10) << 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 + << " * New parent: index " << previousLayer[j]->index << " on layer " << layer - depth << "; energy " << previousLayer[j]->total_energy; ldmx_log(trace) - << " New delta: " << std::setprecision(20) << d; + << " * New delta: " << std::setprecision(10) << d; ldmx_log(trace) - << " New deltaZ: " << std::setprecision(20) << dz; + << " * New deltaZ: " << std::setprecision(10) << dz; currentLayer[i]->delta = d; currentLayer[i]->z_delta = dz; currentLayer[i]->follower_of = previousLayer[j]->index; @@ -450,24 +461,26 @@ std::vector> CLUE::layerSetup() { } } if (layer + depth < nbr_of_layers_) { + ldmx_log(trace) << " In next layer: " << layer + depth; auto& nextLayer = seeds_[layer + depth]; + ldmx_log(trace) << " Got " << nextLayer.size() << " seeds"; 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); + auto dz = std::abs(currentLayer[i]->z - nextLayer[j]->z); ldmx_log(trace) << " Delta to index " << nextLayer[j]->index - << ": " << std::setprecision(20) << d; + << ": " << std::setprecision(10) << d; ldmx_log(trace) << " DeltaZ to index " << nextLayer[j]->index - << ": " << std::setprecision(20) << dz; + << ": " << std::setprecision(10) << 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 + 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; + << " * New delta: " << std::setprecision(10) << d; ldmx_log(trace) - << " New deltaZ: " << std::setprecision(20) << dz; + << " * New deltaZ: " << std::setprecision(10) << dz; currentLayer[i]->delta = d; currentLayer[i]->z_delta = dz; currentLayer[i]->follower_of = nextLayer[j]->index; @@ -475,12 +488,14 @@ std::vector> CLUE::layerSetup() { currentLayer[i]->total_energy) { break; } - } - } + }//over seeds in next layer + }// if not outside layer range + // } while (currentLayer[i]->layerFollowerOf == -1 && (layer - depth >= // 0 || layer + depth < nbrOfLayers_)); + ldmx_log(trace) << " Done setting parents."; densities.push_back(currentLayer[i]); - } + }//over seeds in current layer layer_rho_c_.push_back(highestEnergy / 2); } return densities; @@ -540,14 +555,14 @@ void CLUE::cluster(const std::vector& unsorted_hits, double dc, return a->getZPos() < b->getZPos(); }); - seeds_.resize(nbrOfLayers); + seeds_.resize(nbr_of_layers_); 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 + // convertToIntermediateClusters(clusters); // uncomment for layer //LK try uncommenting // clustering without 3D } // Below for CLUE3D, comment for just layer clustering diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index d93802b38a..34219c6345 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -48,7 +48,10 @@ void PileupFinder::produce(framework::Event& event) { // get the list of hits associated with pileup candidates // if a rechit is not on that list, add to output collection. std::vector output_hits; + std::vector pileup_hitIDs; + //this needs to be a two-step procedure: loop over all clusters deemed to be pileup to find all their associated hits + // then loop over that list to make a collection of output hits that doesn't contain them for (const auto& pf_cand : pf_cands){ if(pf_cand.getPID()==3 || pf_cand.getPID()==7){ //we have both ecal cluster and track @@ -61,32 +64,25 @@ void PileupFinder::produce(framework::Event& event) { ldmx_log(trace) << "Got pileup candidate with PID = " << pf_cand.getPID() << " and momentum = " << mom << " MeV." ; // now! use the hit-candidate association to get the associated ecal hits. - - - // TODO: - // write a header for this file --> DONE - // clean up all the methods below here --> DONE - // compile, possibly increment particleflow candidate class imp nb --> DONE - // try out hit association (print list?) int pf_cl_idx = pf_cand.getEcalIndex(); ldmx_log(trace) << "Got Ecal cluster with index " << pf_cl_idx << " while cluster array length is " << clusters.size(); if (pf_cl_idx < 0 ) // was never set continue; auto cl = clusters[pf_cl_idx]; auto hitIDs = cl.getHitIDs(); - // make a collection without pileup hits - for (auto hit : ecal_hits ) { - - auto foundIndex = std::find(std::begin(hitIDs), std::end(hitIDs), hit.getID()); - // When the element is not found, std::find returns the end of the range - if (foundIndex == std::end(hitIDs)) { //hit not found in the pileup cluster - output_hits.emplace_back(hit); //keep it - ldmx_log(trace) << "Got no-pileup hit! "; - } - }// over hits + // add to collection of pileup hits + pileup_hitIDs.insert(pileup_hitIDs.end(), hitIDs.begin(), hitIDs.end()); }// if trk/ecal matched }// over PF objects - + + for (auto hit : ecal_hits ) { + auto foundIndex = std::find(std::begin(pileup_hitIDs), std::end(pileup_hitIDs), hit.getID()); + // When the element is not found, std::find returns the end of the range + if (foundIndex == std::end(pileup_hitIDs)) { //hit not part of any pileup cluster + output_hits.emplace_back(hit); //keep it + ldmx_log(trace) << "Got no-pileup hit! "; + } + } event.add(output_rec_hit_coll_name_, output_hits); } diff --git a/TrigScint/src/TrigScint/TruthHitProducer.cxx b/TrigScint/src/TrigScint/TruthHitProducer.cxx index 0689aac192..2dab377ee4 100644 --- a/TrigScint/src/TrigScint/TruthHitProducer.cxx +++ b/TrigScint/src/TrigScint/TruthHitProducer.cxx @@ -25,9 +25,15 @@ void TruthHitProducer::produce(framework::Event &event) { // Check if the collection exists. If not, don't bother processing the event. if (!event.exists(input_collection_, input_pass_name_)) { ldmx_log(error) << "No input collection called " << input_collection_ + << "_" << input_pass_name_ << " found; skipping!"; + return; + } + if (!event.exists("SimParticles",sim_particles_pass_name_)) { + ldmx_log(error) << "No input SimParticle collection with pass " << sim_particles_pass_name_ << " found; skipping!"; return; } + // looper over sim hits and aggregate energy depositions for each detID const auto simHits{event.getCollection( input_collection_, input_pass_name_)}; From 7284573361e5afc8c674564d085f221ca51ddd11 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Thu, 18 Sep 2025 09:42:50 +0200 Subject: [PATCH 20/47] resolve comment merge conflict --- Recon/include/Recon/Event/CaloCluster.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Recon/include/Recon/Event/CaloCluster.h b/Recon/include/Recon/Event/CaloCluster.h index 3056c5becf..fa7aeca8a1 100644 --- a/Recon/include/Recon/Event/CaloCluster.h +++ b/Recon/include/Recon/Event/CaloCluster.h @@ -135,8 +135,8 @@ class CaloCluster { /// @brief Delta unc on unc in y-z plane double getEDYDZ() const { return err_dydz_; } - // get hit rawIDs (unused) - const std::vector& getHitIDs() const { return hit_ids_; } + // get hit rawIDs + const std::vector& getHitIDs() const { return hitIDs_; } // ability to store limited hit info const std::vector& getHitX() const { return hit_x_; } From f0b3c9bb74c5d36c56fc5bd670e9cc3fbdc6b82c Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Thu, 18 Sep 2025 09:47:50 +0200 Subject: [PATCH 21/47] resolve merge conflict --- Recon/include/Recon/Event/PFCandidate.h | 46 ++++++++++++++++++++++++- Recon/src/Recon/ParticleFlow.cxx | 3 ++ 2 files changed, 48 insertions(+), 1 deletion(-) diff --git a/Recon/include/Recon/Event/PFCandidate.h b/Recon/include/Recon/Event/PFCandidate.h index ce01f1871a..12539ea608 100644 --- a/Recon/include/Recon/Event/PFCandidate.h +++ b/Recon/include/Recon/Event/PFCandidate.h @@ -9,6 +9,9 @@ // ROOT #include "TObject.h" //For ClassDef +// ldmx-sw objects +//#include "Ecal/Event/EcalHit.h" +//#include "Hcal/Event/HcalHit.h" namespace ldmx { @@ -54,7 +57,11 @@ class PFCandidate { track_px_ = x_; track_py_ = y_; track_pz_ = z_; - } + + // associate component indices to the pf candidate + void setTrackIndex(int x) { track_idx_ = x; } + void setEcalIndex(int x) { ecal_idx_ = x; } + void setHcalIndex(int x) { hcal_idx_ = x; } void setEcalEnergy(float x_) { ecal_energy_ = x_; } void setEcalRawEnergy(float x_) { ecal_raw_energy_ = x_; } @@ -104,6 +111,19 @@ class PFCandidate { void setTruthEnergy(double x_) { truth_energy_ = x_; } void setTruthPdgId(int x_) { truth_pdg_id_ = x_; } + /** + * Take in the ecal hits that make up the candidate. + * @param hit The digi hit's entry number in the events digi + * collection. + */ + //void setEcalHits(const std::vector hits) { ecal_hits_ = hits; } + + /** + * Take in the hcal hits that make up the candidate. + * @param hit The digi hit's entry number in the events digi + * collection. + */ + // void setHcalHits(const std::vector hits) { hcal_hits_ = hits; } /* Getters */ @@ -120,6 +140,10 @@ class PFCandidate { std::vector getHcalPositionXYZ() const { return {pos_hcal_x_, pos_hcal_y_, pos_hcal_z_}; } + // associate component indices to the pf candidate + int getTrackIndex() { return track_idx_; } + int getEcalIndex() { return ecal_idx_ ; } + int getHcalIndex() { return hcal_idx_ ; } std::vector getTrackPxPyPz() const { return {track_px_, track_py_, track_pz_}; @@ -161,6 +185,21 @@ class PFCandidate { double getTruthEnergy() { return truth_energy_; } int getTruthPdgId() { return truth_pdg_id_; } + /** + * Take in the ecal hits that make up the candidate. + * @param hit The digi hit's entry number in the events digi + * collection. + */ + //std::vector getEcalHits() { return ecal_hits_; } + + /** + * Take in the hcal hits that make up the candidate. + * @param hit The digi hit's entry number in the events digi + * collection. + */ + // std::vector getHcalHits() { return hcal_hits_; } + + private: /* Particle ID enum */ int pid_{0}; @@ -224,6 +263,11 @@ class PFCandidate { double truth_energy_{0}; int truth_pdg_id_{0}; + /* Indices of the components making up the PFlow object */ + int track_idx_{-1}; + int ecal_idx_{-1}; + int hcal_idx_{-1}; + /* The ROOT class definition. */ ClassDef(PFCandidate, 2); }; diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx index aec3a2ea72..9ef8253e59 100644 --- a/Recon/src/Recon/ParticleFlow.cxx +++ b/Recon/src/Recon/ParticleFlow.cxx @@ -313,6 +313,7 @@ void ParticleFlow::produce(framework::Event& event) { ldmx::PFCandidate cand; fillCandTrack(cand, tracks[i]); // append track info to candidate + cand.setTrackIndex(i); if (!tk_is_em_linked[i]) { // chargedUnmatch.push_back(cand); } else { // if track is linked with ECal cluster @@ -334,6 +335,7 @@ void ParticleFlow::produce(framework::Event& event) { ldmx::PFCandidate cand; fillCandEMCalo(cand, ecal_clusters[i]); + cand.setEcalIndex(i); if (em_is_had_linked[tk_em_pairs[i]]) { fillCandHadCalo(cand, hcal_clusters[em_had_pairs[i]]); // emMatch.push_back(cand); @@ -347,6 +349,7 @@ void ParticleFlow::produce(framework::Event& event) { if (had_is_em_linked[i]) continue; ldmx::PFCandidate cand; fillCandHadCalo(cand, hcal_clusters[i]); + cand.setHcalIndex(i); // hadOnly.push_back(cand); pf_cands.push_back(cand); } From 1067ee2bafe727d29ef48c9fc0e3e8af39eae56f Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 15:16:03 +0200 Subject: [PATCH 22/47] Add first attempts of pileup hit finding --- Recon/include/Recon/PileupFinder.h | 55 +++++++++++++++++ Recon/src/Recon/PileupFinder.cxx | 97 ++++++++++++++++++++++++++++++ 2 files changed, 152 insertions(+) create mode 100644 Recon/include/Recon/PileupFinder.h create mode 100644 Recon/src/Recon/PileupFinder.cxx diff --git a/Recon/include/Recon/PileupFinder.h b/Recon/include/Recon/PileupFinder.h new file mode 100644 index 0000000000..657f22bd37 --- /dev/null +++ b/Recon/include/Recon/PileupFinder.h @@ -0,0 +1,55 @@ +/** + * @file PileupFinder.h + * @brief Simple PFlow algorithm + * @author Christian Herwig, Fermilab + */ + +#ifndef PARTICLEFLOW_H +#define PARTICLEFLOW_H + +// LDMX Framework +#include "Ecal/Event/EcalCluster.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 +#include "Hcal/Event/HcalCluster.h" +#include "Recon/Event/CaloCluster.h" +#include "Recon/Event/PFCandidate.h" +#include "SimCore/Event/SimParticle.h" +#include "SimCore/Event/SimTrackerHit.h" +#include "TGraph.h" + +namespace recon { + +/** + * @class PileupFinder + * @brief + */ +class PileupFinder : public framework::Producer { + public: + PileupFinder(const std::string& name, framework::Process& process) + : framework::Producer(name, process) {} + + virtual void configure(framework::config::Parameters& ps); + + virtual void produce(framework::Event& event); + + virtual void onProcessEnd(); + + private: + + // name of collections for PF input object to be passed + std::string rec_hit_coll_name_; + std::string rec_hit_pass_name_; + std::string pf_cand_coll_name_; + std::string pf_cand_pass_name_; + std::string cluster_coll_name_; + std::string cluster_pass_name_; + // name of collection for pileup-free output hit coll + std::string output_rec_hit_coll_name_; + // configuration + +}; +} // namespace recon + +#endif /* PILEUPFINDER_H */ diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx new file mode 100644 index 0000000000..554c162f64 --- /dev/null +++ b/Recon/src/Recon/PileupFinder.cxx @@ -0,0 +1,97 @@ +#include "Recon/PileupFinder.h" + +#include + +namespace recon { + +void PileupFinder::configure(framework::config::Parameters& ps) { + // I/O + rec_hit_coll_name_ = ps.getParameter("rec_hit_coll_name"); + rec_hit_pass_name_ = ps.getParameter("rec_hit_pass_name"); + pf_cand_coll_name_ = ps.getParameter("pf_cand_coll_name"); + pf_cand_pass_name_ = ps.getParameter("pf_cand_pass_name"); + cluster_coll_name_ = ps.getParameter("cluster_coll_name"); + cluster_pass_name_ = ps.getParameter("cluster_pass_name"); + output_rec_hit_coll_name_ = ps.getParameter("output_rec_hit_coll_name"); + // Algorithm configuration + min_mom_ = ps.getParameter("min_momentum"); +} + +// get pileup candidates from PFlow and make a cleaned-up hit collection +void PileupFinder::produce(framework::Event& event) { + if (!event.exists(rec_hit_coll_name_,rec_hit_pass_name_)) { //ecal rechits + ldmx_log(error) << "Unable to find (one) collection named " << rec_hit_coll_name << + "_" << rec_hit_pass_name; + return; + } + if (!event.exists(pf_cand_coll_name_, pf_cand_pass_name_ )) { + ldmx_log(error) << "Unable to find (one) collection named " << pf_cand_coll_name << + "_" << pf_cand_pass_name; + return; + } + if (!event.exists(cluster_coll_name_, cluster_pass_name_ )) { + ldmx_log(error) << "Unable to find (one) collection named " << cluster_coll_name << + "_" << cluster_pass_name; + return; + } + + const auto& ecal_hits{event.getCollection(rec_hit_coll_name_, + rec_hit_pass_name_)}; + + const auto& pf_cands{event.getCollection(pf_cand_coll_name_, + pf_cand_pass_name_)}; + + const auto& clusters{event.getCollection(cluster_coll_name_, + cluster_pass_name_)}; + // get PID 3 and 7 -- the ones with track and ecal matching + // get the high-momentum track ones from there -- pileup candidate! + // get the list of hits associated with pileup candidates + // if a rechit is not on that list, add to output collection. + std::vector output_hits; + + for (const auto& pf_cand : pfCandidates){ + if(pf_cand.getPID()==3 || pf_cand.getPID()==7){ + + double mom_vec = pf_cand.getTrackPxPyPz(); + float mom = mom_vec[0]*mom_vec[0]+mom_vec[1]*mom_vec[1]+mom_vec[2]*mom_vec[2]; + mom = sqrt(mom); + + if (mom < min_mom_) + continue; + ldmx_log(trace) << "Got pileup candidate with PID = " << pf_cand.getPID() << " and momentum = " << mom << " MeV." + + // now! use the hit-candidate association to get the associated ecal hits. + + + // TODO: + // write a header for this file --> DONE + // clean up all the methods below here --> DONE + // compile, possibly incrememnt particleflow candidate class imp nb + // try out hit association (print list?) + int pf_cl_idx = pf_cand.getEcalIndex(); + auto cl = clusters[pf_cl_idx]; + auto hitIds = cl.getHitIDs(); + // make a collection without pileup hits + for (auto hit : ecal_hits ) { + + int *foundIndex = std::find(std::begin(hitIDs), std::end(hitIDs), hit.id()); + // When the element is not found, std::find returns the end of the range + if (foundIndex == std::end(hitIDs)) { //hit not found in the pileup cluster + output_hits.emplace_back(hit); //keep it + } + }// if trk/ecal matched + }// over PF objects + + } + +void PileupFinder::onProcessEnd() { + ldmx_log(debug) << "Process ends!"; + delete eCorr_; + delete hCorr_; + + return; +} + +} // namespace recon + +DECLARE_PRODUCER(recon::PileupFinder); From 8fb45f9d41db74c87680ed76885f6c714b6b4dde Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 15:26:38 +0200 Subject: [PATCH 23/47] Fix compilation problems --- Recon/include/Recon/PileupFinder.h | 1 + Recon/src/Recon/PileupFinder.cxx | 30 +++++++++++++++--------------- 2 files changed, 16 insertions(+), 15 deletions(-) diff --git a/Recon/include/Recon/PileupFinder.h b/Recon/include/Recon/PileupFinder.h index 657f22bd37..120b1dada4 100644 --- a/Recon/include/Recon/PileupFinder.h +++ b/Recon/include/Recon/PileupFinder.h @@ -49,6 +49,7 @@ class PileupFinder : public framework::Producer { std::string output_rec_hit_coll_name_; // configuration + float min_mom_{0.}; //MeV }; } // namespace recon diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index 554c162f64..5ab4f16f54 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -20,18 +20,18 @@ void PileupFinder::configure(framework::config::Parameters& ps) { // get pileup candidates from PFlow and make a cleaned-up hit collection void PileupFinder::produce(framework::Event& event) { if (!event.exists(rec_hit_coll_name_,rec_hit_pass_name_)) { //ecal rechits - ldmx_log(error) << "Unable to find (one) collection named " << rec_hit_coll_name << - "_" << rec_hit_pass_name; + ldmx_log(error) << "Unable to find (one) collection named " << rec_hit_coll_name_ << + "_" << rec_hit_pass_name_; return; } if (!event.exists(pf_cand_coll_name_, pf_cand_pass_name_ )) { - ldmx_log(error) << "Unable to find (one) collection named " << pf_cand_coll_name << - "_" << pf_cand_pass_name; + ldmx_log(error) << "Unable to find (one) collection named " << pf_cand_coll_name_ << + "_" << pf_cand_pass_name_; return; } if (!event.exists(cluster_coll_name_, cluster_pass_name_ )) { - ldmx_log(error) << "Unable to find (one) collection named " << cluster_coll_name << - "_" << cluster_pass_name; + ldmx_log(error) << "Unable to find (one) collection named " << cluster_coll_name_ << + "_" << cluster_pass_name_; return; } @@ -49,7 +49,7 @@ void PileupFinder::produce(framework::Event& event) { // if a rechit is not on that list, add to output collection. std::vector output_hits; - for (const auto& pf_cand : pfCandidates){ + for (const auto& pf_cand : pf_cands){ if(pf_cand.getPID()==3 || pf_cand.getPID()==7){ double mom_vec = pf_cand.getTrackPxPyPz(); @@ -70,20 +70,20 @@ void PileupFinder::produce(framework::Event& event) { // try out hit association (print list?) int pf_cl_idx = pf_cand.getEcalIndex(); auto cl = clusters[pf_cl_idx]; - auto hitIds = cl.getHitIDs(); + auto hitIDs = cl.getHitIDs(); // make a collection without pileup hits for (auto hit : ecal_hits ) { - - int *foundIndex = std::find(std::begin(hitIDs), std::end(hitIDs), hit.id()); + + int *foundIndex = std::find(std::begin(hitIDs), std::end(hitIDs), hit.getID()); // When the element is not found, std::find returns the end of the range if (foundIndex == std::end(hitIDs)) { //hit not found in the pileup cluster output_hits.emplace_back(hit); //keep it - } - }// if trk/ecal matched + } + }// over hits + }// if trk/ecal matched }// over PF objects - - } - + +} void PileupFinder::onProcessEnd() { ldmx_log(debug) << "Process ends!"; delete eCorr_; From 94b16e55915038d30fd1da0da7c6646a6faf32c0 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 15:29:05 +0200 Subject: [PATCH 24/47] Fix more compilation problems --- Recon/src/Recon/PileupFinder.cxx | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index 5ab4f16f54..c4468e9435 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -52,23 +52,23 @@ void PileupFinder::produce(framework::Event& event) { for (const auto& pf_cand : pf_cands){ if(pf_cand.getPID()==3 || pf_cand.getPID()==7){ - double mom_vec = pf_cand.getTrackPxPyPz(); + std::vector mom_vec = pf_cand.getTrackPxPyPz(); float mom = mom_vec[0]*mom_vec[0]+mom_vec[1]*mom_vec[1]+mom_vec[2]*mom_vec[2]; mom = sqrt(mom); if (mom < min_mom_) continue; - ldmx_log(trace) << "Got pileup candidate with PID = " << pf_cand.getPID() << " and momentum = " << mom << " MeV." - - // now! use the hit-candidate association to get the associated ecal hits. - - - // TODO: - // write a header for this file --> DONE - // clean up all the methods below here --> DONE - // compile, possibly incrememnt particleflow candidate class imp nb - // try out hit association (print list?) - int pf_cl_idx = pf_cand.getEcalIndex(); + ldmx_log(trace) << "Got pileup candidate with PID = " << pf_cand.getPID() << " and momentum = " << mom << " MeV." ; + + // now! use the hit-candidate association to get the associated ecal hits. + + + // TODO: + // write a header for this file --> DONE + // clean up all the methods below here --> DONE + // compile, possibly incrememnt particleflow candidate class imp nb + // try out hit association (print list?) + int pf_cl_idx = pf_cand.getEcalIndex(); auto cl = clusters[pf_cl_idx]; auto hitIDs = cl.getHitIDs(); // make a collection without pileup hits @@ -86,8 +86,6 @@ void PileupFinder::produce(framework::Event& event) { } void PileupFinder::onProcessEnd() { ldmx_log(debug) << "Process ends!"; - delete eCorr_; - delete hCorr_; return; } From 3a075b57ac51031d7c2c4b8dcaf0cc1c9ae968d8 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 15:30:53 +0200 Subject: [PATCH 25/47] Fix more compilation problems --- Recon/src/Recon/PileupFinder.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index c4468e9435..0e93ee7eb5 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -74,7 +74,7 @@ void PileupFinder::produce(framework::Event& event) { // make a collection without pileup hits for (auto hit : ecal_hits ) { - int *foundIndex = std::find(std::begin(hitIDs), std::end(hitIDs), hit.getID()); + auto foundIndex = std::find(std::begin(hitIDs), std::end(hitIDs), hit.getID()); // When the element is not found, std::find returns the end of the range if (foundIndex == std::end(hitIDs)) { //hit not found in the pileup cluster output_hits.emplace_back(hit); //keep it From 19a2b99d492c08029b3165ceb232977a78dc4aa2 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 16:14:39 +0200 Subject: [PATCH 26/47] Make idx getters const --- Recon/include/Recon/Event/PFCandidate.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Recon/include/Recon/Event/PFCandidate.h b/Recon/include/Recon/Event/PFCandidate.h index 12539ea608..61e533b668 100644 --- a/Recon/include/Recon/Event/PFCandidate.h +++ b/Recon/include/Recon/Event/PFCandidate.h @@ -141,9 +141,9 @@ class PFCandidate { return {pos_hcal_x_, pos_hcal_y_, pos_hcal_z_}; } // associate component indices to the pf candidate - int getTrackIndex() { return track_idx_; } - int getEcalIndex() { return ecal_idx_ ; } - int getHcalIndex() { return hcal_idx_ ; } + int getTrackIndex() const { return track_idx_; } + int getEcalIndex() const { return ecal_idx_ ; } + int getHcalIndex() const { return hcal_idx_ ; } std::vector getTrackPxPyPz() const { return {track_px_, track_py_, track_pz_}; From 00ec9f858db7aac9c8924e41cecd6b64f9164df9 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 16:25:47 +0200 Subject: [PATCH 27/47] Add python config --- Recon/python/pileupFinder.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 Recon/python/pileupFinder.py diff --git a/Recon/python/pileupFinder.py b/Recon/python/pileupFinder.py new file mode 100644 index 0000000000..db4e8e9318 --- /dev/null +++ b/Recon/python/pileupFinder.py @@ -0,0 +1,23 @@ +"""Configuration for pfReco producers + +Sets all parameters to reasonable defaults. + +Examples +-------- + from LDMX.Recon.pfReco import pileupFinder + p.sequence.append( pileupFinder ) +""" + +from LDMX.Framework import ldmxcfg + +class pileupFinder(ldmxcfg.Producer) : + """Configuration for pileup finding from particle flow objects""" + def __init__(self, name='PileupFinder') : + super().__init__(name, 'recon::PileupFinder','Recon') + self.rec_hit_coll_name = 'EcalRecHits' + self.rec_hit_pass_name = '' + self.cluster_coll_name = 'PFEcalClusters' + self.pf_cand_coll_name = 'PFCandidates' + self.pf_cand_pass_name = '' + self.output_rec_hit_coll_name = 'EcalRecHitsNoPileup' + self.min_mom = 4000. From 4eeb758bb20d03a0713c8b95666421bd1188a552 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 19:43:25 +0200 Subject: [PATCH 28/47] Add missing parameter in python config --- Recon/python/pileupFinder.py | 1 + 1 file changed, 1 insertion(+) diff --git a/Recon/python/pileupFinder.py b/Recon/python/pileupFinder.py index db4e8e9318..b125dc3315 100644 --- a/Recon/python/pileupFinder.py +++ b/Recon/python/pileupFinder.py @@ -17,6 +17,7 @@ def __init__(self, name='PileupFinder') : self.rec_hit_coll_name = 'EcalRecHits' self.rec_hit_pass_name = '' self.cluster_coll_name = 'PFEcalClusters' + self.cluster_pass_name = '' self.pf_cand_coll_name = 'PFCandidates' self.pf_cand_pass_name = '' self.output_rec_hit_coll_name = 'EcalRecHitsNoPileup' From 1a26ed1a1387eec919941befa690d91ffa6414cc Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 19:49:58 +0200 Subject: [PATCH 29/47] Another missing parameter in python config --- Recon/python/pileupFinder.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Recon/python/pileupFinder.py b/Recon/python/pileupFinder.py index b125dc3315..e9d28dc873 100644 --- a/Recon/python/pileupFinder.py +++ b/Recon/python/pileupFinder.py @@ -21,4 +21,4 @@ def __init__(self, name='PileupFinder') : self.pf_cand_coll_name = 'PFCandidates' self.pf_cand_pass_name = '' self.output_rec_hit_coll_name = 'EcalRecHitsNoPileup' - self.min_mom = 4000. + self.min_momentum = 4000. From d71dce6274eff356c8963d53d8880ebe1602b599 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 19:57:13 +0200 Subject: [PATCH 30/47] Fix python double vs float nuisance --- Recon/include/Recon/PileupFinder.h | 2 +- Recon/src/Recon/PileupFinder.cxx | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Recon/include/Recon/PileupFinder.h b/Recon/include/Recon/PileupFinder.h index 120b1dada4..745df7fdc0 100644 --- a/Recon/include/Recon/PileupFinder.h +++ b/Recon/include/Recon/PileupFinder.h @@ -49,7 +49,7 @@ class PileupFinder : public framework::Producer { std::string output_rec_hit_coll_name_; // configuration - float min_mom_{0.}; //MeV + double min_mom_{0.}; //MeV }; } // namespace recon diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index 0e93ee7eb5..7984f8f53b 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -14,7 +14,7 @@ void PileupFinder::configure(framework::config::Parameters& ps) { cluster_pass_name_ = ps.getParameter("cluster_pass_name"); output_rec_hit_coll_name_ = ps.getParameter("output_rec_hit_coll_name"); // Algorithm configuration - min_mom_ = ps.getParameter("min_momentum"); + min_mom_ = ps.getParameter("min_momentum"); } // get pileup candidates from PFlow and make a cleaned-up hit collection From e190e8c33079898bb17c7596215d7c2593d4cb6b Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 7 Jul 2025 20:57:40 +0200 Subject: [PATCH 31/47] Add output collection to the bus --- Recon/src/Recon/PileupFinder.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index 7984f8f53b..b086d0a3c2 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -51,7 +51,7 @@ void PileupFinder::produce(framework::Event& event) { for (const auto& pf_cand : pf_cands){ if(pf_cand.getPID()==3 || pf_cand.getPID()==7){ - + //we have both ecal cluster and track std::vector mom_vec = pf_cand.getTrackPxPyPz(); float mom = mom_vec[0]*mom_vec[0]+mom_vec[1]*mom_vec[1]+mom_vec[2]*mom_vec[2]; mom = sqrt(mom); @@ -83,6 +83,8 @@ void PileupFinder::produce(framework::Event& event) { }// if trk/ecal matched }// over PF objects + event.add(output_rec_hit_coll_name_, output_hits); + } void PileupFinder::onProcessEnd() { ldmx_log(debug) << "Process ends!"; From ba37289e1b6914f01b6434e967e857745843b67e Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Wed, 9 Jul 2025 11:40:23 +0200 Subject: [PATCH 32/47] Check that index was set before doing lookup --- Recon/src/Recon/PileupFinder.cxx | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index b086d0a3c2..d93802b38a 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -66,9 +66,12 @@ void PileupFinder::produce(framework::Event& event) { // TODO: // write a header for this file --> DONE // clean up all the methods below here --> DONE - // compile, possibly incrememnt particleflow candidate class imp nb + // compile, possibly increment particleflow candidate class imp nb --> DONE // try out hit association (print list?) int pf_cl_idx = pf_cand.getEcalIndex(); + ldmx_log(trace) << "Got Ecal cluster with index " << pf_cl_idx << " while cluster array length is " << clusters.size(); + if (pf_cl_idx < 0 ) // was never set + continue; auto cl = clusters[pf_cl_idx]; auto hitIDs = cl.getHitIDs(); // make a collection without pileup hits @@ -78,6 +81,7 @@ void PileupFinder::produce(framework::Event& event) { // When the element is not found, std::find returns the end of the range if (foundIndex == std::end(hitIDs)) { //hit not found in the pileup cluster output_hits.emplace_back(hit); //keep it + ldmx_log(trace) << "Got no-pileup hit! "; } }// over hits }// if trk/ecal matched From 74668fcc80fcd317cc3a63e8912986aaf255fc41 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Thu, 18 Sep 2025 09:49:31 +0200 Subject: [PATCH 33/47] resolve merge conflict --- Recon/src/Recon/ParticleFlow.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx index 9ef8253e59..1dc3c509c4 100644 --- a/Recon/src/Recon/ParticleFlow.cxx +++ b/Recon/src/Recon/ParticleFlow.cxx @@ -275,6 +275,7 @@ void ParticleFlow::produce(framework::Event& event) { em_is_tk_linked[em_idx] = true; tk_is_em_linked[i] = true; tk_em_pairs[i] = em_idx; + cand.setEcalIndex(em_idx); break; } } From 1fb3d04f8173c3e2f3f13e8ac7d434ee3042b96c Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Thu, 18 Sep 2025 09:58:42 +0200 Subject: [PATCH 34/47] resolve merge conflict --- Recon/src/Recon/ParticleFlow.cxx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx index 1dc3c509c4..54c91fbd0d 100644 --- a/Recon/src/Recon/ParticleFlow.cxx +++ b/Recon/src/Recon/ParticleFlow.cxx @@ -275,7 +275,6 @@ void ParticleFlow::produce(framework::Event& event) { em_is_tk_linked[em_idx] = true; tk_is_em_linked[i] = true; tk_em_pairs[i] = em_idx; - cand.setEcalIndex(em_idx); break; } } @@ -319,6 +318,7 @@ void ParticleFlow::produce(framework::Event& event) { // chargedUnmatch.push_back(cand); } else { // if track is linked with ECal cluster fillCandEMCalo(cand, ecal_clusters[tk_em_pairs[i]]); + cand.setEcalIndex(tkEMPairs[i]); if (em_is_had_linked[tk_em_pairs[i]]) { // if ECal is linked with HCal // cluster fillCandHadCalo(cand, hcal_clusters[em_had_pairs[tk_em_pairs[i]]]); @@ -333,7 +333,6 @@ void ParticleFlow::produce(framework::Event& event) { for (int i = 0; i < ecal_clusters.size(); i++) { // already linked with ECal in the previous step if (em_is_tk_linked[i]) continue; - ldmx::PFCandidate cand; fillCandEMCalo(cand, ecal_clusters[i]); cand.setEcalIndex(i); From 656f7b23fd99fbb7ff205241ad0ee89ca319edf8 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Thu, 18 Sep 2025 10:06:08 +0200 Subject: [PATCH 35/47] resolve merge conflict --- .../include/TrigScint/TruthHitProducer.h | 9 +-- TrigScint/python/truthHits.py | 5 +- TrigScint/src/TrigScint/TruthHitProducer.cxx | 60 +++++++++---------- 3 files changed, 31 insertions(+), 43 deletions(-) diff --git a/TrigScint/include/TrigScint/TruthHitProducer.h b/TrigScint/include/TrigScint/TruthHitProducer.h index e900d9de32..4167bc8ca9 100644 --- a/TrigScint/include/TrigScint/TruthHitProducer.h +++ b/TrigScint/include/TrigScint/TruthHitProducer.h @@ -52,9 +52,7 @@ class TruthHitProducer : public framework::Producer { */ void produce(framework::Event &event) override; - /// Class to set the verbosity level. - // TODO: Make use of the global verbose parameter. - bool verbose_{false}; + private: /// Name of the input collection containing the sim hits std::string input_collection_; @@ -62,14 +60,13 @@ class TruthHitProducer : public framework::Producer { /// Name of the pass that the input collection is on (empty string means take /// any pass) std::string input_pass_name_; + /// Name of the pass that the input simparticles collection is on + std::string sim_particles_pass_name_; /// Name of the output collection that will be used to store the /// selected sim hits std::string output_collection_; - private: - std::string sim_particles_passname_; - std::string input_collection_events_passname_; }; // TruthHitProducer diff --git a/TrigScint/python/truthHits.py b/TrigScint/python/truthHits.py index e9d9023493..933f9dd355 100644 --- a/TrigScint/python/truthHits.py +++ b/TrigScint/python/truthHits.py @@ -20,10 +20,7 @@ def __init__(self,name) : self.input_collection="TriggerPadUpSimHits" self.input_pass_name="" #take any pass self.output_collection="truthBeamElectronsUp" - self.verbose = False - - self.sim_particles_passname = "" - self.input_collection_events_passname = "" + self.sim_particles_pass_name = "" def up() : """Get the beam electron truth hits for the trigger pad upstream of target""" diff --git a/TrigScint/src/TrigScint/TruthHitProducer.cxx b/TrigScint/src/TrigScint/TruthHitProducer.cxx index dbff995750..22f60b6c99 100644 --- a/TrigScint/src/TrigScint/TruthHitProducer.cxx +++ b/TrigScint/src/TrigScint/TruthHitProducer.cxx @@ -8,38 +8,33 @@ TruthHitProducer::TruthHitProducer(const std::string &name, : Producer(name, process) {} void TruthHitProducer::configure(framework::config::Parameters ¶meters) { - input_collection_ = parameters.get("input_collection"); - input_pass_name_ = parameters.get("input_pass_name"); - output_collection_ = parameters.get("output_collection"); - sim_particles_passname_ = - parameters.get("sim_particles_passname"); - input_collection_events_passname_ = - parameters.get("input_collection_events_passname"); - verbose_ = parameters.get("verbose"); + input_collection_ = parameters.getParameter("input_collection"); + input_pass_name_ = parameters.getParameter("input_pass_name"); + output_collection_ = parameters.getParameter("output_collection"); + sim_particles_pass_name_ = + parameters.getParameter("sim_particles_pass_name"); + + ldmx_log(info) << "In TruthHitProducer: configure done!"; + ldmx_log(info) << "Got parameters: " << "\nInput collection: " + << input_collection_ + << "\nInput pass name: " << input_pass_name_ + << "\nOutput collection: " << output_collection_; - if (verbose_) { - ldmx_log(info) << "In TruthHitProducer: configure done!"; - ldmx_log(info) << "Got parameters: " << "\nInput collection: " - << input_collection_ - << "\nInput pass name: " << input_pass_name_ - << "\nOutput collection: " << output_collection_ - << "\nVerbose: " << verbose_; - } } void TruthHitProducer::produce(framework::Event &event) { // Check if the collection exists. If not, don't bother processing the event. - if (!event.exists(input_collection_, input_collection_events_passname_)) { + if (!event.exists(input_collection_, input_pass_name_)) { ldmx_log(error) << "No input collection called " << input_collection_ << " found; skipping!"; return; } // looper over sim hits and aggregate energy depositions for each detID - const auto sim_hits{event.getCollection( + const auto simHits{event.getCollection( input_collection_, input_pass_name_)}; - auto particle_map{event.getMap( - "SimParticles", sim_particles_passname_)}; + auto particleMap{event.getMap( + "SimParticles", sim_particles_pass_name_)}; std::vector truth_beam_electrons; @@ -48,16 +43,15 @@ void TruthHitProducer::produce(framework::Event &event) { bool keep{false}; // check if hit is from beam electron and, if so, add to output collection for (int i = 0; i < sim_hit.getNumberOfContribs(); i++) { - auto contrib = sim_hit.getContrib(i); - if (verbose_) { - ldmx_log(debug) << "contrib " << i << " trackID: " << contrib.track_id_ - << " pdgID: " << contrib.pdg_code_ - << " edep: " << contrib.edep_; - ldmx_log(debug) << "\t particle id: " - << particle_map[contrib.track_id_].getPdgID() - << " particle status: " - << particle_map[contrib.track_id_].getGenStatus(); - } + auto contrib = simHit.getContrib(i); + ldmx_log(trace) << "contrib " << i << " trackID: " << contrib.trackID + << " pdgID: " << contrib.pdgCode + << " edep: " << contrib.edep; + ldmx_log(trace) << "\t particle id: " + << particleMap[contrib.trackID].getPdgID() + << " particle status: " + << particleMap[contrib.trackID].getGenStatus(); + // if the trackID is in the map if (particle_map.find(contrib.track_id_) != particle_map.end()) { // beam electron (PDGID = 11, genStatus == 1) @@ -66,9 +60,9 @@ void TruthHitProducer::produce(framework::Event &event) { keep = true; } } - if (keep) truth_beam_electrons.push_back(sim_hit); - } - } + if (keep) truthBeamElectrons.push_back(sim_hit); + }//over simhit contribs + }//over simhits event.add(output_collection_, truth_beam_electrons); } } // namespace trigscint From 6f0a26ea5eb22518e83a29f0ed006894d00ac13f Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Thu, 18 Sep 2025 11:05:38 +0200 Subject: [PATCH 36/47] resolve merge conflict --- DQM/include/DQM/EcalClusterAnalyzer.h | 9 ++- DQM/python/dqm.py | 25 +++++--- DQM/src/DQM/EcalClusterAnalyzer.cxx | 85 +++++++++++++++++++++++---- 3 files changed, 100 insertions(+), 19 deletions(-) diff --git a/DQM/include/DQM/EcalClusterAnalyzer.h b/DQM/include/DQM/EcalClusterAnalyzer.h index eb3790cade..5cc46c0ea3 100644 --- a/DQM/include/DQM/EcalClusterAnalyzer.h +++ b/DQM/include/DQM/EcalClusterAnalyzer.h @@ -64,7 +64,14 @@ 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 diff --git a/DQM/python/dqm.py b/DQM/python/dqm.py index 6cefca5b03..df7a92d6a9 100644 --- a/DQM/python/dqm.py +++ b/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,29 @@ 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.build1DHistogram("correctly_predicted_events", "Correctly predicted events", 3, 0., 3.) - 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(), diff --git a/DQM/src/DQM/EcalClusterAnalyzer.cxx b/DQM/src/DQM/EcalClusterAnalyzer.cxx index 766b4df00c..d5c717fa73 100644 --- a/DQM/src/DQM/EcalClusterAnalyzer.cxx +++ b/DQM/src/DQM/EcalClusterAnalyzer.cxx @@ -16,7 +16,12 @@ void EcalClusterAnalyzer::configure(framework::config::Parameters& ps) { cluster_coll_name_ = ps.get("cluster_coll_name"); cluster_pass_name_ = ps.get("cluster_pass_name"); - ecal_sp_hits_passname_ = ps.get("ecal_sp_hits_passname"); + ecal_sp_hits_coll_name_ = + ps.getParameter("ecal_sp_hits_coll_name"); + ecal_sp_hits_pass_name_ = + ps.getParameter("ecal_sp_hits_pass_name"); + mixed_hit_cutoff_ = ps.getParameter("mixed_hit_cutoff"); + return; } @@ -42,6 +47,7 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { layer_cluster_count[layer]++; } + int total_clusters = 0; for (const auto& [layer, count] : layer_cluster_count) { total_clusters += count; @@ -109,6 +115,8 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { << ", TS counted electrons: " << nbr_of_electrons << ", SP electrons: " << sp_electron_positions.size(); + std::map true_energy; + std::map 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 @@ -119,12 +127,20 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { pos2 = sp_electron_positions[1]; 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 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 +149,53 @@ 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 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 " << c.edep; + tot_event_energy+=c.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 nMixed++; histograms_.fill("ancestors", tag); hit_info.insert({hit.getID(), std::make_pair(tag, edep)}); } // end if simhit found @@ -166,6 +204,16 @@ 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",dist, (float)nMixed/ecal_rec_hits.size()); + ldmx_log(debug) << "Got " << nMixed << " mixed hits, a fraction of " << (float)nMixed/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 +260,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 +298,10 @@ 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 +324,23 @@ 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 - - histograms_.fill("clusterless_hits", (ecal_rec_hits.size() - clustered_hits)); + 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("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()); } From a57c4d266e67afb9c42c4997425160192c724be5 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Thu, 18 Sep 2025 11:28:37 +0200 Subject: [PATCH 37/47] resolve merge conflict --- Ecal/include/Ecal/CLUE.h | 3 + Ecal/src/Ecal/CLUE.cxx | 145 +++++++++++-------- Recon/src/Recon/PileupFinder.cxx | 32 ++-- TrigScint/src/TrigScint/TruthHitProducer.cxx | 6 + 4 files changed, 106 insertions(+), 80 deletions(-) diff --git a/Ecal/include/Ecal/CLUE.h b/Ecal/include/Ecal/CLUE.h index 3dc7b14ae8..8fa9290478 100644 --- a/Ecal/include/Ecal/CLUE.h +++ b/Ecal/include/Ecal/CLUE.h @@ -79,6 +79,9 @@ class CLUE { std::vector> setup( const std::vector& hits); + // get distance between clusters in the first layers, proxy for electron sep. + void electronSeparation(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 diff --git a/Ecal/src/Ecal/CLUE.cxx b/Ecal/src/Ecal/CLUE.cxx index 169a9001d6..03f4942fd4 100644 --- a/Ecal/src/Ecal/CLUE.cxx +++ b/Ecal/src/Ecal/CLUE.cxx @@ -34,64 +34,68 @@ T CLUE::dist(T x1, T y1, T z1, T x2, T y2, T z2) { 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 layer_index = 0; -// double layerZ = hits_[0].getZPos(); -// for (const auto& hit : hits_) { -// if (hit.getZPos() > layerZ + layerThickness[layer_index] + air) { -// layer_index++; -// // if (layer_index > limit) break; -// break; -// } -// firstLayers.push_back(hit); -// firstLayerClusters.push_back(IntermediateEcalCluster(hit, layer_index)); - -// } -// bool merge = false; -// do { -// merge = false; -// for (int i = 0; i < firstLayerClusters.size(); i++) { -// if (firstLayerClusters[i].empty()) continue; -// // if (firstLayerClusters[i].centroid().E() >= seed_threshold_) { -// 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; -// } -// } -// } + void CLUE::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.; + // sort hits in z + 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(IntermediateCluster(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; + } + } +} // Function to create layers from hits based on their layer number std::vector> CLUE::createLayers( @@ -99,6 +103,7 @@ std::vector> CLUE::createLayers( ldmx_log(trace) << "--- LAYER CREATION ---"; ldmx_log(trace) << "Number of layers: " << nbr_of_layers_; + // vector of layers, each layer having a vector of hits // initialize with nbr_of_layers_ empty vectors std::vector> hits_per_layer(nbr_of_layers_); @@ -162,7 +167,8 @@ std::vector> CLUE::setup( // collapse z dimension float x = roundToDecimal(hit->getXPos(), 4); float y = roundToDecimal(hit->getYPos(), 4); - ldmx_log(trace) << " New hit { x: " << x << " y: " << y << "}"; + float z = roundToDecimal(hit->getZPos(), 4); + ldmx_log(trace) << " New hit { x: " << x << " y: " << y << "}" << " (and z: " << z << ")"; std::pair coords; if (dc_ != 0 && nbr_of_layers_ > 1) { // if more than one layer, divide hit into densities with side dc @@ -188,11 +194,12 @@ std::vector> CLUE::setup( // density coords = {x, y}; } + if (density_map.find(coords) == density_map.end()) { density_map.emplace(coords, std::make_shared(x, y)); - ldmx_log(trace) << " New density created"; + ldmx_log(trace) << " * New density created"; } else { - ldmx_log(trace) << " Found density with x: " << density_map[coords]->x_ + ldmx_log(trace) << " --> Found density with x: " << density_map[coords]->x_ << " y: " << density_map[coords]->y_; } density_map[coords]->hits_.push_back(hit); @@ -334,6 +341,11 @@ std::vector> CLUE::clustering( density->total_energy_ > rhoc_ && density->delta_ > delta_c_mod; } else { is_seed = density->total_energy_ > rhoc_ && density->delta_ > deltac_; + if (is_seed) { + ldmx_log(trace) << " Distance to event centroid: " + << floatDist(density->x_, density->y_, + event_centroid_.centroid().x(), + event_centroid_.centroid().y()); } bool is_outlier = @@ -359,6 +371,8 @@ std::vector> CLUE::clustering( int& parent_index = density->follower_of_; if (parent_index != -1) followers[parent_index].push_back(density->index_); + else + ldmx_log(error) << " Somehow found a follower with parent index -1: id = " << density->index_; } else { ldmx_log(trace) << " This is an Outlier"; } @@ -447,6 +461,8 @@ std::vector> CLUE::setupForClue3D() { highest_energy = current_seed->total_energy_; ldmx_log(trace) << " Density with index " << current_seed->index_ << ", energy: " << current_seed->total_energy_; + << " position (x,y)= {" << current_seed->x_ << "," + << current_seed->y_ << ")"; int depth = 1; // decide delta_ and followerof from seeds in previous and next layer_ // do { @@ -454,7 +470,9 @@ std::vector> CLUE::setupForClue3D() { if ((layer - depth >= 0) && (layer - depth < seeds_.size())) { ldmx_log(trace) << " Looking at pre-layer: " << layer - depth; // look at previous layer + ldmx_log(trace) << " In previous layer... "; auto& previous_layer = seeds_[layer - depth]; + ldmx_log(trace) << " Got " << previous_layer.size() << " seeds"; for (const auto& prev_seed : previous_layer) { // for each seed in previous layer auto distance_2d_prev = dist(current_seed->x_, current_seed->y_, @@ -484,9 +502,11 @@ std::vector> CLUE::setupForClue3D() { } } } + if (layer + depth < nbr_of_layers_ && layer + depth < seeds_.size()) { ldmx_log(trace) << " Looking at post-layer: " << layer + depth; auto& next_layer = seeds_[layer + depth]; + ldmx_log(trace) << " Got " << next_layer.size() << " seeds"; for (const auto& next_seed : next_layer) { auto distance_2d_next = dist(current_seed->x_, current_seed->y_, next_seed->x_, next_seed->y_); @@ -517,6 +537,7 @@ std::vector> CLUE::setupForClue3D() { } // end of looking at next layer // } while (currentLayer[i]->layerFollowerOf == -1 && (layer - depth >= // 0 || layer + depth < nbr_of_layers_)); + ldmx_log(trace) << " Done setting parents."; densities.push_back(current_seed); } // TODO: This 2 needs to be configurable @@ -601,7 +622,7 @@ void CLUE::cluster(const std::vector& unsorted_hits, double dc, return a->getZPos() < b->getZPos(); }); - seeds_.resize(nbr_of_layers); + seeds_.resize(nbr_of_layers_); if (nbr_of_layers_ > 1) { ldmx_log(debug) << "Creating layers"; diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index d93802b38a..34219c6345 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -48,7 +48,10 @@ void PileupFinder::produce(framework::Event& event) { // get the list of hits associated with pileup candidates // if a rechit is not on that list, add to output collection. std::vector output_hits; + std::vector pileup_hitIDs; + //this needs to be a two-step procedure: loop over all clusters deemed to be pileup to find all their associated hits + // then loop over that list to make a collection of output hits that doesn't contain them for (const auto& pf_cand : pf_cands){ if(pf_cand.getPID()==3 || pf_cand.getPID()==7){ //we have both ecal cluster and track @@ -61,32 +64,25 @@ void PileupFinder::produce(framework::Event& event) { ldmx_log(trace) << "Got pileup candidate with PID = " << pf_cand.getPID() << " and momentum = " << mom << " MeV." ; // now! use the hit-candidate association to get the associated ecal hits. - - - // TODO: - // write a header for this file --> DONE - // clean up all the methods below here --> DONE - // compile, possibly increment particleflow candidate class imp nb --> DONE - // try out hit association (print list?) int pf_cl_idx = pf_cand.getEcalIndex(); ldmx_log(trace) << "Got Ecal cluster with index " << pf_cl_idx << " while cluster array length is " << clusters.size(); if (pf_cl_idx < 0 ) // was never set continue; auto cl = clusters[pf_cl_idx]; auto hitIDs = cl.getHitIDs(); - // make a collection without pileup hits - for (auto hit : ecal_hits ) { - - auto foundIndex = std::find(std::begin(hitIDs), std::end(hitIDs), hit.getID()); - // When the element is not found, std::find returns the end of the range - if (foundIndex == std::end(hitIDs)) { //hit not found in the pileup cluster - output_hits.emplace_back(hit); //keep it - ldmx_log(trace) << "Got no-pileup hit! "; - } - }// over hits + // add to collection of pileup hits + pileup_hitIDs.insert(pileup_hitIDs.end(), hitIDs.begin(), hitIDs.end()); }// if trk/ecal matched }// over PF objects - + + for (auto hit : ecal_hits ) { + auto foundIndex = std::find(std::begin(pileup_hitIDs), std::end(pileup_hitIDs), hit.getID()); + // When the element is not found, std::find returns the end of the range + if (foundIndex == std::end(pileup_hitIDs)) { //hit not part of any pileup cluster + output_hits.emplace_back(hit); //keep it + ldmx_log(trace) << "Got no-pileup hit! "; + } + } event.add(output_rec_hit_coll_name_, output_hits); } diff --git a/TrigScint/src/TrigScint/TruthHitProducer.cxx b/TrigScint/src/TrigScint/TruthHitProducer.cxx index 22f60b6c99..56a64170e4 100644 --- a/TrigScint/src/TrigScint/TruthHitProducer.cxx +++ b/TrigScint/src/TrigScint/TruthHitProducer.cxx @@ -27,9 +27,15 @@ void TruthHitProducer::produce(framework::Event &event) { // Check if the collection exists. If not, don't bother processing the event. if (!event.exists(input_collection_, input_pass_name_)) { ldmx_log(error) << "No input collection called " << input_collection_ + << "_" << input_pass_name_ << " found; skipping!"; + return; + } + if (!event.exists("SimParticles",sim_particles_pass_name_)) { + ldmx_log(error) << "No input SimParticle collection with pass " << sim_particles_pass_name_ << " found; skipping!"; return; } + // looper over sim hits and aggregate energy depositions for each detID const auto simHits{event.getCollection( input_collection_, input_pass_name_)}; From da89764aba057bb1cac3554617b468d4c33b2b4d Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Thu, 18 Sep 2025 13:21:44 +0200 Subject: [PATCH 38/47] harmonize variables after merge, now compiles --- DQM/src/DQM/EcalClusterAnalyzer.cxx | 12 ++++++------ Ecal/src/Ecal/CLUE.cxx | 14 +++++++------- Recon/include/Recon/Event/CaloCluster.h | 2 +- Recon/include/Recon/Event/PFCandidate.h | 3 +-- Recon/src/Recon/ParticleFlow.cxx | 2 +- TrigScint/src/TrigScint/TruthHitProducer.cxx | 18 +++++++++--------- 6 files changed, 25 insertions(+), 26 deletions(-) diff --git a/DQM/src/DQM/EcalClusterAnalyzer.cxx b/DQM/src/DQM/EcalClusterAnalyzer.cxx index d5c717fa73..22b7c9164a 100644 --- a/DQM/src/DQM/EcalClusterAnalyzer.cxx +++ b/DQM/src/DQM/EcalClusterAnalyzer.cxx @@ -80,7 +80,7 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { // Determine the truth information for the recoil electron std::vector> sp_electron_positions; const auto& ecal_sp_hits{event.getCollection( - "EcalScoringPlaneHits", ecal_sp_hits_passname_)}; + "EcalScoringPlaneHits", ecal_sp_hits_pass_name_)}; std::vector sorted_sp_hits = ecal_sp_hits; std::sort(sorted_sp_hits.begin(), sorted_sp_hits.end(), @@ -163,8 +163,8 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { const auto& contrib = it->getContrib(i); // get origin electron ID ancestor = contrib.origin_id_; - ldmx_log(trace) << "\t\t\tAncestor ID " << ancestor << " with edep " << c.edep; - tot_event_energy+=c.edep; + 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_; @@ -195,7 +195,7 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { // if not tagged, hit was from a single electron (within acceptable purity) tag = ancestor; //prev_ancestor; } - else nMixed++; + else n_mixed++; histograms_.fill("ancestors", tag); hit_info.insert({hit.getID(), std::make_pair(tag, edep)}); } // end if simhit found @@ -204,8 +204,8 @@ 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",dist, (float)nMixed/ecal_rec_hits.size()); - ldmx_log(debug) << "Got " << nMixed << " mixed hits, a fraction of " << (float)nMixed/ecal_rec_hits.size(); + 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) + diff --git a/Ecal/src/Ecal/CLUE.cxx b/Ecal/src/Ecal/CLUE.cxx index 03f4942fd4..9bae87e775 100644 --- a/Ecal/src/Ecal/CLUE.cxx +++ b/Ecal/src/Ecal/CLUE.cxx @@ -341,13 +341,13 @@ std::vector> CLUE::clustering( density->total_energy_ > rhoc_ && density->delta_ > delta_c_mod; } else { is_seed = density->total_energy_ > rhoc_ && density->delta_ > deltac_; - if (is_seed) { - ldmx_log(trace) << " Distance to event centroid: " - << floatDist(density->x_, density->y_, - event_centroid_.centroid().x(), - event_centroid_.centroid().y()); + if (is_seed) { + ldmx_log(trace) << " Distance to event centroid: " + << dist(density->x_, density->y_, + event_centroid_.centroid().x(), + event_centroid_.centroid().y()); + } } - bool is_outlier = (density->total_energy_ < rhoc_) && (density->delta_ > deltao_); density->cluster_id_ = -1; @@ -460,7 +460,7 @@ std::vector> CLUE::setupForClue3D() { if (current_seed->total_energy_ > highest_energy) highest_energy = current_seed->total_energy_; ldmx_log(trace) << " Density with index " << current_seed->index_ - << ", energy: " << current_seed->total_energy_; + << ", energy: " << current_seed->total_energy_ << " position (x,y)= {" << current_seed->x_ << "," << current_seed->y_ << ")"; int depth = 1; diff --git a/Recon/include/Recon/Event/CaloCluster.h b/Recon/include/Recon/Event/CaloCluster.h index fa7aeca8a1..cbb85495a5 100644 --- a/Recon/include/Recon/Event/CaloCluster.h +++ b/Recon/include/Recon/Event/CaloCluster.h @@ -136,7 +136,7 @@ class CaloCluster { double getEDYDZ() const { return err_dydz_; } // get hit rawIDs - const std::vector& getHitIDs() const { return hitIDs_; } + const std::vector& getHitIDs() const { return hit_ids_; } // ability to store limited hit info const std::vector& getHitX() const { return hit_x_; } diff --git a/Recon/include/Recon/Event/PFCandidate.h b/Recon/include/Recon/Event/PFCandidate.h index 61e533b668..2a40cbff4e 100644 --- a/Recon/include/Recon/Event/PFCandidate.h +++ b/Recon/include/Recon/Event/PFCandidate.h @@ -52,12 +52,11 @@ class PFCandidate { pos_hcal_y_ = y_; pos_hcal_z_ = z_; } - void setTrackPxPyPz(float x_, float y_, float z_) { track_px_ = x_; track_py_ = y_; track_pz_ = z_; - + } // associate component indices to the pf candidate void setTrackIndex(int x) { track_idx_ = x; } void setEcalIndex(int x) { ecal_idx_ = x; } diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx index 54c91fbd0d..fb37f0c836 100644 --- a/Recon/src/Recon/ParticleFlow.cxx +++ b/Recon/src/Recon/ParticleFlow.cxx @@ -318,7 +318,7 @@ void ParticleFlow::produce(framework::Event& event) { // chargedUnmatch.push_back(cand); } else { // if track is linked with ECal cluster fillCandEMCalo(cand, ecal_clusters[tk_em_pairs[i]]); - cand.setEcalIndex(tkEMPairs[i]); + cand.setEcalIndex(tk_em_pairs[i]); if (em_is_had_linked[tk_em_pairs[i]]) { // if ECal is linked with HCal // cluster fillCandHadCalo(cand, hcal_clusters[em_had_pairs[tk_em_pairs[i]]]); diff --git a/TrigScint/src/TrigScint/TruthHitProducer.cxx b/TrigScint/src/TrigScint/TruthHitProducer.cxx index 56a64170e4..d73e7e823b 100644 --- a/TrigScint/src/TrigScint/TruthHitProducer.cxx +++ b/TrigScint/src/TrigScint/TruthHitProducer.cxx @@ -37,9 +37,9 @@ void TruthHitProducer::produce(framework::Event &event) { } // looper over sim hits and aggregate energy depositions for each detID - const auto simHits{event.getCollection( + const auto sim_hits{event.getCollection( input_collection_, input_pass_name_)}; - auto particleMap{event.getMap( + auto particle_map{event.getMap( "SimParticles", sim_particles_pass_name_)}; std::vector truth_beam_electrons; @@ -49,14 +49,14 @@ void TruthHitProducer::produce(framework::Event &event) { bool keep{false}; // check if hit is from beam electron and, if so, add to output collection for (int i = 0; i < sim_hit.getNumberOfContribs(); i++) { - auto contrib = simHit.getContrib(i); - ldmx_log(trace) << "contrib " << i << " trackID: " << contrib.trackID - << " pdgID: " << contrib.pdgCode - << " edep: " << contrib.edep; + auto contrib = sim_hit.getContrib(i); + ldmx_log(trace) << "contrib " << i << " track_id: " << contrib.track_id_ + << " pdgID: " << contrib.pdg_code_ + << " edep: " << contrib.edep_; ldmx_log(trace) << "\t particle id: " - << particleMap[contrib.trackID].getPdgID() + << particle_map[contrib.track_id_].getPdgID() << " particle status: " - << particleMap[contrib.trackID].getGenStatus(); + << particle_map[contrib.track_id_].getGenStatus(); // if the trackID is in the map if (particle_map.find(contrib.track_id_) != particle_map.end()) { @@ -66,7 +66,7 @@ void TruthHitProducer::produce(framework::Event &event) { keep = true; } } - if (keep) truthBeamElectrons.push_back(sim_hit); + if (keep) truth_beam_electrons.push_back(sim_hit); }//over simhit contribs }//over simhits event.add(output_collection_, truth_beam_electrons); From 864fdddcc02f564f45455cc20e091b6d94c88792 Mon Sep 17 00:00:00 2001 From: Lene Kristian Bryngemark Date: Mon, 29 Sep 2025 10:15:36 +0200 Subject: [PATCH 39/47] fix float-->int number of bins + double hist def --- DQM/python/dqm.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/DQM/python/dqm.py b/DQM/python/dqm.py index df7a92d6a9..5fcc87d861 100644 --- a/DQM/python/dqm.py +++ b/DQM/python/dqm.py @@ -1079,11 +1079,9 @@ def __init__(self,name='EcalClusterAnalyzer') : 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.build1DHistogram("correctly_predicted_events", "Correctly predicted events", 3, 0., 3.) - 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("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) From 608d622e242f1ae775ece14a089e3e9c99ba84a9 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Thu, 2 Oct 2025 11:37:25 +0000 Subject: [PATCH 40/47] Apply clang-format --- DQM/include/DQM/EcalClusterAnalyzer.h | 5 +- DQM/src/DQM/EcalClusterAnalyzer.cxx | 111 ++++++++++-------- Ecal/include/Ecal/CLUE.h | 2 +- Ecal/src/Ecal/CLUE.cxx | 77 ++++++------ Recon/include/Recon/Event/CaloCluster.h | 2 +- Recon/include/Recon/Event/PFCandidate.h | 23 ++-- Recon/include/Recon/PileupFinder.h | 3 +- Recon/src/Recon/ParticleFlow.cxx | 2 +- Recon/src/Recon/PileupFinder.cxx | 88 +++++++------- .../include/TrigScint/TruthHitProducer.h | 4 +- TrigScint/src/TrigScint/TruthHitProducer.cxx | 37 +++--- 11 files changed, 188 insertions(+), 166 deletions(-) diff --git a/DQM/include/DQM/EcalClusterAnalyzer.h b/DQM/include/DQM/EcalClusterAnalyzer.h index 5cc46c0ea3..05878b069f 100644 --- a/DQM/include/DQM/EcalClusterAnalyzer.h +++ b/DQM/include/DQM/EcalClusterAnalyzer.h @@ -66,12 +66,11 @@ class EcalClusterAnalyzer : public framework::Analyzer { // Collection Name for Scoring Plane hits std::string ecal_sp_hits_coll_name_; - // Pass Name for Scoring Plane hits + // 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 diff --git a/DQM/src/DQM/EcalClusterAnalyzer.cxx b/DQM/src/DQM/EcalClusterAnalyzer.cxx index 22b7c9164a..ac55e88de2 100644 --- a/DQM/src/DQM/EcalClusterAnalyzer.cxx +++ b/DQM/src/DQM/EcalClusterAnalyzer.cxx @@ -21,7 +21,7 @@ void EcalClusterAnalyzer::configure(framework::config::Parameters& ps) { ecal_sp_hits_pass_name_ = ps.getParameter("ecal_sp_hits_pass_name"); mixed_hit_cutoff_ = ps.getParameter("mixed_hit_cutoff"); - + return; } @@ -47,7 +47,6 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { layer_cluster_count[layer]++; } - int total_clusters = 0; for (const auto& [layer, count] : layer_cluster_count) { total_clusters += count; @@ -115,8 +114,8 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { << ", TS counted electrons: " << nbr_of_electrons << ", SP electrons: " << sp_electron_positions.size(); - std::map true_energy; - std::map delta_energy; + std::map true_energy; + std::map 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 @@ -127,20 +126,19 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { pos2 = sp_electron_positions[1]; 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; + double tot_event_energy = 0; std::vector 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) { @@ -149,53 +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(); + 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 edep; 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. "; + 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_; + 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_; - tot_origin_edep[ancestor] += contrib.edep_; - eTot+= contrib.edep_; - } + 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; + 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. + } // 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 (within acceptable purity) - tag = ancestor; //prev_ancestor; + 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; + } + } } - else n_mixed++; + if (!tagged) { + // 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 @@ -204,12 +214,18 @@ 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(); + 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) ); + 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; } @@ -299,9 +315,10 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { 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 ; - + 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)); @@ -324,18 +341,20 @@ 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())); + // 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[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 ; - + 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("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()); diff --git a/Ecal/include/Ecal/CLUE.h b/Ecal/include/Ecal/CLUE.h index 8fa9290478..4c2ab6ac81 100644 --- a/Ecal/include/Ecal/CLUE.h +++ b/Ecal/include/Ecal/CLUE.h @@ -81,7 +81,7 @@ class CLUE { // get distance between clusters in the first layers, proxy for electron sep. void electronSeparation(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 diff --git a/Ecal/src/Ecal/CLUE.cxx b/Ecal/src/Ecal/CLUE.cxx index 9bae87e775..d873da4465 100644 --- a/Ecal/src/Ecal/CLUE.cxx +++ b/Ecal/src/Ecal/CLUE.cxx @@ -34,17 +34,17 @@ T CLUE::dist(T x1, T y1, T z1, T x2, T y2, T z2) { likely merged => recluster Did not quite work and I don't remember the idea anymore but leaving the code here for inspo */ - void CLUE::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 }; +void CLUE::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.; - // sort hits in z - std::sort(hits.begin(), hits.end(), [](const - ldmx::EcalHit& a, const ldmx::EcalHit& b) { - return a.getZPos() < b.getZPos(); - }); - + // sort hits in z + 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; @@ -57,7 +57,6 @@ T CLUE::dist(T x1, T y1, T z1, T x2, T y2, T z2) { } firstLayers.push_back(hit); firstLayerClusters.push_back(IntermediateCluster(hit, layerTag)); - } bool merge = false; do { @@ -68,9 +67,9 @@ T CLUE::dist(T x1, T y1, T z1, T x2, T y2, T z2) { 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].centroid().Py(), + firstLayerClusters[j].centroid().Px(), + firstLayerClusters[j].centroid().Py()) < 8.) { firstLayerClusters[i].add(firstLayerClusters[j]); firstLayerClusters[j].clear(); merge = true; @@ -82,16 +81,15 @@ T CLUE::dist(T x1, T y1, T z1, T x2, T y2, T z2) { 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++) { + 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()); + firstLayerClusters[i].centroid().Py(), + firstLayerClusters[j].centroid().Px(), + firstLayerClusters[j].centroid().Py()); ldmx_log(trace) << "Dist to cluster " << j << ": " << d; } } @@ -103,7 +101,6 @@ std::vector> CLUE::createLayers( ldmx_log(trace) << "--- LAYER CREATION ---"; ldmx_log(trace) << "Number of layers: " << nbr_of_layers_; - // vector of layers, each layer having a vector of hits // initialize with nbr_of_layers_ empty vectors std::vector> hits_per_layer(nbr_of_layers_); @@ -168,7 +165,8 @@ std::vector> CLUE::setup( float x = roundToDecimal(hit->getXPos(), 4); float y = roundToDecimal(hit->getYPos(), 4); float z = roundToDecimal(hit->getZPos(), 4); - ldmx_log(trace) << " New hit { x: " << x << " y: " << y << "}" << " (and z: " << z << ")"; + ldmx_log(trace) << " New hit { x: " << x << " y: " << y << "}" + << " (and z: " << z << ")"; std::pair coords; if (dc_ != 0 && nbr_of_layers_ > 1) { // if more than one layer, divide hit into densities with side dc @@ -199,7 +197,8 @@ std::vector> CLUE::setup( density_map.emplace(coords, std::make_shared(x, y)); ldmx_log(trace) << " * New density created"; } else { - ldmx_log(trace) << " --> Found density with x: " << density_map[coords]->x_ + ldmx_log(trace) << " --> Found density with x: " + << density_map[coords]->x_ << " y: " << density_map[coords]->y_; } density_map[coords]->hits_.push_back(hit); @@ -341,12 +340,12 @@ std::vector> CLUE::clustering( density->total_energy_ > rhoc_ && density->delta_ > delta_c_mod; } else { is_seed = density->total_energy_ > rhoc_ && density->delta_ > deltac_; - if (is_seed) { - ldmx_log(trace) << " Distance to event centroid: " - << dist(density->x_, density->y_, - event_centroid_.centroid().x(), - event_centroid_.centroid().y()); - } + if (is_seed) { + ldmx_log(trace) << " Distance to event centroid: " + << dist(density->x_, density->y_, + event_centroid_.centroid().x(), + event_centroid_.centroid().y()); + } } bool is_outlier = (density->total_energy_ < rhoc_) && (density->delta_ > deltao_); @@ -371,8 +370,10 @@ std::vector> CLUE::clustering( int& parent_index = density->follower_of_; if (parent_index != -1) followers[parent_index].push_back(density->index_); - else - ldmx_log(error) << " Somehow found a follower with parent index -1: id = " << density->index_; + else + ldmx_log(error) + << " Somehow found a follower with parent index -1: id = " + << density->index_; } else { ldmx_log(trace) << " This is an Outlier"; } @@ -461,8 +462,8 @@ std::vector> CLUE::setupForClue3D() { highest_energy = current_seed->total_energy_; ldmx_log(trace) << " Density with index " << current_seed->index_ << ", energy: " << current_seed->total_energy_ - << " position (x,y)= {" << current_seed->x_ << "," - << current_seed->y_ << ")"; + << " position (x,y)= {" << current_seed->x_ << "," + << current_seed->y_ << ")"; int depth = 1; // decide delta_ and followerof from seeds in previous and next layer_ // do { @@ -470,9 +471,9 @@ std::vector> CLUE::setupForClue3D() { if ((layer - depth >= 0) && (layer - depth < seeds_.size())) { ldmx_log(trace) << " Looking at pre-layer: " << layer - depth; // look at previous layer - ldmx_log(trace) << " In previous layer... "; + ldmx_log(trace) << " In previous layer... "; auto& previous_layer = seeds_[layer - depth]; - ldmx_log(trace) << " Got " << previous_layer.size() << " seeds"; + ldmx_log(trace) << " Got " << previous_layer.size() << " seeds"; for (const auto& prev_seed : previous_layer) { // for each seed in previous layer auto distance_2d_prev = dist(current_seed->x_, current_seed->y_, @@ -506,7 +507,7 @@ std::vector> CLUE::setupForClue3D() { if (layer + depth < nbr_of_layers_ && layer + depth < seeds_.size()) { ldmx_log(trace) << " Looking at post-layer: " << layer + depth; auto& next_layer = seeds_[layer + depth]; - ldmx_log(trace) << " Got " << next_layer.size() << " seeds"; + ldmx_log(trace) << " Got " << next_layer.size() << " seeds"; for (const auto& next_seed : next_layer) { auto distance_2d_next = dist(current_seed->x_, current_seed->y_, next_seed->x_, next_seed->y_); diff --git a/Recon/include/Recon/Event/CaloCluster.h b/Recon/include/Recon/Event/CaloCluster.h index cbb85495a5..648b7ae30d 100644 --- a/Recon/include/Recon/Event/CaloCluster.h +++ b/Recon/include/Recon/Event/CaloCluster.h @@ -135,7 +135,7 @@ class CaloCluster { /// @brief Delta unc on unc in y-z plane double getEDYDZ() const { return err_dydz_; } - // get hit rawIDs + // get hit rawIDs const std::vector& getHitIDs() const { return hit_ids_; } // ability to store limited hit info diff --git a/Recon/include/Recon/Event/PFCandidate.h b/Recon/include/Recon/Event/PFCandidate.h index 2a40cbff4e..e678fd6a61 100644 --- a/Recon/include/Recon/Event/PFCandidate.h +++ b/Recon/include/Recon/Event/PFCandidate.h @@ -10,8 +10,8 @@ // ROOT #include "TObject.h" //For ClassDef // ldmx-sw objects -//#include "Ecal/Event/EcalHit.h" -//#include "Hcal/Event/HcalHit.h" +// #include "Ecal/Event/EcalHit.h" +// #include "Hcal/Event/HcalHit.h" namespace ldmx { @@ -115,14 +115,16 @@ class PFCandidate { * @param hit The digi hit's entry number in the events digi * collection. */ - //void setEcalHits(const std::vector hits) { ecal_hits_ = hits; } + // void setEcalHits(const std::vector hits) { ecal_hits_ + // = hits; } /** * Take in the hcal hits that make up the candidate. * @param hit The digi hit's entry number in the events digi * collection. */ - // void setHcalHits(const std::vector hits) { hcal_hits_ = hits; } + // void setHcalHits(const std::vector hits) { + // hcal_hits_ = hits; } /* Getters */ @@ -141,8 +143,8 @@ class PFCandidate { } // associate component indices to the pf candidate int getTrackIndex() const { return track_idx_; } - int getEcalIndex() const { return ecal_idx_ ; } - int getHcalIndex() const { return hcal_idx_ ; } + int getEcalIndex() const { return ecal_idx_; } + int getHcalIndex() const { return hcal_idx_; } std::vector getTrackPxPyPz() const { return {track_px_, track_py_, track_pz_}; @@ -184,20 +186,19 @@ class PFCandidate { double getTruthEnergy() { return truth_energy_; } int getTruthPdgId() { return truth_pdg_id_; } - /** + /** * Take in the ecal hits that make up the candidate. * @param hit The digi hit's entry number in the events digi * collection. */ - //std::vector getEcalHits() { return ecal_hits_; } + // std::vector getEcalHits() { return ecal_hits_; } /** * Take in the hcal hits that make up the candidate. * @param hit The digi hit's entry number in the events digi * collection. */ - // std::vector getHcalHits() { return hcal_hits_; } - + // std::vector getHcalHits() { return hcal_hits_; } private: /* Particle ID enum */ @@ -266,7 +267,7 @@ class PFCandidate { int track_idx_{-1}; int ecal_idx_{-1}; int hcal_idx_{-1}; - + /* The ROOT class definition. */ ClassDef(PFCandidate, 2); }; diff --git a/Recon/include/Recon/PileupFinder.h b/Recon/include/Recon/PileupFinder.h index 745df7fdc0..1fae41151a 100644 --- a/Recon/include/Recon/PileupFinder.h +++ b/Recon/include/Recon/PileupFinder.h @@ -37,7 +37,6 @@ class PileupFinder : public framework::Producer { virtual void onProcessEnd(); private: - // name of collections for PF input object to be passed std::string rec_hit_coll_name_; std::string rec_hit_pass_name_; @@ -49,7 +48,7 @@ class PileupFinder : public framework::Producer { std::string output_rec_hit_coll_name_; // configuration - double min_mom_{0.}; //MeV + double min_mom_{0.}; // MeV }; } // namespace recon diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx index fb37f0c836..b5cb7c0495 100644 --- a/Recon/src/Recon/ParticleFlow.cxx +++ b/Recon/src/Recon/ParticleFlow.cxx @@ -318,7 +318,7 @@ void ParticleFlow::produce(framework::Event& event) { // chargedUnmatch.push_back(cand); } else { // if track is linked with ECal cluster fillCandEMCalo(cand, ecal_clusters[tk_em_pairs[i]]); - cand.setEcalIndex(tk_em_pairs[i]); + cand.setEcalIndex(tk_em_pairs[i]); if (em_is_had_linked[tk_em_pairs[i]]) { // if ECal is linked with HCal // cluster fillCandHadCalo(cand, hcal_clusters[em_had_pairs[tk_em_pairs[i]]]); diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index 34219c6345..87dc4d8498 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -12,37 +12,38 @@ void PileupFinder::configure(framework::config::Parameters& ps) { pf_cand_pass_name_ = ps.getParameter("pf_cand_pass_name"); cluster_coll_name_ = ps.getParameter("cluster_coll_name"); cluster_pass_name_ = ps.getParameter("cluster_pass_name"); - output_rec_hit_coll_name_ = ps.getParameter("output_rec_hit_coll_name"); + output_rec_hit_coll_name_ = + ps.getParameter("output_rec_hit_coll_name"); // Algorithm configuration - min_mom_ = ps.getParameter("min_momentum"); + min_mom_ = ps.getParameter("min_momentum"); } // get pileup candidates from PFlow and make a cleaned-up hit collection void PileupFinder::produce(framework::Event& event) { - if (!event.exists(rec_hit_coll_name_,rec_hit_pass_name_)) { //ecal rechits - ldmx_log(error) << "Unable to find (one) collection named " << rec_hit_coll_name_ << - "_" << rec_hit_pass_name_; - return; + if (!event.exists(rec_hit_coll_name_, rec_hit_pass_name_)) { // ecal rechits + ldmx_log(error) << "Unable to find (one) collection named " + << rec_hit_coll_name_ << "_" << rec_hit_pass_name_; + return; } - if (!event.exists(pf_cand_coll_name_, pf_cand_pass_name_ )) { - ldmx_log(error) << "Unable to find (one) collection named " << pf_cand_coll_name_ << - "_" << pf_cand_pass_name_; - return; + if (!event.exists(pf_cand_coll_name_, pf_cand_pass_name_)) { + ldmx_log(error) << "Unable to find (one) collection named " + << pf_cand_coll_name_ << "_" << pf_cand_pass_name_; + return; } - if (!event.exists(cluster_coll_name_, cluster_pass_name_ )) { - ldmx_log(error) << "Unable to find (one) collection named " << cluster_coll_name_ << - "_" << cluster_pass_name_; - return; + if (!event.exists(cluster_coll_name_, cluster_pass_name_)) { + ldmx_log(error) << "Unable to find (one) collection named " + << cluster_coll_name_ << "_" << cluster_pass_name_; + return; } const auto& ecal_hits{event.getCollection(rec_hit_coll_name_, rec_hit_pass_name_)}; - const auto& pf_cands{event.getCollection(pf_cand_coll_name_, - pf_cand_pass_name_)}; + const auto& pf_cands{event.getCollection( + pf_cand_coll_name_, pf_cand_pass_name_)}; - const auto& clusters{event.getCollection(cluster_coll_name_, - cluster_pass_name_)}; + const auto& clusters{event.getCollection( + cluster_coll_name_, cluster_pass_name_)}; // get PID 3 and 7 -- the ones with track and ecal matching // get the high-momentum track ones from there -- pileup candidate! // get the list of hits associated with pileup candidates @@ -50,41 +51,46 @@ void PileupFinder::produce(framework::Event& event) { std::vector output_hits; std::vector pileup_hitIDs; - //this needs to be a two-step procedure: loop over all clusters deemed to be pileup to find all their associated hits - // then loop over that list to make a collection of output hits that doesn't contain them - for (const auto& pf_cand : pf_cands){ - if(pf_cand.getPID()==3 || pf_cand.getPID()==7){ - //we have both ecal cluster and track + // this needs to be a two-step procedure: loop over all clusters deemed to be + // pileup to find all their associated hits + // then loop over that list to make a collection of output hits that doesn't + // contain them + for (const auto& pf_cand : pf_cands) { + if (pf_cand.getPID() == 3 || pf_cand.getPID() == 7) { + // we have both ecal cluster and track std::vector mom_vec = pf_cand.getTrackPxPyPz(); - float mom = mom_vec[0]*mom_vec[0]+mom_vec[1]*mom_vec[1]+mom_vec[2]*mom_vec[2]; + float mom = mom_vec[0] * mom_vec[0] + mom_vec[1] * mom_vec[1] + + mom_vec[2] * mom_vec[2]; mom = sqrt(mom); - if (mom < min_mom_) - continue; - ldmx_log(trace) << "Got pileup candidate with PID = " << pf_cand.getPID() << " and momentum = " << mom << " MeV." ; - - // now! use the hit-candidate association to get the associated ecal hits. + if (mom < min_mom_) continue; + ldmx_log(trace) << "Got pileup candidate with PID = " << pf_cand.getPID() + << " and momentum = " << mom << " MeV."; + + // now! use the hit-candidate association to get the associated ecal hits. int pf_cl_idx = pf_cand.getEcalIndex(); - ldmx_log(trace) << "Got Ecal cluster with index " << pf_cl_idx << " while cluster array length is " << clusters.size(); - if (pf_cl_idx < 0 ) // was never set - continue; + ldmx_log(trace) << "Got Ecal cluster with index " << pf_cl_idx + << " while cluster array length is " << clusters.size(); + if (pf_cl_idx < 0) // was never set + continue; auto cl = clusters[pf_cl_idx]; auto hitIDs = cl.getHitIDs(); // add to collection of pileup hits pileup_hitIDs.insert(pileup_hitIDs.end(), hitIDs.begin(), hitIDs.end()); - }// if trk/ecal matched - }// over PF objects - - for (auto hit : ecal_hits ) { - auto foundIndex = std::find(std::begin(pileup_hitIDs), std::end(pileup_hitIDs), hit.getID()); - // When the element is not found, std::find returns the end of the range - if (foundIndex == std::end(pileup_hitIDs)) { //hit not part of any pileup cluster - output_hits.emplace_back(hit); //keep it + } // if trk/ecal matched + } // over PF objects + + for (auto hit : ecal_hits) { + auto foundIndex = std::find(std::begin(pileup_hitIDs), + std::end(pileup_hitIDs), hit.getID()); + // When the element is not found, std::find returns the end of the range + if (foundIndex == + std::end(pileup_hitIDs)) { // hit not part of any pileup cluster + output_hits.emplace_back(hit); // keep it ldmx_log(trace) << "Got no-pileup hit! "; } } event.add(output_rec_hit_coll_name_, output_hits); - } void PileupFinder::onProcessEnd() { ldmx_log(debug) << "Process ends!"; diff --git a/TrigScint/include/TrigScint/TruthHitProducer.h b/TrigScint/include/TrigScint/TruthHitProducer.h index 4167bc8ca9..e7cb6b1090 100644 --- a/TrigScint/include/TrigScint/TruthHitProducer.h +++ b/TrigScint/include/TrigScint/TruthHitProducer.h @@ -52,8 +52,7 @@ class TruthHitProducer : public framework::Producer { */ void produce(framework::Event &event) override; - private: - + private: /// Name of the input collection containing the sim hits std::string input_collection_; @@ -67,7 +66,6 @@ class TruthHitProducer : public framework::Producer { /// selected sim hits std::string output_collection_; - }; // TruthHitProducer } // namespace trigscint diff --git a/TrigScint/src/TrigScint/TruthHitProducer.cxx b/TrigScint/src/TrigScint/TruthHitProducer.cxx index d73e7e823b..80ead80902 100644 --- a/TrigScint/src/TrigScint/TruthHitProducer.cxx +++ b/TrigScint/src/TrigScint/TruthHitProducer.cxx @@ -8,31 +8,30 @@ TruthHitProducer::TruthHitProducer(const std::string &name, : Producer(name, process) {} void TruthHitProducer::configure(framework::config::Parameters ¶meters) { - input_collection_ = parameters.getParameter("input_collection"); input_pass_name_ = parameters.getParameter("input_pass_name"); - output_collection_ = parameters.getParameter("output_collection"); + output_collection_ = + parameters.getParameter("output_collection"); sim_particles_pass_name_ = parameters.getParameter("sim_particles_pass_name"); ldmx_log(info) << "In TruthHitProducer: configure done!"; ldmx_log(info) << "Got parameters: " << "\nInput collection: " - << input_collection_ - << "\nInput pass name: " << input_pass_name_ - << "\nOutput collection: " << output_collection_; - + << input_collection_ + << "\nInput pass name: " << input_pass_name_ + << "\nOutput collection: " << output_collection_; } void TruthHitProducer::produce(framework::Event &event) { // Check if the collection exists. If not, don't bother processing the event. if (!event.exists(input_collection_, input_pass_name_)) { - ldmx_log(error) << "No input collection called " << input_collection_ - << "_" << input_pass_name_ << " found; skipping!"; + ldmx_log(error) << "No input collection called " << input_collection_ << "_" + << input_pass_name_ << " found; skipping!"; return; } - if (!event.exists("SimParticles",sim_particles_pass_name_)) { - ldmx_log(error) << "No input SimParticle collection with pass " << sim_particles_pass_name_ - << " found; skipping!"; + if (!event.exists("SimParticles", sim_particles_pass_name_)) { + ldmx_log(error) << "No input SimParticle collection with pass " + << sim_particles_pass_name_ << " found; skipping!"; return; } @@ -51,13 +50,13 @@ void TruthHitProducer::produce(framework::Event &event) { for (int i = 0; i < sim_hit.getNumberOfContribs(); i++) { auto contrib = sim_hit.getContrib(i); ldmx_log(trace) << "contrib " << i << " track_id: " << contrib.track_id_ - << " pdgID: " << contrib.pdg_code_ - << " edep: " << contrib.edep_; + << " pdgID: " << contrib.pdg_code_ + << " edep: " << contrib.edep_; ldmx_log(trace) << "\t particle id: " - << particle_map[contrib.track_id_].getPdgID() - << " particle status: " - << particle_map[contrib.track_id_].getGenStatus(); - + << particle_map[contrib.track_id_].getPdgID() + << " particle status: " + << particle_map[contrib.track_id_].getGenStatus(); + // if the trackID is in the map if (particle_map.find(contrib.track_id_) != particle_map.end()) { // beam electron (PDGID = 11, genStatus == 1) @@ -67,8 +66,8 @@ void TruthHitProducer::produce(framework::Event &event) { } } if (keep) truth_beam_electrons.push_back(sim_hit); - }//over simhit contribs - }//over simhits + } // over simhit contribs + } // over simhits event.add(output_collection_, truth_beam_electrons); } } // namespace trigscint From 9a371408f91a8d1bfa01c360c94c1075814416d2 Mon Sep 17 00:00:00 2001 From: bryngemark Date: Thu, 2 Oct 2025 16:32:48 +0200 Subject: [PATCH 41/47] Update Recon/include/Recon/Event/PFCandidate.h Remove remnant of attempt to set hits directly Co-authored-by: Tom Eichlersmith <31970302+tomeichlersmith@users.noreply.github.com> --- Recon/include/Recon/Event/PFCandidate.h | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/Recon/include/Recon/Event/PFCandidate.h b/Recon/include/Recon/Event/PFCandidate.h index e678fd6a61..a21846cea1 100644 --- a/Recon/include/Recon/Event/PFCandidate.h +++ b/Recon/include/Recon/Event/PFCandidate.h @@ -110,21 +110,6 @@ class PFCandidate { void setTruthEnergy(double x_) { truth_energy_ = x_; } void setTruthPdgId(int x_) { truth_pdg_id_ = x_; } - /** - * Take in the ecal hits that make up the candidate. - * @param hit The digi hit's entry number in the events digi - * collection. - */ - // void setEcalHits(const std::vector hits) { ecal_hits_ - // = hits; } - - /** - * Take in the hcal hits that make up the candidate. - * @param hit The digi hit's entry number in the events digi - * collection. - */ - // void setHcalHits(const std::vector hits) { - // hcal_hits_ = hits; } /* Getters */ From 32fb4db81dd6ae632846be4832682fdf7cad14e4 Mon Sep 17 00:00:00 2001 From: bryngemark Date: Thu, 2 Oct 2025 16:32:58 +0200 Subject: [PATCH 42/47] Update Recon/include/Recon/Event/PFCandidate.h Co-authored-by: Tom Eichlersmith <31970302+tomeichlersmith@users.noreply.github.com> --- Recon/include/Recon/Event/PFCandidate.h | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/Recon/include/Recon/Event/PFCandidate.h b/Recon/include/Recon/Event/PFCandidate.h index a21846cea1..e6bd1af43d 100644 --- a/Recon/include/Recon/Event/PFCandidate.h +++ b/Recon/include/Recon/Event/PFCandidate.h @@ -171,20 +171,6 @@ class PFCandidate { double getTruthEnergy() { return truth_energy_; } int getTruthPdgId() { return truth_pdg_id_; } - /** - * Take in the ecal hits that make up the candidate. - * @param hit The digi hit's entry number in the events digi - * collection. - */ - // std::vector getEcalHits() { return ecal_hits_; } - - /** - * Take in the hcal hits that make up the candidate. - * @param hit The digi hit's entry number in the events digi - * collection. - */ - // std::vector getHcalHits() { return hcal_hits_; } - private: /* Particle ID enum */ int pid_{0}; From d04e06dd5686f5106ec1eb0e683c4e892ac385c9 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Thu, 2 Oct 2025 14:45:43 +0000 Subject: [PATCH 43/47] Apply clang-tidy --- DQM/src/DQM/EcalClusterAnalyzer.cxx | 12 +++--- Ecal/src/Ecal/CLUE.cxx | 58 ++++++++++++++--------------- Recon/src/Recon/PileupFinder.cxx | 14 +++---- 3 files changed, 42 insertions(+), 42 deletions(-) diff --git a/DQM/src/DQM/EcalClusterAnalyzer.cxx b/DQM/src/DQM/EcalClusterAnalyzer.cxx index ac55e88de2..676248c986 100644 --- a/DQM/src/DQM/EcalClusterAnalyzer.cxx +++ b/DQM/src/DQM/EcalClusterAnalyzer.cxx @@ -155,7 +155,7 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { int tag = 0; std::vector edep; edep.resize(nbr_of_electrons_ + 1); - double eTot = 0; // keep track of total from all counted ancestors + double e_tot = 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++) { @@ -170,7 +170,7 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { if (ancestor <= nbr_of_electrons) { edep[ancestor] += contrib.edep_; tot_origin_edep[ancestor] += contrib.edep_; - eTot += contrib.edep_; + e_tot += contrib.edep_; } if (!tagged && i != 0 && prev_ancestor != ancestor) { // if origin electron ID does not match previous origin electron ID @@ -186,7 +186,7 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { // from a second electron. if (tagged) { for (int i = 1; i < nbr_of_electrons_ + 1; i++) { - if (edep[i] / eTot > + if (edep[i] / e_tot > 1 - mixed_hit_cutoff_) { // one ancestor contributes at least the // complement to the allowed mixing // fraction @@ -220,14 +220,14 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { << (float)n_mixed / ecal_rec_hits.size(); if (ecal_clusters.size() >= 2) { - float dR = std::sqrt(std::pow((ecal_clusters[0].getCentroidX() - + float d_r = 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; + histograms_.fill("cluster_distance", d_r); + ldmx_log(trace) << "Gt cluster distance (0,1) = " << d_r; } for (const auto& cl : ecal_clusters) { diff --git a/Ecal/src/Ecal/CLUE.cxx b/Ecal/src/Ecal/CLUE.cxx index d873da4465..6718fb7d0c 100644 --- a/Ecal/src/Ecal/CLUE.cxx +++ b/Ecal/src/Ecal/CLUE.cxx @@ -35,7 +35,7 @@ T CLUE::dist(T x1, T y1, T z1, T x2, T y2, T z2) { anymore but leaving the code here for inspo */ void CLUE::electronSeparation(std::vector hits) { - std::vector layerThickness = {2., 3.5, 5.3, 5.3, 5.3, 5.3, + 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}; double air = 10.; @@ -45,33 +45,33 @@ void CLUE::electronSeparation(std::vector hits) { return a.getZPos() < b.getZPos(); }); - std::vector firstLayers; - std::vector firstLayerClusters; - int layerTag = 0; - double layerZ = hits[0].getZPos(); + std::vector first_layers; + std::vector first_layer_clusters; + int layer_tag = 0; + double layer_z = hits[0].getZPos(); for (const auto& hit : hits) { - if (hit.getZPos() > layerZ + layerThickness[layerTag] + air) { - layerTag++; + if (hit.getZPos() > layer_z + layer_thickness[layer_tag] + air) { + layer_tag++; // if (layerTag > limit) break; break; } - firstLayers.push_back(hit); - firstLayerClusters.push_back(IntermediateCluster(hit, layerTag)); + first_layers.push_back(hit); + first_layer_clusters.push_back(IntermediateCluster(hit, layer_tag)); } bool merge = false; do { merge = false; - for (int i = 0; i < firstLayerClusters.size(); i++) { - if (firstLayerClusters[i].empty()) continue; + for (int i = 0; i < first_layer_clusters.size(); i++) { + if (first_layer_clusters[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(); + for (int j = i + 1; j < first_layer_clusters.size(); j++) { + if (first_layer_clusters[j].empty()) continue; + if (dist(first_layer_clusters[i].centroid().Px(), + first_layer_clusters[i].centroid().Py(), + first_layer_clusters[j].centroid().Px(), + first_layer_clusters[j].centroid().Py()) < 8.) { + first_layer_clusters[i].add(first_layer_clusters[j]); + first_layer_clusters[j].clear(); merge = true; } } @@ -79,17 +79,17 @@ void CLUE::electronSeparation(std::vector hits) { } } while (merge); ldmx_log(trace) << "--- ELECTRON SEPARATION ---"; - for (int i = 0; i < firstLayerClusters.size(); i++) { - if (firstLayerClusters[i].empty()) continue; + for (int i = 0; i < first_layer_clusters.size(); i++) { + if (first_layer_clusters[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()); + << " x: " << first_layer_clusters[i].centroid().Px() + << " y: " << first_layer_clusters[i].centroid().Py(); + for (int j = i + 1; j < first_layer_clusters.size(); j++) { + if (first_layer_clusters[j].empty()) continue; + auto d = dist(first_layer_clusters[i].centroid().Px(), + first_layer_clusters[i].centroid().Py(), + first_layer_clusters[j].centroid().Px(), + first_layer_clusters[j].centroid().Py()); ldmx_log(trace) << "Dist to cluster " << j << ": " << d; } } diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index 87dc4d8498..7d257d493e 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -49,7 +49,7 @@ void PileupFinder::produce(framework::Event& event) { // get the list of hits associated with pileup candidates // if a rechit is not on that list, add to output collection. std::vector output_hits; - std::vector pileup_hitIDs; + std::vector pileup_hit_i_ds; // this needs to be a two-step procedure: loop over all clusters deemed to be // pileup to find all their associated hits @@ -74,18 +74,18 @@ void PileupFinder::produce(framework::Event& event) { if (pf_cl_idx < 0) // was never set continue; auto cl = clusters[pf_cl_idx]; - auto hitIDs = cl.getHitIDs(); + auto hit_i_ds = cl.getHitIDs(); // add to collection of pileup hits - pileup_hitIDs.insert(pileup_hitIDs.end(), hitIDs.begin(), hitIDs.end()); + pileup_hit_i_ds.insert(pileup_hit_i_ds.end(), hit_i_ds.begin(), hit_i_ds.end()); } // if trk/ecal matched } // over PF objects for (auto hit : ecal_hits) { - auto foundIndex = std::find(std::begin(pileup_hitIDs), - std::end(pileup_hitIDs), hit.getID()); + auto found_index = std::find(std::begin(pileup_hit_i_ds), + std::end(pileup_hit_i_ds), hit.getID()); // When the element is not found, std::find returns the end of the range - if (foundIndex == - std::end(pileup_hitIDs)) { // hit not part of any pileup cluster + if (found_index == + std::end(pileup_hit_i_ds)) { // hit not part of any pileup cluster output_hits.emplace_back(hit); // keep it ldmx_log(trace) << "Got no-pileup hit! "; } From ca1b4cf05276bf0c2d2a0ec46925edf3a3be7a9b Mon Sep 17 00:00:00 2001 From: bryngemark Date: Thu, 13 Nov 2025 10:17:15 +0100 Subject: [PATCH 44/47] substitute x*x for pow(x,2) --- DQM/src/DQM/EcalClusterAnalyzer.cxx | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/DQM/src/DQM/EcalClusterAnalyzer.cxx b/DQM/src/DQM/EcalClusterAnalyzer.cxx index 676248c986..16fdcb9d51 100644 --- a/DQM/src/DQM/EcalClusterAnalyzer.cxx +++ b/DQM/src/DQM/EcalClusterAnalyzer.cxx @@ -220,12 +220,11 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { << (float)n_mixed / ecal_rec_hits.size(); if (ecal_clusters.size() >= 2) { - float d_r = std::sqrt(std::pow((ecal_clusters[0].getCentroidX() - - ecal_clusters[1].getCentroidX()), - 2) + - std::pow((ecal_clusters[0].getCentroidY() - - ecal_clusters[1].getCentroidY()), - 2)); + float d_x = ecal_clusters[0].getCentroidX() - + ecal_clusters[1].getCentroidX(); + float d_y = ecal_clusters[0].getCentroidY() - + ecal_clusters[1].getCentroidY(); + float d_r = std::sqrt( d_x*d_x + d_y*d_y ); histograms_.fill("cluster_distance", d_r); ldmx_log(trace) << "Gt cluster distance (0,1) = " << d_r; } From b2eedbea9140dc35c0eeb4475be392c77f24e13a Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Tue, 25 Nov 2025 08:28:51 +0000 Subject: [PATCH 45/47] apply clang-format and clang-tidy --- DQM/src/DQM/EcalClusterAnalyzer.cxx | 10 +++++----- Ecal/src/Ecal/CLUE.cxx | 4 ++-- Recon/src/Recon/PileupFinder.cxx | 7 ++++--- 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/DQM/src/DQM/EcalClusterAnalyzer.cxx b/DQM/src/DQM/EcalClusterAnalyzer.cxx index ae0a903e99..3a64e54dd2 100644 --- a/DQM/src/DQM/EcalClusterAnalyzer.cxx +++ b/DQM/src/DQM/EcalClusterAnalyzer.cxx @@ -226,11 +226,11 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { << (float)n_mixed / ecal_rec_hits.size(); if (ecal_clusters.size() >= 2) { - float d_x = ecal_clusters[0].getCentroidX() - - ecal_clusters[1].getCentroidX(); - float d_y = ecal_clusters[0].getCentroidY() - - ecal_clusters[1].getCentroidY(); - float d_r = std::sqrt( d_x*d_x + d_y*d_y ); + float d_x = + ecal_clusters[0].getCentroidX() - ecal_clusters[1].getCentroidX(); + float d_y = + ecal_clusters[0].getCentroidY() - ecal_clusters[1].getCentroidY(); + float d_r = std::sqrt(d_x * d_x + d_y * d_y); histograms_.fill("cluster_distance", d_r); ldmx_log(trace) << "Gt cluster distance (0,1) = " << d_r; } diff --git a/Ecal/src/Ecal/CLUE.cxx b/Ecal/src/Ecal/CLUE.cxx index 6718fb7d0c..60c4849a2c 100644 --- a/Ecal/src/Ecal/CLUE.cxx +++ b/Ecal/src/Ecal/CLUE.cxx @@ -36,8 +36,8 @@ T CLUE::dist(T x1, T y1, T z1, T x2, T y2, T z2) { void CLUE::electronSeparation(std::vector hits) { 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}; + 5.3, 5.3, 5.3, 5.3, 5.3, 10.5, + 10.5, 10.5, 10.5, 10.5}; double air = 10.; // sort hits in z std::sort(hits.begin(), hits.end(), diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx index 7d257d493e..66f480b651 100644 --- a/Recon/src/Recon/PileupFinder.cxx +++ b/Recon/src/Recon/PileupFinder.cxx @@ -76,16 +76,17 @@ void PileupFinder::produce(framework::Event& event) { auto cl = clusters[pf_cl_idx]; auto hit_i_ds = cl.getHitIDs(); // add to collection of pileup hits - pileup_hit_i_ds.insert(pileup_hit_i_ds.end(), hit_i_ds.begin(), hit_i_ds.end()); + pileup_hit_i_ds.insert(pileup_hit_i_ds.end(), hit_i_ds.begin(), + hit_i_ds.end()); } // if trk/ecal matched } // over PF objects for (auto hit : ecal_hits) { auto found_index = std::find(std::begin(pileup_hit_i_ds), - std::end(pileup_hit_i_ds), hit.getID()); + std::end(pileup_hit_i_ds), hit.getID()); // When the element is not found, std::find returns the end of the range if (found_index == - std::end(pileup_hit_i_ds)) { // hit not part of any pileup cluster + std::end(pileup_hit_i_ds)) { // hit not part of any pileup cluster output_hits.emplace_back(hit); // keep it ldmx_log(trace) << "Got no-pileup hit! "; } From 801c25f7e07b570640566f675505a0ff64cba16b Mon Sep 17 00:00:00 2001 From: bryngemark Date: Tue, 25 Nov 2025 09:55:05 +0100 Subject: [PATCH 46/47] Add PFlow + pileup finding to 2e CI --- .../validation_samples/it_pileup/config.py | 54 +++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/.github/validation_samples/it_pileup/config.py b/.github/validation_samples/it_pileup/config.py index 07c840e0e0..6107ccbde0 100644 --- a/.github/validation_samples/it_pileup/config.py +++ b/.github/validation_samples/it_pileup/config.py @@ -104,6 +104,51 @@ hcal_veto = hcal.HcalVetoProcessor() hcal_veto.input_hit_pass_name = thisPassName +# Load and configure particle flow sequence. +# Here we use PF "tracking" and CLUE Ecal clustering +from LDMX.Recon import pfReco +trackPF = pfReco.pfTrackProducer() +trackPF.inputTrackCollName=trackPF.inputTrackCollName+overlayStr #"EcalScoringPlaneHitsOverlay" # +trackPF.input_pass_name=inputPassName +trackPF.doElectronTracking=True +# reference info +truthPF = pfReco.pfTruthProducer() + +# CLUE +import LDMX.Ecal.ecalClusters as cl +cluster = cl.EcalClusterProducer() +cluster.seedThreshold = 350. +cluster.cutoff = 10. +# CLUE specific +cluster.CLUE = True #CLUE +cluster.dc = 0.3 +cluster.rhoc = 550. +cluster.deltac = 10. +cluster.deltao = 40. +cluster.debug = False +cluster.nbrOfLayers = 1 +cluster.nbr_of_layers = cluster.nbrOfLayers +cluster.reclustering = True +cluster.rec_hit_pass_name=thisPassName #run on process+pileup + +# particle flow: +pfComb=pfReco.pfProducer() +pfComb.inputEcalCollName = cluster.cluster_coll_name # use CLUE +pfComb.input_ecal_passname = thisPassName +# trigger recasting existing CLUE to caloclusters +pfComb.use_existing_ecal_clusters = True + +from LDMX.Recon import pileupFinder +puFinder = pileupFinder.pileupFinder() +puFinder.rec_hit_pass_name=thisPassName +#needs recast caloclusters, not (CLUE) ecalclusters +puFinder.cluster_coll_name=pfComb.inputEcalCollName+"Cast" +puFinder.pf_cand_coll_name=pfComb.outputCollName +puFinder.min_momentum=3000. + +# Load pileup finder +from LDMX.Recon.pfReco import pileupFinder + # Load the DQM modules from LDMX.DQM import dqm @@ -204,6 +249,15 @@ p.sequence.extend(dqm_with_overlay) +# Add PFlow + pileup finding sequence +p.sequence.extend([ + cluster, + trackPF, + truthPF, + pfComb, + puFinder +)] + p.inputFiles = ['ecal_pn.root'] p.outputFiles= ['events.root'] p.histogramFile = 'hist.root' From df608cb15853f23c45e4ab0a1295a515a0f5e77a Mon Sep 17 00:00:00 2001 From: bryngemark Date: Thu, 27 Nov 2025 16:08:04 +0100 Subject: [PATCH 47/47] Remove parameters set to default and double CLUE --- .github/validation_samples/it_pileup/config.py | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/.github/validation_samples/it_pileup/config.py b/.github/validation_samples/it_pileup/config.py index 6107ccbde0..d85645d242 100644 --- a/.github/validation_samples/it_pileup/config.py +++ b/.github/validation_samples/it_pileup/config.py @@ -117,17 +117,9 @@ # CLUE import LDMX.Ecal.ecalClusters as cl cluster = cl.EcalClusterProducer() -cluster.seedThreshold = 350. -cluster.cutoff = 10. -# CLUE specific -cluster.CLUE = True #CLUE -cluster.dc = 0.3 -cluster.rhoc = 550. -cluster.deltac = 10. -cluster.deltao = 40. -cluster.debug = False +cluster.seedThreshold = 350. cluster.dc = 0.3 cluster.nbrOfLayers = 1 -cluster.nbr_of_layers = cluster.nbrOfLayers +cluster.nbr_of_layers = cluster.nbrOfLayers #just here to highlight confusion cluster.reclustering = True cluster.rec_hit_pass_name=thisPassName #run on process+pileup @@ -236,7 +228,6 @@ ecalVeto, ecalMip, ecal_veto_pnet, - ecal_cluster.EcalClusterProducer(), hcal_digi_reco, hcal_veto, *ts_digis,