diff --git a/Tools/include/Tools/ClusterTripletMaker.h b/Tools/include/Tools/ClusterTripletMaker.h new file mode 100644 index 0000000000..f8cda6dda5 --- /dev/null +++ b/Tools/include/Tools/ClusterTripletMaker.h @@ -0,0 +1,62 @@ +/** + * @file ClusterTripletMaker.h + * @brief Writes cluster combinations to text file + * @author Lucia Kvarnstrom, Lund University + */ + +#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 +#include +#include + +#include "TrigScint/Event/TrigScintCluster.h" + +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_collection_{"clusters.txt"}; + + // output stream + std::ofstream output_stream_; +}; + +} // namespace trigscint + +#endif // TRIGSCINT_CLUSTERTRIPLETMAKER_H diff --git a/Tools/include/Tools/PatternLUTMaker.h b/Tools/include/Tools/PatternLUTMaker.h new file mode 100644 index 0000000000..49d780776c --- /dev/null +++ b/Tools/include/Tools/PatternLUTMaker.h @@ -0,0 +1,64 @@ +/** + * @file PatternLUTMaker.h + * @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 +#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_; + }; + + 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_collection_; + + // output LUT file name + std::string output_collection_; + + // 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 total_lines_{0}; +}; + +} // namespace trigscint + +#endif /* TRIGSCINT_PATTERNLUTMAKER_H */ diff --git a/Tools/python/lut_making.py b/Tools/python/lut_making.py new file mode 100644 index 0000000000..1c9125df53 --- /dev/null +++ b/Tools/python/lut_making.py @@ -0,0 +1,23 @@ +"""Configuration for Cluster Triplet text file Maker and Pattern LUT Maker""" + +from LDMX.Framework import parameter_set + +@parameter_set +class ClusterTripletMaker : + """Configuration for cluster text file maker for Trigger Scintillators""" + cluster_input_collections: list[str] = ["TriggerPad1Clusters", + "TriggerPad2Clusters","TriggerPad3Clusters"] + pass_name: str = "" + output_collection: str = "clusters.txt" + verbosity: int = 0 + + +@parameter_set +class PatternLUTMaker : + """Configuration for track-pattern LUT-writing analyzer for Trigger Scintillators""" + + input_collection: str ="clusters.txt" + input_pass_name: str = "" + output_collection: str ="LUT.txt" + lut_threshold: float = 0.0008 + verbosity: int = 0 \ No newline at end of file diff --git a/Tools/src/Tools/ClusterTripletMaker.cxx b/Tools/src/Tools/ClusterTripletMaker.cxx new file mode 100644 index 0000000000..02d25f8c5b --- /dev/null +++ b/Tools/src/Tools/ClusterTripletMaker.cxx @@ -0,0 +1,116 @@ +#include "TrigScint/ClusterTripletMaker.h" + +#include + +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_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_; + } + + return; +} + +void ClusterTripletMaker::onProcessStart() { + output_stream_.open(output_collection_); + + if (verbose_) { + ldmx_log(info) << "Opened output file " << output_collection_; + } + + 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_collection_; + } + + return; +} + +} // namespace trigscint + +DECLARE_ANALYZER(trigscint::ClusterTripletMaker); diff --git a/Tools/src/Tools/PatternLUTMaker.cxx b/Tools/src/Tools/PatternLUTMaker.cxx new file mode 100644 index 0000000000..c68fd7e8a1 --- /dev/null +++ b/Tools/src/Tools/PatternLUTMaker.cxx @@ -0,0 +1,115 @@ +// 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_collection_ = ps.get("input_collection"); + 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. + + 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_ + << "\nVerbosity: " << verbose_; + } + return; +} + +void PatternLUTMaker::onProcessStart() { + infile.open(input_collection_); + outfile.open(output_collection_); + 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/TrigScint/exampleConfigs/runClusterstxt.py b/TrigScint/exampleConfigs/runClusterstxt.py new file mode 100644 index 0000000000..778dc1157c --- /dev/null +++ b/TrigScint/exampleConfigs/runClusterstxt.py @@ -0,0 +1,49 @@ +#example config which makes digis, clusters, and writes clusters +#to .txt file using ClusterTripletMaker + +from LDMX.Framework import ldmxcfg + +p = ldmxcfg.Process('nontruthclusters') +p.input_files = ["SimSamples.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 + +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 + +from LDMX.TrigScint.trig_scint import ClusterTripletMaker + +triplets = ClusterTripletMaker("tripletmaker") +triplets.output_collection = "clusters.txt" + +p.sequence = [ + #*truth_hits, + *digis, + *clusters, + triplets + ] + + diff --git a/TrigScint/exampleConfigs/runLUTana.py b/TrigScint/exampleConfigs/runLUTana.py new file mode 100644 index 0000000000..468a4eab7c --- /dev/null +++ b/TrigScint/exampleConfigs/runLUTana.py @@ -0,0 +1,15 @@ +#config file to make LUT for tracking using PatternLUTMaker + +from LDMX.Framework import ldmxcfg + +p = ldmxcfg.Process("LUTmaker") +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_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 new file mode 100644 index 0000000000..68bab3cf6a --- /dev/null +++ b/TrigScint/exampleConfigs/runLUTtracking.py @@ -0,0 +1,27 @@ +#Performs tracking with TrigScintTrackProducer, with and without LUT method for comparison +from LDMX.Framework import ldmxcfg + +p = ldmxcfg.Process('tracking') +p.input_files = ["Clusters.root"] +p.output_files = ["Tracks.root"] + +from LDMX.TrigScint.trig_scint import TrigScintTrackProducer + +tstp_tracks = TrigScintTrackProducer("tstp") +lut_tracks = TrigScintTrackProducer("lut") + +tstp_tracks.delta_max = 1.0 +tstp_tracks.verbosity = 1 +tstp_tracks.output_collection = "tstpTracks" + +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 + +p.sequence = [tstp_tracks, + lut_tracks + ] diff --git a/TrigScint/include/TrigScint/TrigScintTrackProducer.h b/TrigScint/include/TrigScint/TrigScintTrackProducer.h index fd608d2086..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 @@ -44,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}; @@ -65,6 +67,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}; @@ -88,6 +93,23 @@ 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 float bar_width_x_{3.}; // mm diff --git a/TrigScint/python/trig_scint.py b/TrigScint/python/trig_scint.py index a5b5ed7bb9..f8967d3408 100644 --- a/TrigScint/python/trig_scint.py +++ b/TrigScint/python/trig_scint.py @@ -391,7 +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") @@ -506,3 +508,4 @@ class TestBeamClusterAnalyzer(Processor): input_hit_pass_name: str = "" start_sample: int = 2 dead_channels: list[int] = [8] + diff --git a/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx b/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx index fe074335a5..0f3f7f7104 100644 --- a/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx +++ b/TrigScint/src/TrigScint/TrigScintTrackProducer.cxx @@ -1,16 +1,19 @@ #include "TrigScint/TrigScintTrackProducer.h" -//Ricardo,05-26 +// Ricardo,05-26 +#include +#include #include // std::next -#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>( @@ -27,7 +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 lenght of the horizontal bars + 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 @@ -40,6 +47,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_ @@ -59,6 +68,20 @@ void TrigScintTrackProducer::configure(framework::config::Parameters &ps) { // 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; } @@ -153,107 +176,136 @@ 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_) || - (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() << ")"; + 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 { + 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 - std::vector three_cluster_vec = { - seed, cluster1, cluster2}; + 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; + } - /* - // 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); - } + // 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}; - } // if possible (x,)y match in pad1 /* - //same here - if (madeTrack) - break; - */ + // 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); + } - } // over clusters in pad1 + } // if possible (x,)y match in pad1 + /* +//same here +if (madeTrack) +break; +*/ + + } // over clusters in pad1 + } // continue to next seed if 0 track candidates if (track_candidates.size() == 0) continue; @@ -333,10 +385,15 @@ void TrigScintTrackProducer::produce(framework::Event &event) { // 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 = @@ -350,12 +407,15 @@ void TrigScintTrackProducer::produce(framework::Event &event) { // 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 @@ -364,11 +424,19 @@ void TrigScintTrackProducer::produce(framework::Event &event) { 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; @@ -537,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); @@ -590,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; } } @@ -624,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 @@ -639,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 @@ -683,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"; @@ -694,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: @@ -716,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); @@ -725,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 @@ -739,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 @@ -750,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); @@ -771,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