From fc6aadda22f63d038bb498a82470d967dffee46b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Wed, 11 Mar 2026 16:27:54 +0100 Subject: [PATCH 01/20] Experimental commit to use ClusterProducer with amplitude weighting and half-integer values without weighting. --- ClusterViewerAnalyzer.cxx | 54 +++++++++++++++++ ClusterViewerFig.py | 7 +++ NewEventcreator.py | 24 ++++++++ .../TrigScint/Event/TrigScintCluster.h | 6 +- .../TrigScint/TrigScintClusterProducer.h | 6 ++ .../TrigScint/TrigScintClusterProducer.cxx | 23 +++++-- extruthfig.py | 60 +++++++++++++++++++ 7 files changed, 172 insertions(+), 8 deletions(-) create mode 100644 ClusterViewerAnalyzer.cxx create mode 100644 ClusterViewerFig.py create mode 100644 NewEventcreator.py create mode 100644 extruthfig.py diff --git a/ClusterViewerAnalyzer.cxx b/ClusterViewerAnalyzer.cxx new file mode 100644 index 0000000000..466dc6c329 --- /dev/null +++ b/ClusterViewerAnalyzer.cxx @@ -0,0 +1,54 @@ +//ANALYZER TO VIEW CLUSTERS + +#include "Framework/EventProcessor.h" +#include "Framework/Event.h" +#include "TrigScint/Event/TrigScintCluster.h" +#include +#include + + +class ClusterViewerAnalyzer : public framework::Analyzer { + public: + + std::ofstream outfile; + + ClusterViewerAnalyzer(const std::string& name, framework::Process& p) + : framework::Analyzer(name, p) { + + outfile.open("clusters.txt"); + } + + void analyze(const framework::Event& event) override { + + const auto& clusters1 = + event.getCollection("TriggerPad1Clusters","truthisolater"); + + const auto& clusters2 = + event.getCollection("TriggerPad2Clusters","truthisolater"); + + const auto& clusters3 = + event.getCollection("TriggerPad3Clusters","truthisolater"); + + std::cout << "Pad1 clusters: " << clusters1.size() << std::endl; + std::cout << "Pad2 clusters: " << clusters2.size() << std::endl; + std::cout << "Pad3 clusters: " << clusters3.size() << std::endl; + + size_t Clusters = std::min({clusters1.size(), clusters2.size(), clusters3.size()}); + + for (size_t i=0; i < Clusters; i++) { + + float c1 = clusters1[i].getCentroid(); + float c2 = clusters2[i].getCentroid(); + float c3 = clusters3[i].getCentroid(); + std::cout << "printed " << c1 << std::endl; + std::cout << "printed " << c2 << std::endl; + std::cout << "printed " << c3 << std::endl; + outfile << c1 << " " + << c2 << " " + << c3 << "\n"; + } + } + void onProcessEnd() override {outfile.close();} + }; + +DECLARE_ANALYZER(ClusterViewerAnalyzer) diff --git a/ClusterViewerFig.py b/ClusterViewerFig.py new file mode 100644 index 0000000000..b003c6dc47 --- /dev/null +++ b/ClusterViewerFig.py @@ -0,0 +1,7 @@ +#CONFIG FILE TO VIEW CLUSTERS + +from LDMX.Framework import ldmxcfg +p = ldmxcfg.Process('clusterviewer') +p.input_files = ["TruthEvents10.root"] #change for new naming conventions +#p.ampl_weighting = False +p.sequence = [ldmxcfg.Analyzer.from_file('ClusterViewerAnalyzer.cxx')] diff --git a/NewEventcreator.py b/NewEventcreator.py new file mode 100644 index 0000000000..8e3640c546 --- /dev/null +++ b/NewEventcreator.py @@ -0,0 +1,24 @@ +#Creates simulation sample .root file with TS info + +from LDMX.Framework import ldmxcfg +p = ldmxcfg.Process('test') + +from LDMX.SimCore import simulator as sim +from LDMX.SimCore import generators as gen + +mySim = sim.simulator( "mySim" ) +mySim.setDetector( 'ldmx-det-v14-8gev', True ) +mySim.generators = [ gen.single_8gev_e_upstream_tagger() ] +mySim.beamSpotSmear = [20.,80.,0.] +mySim.description = 'Basic test Simulation' + +p.sequence = [ mySim ] +p.run = 1 +p.max_events = 10 #new number of electrons +p.output_files = [ 'NewEvents10.root' ] #new output file name + +import LDMX.Ecal.ecal_geometry +import LDMX.Ecal.ecal_hardcoded_conditions +import LDMX.Hcal.hcal_geometry + + diff --git a/TrigScint/include/TrigScint/Event/TrigScintCluster.h b/TrigScint/include/TrigScint/Event/TrigScintCluster.h index bfedb52210..a6839ae3eb 100644 --- a/TrigScint/include/TrigScint/Event/TrigScintCluster.h +++ b/TrigScint/include/TrigScint/Event/TrigScintCluster.h @@ -93,7 +93,7 @@ class TrigScintCluster { /** * @param centroid The channel ID centroid */ - void setCentroid(double centroid) { centroid_ = centroid; } + void setCentroid(float centroid) { centroid_ = centroid; } /** Set time of hit. */ void setTime(float t) { time_ = t; } @@ -132,7 +132,7 @@ class TrigScintCluster { const std::vector &getHitIDs() const { return hit_ids_; } /** Get the cluster centroid in units of channel nb */ - double getCentroid() const { return centroid_; } + float getCentroid() const { return centroid_; } bool operator<(const TrigScintCluster &rhs) const { return this->getEnergy() < rhs.getEnergy(); @@ -156,7 +156,7 @@ class TrigScintCluster { // hit centroid in units of channel nb: energy weighted average of the IDs of // the hits forming the cluster - double centroid_{-1}; + float centroid_{-1}; // hit centroid in x [mm] (not implemented) double centroid_x_{0}; diff --git a/TrigScint/include/TrigScint/TrigScintClusterProducer.h b/TrigScint/include/TrigScint/TrigScintClusterProducer.h index 90411f79e5..f58f18faee 100644 --- a/TrigScint/include/TrigScint/TrigScintClusterProducer.h +++ b/TrigScint/include/TrigScint/TrigScintClusterProducer.h @@ -45,6 +45,9 @@ class TrigScintClusterProducer : public framework::Producer { // collection of clusters produced std::vector clusters_; + //amplitude weighting for centroid values + bool ampl_weighting{true}; + // cluster seeding threshold double seed_{0.}; @@ -89,6 +92,9 @@ class TrigScintClusterProducer : public framework::Producer { // edep content, only; leave val_ for PE float val_e_{0.}; + + // sum of hit cluster weights + float sumw_{0.}; // book keep which channels have already been added to the cluster at hand std::vector v_added_indices_; diff --git a/TrigScint/src/TrigScint/TrigScintClusterProducer.cxx b/TrigScint/src/TrigScint/TrigScintClusterProducer.cxx index da617658d2..93dfc88cac 100644 --- a/TrigScint/src/TrigScint/TrigScintClusterProducer.cxx +++ b/TrigScint/src/TrigScint/TrigScintClusterProducer.cxx @@ -7,6 +7,7 @@ namespace trigscint { void TrigScintClusterProducer::configure(framework::config::Parameters &ps) { + ampl_weighting = ps.get("ampl_weighting"); seed_ = ps.get("seed_threshold"); min_thr_ = ps.get("clustering_threshold"); max_width_ = ps.get("max_cluster_width"); @@ -402,7 +403,9 @@ void TrigScintClusterProducer::produce(framework::Event &event) { } // if adding another hit, going forward, was allowed // done adding hits to cluster. calculate centroid - centroid_ /= val_; // final weighting step: divide by total + centroid_ /= sumw_ ; // final weighting step: divide by total amplitude sum + + centroid_ -= 1; // shift back to actual channel center ldmx::TrigScintCluster cluster; @@ -445,6 +448,7 @@ void TrigScintClusterProducer::produce(framework::Event &event) { val_e_ = 0; beam_e_ = 0; time_ = 0; + sumw_ = 0; // book keep which channels have already been added to a cluster v_added_indices_.resize(0); @@ -476,13 +480,22 @@ void TrigScintClusterProducer::produce(framework::Event &event) { void TrigScintClusterProducer::addHit(uint idx, ldmx::TrigScintHit hit) { float ampl = hit.getPE(); - val_ += ampl; + float w = 1; + if (ampl_weighting) { + w = ampl; + } + + float energy = hit.getEnergy(); val_e_ += energy; - centroid_ += (idx + 1) * ampl; // need non-zero weight of channel 0. shifting - // centroid back by 1 in the end - // this number gets divided by val at the end + val_ += ampl; + centroid_ += (idx + 1) * w; // need non-zero weight of channel 0. shifting + // centroid back by 1 in the end + // this number gets divided by val at the end + + sumw_ += w; + v_added_indices_.push_back(idx); beam_e_ += hit.getBeamEfrac() * energy; diff --git a/extruthfig.py b/extruthfig.py new file mode 100644 index 0000000000..883fafe7e6 --- /dev/null +++ b/extruthfig.py @@ -0,0 +1,60 @@ +#example config for truth hit producer + +#Uses TruthHitProducer to isolate true hits, +#(then turns them to digis, then clusters) + +from LDMX.Framework import ldmxcfg + +p = ldmxcfg.Process('truthisolater') +p.input_files = ["NewEvents10.root"] +p.output_files = ["TruthEvents10.root"] + +from LDMX.TrigScint.truth_hits import TruthHitProducer +from LDMX.TrigScint.trig_scint import TrigScintDigiProducer + +truth_hits = [TruthHitProducer('beamElectronsPad1'), + TruthHitProducer('beamElectronsPad2'), + TruthHitProducer('beamElectronsPad3') + ] + +pad_num = 1 + +for hits in truth_hits: + hits.input_collection=f"TriggerPad{pad_num}SimHits" + hits.output_collection=f"truthBeamElectronsPad{pad_num}" + pad_num+=1 + +truth_digis = [TrigScintDigiProducer.pad1(), + TrigScintDigiProducer.pad2(), + TrigScintDigiProducer.pad3() + ] + +for digi,hits in zip(truth_digis, truth_hits): + digi.input_collection = hits.output_collection + +from LDMX.TrigScint.trig_scint import TrigScintClusterProducer + +truth_clusters = [ + TrigScintClusterProducer.pad1(), + TrigScintClusterProducer.pad2(), + TrigScintClusterProducer.pad3(), + ] + +for cluster, digi in zip(truth_clusters, truth_digis): + cluster.input_collection = digi.output_collection + cluster.ampl_weighting = True + + +p.sequence = [*truth_hits, + *truth_digis, + *truth_clusters, + #trigScintTrack, + #count, +# TriggerProcessor('trigger', 800.), + # ldmxcfg.Producer.from_file('TruthHitProducer.cxx') + ] + + + + + From 4f156aafa017b17b172df1b7ba76f40deb3d9fa2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Mon, 20 Apr 2026 14:43:22 +0200 Subject: [PATCH 02/20] Add LUT tracking method --- .../include/TrigScint/TrigScintLUTTracker.h | 125 +++ .../src/TrigScint/TrigScintLUTTracker.cxx | 844 ++++++++++++++++++ 2 files changed, 969 insertions(+) create mode 100644 TrigScint/include/TrigScint/TrigScintLUTTracker.h create mode 100644 TrigScint/src/TrigScint/TrigScintLUTTracker.cxx diff --git a/TrigScint/include/TrigScint/TrigScintLUTTracker.h b/TrigScint/include/TrigScint/TrigScintLUTTracker.h new file mode 100644 index 0000000000..2102fd5f83 --- /dev/null +++ b/TrigScint/include/TrigScint/TrigScintLUTTracker.h @@ -0,0 +1,125 @@ + +#ifndef TRIGSCINT_TRIGSCINTTRACKPRODUCER_H +#define TRIGSCINT_TRIGSCINTTRACKPRODUCER_H + +// LDMX Framework +#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 "Recon/Event/EventConstants.h" +#include "TrigScint/Event/TrigScintCluster.h" +#include "TrigScint/Event/TrigScintTrack.h" + +#include + +namespace trigscint { + +/** + * @class TrigScintTrackProducer + * @brief making tracks from trigger scintillator clusters + */ +class TrigScintLUTTracker : public framework::Producer { + public: + TrigScintLUTTracker(const std::string &name, framework::Process &process) + : Producer(name, process) {} + + void configure(framework::config::Parameters &ps) override; + + void produce(framework::Event &event) override; + + void onProcessStart() override; + void onProcessEnd() override; + + private: + // collection of produced tracks + std::vector tracks_; + + // add a cluster to a track + ldmx::TrigScintTrack makeTrack(std::vector clusters); + + // match x, y tracks and set their x,y spatial coordinates + void matchXYTracks(std::vector &tracks); + // std::vector matchXYTracks( + // std::vector &tracks); + + // maximum difference (in channel number space) between track seed and cluster + // in the next pad tolerated to form a track + double max_delta_{0.}; + + // producer specific verbosity + int verbose_{0}; + + // collection used to seed the tracks + std::string seeding_collection_; + + // other cluster collections used in track making + std::vector input_collections_; + + // output collection (tracks) + std::string output_collection_; + + // specific pass name to use for track making + std::string pass_name_{""}; + + // allow forming tracks without match in the last collection + bool skip_last_{false}; + + // vertical bar start index + int vert_bar_start_idx_{52}; + + // number of horizontal bars (one layer) in geometry + int n_bars_y_{16}; + + // number of vertical bars (one row) in geometry + int n_bars_x_{8}; + + // track centroid in units of channel nb (will not be content weighted) + // float centroid_{0.}; + + // track horizontal centroid in units of channel nb (will not be content + // weighted) + // float centroidX_{-1}; + + // track vertical centroid in units of channel nb (will not be content + // weighted) + // float centroidY_{-1}; + + // track residual in units of channel nb (will not be content weighted) + // float residual_{0.}; + + // LUT part starts here + + struct LUTKey { + float p1, p2, p3; + + bool operator==(const LUTKey& other) const { + return p1 == other.p1 && p2 == other.p2 && p3 == other.p3; + } + }; + + struct LUTKeyHash { + size_t operator()(const LUTKey& k) const { + return std::hash()(k.p1) ^ + (std::hash()(k.p2) << 1) ^ + (std::hash()(k.p3) << 2); + } + }; + + std::unordered_set lut_; + + // LUT ends here + + float bar_width_y_{3.}; // mm + float bar_gap_y_{2.1}; // mm + float bar_width_x_{3.}; // mm + float bar_gap_x_{0.1}; // mm + + float x_conv_factor_; // geometry conversion factors + float x_start_; + float y_conv_factor_; + float y_start_; + }; + +} // namespace trigscint + +#endif // TRIGSCINT_TRIGSCINTLUTTRACKER_H diff --git a/TrigScint/src/TrigScint/TrigScintLUTTracker.cxx b/TrigScint/src/TrigScint/TrigScintLUTTracker.cxx new file mode 100644 index 0000000000..23ec09f333 --- /dev/null +++ b/TrigScint/src/TrigScint/TrigScintLUTTracker.cxx @@ -0,0 +1,844 @@ +#include "TrigScint/TrigScintLUTTracker.h" + +#include // std::next +#include +#include +#include + +namespace trigscint { + +void TrigScintLUTTracker::configure(framework::config::Parameters &ps) { + max_delta_ = ps.get( + "delta_max"); // max distance to consider adding in a cluster to track + seeding_collection_ = ps.get( + "seeding_collection"); // probably tagger pad, "TriggerPadTagClusters" + input_collections_ = ps.get>( + "further_input_collections"); // {"TriggerPadUpClusters" , + // "TriggerPadDnClusters" } + output_collection_ = ps.get("output_collection"); + pass_name_ = ps.get("input_pass_name"); + verbose_ = ps.get("verbosity"); + vert_bar_start_idx_ = ps.get("vertical_bar_start_index"); + n_bars_y_ = ps.get("number_horizontal_bars"); + bar_width_y_ = ps.get("horizontal_bar_width"); + bar_gap_y_ = ps.get("horizontal_bar_gap"); + n_bars_x_ = ps.get("number_vertical_bars"); + bar_width_x_ = ps.get("vertical_bar_width"); + bar_gap_x_ = ps.get("vertical_bar_gap"); + skip_last_ = ps.get("allow_skip_last_collection"); + std::string lut_file_ = ps.get("lut_file"); //LUT ! new! + + // TO DO: allow any number of input collections + + if (verbose_) { + ldmx_log(info) << "In TrigScintLUTTracker: configure done!" << std::endl; + ldmx_log(info) << "Got parameters: \nSeeding: " << seeding_collection_ + << "\nTolerance: " << max_delta_ + << "\nInput: " << input_collections_.at(0) << " and " + << input_collections_.at(1) + << "\nInput pass name: " << pass_name_ + << "\nAllow tracks with no hit in last collection: " + << skip_last_ + << "\nVertical bar start index: " << vert_bar_start_idx_ + << "\nNumber of horizontal bars: " << n_bars_y_ + << "\nHorizontal bar width: " << bar_width_y_ + << "\nHorizontal bar gap: " << bar_gap_y_ + << "\nNumber of vertical bars: " << n_bars_x_ + << "\nVertical bar width: " << bar_width_x_ + << "\nVertical bar gap: " << bar_gap_x_ + << "\nOutput: " << output_collection_ + << "\nVerbosity: " << verbose_ + << "\nTracks (LUT): " << lut_file_; //is this necessary too? + } + // each bar only goes half this distance up (overlap/zig-zag) + y_conv_factor_ = (bar_width_y_ + bar_gap_y_) / 2.; + // half height of pad + y_start_ = -(n_bars_y_ * (bar_width_y_ + bar_gap_y_) - bar_gap_y_) / 2.; + // each bar goes entire distance sideways (no overlap) + x_conv_factor_ = bar_width_x_ + bar_gap_x_; + // half width of pad + x_start_ = -(n_bars_x_ * (bar_width_x_ + bar_gap_x_) - bar_gap_x_) / 2.; + + //also LUT related here + std::ifstream file(lut_file_); //open file + float a, b, c; + + while (file >> a >> b >> c) { + float p1 = a; + float p2 = b; + float p3 = c; + lut_.insert({p1, p2, p3}); + } + + ldmx_log(info) << "Loaded LUT with size: " << lut_.size(); + + return; +} //end of configure + + + +void TrigScintLUTTracker::produce(framework::Event &event) { + // parameters. + // one pad cluster collection to use as seed + // a vector with the other two + // a maximum distance between seed centroid and other pad clusters + // allowed to belong to the same track + // an output collection name a verbosity controller + + if (verbose_) { + ldmx_log(debug) + << "TrigScintLUTTracker: produce() starts! Event number: " + << event.getEventHeader().getEventNumber(); + } + + if (!event.exists(seeding_collection_, pass_name_)) { + ldmx_log(info) << "No collection called " << seeding_collection_ + << "; skipping event"; + // << "; still, not skipping event"; + return; + } + + if (!event.exists(seeding_collection_, pass_name_)) { + ldmx_log(info) << "No collection called " << seeding_collection_ + << "; skipping event"; + return; + } + + const auto seeds{event.getCollection( + seeding_collection_, pass_name_)}; + uint num_seeds = seeds.size(); + + if (verbose_) { + ldmx_log(debug) << "Got track seeding cluster collection " + << seeding_collection_ << " with " << num_seeds + << " entries "; + } + + if (!event.exists(input_collections_.at(0), pass_name_)) { + ldmx_log(info) << "No collection called " << input_collections_.at(0) + << "; skipping event"; + // << "; still, not skipping event"; + + return; + } + + const auto clusters_pad1{event.getCollection( + input_collections_.at(0), pass_name_)}; + + if (!event.exists(input_collections_.at(1), pass_name_)) { + ldmx_log(info) << "No collection called " + << input_collections_.at(1) + // << "; still, not skipping event"; + << "; skipping event"; + std::vector empty{}; + event.add(output_collection_, empty); + return; + } + + const auto clusters_pad2{event.getCollection( + input_collections_.at(1), pass_name_)}; + + if (verbose_) { + ldmx_log(debug) << "Got the other two pad collections:" + << input_collections_.at(0) << " with " + << clusters_pad1.size() << " entries, and " + << input_collections_.at(1) << " with " + << clusters_pad2.size() << " entries."; + } + + std::vector cleaned_tracks; + std::vector cleaned_tracks_y; + std::vector cleaned_tracks_x; + + + + +//----------------------------------------------------------------------------- + + //this is what 130-250 have been replaced with: + + // loop over the clusters in the seeding pad collection, if there are clusters + // in all pads + // bool skipDn = false; + if (num_seeds && clusters_pad1.size()) { //ends on 407, if there are clusters in all pads + // could check this explicitly here: and then just get out of all checks on + // the dn pad immediately + // if (! clusters_pad2.size()) + // skipDn = true ; + for (const auto &seed : seeds) { //ends on line 247, over all seeds + // for each seed, search through the other two pads to match all clusters + // with centroids within tolerance to tracks + float centroid = seed.getCentroid(); + + std::vector track_candidates; + + if (verbose_ > 1) { + ldmx_log(debug) << "Got seed with centroid " << centroid; + } + + // reset for each seed + // bool madeTrack = false; + + + //LUT part here! + for (const auto &cluster1 : clusters_pad1) { + for (const auto &cluster2 : clusters_pad2) { + + float seed_bin = seed.getCentroid(); + float pad1_bin = cluster1.getCentroid(); + float pad2_bin = cluster2.getCentroid(); + + LUTKey key{seed_bin, pad1_bin, pad2_bin}; //LUTKey defined in header file + + if (lut_.find(key) != lut_.end()) { + + std::vector three_cluster_vec = { + seed, cluster1, cluster2}; + + ldmx::TrigScintTrack track = makeTrack(three_cluster_vec); + track_candidates.push_back(track); + } + } + } // over clusters in pad1, //LUT ends here + +//----------------------------------------------------------------------------- + + + // continue to next seed if 0 track candidates + if (track_candidates.size() == 0) continue; + + int keep_idx = 0; + float min_residual = 1000; // some large number + + // no need to choose between only one candidate track + if (track_candidates.size() > 1) { + // now for each seed, pick only the track with the smallest residual. + + if (verbose_) { + ldmx_log(debug) << "Got " << track_candidates.size() + << " tracks to check."; + } + + for (uint idx = 0; idx < track_candidates.size(); idx++) { + if ((track_candidates.at(idx)).getResidual() < min_residual) { + keep_idx = (int)idx; + min_residual = + (track_candidates.at(idx)).getResidual(); // update minimum + + if (verbose_ > 1) { + ldmx_log(debug) + << "Track at index " << idx + << " has smallest residual so far: " << min_residual; + } + + } // finding min residual + } // over track candidates + } // if more than 1 to choose from + + // store the track at keepIdx, if there was one we made it this far and + // keepIdx is 0 or has been updated to the smallest residual track idx + // if (keepIdx >= 0) { + tracks_.push_back(track_candidates.at(keep_idx)); + if (verbose_) { + ldmx_log(debug) << "Kept track at index " << keep_idx; + ldmx_log(trace) << track_candidates.at(keep_idx); + } + //} + } // over seeds + + // done here if there were no tracks found + if (tracks_.size() == 0) { + if (verbose_) { + ldmx_log(debug) << "No tracks found!"; + } + std::vector empty{}; + event.add(output_collection_, empty); + return; + } + // now, if there are multiple seeds sharing the same downstream hits, this + // should also be remedied with a selection on min residual. + + // The logic of this loop kind of assumes I can remove tracks immediately -- + // that way I can do pairwise checks between more tracks within a single + // loop. But for now I haven't figured out how to erase elements in a fool + // proof way. So I iterate over a vector... + + std::vector keep_indices(tracks_.size(), 1); + if (verbose_ > 1) + ldmx_log(debug) << "vector of indices to keep has size " + << keep_indices.size(); + + for (uint idx = tracks_.size() - 1; idx > 0; idx--) { + // since we start in one end, we only have to check matches in one + // direction + ldmx::TrigScintTrack track = tracks_.at(idx); + for (int idx_comp = idx - 1; idx_comp >= 0; idx_comp--) { + if (verbose_ > 1) + ldmx_log(debug) << "In track disambiguation loop, idx points at " + << idx << " and prev idx points at " << idx_comp; + + ldmx::TrigScintTrack next_track = tracks_.at(idx_comp); + + // no need to start pulling constituents from tracks that are + // ridiculously far apart + if (fabs(track.getCentroid() - next_track.getCentroid()) < + 3 * max_delta_) { + std::vector consts_1 = + track.getConstituents(); + std::vector consts_2 = + next_track.getConstituents(); + if (verbose_ > 1) + ldmx_log(debug) + << "In track disambiguation loop, got the two tracks, " + "with nConstituents " + << consts_1.size() << " and " << consts_2.size() + << ", respectively. "; + // let's do "if either cluster is shared" right now... but could also + // have it settable to use a stricter cut: an AND + if (consts_1[1].getCentroid() == consts_2[1].getCentroid() || + (consts_1.size() > 2 && consts_2.size() > 2 && + consts_1[2].getCentroid() == + consts_2[2] + .getCentroid())) { // we have overlap downstream of the + // seeding pad. probably, one cluster + // in seeding pad is noise + + if (verbose_ > 1) { + ldmx_log(debug) << "Found overlap! Tracks at index " << idx + << " and " << idx_comp; + ldmx_log(trace) << tracks_.at(idx); + ldmx_log(trace) << tracks_.at(idx_comp); + } + + if ((tracks_.at(idx)).getResidual() < + (tracks_.at(idx_comp)).getResidual()) { + // next track (lower index) is a worse choice, remove its flag for + // keeping + keep_indices.at(idx_comp) = 0; + } else // prefer next track over current. remove current track's + // keep + // flag + keep_indices.at(idx) = 0; + /*} + else { + tracks_.erase(itNext); + // removeIdx.push_back(idx+1); + // we might see the same index two times in the loop in this + case, if there are three seeds sharing the same clusters + downstream. + // then the third only gets removed if it's even worse than + the second. + // one could deal with this with an extra overlap check. not + sure we will be in this situation any time soon though. + }*/ + } // over matching/overlapping tracks + } // over tracks close enough to share constituents + } // over constructed tracks at other indices, to match + } // over constructed tracks + + for (uint idx = 0; idx < tracks_.size(); idx++) { + if (verbose_ > 1) { + ldmx_log(debug) << "keep flag for idx " << idx << " is " + << keep_indices.at(idx); + } + if (keep_indices.at(idx)) { // this hasn't been flagged for removal + + cleaned_tracks.push_back(tracks_.at(idx)); + + if (verbose_) { + ldmx_log(debug) << "After cleaning, keeping track at index " << idx + << ": Centroid = " << (tracks_.at(idx)).getCentroid() + << "; CentroidX = " + << (tracks_.at(idx)).getCentroidX() + << "; CentroidY = " + << (tracks_.at(idx)).getCentroidY() + << "; track PE = " << (tracks_.at(idx)).getPE() + << tracks_.at(idx); + } + } // if index flagged for keeping + } // over all (uniquely seeded) tracks in the event + + if (verbose_) { + for (uint idx = 0; idx < tracks_.size(); idx++) { + ldmx_log(debug) << "Keeping track at index " << idx << ":" + << tracks_.at(idx); + } + } + + if (verbose_) { + ldmx_log(debug) << "Running track x,y matching "; + } + + if (cleaned_tracks.size() > 0) { + matchXYTracks(cleaned_tracks); + std::vector matched_tracks = + cleaned_tracks; // don't know why this copying needs to happen but it + // does + // std::vector matchXYTracks( cleanedTracks + //); std::vector matchedTracks = + // matchXYTracks( cleanedTracks ); + for (auto trk : matched_tracks) { + /* for (uint idx = 0; idx < tracks_.size(); idx++) { + if (verbose_ > 1) { + ldmx_log(debug) << "keep flag for idx " << idx << " is " + << keepIndices.at(idx); + } + if (keepIndices.at(idx)) { // this hasn't + been flagged for removal + //check if channel nb is above that of horizontal bars + if (tracks_.at(idx).getCentroid() >= vert_bar_start_idx_) + */ + if (trk.getCentroid() >= vert_bar_start_idx_) + cleaned_tracks_x.push_back(trk); // acks_.at(idx)); + else + cleaned_tracks_y.push_back(trk); // acks_.at(idx)); + // cleanedTracksY.push_back(trk); + if (verbose_ > 1) { + float centr = trk.getCentroid(); // tracks_.at(idx).getCentroid(); // + std::string coll_str = centr >= vert_bar_start_idx_ ? "X" : "Y"; + coll_str = output_collection_ + coll_str; + ldmx_log(debug) << "saving track with centroid " << centr + << " to output track collection " << coll_str; + } + // } + } + } + + } // if there are clusters in all pads + + + + else if (verbose_) { //if theyre arent clusters in all pads + ldmx_log(info) + << "Not all pads had clusters; (maybe) skipping tracking attempt"; + } + + if (verbose_) { + ldmx_log(debug) << "Done with tracking step. "; + } + + event.add(output_collection_, cleaned_tracks); + // event.add(output_collection_, matchedTracks); + + event.add(output_collection_ + "Y", cleaned_tracks_y); + event.add(output_collection_ + "X", cleaned_tracks_x); + + tracks_.resize(0); + + return; +} //produce ends + + +//---------------------------------------------------------------------------- + + + +ldmx::TrigScintTrack TrigScintLUTTracker::makeTrack( + std::vector clusters) { + // for now let's keep a straight, unweighted centroid + // consider the possibility that at least one cluster has a centroid + // identically == 0. then we need to shift them by 1 if we want to do energy + // weighted track centroid later. but no need now + ldmx::TrigScintTrack tr; + float centroid = 0; + float centroid_x = 0; + float centroid_y = 0; + float beam_efrac = 0; + float pe = 0; + for (uint i = 0; i < clusters.size(); i++) { + centroid += (clusters.at(i)).getCentroid(); + centroid_x += (clusters.at(i)).getCentroidX(); + centroid_y += (clusters.at(i)).getCentroidY(); + tr.addConstituent(clusters.at(i)); + beam_efrac += (clusters.at(i)).getBeamEfrac(); + pe += (clusters.at(i)).getPE(); + } + centroid /= clusters.size(); + centroid_x /= clusters.size(); + if (centroid >= vert_bar_start_idx_) { + if (verbose_) { + ldmx_log(debug) + << " -- In makeTrack made vertical bar track with centroid " + << centroid << " and y flag sum " << centroid_y; + // try commenting this to check if that helps with an out-of-bounds + // problem + // << " from clusters with y centroids"; + // for (uint i = 0; i < clusters.size(); i++) + // ldmx_log(debug) << "\tpad " << i << ": centroidY " + // << (clusters.at(i)).getCentroidY(); + } + // then the sum of centroid y is 0, 2, 4 or 6 + // we have 4 divisions, so, the center of it should be divNb/8 + // (or rather, that's where channel nBars/8 begins) + // and then a factor 2 for the zig-zag pattern + centroid_y = (centroid_y + 1) * 2 * n_bars_y_ / 8.; + // TODO: here we could instead just use quadrant indices 0-3 by dividing by + // 2 but that would mean that in the raw, x and y track centroidY would mean + // different things + if (verbose_) ldmx_log(debug) << " -- new centroidY = " << centroid_y; + } else + centroid_y /= clusters.size(); + + beam_efrac /= clusters.size(); + pe /= clusters.size(); + + float residual = 0; + for (uint i = 0; i < clusters.size(); i++) + residual += ((clusters.at(i)).getCentroid() - centroid) * + ((clusters.at(i)).getCentroid() - centroid); + residual = sqrt(residual / clusters.size()); + + tr.setCentroid(centroid); + tr.setCentroidX(centroid_x); + tr.setCentroidY(centroid_y); + tr.setResidual(residual); + tr.setBeamEfrac(beam_efrac); + tr.setPE(pe); + + if (verbose_) { + ldmx_log(debug) << " -- In makeTrack made track with centroid " + << centroid << " and residual " << residual << " and pe " + << pe << " from clusters with centroids"; + for (uint i = 0; i < clusters.size(); i++) + ldmx_log(debug) << "\tpad " << i << ": centroid " + << (clusters.at(i)).getCentroid(); + } + + return tr; +} + + +// std::vector TrigScintTrackProducer::matchXYTracks( +void TrigScintLUTTracker::matchXYTracks( + std::vector &tracks) { + // map quadrant nb to track (can be multiple per quadrant) + std::multimap + y_idx_quad_map; // key = quad, val = track index in collection + std::multimap x_idx_quad_map; + + std::multimap y_quad_map; + std::multimap x_quad_map; + // map track in quadrant back to index in entire track collection + // used for updating collection track variables + std::map y_track_map; + std::map x_track_map; + + uint trk_idx = -1; + for (auto trk : tracks) { + trk_idx++; + // 1. get the y bar tracks with centroidX = -1 + if (trk.getCentroidX() == -1) { + if (verbose_) + ldmx_log(debug) << " -- In matchXYTracks found y track at " + << trk.getCentroidY() << "; mapping to quad " + << (int)trk.getCentroidY() / 8 << " with trk index " + << trk_idx; + // 2. order them... or map them to quadrants. note that there are 2 layers + // so 2*n_bars_y_/4 channels per quadrant + y_quad_map.insert(std::make_pair((int)(trk.getCentroidY() / 8), trk)); + y_track_map[trk] = trk_idx; + y_idx_quad_map.insert( + std::make_pair((int)(trk.getCentroidY() / 8), trk_idx)); + + } else { // 3. get the remaining tracks (from vertical bars) and map them + // (back) to (middle of) quadrants + x_quad_map.insert(std::make_pair((int)(trk.getCentroidY() / 8), trk)); + x_track_map[trk] = trk_idx; + x_idx_quad_map.insert( + std::make_pair((int)(trk.getCentroidY() / 8), trk_idx)); + if (verbose_) + ldmx_log(debug) << " -- In matchXYTracks found x track at (x,y) = (" + << trk.getCentroidX() << ", " << trk.getCentroidY() + << "); mapping to quad " << (int)trk.getCentroidY() / 8 + << " with trk index " << trk_idx; + } + } + + // 4a + // + // 1) here use the geometry? if we can assume perfect alignment we can take + // width and nBars and take nBars/2 as origin + // --- now do the matching --- + + // if there is no useful matching to be done: these are the pad width wide + // numbers + float x0 = 0; + // this should be half the pad... could also set + // it to full beam spot width + float sx0 = fabs(x_start_); + + // y_start_ is half the pad, so this should be half a quadrant + float sy0 = fabs(y_start_) / 4.; + + // assume at least one y track. will have to figure out if there is ever a + // reason to use an isolated x track in its place. + for (auto yitr = y_quad_map.begin(); yitr != y_quad_map.end(); ++yitr) { + int n_yin_quad = y_quad_map.count((*yitr).first); + int n_xin_quad = x_quad_map.count((*yitr).first); + float y{-9999.}, sy{-9999.}, x{-9999.}, x1{-9999.}, x2{-9999.}, sx1{-9999.}, + sx2{-9999.}, y1{-9999.}, y2{-9999.}, sy1{-9999.}, sy2{-9999.}; + // quad midpoint: + float y0 = (*yitr).first * 8 + sy0; + float sx = 1. / 2 * + x_conv_factor_; // rely on x precision being one single bar + // width; always used unless x is undeterminable + + // check all x first + // do the easiest first: + if (n_xin_quad == 0) { // then there's no hope of setting a better x here + // just use the beam spot width... and center of pad + x = x0; + sx = sx0; + if (verbose_) + ldmx_log(debug) << "\t\t\t no x info in quad " << (*yitr).first + << "; will set x to middle of pad, pad half-width as " + "precision: set (x, sx)=(" + << x << ", " << sx << ")"; + } // 0 x tracks in quadrant + else if (n_xin_quad == + 1) { // slightly harder: 1 x track -- might be easy if + // it's just one y track; if several, need to + // think about overlaps. but in overlap case, just + // revert to setting x0 and sx0, when we know + auto xitr = x_quad_map.find((*yitr).first); + x = ((*xitr).second).getCentroidX() * x_conv_factor_ + x_start_; + + if (verbose_) + ldmx_log(debug) << "\t\t\t 1 x in quad " << (*yitr).first + << ", getting (x, sx)=(" << x << ", " << sx << ")"; + } // 1 x track in quadrant + else if (n_xin_quad == 2) { // finally if we have two tracks, get x1 and x2 + // and decide later how to use them + // don't think we want to experiment with discerning three overlapping + // tracks, so not >= 2 + // continue; //debugging: skip for now -- didn't help + auto xitr1 = x_quad_map.lower_bound((*yitr).first); + auto xitr2 = x_quad_map.upper_bound((*yitr).first); + xitr2--; // upper_bound points to next element + + if (xitr1 != xitr2) { // should be true already but... + x1 = ((*xitr1).second).getCentroidX() * x_conv_factor_ + x_start_; + x2 = ((*xitr2).second).getCentroidX() * x_conv_factor_ + x_start_; + sx1 = x_conv_factor_ / 2.; // 1 bar width + sx2 = sx1; + x = (x1 + x2) / 2.; + sx = fabs(x1 - x2) / 2 * + x_conv_factor_; // rely on x precision being one single pad width + if (verbose_) + ldmx_log(debug) << "\t\t -- 2 x in quad: setting y track x " + "coordinate to midpoint"; + } + } // if 2 x tracks in quad + + // ok! over y: + // can skip 0 y case by construction + if (n_yin_quad == 1) { // we can already now tell what the y coordinate and + // its precision is + y = ((*yitr).second).getCentroidY() * y_conv_factor_ + y_start_; + sy = ((*yitr).second).getResidual() * y_conv_factor_; + // if all clusters lined up, assign + // precision of 1 bar width + if (sy == 0) sy = 1. / 2 * y_conv_factor_; + + if (n_xin_quad <= 1) { + // 4. every quadrant which just has one of each --> done ; + // b) set the sx, sy of the x track now, using the residuals from the + // other b1) special case: no x tracks; then x, sx have been set above + auto xidx = x_idx_quad_map.find((*yitr).first); + auto yidx = y_idx_quad_map.find((*yitr).first); + tracks.at((*xidx).second).setPosition(x, y); + tracks.at((*xidx).second).setSigmaXY(sx, sy); + tracks.at((*yidx).second).setPosition(x, y); + tracks.at((*yidx).second).setSigmaXY(sx, sy); + if (verbose_) + ldmx_log(debug) << "\t\t\t in quad " << (*yitr).first + << ", set (x, y) = (" << x << ", " << y + << ") and (sx, sy) = " << sx << ", " << sy << ")"; + continue; + } + } // 1 y, 0 or 1 x track in quadrant + + if (verbose_) + ldmx_log(debug) << "\t\t in quad " << (*yitr).first + << ", not single x,y tracks: " << n_xin_quad + << " of x and " << n_yin_quad << " of y"; + + if (n_yin_quad == 2) { // let's start here and see if we can do >= 2 later + // here one could do sth to avoid checking the other y track again in the + // outermost loop over y + auto yitr1 = y_quad_map.lower_bound((*yitr).first); + auto yitr2 = y_quad_map.upper_bound((*yitr).first); + yitr2--; // back up once + y1 = ((*yitr1).second).getCentroidY() * y_conv_factor_ + y_start_; + y2 = ((*yitr2).second).getCentroidY() * y_conv_factor_ + y_start_; + sy1 = ((*yitr1).second).getResidual() * y_conv_factor_; + sy2 = ((*yitr2).second).getResidual() * y_conv_factor_; + y = (y1 + y2) / 2.; + sy = fabs(y1 - y2) / 2 * y_conv_factor_; + if (verbose_) + ldmx_log(debug) + << "\t\t -- 2 y in quad: setting x track y coordinate to midpoint"; + } // 2y in quad + + if (n_yin_quad == 1 && + n_xin_quad == 2) { // don't think we want to experiment with discerning + // three overlapping tracks, so not >= 2 + + // first: set the y track coordinates to x = the mid of x tracks, y = y + // of y track + auto yidx = y_idx_quad_map.find((*yitr).first); + tracks.at((*yidx).second).setPosition(x, y); + tracks.at((*yidx).second).setSigmaXY(sx, sy); + + int min_overlap_pe = 250; + if (((*yitr).second).getPE() < min_overlap_pe) { + // can't tell, really, that either of these belong to the y track. so. + // let them keep their own x coordinate but set y to quadrant midpoint, + // with uncertainty +/- half quadrant width (1/8 of pad height) + y = y0; + sy = sy0; + if (verbose_) + ldmx_log(debug) << "\t\t -- Can't tell which x track should be " + "matched to single y track. Setting both x track " + "coordinates to y quadrant value:"; + } // if can't assume overlap + else if (verbose_) + ldmx_log(debug) << "\t\t -- Found large PE count (" + << ((*yitr).second).getPE() << " > " << min_overlap_pe + << "), suggesting overlap! Setting both x track " + "coordinates to y track value:"; + + // consider making two x tracks out if this one, and, anyway have to set + // their average as the y track x cocordinate + // EXPERIMENTAL : apply only to x tracks, which can be disregarded for + // electron counting + if (verbose_) + ldmx_log(debug) << "\t\t -- (x1, x2, y) = (" << x1 << ", " << x2 + << ", " << y << ") and (sx1, sx2, sy) = " << sx1 << ", " + << sx2 << ", " << sy << ")"; + + // now set x track coordinates according to overlap check result + auto xidx1 = x_idx_quad_map.lower_bound((*yitr).first); + auto xidx2 = x_idx_quad_map.upper_bound((*yitr).first); + xidx2--; // upper_bound points to (last+1) element + tracks.at((*xidx1).second).setPosition(x1, y); + tracks.at((*xidx1).second).setSigmaXY(sx1, sy); + tracks.at((*xidx2).second).setPosition(x2, y); + tracks.at((*xidx2).second).setSigmaXY(sx2, sy); + + } // 1 y, 2 x tracks in the quadrant + else if (n_yin_quad == 2 && n_xin_quad == 1) { + // 5b) if there are more y than x: could be an overlap + + // first: set the x track coordinates to x = x of x track, y = the mid of + // y tracks + auto xidx = x_idx_quad_map.find((*yitr).first); + tracks.at((*xidx).second).setPosition(x, y); + tracks.at((*xidx).second).setSigmaXY(sx, sy); + + auto xitr = x_quad_map.lower_bound((*yitr).first); + int min_overlap_pe = 300; + if (((*xitr).second).getPE() < min_overlap_pe) { + if (verbose_) + ldmx_log(debug) + << "\t\t just 1 x track with not-unusual PE in the quad -- can't " + "match; setting mid-point values for x "; + x = x0; + sx = sx0; + } // if can't assume overlap + else { + // consider making two x tracks out if this one, and, anyway have to set + // their average as the y track x cocordinate + // EXPERIMENTAL : apply only to x tracks, which can be disregarded for + // electron counting + if (verbose_) + ldmx_log(debug) << "\t\t -- Found large PE count (" + << ((*xitr).second).getPE() << " > " << min_overlap_pe + << ") in x track, suggesting overlap! Setting both y " + "track coordinates to x track value:"; + } // if can assume overlap + if (verbose_) + ldmx_log(debug) << "\t\t -- (x, y1, y2) = (" << x << ", " << y1 << ", " + << y2 << ") and (sx, sy1, sy2) = " << sx << ", " << sy1 + << ", " << sy2 << ")"; + + auto yidx1 = y_idx_quad_map.lower_bound((*yitr).first); + auto yidx2 = y_idx_quad_map.upper_bound((*yitr).first); + yidx2--; // upper_bound points to next element + tracks.at((*yidx1).second).setPosition(x, y1); + tracks.at((*yidx1).second).setSigmaXY(sx, sy1); + tracks.at((*yidx2).second).setPosition(x, y2); + tracks.at((*yidx2).second).setSigmaXY(sx, sy2); + + } // 2 y and 1 x track in quad + else if (n_yin_quad == 2 && n_xin_quad == 2) { + // MIDPONTS ALL OVER! + auto xidx1 = x_idx_quad_map.lower_bound((*yitr).first); + auto xidx2 = x_idx_quad_map.upper_bound((*yitr).first); + xidx2--; + auto yidx1 = y_idx_quad_map.lower_bound((*yitr).first); + auto yidx2 = y_idx_quad_map.upper_bound((*yitr).first); + yidx2--; + + if (y_idx_quad_map.find((*yitr).first) == y_idx_quad_map.end()) + ldmx_log(error) << "The two y tracks in the same quadrant at " + << (*yitr).first + << " appear to not be found in the y track map! " + "investigate. Note that yidx1.first = " + << (*yidx1).first + << " and yidx2.first = " << (*yidx2).first; + else { + tracks.at((*xidx1).second).setPosition(x1, y); + tracks.at((*xidx1).second).setSigmaXY(sx1, sy); + tracks.at((*xidx2).second).setPosition(x2, y); + tracks.at((*xidx2).second).setSigmaXY(sx2, sy); + + tracks.at((*yidx1).second).setPosition(x, y1); + tracks.at((*yidx1).second).setSigmaXY(sx, sy1); + tracks.at((*yidx2).second).setPosition(x, y2); + tracks.at((*yidx2).second).setSigmaXY(sx, sy2); + + if (verbose_) + ldmx_log(debug) << "\t\t -- in a 2 x 2 situaiton; midpoint y: " << y + << " for both x tracks, midpoint x: " << x + << " for both y tracks"; + } + } // if 2 y, 2 x tracks + + if (n_xin_quad > 2) { + if (verbose_) + ldmx_log(debug) << "\t\t -*-*-*- more than 2 x tracks in the same quad " + "-- nothing done about the x,y coordinates in this " + "situation -- implement if needed!!"; + } + if (n_yin_quad > 2) { + if (verbose_) + ldmx_log(debug) << "\t\t -*-*-*- more than 2 y tracks in the same quad " + "-- nothing done about the x,y coordinates in this " + "situation -- implement if needed!!"; + } + + } // over y tracks + + y_quad_map.clear(); + x_quad_map.clear(); + + // return tracks; +} + +void TrigScintLUTTracker::onProcessStart() { + ldmx_log(debug) << "Process starts!"; + + return; +} + +void TrigScintLUTTracker::onProcessEnd() { + ldmx_log(debug) << "Process ends!"; + + return; +} + +} // namespace trigscint + +DECLARE_PRODUCER(trigscint::TrigScintLUTTracker); + From f17ac9ec1a73e628a169bb99c2960342b176c0a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Mon, 20 Apr 2026 16:51:02 +0200 Subject: [PATCH 03/20] Add LUT method tracking (continued) --- ClusterViewerAnalyzer.cxx | 49 ++++++++++++------ ClusterViewerFig.py | 6 ++- LUTAnalyzer.cxx | 101 ++++++++++++++++++++++++++++++++++++++ LUTconfig.py | 10 ++++ NewEventcreator.py | 8 +-- TrackFig.py | 63 ++++++++++++++++++++++++ extruthfig.py | 65 +++++++++++------------- notruthfig.py | 53 ++++++++++++++++++++ 8 files changed, 297 insertions(+), 58 deletions(-) create mode 100644 LUTAnalyzer.cxx create mode 100644 LUTconfig.py create mode 100644 TrackFig.py create mode 100644 notruthfig.py diff --git a/ClusterViewerAnalyzer.cxx b/ClusterViewerAnalyzer.cxx index 466dc6c329..7e6d344bd4 100644 --- a/ClusterViewerAnalyzer.cxx +++ b/ClusterViewerAnalyzer.cxx @@ -17,9 +17,9 @@ class ClusterViewerAnalyzer : public framework::Analyzer { outfile.open("clusters.txt"); } - + void analyze(const framework::Event& event) override { - + const auto& clusters1 = event.getCollection("TriggerPad1Clusters","truthisolater"); @@ -29,24 +29,43 @@ class ClusterViewerAnalyzer : public framework::Analyzer { const auto& clusters3 = event.getCollection("TriggerPad3Clusters","truthisolater"); - std::cout << "Pad1 clusters: " << clusters1.size() << std::endl; - std::cout << "Pad2 clusters: " << clusters2.size() << std::endl; - std::cout << "Pad3 clusters: " << clusters3.size() << std::endl; - size_t Clusters = std::min({clusters1.size(), clusters2.size(), clusters3.size()}); - - for (size_t i=0; i < Clusters; i++) { + for (size_t i = 0; i < Clusters; i++) { + + int eventNumber = event.getEventNumber(); float c1 = clusters1[i].getCentroid(); float c2 = clusters2[i].getCentroid(); float c3 = clusters3[i].getCentroid(); - std::cout << "printed " << c1 << std::endl; - std::cout << "printed " << c2 << std::endl; - std::cout << "printed " << c3 << std::endl; - outfile << c1 << " " - << c2 << " " - << c3 << "\n"; - } + float c3_alt = c3; + + if (clusters3.size() > i + 1) { + c3_alt = clusters3[i + 1].getCentroid(); + } + + outfile << eventNumber << " " + << c1 << " " + << c2 << " " + << c3 << " " + << c3_alt << "\n"; + } + + + for (size_t i = 1; i < clusters1.size(); i++){ + std::cout << "In event " << event.getEventNumber() <<", extra cluster in pad 1: " + << clusters1[i].getCentroid() << "\n"; + } + + for (size_t i = 1; i < clusters2.size(); i++){ + std::cout << "In event " << event.getEventNumber() <<", extra cluster in pad 2: " + << clusters2[i].getCentroid() << "\n"; + } + + for (size_t i = 1; i < clusters3.size(); i++){ + std::cout << "In event " << event.getEventNumber() <<", extra cluster in pad 3: " + << clusters3[i].getCentroid() << "\n"; + } + } void onProcessEnd() override {outfile.close();} }; diff --git a/ClusterViewerFig.py b/ClusterViewerFig.py index b003c6dc47..e7f2861c2a 100644 --- a/ClusterViewerFig.py +++ b/ClusterViewerFig.py @@ -2,6 +2,8 @@ from LDMX.Framework import ldmxcfg p = ldmxcfg.Process('clusterviewer') -p.input_files = ["TruthEvents10.root"] #change for new naming conventions +#p.max_events = 10 +p.input_files = ["TruthEvents1000.root"] #change for new naming conventions #p.ampl_weighting = False -p.sequence = [ldmxcfg.Analyzer.from_file('ClusterViewerAnalyzer.cxx')] +p.sequence = [ldmxcfg.Analyzer.from_file('ClusterViewerAnalyzer.cxx', + needs = ['TrigScint_Event', 'SimCore_Event'])] diff --git a/LUTAnalyzer.cxx b/LUTAnalyzer.cxx new file mode 100644 index 0000000000..63040970ab --- /dev/null +++ b/LUTAnalyzer.cxx @@ -0,0 +1,101 @@ + +#include "Framework/EventProcessor.h" +#include "Framework/Event.h" +#include "TrigScint/Event/TrigScintTrack.h" +#include +#include +#include +#include + +class LUTAnalyzer : public framework::Analyzer { + public: + + struct Line {int event; float p1, p2, p3, p3alt; + }; + + std::ifstream infile; + std::ofstream outfile; + + std::map, std::vector> groups; + + int totalLines = 0; + + LUTAnalyzer(const std::string& name, framework::Process& p) + : framework::Analyzer(name,p) { + + infile.open("clusters.txt"); + outfile.open("LUT.txt"); + } + + void analyze(const framework::Event& event) override { + if (totalLines > 0) return; + + int ev; + float p1, p2, p3, p3alt; + + while (infile >> ev >> p1 >> p2 >> p3 >> p3alt) { + + float p12 = p2-p1; + float p23 = p3-p2; + + groups[{p12,p23}].push_back({ev,p1,p2,p3,p3alt}); + totalLines++; + } + } + + void onProcessEnd() override { + + int combs = 0; + int tracks = 0; + + std::cout << "Total number of track candidates: " << totalLines << "\n"; + std::cout << "Number of track candidate types: " << groups.size() << "\n"; + + double threshold = 1.0/1000.0; + std::cout << "Threshold: " << threshold << "\n"; + + for (auto& g : groups) { + + float p12 = g.first.first; + float p23 = g.first.second; + int count = g.second.size(); + + double frac = (double)count / totalLines; + + std::cout << "(" << p12 << "," << p23 << ") appears " + << count << " times, (" << frac << " %)" << "\n"; + } + + for (auto& g : groups) { + + int count = g.second.size(); + double frac = (double)count / totalLines; + + if (frac > threshold) { + + combs++; + + for (auto& line : g.second) { + + tracks++; + + outfile << line.p1 << " " + << line.p2 << " " + << line.p3 << "\n"; + } + } + } + std::cout << "\nLUT textfile written."; + std::cout << "\n" << combs << " combinations (" << tracks << " tracks) written to LUT." << "\n"; + } +}; + +DECLARE_ANALYZER(LUTAnalyzer) + + + + + + + + diff --git a/LUTconfig.py b/LUTconfig.py new file mode 100644 index 0000000000..166275734f --- /dev/null +++ b/LUTconfig.py @@ -0,0 +1,10 @@ +#config file to make LUT for tracking + +from LDMX.Framework import ldmxcfg + +p = ldmxcfg.Process("LUTmaker") +p.input_files = ["TruthSamples.root"] #necessary to define but it doesnt use it +#p.output_files = ["Tracks.root"] + +p.sequence = [ldmxcfg.Analyzer.from_file('LUTAnalyzer.cxx', needs = ['TrigScint_Event', 'SimCore_Event']) + ] diff --git a/NewEventcreator.py b/NewEventcreator.py index 8e3640c546..8c6d9e50ed 100644 --- a/NewEventcreator.py +++ b/NewEventcreator.py @@ -1,7 +1,7 @@ #Creates simulation sample .root file with TS info from LDMX.Framework import ldmxcfg -p = ldmxcfg.Process('test') +p = ldmxcfg.Process('simulation') from LDMX.SimCore import simulator as sim from LDMX.SimCore import generators as gen @@ -13,9 +13,9 @@ mySim.description = 'Basic test Simulation' p.sequence = [ mySim ] -p.run = 1 -p.max_events = 10 #new number of electrons -p.output_files = [ 'NewEvents10.root' ] #new output file name +p.run = 1 #different for each simulation +p.max_events = 10000 +p.output_files = [ 'SimSamples.root' ] #new output file name import LDMX.Ecal.ecal_geometry import LDMX.Ecal.ecal_hardcoded_conditions diff --git a/TrackFig.py b/TrackFig.py new file mode 100644 index 0000000000..d1384ee13a --- /dev/null +++ b/TrackFig.py @@ -0,0 +1,63 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Apr 13 12:17:11 2026 + +@author: kvarn +""" + +from LDMX.Framework import ldmxcfg + +p = ldmxcfg.Process('tracking') +p.input_files = ["TruthSamples.root"] +p.output_files = ["Tracks.root"] + +#from LDMX.TrigScint.trig_scint import TrigScintTrackProducer + +prod_tracks = ldmxcfg.Producer("lut_track", "trigscint::TrigScintTrackProducer", + "TrigScint") + +LUT_tracks = ldmxcfg.Producer("prod_track", "trigscint::TrigScintLUTTracker", + "TrigScint") + +prod_tracks.seeding_collection = "TriggerPad1Clusters" +prod_tracks.further_input_collections = ["TriggerPad2Clusters", "TriggerPad3Clusters"] +prod_tracks.delta_max = 1.0 +prod_tracks.input_pass_name = "" +prod_tracks.verbosity = 1 +prod_tracks.vertical_bar_start_index = 100 +prod_tracks.number_horizontal_bars = 16 +prod_tracks.horizontal_bar_width = 10.0 +prod_tracks.horizontal_bar_gap = 1.0 +prod_tracks.number_vertical_bars = 16 +prod_tracks.vertical_bar_width = 10.0 +prod_tracks.vertical_bar_gap = 1.0 +prod_tracks.allow_skip_last_collection = False + +p.logger.term_level = 1 + +LUT_tracks.seeding_collection = "TriggerPad1Clusters" +LUT_tracks.further_input_collections = ["TriggerPad2Clusters", "TriggerPad3Clusters"] +LUT_tracks.delta_max = 1.0 +LUT_tracks.input_pass_name = "" +LUT_tracks.verbosity = 1 +LUT_tracks.vertical_bar_start_index = 100 +LUT_tracks.number_horizontal_bars = 16 +LUT_tracks.horizontal_bar_width = 10.0 +LUT_tracks.horizontal_bar_gap = 1.0 +LUT_tracks.number_vertical_bars = 16 +LUT_tracks.vertical_bar_width = 10.0 +LUT_tracks.vertical_bar_gap = 1.0 +LUT_tracks.allow_skip_last_collection = False +LUT_tracks.lut_file = "LUT.txt" + + +prod_tracks.output_collection = "TrigScintTrackProdTracks" +LUT_tracks.output_collection = "TrigScintLUTTracks" + + +p.sequence = [#ldmxcfg.Analyzer.from_file('LUTAnalyzer.cxx', + #needs = ['TrigScint_Event', 'SimCore_Event']), + prod_tracks, + #LUT_tracks + ] \ No newline at end of file diff --git a/extruthfig.py b/extruthfig.py index 883fafe7e6..9bac416ffd 100644 --- a/extruthfig.py +++ b/extruthfig.py @@ -1,21 +1,20 @@ -#example config for truth hit producer - -#Uses TruthHitProducer to isolate true hits, -#(then turns them to digis, then clusters) +#Uses TruthHitProducer to isolate true hits, +#then turns them to digis, then clusters, then makes clusters.txt from LDMX.Framework import ldmxcfg p = ldmxcfg.Process('truthisolater') -p.input_files = ["NewEvents10.root"] -p.output_files = ["TruthEvents10.root"] - +p.input_files = ["SimSamples.root"] +p.output_files = ["TruthSamples.root"] +#an additional output file called clusters.txt will be created as well + from LDMX.TrigScint.truth_hits import TruthHitProducer from LDMX.TrigScint.trig_scint import TrigScintDigiProducer truth_hits = [TruthHitProducer('beamElectronsPad1'), TruthHitProducer('beamElectronsPad2'), TruthHitProducer('beamElectronsPad3') - ] + ] pad_num = 1 @@ -23,38 +22,30 @@ hits.input_collection=f"TriggerPad{pad_num}SimHits" hits.output_collection=f"truthBeamElectronsPad{pad_num}" pad_num+=1 - -truth_digis = [TrigScintDigiProducer.pad1(), - TrigScintDigiProducer.pad2(), - TrigScintDigiProducer.pad3() - ] - -for digi,hits in zip(truth_digis, truth_hits): - digi.input_collection = hits.output_collection - -from LDMX.TrigScint.trig_scint import TrigScintClusterProducer - -truth_clusters = [ - TrigScintClusterProducer.pad1(), - TrigScintClusterProducer.pad2(), - TrigScintClusterProducer.pad3(), - ] - -for cluster, digi in zip(truth_clusters, truth_digis): - cluster.input_collection = digi.output_collection - cluster.ampl_weighting = True +digis = [TrigScintDigiProducer.pad1(), + TrigScintDigiProducer.pad2(), + TrigScintDigiProducer.pad3() + ] -p.sequence = [*truth_hits, - *truth_digis, - *truth_clusters, - #trigScintTrack, - #count, -# TriggerProcessor('trigger', 800.), - # ldmxcfg.Producer.from_file('TruthHitProducer.cxx') - ] +for digi,hits in zip(digis, truth_hits): + digi.input_collection = hits.output_collection +from LDMX.TrigScint.trig_scint import TrigScintClusterProducer +clusters = [TrigScintClusterProducer.pad1(), + TrigScintClusterProducer.pad2(), + TrigScintClusterProducer.pad3(), + ] +for cluster, digi in zip(clusters, digis): + cluster.input_collection = digi.output_collection + cluster.ampl_weighting = False + cluster.clustering_threshold = 3.0 - +p.sequence = [*truth_hits, + *digis, + *clusters, + ldmxcfg.Analyzer.from_file('ClusterViewerAnalyzer.cxx', + needs = ['TrigScint_Event', 'SimCore_Event']), + ] diff --git a/notruthfig.py b/notruthfig.py new file mode 100644 index 0000000000..442824e079 --- /dev/null +++ b/notruthfig.py @@ -0,0 +1,53 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Tue Apr 14 21:58:25 2026 + +@author: kvarn +""" + +#digis, then clusters, then makes clusters.txt + +from LDMX.Framework import ldmxcfg + +p = ldmxcfg.Process('truthisolater') +p.input_files = ["SimSamples.root"] +p.output_files = ["AllSamples.root"] +#an additional output file called clusters.txt will be created as well + +from LDMX.TrigScint.trig_scint import TrigScintDigiProducer + +pad_num = 1 + +digis = [TrigScintDigiProducer.pad1(), + TrigScintDigiProducer.pad2(), + TrigScintDigiProducer.pad3() + ] + +for digi in digis: + digi.input_collection = f"TriggerPad{pad_num}SimHits" + pad_num+=1 + +from LDMX.TrigScint.trig_scint import TrigScintClusterProducer + +clusters = [ + TrigScintClusterProducer.pad1(), + TrigScintClusterProducer.pad2(), + TrigScintClusterProducer.pad3(), + ] + +for cluster, digi in zip(clusters, digis): + cluster.input_collection = digi.output_collection + cluster.ampl_weighting = False + cluster.clustering_threshold = 3.0 + +p.sequence = [ + #*truth_hits, + *digis, + *clusters, + ldmxcfg.Analyzer.from_file('ClusterViewerAnalyzer.cxx', + needs = ['TrigScint_Event', 'SimCore_Event']), + #tracks + ] + + From f567ba886ab29ccc88b14e890dde25fe1103f030 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Fri, 24 Apr 2026 21:43:24 +0200 Subject: [PATCH 04/20] Move away from amplitude weighted TS cluster centroids --- TrigScint/include/TrigScint/TrigScintClusterProducer.h | 6 +++--- TrigScint/python/trig_scint.py | 1 + TrigScint/src/TrigScint/TrigScintClusterProducer.cxx | 5 +++-- 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/TrigScint/include/TrigScint/TrigScintClusterProducer.h b/TrigScint/include/TrigScint/TrigScintClusterProducer.h index f58f18faee..4426c76ee9 100644 --- a/TrigScint/include/TrigScint/TrigScintClusterProducer.h +++ b/TrigScint/include/TrigScint/TrigScintClusterProducer.h @@ -45,9 +45,6 @@ class TrigScintClusterProducer : public framework::Producer { // collection of clusters produced std::vector clusters_; - //amplitude weighting for centroid values - bool ampl_weighting{true}; - // cluster seeding threshold double seed_{0.}; @@ -57,6 +54,9 @@ class TrigScintClusterProducer : public framework::Producer { // max number of neighboring hits to combine when forming a cluster int max_width_{2}; + //PE-amplitude weighting for centroid values + bool ampl_weighting_{true}; + // specific verbosity of this producer int verbose_{0}; diff --git a/TrigScint/python/trig_scint.py b/TrigScint/python/trig_scint.py index f87bbbf05e..57c5cc4255 100644 --- a/TrigScint/python/trig_scint.py +++ b/TrigScint/python/trig_scint.py @@ -280,6 +280,7 @@ def __init__(self,name) : super().__init__(name,'trigscint::TrigScintClusterProducer','TrigScint') self.max_cluster_width = 2 + self.ampl_weighting = True self.clustering_threshold = 0. #to add in neighboring channels self.seed_threshold = 30. self.pad_time = 0. diff --git a/TrigScint/src/TrigScint/TrigScintClusterProducer.cxx b/TrigScint/src/TrigScint/TrigScintClusterProducer.cxx index 93dfc88cac..af0bb294db 100644 --- a/TrigScint/src/TrigScint/TrigScintClusterProducer.cxx +++ b/TrigScint/src/TrigScint/TrigScintClusterProducer.cxx @@ -7,10 +7,10 @@ namespace trigscint { void TrigScintClusterProducer::configure(framework::config::Parameters &ps) { - ampl_weighting = ps.get("ampl_weighting"); seed_ = ps.get("seed_threshold"); min_thr_ = ps.get("clustering_threshold"); max_width_ = ps.get("max_cluster_width"); + ampl_weighting_ = ps.get("ampl_weighting"); input_collection_ = ps.get("input_collection"); pass_name_ = ps.get("input_pass_name"); output_collection_ = ps.get("output_collection"); @@ -23,6 +23,7 @@ void TrigScintClusterProducer::configure(framework::config::Parameters &ps) { ldmx_log(info) << "Got parameters: \nSeed threshold: " << seed_ << "\nClustering threshold: " << min_thr_ << "\nMax cluster width: " << max_width_ + << "\nAmplitude weighting: " << ampl_weighting_ << "\nExpected pad hit time: " << pad_time_ << "\nMax hit time delay: " << time_tolerance_ << "\nVertical bar start index: " << vert_bar_start_idx_ @@ -481,7 +482,7 @@ void TrigScintClusterProducer::produce(framework::Event &event) { void TrigScintClusterProducer::addHit(uint idx, ldmx::TrigScintHit hit) { float ampl = hit.getPE(); float w = 1; - if (ampl_weighting) { + if (ampl_weighting_) { //if choosing to PE-weight centroid positions w = ampl; } From 9fa0fda99a106c28159e5c0e1bcf7f40049ef305 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Sun, 26 Apr 2026 14:32:22 +0200 Subject: [PATCH 05/20] Add LUT tracking as variable in TSTrackProducer --- .../TrigScint/TrigScintTrackProducer.h | 23 ++ TrigScint/python/trig_scint.py | 24 +- .../src/TrigScint/TrigScintTrackProducer.cxx | 210 +++++++++++------- 3 files changed, 155 insertions(+), 102 deletions(-) diff --git a/TrigScint/include/TrigScint/TrigScintTrackProducer.h b/TrigScint/include/TrigScint/TrigScintTrackProducer.h index 42667a3e30..4382ffc6e2 100644 --- a/TrigScint/include/TrigScint/TrigScintTrackProducer.h +++ b/TrigScint/include/TrigScint/TrigScintTrackProducer.h @@ -10,6 +10,8 @@ #include "TrigScint/Event/TrigScintCluster.h" #include "TrigScint/Event/TrigScintTrack.h" +#include + namespace trigscint { /** @@ -61,6 +63,9 @@ class TrigScintTrackProducer : public framework::Producer { // allow forming tracks without match in the last collection bool skip_last_{false}; + + // do tracking using LUT method instead of with max_delta + bool lut_tracking_{false}; // vertical bar start index int vert_bar_start_idx_{52}; @@ -84,6 +89,24 @@ class TrigScintTrackProducer : public framework::Producer { // track residual in units of channel nb (will not be content weighted) // float residual_{0.}; + + struct LUTKey { + float p1, p2, p3; + + bool operator==(const LUTKey& other) const { + return p1 == other.p1 && p2 == other.p2 && p3 == other.p3; + } + }; + + struct LUTKeyHash { + size_t operator()(const LUTKey& k) const { + return std::hash()(k.p1) ^ + (std::hash()(k.p2) << 1) ^ + (std::hash()(k.p3) << 2); + } + }; + + std::unordered_set lut_; float bar_width_y_{3.}; // mm float bar_gap_y_{2.1}; // mm diff --git a/TrigScint/python/trig_scint.py b/TrigScint/python/trig_scint.py index 57c5cc4255..393fcd8e9d 100644 --- a/TrigScint/python/trig_scint.py +++ b/TrigScint/python/trig_scint.py @@ -335,6 +335,8 @@ def __init__(self,name) : self.seeding_collection = "TriggerPad1Clusters" self.further_input_collections = ["TriggerPad2Clusters","TriggerPad3Clusters"] self.allow_skip_last_collection = False + self.lut_tracking = False + self.lut_file = "LUT.txt" self.vertical_bar_start_index = 52 self.number_horizontal_bars = 24 #16 for x,y segmented geometry only self.number_vertical_bars = 0 #8 for x,y segmented geometry only @@ -348,28 +350,6 @@ def __init__(self,name) : trig_scint_track = TrigScintTrackProducer( "trig_scint_track" ) -class TrigScintTrackProducer(ldmxcfg.Producer) : - """Configuration for track producer for Trigger Scintillators""" - - def __init__(self,name) : - super().__init__(name,'trigscint::','TrigScint') - - self.delta_max = 0.75 - self.tracking_threshold = 0. #to add in neighboring channels - self.seeding_collection = "TriggerPad1Clusters" - self.further_input_collections = ["TriggerPad2Clusters","TriggerPad3Clusters"] - self.allow_skip_last_collection = False - self.vertical_bar_start_index = 52 - self.number_horizontal_bars = 24 #16 for x,y segmented geometry only - self.number_vertical_bars = 0 #8 for x,y segmented geometry only - self.horizontal_bar_width = 3. - self.horizontal_bar_gap = 0.3 - self.vertical_bar_width = 3. - self.vertical_bar_gap = 0.3 - self.input_pass_name="" #take any pass - self.output_collection="TriggerPadTracks" - self.verbosity = 0 - class TrigScintFirmwareTracker(ldmxcfg.Producer) : """Configuration for the track producer from the Firmware Tracker""" def __init__(self,name) : diff --git a/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx b/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx index b17d2b7613..069ecdf58f 100644 --- a/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx +++ b/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx @@ -2,6 +2,8 @@ #include // std::next #include +#include +#include namespace trigscint { @@ -14,6 +16,7 @@ void TrigScintTrackProducer::configure(framework::config::Parameters &ps) { "further_input_collections"); // {"TriggerPadUpClusters" , // "TriggerPadDnClusters" } output_collection_ = ps.get("output_collection"); + pass_name_ = ps.get("input_pass_name"); verbose_ = ps.get("verbosity"); vert_bar_start_idx_ = ps.get("vertical_bar_start_index"); @@ -24,6 +27,8 @@ void TrigScintTrackProducer::configure(framework::config::Parameters &ps) { bar_width_x_ = ps.get("vertical_bar_width"); bar_gap_x_ = ps.get("vertical_bar_gap"); skip_last_ = ps.get("allow_skip_last_collection"); + lut_tracking_ = ps.get("lut_tracking"); //if choosing LUT tracking + std::string lut_file_ = ps.get("lut_file"); //from LUTAnalyzer // TO DO: allow any number of input collections @@ -36,6 +41,8 @@ void TrigScintTrackProducer::configure(framework::config::Parameters &ps) { << "\nInput pass name: " << pass_name_ << "\nAllow tracks with no hit in last collection: " << skip_last_ + << "\nUsing LUT Tracking Method: " << lut_tracking_ + << "\nIf using LUT Method, LUT from: " << lut_file_ << "\nVertical bar start index: " << vert_bar_start_idx_ << "\nNumber of horizontal bars: " << n_bars_y_ << "\nHorizontal bar width: " << bar_width_y_ @@ -54,7 +61,23 @@ void TrigScintTrackProducer::configure(framework::config::Parameters &ps) { x_conv_factor_ = bar_width_x_ + bar_gap_x_; // half width of pad x_start_ = -(n_bars_x_ * (bar_width_x_ + bar_gap_x_) - bar_gap_x_) / 2.; - + + if (lut_tracking_) { + + std::ifstream file(lut_file_); + float a, b, c; + + while (file >> a >> b >> c) { + float p1 = a; + float p2 = b; + float p3 = c; + lut_.insert({p1, p2, p3}); + } + + ldmx_log(info) << "Loaded LUT with size: " << lut_.size(); + + } + return; } @@ -149,107 +172,134 @@ void TrigScintTrackProducer::produce(framework::Event &event) { // reset for each seed // bool madeTrack = false; - - for (const auto &cluster1 : clusters_pad1) { - if (verbose_ > 1) { - ldmx_log(debug) << "\tGot pad1 cluster with centroid " - << cluster1.getCentroid(); - } - if (fabs(cluster1.getCentroid() - centroid) < max_delta_ || - (centroid >= vert_bar_start_idx_ // then in vertical bars - && seed.getCentroidX() == cluster1.getCentroidX())) { - // use geometry y overlap scheme to see if this is really a match in x - // should be done in a map - - if (centroid >= vert_bar_start_idx_ && - seed.getCentroidY() < cluster1.getCentroidY()) { - // impossible combination - if (verbose_ > 1) { - ldmx_log(debug) << "\tSkipping impossible x cluster combination " - "with y flags (tag up) (" - << seed.getCentroidY() << " " - << cluster1.getCentroidY() << ")"; + + if (lut_tracking_) { //if using LUT method + for (const auto &cluster1 : clusters_pad1) { + for (const auto &cluster2 : clusters_pad2) { + + float seed_bin = seed.getCentroid(); + float pad1_bin = cluster1.getCentroid(); + float pad2_bin = cluster2.getCentroid(); + + LUTKey key{seed_bin, pad1_bin, pad2_bin}; //LUTKey defined in header file + + if (lut_.find(key) != lut_.end()) { + + std::vector three_cluster_vec = { + seed, cluster1, cluster2}; + + ldmx::TrigScintTrack track = makeTrack(three_cluster_vec); + track_candidates.push_back(track); } - continue; } - - // else: first (possible) match! loop through next pad too - + } + + } + + else { //ends line 283, normal tstp tracking + + for (const auto &cluster1 : clusters_pad1) { if (verbose_ > 1) { - ldmx_log(debug) << "\t\tIt is close enough!. Check pad2"; + ldmx_log(debug) << "\tGot pad1 cluster with centroid " + << cluster1.getCentroid(); } + if (fabs(cluster1.getCentroid() - centroid) < max_delta_ || + (centroid >= vert_bar_start_idx_ // then in vertical bars + && seed.getCentroidX() == cluster1.getCentroidX())) { + // use geometry y overlap scheme to see if this is really a match in x + // should be done in a map + + if (centroid >= vert_bar_start_idx_ && + seed.getCentroidY() < cluster1.getCentroidY()) { + // impossible combination + if (verbose_ > 1) { + ldmx_log(debug) << "\tSkipping impossible x cluster combination " + "with y flags (tag up) (" + << seed.getCentroidY() << " " + << cluster1.getCentroidY() << ")"; + } + continue; + } - // try making third pad clusters an optional part of track - - std::vector cluster_vec = {seed, cluster1}; - - bool has_match_dn = false; + // else: first (possible) match! loop through next pad too - for (const auto &cluster2 : clusters_pad2) { if (verbose_ > 1) { - ldmx_log(debug) << "\tGot pad2 cluster with centroid " - << cluster2.getCentroid(); + ldmx_log(debug) << "\t\tIt is close enough!. Check pad2"; } - if (fabs(cluster2.getCentroid() - centroid) < max_delta_ || - (centroid >= vert_bar_start_idx_ // then in vertical bars - && seed.getCentroidX() == cluster2.getCentroidX())) { - // use geometry y overlap scheme to see if this is really a match - // in x + // try making third pad clusters an optional part of track - if (centroid >= vert_bar_start_idx_ && - (seed.getCentroidY() < cluster2.getCentroidY() || - cluster1.getCentroidY() > - cluster2.getCentroidY())) { // impossible - if (verbose_ > 1) { - ldmx_log(debug) - << "\tSkipping impossible x cluster combination with y " - "flags (tag up dn) (" - << seed.getCentroidY() << " " << cluster1.getCentroidY() - << " " << cluster2.getCentroidY() << ")"; - } - continue; - } + std::vector cluster_vec = {seed, cluster1}; - // first match! loop through next pad too + bool has_match_dn = false; + for (const auto &cluster2 : clusters_pad2) { if (verbose_ > 1) { - ldmx_log(debug) << "\t\tIt is close enough!. Make a track"; + ldmx_log(debug) << "\tGot pad2 cluster with centroid " + << cluster2.getCentroid(); } - // only make this vector now! this ensures against hanging - // clusters with indices from earlier in the loop - std::vector three_cluster_vec = { - seed, cluster1, cluster2}; - - /* - // here we could break if we didn't want to allow all possible - combinations madeTrack=true; break; //we're done with this - iteration once there's a track made - */ - // make a track - ldmx::TrigScintTrack track = makeTrack(three_cluster_vec); + if (fabs(cluster2.getCentroid() - centroid) < max_delta_ || + (centroid >= vert_bar_start_idx_ // then in vertical bars + && seed.getCentroidX() == cluster2.getCentroidX())) { + // use geometry y overlap scheme to see if this is really a match + // in x + + if (centroid >= vert_bar_start_idx_ && + (seed.getCentroidY() < cluster2.getCentroidY() || + cluster1.getCentroidY() > + cluster2.getCentroidY())) { // impossible + if (verbose_ > 1) { + ldmx_log(debug) + << "\tSkipping impossible x cluster combination with y " + "flags (tag up dn) (" + << seed.getCentroidY() << " " << cluster1.getCentroidY() + << " " << cluster2.getCentroidY() << ")"; + } + continue; + } + + // first match! loop through next pad too + + if (verbose_ > 1) { + ldmx_log(debug) << "\t\tIt is close enough!. Make a track"; + } + + // only make this vector now! this ensures against hanging + // clusters with indices from earlier in the loop + std::vector three_cluster_vec = { + seed, cluster1, cluster2}; + + /* + // here we could break if we didn't want to allow all possible + combinations madeTrack=true; break; //we're done with this + iteration once there's a track made + */ + // make a track + ldmx::TrigScintTrack track = makeTrack(three_cluster_vec); + track_candidates.push_back(track); + has_match_dn = true; + } // if match in pad2 + } // over clusters in pad2 + // if there was no match to this in pad 2, make a track with just + // these two clusters + if (!has_match_dn && skip_last_) { + // we allow skipping last pad if needed + ldmx::TrigScintTrack track = makeTrack(cluster_vec); track_candidates.push_back(track); - has_match_dn = true; - } // if match in pad2 - } // over clusters in pad2 - // if there was no match to this in pad 2, make a track with just - // these two clusters - if (!has_match_dn && skip_last_) { - // we allow skipping last pad if needed - ldmx::TrigScintTrack track = makeTrack(cluster_vec); - track_candidates.push_back(track); - } + } - } // if possible (x,)y match in pad1 - /* + } // if possible (x,)y match in pad1 + /* //same here if (madeTrack) break; */ - } // over clusters in pad1 - + } // over clusters in pad1 + + } + // continue to next seed if 0 track candidates if (track_candidates.size() == 0) continue; From e5b95d0a4b3d54b4b2b56c4328e08e03a2e58537 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Mon, 27 Apr 2026 10:36:53 +0200 Subject: [PATCH 06/20] Remove previous LUTmethod attempt, update configs --- ClusterViewerAnalyzer.cxx | 14 +- ClusterViewerFig.py | 9 - LUTAnalyzer.cxx | 49 +- LUTconfig.py | 10 +- .../include/TrigScint/TrigScintLUTTracker.h | 125 --- .../src/TrigScint/TrigScintLUTTracker.cxx | 844 ------------------ extruthfig.py | 8 +- notruthfig.py | 11 +- 8 files changed, 59 insertions(+), 1011 deletions(-) delete mode 100644 ClusterViewerFig.py delete mode 100644 TrigScint/include/TrigScint/TrigScintLUTTracker.h delete mode 100644 TrigScint/src/TrigScint/TrigScintLUTTracker.cxx diff --git a/ClusterViewerAnalyzer.cxx b/ClusterViewerAnalyzer.cxx index 7e6d344bd4..2ece585615 100644 --- a/ClusterViewerAnalyzer.cxx +++ b/ClusterViewerAnalyzer.cxx @@ -11,7 +11,13 @@ class ClusterViewerAnalyzer : public framework::Analyzer { public: std::ofstream outfile; - + + void configure(framework::config::Parameters& ps) override { + pass_name_ = ps.get("pass_name"); + } + + std::string pass_name_; + ClusterViewerAnalyzer(const std::string& name, framework::Process& p) : framework::Analyzer(name, p) { @@ -21,13 +27,13 @@ class ClusterViewerAnalyzer : public framework::Analyzer { void analyze(const framework::Event& event) override { const auto& clusters1 = - event.getCollection("TriggerPad1Clusters","truthisolater"); + event.getCollection("TriggerPad1Clusters", pass_name_); const auto& clusters2 = - event.getCollection("TriggerPad2Clusters","truthisolater"); + event.getCollection("TriggerPad2Clusters", pass_name_); const auto& clusters3 = - event.getCollection("TriggerPad3Clusters","truthisolater"); + event.getCollection("TriggerPad3Clusters", pass_name_); size_t Clusters = std::min({clusters1.size(), clusters2.size(), clusters3.size()}); diff --git a/ClusterViewerFig.py b/ClusterViewerFig.py deleted file mode 100644 index e7f2861c2a..0000000000 --- a/ClusterViewerFig.py +++ /dev/null @@ -1,9 +0,0 @@ -#CONFIG FILE TO VIEW CLUSTERS - -from LDMX.Framework import ldmxcfg -p = ldmxcfg.Process('clusterviewer') -#p.max_events = 10 -p.input_files = ["TruthEvents1000.root"] #change for new naming conventions -#p.ampl_weighting = False -p.sequence = [ldmxcfg.Analyzer.from_file('ClusterViewerAnalyzer.cxx', - needs = ['TrigScint_Event', 'SimCore_Event'])] diff --git a/LUTAnalyzer.cxx b/LUTAnalyzer.cxx index 63040970ab..07c1bed750 100644 --- a/LUTAnalyzer.cxx +++ b/LUTAnalyzer.cxx @@ -1,4 +1,4 @@ - +//WRITES LUT BASED ON FREQUENCY OF TRACK CANDIDATE PROPAGATION COMBINATIONS #include "Framework/EventProcessor.h" #include "Framework/Event.h" #include "TrigScint/Event/TrigScintTrack.h" @@ -13,20 +13,39 @@ class LUTAnalyzer : public framework::Analyzer { struct Line {int event; float p1, p2, p3, p3alt; }; - std::ifstream infile; - std::ofstream outfile; + void configure(framework::config::Parameters& ps) override { + input_file_ = ps.get("input_file"); + output_file_ = ps.get("output_file"); + lut_threshold_ = ps.get("lut_threshold"); + //The LUT threshold is the minimum percentage of times that a track + //must appear in a pool of events to be written to the LUT. + //Example --- straight tracks whose pad 1-2 and pad 2-3 delta values + // are both 0 make up ~85% of tracks out of 10.000 events, + // while a track with deltas (+22,-22), meaning pad 1 cluster in + // e.g. bar 4, pad 2 cluster in bar 26, and pad 3 cluster in bar 4 + // only appears once in 10.000 events (0.01%). The straight + // tracks are written to the LUT and the single "anomaly" is not + + } + std::ifstream infile; + std::ofstream outfile; + std::string input_file_; + std::string output_file_; + double lut_threshold_; + std::map, std::vector> groups; int totalLines = 0; LUTAnalyzer(const std::string& name, framework::Process& p) : framework::Analyzer(name,p) { - - infile.open("clusters.txt"); - outfile.open("LUT.txt"); } - + void onProcessStart() override { + infile.open(input_file_); + outfile.open(output_file_); + } + void analyze(const framework::Event& event) override { if (totalLines > 0) return; @@ -50,9 +69,7 @@ class LUTAnalyzer : public framework::Analyzer { std::cout << "Total number of track candidates: " << totalLines << "\n"; std::cout << "Number of track candidate types: " << groups.size() << "\n"; - - double threshold = 1.0/1000.0; - std::cout << "Threshold: " << threshold << "\n"; + std::cout << "LUT Threshold: " << lut_threshold_ * 100 << "%\n"; for (auto& g : groups) { @@ -63,7 +80,7 @@ class LUTAnalyzer : public framework::Analyzer { double frac = (double)count / totalLines; std::cout << "(" << p12 << "," << p23 << ") appears " - << count << " times, (" << frac << " %)" << "\n"; + << count << " times, (" << frac * 100 << " %)" << "\n"; } for (auto& g : groups) { @@ -71,7 +88,7 @@ class LUTAnalyzer : public framework::Analyzer { int count = g.second.size(); double frac = (double)count / totalLines; - if (frac > threshold) { + if (frac > lut_threshold_) { combs++; @@ -91,11 +108,3 @@ class LUTAnalyzer : public framework::Analyzer { }; DECLARE_ANALYZER(LUTAnalyzer) - - - - - - - - diff --git a/LUTconfig.py b/LUTconfig.py index 166275734f..fe2dfd1e3d 100644 --- a/LUTconfig.py +++ b/LUTconfig.py @@ -4,7 +4,11 @@ p = ldmxcfg.Process("LUTmaker") p.input_files = ["TruthSamples.root"] #necessary to define but it doesnt use it -#p.output_files = ["Tracks.root"] -p.sequence = [ldmxcfg.Analyzer.from_file('LUTAnalyzer.cxx', needs = ['TrigScint_Event', 'SimCore_Event']) - ] +make_lut = ldmxcfg.Analyzer.from_file('LUTAnalyzer.cxx', + needs = ['TrigScint_Event', 'SimCore_Event']) +make_lut.lut_threshold = 1.0/1000 +make_lut.input_file = "clusters.txt" +make_lut.output_file = "LUT.txt" + +p.sequence = [make_lut] diff --git a/TrigScint/include/TrigScint/TrigScintLUTTracker.h b/TrigScint/include/TrigScint/TrigScintLUTTracker.h deleted file mode 100644 index 2102fd5f83..0000000000 --- a/TrigScint/include/TrigScint/TrigScintLUTTracker.h +++ /dev/null @@ -1,125 +0,0 @@ - -#ifndef TRIGSCINT_TRIGSCINTTRACKPRODUCER_H -#define TRIGSCINT_TRIGSCINTTRACKPRODUCER_H - -// LDMX Framework -#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 "Recon/Event/EventConstants.h" -#include "TrigScint/Event/TrigScintCluster.h" -#include "TrigScint/Event/TrigScintTrack.h" - -#include - -namespace trigscint { - -/** - * @class TrigScintTrackProducer - * @brief making tracks from trigger scintillator clusters - */ -class TrigScintLUTTracker : public framework::Producer { - public: - TrigScintLUTTracker(const std::string &name, framework::Process &process) - : Producer(name, process) {} - - void configure(framework::config::Parameters &ps) override; - - void produce(framework::Event &event) override; - - void onProcessStart() override; - void onProcessEnd() override; - - private: - // collection of produced tracks - std::vector tracks_; - - // add a cluster to a track - ldmx::TrigScintTrack makeTrack(std::vector clusters); - - // match x, y tracks and set their x,y spatial coordinates - void matchXYTracks(std::vector &tracks); - // std::vector matchXYTracks( - // std::vector &tracks); - - // maximum difference (in channel number space) between track seed and cluster - // in the next pad tolerated to form a track - double max_delta_{0.}; - - // producer specific verbosity - int verbose_{0}; - - // collection used to seed the tracks - std::string seeding_collection_; - - // other cluster collections used in track making - std::vector input_collections_; - - // output collection (tracks) - std::string output_collection_; - - // specific pass name to use for track making - std::string pass_name_{""}; - - // allow forming tracks without match in the last collection - bool skip_last_{false}; - - // vertical bar start index - int vert_bar_start_idx_{52}; - - // number of horizontal bars (one layer) in geometry - int n_bars_y_{16}; - - // number of vertical bars (one row) in geometry - int n_bars_x_{8}; - - // track centroid in units of channel nb (will not be content weighted) - // float centroid_{0.}; - - // track horizontal centroid in units of channel nb (will not be content - // weighted) - // float centroidX_{-1}; - - // track vertical centroid in units of channel nb (will not be content - // weighted) - // float centroidY_{-1}; - - // track residual in units of channel nb (will not be content weighted) - // float residual_{0.}; - - // LUT part starts here - - struct LUTKey { - float p1, p2, p3; - - bool operator==(const LUTKey& other) const { - return p1 == other.p1 && p2 == other.p2 && p3 == other.p3; - } - }; - - struct LUTKeyHash { - size_t operator()(const LUTKey& k) const { - return std::hash()(k.p1) ^ - (std::hash()(k.p2) << 1) ^ - (std::hash()(k.p3) << 2); - } - }; - - std::unordered_set lut_; - - // LUT ends here - - float bar_width_y_{3.}; // mm - float bar_gap_y_{2.1}; // mm - float bar_width_x_{3.}; // mm - float bar_gap_x_{0.1}; // mm - - float x_conv_factor_; // geometry conversion factors - float x_start_; - float y_conv_factor_; - float y_start_; - }; - -} // namespace trigscint - -#endif // TRIGSCINT_TRIGSCINTLUTTRACKER_H diff --git a/TrigScint/src/TrigScint/TrigScintLUTTracker.cxx b/TrigScint/src/TrigScint/TrigScintLUTTracker.cxx deleted file mode 100644 index 23ec09f333..0000000000 --- a/TrigScint/src/TrigScint/TrigScintLUTTracker.cxx +++ /dev/null @@ -1,844 +0,0 @@ -#include "TrigScint/TrigScintLUTTracker.h" - -#include // std::next -#include -#include -#include - -namespace trigscint { - -void TrigScintLUTTracker::configure(framework::config::Parameters &ps) { - max_delta_ = ps.get( - "delta_max"); // max distance to consider adding in a cluster to track - seeding_collection_ = ps.get( - "seeding_collection"); // probably tagger pad, "TriggerPadTagClusters" - input_collections_ = ps.get>( - "further_input_collections"); // {"TriggerPadUpClusters" , - // "TriggerPadDnClusters" } - output_collection_ = ps.get("output_collection"); - pass_name_ = ps.get("input_pass_name"); - verbose_ = ps.get("verbosity"); - vert_bar_start_idx_ = ps.get("vertical_bar_start_index"); - n_bars_y_ = ps.get("number_horizontal_bars"); - bar_width_y_ = ps.get("horizontal_bar_width"); - bar_gap_y_ = ps.get("horizontal_bar_gap"); - n_bars_x_ = ps.get("number_vertical_bars"); - bar_width_x_ = ps.get("vertical_bar_width"); - bar_gap_x_ = ps.get("vertical_bar_gap"); - skip_last_ = ps.get("allow_skip_last_collection"); - std::string lut_file_ = ps.get("lut_file"); //LUT ! new! - - // TO DO: allow any number of input collections - - if (verbose_) { - ldmx_log(info) << "In TrigScintLUTTracker: configure done!" << std::endl; - ldmx_log(info) << "Got parameters: \nSeeding: " << seeding_collection_ - << "\nTolerance: " << max_delta_ - << "\nInput: " << input_collections_.at(0) << " and " - << input_collections_.at(1) - << "\nInput pass name: " << pass_name_ - << "\nAllow tracks with no hit in last collection: " - << skip_last_ - << "\nVertical bar start index: " << vert_bar_start_idx_ - << "\nNumber of horizontal bars: " << n_bars_y_ - << "\nHorizontal bar width: " << bar_width_y_ - << "\nHorizontal bar gap: " << bar_gap_y_ - << "\nNumber of vertical bars: " << n_bars_x_ - << "\nVertical bar width: " << bar_width_x_ - << "\nVertical bar gap: " << bar_gap_x_ - << "\nOutput: " << output_collection_ - << "\nVerbosity: " << verbose_ - << "\nTracks (LUT): " << lut_file_; //is this necessary too? - } - // each bar only goes half this distance up (overlap/zig-zag) - y_conv_factor_ = (bar_width_y_ + bar_gap_y_) / 2.; - // half height of pad - y_start_ = -(n_bars_y_ * (bar_width_y_ + bar_gap_y_) - bar_gap_y_) / 2.; - // each bar goes entire distance sideways (no overlap) - x_conv_factor_ = bar_width_x_ + bar_gap_x_; - // half width of pad - x_start_ = -(n_bars_x_ * (bar_width_x_ + bar_gap_x_) - bar_gap_x_) / 2.; - - //also LUT related here - std::ifstream file(lut_file_); //open file - float a, b, c; - - while (file >> a >> b >> c) { - float p1 = a; - float p2 = b; - float p3 = c; - lut_.insert({p1, p2, p3}); - } - - ldmx_log(info) << "Loaded LUT with size: " << lut_.size(); - - return; -} //end of configure - - - -void TrigScintLUTTracker::produce(framework::Event &event) { - // parameters. - // one pad cluster collection to use as seed - // a vector with the other two - // a maximum distance between seed centroid and other pad clusters - // allowed to belong to the same track - // an output collection name a verbosity controller - - if (verbose_) { - ldmx_log(debug) - << "TrigScintLUTTracker: produce() starts! Event number: " - << event.getEventHeader().getEventNumber(); - } - - if (!event.exists(seeding_collection_, pass_name_)) { - ldmx_log(info) << "No collection called " << seeding_collection_ - << "; skipping event"; - // << "; still, not skipping event"; - return; - } - - if (!event.exists(seeding_collection_, pass_name_)) { - ldmx_log(info) << "No collection called " << seeding_collection_ - << "; skipping event"; - return; - } - - const auto seeds{event.getCollection( - seeding_collection_, pass_name_)}; - uint num_seeds = seeds.size(); - - if (verbose_) { - ldmx_log(debug) << "Got track seeding cluster collection " - << seeding_collection_ << " with " << num_seeds - << " entries "; - } - - if (!event.exists(input_collections_.at(0), pass_name_)) { - ldmx_log(info) << "No collection called " << input_collections_.at(0) - << "; skipping event"; - // << "; still, not skipping event"; - - return; - } - - const auto clusters_pad1{event.getCollection( - input_collections_.at(0), pass_name_)}; - - if (!event.exists(input_collections_.at(1), pass_name_)) { - ldmx_log(info) << "No collection called " - << input_collections_.at(1) - // << "; still, not skipping event"; - << "; skipping event"; - std::vector empty{}; - event.add(output_collection_, empty); - return; - } - - const auto clusters_pad2{event.getCollection( - input_collections_.at(1), pass_name_)}; - - if (verbose_) { - ldmx_log(debug) << "Got the other two pad collections:" - << input_collections_.at(0) << " with " - << clusters_pad1.size() << " entries, and " - << input_collections_.at(1) << " with " - << clusters_pad2.size() << " entries."; - } - - std::vector cleaned_tracks; - std::vector cleaned_tracks_y; - std::vector cleaned_tracks_x; - - - - -//----------------------------------------------------------------------------- - - //this is what 130-250 have been replaced with: - - // loop over the clusters in the seeding pad collection, if there are clusters - // in all pads - // bool skipDn = false; - if (num_seeds && clusters_pad1.size()) { //ends on 407, if there are clusters in all pads - // could check this explicitly here: and then just get out of all checks on - // the dn pad immediately - // if (! clusters_pad2.size()) - // skipDn = true ; - for (const auto &seed : seeds) { //ends on line 247, over all seeds - // for each seed, search through the other two pads to match all clusters - // with centroids within tolerance to tracks - float centroid = seed.getCentroid(); - - std::vector track_candidates; - - if (verbose_ > 1) { - ldmx_log(debug) << "Got seed with centroid " << centroid; - } - - // reset for each seed - // bool madeTrack = false; - - - //LUT part here! - for (const auto &cluster1 : clusters_pad1) { - for (const auto &cluster2 : clusters_pad2) { - - float seed_bin = seed.getCentroid(); - float pad1_bin = cluster1.getCentroid(); - float pad2_bin = cluster2.getCentroid(); - - LUTKey key{seed_bin, pad1_bin, pad2_bin}; //LUTKey defined in header file - - if (lut_.find(key) != lut_.end()) { - - std::vector three_cluster_vec = { - seed, cluster1, cluster2}; - - ldmx::TrigScintTrack track = makeTrack(three_cluster_vec); - track_candidates.push_back(track); - } - } - } // over clusters in pad1, //LUT ends here - -//----------------------------------------------------------------------------- - - - // continue to next seed if 0 track candidates - if (track_candidates.size() == 0) continue; - - int keep_idx = 0; - float min_residual = 1000; // some large number - - // no need to choose between only one candidate track - if (track_candidates.size() > 1) { - // now for each seed, pick only the track with the smallest residual. - - if (verbose_) { - ldmx_log(debug) << "Got " << track_candidates.size() - << " tracks to check."; - } - - for (uint idx = 0; idx < track_candidates.size(); idx++) { - if ((track_candidates.at(idx)).getResidual() < min_residual) { - keep_idx = (int)idx; - min_residual = - (track_candidates.at(idx)).getResidual(); // update minimum - - if (verbose_ > 1) { - ldmx_log(debug) - << "Track at index " << idx - << " has smallest residual so far: " << min_residual; - } - - } // finding min residual - } // over track candidates - } // if more than 1 to choose from - - // store the track at keepIdx, if there was one we made it this far and - // keepIdx is 0 or has been updated to the smallest residual track idx - // if (keepIdx >= 0) { - tracks_.push_back(track_candidates.at(keep_idx)); - if (verbose_) { - ldmx_log(debug) << "Kept track at index " << keep_idx; - ldmx_log(trace) << track_candidates.at(keep_idx); - } - //} - } // over seeds - - // done here if there were no tracks found - if (tracks_.size() == 0) { - if (verbose_) { - ldmx_log(debug) << "No tracks found!"; - } - std::vector empty{}; - event.add(output_collection_, empty); - return; - } - // now, if there are multiple seeds sharing the same downstream hits, this - // should also be remedied with a selection on min residual. - - // The logic of this loop kind of assumes I can remove tracks immediately -- - // that way I can do pairwise checks between more tracks within a single - // loop. But for now I haven't figured out how to erase elements in a fool - // proof way. So I iterate over a vector... - - std::vector keep_indices(tracks_.size(), 1); - if (verbose_ > 1) - ldmx_log(debug) << "vector of indices to keep has size " - << keep_indices.size(); - - for (uint idx = tracks_.size() - 1; idx > 0; idx--) { - // since we start in one end, we only have to check matches in one - // direction - ldmx::TrigScintTrack track = tracks_.at(idx); - for (int idx_comp = idx - 1; idx_comp >= 0; idx_comp--) { - if (verbose_ > 1) - ldmx_log(debug) << "In track disambiguation loop, idx points at " - << idx << " and prev idx points at " << idx_comp; - - ldmx::TrigScintTrack next_track = tracks_.at(idx_comp); - - // no need to start pulling constituents from tracks that are - // ridiculously far apart - if (fabs(track.getCentroid() - next_track.getCentroid()) < - 3 * max_delta_) { - std::vector consts_1 = - track.getConstituents(); - std::vector consts_2 = - next_track.getConstituents(); - if (verbose_ > 1) - ldmx_log(debug) - << "In track disambiguation loop, got the two tracks, " - "with nConstituents " - << consts_1.size() << " and " << consts_2.size() - << ", respectively. "; - // let's do "if either cluster is shared" right now... but could also - // have it settable to use a stricter cut: an AND - if (consts_1[1].getCentroid() == consts_2[1].getCentroid() || - (consts_1.size() > 2 && consts_2.size() > 2 && - consts_1[2].getCentroid() == - consts_2[2] - .getCentroid())) { // we have overlap downstream of the - // seeding pad. probably, one cluster - // in seeding pad is noise - - if (verbose_ > 1) { - ldmx_log(debug) << "Found overlap! Tracks at index " << idx - << " and " << idx_comp; - ldmx_log(trace) << tracks_.at(idx); - ldmx_log(trace) << tracks_.at(idx_comp); - } - - if ((tracks_.at(idx)).getResidual() < - (tracks_.at(idx_comp)).getResidual()) { - // next track (lower index) is a worse choice, remove its flag for - // keeping - keep_indices.at(idx_comp) = 0; - } else // prefer next track over current. remove current track's - // keep - // flag - keep_indices.at(idx) = 0; - /*} - else { - tracks_.erase(itNext); - // removeIdx.push_back(idx+1); - // we might see the same index two times in the loop in this - case, if there are three seeds sharing the same clusters - downstream. - // then the third only gets removed if it's even worse than - the second. - // one could deal with this with an extra overlap check. not - sure we will be in this situation any time soon though. - }*/ - } // over matching/overlapping tracks - } // over tracks close enough to share constituents - } // over constructed tracks at other indices, to match - } // over constructed tracks - - for (uint idx = 0; idx < tracks_.size(); idx++) { - if (verbose_ > 1) { - ldmx_log(debug) << "keep flag for idx " << idx << " is " - << keep_indices.at(idx); - } - if (keep_indices.at(idx)) { // this hasn't been flagged for removal - - cleaned_tracks.push_back(tracks_.at(idx)); - - if (verbose_) { - ldmx_log(debug) << "After cleaning, keeping track at index " << idx - << ": Centroid = " << (tracks_.at(idx)).getCentroid() - << "; CentroidX = " - << (tracks_.at(idx)).getCentroidX() - << "; CentroidY = " - << (tracks_.at(idx)).getCentroidY() - << "; track PE = " << (tracks_.at(idx)).getPE() - << tracks_.at(idx); - } - } // if index flagged for keeping - } // over all (uniquely seeded) tracks in the event - - if (verbose_) { - for (uint idx = 0; idx < tracks_.size(); idx++) { - ldmx_log(debug) << "Keeping track at index " << idx << ":" - << tracks_.at(idx); - } - } - - if (verbose_) { - ldmx_log(debug) << "Running track x,y matching "; - } - - if (cleaned_tracks.size() > 0) { - matchXYTracks(cleaned_tracks); - std::vector matched_tracks = - cleaned_tracks; // don't know why this copying needs to happen but it - // does - // std::vector matchXYTracks( cleanedTracks - //); std::vector matchedTracks = - // matchXYTracks( cleanedTracks ); - for (auto trk : matched_tracks) { - /* for (uint idx = 0; idx < tracks_.size(); idx++) { - if (verbose_ > 1) { - ldmx_log(debug) << "keep flag for idx " << idx << " is " - << keepIndices.at(idx); - } - if (keepIndices.at(idx)) { // this hasn't - been flagged for removal - //check if channel nb is above that of horizontal bars - if (tracks_.at(idx).getCentroid() >= vert_bar_start_idx_) - */ - if (trk.getCentroid() >= vert_bar_start_idx_) - cleaned_tracks_x.push_back(trk); // acks_.at(idx)); - else - cleaned_tracks_y.push_back(trk); // acks_.at(idx)); - // cleanedTracksY.push_back(trk); - if (verbose_ > 1) { - float centr = trk.getCentroid(); // tracks_.at(idx).getCentroid(); // - std::string coll_str = centr >= vert_bar_start_idx_ ? "X" : "Y"; - coll_str = output_collection_ + coll_str; - ldmx_log(debug) << "saving track with centroid " << centr - << " to output track collection " << coll_str; - } - // } - } - } - - } // if there are clusters in all pads - - - - else if (verbose_) { //if theyre arent clusters in all pads - ldmx_log(info) - << "Not all pads had clusters; (maybe) skipping tracking attempt"; - } - - if (verbose_) { - ldmx_log(debug) << "Done with tracking step. "; - } - - event.add(output_collection_, cleaned_tracks); - // event.add(output_collection_, matchedTracks); - - event.add(output_collection_ + "Y", cleaned_tracks_y); - event.add(output_collection_ + "X", cleaned_tracks_x); - - tracks_.resize(0); - - return; -} //produce ends - - -//---------------------------------------------------------------------------- - - - -ldmx::TrigScintTrack TrigScintLUTTracker::makeTrack( - std::vector clusters) { - // for now let's keep a straight, unweighted centroid - // consider the possibility that at least one cluster has a centroid - // identically == 0. then we need to shift them by 1 if we want to do energy - // weighted track centroid later. but no need now - ldmx::TrigScintTrack tr; - float centroid = 0; - float centroid_x = 0; - float centroid_y = 0; - float beam_efrac = 0; - float pe = 0; - for (uint i = 0; i < clusters.size(); i++) { - centroid += (clusters.at(i)).getCentroid(); - centroid_x += (clusters.at(i)).getCentroidX(); - centroid_y += (clusters.at(i)).getCentroidY(); - tr.addConstituent(clusters.at(i)); - beam_efrac += (clusters.at(i)).getBeamEfrac(); - pe += (clusters.at(i)).getPE(); - } - centroid /= clusters.size(); - centroid_x /= clusters.size(); - if (centroid >= vert_bar_start_idx_) { - if (verbose_) { - ldmx_log(debug) - << " -- In makeTrack made vertical bar track with centroid " - << centroid << " and y flag sum " << centroid_y; - // try commenting this to check if that helps with an out-of-bounds - // problem - // << " from clusters with y centroids"; - // for (uint i = 0; i < clusters.size(); i++) - // ldmx_log(debug) << "\tpad " << i << ": centroidY " - // << (clusters.at(i)).getCentroidY(); - } - // then the sum of centroid y is 0, 2, 4 or 6 - // we have 4 divisions, so, the center of it should be divNb/8 - // (or rather, that's where channel nBars/8 begins) - // and then a factor 2 for the zig-zag pattern - centroid_y = (centroid_y + 1) * 2 * n_bars_y_ / 8.; - // TODO: here we could instead just use quadrant indices 0-3 by dividing by - // 2 but that would mean that in the raw, x and y track centroidY would mean - // different things - if (verbose_) ldmx_log(debug) << " -- new centroidY = " << centroid_y; - } else - centroid_y /= clusters.size(); - - beam_efrac /= clusters.size(); - pe /= clusters.size(); - - float residual = 0; - for (uint i = 0; i < clusters.size(); i++) - residual += ((clusters.at(i)).getCentroid() - centroid) * - ((clusters.at(i)).getCentroid() - centroid); - residual = sqrt(residual / clusters.size()); - - tr.setCentroid(centroid); - tr.setCentroidX(centroid_x); - tr.setCentroidY(centroid_y); - tr.setResidual(residual); - tr.setBeamEfrac(beam_efrac); - tr.setPE(pe); - - if (verbose_) { - ldmx_log(debug) << " -- In makeTrack made track with centroid " - << centroid << " and residual " << residual << " and pe " - << pe << " from clusters with centroids"; - for (uint i = 0; i < clusters.size(); i++) - ldmx_log(debug) << "\tpad " << i << ": centroid " - << (clusters.at(i)).getCentroid(); - } - - return tr; -} - - -// std::vector TrigScintTrackProducer::matchXYTracks( -void TrigScintLUTTracker::matchXYTracks( - std::vector &tracks) { - // map quadrant nb to track (can be multiple per quadrant) - std::multimap - y_idx_quad_map; // key = quad, val = track index in collection - std::multimap x_idx_quad_map; - - std::multimap y_quad_map; - std::multimap x_quad_map; - // map track in quadrant back to index in entire track collection - // used for updating collection track variables - std::map y_track_map; - std::map x_track_map; - - uint trk_idx = -1; - for (auto trk : tracks) { - trk_idx++; - // 1. get the y bar tracks with centroidX = -1 - if (trk.getCentroidX() == -1) { - if (verbose_) - ldmx_log(debug) << " -- In matchXYTracks found y track at " - << trk.getCentroidY() << "; mapping to quad " - << (int)trk.getCentroidY() / 8 << " with trk index " - << trk_idx; - // 2. order them... or map them to quadrants. note that there are 2 layers - // so 2*n_bars_y_/4 channels per quadrant - y_quad_map.insert(std::make_pair((int)(trk.getCentroidY() / 8), trk)); - y_track_map[trk] = trk_idx; - y_idx_quad_map.insert( - std::make_pair((int)(trk.getCentroidY() / 8), trk_idx)); - - } else { // 3. get the remaining tracks (from vertical bars) and map them - // (back) to (middle of) quadrants - x_quad_map.insert(std::make_pair((int)(trk.getCentroidY() / 8), trk)); - x_track_map[trk] = trk_idx; - x_idx_quad_map.insert( - std::make_pair((int)(trk.getCentroidY() / 8), trk_idx)); - if (verbose_) - ldmx_log(debug) << " -- In matchXYTracks found x track at (x,y) = (" - << trk.getCentroidX() << ", " << trk.getCentroidY() - << "); mapping to quad " << (int)trk.getCentroidY() / 8 - << " with trk index " << trk_idx; - } - } - - // 4a - // - // 1) here use the geometry? if we can assume perfect alignment we can take - // width and nBars and take nBars/2 as origin - // --- now do the matching --- - - // if there is no useful matching to be done: these are the pad width wide - // numbers - float x0 = 0; - // this should be half the pad... could also set - // it to full beam spot width - float sx0 = fabs(x_start_); - - // y_start_ is half the pad, so this should be half a quadrant - float sy0 = fabs(y_start_) / 4.; - - // assume at least one y track. will have to figure out if there is ever a - // reason to use an isolated x track in its place. - for (auto yitr = y_quad_map.begin(); yitr != y_quad_map.end(); ++yitr) { - int n_yin_quad = y_quad_map.count((*yitr).first); - int n_xin_quad = x_quad_map.count((*yitr).first); - float y{-9999.}, sy{-9999.}, x{-9999.}, x1{-9999.}, x2{-9999.}, sx1{-9999.}, - sx2{-9999.}, y1{-9999.}, y2{-9999.}, sy1{-9999.}, sy2{-9999.}; - // quad midpoint: - float y0 = (*yitr).first * 8 + sy0; - float sx = 1. / 2 * - x_conv_factor_; // rely on x precision being one single bar - // width; always used unless x is undeterminable - - // check all x first - // do the easiest first: - if (n_xin_quad == 0) { // then there's no hope of setting a better x here - // just use the beam spot width... and center of pad - x = x0; - sx = sx0; - if (verbose_) - ldmx_log(debug) << "\t\t\t no x info in quad " << (*yitr).first - << "; will set x to middle of pad, pad half-width as " - "precision: set (x, sx)=(" - << x << ", " << sx << ")"; - } // 0 x tracks in quadrant - else if (n_xin_quad == - 1) { // slightly harder: 1 x track -- might be easy if - // it's just one y track; if several, need to - // think about overlaps. but in overlap case, just - // revert to setting x0 and sx0, when we know - auto xitr = x_quad_map.find((*yitr).first); - x = ((*xitr).second).getCentroidX() * x_conv_factor_ + x_start_; - - if (verbose_) - ldmx_log(debug) << "\t\t\t 1 x in quad " << (*yitr).first - << ", getting (x, sx)=(" << x << ", " << sx << ")"; - } // 1 x track in quadrant - else if (n_xin_quad == 2) { // finally if we have two tracks, get x1 and x2 - // and decide later how to use them - // don't think we want to experiment with discerning three overlapping - // tracks, so not >= 2 - // continue; //debugging: skip for now -- didn't help - auto xitr1 = x_quad_map.lower_bound((*yitr).first); - auto xitr2 = x_quad_map.upper_bound((*yitr).first); - xitr2--; // upper_bound points to next element - - if (xitr1 != xitr2) { // should be true already but... - x1 = ((*xitr1).second).getCentroidX() * x_conv_factor_ + x_start_; - x2 = ((*xitr2).second).getCentroidX() * x_conv_factor_ + x_start_; - sx1 = x_conv_factor_ / 2.; // 1 bar width - sx2 = sx1; - x = (x1 + x2) / 2.; - sx = fabs(x1 - x2) / 2 * - x_conv_factor_; // rely on x precision being one single pad width - if (verbose_) - ldmx_log(debug) << "\t\t -- 2 x in quad: setting y track x " - "coordinate to midpoint"; - } - } // if 2 x tracks in quad - - // ok! over y: - // can skip 0 y case by construction - if (n_yin_quad == 1) { // we can already now tell what the y coordinate and - // its precision is - y = ((*yitr).second).getCentroidY() * y_conv_factor_ + y_start_; - sy = ((*yitr).second).getResidual() * y_conv_factor_; - // if all clusters lined up, assign - // precision of 1 bar width - if (sy == 0) sy = 1. / 2 * y_conv_factor_; - - if (n_xin_quad <= 1) { - // 4. every quadrant which just has one of each --> done ; - // b) set the sx, sy of the x track now, using the residuals from the - // other b1) special case: no x tracks; then x, sx have been set above - auto xidx = x_idx_quad_map.find((*yitr).first); - auto yidx = y_idx_quad_map.find((*yitr).first); - tracks.at((*xidx).second).setPosition(x, y); - tracks.at((*xidx).second).setSigmaXY(sx, sy); - tracks.at((*yidx).second).setPosition(x, y); - tracks.at((*yidx).second).setSigmaXY(sx, sy); - if (verbose_) - ldmx_log(debug) << "\t\t\t in quad " << (*yitr).first - << ", set (x, y) = (" << x << ", " << y - << ") and (sx, sy) = " << sx << ", " << sy << ")"; - continue; - } - } // 1 y, 0 or 1 x track in quadrant - - if (verbose_) - ldmx_log(debug) << "\t\t in quad " << (*yitr).first - << ", not single x,y tracks: " << n_xin_quad - << " of x and " << n_yin_quad << " of y"; - - if (n_yin_quad == 2) { // let's start here and see if we can do >= 2 later - // here one could do sth to avoid checking the other y track again in the - // outermost loop over y - auto yitr1 = y_quad_map.lower_bound((*yitr).first); - auto yitr2 = y_quad_map.upper_bound((*yitr).first); - yitr2--; // back up once - y1 = ((*yitr1).second).getCentroidY() * y_conv_factor_ + y_start_; - y2 = ((*yitr2).second).getCentroidY() * y_conv_factor_ + y_start_; - sy1 = ((*yitr1).second).getResidual() * y_conv_factor_; - sy2 = ((*yitr2).second).getResidual() * y_conv_factor_; - y = (y1 + y2) / 2.; - sy = fabs(y1 - y2) / 2 * y_conv_factor_; - if (verbose_) - ldmx_log(debug) - << "\t\t -- 2 y in quad: setting x track y coordinate to midpoint"; - } // 2y in quad - - if (n_yin_quad == 1 && - n_xin_quad == 2) { // don't think we want to experiment with discerning - // three overlapping tracks, so not >= 2 - - // first: set the y track coordinates to x = the mid of x tracks, y = y - // of y track - auto yidx = y_idx_quad_map.find((*yitr).first); - tracks.at((*yidx).second).setPosition(x, y); - tracks.at((*yidx).second).setSigmaXY(sx, sy); - - int min_overlap_pe = 250; - if (((*yitr).second).getPE() < min_overlap_pe) { - // can't tell, really, that either of these belong to the y track. so. - // let them keep their own x coordinate but set y to quadrant midpoint, - // with uncertainty +/- half quadrant width (1/8 of pad height) - y = y0; - sy = sy0; - if (verbose_) - ldmx_log(debug) << "\t\t -- Can't tell which x track should be " - "matched to single y track. Setting both x track " - "coordinates to y quadrant value:"; - } // if can't assume overlap - else if (verbose_) - ldmx_log(debug) << "\t\t -- Found large PE count (" - << ((*yitr).second).getPE() << " > " << min_overlap_pe - << "), suggesting overlap! Setting both x track " - "coordinates to y track value:"; - - // consider making two x tracks out if this one, and, anyway have to set - // their average as the y track x cocordinate - // EXPERIMENTAL : apply only to x tracks, which can be disregarded for - // electron counting - if (verbose_) - ldmx_log(debug) << "\t\t -- (x1, x2, y) = (" << x1 << ", " << x2 - << ", " << y << ") and (sx1, sx2, sy) = " << sx1 << ", " - << sx2 << ", " << sy << ")"; - - // now set x track coordinates according to overlap check result - auto xidx1 = x_idx_quad_map.lower_bound((*yitr).first); - auto xidx2 = x_idx_quad_map.upper_bound((*yitr).first); - xidx2--; // upper_bound points to (last+1) element - tracks.at((*xidx1).second).setPosition(x1, y); - tracks.at((*xidx1).second).setSigmaXY(sx1, sy); - tracks.at((*xidx2).second).setPosition(x2, y); - tracks.at((*xidx2).second).setSigmaXY(sx2, sy); - - } // 1 y, 2 x tracks in the quadrant - else if (n_yin_quad == 2 && n_xin_quad == 1) { - // 5b) if there are more y than x: could be an overlap - - // first: set the x track coordinates to x = x of x track, y = the mid of - // y tracks - auto xidx = x_idx_quad_map.find((*yitr).first); - tracks.at((*xidx).second).setPosition(x, y); - tracks.at((*xidx).second).setSigmaXY(sx, sy); - - auto xitr = x_quad_map.lower_bound((*yitr).first); - int min_overlap_pe = 300; - if (((*xitr).second).getPE() < min_overlap_pe) { - if (verbose_) - ldmx_log(debug) - << "\t\t just 1 x track with not-unusual PE in the quad -- can't " - "match; setting mid-point values for x "; - x = x0; - sx = sx0; - } // if can't assume overlap - else { - // consider making two x tracks out if this one, and, anyway have to set - // their average as the y track x cocordinate - // EXPERIMENTAL : apply only to x tracks, which can be disregarded for - // electron counting - if (verbose_) - ldmx_log(debug) << "\t\t -- Found large PE count (" - << ((*xitr).second).getPE() << " > " << min_overlap_pe - << ") in x track, suggesting overlap! Setting both y " - "track coordinates to x track value:"; - } // if can assume overlap - if (verbose_) - ldmx_log(debug) << "\t\t -- (x, y1, y2) = (" << x << ", " << y1 << ", " - << y2 << ") and (sx, sy1, sy2) = " << sx << ", " << sy1 - << ", " << sy2 << ")"; - - auto yidx1 = y_idx_quad_map.lower_bound((*yitr).first); - auto yidx2 = y_idx_quad_map.upper_bound((*yitr).first); - yidx2--; // upper_bound points to next element - tracks.at((*yidx1).second).setPosition(x, y1); - tracks.at((*yidx1).second).setSigmaXY(sx, sy1); - tracks.at((*yidx2).second).setPosition(x, y2); - tracks.at((*yidx2).second).setSigmaXY(sx, sy2); - - } // 2 y and 1 x track in quad - else if (n_yin_quad == 2 && n_xin_quad == 2) { - // MIDPONTS ALL OVER! - auto xidx1 = x_idx_quad_map.lower_bound((*yitr).first); - auto xidx2 = x_idx_quad_map.upper_bound((*yitr).first); - xidx2--; - auto yidx1 = y_idx_quad_map.lower_bound((*yitr).first); - auto yidx2 = y_idx_quad_map.upper_bound((*yitr).first); - yidx2--; - - if (y_idx_quad_map.find((*yitr).first) == y_idx_quad_map.end()) - ldmx_log(error) << "The two y tracks in the same quadrant at " - << (*yitr).first - << " appear to not be found in the y track map! " - "investigate. Note that yidx1.first = " - << (*yidx1).first - << " and yidx2.first = " << (*yidx2).first; - else { - tracks.at((*xidx1).second).setPosition(x1, y); - tracks.at((*xidx1).second).setSigmaXY(sx1, sy); - tracks.at((*xidx2).second).setPosition(x2, y); - tracks.at((*xidx2).second).setSigmaXY(sx2, sy); - - tracks.at((*yidx1).second).setPosition(x, y1); - tracks.at((*yidx1).second).setSigmaXY(sx, sy1); - tracks.at((*yidx2).second).setPosition(x, y2); - tracks.at((*yidx2).second).setSigmaXY(sx, sy2); - - if (verbose_) - ldmx_log(debug) << "\t\t -- in a 2 x 2 situaiton; midpoint y: " << y - << " for both x tracks, midpoint x: " << x - << " for both y tracks"; - } - } // if 2 y, 2 x tracks - - if (n_xin_quad > 2) { - if (verbose_) - ldmx_log(debug) << "\t\t -*-*-*- more than 2 x tracks in the same quad " - "-- nothing done about the x,y coordinates in this " - "situation -- implement if needed!!"; - } - if (n_yin_quad > 2) { - if (verbose_) - ldmx_log(debug) << "\t\t -*-*-*- more than 2 y tracks in the same quad " - "-- nothing done about the x,y coordinates in this " - "situation -- implement if needed!!"; - } - - } // over y tracks - - y_quad_map.clear(); - x_quad_map.clear(); - - // return tracks; -} - -void TrigScintLUTTracker::onProcessStart() { - ldmx_log(debug) << "Process starts!"; - - return; -} - -void TrigScintLUTTracker::onProcessEnd() { - ldmx_log(debug) << "Process ends!"; - - return; -} - -} // namespace trigscint - -DECLARE_PRODUCER(trigscint::TrigScintLUTTracker); - diff --git a/extruthfig.py b/extruthfig.py index 9bac416ffd..d08f02bc7f 100644 --- a/extruthfig.py +++ b/extruthfig.py @@ -42,10 +42,14 @@ cluster.input_collection = digi.output_collection cluster.ampl_weighting = False cluster.clustering_threshold = 3.0 + + +CVA = ldmxcfg.Analyzer.from_file('ClusterViewerAnalyzer.cxx', + needs = ['TrigScint_Event', 'SimCore_Event']) +CVA.pass_name = "truthisolater" p.sequence = [*truth_hits, *digis, *clusters, - ldmxcfg.Analyzer.from_file('ClusterViewerAnalyzer.cxx', - needs = ['TrigScint_Event', 'SimCore_Event']), + CVA ] diff --git a/notruthfig.py b/notruthfig.py index 442824e079..f6bb417d89 100644 --- a/notruthfig.py +++ b/notruthfig.py @@ -10,7 +10,7 @@ from LDMX.Framework import ldmxcfg -p = ldmxcfg.Process('truthisolater') +p = ldmxcfg.Process('nontruthclusters') p.input_files = ["SimSamples.root"] p.output_files = ["AllSamples.root"] #an additional output file called clusters.txt will be created as well @@ -41,13 +41,16 @@ cluster.ampl_weighting = False cluster.clustering_threshold = 3.0 +CVA= ldmxcfg.Analyzer.from_file('ClusterViewerAnalyzer.cxx', + needs = ['TrigScint_Event', 'SimCore_Event']) + +CVA.pass_name = "nontruthclusters" + p.sequence = [ #*truth_hits, *digis, *clusters, - ldmxcfg.Analyzer.from_file('ClusterViewerAnalyzer.cxx', - needs = ['TrigScint_Event', 'SimCore_Event']), - #tracks + CVA ] From f61d6526b15c9224fad656366a60e2bb5a870869 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Mon, 27 Apr 2026 10:54:04 +0200 Subject: [PATCH 07/20] update configs --- TrackFig.py | 66 ++++++++++++----------------------------------------- 1 file changed, 15 insertions(+), 51 deletions(-) diff --git a/TrackFig.py b/TrackFig.py index d1384ee13a..3f9d8c9b3b 100644 --- a/TrackFig.py +++ b/TrackFig.py @@ -1,63 +1,27 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Mon Apr 13 12:17:11 2026 - -@author: kvarn -""" - +#Performs tracking with TrigScintTrackProducer, with and without LUT method for comparison from LDMX.Framework import ldmxcfg p = ldmxcfg.Process('tracking') p.input_files = ["TruthSamples.root"] p.output_files = ["Tracks.root"] -#from LDMX.TrigScint.trig_scint import TrigScintTrackProducer +from LDMX.TrigScint.trig_scint import TrigScintTrackProducer -prod_tracks = ldmxcfg.Producer("lut_track", "trigscint::TrigScintTrackProducer", - "TrigScint") +tstp_tracks = TrigScintTrackProducer("tstp") +lut_tracks = TrigScintTrackProducer("lut") -LUT_tracks = ldmxcfg.Producer("prod_track", "trigscint::TrigScintLUTTracker", - "TrigScint") +tstp_tracks.delta_max = 1.0 +tstp_tracks.verbosity = 1 +tstp_tracks.output_collection = "tstpTracks" -prod_tracks.seeding_collection = "TriggerPad1Clusters" -prod_tracks.further_input_collections = ["TriggerPad2Clusters", "TriggerPad3Clusters"] -prod_tracks.delta_max = 1.0 -prod_tracks.input_pass_name = "" -prod_tracks.verbosity = 1 -prod_tracks.vertical_bar_start_index = 100 -prod_tracks.number_horizontal_bars = 16 -prod_tracks.horizontal_bar_width = 10.0 -prod_tracks.horizontal_bar_gap = 1.0 -prod_tracks.number_vertical_bars = 16 -prod_tracks.vertical_bar_width = 10.0 -prod_tracks.vertical_bar_gap = 1.0 -prod_tracks.allow_skip_last_collection = False +lut_tracks.delta_max = 1.0 +lut_tracks.verbosity = 1 +lut_tracks.lut_tracking = True +lut_tracks.lut_file = "LUT.txt" +lut_tracks.output_collection = "lutTracks" p.logger.term_level = 1 -LUT_tracks.seeding_collection = "TriggerPad1Clusters" -LUT_tracks.further_input_collections = ["TriggerPad2Clusters", "TriggerPad3Clusters"] -LUT_tracks.delta_max = 1.0 -LUT_tracks.input_pass_name = "" -LUT_tracks.verbosity = 1 -LUT_tracks.vertical_bar_start_index = 100 -LUT_tracks.number_horizontal_bars = 16 -LUT_tracks.horizontal_bar_width = 10.0 -LUT_tracks.horizontal_bar_gap = 1.0 -LUT_tracks.number_vertical_bars = 16 -LUT_tracks.vertical_bar_width = 10.0 -LUT_tracks.vertical_bar_gap = 1.0 -LUT_tracks.allow_skip_last_collection = False -LUT_tracks.lut_file = "LUT.txt" - - -prod_tracks.output_collection = "TrigScintTrackProdTracks" -LUT_tracks.output_collection = "TrigScintLUTTracks" - - -p.sequence = [#ldmxcfg.Analyzer.from_file('LUTAnalyzer.cxx', - #needs = ['TrigScint_Event', 'SimCore_Event']), - prod_tracks, - #LUT_tracks - ] \ No newline at end of file +p.sequence = [tstp_tracks, + lut_tracks + ] From fe9106ac7a2579965b6633965485797fe5a636c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Thu, 4 Jun 2026 15:57:04 +0200 Subject: [PATCH 08/20] Rename files for PR --- LUTAnalyzer.cxx | 110 ---------------- NewEventcreator.py | 24 ---- .../include/TrigScint/ClusterTripletMaker.h | 56 ++++++++ TrigScint/include/TrigScint/PatternLUTMaker.h | 65 +++++++++ .../src/TrigScint/ClusterTripletMaker.cxx | 123 ++++++++++++++++++ TrigScint/src/TrigScint/PatternLUTMaker.cxx | 119 +++++++++++++++++ extruthfig.py | 55 -------- 7 files changed, 363 insertions(+), 189 deletions(-) delete mode 100644 LUTAnalyzer.cxx delete mode 100644 NewEventcreator.py create mode 100644 TrigScint/include/TrigScint/ClusterTripletMaker.h create mode 100644 TrigScint/include/TrigScint/PatternLUTMaker.h create mode 100644 TrigScint/src/TrigScint/ClusterTripletMaker.cxx create mode 100644 TrigScint/src/TrigScint/PatternLUTMaker.cxx delete mode 100644 extruthfig.py diff --git a/LUTAnalyzer.cxx b/LUTAnalyzer.cxx deleted file mode 100644 index 07c1bed750..0000000000 --- a/LUTAnalyzer.cxx +++ /dev/null @@ -1,110 +0,0 @@ -//WRITES LUT BASED ON FREQUENCY OF TRACK CANDIDATE PROPAGATION COMBINATIONS -#include "Framework/EventProcessor.h" -#include "Framework/Event.h" -#include "TrigScint/Event/TrigScintTrack.h" -#include -#include -#include -#include - -class LUTAnalyzer : public framework::Analyzer { - public: - - struct Line {int event; float p1, p2, p3, p3alt; - }; - - void configure(framework::config::Parameters& ps) override { - input_file_ = ps.get("input_file"); - output_file_ = ps.get("output_file"); - lut_threshold_ = ps.get("lut_threshold"); - //The LUT threshold is the minimum percentage of times that a track - //must appear in a pool of events to be written to the LUT. - //Example --- straight tracks whose pad 1-2 and pad 2-3 delta values - // are both 0 make up ~85% of tracks out of 10.000 events, - // while a track with deltas (+22,-22), meaning pad 1 cluster in - // e.g. bar 4, pad 2 cluster in bar 26, and pad 3 cluster in bar 4 - // only appears once in 10.000 events (0.01%). The straight - // tracks are written to the LUT and the single "anomaly" is not - - } - - std::ifstream infile; - std::ofstream outfile; - std::string input_file_; - std::string output_file_; - double lut_threshold_; - - std::map, std::vector> groups; - - int totalLines = 0; - - LUTAnalyzer(const std::string& name, framework::Process& p) - : framework::Analyzer(name,p) { - } - void onProcessStart() override { - infile.open(input_file_); - outfile.open(output_file_); - } - - void analyze(const framework::Event& event) override { - if (totalLines > 0) return; - - int ev; - float p1, p2, p3, p3alt; - - while (infile >> ev >> p1 >> p2 >> p3 >> p3alt) { - - float p12 = p2-p1; - float p23 = p3-p2; - - groups[{p12,p23}].push_back({ev,p1,p2,p3,p3alt}); - totalLines++; - } - } - - void onProcessEnd() override { - - int combs = 0; - int tracks = 0; - - std::cout << "Total number of track candidates: " << totalLines << "\n"; - std::cout << "Number of track candidate types: " << groups.size() << "\n"; - std::cout << "LUT Threshold: " << lut_threshold_ * 100 << "%\n"; - - for (auto& g : groups) { - - float p12 = g.first.first; - float p23 = g.first.second; - int count = g.second.size(); - - double frac = (double)count / totalLines; - - std::cout << "(" << p12 << "," << p23 << ") appears " - << count << " times, (" << frac * 100 << " %)" << "\n"; - } - - for (auto& g : groups) { - - int count = g.second.size(); - double frac = (double)count / totalLines; - - if (frac > lut_threshold_) { - - combs++; - - for (auto& line : g.second) { - - tracks++; - - outfile << line.p1 << " " - << line.p2 << " " - << line.p3 << "\n"; - } - } - } - std::cout << "\nLUT textfile written."; - std::cout << "\n" << combs << " combinations (" << tracks << " tracks) written to LUT." << "\n"; - } -}; - -DECLARE_ANALYZER(LUTAnalyzer) diff --git a/NewEventcreator.py b/NewEventcreator.py deleted file mode 100644 index 8c6d9e50ed..0000000000 --- a/NewEventcreator.py +++ /dev/null @@ -1,24 +0,0 @@ -#Creates simulation sample .root file with TS info - -from LDMX.Framework import ldmxcfg -p = ldmxcfg.Process('simulation') - -from LDMX.SimCore import simulator as sim -from LDMX.SimCore import generators as gen - -mySim = sim.simulator( "mySim" ) -mySim.setDetector( 'ldmx-det-v14-8gev', True ) -mySim.generators = [ gen.single_8gev_e_upstream_tagger() ] -mySim.beamSpotSmear = [20.,80.,0.] -mySim.description = 'Basic test Simulation' - -p.sequence = [ mySim ] -p.run = 1 #different for each simulation -p.max_events = 10000 -p.output_files = [ 'SimSamples.root' ] #new output file name - -import LDMX.Ecal.ecal_geometry -import LDMX.Ecal.ecal_hardcoded_conditions -import LDMX.Hcal.hcal_geometry - - diff --git a/TrigScint/include/TrigScint/ClusterTripletMaker.h b/TrigScint/include/TrigScint/ClusterTripletMaker.h new file mode 100644 index 0000000000..8db952139a --- /dev/null +++ b/TrigScint/include/TrigScint/ClusterTripletMaker.h @@ -0,0 +1,56 @@ +#ifndef TRIGSCINT_CLUSTERTRIPLETMAKER_H +#define TRIGSCINT_CLUSTERTRIPLETMAKER_H + +// LDMX Framework +#include "Framework/Configure/Parameters.h" +#include "Framework/Event.h" +#include "Framework/EventProcessor.h" + +// TrigScint +#include "TrigScint/Event/TrigScintCluster.h" + +#include +#include +#include + +namespace trigscint { + +/** + * @class ClusterTripletMaker + * @brief Write trigger scintillator cluster triplets to a text file + */ + +class ClusterTripletMaker : public framework::Analyzer { + public: + ClusterTripletMaker(const std::string& name, framework::Process& process); + + virtual ~ClusterTripletMaker() = default; + + void configure(framework::config::Parameters& ps) override; + + void analyze(const framework::Event& event) override; + + void onProcessStart() override; + + void onProcessEnd() override; + + private: + // verbosity + int verbose_{0}; + + // specific pass name + std::string pass_name_{""}; + + // input cluster collections + std::vector cluster_input_collections_; + + // output text file + std::string output_file_{"clusters.txt"}; + + // output stream + std::ofstream output_stream_; +}; + +} // namespace trigscint + +#endif // TRIGSCINT_CLUSTERTRIPLETMAKER_H diff --git a/TrigScint/include/TrigScint/PatternLUTMaker.h b/TrigScint/include/TrigScint/PatternLUTMaker.h new file mode 100644 index 0000000000..87f72269f9 --- /dev/null +++ b/TrigScint/include/TrigScint/PatternLUTMaker.h @@ -0,0 +1,65 @@ +/** + * @file LUTAnalyzer.h + * @brief Outputs a LUT based on frequency of track propagation patterns + * (PATTERNS, NOT INDIVIDUAL BAR ID COMBINATIONS), which can then be used + * as input in TrigScintTrackProducer for LUT-method track-making. + * @author Lucia Kvarnstrom, Lund University + */ + +#ifndef TRIGSCINT_PATTERNLUTMAKER_H +#define TRIGSCINT_PATTERNLUTMAKER_H +#include "Framework/Configure/Parameters.h" +#include "Framework/EventProcessor.h" + +#include +#include +#include +#include +#include + +namespace trigscint { + +class PatternLUTMaker : public framework::Analyzer { + public: + struct Line { + int event; + float p1, p2, p3; + }; + + PatternLUTMaker(const std::string& name, framework::Process& process); + + virtual ~PatternLUTMaker() = default; + + void configure(framework::config::Parameters& parameters) override; + + void analyze(const framework::Event& event) override; + + void onProcessStart() override; + + void onProcessEnd() override; + + private: + // verbosity + int verbose_{0}; + + // input text file of clusters + std::string input_file_; + + //output LUT file name + std::string output_file_; + + //minimum frequency of a specific pattern to be written to the LUT + double lut_threshold_{0.0008}; + + std::ifstream infile; + std::ofstream outfile; + + // to group cluster combinations by track pattern + std::map, std::vector> groups; + + int totalLines{0}; +}; + +} // namespace trigscint + +#endif /* TRIGSCINT_PATTERNLUTMAKER_H */ diff --git a/TrigScint/src/TrigScint/ClusterTripletMaker.cxx b/TrigScint/src/TrigScint/ClusterTripletMaker.cxx new file mode 100644 index 0000000000..4a1d42d238 --- /dev/null +++ b/TrigScint/src/TrigScint/ClusterTripletMaker.cxx @@ -0,0 +1,123 @@ +#include "TrigScint/ClusterTripletMaker.h" + +#include + +namespace trigscint { + +void ClusterTripletMaker::configure(framework::config::Parameters& ps) { + pass_name_ = ps.get("pass_name"); + cluster_input_collections_ = ps.get>("cluster_input_collections"); + output_file_ = ps.get("output_file"); + verbose_ = ps.get("verbosity"); + + if (verbose_) { + ldmx_log(info) + << "In ClusterTripletMaker: configure done!" + << "\nInput collection 1: " + << cluster_input_collections_.at(0) + << "\nInput collection 2: " + << cluster_input_collections_.at(1) + << "\nInput collection 3: " + << cluster_input_collections_.at(2) + << "\nPass name: " << pass_name_ + << "\nOutput file: " << output_file_ + << "\nVerbosity: " << verbose_; + } + + return; +} + +void ClusterTripletMaker::onProcessStart() { + output_stream_.open(output_file_); + + if (verbose_) { + ldmx_log(info) << "Opened output file " << output_file_; + } + + return; +} + +void ClusterTripletMaker::analyze( + const framework::Event& event) { + const auto clusters_pad1{event.getCollection( + cluster_input_collections_.at(0), pass_name_)}; + + const auto clusters_pad2{event.getCollection( + cluster_input_collections_.at(1), pass_name_)}; + + const auto clusters_pad3{event.getCollection( + cluster_input_collections_.at(2), pass_name_)}; + + const size_t num_clusters =std::min( + {clusters_pad1.size(), clusters_pad2.size(), clusters_pad3.size()}); + + if (verbose_ > 1) { + ldmx_log(debug) + << "Event " << event.getEventNumber() + << " has " << clusters_pad1.size() + << ", " << clusters_pad2.size() << ", " + << clusters_pad3.size() + << " clusters in pads 1, 2, and 3."; + } + + for (size_t cluster_index = 0; cluster_index < num_clusters; ++cluster_index) { + const int event_number = event.getEventNumber(); + + const float pad1_centroid = clusters_pad1.at(cluster_index).getCentroid(); + const float pad2_centroid = clusters_pad2.at(cluster_index).getCentroid(); + const float pad3_centroid = clusters_pad3.at(cluster_index).getCentroid(); + + output_stream_ + << event_number << " " + << pad1_centroid << " " + << pad2_centroid << " " + << pad3_centroid << " " + << "\n"; + } + + if (verbose_) { + for (size_t cluster_index = 1; + cluster_index < clusters_pad1.size(); ++cluster_index) { + ldmx_log(info) + << "Event " << event.getEventNumber() + << ", extra cluster in pad 1: " + << clusters_pad1.at(cluster_index).getCentroid(); + } + + for (size_t cluster_index = 1; + cluster_index < clusters_pad2.size(); ++cluster_index) { + ldmx_log(info) + << "Event " << event.getEventNumber() + << ", extra cluster in pad 2: " + << clusters_pad2.at(cluster_index).getCentroid(); + } + + for (size_t cluster_index = 1; + cluster_index < clusters_pad3.size(); ++cluster_index) { + ldmx_log(info) + << "Event " << event.getEventNumber() + << ", extra cluster in pad 3: " + << clusters_pad3.at(cluster_index).getCentroid(); + } + } + + return; +} + +void ClusterTripletMaker::onProcessEnd() { + if (output_stream_.is_open()) { + output_stream_.close(); + } + + if (verbose_) { + ldmx_log(info) + << "Closed output file " + << output_file_; + } + + return; +} + +} // namespace trigscint + +DECLARE_ANALYZER(trigscint::ClusterTripletMaker); diff --git a/TrigScint/src/TrigScint/PatternLUTMaker.cxx b/TrigScint/src/TrigScint/PatternLUTMaker.cxx new file mode 100644 index 0000000000..27d17cc538 --- /dev/null +++ b/TrigScint/src/TrigScint/PatternLUTMaker.cxx @@ -0,0 +1,119 @@ +//Analyzer file which writes a LUT text file from input cluster text file (made in +//ClusterTripletMaker) containing only those cluster combination entries which +//occur with a frequency above the lut_threshold parameter + +#include "TrigScint/PatternLUTMaker.h" + +namespace trigscint { + +PatternLUTMaker::PatternLUTMaker(const std::string& name, + framework::Process& process) + : Analyzer(name, process) {} + +void PatternLUTMaker::configure(framework::config::Parameters& ps) { + input_file_ = ps.get("input_file"); + output_file_ = ps.get("output_file"); + lut_threshold_ = ps.get("lut_threshold"); + verbose_ = ps.get("verbosity"); + + //The LUT threshold is the minimum percentage of times that a track + //must appear in a pool of events to be written to the LUT. + //Example --- horizontal tracks whose pad 1-2 and pad 2-3 delta values + // are both 0 make up ~80% of tracks out of 10.000 events, + // while a track with deltas (+22,-22), meaning pad 1 cluster in + // e.g. bar 4, pad 2 cluster in bar 26, and pad 3 cluster in bar 4 + // only appears once in 10.000 events (0.01%). The straight + // tracks are written to the LUT and the single "anomaly" is not. + //Primary motivation here is for use in the case of an unknown TS misalignment. + + if (verbose_) { + ldmx_log(info) << "In PatternLUTMaker: configure done!" << std::endl; + ldmx_log(info) << "Got parameters: \nInput file: " << input_file_ + << "\nOutput file: " << output_file_ + << "\nLUT threshold: " << lut_threshold_ + << "\nVerbosity: " << verbose_; + } + return; +} + +void PatternLUTMaker::onProcessStart() { + infile.open(input_file_); + outfile.open(output_file_); + return; +} + +void PatternLUTMaker::analyze(const framework::Event& event) { + if (totalLines > 0) return; + + int ev; + float p1, p2, p3; + + while (infile >> ev >> p1 >> p2 >> p3) { + float p12 = p2 - p1; + float p23 = p3 - p2; + + groups[{p12, p23}].push_back({ev, p1, p2, p3}); + + totalLines++; + } + + return; +} + +void PatternLUTMaker::onProcessEnd() { + int combs = 0; + int tracks = 0; + + if(verbose_){ + ldmx_log(info) << "Total number of track candidates: " + << totalLines << "\n" + << "Number of track candidate types: " + << groups.size() << "\n" + << "LUT Threshold: " + << lut_threshold_ * 100 << "%\n"; + } + + for (auto& g : groups) { + float p12 = g.first.first; + float p23 = g.first.second; + int count = g.second.size(); + + double frac = static_cast(count) / totalLines; + + if (verbose_) { + ldmx_log(info) << "(" << p12 << "," << p23 << ") appears " + << count << " times, (" << frac * 100 << " %)" + << "\n"; + } + } + + for (auto& g : groups) { + int count = g.second.size(); + + double frac = static_cast(count) / totalLines; + + if (frac > lut_threshold_) { //write to outfile if over threshold + combs++; + + for (auto& line : g.second) { + tracks++; + + outfile << line.p1 << " " + << line.p2 << " " + << line.p3 << "\n"; + } + } + } + + if (verbose_) { + ldmx_log(info) << "\nLUT textfile written." + << "\n" << combs << " combinations (" + << tracks << " tracks) written to LUT.\n"; + } + + return; +} + +} // namespace trigscint + +DECLARE_ANALYZER(trigscint::PatternLUTMaker) diff --git a/extruthfig.py b/extruthfig.py deleted file mode 100644 index d08f02bc7f..0000000000 --- a/extruthfig.py +++ /dev/null @@ -1,55 +0,0 @@ -#Uses TruthHitProducer to isolate true hits, -#then turns them to digis, then clusters, then makes clusters.txt - -from LDMX.Framework import ldmxcfg - -p = ldmxcfg.Process('truthisolater') -p.input_files = ["SimSamples.root"] -p.output_files = ["TruthSamples.root"] -#an additional output file called clusters.txt will be created as well - -from LDMX.TrigScint.truth_hits import TruthHitProducer -from LDMX.TrigScint.trig_scint import TrigScintDigiProducer - -truth_hits = [TruthHitProducer('beamElectronsPad1'), - TruthHitProducer('beamElectronsPad2'), - TruthHitProducer('beamElectronsPad3') - ] - -pad_num = 1 - -for hits in truth_hits: - hits.input_collection=f"TriggerPad{pad_num}SimHits" - hits.output_collection=f"truthBeamElectronsPad{pad_num}" - pad_num+=1 - -digis = [TrigScintDigiProducer.pad1(), - TrigScintDigiProducer.pad2(), - TrigScintDigiProducer.pad3() - ] - -for digi,hits in zip(digis, truth_hits): - digi.input_collection = hits.output_collection - -from LDMX.TrigScint.trig_scint import TrigScintClusterProducer - -clusters = [TrigScintClusterProducer.pad1(), - TrigScintClusterProducer.pad2(), - TrigScintClusterProducer.pad3(), - ] - -for cluster, digi in zip(clusters, digis): - cluster.input_collection = digi.output_collection - cluster.ampl_weighting = False - cluster.clustering_threshold = 3.0 - - -CVA = ldmxcfg.Analyzer.from_file('ClusterViewerAnalyzer.cxx', - needs = ['TrigScint_Event', 'SimCore_Event']) -CVA.pass_name = "truthisolater" - -p.sequence = [*truth_hits, - *digis, - *clusters, - CVA - ] From 823b30f478f983e9a5adc7e4308a49c97c514d6c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Thu, 4 Jun 2026 16:08:45 +0200 Subject: [PATCH 09/20] Add new analyzers and rm extra TSTrackProd --- ClusterViewerAnalyzer.cxx | 79 ---------------------------------- TrigScint/python/trig_scint.py | 22 ++++++++++ 2 files changed, 22 insertions(+), 79 deletions(-) delete mode 100644 ClusterViewerAnalyzer.cxx diff --git a/ClusterViewerAnalyzer.cxx b/ClusterViewerAnalyzer.cxx deleted file mode 100644 index 2ece585615..0000000000 --- a/ClusterViewerAnalyzer.cxx +++ /dev/null @@ -1,79 +0,0 @@ -//ANALYZER TO VIEW CLUSTERS - -#include "Framework/EventProcessor.h" -#include "Framework/Event.h" -#include "TrigScint/Event/TrigScintCluster.h" -#include -#include - - -class ClusterViewerAnalyzer : public framework::Analyzer { - public: - - std::ofstream outfile; - - void configure(framework::config::Parameters& ps) override { - pass_name_ = ps.get("pass_name"); - } - - std::string pass_name_; - - ClusterViewerAnalyzer(const std::string& name, framework::Process& p) - : framework::Analyzer(name, p) { - - outfile.open("clusters.txt"); - } - - void analyze(const framework::Event& event) override { - - const auto& clusters1 = - event.getCollection("TriggerPad1Clusters", pass_name_); - - const auto& clusters2 = - event.getCollection("TriggerPad2Clusters", pass_name_); - - const auto& clusters3 = - event.getCollection("TriggerPad3Clusters", pass_name_); - - size_t Clusters = std::min({clusters1.size(), clusters2.size(), clusters3.size()}); - - for (size_t i = 0; i < Clusters; i++) { - - int eventNumber = event.getEventNumber(); - float c1 = clusters1[i].getCentroid(); - float c2 = clusters2[i].getCentroid(); - float c3 = clusters3[i].getCentroid(); - float c3_alt = c3; - - if (clusters3.size() > i + 1) { - c3_alt = clusters3[i + 1].getCentroid(); - } - - outfile << eventNumber << " " - << c1 << " " - << c2 << " " - << c3 << " " - << c3_alt << "\n"; - } - - - for (size_t i = 1; i < clusters1.size(); i++){ - std::cout << "In event " << event.getEventNumber() <<", extra cluster in pad 1: " - << clusters1[i].getCentroid() << "\n"; - } - - for (size_t i = 1; i < clusters2.size(); i++){ - std::cout << "In event " << event.getEventNumber() <<", extra cluster in pad 2: " - << clusters2[i].getCentroid() << "\n"; - } - - for (size_t i = 1; i < clusters3.size(); i++){ - std::cout << "In event " << event.getEventNumber() <<", extra cluster in pad 3: " - << clusters3[i].getCentroid() << "\n"; - } - - } - void onProcessEnd() override {outfile.close();} - }; - -DECLARE_ANALYZER(ClusterViewerAnalyzer) diff --git a/TrigScint/python/trig_scint.py b/TrigScint/python/trig_scint.py index 393fcd8e9d..0ce0380916 100644 --- a/TrigScint/python/trig_scint.py +++ b/TrigScint/python/trig_scint.py @@ -467,3 +467,25 @@ def __init__(self,name) : self.start_sample=2 #first time sample included in reformatting self.deadChannels=[ 8 ] + +class ClusterTripletMaker(ldmxcfg.Analyzer) : + """Configuration for cluster text file maker for Trigger Scintillators""" + + def __init__(self,name) : + super().__init__(name,'trigscint::ClusterTripletMaker','TrigScint') + + self.cluster_input_collections = ["TriggerPad1Clusters", "TriggerPad2Clusters","TriggerPad3Clusters"] + self.input_pass_name="" #take any pass + self.output_collection="clusters.txt" + + +class PatternLUTMaker(ldmxcfg.Analyzer) : + """Configuration for track-pattern LUT-writing analyzer for Trigger Scintillators""" + + def __init__(self,name) : + super().__init__(name,'trigscint::PatternLUTMaker','TrigScint') + + self.input_collection="clusters.txt" + self.input_pass_name="" #take any pass + self.output_collection="LUT.txt" + self.lut_threshold=0.0008 From 34378da1037a334c1ab0a82738fc066214f35764 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Thu, 4 Jun 2026 16:33:31 +0200 Subject: [PATCH 10/20] Update configs --- notruthfig.py => TrigScint/exampleConfigs/runClusterstxt.py | 0 LUTconfig.py => TrigScint/exampleConfigs/runLUTana.py | 0 TrackFig.py => TrigScint/exampleConfigs/runLUTtracking.py | 0 3 files changed, 0 insertions(+), 0 deletions(-) rename notruthfig.py => TrigScint/exampleConfigs/runClusterstxt.py (100%) rename LUTconfig.py => TrigScint/exampleConfigs/runLUTana.py (100%) rename TrackFig.py => TrigScint/exampleConfigs/runLUTtracking.py (100%) diff --git a/notruthfig.py b/TrigScint/exampleConfigs/runClusterstxt.py similarity index 100% rename from notruthfig.py rename to TrigScint/exampleConfigs/runClusterstxt.py diff --git a/LUTconfig.py b/TrigScint/exampleConfigs/runLUTana.py similarity index 100% rename from LUTconfig.py rename to TrigScint/exampleConfigs/runLUTana.py diff --git a/TrackFig.py b/TrigScint/exampleConfigs/runLUTtracking.py similarity index 100% rename from TrackFig.py rename to TrigScint/exampleConfigs/runLUTtracking.py From 9c4667f8fc381823e378aa9094eaf449294ea827 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Thu, 4 Jun 2026 18:55:20 +0200 Subject: [PATCH 11/20] formatting fixes --- TrigScint/include/TrigScint/Event/TrigScintCluster.h | 6 +++--- TrigScint/include/TrigScint/TrigScintClusterProducer.h | 3 --- TrigScint/src/TrigScint/TrigScintClusterProducer.cxx | 3 --- 3 files changed, 3 insertions(+), 9 deletions(-) diff --git a/TrigScint/include/TrigScint/Event/TrigScintCluster.h b/TrigScint/include/TrigScint/Event/TrigScintCluster.h index a6839ae3eb..bfedb52210 100644 --- a/TrigScint/include/TrigScint/Event/TrigScintCluster.h +++ b/TrigScint/include/TrigScint/Event/TrigScintCluster.h @@ -93,7 +93,7 @@ class TrigScintCluster { /** * @param centroid The channel ID centroid */ - void setCentroid(float centroid) { centroid_ = centroid; } + void setCentroid(double centroid) { centroid_ = centroid; } /** Set time of hit. */ void setTime(float t) { time_ = t; } @@ -132,7 +132,7 @@ class TrigScintCluster { const std::vector &getHitIDs() const { return hit_ids_; } /** Get the cluster centroid in units of channel nb */ - float getCentroid() const { return centroid_; } + double getCentroid() const { return centroid_; } bool operator<(const TrigScintCluster &rhs) const { return this->getEnergy() < rhs.getEnergy(); @@ -156,7 +156,7 @@ class TrigScintCluster { // hit centroid in units of channel nb: energy weighted average of the IDs of // the hits forming the cluster - float centroid_{-1}; + double centroid_{-1}; // hit centroid in x [mm] (not implemented) double centroid_x_{0}; diff --git a/TrigScint/include/TrigScint/TrigScintClusterProducer.h b/TrigScint/include/TrigScint/TrigScintClusterProducer.h index 92b2149f4d..b09eeb3bd6 100644 --- a/TrigScint/include/TrigScint/TrigScintClusterProducer.h +++ b/TrigScint/include/TrigScint/TrigScintClusterProducer.h @@ -96,9 +96,6 @@ class TrigScintClusterProducer : public framework::Producer { // sum of hit cluster weights float sumw_{0.}; - // sum of hit cluster weights - float sumw_{0.}; - // book keep which channels have already been added to the cluster at hand std::vector v_added_indices_; diff --git a/TrigScint/src/TrigScint/TrigScintClusterProducer.cxx b/TrigScint/src/TrigScint/TrigScintClusterProducer.cxx index ae52966acf..58a7f8557f 100644 --- a/TrigScint/src/TrigScint/TrigScintClusterProducer.cxx +++ b/TrigScint/src/TrigScint/TrigScintClusterProducer.cxx @@ -404,7 +404,6 @@ void TrigScintClusterProducer::produce(framework::Event &event) { } // if adding another hit, going forward, was allowed // done adding hits to cluster. calculate centroid - centroid_ /= sumw_; // final weighting step: divide by total amplitude sum @@ -483,7 +482,6 @@ void TrigScintClusterProducer::produce(framework::Event &event) { void TrigScintClusterProducer::addHit(uint idx, ldmx::TrigScintHit hit) { float ampl = hit.getPE(); float w = 1; - if (ampl_weighting_) { // if choosing to PE-weight centroid positions w = ampl; } @@ -492,7 +490,6 @@ void TrigScintClusterProducer::addHit(uint idx, ldmx::TrigScintHit hit) { val_e_ += energy; val_ += ampl; - centroid_ += (idx + 1) * w; // need non-zero weight of channel 0. shifting // centroid back by 1 in the end // this number gets divided by val at the end From 9bd018f40c58d3ced39ccc95c246f639a135ca4d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Thu, 4 Jun 2026 18:59:45 +0200 Subject: [PATCH 12/20] formatting fixes --- TrigScint/include/TrigScint/PatternLUTMaker.h | 2 +- TrigScint/include/TrigScint/TrigScintClusterProducer.h | 2 +- TrigScint/src/TrigScint/TrigScintTrackProducer.cxx | 2 -- 3 files changed, 2 insertions(+), 4 deletions(-) diff --git a/TrigScint/include/TrigScint/PatternLUTMaker.h b/TrigScint/include/TrigScint/PatternLUTMaker.h index 87f72269f9..9f5c064fa4 100644 --- a/TrigScint/include/TrigScint/PatternLUTMaker.h +++ b/TrigScint/include/TrigScint/PatternLUTMaker.h @@ -1,5 +1,5 @@ /** - * @file LUTAnalyzer.h + * @file PatternLUTMaker.h * @brief Outputs a LUT based on frequency of track propagation patterns * (PATTERNS, NOT INDIVIDUAL BAR ID COMBINATIONS), which can then be used * as input in TrigScintTrackProducer for LUT-method track-making. diff --git a/TrigScint/include/TrigScint/TrigScintClusterProducer.h b/TrigScint/include/TrigScint/TrigScintClusterProducer.h index b09eeb3bd6..a0ed9a242e 100644 --- a/TrigScint/include/TrigScint/TrigScintClusterProducer.h +++ b/TrigScint/include/TrigScint/TrigScintClusterProducer.h @@ -92,7 +92,7 @@ class TrigScintClusterProducer : public framework::Producer { // edep content, only; leave val_ for PE float val_e_{0.}; - + // sum of hit cluster weights float sumw_{0.}; diff --git a/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx b/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx index 6ad0e0eb77..ac0db5dcae 100644 --- a/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx +++ b/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx @@ -19,7 +19,6 @@ void TrigScintTrackProducer::configure(framework::config::Parameters &ps) { "further_input_collections"); // {"TriggerPadUpClusters" , // "TriggerPadDnClusters" } output_collection_ = ps.get("output_collection"); - pass_name_ = ps.get("input_pass_name"); verbose_ = ps.get("verbosity"); vert_bar_start_idx_ = ps.get("vertical_bar_start_index"); @@ -175,7 +174,6 @@ void TrigScintTrackProducer::produce(framework::Event &event) { // reset for each seed // bool madeTrack = false; -<<<<<<< HEAD if (lut_tracking_) { //if using LUT method for (const auto &cluster1 : clusters_pad1) { From 84d4ad9310b910cdaa2d102b75e8f8ae0d0d009b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Thu, 4 Jun 2026 20:39:56 +0200 Subject: [PATCH 13/20] separate LUT and deltamax tracking --- .../src/TrigScint/TrigScintTrackProducer.cxx | 69 ------------------- 1 file changed, 69 deletions(-) diff --git a/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx b/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx index ac0db5dcae..799c825e8e 100644 --- a/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx +++ b/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx @@ -270,75 +270,6 @@ void TrigScintTrackProducer::produce(framework::Event &event) { // only make this vector now! this ensures against hanging // clusters with indices from earlier in the loop - for (const auto &cluster1 : clusters_pad1) { - if (verbose_ > 1) { - ldmx_log(debug) << "\tGot pad1 cluster with centroid " - << cluster1.getCentroid(); - } - if (fabs(cluster1.getCentroid() - centroid) < max_delta_ || - (centroid >= vert_bar_start_idx_ // then in vertical bars - && seed.getCentroidX() == cluster1.getCentroidX())) { - // use geometry y overlap scheme to see if this is really a match in x - // should be done in a map - - if (centroid >= vert_bar_start_idx_ && - seed.getCentroidY() < cluster1.getCentroidY()) { - // impossible combination - if (verbose_ > 1) { - ldmx_log(debug) << "\tSkipping impossible x cluster combination " - "with y flags (tag up) (" - << seed.getCentroidY() << " " - << cluster1.getCentroidY() << ")"; - } - continue; - } - - // else: first (possible) match! loop through next pad too - - if (verbose_ > 1) { - ldmx_log(debug) << "\t\tIt is close enough!. Check pad2"; - } - - // try making third pad clusters an optional part of track - - std::vector cluster_vec = {seed, cluster1}; - - bool has_match_dn = false; - - for (const auto &cluster2 : clusters_pad2) { - if (verbose_ > 1) { - ldmx_log(debug) << "\tGot pad2 cluster with centroid " - << cluster2.getCentroid(); - } - - if (fabs(cluster2.getCentroid() - centroid) < max_delta_ || - (centroid >= vert_bar_start_idx_ // then in vertical bars - && seed.getCentroidX() == cluster2.getCentroidX())) { - // use geometry y overlap scheme to see if this is really a match - // in x - - if (centroid >= vert_bar_start_idx_ && - (seed.getCentroidY() < cluster2.getCentroidY() || - cluster1.getCentroidY() > - cluster2.getCentroidY())) { // impossible - if (verbose_ > 1) { - ldmx_log(debug) - << "\tSkipping impossible x cluster combination with y " - "flags (tag up dn) (" - << seed.getCentroidY() << " " << cluster1.getCentroidY() - << " " << cluster2.getCentroidY() << ")"; - } - continue; - } - - // first match! loop through next pad too - - if (verbose_ > 1) { - ldmx_log(debug) << "\t\tIt is close enough!. Make a track"; - } - - // only make this vector now! this ensures against hanging - // clusters with indices from earlier in the loop std::vector three_cluster_vec = { seed, cluster1, cluster2}; From 8981aeb2b3b738f2434ac5a10ea059dba5b2e25d Mon Sep 17 00:00:00 2001 From: kvarnla Date: Thu, 4 Jun 2026 22:57:31 +0200 Subject: [PATCH 14/20] Update runClusterstxt.py --- TrigScint/exampleConfigs/runClusterstxt.py | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/TrigScint/exampleConfigs/runClusterstxt.py b/TrigScint/exampleConfigs/runClusterstxt.py index f6bb417d89..da82ccfe7e 100644 --- a/TrigScint/exampleConfigs/runClusterstxt.py +++ b/TrigScint/exampleConfigs/runClusterstxt.py @@ -1,12 +1,5 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Apr 14 21:58:25 2026 - -@author: kvarn -""" - -#digis, then clusters, then makes clusters.txt +#example config which makes digis, clusters, and writes clusters +#to .txt file using ClusterTripletMaker from LDMX.Framework import ldmxcfg @@ -41,16 +34,15 @@ cluster.ampl_weighting = False cluster.clustering_threshold = 3.0 -CVA= ldmxcfg.Analyzer.from_file('ClusterViewerAnalyzer.cxx', - needs = ['TrigScint_Event', 'SimCore_Event']) +from LDMX.TrigScint.trig_scint import ClusterTripletMaker -CVA.pass_name = "nontruthclusters" +triplets = ClusterTripletMaker("tripletmaker") p.sequence = [ #*truth_hits, *digis, *clusters, - CVA + triplets ] From 6a739f9426e3dc244a80d68b8199712324ff8c8a Mon Sep 17 00:00:00 2001 From: kvarnla Date: Thu, 4 Jun 2026 23:01:13 +0200 Subject: [PATCH 15/20] Update runLUTana.py --- TrigScint/exampleConfigs/runLUTana.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/TrigScint/exampleConfigs/runLUTana.py b/TrigScint/exampleConfigs/runLUTana.py index fe2dfd1e3d..9a5f970a1f 100644 --- a/TrigScint/exampleConfigs/runLUTana.py +++ b/TrigScint/exampleConfigs/runLUTana.py @@ -1,14 +1,15 @@ -#config file to make LUT for tracking +#config file to make LUT for tracking using PatternLUTMaker from LDMX.Framework import ldmxcfg p = ldmxcfg.Process("LUTmaker") -p.input_files = ["TruthSamples.root"] #necessary to define but it doesnt use it +p.input_files = ["TruthSamples.root"] #necessary to define but it doesnt use it, the real input file is defined below -make_lut = ldmxcfg.Analyzer.from_file('LUTAnalyzer.cxx', - needs = ['TrigScint_Event', 'SimCore_Event']) -make_lut.lut_threshold = 1.0/1000 -make_lut.input_file = "clusters.txt" +from LDMX.TrigScint.trig_scint import PatternLUTMaker + +make_lut = PatternLUTMaker("LUTmaker") +make_lut.lut_threshold = 1.0/1000 #or whichever threshold you like +make_lut.input_file = "clusters.txt" #here! input list of cluster triplets make_lut.output_file = "LUT.txt" p.sequence = [make_lut] From f08c3cd350534c2e15d58c48c5a7374d94329538 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Fri, 5 Jun 2026 09:32:33 +0200 Subject: [PATCH 16/20] update analyzer names --- TrigScint/python/trig_scint.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/TrigScint/python/trig_scint.py b/TrigScint/python/trig_scint.py index 2b90118501..87adecd633 100644 --- a/TrigScint/python/trig_scint.py +++ b/TrigScint/python/trig_scint.py @@ -508,7 +508,7 @@ class TestBeamClusterAnalyzer(Processor): @processor("trigscint::ClusterTripletMaker", "TrigScint") -class ClusterTripletMaker(ldmxcfg.Analyzer) : +class ClusterTripletMaker(Processor) : """Configuration for cluster text file maker for Trigger Scintillators""" cluster_input_collections: list[str] = ["TriggerPad1Clusters", "TriggerPad2Clusters","TriggerPad3Clusters"] @@ -517,7 +517,7 @@ class ClusterTripletMaker(ldmxcfg.Analyzer) : @processor("trigscint::PatternLUTMaker", "TrigScint") -class PatternLUTMaker(ldmxcfg.Analyzer) : +class PatternLUTMaker(Processor) : """Configuration for track-pattern LUT-writing analyzer for Trigger Scintillators""" input_collection: str ="clusters.txt" From 34a5045f65e72e639cae0a672532511e1c24179f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lucia=20Kvarnstr=C3=B6m?= Date: Sat, 6 Jun 2026 19:57:26 +0200 Subject: [PATCH 17/20] fix parameter names and library errors --- TrigScint/exampleConfigs/runClusterstxt.py | 11 ++++++----- TrigScint/exampleConfigs/runLUTana.py | 6 +++--- TrigScint/exampleConfigs/runLUTtracking.py | 2 +- TrigScint/include/TrigScint/ClusterTripletMaker.h | 2 +- TrigScint/include/TrigScint/PatternLUTMaker.h | 8 ++++---- TrigScint/python/trig_scint.py | 7 ++++++- TrigScint/src/TrigScint/ClusterTripletMaker.cxx | 14 +++++++++----- TrigScint/src/TrigScint/PatternLUTMaker.cxx | 12 ++++++------ 8 files changed, 36 insertions(+), 26 deletions(-) diff --git a/TrigScint/exampleConfigs/runClusterstxt.py b/TrigScint/exampleConfigs/runClusterstxt.py index da82ccfe7e..778dc1157c 100644 --- a/TrigScint/exampleConfigs/runClusterstxt.py +++ b/TrigScint/exampleConfigs/runClusterstxt.py @@ -5,7 +5,7 @@ p = ldmxcfg.Process('nontruthclusters') p.input_files = ["SimSamples.root"] -p.output_files = ["AllSamples.root"] +p.output_files = ["Clusters.root"] #an additional output file called clusters.txt will be created as well from LDMX.TrigScint.trig_scint import TrigScintDigiProducer @@ -20,9 +20,9 @@ for digi in digis: digi.input_collection = f"TriggerPad{pad_num}SimHits" pad_num+=1 - + from LDMX.TrigScint.trig_scint import TrigScintClusterProducer - + clusters = [ TrigScintClusterProducer.pad1(), TrigScintClusterProducer.pad2(), @@ -37,12 +37,13 @@ from LDMX.TrigScint.trig_scint import ClusterTripletMaker triplets = ClusterTripletMaker("tripletmaker") +triplets.output_collection = "clusters.txt" p.sequence = [ #*truth_hits, *digis, - *clusters, + *clusters, triplets ] - + diff --git a/TrigScint/exampleConfigs/runLUTana.py b/TrigScint/exampleConfigs/runLUTana.py index 9a5f970a1f..468a4eab7c 100644 --- a/TrigScint/exampleConfigs/runLUTana.py +++ b/TrigScint/exampleConfigs/runLUTana.py @@ -3,13 +3,13 @@ from LDMX.Framework import ldmxcfg p = ldmxcfg.Process("LUTmaker") -p.input_files = ["TruthSamples.root"] #necessary to define but it doesnt use it, the real input file is defined below +p.input_files = ["Clusters.root"] #necessary to define but unused, the real input file is defined below from LDMX.TrigScint.trig_scint import PatternLUTMaker make_lut = PatternLUTMaker("LUTmaker") make_lut.lut_threshold = 1.0/1000 #or whichever threshold you like -make_lut.input_file = "clusters.txt" #here! input list of cluster triplets -make_lut.output_file = "LUT.txt" +make_lut.input_collection = "clusters.txt" #here! input list of cluster triplets as produced by ClusterTripletMaker +make_lut.output_collection = "LUT.txt" p.sequence = [make_lut] diff --git a/TrigScint/exampleConfigs/runLUTtracking.py b/TrigScint/exampleConfigs/runLUTtracking.py index 3f9d8c9b3b..68bab3cf6a 100644 --- a/TrigScint/exampleConfigs/runLUTtracking.py +++ b/TrigScint/exampleConfigs/runLUTtracking.py @@ -2,7 +2,7 @@ from LDMX.Framework import ldmxcfg p = ldmxcfg.Process('tracking') -p.input_files = ["TruthSamples.root"] +p.input_files = ["Clusters.root"] p.output_files = ["Tracks.root"] from LDMX.TrigScint.trig_scint import TrigScintTrackProducer diff --git a/TrigScint/include/TrigScint/ClusterTripletMaker.h b/TrigScint/include/TrigScint/ClusterTripletMaker.h index 8db952139a..b06d6753bc 100644 --- a/TrigScint/include/TrigScint/ClusterTripletMaker.h +++ b/TrigScint/include/TrigScint/ClusterTripletMaker.h @@ -45,7 +45,7 @@ class ClusterTripletMaker : public framework::Analyzer { std::vector cluster_input_collections_; // output text file - std::string output_file_{"clusters.txt"}; + std::string output_collection_{"clusters.txt"}; // output stream std::ofstream output_stream_; diff --git a/TrigScint/include/TrigScint/PatternLUTMaker.h b/TrigScint/include/TrigScint/PatternLUTMaker.h index 9f5c064fa4..31f64f91a6 100644 --- a/TrigScint/include/TrigScint/PatternLUTMaker.h +++ b/TrigScint/include/TrigScint/PatternLUTMaker.h @@ -1,8 +1,8 @@ /** * @file PatternLUTMaker.h * @brief Outputs a LUT based on frequency of track propagation patterns - * (PATTERNS, NOT INDIVIDUAL BAR ID COMBINATIONS), which can then be used - * as input in TrigScintTrackProducer for LUT-method track-making. + * which can then be used as input in TrigScintTrackProducer for + * LUT-method track-making. * @author Lucia Kvarnstrom, Lund University */ @@ -43,10 +43,10 @@ class PatternLUTMaker : public framework::Analyzer { int verbose_{0}; // input text file of clusters - std::string input_file_; + std::string input_collection_; //output LUT file name - std::string output_file_; + std::string output_collection_; //minimum frequency of a specific pattern to be written to the LUT double lut_threshold_{0.0008}; diff --git a/TrigScint/python/trig_scint.py b/TrigScint/python/trig_scint.py index 87adecd633..f5b939b09c 100644 --- a/TrigScint/python/trig_scint.py +++ b/TrigScint/python/trig_scint.py @@ -391,6 +391,9 @@ class TrigScintTrackProducer(Processor): input_pass_name: str = "" output_collection: str = "TriggerPadTracks" verbosity: int = 0 + lut_tracking: bool = False + lut_file: str = "LUT.txt" + trig_scint_track = TrigScintTrackProducer(instance_name="trig_scint_track") @@ -512,8 +515,9 @@ class ClusterTripletMaker(Processor) : """Configuration for cluster text file maker for Trigger Scintillators""" cluster_input_collections: list[str] = ["TriggerPad1Clusters", "TriggerPad2Clusters","TriggerPad3Clusters"] - input_pass_name: str = "" + pass_name: str = "" output_collection: str = "clusters.txt" + verbosity: int = 0 @processor("trigscint::PatternLUTMaker", "TrigScint") @@ -524,3 +528,4 @@ class PatternLUTMaker(Processor) : input_pass_name: str = "" output_collection: str ="LUT.txt" lut_threshold: float = 0.0008 + verbosity: int = 0 diff --git a/TrigScint/src/TrigScint/ClusterTripletMaker.cxx b/TrigScint/src/TrigScint/ClusterTripletMaker.cxx index 4a1d42d238..2c2dffc514 100644 --- a/TrigScint/src/TrigScint/ClusterTripletMaker.cxx +++ b/TrigScint/src/TrigScint/ClusterTripletMaker.cxx @@ -4,10 +4,14 @@ namespace trigscint { +ClusterTripletMaker::ClusterTripletMaker(const std::string& name, + framework::Process& process) + : Analyzer(name, process) {} + void ClusterTripletMaker::configure(framework::config::Parameters& ps) { pass_name_ = ps.get("pass_name"); cluster_input_collections_ = ps.get>("cluster_input_collections"); - output_file_ = ps.get("output_file"); + output_collection_ = ps.get("output_collection"); verbose_ = ps.get("verbosity"); if (verbose_) { @@ -20,7 +24,7 @@ void ClusterTripletMaker::configure(framework::config::Parameters& ps) { << "\nInput collection 3: " << cluster_input_collections_.at(2) << "\nPass name: " << pass_name_ - << "\nOutput file: " << output_file_ + << "\nOutput file: " << output_collection_ << "\nVerbosity: " << verbose_; } @@ -28,10 +32,10 @@ void ClusterTripletMaker::configure(framework::config::Parameters& ps) { } void ClusterTripletMaker::onProcessStart() { - output_stream_.open(output_file_); + output_stream_.open(output_collection_); if (verbose_) { - ldmx_log(info) << "Opened output file " << output_file_; + ldmx_log(info) << "Opened output file " << output_collection_; } return; @@ -112,7 +116,7 @@ void ClusterTripletMaker::onProcessEnd() { if (verbose_) { ldmx_log(info) << "Closed output file " - << output_file_; + << output_collection_; } return; diff --git a/TrigScint/src/TrigScint/PatternLUTMaker.cxx b/TrigScint/src/TrigScint/PatternLUTMaker.cxx index 27d17cc538..07e0bc33a1 100644 --- a/TrigScint/src/TrigScint/PatternLUTMaker.cxx +++ b/TrigScint/src/TrigScint/PatternLUTMaker.cxx @@ -11,8 +11,8 @@ PatternLUTMaker::PatternLUTMaker(const std::string& name, : Analyzer(name, process) {} void PatternLUTMaker::configure(framework::config::Parameters& ps) { - input_file_ = ps.get("input_file"); - output_file_ = ps.get("output_file"); + input_collection_ = ps.get("input_collection"); + output_collection_ = ps.get("output_collection"); lut_threshold_ = ps.get("lut_threshold"); verbose_ = ps.get("verbosity"); @@ -28,8 +28,8 @@ void PatternLUTMaker::configure(framework::config::Parameters& ps) { if (verbose_) { ldmx_log(info) << "In PatternLUTMaker: configure done!" << std::endl; - ldmx_log(info) << "Got parameters: \nInput file: " << input_file_ - << "\nOutput file: " << output_file_ + ldmx_log(info) << "Got parameters: \nInput file: " << input_collection_ + << "\nOutput file: " << output_collection_ << "\nLUT threshold: " << lut_threshold_ << "\nVerbosity: " << verbose_; } @@ -37,8 +37,8 @@ void PatternLUTMaker::configure(framework::config::Parameters& ps) { } void PatternLUTMaker::onProcessStart() { - infile.open(input_file_); - outfile.open(output_file_); + infile.open(input_collection_); + outfile.open(output_collection_); return; } From 7e96e3936f20c1cf432f9040acb2da8dd27108d6 Mon Sep 17 00:00:00 2001 From: kvarnla Date: Sat, 6 Jun 2026 20:59:25 +0200 Subject: [PATCH 18/20] Update ClusterTripletMaker.h --- TrigScint/include/TrigScint/ClusterTripletMaker.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/TrigScint/include/TrigScint/ClusterTripletMaker.h b/TrigScint/include/TrigScint/ClusterTripletMaker.h index b06d6753bc..c3b9e09485 100644 --- a/TrigScint/include/TrigScint/ClusterTripletMaker.h +++ b/TrigScint/include/TrigScint/ClusterTripletMaker.h @@ -1,3 +1,9 @@ +/** + * @file ClusterTripletMaker.h + * @brief Writes cluster combinations to text file + * @author Lucia Kvarnstrom, Lund University + */ + #ifndef TRIGSCINT_CLUSTERTRIPLETMAKER_H #define TRIGSCINT_CLUSTERTRIPLETMAKER_H From f3b5745dcd9ac250c373bf32b679d91621915433 Mon Sep 17 00:00:00 2001 From: kvarnla Date: Sat, 6 Jun 2026 21:00:14 +0200 Subject: [PATCH 19/20] Update PatternLUTMaker.h --- TrigScint/include/TrigScint/PatternLUTMaker.h | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/TrigScint/include/TrigScint/PatternLUTMaker.h b/TrigScint/include/TrigScint/PatternLUTMaker.h index 31f64f91a6..622f3d0473 100644 --- a/TrigScint/include/TrigScint/PatternLUTMaker.h +++ b/TrigScint/include/TrigScint/PatternLUTMaker.h @@ -1,8 +1,6 @@ /** * @file PatternLUTMaker.h - * @brief Outputs a LUT based on frequency of track propagation patterns - * which can then be used as input in TrigScintTrackProducer for - * LUT-method track-making. + * @brief Writes LUT based on frequency of track propagation patterns for LUT-based TS tracking. * @author Lucia Kvarnstrom, Lund University */ From c2775592a49c7a0d70535b097447d512d71f43bf Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Sat, 6 Jun 2026 19:44:51 +0000 Subject: [PATCH 20/20] Apply clang-format and clang-tidy --- .../include/TrigScint/ClusterTripletMaker.h | 4 +- TrigScint/include/TrigScint/PatternLUTMaker.h | 31 +- .../TrigScint/TrigScintTrackProducer.h | 31 +- .../src/TrigScint/ClusterTripletMaker.cxx | 103 +++--- TrigScint/src/TrigScint/PatternLUTMaker.cxx | 62 ++-- .../src/TrigScint/TrigScintTrackProducer.cxx | 304 ++++++++++-------- 6 files changed, 272 insertions(+), 263 deletions(-) diff --git a/TrigScint/include/TrigScint/ClusterTripletMaker.h b/TrigScint/include/TrigScint/ClusterTripletMaker.h index c3b9e09485..f8cda6dda5 100644 --- a/TrigScint/include/TrigScint/ClusterTripletMaker.h +++ b/TrigScint/include/TrigScint/ClusterTripletMaker.h @@ -13,12 +13,12 @@ #include "Framework/EventProcessor.h" // TrigScint -#include "TrigScint/Event/TrigScintCluster.h" - #include #include #include +#include "TrigScint/Event/TrigScintCluster.h" + namespace trigscint { /** diff --git a/TrigScint/include/TrigScint/PatternLUTMaker.h b/TrigScint/include/TrigScint/PatternLUTMaker.h index 622f3d0473..49d780776c 100644 --- a/TrigScint/include/TrigScint/PatternLUTMaker.h +++ b/TrigScint/include/TrigScint/PatternLUTMaker.h @@ -1,27 +1,28 @@ /** * @file PatternLUTMaker.h - * @brief Writes LUT based on frequency of track propagation patterns for LUT-based TS tracking. + * @brief Writes LUT based on frequency of track propagation patterns for + * LUT-based TS tracking. * @author Lucia Kvarnstrom, Lund University */ #ifndef TRIGSCINT_PATTERNLUTMAKER_H #define TRIGSCINT_PATTERNLUTMAKER_H -#include "Framework/Configure/Parameters.h" -#include "Framework/EventProcessor.h" - #include #include #include #include #include +#include "Framework/Configure/Parameters.h" +#include "Framework/EventProcessor.h" + namespace trigscint { class PatternLUTMaker : public framework::Analyzer { public: struct Line { - int event; - float p1, p2, p3; + int event_; + float p1_, p2_, p3_; }; PatternLUTMaker(const std::string& name, framework::Process& process); @@ -39,23 +40,23 @@ class PatternLUTMaker : public framework::Analyzer { private: // verbosity int verbose_{0}; - + // input text file of clusters std::string input_collection_; - - //output LUT file name + + // output LUT file name std::string output_collection_; - - //minimum frequency of a specific pattern to be written to the LUT + + // minimum frequency of a specific pattern to be written to the LUT double lut_threshold_{0.0008}; - std::ifstream infile; - std::ofstream outfile; + std::ifstream infile_; + std::ofstream outfile_; // to group cluster combinations by track pattern - std::map, std::vector> groups; + std::map, std::vector> groups_; - int totalLines{0}; + int total_lines_{0}; }; } // namespace trigscint diff --git a/TrigScint/include/TrigScint/TrigScintTrackProducer.h b/TrigScint/include/TrigScint/TrigScintTrackProducer.h index 390a9ebfc0..d82e09b972 100644 --- a/TrigScint/include/TrigScint/TrigScintTrackProducer.h +++ b/TrigScint/include/TrigScint/TrigScintTrackProducer.h @@ -3,6 +3,8 @@ #define TRIGSCINT_TRIGSCINTTRACKPRODUCER_H // LDMX Framework +#include + #include "Framework/Configure/Parameters.h" // Needed to import parameters from configuration file #include "Framework/Event.h" #include "Framework/EventProcessor.h" //Needed to declare processor @@ -10,8 +12,6 @@ #include "TrigScint/Event/TrigScintCluster.h" #include "TrigScint/Event/TrigScintTrack.h" -#include - namespace trigscint { /** @@ -46,8 +46,8 @@ class TrigScintTrackProducer : public framework::Producer { // in the next pad tolerated to form a track double max_delta_{0.}; - double max_delta_vert{0.}; - double bar_length_y_{30.}; + double max_delta_vert_{0.}; + double bar_length_y_{30.}; // producer specific verbosity int verbose_{0}; @@ -66,7 +66,7 @@ class TrigScintTrackProducer : public framework::Producer { // allow forming tracks without match in the last collection bool skip_last_{false}; - + // do tracking using LUT method instead of with max_delta bool lut_tracking_{false}; @@ -92,23 +92,22 @@ class TrigScintTrackProducer : public framework::Producer { // track residual in units of channel nb (will not be content weighted) // float residual_{0.}; - - struct LUTKey { - float p1, p2, p3; - bool operator==(const LUTKey& other) const { - return p1 == other.p1 && p2 == other.p2 && p3 == other.p3; + struct LUTKey { + float p1_, p2_, p3_; + + bool operator==(const LUTKey &other) const { + return p1_ == other.p1_ && p2_ == other.p2_ && p3_ == other.p3_; } }; - + struct LUTKeyHash { - size_t operator()(const LUTKey& k) const { - return std::hash()(k.p1) ^ - (std::hash()(k.p2) << 1) ^ - (std::hash()(k.p3) << 2); + size_t operator()(const LUTKey &k) const { + return std::hash()(k.p1_) ^ (std::hash()(k.p2_) << 1) ^ + (std::hash()(k.p3_) << 2); } }; - + std::unordered_set lut_; float bar_width_y_{3.}; // mm diff --git a/TrigScint/src/TrigScint/ClusterTripletMaker.cxx b/TrigScint/src/TrigScint/ClusterTripletMaker.cxx index 2c2dffc514..02d25f8c5b 100644 --- a/TrigScint/src/TrigScint/ClusterTripletMaker.cxx +++ b/TrigScint/src/TrigScint/ClusterTripletMaker.cxx @@ -5,27 +5,27 @@ namespace trigscint { ClusterTripletMaker::ClusterTripletMaker(const std::string& name, - framework::Process& process) + framework::Process& process) : Analyzer(name, process) {} void ClusterTripletMaker::configure(framework::config::Parameters& ps) { - pass_name_ = ps.get("pass_name"); - cluster_input_collections_ = ps.get>("cluster_input_collections"); - output_collection_ = ps.get("output_collection"); - verbose_ = ps.get("verbosity"); + pass_name_ = ps.get("pass_name"); + cluster_input_collections_ = + ps.get>("cluster_input_collections"); + output_collection_ = ps.get("output_collection"); + verbose_ = ps.get("verbosity"); if (verbose_) { - ldmx_log(info) - << "In ClusterTripletMaker: configure done!" - << "\nInput collection 1: " - << cluster_input_collections_.at(0) - << "\nInput collection 2: " - << cluster_input_collections_.at(1) - << "\nInput collection 3: " - << cluster_input_collections_.at(2) - << "\nPass name: " << pass_name_ - << "\nOutput file: " << output_collection_ - << "\nVerbosity: " << verbose_; + ldmx_log(info) << "In ClusterTripletMaker: configure done!" + << "\nInput collection 1: " + << cluster_input_collections_.at(0) + << "\nInput collection 2: " + << cluster_input_collections_.at(1) + << "\nInput collection 3: " + << cluster_input_collections_.at(2) + << "\nPass name: " << pass_name_ + << "\nOutput file: " << output_collection_ + << "\nVerbosity: " << verbose_; } return; @@ -41,67 +41,58 @@ void ClusterTripletMaker::onProcessStart() { return; } -void ClusterTripletMaker::analyze( - const framework::Event& event) { +void ClusterTripletMaker::analyze(const framework::Event& event) { const auto clusters_pad1{event.getCollection( - cluster_input_collections_.at(0), pass_name_)}; + cluster_input_collections_.at(0), pass_name_)}; const auto clusters_pad2{event.getCollection( - cluster_input_collections_.at(1), pass_name_)}; + cluster_input_collections_.at(1), pass_name_)}; const auto clusters_pad3{event.getCollection( - cluster_input_collections_.at(2), pass_name_)}; + cluster_input_collections_.at(2), pass_name_)}; - const size_t num_clusters =std::min( - {clusters_pad1.size(), clusters_pad2.size(), clusters_pad3.size()}); + const size_t num_clusters = std::min( + {clusters_pad1.size(), clusters_pad2.size(), clusters_pad3.size()}); if (verbose_ > 1) { - ldmx_log(debug) - << "Event " << event.getEventNumber() - << " has " << clusters_pad1.size() - << ", " << clusters_pad2.size() << ", " - << clusters_pad3.size() - << " clusters in pads 1, 2, and 3."; + ldmx_log(debug) << "Event " << event.getEventNumber() << " has " + << clusters_pad1.size() << ", " << clusters_pad2.size() + << ", " << clusters_pad3.size() + << " clusters in pads 1, 2, and 3."; } - for (size_t cluster_index = 0; cluster_index < num_clusters; ++cluster_index) { + for (size_t cluster_index = 0; cluster_index < num_clusters; + ++cluster_index) { const int event_number = event.getEventNumber(); const float pad1_centroid = clusters_pad1.at(cluster_index).getCentroid(); const float pad2_centroid = clusters_pad2.at(cluster_index).getCentroid(); const float pad3_centroid = clusters_pad3.at(cluster_index).getCentroid(); - output_stream_ - << event_number << " " - << pad1_centroid << " " - << pad2_centroid << " " - << pad3_centroid << " " - << "\n"; + output_stream_ << event_number << " " << pad1_centroid << " " + << pad2_centroid << " " << pad3_centroid << " " << "\n"; } if (verbose_) { - for (size_t cluster_index = 1; - cluster_index < clusters_pad1.size(); ++cluster_index) { - ldmx_log(info) - << "Event " << event.getEventNumber() - << ", extra cluster in pad 1: " - << clusters_pad1.at(cluster_index).getCentroid(); + for (size_t cluster_index = 1; cluster_index < clusters_pad1.size(); + ++cluster_index) { + ldmx_log(info) << "Event " << event.getEventNumber() + << ", extra cluster in pad 1: " + << clusters_pad1.at(cluster_index).getCentroid(); } - for (size_t cluster_index = 1; - cluster_index < clusters_pad2.size(); ++cluster_index) { - ldmx_log(info) - << "Event " << event.getEventNumber() - << ", extra cluster in pad 2: " - << clusters_pad2.at(cluster_index).getCentroid(); + for (size_t cluster_index = 1; cluster_index < clusters_pad2.size(); + ++cluster_index) { + ldmx_log(info) << "Event " << event.getEventNumber() + << ", extra cluster in pad 2: " + << clusters_pad2.at(cluster_index).getCentroid(); } - for (size_t cluster_index = 1; - cluster_index < clusters_pad3.size(); ++cluster_index) { - ldmx_log(info) - << "Event " << event.getEventNumber() - << ", extra cluster in pad 3: " - << clusters_pad3.at(cluster_index).getCentroid(); + for (size_t cluster_index = 1; cluster_index < clusters_pad3.size(); + ++cluster_index) { + ldmx_log(info) << "Event " << event.getEventNumber() + << ", extra cluster in pad 3: " + << clusters_pad3.at(cluster_index).getCentroid(); } } @@ -114,9 +105,7 @@ void ClusterTripletMaker::onProcessEnd() { } if (verbose_) { - ldmx_log(info) - << "Closed output file " - << output_collection_; + ldmx_log(info) << "Closed output file " << output_collection_; } return; diff --git a/TrigScint/src/TrigScint/PatternLUTMaker.cxx b/TrigScint/src/TrigScint/PatternLUTMaker.cxx index 07e0bc33a1..c68fd7e8a1 100644 --- a/TrigScint/src/TrigScint/PatternLUTMaker.cxx +++ b/TrigScint/src/TrigScint/PatternLUTMaker.cxx @@ -1,13 +1,13 @@ -//Analyzer file which writes a LUT text file from input cluster text file (made in -//ClusterTripletMaker) containing only those cluster combination entries which -//occur with a frequency above the lut_threshold parameter +// Analyzer file which writes a LUT text file from input cluster text file (made +// in ClusterTripletMaker) containing only those cluster combination entries +// which occur with a frequency above the lut_threshold parameter #include "TrigScint/PatternLUTMaker.h" namespace trigscint { PatternLUTMaker::PatternLUTMaker(const std::string& name, - framework::Process& process) + framework::Process& process) : Analyzer(name, process) {} void PatternLUTMaker::configure(framework::config::Parameters& ps) { @@ -15,22 +15,23 @@ void PatternLUTMaker::configure(framework::config::Parameters& ps) { output_collection_ = ps.get("output_collection"); lut_threshold_ = ps.get("lut_threshold"); verbose_ = ps.get("verbosity"); - - //The LUT threshold is the minimum percentage of times that a track - //must appear in a pool of events to be written to the LUT. - //Example --- horizontal tracks whose pad 1-2 and pad 2-3 delta values - // are both 0 make up ~80% of tracks out of 10.000 events, - // while a track with deltas (+22,-22), meaning pad 1 cluster in - // e.g. bar 4, pad 2 cluster in bar 26, and pad 3 cluster in bar 4 - // only appears once in 10.000 events (0.01%). The straight - // tracks are written to the LUT and the single "anomaly" is not. - //Primary motivation here is for use in the case of an unknown TS misalignment. - + + // The LUT threshold is the minimum percentage of times that a track + // must appear in a pool of events to be written to the LUT. + // Example --- horizontal tracks whose pad 1-2 and pad 2-3 delta values + // are both 0 make up ~80% of tracks out of 10.000 events, + // while a track with deltas (+22,-22), meaning pad 1 cluster in + // e.g. bar 4, pad 2 cluster in bar 26, and pad 3 cluster in bar 4 + // only appears once in 10.000 events (0.01%). The straight + // tracks are written to the LUT and the single "anomaly" is not. + // Primary motivation here is for use in the case of an unknown TS + // misalignment. + if (verbose_) { ldmx_log(info) << "In PatternLUTMaker: configure done!" << std::endl; ldmx_log(info) << "Got parameters: \nInput file: " << input_collection_ << "\nOutput file: " << output_collection_ - << "\nLUT threshold: " << lut_threshold_ + << "\nLUT threshold: " << lut_threshold_ << "\nVerbosity: " << verbose_; } return; @@ -64,13 +65,11 @@ void PatternLUTMaker::onProcessEnd() { int combs = 0; int tracks = 0; - if(verbose_){ - ldmx_log(info) << "Total number of track candidates: " - << totalLines << "\n" - << "Number of track candidate types: " - << groups.size() << "\n" - << "LUT Threshold: " - << lut_threshold_ * 100 << "%\n"; + if (verbose_) { + ldmx_log(info) << "Total number of track candidates: " << totalLines << "\n" + << "Number of track candidate types: " << groups.size() + << "\n" + << "LUT Threshold: " << lut_threshold_ * 100 << "%\n"; } for (auto& g : groups) { @@ -81,9 +80,8 @@ void PatternLUTMaker::onProcessEnd() { double frac = static_cast(count) / totalLines; if (verbose_) { - ldmx_log(info) << "(" << p12 << "," << p23 << ") appears " - << count << " times, (" << frac * 100 << " %)" - << "\n"; + ldmx_log(info) << "(" << p12 << "," << p23 << ") appears " << count + << " times, (" << frac * 100 << " %)" << "\n"; } } @@ -92,23 +90,21 @@ void PatternLUTMaker::onProcessEnd() { double frac = static_cast(count) / totalLines; - if (frac > lut_threshold_) { //write to outfile if over threshold + if (frac > lut_threshold_) { // write to outfile if over threshold combs++; for (auto& line : g.second) { tracks++; - outfile << line.p1 << " " - << line.p2 << " " - << line.p3 << "\n"; + outfile << line.p1 << " " << line.p2 << " " << line.p3 << "\n"; } } } if (verbose_) { - ldmx_log(info) << "\nLUT textfile written." - << "\n" << combs << " combinations (" - << tracks << " tracks) written to LUT.\n"; + ldmx_log(info) << "\nLUT textfile written." << "\n" + << combs << " combinations (" << tracks + << " tracks) written to LUT.\n"; } return; diff --git a/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx b/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx index 799c825e8e..0f3f7f7104 100644 --- a/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx +++ b/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx @@ -1,18 +1,19 @@ #include "TrigScint/TrigScintTrackProducer.h" -//Ricardo,05-26 +// Ricardo,05-26 +#include +#include #include // std::next #include -#include -#include namespace trigscint { void TrigScintTrackProducer::configure(framework::config::Parameters &ps) { max_delta_ = ps.get( - "delta_max"); // max distance to consider adding in a cluster to track - max_delta_vert= ps.get("delta_vert_max"); //max distance between pad 1/2 and 3 along the x axis - //to consider make a track using the vertical bars + "delta_max"); // max distance to consider adding in a cluster to track + max_delta_vert = ps.get( + "delta_vert_max"); // max distance between pad 1/2 and 3 along the x axis + // to consider make a track using the vertical bars seeding_collection_ = ps.get( "seeding_collection"); // probably tagger pad, "TriggerPadTagClusters" input_collections_ = ps.get>( @@ -29,9 +30,11 @@ void TrigScintTrackProducer::configure(framework::config::Parameters &ps) { bar_width_x_ = ps.get("vertical_bar_width"); bar_gap_x_ = ps.get("vertical_bar_gap"); skip_last_ = ps.get("allow_skip_last_collection"); - bar_length_y_=ps.get("horizontal_bar_length"); //bar length of the horizontal bars - lut_tracking_ = ps.get("lut_tracking"); - std::string lut_file_ = ps.get("lut_file"); //from PatternLUTMaker + bar_length_y_ = ps.get( + "horizontal_bar_length"); // bar length of the horizontal bars + lut_tracking_ = ps.get("lut_tracking"); + std::string lut_file = + ps.get("lut_file"); // from PatternLUTMaker // TO DO: allow any number of input collections @@ -45,7 +48,7 @@ void TrigScintTrackProducer::configure(framework::config::Parameters &ps) { << "\nAllow tracks with no hit in last collection: " << skip_last_ << "\nUsing LUT Tracking Method: " << lut_tracking_ - << "\nIf using LUT Method, LUT from: " << lut_file_ + << "\nIf using LUT Method, LUT from: " << lut_file << "\nVertical bar start index: " << vert_bar_start_idx_ << "\nNumber of horizontal bars: " << n_bars_y_ << "\nHorizontal bar width: " << bar_width_y_ @@ -64,23 +67,21 @@ void TrigScintTrackProducer::configure(framework::config::Parameters &ps) { x_conv_factor_ = bar_width_x_ + bar_gap_x_; // half width of pad x_start_ = -(n_bars_x_ * (bar_width_x_ + bar_gap_x_) - bar_gap_x_) / 2.; - - if (lut_tracking_) { - - std::ifstream file(lut_file_); + + if (lut_tracking_) { + std::ifstream file(lut_file); float a, b, c; - + while (file >> a >> b >> c) { float p1 = a; float p2 = b; float p3 = c; lut_.insert({p1, p2, p3}); } - + ldmx_log(info) << "Loaded LUT with size: " << lut_.size(); - } - + return; } @@ -174,102 +175,106 @@ void TrigScintTrackProducer::produce(framework::Event &event) { // reset for each seed // bool madeTrack = false; - - if (lut_tracking_) { //if using LUT method + + if (lut_tracking_) { // if using LUT method for (const auto &cluster1 : clusters_pad1) { for (const auto &cluster2 : clusters_pad2) { - float seed_bin = seed.getCentroid(); float pad1_bin = cluster1.getCentroid(); float pad2_bin = cluster2.getCentroid(); - - LUTKey key{seed_bin, pad1_bin, pad2_bin}; //LUTKey defined in header file - + + LUTKey key{seed_bin, pad1_bin, + pad2_bin}; // LUTKey defined in header file + if (lut_.find(key) != lut_.end()) { std::vector three_cluster_vec = { - seed, cluster1, cluster2}; - + seed, cluster1, cluster2}; + ldmx::TrigScintTrack track = makeTrack(three_cluster_vec); track_candidates.push_back(track); } } - } - - } - - else { - - for (const auto &cluster1 : clusters_pad1) { - if (verbose_ > 1) { - ldmx_log(debug) << "\tGot pad1 cluster with centroid " - << cluster1.getCentroid(); } - if ((fabs(cluster1.getCentroid() - centroid) < max_delta_ && - centroid < vert_bar_start_idx_) || - (centroid >= vert_bar_start_idx_ && cluster1.getCentroid() >= vert_bar_start_idx_ && - seed.getCentroidX() == cluster1.getCentroidX())) { - // use geometry y overlap scheme to see if this is really a match in x - // should be done in a map - - if (centroid >= vert_bar_start_idx_ && - seed.getCentroidY() < cluster1.getCentroidY()) { - // impossible combination - if (verbose_ > 1) { - ldmx_log(debug) << "\tSkipping impossible x cluster combination " - "with y flags (tag up) (" - << seed.getCentroidY() << " " - << cluster1.getCentroidY() << ")"; - } - continue; - } - // else: first (possible) match! loop through next pad too + } + else { + for (const auto &cluster1 : clusters_pad1) { if (verbose_ > 1) { - ldmx_log(debug) << "\t\tIt is close enough!. Check pad2"; + ldmx_log(debug) << "\tGot pad1 cluster with centroid " + << cluster1.getCentroid(); } + if ((fabs(cluster1.getCentroid() - centroid) < max_delta_ && + centroid < vert_bar_start_idx_) || + (centroid >= vert_bar_start_idx_ && + cluster1.getCentroid() >= vert_bar_start_idx_ && + seed.getCentroidX() == cluster1.getCentroidX())) { + // use geometry y overlap scheme to see if this is really a match in + // x should be done in a map + + if (centroid >= vert_bar_start_idx_ && + seed.getCentroidY() < cluster1.getCentroidY()) { + // impossible combination + if (verbose_ > 1) { + ldmx_log(debug) + << "\tSkipping impossible x cluster combination " + "with y flags (tag up) (" + << seed.getCentroidY() << " " << cluster1.getCentroidY() + << ")"; + } + continue; + } - // try making third pad clusters an optional part of track - - std::vector cluster_vec = {seed, cluster1}; - - bool has_match_dn = false; + // else: first (possible) match! loop through next pad too - for (const auto &cluster2 : clusters_pad2) { if (verbose_ > 1) { - ldmx_log(debug) << "\tGot pad2 cluster with centroid " - << cluster2.getCentroid(); + ldmx_log(debug) << "\t\tIt is close enough!. Check pad2"; } - if ((fabs(cluster2.getCentroid() - centroid) < max_delta_ && - centroid < vert_bar_start_idx_) || - (centroid >= vert_bar_start_idx_ && cluster2.getCentroid() >= vert_bar_start_idx_ && - fabs(seed.getCentroidX() - cluster2.getCentroidX()) <= max_delta_vert)) { - // use geometry y overlap scheme to see if this is really a match - // in x - - if (centroid >= vert_bar_start_idx_ && - (seed.getCentroidY() < cluster2.getCentroidY() || - cluster1.getCentroidY() > - cluster2.getCentroidY())) { // impossible - if (verbose_ > 1) { - ldmx_log(debug) - << "\tSkipping impossible x cluster combination with y " - "flags (tag up dn) (" - << seed.getCentroidY() << " " << cluster1.getCentroidY() - << " " << cluster2.getCentroidY() << ")"; - } - continue; - } + // try making third pad clusters an optional part of track + + std::vector cluster_vec = {seed, cluster1}; - // first match! loop through next pad too + bool has_match_dn = false; + for (const auto &cluster2 : clusters_pad2) { if (verbose_ > 1) { - ldmx_log(debug) << "\t\tIt is close enough!. Make a track"; + ldmx_log(debug) << "\tGot pad2 cluster with centroid " + << cluster2.getCentroid(); } - // only make this vector now! this ensures against hanging - // clusters with indices from earlier in the loop + if ((fabs(cluster2.getCentroid() - centroid) < max_delta_ && + centroid < vert_bar_start_idx_) || + (centroid >= vert_bar_start_idx_ && + cluster2.getCentroid() >= vert_bar_start_idx_ && + fabs(seed.getCentroidX() - cluster2.getCentroidX()) <= + max_delta_vert)) { + // use geometry y overlap scheme to see if this is really a + // match + // in x + + if (centroid >= vert_bar_start_idx_ && + (seed.getCentroidY() < cluster2.getCentroidY() || + cluster1.getCentroidY() > + cluster2.getCentroidY())) { // impossible + if (verbose_ > 1) { + ldmx_log(debug) + << "\tSkipping impossible x cluster combination with y " + "flags (tag up dn) (" + << seed.getCentroidY() << " " << cluster1.getCentroidY() + << " " << cluster2.getCentroidY() << ")"; + } + continue; + } + + // first match! loop through next pad too + + if (verbose_ > 1) { + ldmx_log(debug) << "\t\tIt is close enough!. Make a track"; + } + + // only make this vector now! this ensures against hanging + // clusters with indices from earlier in the loop std::vector three_cluster_vec = { seed, cluster1, cluster2}; @@ -300,9 +305,8 @@ break; */ } // over clusters in pad1 - } - + // continue to next seed if 0 track candidates if (track_candidates.size() == 0) continue; @@ -381,10 +385,15 @@ break; // no need to start pulling constituents from tracks that are // ridiculously far apart if (((fabs(track.getCentroid() - next_track.getCentroid()) < - 3 * max_delta_) && (track.getCentroid()=vert_bar_start_idx_))) { - //and for the vertical bars, check if they are in the same quad and close enough + 3 * max_delta_) && + (track.getCentroid() < + vert_bar_start_idx_)) // for the horizontal bars + || ((fabs(track.getCentroidX() - next_track.getCentroidX()) < + 2 * max_delta_vert) && + (track.getCentroidY() == next_track.getCentroidY()) && + (track.getCentroid() >= vert_bar_start_idx_))) { + // and for the vertical bars, check if they are in the same quad and + // close enough std::vector consts_1 = track.getConstituents(); std::vector consts_2 = @@ -398,12 +407,15 @@ break; // let's do "if either cluster is shared" right now... but could also // have it settable to use a stricter cut: an AND if (((consts_1[1].getCentroid() == consts_2[1].getCentroid() || - ((consts_1.size() > 2) && (consts_2.size() > 2) && - (consts_1[2].getCentroid() == consts_2[2].getCentroid()))) && (track.getCentroid()=vert_bar_start_idx_) && ((consts_1[1].getCentroidX() == consts_2[1].getCentroidX()) || - (consts_1[2].getCentroidX() == consts_2[2].getCentroidX()) - || (consts_1[0].getCentroidX() == consts_2[0].getCentroidX())))) { //and vertical bars + ((consts_1.size() > 2) && (consts_2.size() > 2) && + (consts_1[2].getCentroid() == consts_2[2].getCentroid()))) && + (track.getCentroid() < vert_bar_start_idx_)) || + // horizontal bars + ((track.getCentroid() >= vert_bar_start_idx_) && + ((consts_1[1].getCentroidX() == consts_2[1].getCentroidX()) || + (consts_1[2].getCentroidX() == consts_2[2].getCentroidX()) || + (consts_1[0].getCentroidX() == + consts_2[0].getCentroidX())))) { // and vertical bars if (verbose_ > 1) { ldmx_log(debug) << "Found overlap! Tracks at index " << idx @@ -412,11 +424,19 @@ break; ldmx_log(trace) << tracks_.at(idx_comp); } - if (((fabs((tracks_.at(idx)).getResidualX() -(tracks_.at(idx_comp)).getResidualX()))< 0.01) //it should be equal - && (track.getCentroid()>=vert_bar_start_idx_)){ //specific case for the vertical bars - continue; //currently we can't do more here - } else if (((tracks_.at(idx)).getResidual() <(tracks_.at(idx_comp)).getResidual() && (track.getCentroid()=vert_bar_start_idx_))) { + if (((fabs((tracks_.at(idx)).getResidualX() - + (tracks_.at(idx_comp)).getResidualX())) < + 0.01) // it should be equal + && + (track.getCentroid() >= + vert_bar_start_idx_)) { // specific case for the vertical bars + continue; // currently we can't do more here + } else if (((tracks_.at(idx)).getResidual() < + (tracks_.at(idx_comp)).getResidual() && + (track.getCentroid() < vert_bar_start_idx_)) || + ((tracks_.at(idx)).getResidualX() < + (tracks_.at(idx_comp)).getResidualX() && + (track.getCentroid() >= vert_bar_start_idx_))) { // next track (lower index) is a worse choice, remove its flag for // keeping keep_indices.at(idx_comp) = 0; @@ -585,16 +605,14 @@ ldmx::TrigScintTrack TrigScintTrackProducer::makeTrack( ((clusters.at(i)).getCentroid() - centroid); residual = sqrt(residual / clusters.size()); - - float residual_x = 0; //only for the vertical bars - if (centroid>=vert_bar_start_idx_) { + float residual_x = 0; // only for the vertical bars + if (centroid >= vert_bar_start_idx_) { for (uint i = 0; i < clusters.size(); i++) residual_x += ((clusters.at(i)).getCentroidX() - centroid_x) * - ((clusters.at(i)).getCentroidX() - centroid_x); + ((clusters.at(i)).getCentroidX() - centroid_x); residual_x = sqrt(residual_x / clusters.size()); } - tr.setResidualX(residual_x); tr.setCentroid(centroid); tr.setCentroidX(centroid_x); @@ -638,25 +656,28 @@ void TrigScintTrackProducer::matchXYTracks( if (verbose_) ldmx_log(debug) << " -- In matchXYTracks found y track at " << trk.getCentroidY() << "; mapping to quad " - << (int)trk.getCentroidY() / (n_bars_y_/2) << " with trk index " - << trk_idx; + << (int)trk.getCentroidY() / (n_bars_y_ / 2) + << " with trk index " << trk_idx; // 2. order them... or map them to quadrants. note that there are 2 layers // so 2*n_bars_y_/4 channels per quadrant - y_quad_map.insert(std::make_pair((int)(trk.getCentroidY() / (n_bars_y_/2)), trk)); + y_quad_map.insert( + std::make_pair((int)(trk.getCentroidY() / (n_bars_y_ / 2)), trk)); y_track_map[trk] = trk_idx; y_idx_quad_map.insert( - std::make_pair((int)(trk.getCentroidY() / (n_bars_y_/2)), trk_idx)); + std::make_pair((int)(trk.getCentroidY() / (n_bars_y_ / 2)), trk_idx)); } else { // 3. get the remaining tracks (from vertical bars) and map them // (back) to (middle of) quadrants - x_quad_map.insert(std::make_pair((int)(trk.getCentroidY() / (n_bars_y_/2)), trk)); + x_quad_map.insert( + std::make_pair((int)(trk.getCentroidY() / (n_bars_y_ / 2)), trk)); x_track_map[trk] = trk_idx; x_idx_quad_map.insert( - std::make_pair((int)(trk.getCentroidY() / (n_bars_y_/2)), trk_idx)); + std::make_pair((int)(trk.getCentroidY() / (n_bars_y_ / 2)), trk_idx)); if (verbose_) ldmx_log(debug) << " -- In matchXYTracks found x track at (x,y) = (" << trk.getCentroidX() << ", " << trk.getCentroidY() - << "); mapping to quad " << (int)trk.getCentroidY() / (n_bars_y_/2) + << "); mapping to quad " + << (int)trk.getCentroidY() / (n_bars_y_ / 2) << " with trk index " << trk_idx; } } @@ -672,8 +693,8 @@ void TrigScintTrackProducer::matchXYTracks( float x0 = 0; // this should be half the pad... could also set // it to full beam spot width - float sx0 = fabs(x_start_); - float sx0_vert=fabs(bar_length_y_/2); // When there are no hits + float sx0 = fabs(x_start_); + float sx0_vert = fabs(bar_length_y_ / 2); // When there are no hits // along the vertical bars // y_start_ is half the pad, so this should be half a quadrant @@ -687,7 +708,7 @@ void TrigScintTrackProducer::matchXYTracks( float y{-9999.}, sy{-9999.}, x{-9999.}, x1{-9999.}, x2{-9999.}, sx1{-9999.}, sx2{-9999.}, y1{-9999.}, y2{-9999.}, sy1{-9999.}, sy2{-9999.}; // quad midpoint: - float y0 = (((*yitr).first * 8)*y_conv_factor_ )+y_start_+ sy0; + float y0 = (((*yitr).first * 8) * y_conv_factor_) + y_start_ + sy0; float sx = 1. / 2 * x_conv_factor_; // rely on x precision being one single bar // width; always used unless x is undeterminable @@ -731,7 +752,8 @@ void TrigScintTrackProducer::matchXYTracks( sx1 = x_conv_factor_ / 2.; // 1 bar width sx2 = sx1; x = (x1 + x2) / 2.; - sx = fabs(x1 - x2) / 2; // rely on x precision being one single pad width + sx = fabs(x1 - x2) / + 2; // rely on x precision being one single pad width if (verbose_) ldmx_log(debug) << "\t\t -- 2 x in quad: setting y track x " "coordinate to midpoint"; @@ -742,12 +764,13 @@ void TrigScintTrackProducer::matchXYTracks( x = x0; sx = sx0; if (verbose_) - ldmx_log(debug) << "\t\t\t currently no x info assigned in ambiguous case of " - << n_xin_quad - << "vertical bar track candidates in quad " << (*yitr).first - << "; will set x to middle of pad, pad half-width as " - "precision: set (x, sx)=(" - << x << ", " << sx << ")"; + ldmx_log(debug) + << "\t\t\t currently no x info assigned in ambiguous case of " + << n_xin_quad << "vertical bar track candidates in quad " + << (*yitr).first + << "; will set x to middle of pad, pad half-width as " + "precision: set (x, sx)=(" + << x << ", " << sx << ")"; } // 3 x tracks in quadrant // ok! over y: @@ -764,7 +787,7 @@ void TrigScintTrackProducer::matchXYTracks( // 4. every quadrant which just has one of each --> done ; // b) set the sx, sy of the x track now, using the residuals from the // other b1) special case: no x tracks; then x, sx have been set above - if (n_xin_quad==1){ + if (n_xin_quad == 1) { auto xidx = x_idx_quad_map.find((*yitr).first); tracks.at((*xidx).second).setPosition(x, y); tracks.at((*xidx).second).setSigmaXY(sx, sy); @@ -773,10 +796,10 @@ void TrigScintTrackProducer::matchXYTracks( ldmx_log(debug) << "\t\t\t in quad " << (*yitr).first << ", set (x, y) = (" << x << ", " << y << ") and (sx, sy) = " << sx << ", " << sy << ")"; - auto yidx = y_idx_quad_map.find((*yitr).first); - tracks.at((*yidx).second).setPosition(x, y); - tracks.at((*yidx).second).setSigmaXY(sx, sy); - continue; + auto yidx = y_idx_quad_map.find((*yitr).first); + tracks.at((*yidx).second).setPosition(x, y); + tracks.at((*yidx).second).setSigmaXY(sx, sy); + continue; } } // 1 y, 0 or 1 or 3+ x track in quadrant @@ -787,7 +810,7 @@ void TrigScintTrackProducer::matchXYTracks( if (n_yin_quad == 2) { // let's start here and see if we can do >= 2 later // here one could do sth to avoid checking the other y track again in the - // outermost loop over y + // outermost loop over y auto yitr1 = y_quad_map.lower_bound((*yitr).first); auto yitr2 = y_quad_map.upper_bound((*yitr).first); yitr2--; // back up once @@ -798,17 +821,18 @@ void TrigScintTrackProducer::matchXYTracks( if (sy1 == 0) sy1 = 1. / 2 * y_conv_factor_; if (sy2 == 0) sy2 = 1. / 2 * y_conv_factor_; y = (y1 + y2) / 2.; - sy = fabs(y1 - y2) / 2 ; + sy = fabs(y1 - y2) / 2; if (verbose_) ldmx_log(debug) << "\t\t -- 2 y in quad: setting x track y coordinate to midpoint"; } // 2y in quad - - if ((n_xin_quad == 0 || n_xin_quad >= 3) && (n_yin_quad == 2)) { //not using the X tracks for now for >=3 - if (n_xin_quad == 0){ + + if ((n_xin_quad == 0 || n_xin_quad >= 3) && + (n_yin_quad == 2)) { // not using the X tracks for now for >=3 + if (n_xin_quad == 0) { if (verbose_) - ldmx_log(debug) - << "\t\t -- No x tracks but 2 y tracks in quad: unusual behaviour"; + ldmx_log(debug) << "\t\t -- No x tracks but 2 y tracks in quad: " + "unusual behaviour"; } auto yidx1 = y_idx_quad_map.lower_bound((*yitr).first); auto yidx2 = y_idx_quad_map.upper_bound((*yitr).first); @@ -819,7 +843,7 @@ void TrigScintTrackProducer::matchXYTracks( tracks.at((*yidx2).second).setSigmaXY(sx, sy2); continue; } - + if (n_yin_quad == 1 && n_xin_quad == 2) { // don't think we want to experiment with discerning // three overlapping tracks, so not >= 2