From 327aa485497ffe6dfc341f4bda742222e9d10912 Mon Sep 17 00:00:00 2001 From: Alexandra Moor Date: Fri, 14 Feb 2025 04:04:41 -0600 Subject: [PATCH 01/18] Updates to test crttpc matching with data --- sbndcode/CRT/CRTAna/CRTAnalysis_module.cc | 389 ++++++++++++++++-- sbndcode/CRT/CRTAna/run_crtana_data.fcl | 9 +- .../CRT/CRTEventDisplay/CRTEventDisplayAlg.cc | 100 +++++ .../CRT/CRTEventDisplay/CRTEventDisplayAlg.h | 4 + .../crteventdisplayalg_sbnd.fcl | 1 + .../crttpcmatchingproducers_sbnd.fcl | 12 +- .../CRT/CalibService/run_caloskim_data.fcl | 10 + 7 files changed, 488 insertions(+), 37 deletions(-) create mode 100644 sbndcode/Calibration/CRT/CalibService/run_caloskim_data.fcl diff --git a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc index 718ad3540..9d1f42c06 100644 --- a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc +++ b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc @@ -305,12 +305,33 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer { std::vector _tpc_sp_ts0; std::vector _tpc_sp_ts1; std::vector _tpc_sp_score; + std::vector _tpc_sp_x; + std::vector _tpc_sp_y; + std::vector _tpc_sp_z; + std::vector _tpc_sp_tag; std::vector _tpc_tr_matchable; std::vector _tpc_tr_matched; std::vector _tpc_tr_good_match; std::vector _tpc_tr_ts0; std::vector _tpc_tr_ts1; std::vector _tpc_tr_score; + std::vector _tpc_tr_start_x; + std::vector _tpc_tr_start_y; + std::vector _tpc_tr_start_z; + std::vector _tpc_tr_dir_x; + std::vector _tpc_tr_dir_y; + std::vector _tpc_tr_dir_z; + std::vector _tpc_tr_tagger1; + std::vector _tpc_tr_tagger2; + std::vector _tpc_tr_tagger3; + std::vector _tpc_diff_tagger_1; + std::vector _tpc_diff1_x; + std::vector _tpc_diff1_y; + std::vector _tpc_diff1_z; + std::vector _tpc_diff_tagger_2; + std::vector _tpc_diff2_x; + std::vector _tpc_diff2_y; + std::vector _tpc_diff2_z; std::vector _ptb_hlt_trigger; std::vector _ptb_hlt_timestamp; @@ -560,6 +581,35 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fTree->Branch("tpc_dir_z", "std::vector", &_tpc_dir_z); fTree->Branch("tpc_length", "std::vector", &_tpc_length); fTree->Branch("tpc_track_score", "std::vector", &_tpc_track_score); + fTree->Branch("tpc_sp_matched", "std::vector", &_tpc_sp_matched); + fTree->Branch("tpc_sp_ts0", "std::vector", &_tpc_sp_ts0); + fTree->Branch("tpc_sp_ts1", "std::vector", &_tpc_sp_ts1); + fTree->Branch("tpc_sp_score", "std::vector", &_tpc_sp_score); + fTree->Branch("tpc_sp_x", "std::vector", &_tpc_sp_x); + fTree->Branch("tpc_sp_y", "std::vector", &_tpc_sp_y); + fTree->Branch("tpc_sp_z", "std::vector", &_tpc_sp_z); + fTree->Branch("tpc_sp_tag", "std::vector", &_tpc_sp_tag); + fTree->Branch("tpc_tr_matched", "std::vector", &_tpc_tr_matched); + fTree->Branch("tpc_tr_ts0", "std::vector", &_tpc_tr_ts0); + fTree->Branch("tpc_tr_ts1", "std::vector", &_tpc_tr_ts1); + fTree->Branch("tpc_tr_score", "std::vector", &_tpc_tr_score); + fTree->Branch("tpc_tr_start_x", "std::vector", &_tpc_tr_start_x); + fTree->Branch("tpc_tr_start_y", "std::vector", &_tpc_tr_start_y); + fTree->Branch("tpc_tr_start_z", "std::vector", &_tpc_tr_start_z); + fTree->Branch("tpc_tr_dir_x", "std::vector", &_tpc_tr_dir_x); + fTree->Branch("tpc_tr_dir_y", "std::vector", &_tpc_tr_dir_y); + fTree->Branch("tpc_tr_dir_z", "std::vector", &_tpc_tr_dir_z); + fTree->Branch("tpc_tr_tagger1", "std::vector", &_tpc_tr_tagger1); + fTree->Branch("tpc_tr_tagger2", "std::vector", &_tpc_tr_tagger2); + fTree->Branch("tpc_tr_tagger3", "std::vector", &_tpc_tr_tagger3); + fTree->Branch("tpc_diff_tagger_1", "std::vector", &_tpc_diff_tagger_1); + fTree->Branch("tpc_diff1_x", "std::vector", &_tpc_diff1_x); + fTree->Branch("tpc_diff1_y", "std::vector", &_tpc_diff1_y); + fTree->Branch("tpc_diff1_z", "std::vector", &_tpc_diff1_z); + fTree->Branch("tpc_diff_tagger_2", "std::vector", &_tpc_diff_tagger_2); + fTree->Branch("tpc_diff2_x", "std::vector", &_tpc_diff2_x); + fTree->Branch("tpc_diff2_y", "std::vector", &_tpc_diff2_y); + fTree->Branch("tpc_diff2_z", "std::vector", &_tpc_diff2_z); if(!fDataMode) { fTree->Branch("tpc_truth_trackid", "std::vector", &_tpc_truth_trackid); @@ -567,17 +617,9 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fTree->Branch("tpc_truth_energy", "std::vector", &_tpc_truth_energy); fTree->Branch("tpc_truth_time", "std::vector", &_tpc_truth_time); fTree->Branch("tpc_sp_matchable", "std::vector", &_tpc_sp_matchable); - fTree->Branch("tpc_sp_matched", "std::vector", &_tpc_sp_matched); fTree->Branch("tpc_sp_good_match", "std::vector", &_tpc_sp_good_match); - fTree->Branch("tpc_sp_ts0", "std::vector", &_tpc_sp_ts0); - fTree->Branch("tpc_sp_ts1", "std::vector", &_tpc_sp_ts1); - fTree->Branch("tpc_sp_score", "std::vector", &_tpc_sp_score); fTree->Branch("tpc_tr_matchable", "std::vector", &_tpc_tr_matchable); - fTree->Branch("tpc_tr_matched", "std::vector", &_tpc_tr_matched); fTree->Branch("tpc_tr_good_match", "std::vector", &_tpc_tr_good_match); - fTree->Branch("tpc_tr_ts0", "std::vector", &_tpc_tr_ts0); - fTree->Branch("tpc_tr_ts1", "std::vector", &_tpc_tr_ts1); - fTree->Branch("tpc_tr_score", "std::vector", &_tpc_tr_score); } } @@ -1511,12 +1553,33 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_sp_ts0.resize(nTracks); _tpc_sp_ts1.resize(nTracks); _tpc_sp_score.resize(nTracks); + _tpc_sp_x.resize(nTracks); + _tpc_sp_y.resize(nTracks); + _tpc_sp_z.resize(nTracks); + _tpc_sp_tag.resize(nTracks); _tpc_tr_matchable.resize(nTracks); _tpc_tr_matched.resize(nTracks); _tpc_tr_good_match.resize(nTracks); _tpc_tr_ts0.resize(nTracks); _tpc_tr_ts1.resize(nTracks); _tpc_tr_score.resize(nTracks); + _tpc_tr_start_x.resize(nTracks); + _tpc_tr_start_y.resize(nTracks); + _tpc_tr_start_z.resize(nTracks); + _tpc_tr_dir_x.resize(nTracks); + _tpc_tr_dir_y.resize(nTracks); + _tpc_tr_dir_z.resize(nTracks); + _tpc_tr_tagger1.resize(nTracks); + _tpc_tr_tagger2.resize(nTracks); + _tpc_tr_tagger3.resize(nTracks); + _tpc_diff_tagger_1.resize(nTracks); + _tpc_diff1_x.resize(nTracks); + _tpc_diff1_y.resize(nTracks); + _tpc_diff1_z.resize(nTracks); + _tpc_diff_tagger_2.resize(nTracks); + _tpc_diff2_x.resize(nTracks); + _tpc_diff2_y.resize(nTracks); + _tpc_diff2_z.resize(nTracks); art::FindOneP tracksToPFPs(TPCTrackHandle, e, fTPCTrackModuleLabel); art::FindOneP pfpsToMetadata(PFPHandle, e, fPFPModuleLabel); @@ -1561,6 +1624,273 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: else _tpc_track_score[nActualTracks] = -std::numeric_limits::max(); + const art::Ptr spacepoint = tracksToSPMatches.at(track.key()); + const art::Ptr crttrack = tracksToTrackMatches.at(track.key()); + if(spacepoint.isNonnull()) + { + const anab::T0 spMatch = tracksToSPMatches.data(track.key()).ref(); + const art::Ptr cluster = spsToClusters.at(spacepoint.key()); + std::cout << "SPACE POINT TAG = " << cluster->Tagger() << std::endl; + _tpc_sp_matched[nActualTracks] = true; + _tpc_sp_ts0[nActualTracks] = spacepoint->Ts0(); + _tpc_sp_ts1[nActualTracks] = spacepoint->Ts1(); + _tpc_sp_score[nActualTracks] = spMatch.TriggerConfidence(); + _tpc_sp_x[nActualTracks] = spacepoint->X(); + _tpc_sp_y[nActualTracks] = spacepoint->Y(); + _tpc_sp_z[nActualTracks] = spacepoint->Z(); + _tpc_sp_tag[nActualTracks] = cluster->Tagger(); + } + else + { + _tpc_sp_matched[nActualTracks] = false; + _tpc_sp_ts0[nActualTracks] = -std::numeric_limits::max(); + _tpc_sp_ts1[nActualTracks] = -std::numeric_limits::max(); + _tpc_sp_score[nActualTracks] = -std::numeric_limits::max(); + _tpc_sp_x[nActualTracks] = -std::numeric_limits::max(); + _tpc_sp_y[nActualTracks] = -std::numeric_limits::max(); + _tpc_sp_z[nActualTracks] = -std::numeric_limits::max(); + _tpc_sp_tag[nActualTracks] = -1; + } + if(crttrack.isNonnull()) + { + const anab::T0 trackMatch = tracksToTrackMatches.data(track.key()).ref(); + + geo::Point_t entry_tag1, exit_tag1, entry_tag2, exit_tag2, entry_tag3, exit_tag3; + bool hasTag1 = false; + bool hasTag2 = false; + bool hasTag3 = false; + int taggerNum1 = -1; + int taggerNum2 = -1; + int taggerNum3 = -1; + unsigned tag_i = 0; + for(auto const &tagger : crttrack->Taggers()) + { + if(tag_i == 0) { + _tpc_tr_tagger1[i] = tagger; + + for(auto const &[name, tagger2] : fCRTGeoAlg.GetTaggers()) + { + if (CRTCommonUtils::GetTaggerEnum(name) ==tagger) { + std::cout << "name = " << name << " " << CRTCommonUtils::GetTaggerEnum(name) << std::endl; + taggerNum1 = CRTCommonUtils::GetTaggerEnum(name); + + const geo::Point_t min(tagger2.minX, tagger2.minY, tagger2.minZ); + const geo::Point_t max(tagger2.maxX, tagger2.maxY, tagger2.maxZ); + + std::vector> distances; + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kX, min.X()), kX); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kX, max.X()), kX); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kY, min.Y()), kY); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kY, max.Y()), kY); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kZ, min.Z()), kZ); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kZ, max.Z()), kZ); + std::vector chosen_distances; + for(auto const& [k, plane] : distances) + { + const geo::Point_t intersection = start + k * dir; + if(CRTCommonUtils::IsInsideRectangle(min, max, intersection, plane)) + chosen_distances.push_back(k); + } + if(chosen_distances.size() == 0) + std::cout << "Looks like the TPC track wouldn't intersect this CRT plane. (Starting tag for CRT track))" << std::endl; + else if(chosen_distances.size() == 2) + { + entry_tag1 = start + chosen_distances[0] * dir; + exit_tag1 = start + chosen_distances[1] * dir; + hasTag1 = true; + if(chosen_distances[1] < chosen_distances[0]) + std::swap(entry_tag1, exit_tag1); + } + } //tagger match + } //tagger loop + + std::cout << "tagger = " << tagger <<" TPC has tag 1 in this tagger = " << hasTag1 <> distances; + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kX, min.X()), kX); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kX, max.X()), kX); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kY, min.Y()), kY); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kY, max.Y()), kY); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kZ, min.Z()), kZ); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kZ, max.Z()), kZ); + std::vector chosen_distances; + for(auto const& [k, plane] : distances) + { + const geo::Point_t intersection = start + k * dir; + if(CRTCommonUtils::IsInsideRectangle(min, max, intersection, plane)) + chosen_distances.push_back(k); + } + if(chosen_distances.size() == 0) + std::cout << "Looks like the TPC track wouldn't intersect this CRT plane. (Starting tag for CRT track))" << std::endl; + else if(chosen_distances.size() == 2) + { + entry_tag2 = start + chosen_distances[0] * dir; + exit_tag2 = start + chosen_distances[1] * dir; + hasTag2 = true; + if(chosen_distances[1] < chosen_distances[0]) + std::swap(entry_tag2, exit_tag2); + } + } //tagger match + } //tagger loop + std::cout << "tagger = " << tagger <<" TPC has tag 2 in this tagger = " << hasTag2<> distances; + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kX, min.X()), kX); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kX, max.X()), kX); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kY, min.Y()), kY); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kY, max.Y()), kY); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kZ, min.Z()), kZ); + distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kZ, max.Z()), kZ); + std::vector chosen_distances; + for(auto const& [k, plane] : distances) + { + const geo::Point_t intersection = start + k * dir; + if(CRTCommonUtils::IsInsideRectangle(min, max, intersection, plane)) + chosen_distances.push_back(k); + } + if(chosen_distances.size() == 0) + std::cout << "Looks like the TPC track wouldn't intersect this CRT plane. (Starting tag for CRT track))" << std::endl; + else if(chosen_distances.size() == 2) + { + entry_tag3 = start + chosen_distances[0] * dir; + exit_tag3 = start + chosen_distances[1] * dir; + hasTag3 = true; + if(chosen_distances[1] < chosen_distances[0]) + std::swap(entry_tag3, exit_tag3); + } + } //tagger match + } //tagger loop + std::cout << "tagger = " << tagger <<" TPC has tag 3 in this tagger = " << hasTag3 <Ts0(); + _tpc_tr_ts1[nActualTracks] = crttrack->Ts1(); + _tpc_tr_score[nActualTracks] = trackMatch.TriggerConfidence(); + const geo::Point_t startcrt = crttrack->Start(); + const geo::Point_t endcrt = crttrack->End(); + _tpc_tr_start_x[nActualTracks] = startcrt.X(); + _tpc_tr_start_y[nActualTracks] = startcrt.Y(); + _tpc_tr_start_z[nActualTracks] = startcrt.Z(); + const geo::Vector_t dircrt = crttrack->Direction(); + _tpc_tr_dir_x[nActualTracks] = dircrt.X(); + _tpc_tr_dir_y[nActualTracks] = dircrt.Y(); + _tpc_tr_dir_z[nActualTracks] = dircrt.Z(); + if (hasTag3 && !hasTag2) { + entry_tag2 = entry_tag3; + exit_tag2 = exit_tag3; + taggerNum2 = taggerNum3; + } + std::cout << "So the four points then are..." << std::endl; + std::cout << "TPC tags = " << entry_tag1 << " and " << entry_tag2 << std::endl; + std::cout << "CRT tags = " << startcrt << " and " << endcrt << std::endl; + std::cout << "TPC Rs = " << entry_tag1.R() << " and " << entry_tag2.R() << std::endl; + std::cout << "CRT Rs = " << startcrt.R() << " and " << endcrt.R() << std::endl; + if ((entry_tag1.R() < entry_tag2.R()) && (startcrt.R() > endcrt.R())) { + std::cout << "Swapping..." << std::endl; + std::swap(entry_tag1, entry_tag2); + std::swap(taggerNum1, taggerNum2); + } + else if ((entry_tag1.R() > entry_tag2.R()) && (startcrt.R() < endcrt.R())) { + std::cout << "Swapping 2..." << std::endl; + std::swap(entry_tag1, entry_tag2); + std::swap(taggerNum1, taggerNum2); + } + std::cout << " " << std::endl; + std::cout << "So the four points then are..." << std::endl; + std::cout << "TPC tags = " << entry_tag1 << " and " << entry_tag2 << std::endl; + std::cout << "CRT tags = " << startcrt << " and " << endcrt << std::endl; + std::cout << " " << std::endl; + double diff1X = 0.0; + double diff1Y = 0.0; + double diff1Z = 0.0; + double diff2X = 0.0; + double diff2Y = 0.0; + double diff2Z = 0.0; + if (hasTag1 && (hasTag2 || hasTag3)) { + diff1X = entry_tag1.X() - startcrt.X(); + diff1Y = entry_tag1.Y() - startcrt.Y(); + diff1Z = entry_tag1.Z() - startcrt.Z(); + diff2X = entry_tag2.X() - endcrt.X(); + diff2Y = entry_tag2.Y() - endcrt.Y(); + diff2Z = entry_tag2.Z() - endcrt.Z(); + _tpc_diff_tagger_1[nActualTracks] = taggerNum1; + _tpc_diff1_x[nActualTracks] = diff1X; + _tpc_diff1_y[nActualTracks] = diff1Y; + _tpc_diff1_z[nActualTracks] = diff1Z; + _tpc_diff_tagger_2[nActualTracks] = taggerNum2; + _tpc_diff2_x[nActualTracks] = diff2X; + _tpc_diff2_y[nActualTracks] = diff2Y; + _tpc_diff2_z[nActualTracks] = diff2Z; + } + else { + _tpc_diff_tagger_1[nActualTracks] = -1; + _tpc_diff1_x[nActualTracks] = -std::numeric_limits::max(); + _tpc_diff1_y[nActualTracks] = -std::numeric_limits::max(); + _tpc_diff1_z[nActualTracks] = -std::numeric_limits::max(); + _tpc_diff_tagger_2[nActualTracks] = -1; + _tpc_diff2_x[nActualTracks] = -std::numeric_limits::max(); + _tpc_diff2_y[nActualTracks] = -std::numeric_limits::max(); + _tpc_diff2_z[nActualTracks] = -std::numeric_limits::max(); + } + std::cout << "In the 1st tagger " << taggerNum1 << std::endl;; + std::cout << " X TCP - CRT = " << diff1X << std::endl; + std::cout << " Y TCP - CRT = " << diff1Y << std::endl; + std::cout << " Z TCP - CRT = " << diff1Z << std::endl;; + std::cout << "In the 2nd tagger " << taggerNum2 << std::endl; + std::cout << " X TCP - CRT = " << diff2X << std::endl; + std::cout << " Y TCP - CRT = " << diff2Y << std::endl; + std::cout << " Z TCP - CRT = " << diff2Z << std::endl; + std::cout << "----------------------------------------------------" << std::endl; + } + else + { + _tpc_tr_matched[nActualTracks] = false; + _tpc_tr_ts0[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_ts1[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_score[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_start_x[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_start_y[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_start_z[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_dir_x[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_dir_y[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_dir_z[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_tagger1[i] = -1; + _tpc_tr_tagger2[i] = -1; + _tpc_tr_tagger3[i] = -1; + _tpc_diff_tagger_1[nActualTracks] = -1; + _tpc_diff1_x[nActualTracks] = -std::numeric_limits::max(); + _tpc_diff1_y[nActualTracks] = -std::numeric_limits::max(); + _tpc_diff1_z[nActualTracks] = -std::numeric_limits::max(); + _tpc_diff_tagger_2[nActualTracks] = -1; + _tpc_diff2_x[nActualTracks] = -std::numeric_limits::max(); + _tpc_diff2_y[nActualTracks] = -std::numeric_limits::max(); + _tpc_diff2_z[nActualTracks] = -std::numeric_limits::max(); + } + if(!fDataMode) { const std::vector> trackHits = tracksToHits.at(track.key()); @@ -1587,49 +1917,27 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_sp_matchable[nActualTracks] = sp_matchable; _tpc_tr_matchable[nActualTracks] = tr_matchable; - const art::Ptr spacepoint = tracksToSPMatches.at(track.key()); - if(spacepoint.isNonnull()) { - const anab::T0 spMatch = tracksToSPMatches.data(track.key()).ref(); const art::Ptr cluster = spsToClusters.at(spacepoint.key()); const CRTBackTrackerAlg::TruthMatchMetrics truthMatch = fCRTBackTrackerAlg.TruthMatching(e, cluster); - _tpc_sp_matched[nActualTracks] = true; _tpc_sp_good_match[nActualTracks] = truthMatch.trackid == trackid; - _tpc_sp_ts0[nActualTracks] = spacepoint->Ts0(); - _tpc_sp_ts1[nActualTracks] = spacepoint->Ts1(); - _tpc_sp_score[nActualTracks] = spMatch.TriggerConfidence(); } else { - _tpc_sp_matched[nActualTracks] = false; _tpc_sp_good_match[nActualTracks] = false; - _tpc_sp_ts0[nActualTracks] = -std::numeric_limits::max(); - _tpc_sp_ts1[nActualTracks] = -std::numeric_limits::max(); - _tpc_sp_score[nActualTracks] = -std::numeric_limits::max(); } - const art::Ptr crttrack = tracksToTrackMatches.at(track.key()); - if(crttrack.isNonnull()) { - const anab::T0 trackMatch = tracksToTrackMatches.data(track.key()).ref(); const CRTBackTrackerAlg::TruthMatchMetrics truthMatch = fCRTBackTrackerAlg.TruthMatching(e, crttrack); - _tpc_tr_matched[nActualTracks] = true; _tpc_tr_good_match[nActualTracks] = truthMatch.trackid == trackid; - _tpc_tr_ts0[nActualTracks] = crttrack->Ts0(); - _tpc_tr_ts1[nActualTracks] = crttrack->Ts1(); - _tpc_tr_score[nActualTracks] = trackMatch.TriggerConfidence(); } else { - _tpc_tr_matched[nActualTracks] = false; _tpc_tr_good_match[nActualTracks] = false; - _tpc_tr_ts0[nActualTracks] = -std::numeric_limits::max(); - _tpc_tr_ts1[nActualTracks] = -std::numeric_limits::max(); - _tpc_tr_score[nActualTracks] = -std::numeric_limits::max(); } } @@ -1657,12 +1965,33 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_sp_ts0.resize(nActualTracks); _tpc_sp_ts1.resize(nActualTracks); _tpc_sp_score.resize(nActualTracks); + _tpc_sp_x.resize(nActualTracks); + _tpc_sp_y.resize(nActualTracks); + _tpc_sp_z.resize(nActualTracks); + _tpc_sp_tag.resize(nActualTracks); _tpc_tr_matchable.resize(nActualTracks); _tpc_tr_matched.resize(nActualTracks); _tpc_tr_good_match.resize(nActualTracks); _tpc_tr_ts0.resize(nActualTracks); _tpc_tr_ts1.resize(nActualTracks); _tpc_tr_score.resize(nActualTracks); + _tpc_tr_start_x.resize(nActualTracks); + _tpc_tr_start_y.resize(nActualTracks); + _tpc_tr_start_z.resize(nActualTracks); + _tpc_tr_dir_x.resize(nActualTracks); + _tpc_tr_dir_y.resize(nActualTracks); + _tpc_tr_dir_z.resize(nActualTracks); + _tpc_tr_tagger1.resize(nActualTracks); + _tpc_tr_tagger2.resize(nActualTracks); + _tpc_tr_tagger3.resize(nActualTracks); + _tpc_diff_tagger_1.resize(nActualTracks); + _tpc_diff1_x.resize(nActualTracks); + _tpc_diff1_y.resize(nActualTracks); + _tpc_diff1_z.resize(nActualTracks); + _tpc_diff_tagger_2.resize(nActualTracks); + _tpc_diff2_x.resize(nActualTracks); + _tpc_diff2_y.resize(nActualTracks); + _tpc_diff2_z.resize(nActualTracks); } DEFINE_ART_MODULE(sbnd::crt::CRTAnalysis) diff --git a/sbndcode/CRT/CRTAna/run_crtana_data.fcl b/sbndcode/CRT/CRTAna/run_crtana_data.fcl index b2ade27a9..77e03ee4e 100644 --- a/sbndcode/CRT/CRTAna/run_crtana_data.fcl +++ b/sbndcode/CRT/CRTAna/run_crtana_data.fcl @@ -8,4 +8,11 @@ services.BackTrackerService: @erase services.CRTChannelMapService: @local::crt_channel_map_standard services.CRTCalibService: @local::crt_calib_service -physics.analyzers.crtana: @local::crtana_data_sbnd +physics.analyzers.crtana.FEBDataModuleLabel: "crtdecoder" +physics.analyzers.crtana.DataMode: true +physics.analyzers.crtana.NoTPC: false +physics.analyzers.crtana.HasPTB: true +physics.analyzers.crtana.HasTDC: true +physics.analyzers.crtana.TPCTrackModuleLabel: "pandoraTrack" +physics.analyzers.crtana.CRTSpacePointMatchingModuleLabel: "crtspacepointmatching" +physics.analyzers.crtana.CRTTrackMatchingModuleLabel: "crttrackmatching" diff --git a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc index 4ef671723..bdece4f77 100644 --- a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc +++ b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc @@ -1,4 +1,9 @@ #include "CRTEventDisplayAlg.h" +#include "lardataobj/AnalysisBase/T0.h" +#include "canvas/Persistency/Common/FindManyP.h" +#include "lardataobj/RecoBase/Track.h" +#include "TH1D.h" +#include "TPaveText.h" namespace sbnd::crt { @@ -23,6 +28,7 @@ namespace sbnd::crt { fClusterLabel = config.ClusterLabel(); fSpacePointLabel = config.SpacePointLabel(); fTrackLabel = config.TrackLabel(); + fTrackMatchLabel = config.TrackMatchLabel(); fSaveRoot = config.SaveRoot(); fSaveViews = config.SaveViews(); @@ -463,6 +469,54 @@ namespace sbnd::crt { const art::Ptr spacepoint = spacePointVec[0]; const geo::Point_t pos = spacepoint->Pos(); const geo::Point_t err = spacepoint->Err(); + try { + auto spacePointsHandle = event.getValidHandle>(fSpacePointLabel); + art::FindOneP CRTSPstoTPCTracks(spacePointsHandle, event, "crtspacepointmatching"); + const art::Ptr TPCTrack = CRTSPstoTPCTracks.at(spacepoint.key()); + if(TPCTrack.isNonnull()) { + const anab::T0 t0Match = CRTSPstoTPCTracks.data(spacepoint.key()).ref(); + std::cout << "t0 SP match (confidence) = " << t0Match.TriggerConfidence() << std::endl; + double t0SPMatchConfidence = t0Match.TriggerConfidence(); + const geo::Point_t startTPC = TPCTrack->Start(); + const geo::Vector_t dirTPC = TPCTrack->StartDirection(); + TPolyLine3D *lineTPC = new TPolyLine3D(2); + geo::Point_t aTPC {0,0,0}; + geo::Point_t bTPC {0,0,0}; + int i = 0; + do + { + aTPC = startTPC + i * dirTPC; + ++i; + } + while(IsPointInsideBox(crtLims, aTPC)); + i = 0; + do + { + bTPC = startTPC + i * dirTPC; + --i; + } + while(IsPointInsideBox(crtLims, bTPC)); + lineTPC->SetPoint(0, aTPC.X(), aTPC.Y(), aTPC.Z()); + lineTPC->SetPoint(1, bTPC.X(), bTPC.Y(), bTPC.Z()); + lineTPC->SetLineColor(3); + lineTPC->SetLineWidth(fLineWidth); + lineTPC->Draw(); + TPaveText *pt = new TPaveText(0.05,0.85,0.35,0.65,"NB"); + pt->SetTextSize(0.02); + pt->SetFillStyle(0); + pt->SetLineStyle(0); + pt->SetTextAlign(12); + pt->SetBorderSize(0); + pt->AddText(Form("t0SPMatchConfidence=%g", t0SPMatchConfidence)); + pt->Draw(); + } + else { + continue; + } + } catch(...) { + continue; + } + double rmin[3] = {pos.X() - err.X(), pos.Y() - err.Y(), pos.Z() - err.Z()}; double rmax[3] = {pos.X() + err.X(), pos.Y() + err.Y(), pos.Z() + err.Z()}; @@ -489,6 +543,8 @@ namespace sbnd::crt { auto tracksHandle = event.getValidHandle>(fTrackLabel); std::vector> tracksVec; art::fill_ptr_vector(tracksVec, tracksHandle); + art::FindOneP CRTTrackstoTPCTracks(tracksHandle, event, "crttrackmatching"); //fTrackMatchLabel + for(auto track : tracksVec) { @@ -511,6 +567,49 @@ namespace sbnd::crt { const geo::Point_t start = track->Start(); const geo::Vector_t dir = track->Direction(); + const art::Ptr TPCTrack = CRTTrackstoTPCTracks.at(track.key()); + if(TPCTrack.isNonnull()) { + const anab::T0 t0Match = CRTTrackstoTPCTracks.data(track.key()).ref(); + std::cout << "t0 match (confidence) = " << t0Match.TriggerConfidence() << std::endl; + double t0MatchConfidence = t0Match.TriggerConfidence(); + + const geo::Point_t startTPC = TPCTrack->Start(); + const geo::Vector_t dirTPC = TPCTrack->StartDirection(); + TPolyLine3D *lineTPC = new TPolyLine3D(2); + geo::Point_t aTPC {0,0,0}; + geo::Point_t bTPC {0,0,0}; + int i = 0; + do + { + aTPC = startTPC + i * dirTPC; + ++i; + } + while(IsPointInsideBox(crtLims, aTPC)); + i = 0; + do + { + bTPC = startTPC + i * dirTPC; + --i; + } + while(IsPointInsideBox(crtLims, bTPC)); + lineTPC->SetPoint(0, aTPC.X(), aTPC.Y(), aTPC.Z()); + lineTPC->SetPoint(1, bTPC.X(), bTPC.Y(), bTPC.Z()); + lineTPC->SetLineColor(2); + lineTPC->SetLineWidth(fLineWidth); + lineTPC->Draw(); + TPaveText *pt = new TPaveText(0.05,0.8,0.35,0.6,"NB"); + pt->SetTextSize(0.02); + pt->SetFillStyle(0); + pt->SetLineStyle(0); + pt->SetTextAlign(12); + pt->SetBorderSize(0); + pt->AddText(Form("t0MatchConfidence=%g", t0MatchConfidence)); + pt->Draw(); + } + else { + continue; + } + TPolyLine3D *line = new TPolyLine3D(2); geo::Point_t a {0,0,0}; @@ -551,6 +650,7 @@ namespace sbnd::crt { } } + std::cout << Form("%s", saveName.Data()) << std::endl; if(fSaveRoot) c1->SaveAs(Form("%s.root", saveName.Data())); diff --git a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h index 9429bd033..b1875fb6e 100644 --- a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h +++ b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h @@ -87,6 +87,9 @@ namespace sbnd::crt { fhicl::Atom TrackLabel { Name("TrackLabel") }; + fhicl::Atom TrackMatchLabel { + Name("TrackMatchLabel") + }; fhicl::Atom SaveRoot { Name("SaveRoot"), @@ -247,6 +250,7 @@ namespace sbnd::crt { art::InputTag fClusterLabel; art::InputTag fSpacePointLabel; art::InputTag fTrackLabel; + art::InputTag fTrackMatchLabel; bool fSaveRoot; bool fSaveViews; diff --git a/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl b/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl index 5eb31515e..f73721397 100644 --- a/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl +++ b/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl @@ -17,6 +17,7 @@ crteventdisplayalg_sbnd: ClusterLabel: "crtclustering" SpacePointLabel: "crtspacepoints" TrackLabel: "crttracks" + TrackMatchLabel: "crttrackmatching" SaveRoot: true SaveViews: false diff --git a/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl b/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl index 77f35b4bc..6d7c5fcd2 100644 --- a/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl +++ b/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl @@ -3,7 +3,7 @@ BEGIN_PROLOG crtspacepointmatchalg_sbnd: { TrackDirectionFrac: 0.5 - DCALimit: 70. + DCALimit: 140. MinTPCTrackLength: 10. DirMethod: 1 DCAuseBox: false @@ -12,7 +12,7 @@ crtspacepointmatchalg_sbnd: MaxUncert: 20. TPCTrackLabel: "pandoraTrack" CRTSpacePointLabel: "crtspacepoints" - UseTs0: false + UseTs0: true } crtspacepointmatchproducer_sbnd: @@ -26,13 +26,13 @@ crtspacepointmatchproducer_sbnd: crttrackmatchalg_sbnd: { - MaxAngleDiff: 0.4 - MaxDCA: 70. - MaxScore: 150. + MaxAngleDiff: 0.8 + MaxDCA: 140. + MaxScore: 300. MinTPCTrackLength: 5. SelectionMetric: "score" TPCTrackLabel: "pandoraTrack" - UseTs0: false + UseTs0: true } crttrackmatchproducer_sbnd: diff --git a/sbndcode/Calibration/CRT/CalibService/run_caloskim_data.fcl b/sbndcode/Calibration/CRT/CalibService/run_caloskim_data.fcl new file mode 100644 index 000000000..a0f818f4b --- /dev/null +++ b/sbndcode/Calibration/CRT/CalibService/run_caloskim_data.fcl @@ -0,0 +1,10 @@ +#include "run_calorimetry_skim_sbnd.fcl" + +services.ParticleInventoryService: @erase +services.BackTrackerService: @erase + +physics.analyzers.caloskim.G4producer: "" +physics.analyzers.caloskim.SimChannelproducer: "" + +#outputs.outpid.dataTier: "reconstructed" +#outputs.outpid.fastCloning: true From 8665f8b9ac68af126f19775db5a0f399f45d4f88 Mon Sep 17 00:00:00 2001 From: Alexandra Moor Date: Tue, 25 Feb 2025 08:51:32 -0600 Subject: [PATCH 02/18] clean up to basics --- sbndcode/CRT/CRTAna/CRTAnalysis_module.cc | 319 +----------------- .../CRT/CRTEventDisplay/CRTEventDisplayAlg.cc | 11 +- .../crttpcmatchingproducers_sbnd.fcl | 8 +- .../CRT/CalibService/run_caloskim_data.fcl | 10 - 4 files changed, 11 insertions(+), 337 deletions(-) delete mode 100644 sbndcode/Calibration/CRT/CalibService/run_caloskim_data.fcl diff --git a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc index 9d1f42c06..793d69eef 100644 --- a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc +++ b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc @@ -305,33 +305,12 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer { std::vector _tpc_sp_ts0; std::vector _tpc_sp_ts1; std::vector _tpc_sp_score; - std::vector _tpc_sp_x; - std::vector _tpc_sp_y; - std::vector _tpc_sp_z; - std::vector _tpc_sp_tag; std::vector _tpc_tr_matchable; std::vector _tpc_tr_matched; std::vector _tpc_tr_good_match; std::vector _tpc_tr_ts0; std::vector _tpc_tr_ts1; std::vector _tpc_tr_score; - std::vector _tpc_tr_start_x; - std::vector _tpc_tr_start_y; - std::vector _tpc_tr_start_z; - std::vector _tpc_tr_dir_x; - std::vector _tpc_tr_dir_y; - std::vector _tpc_tr_dir_z; - std::vector _tpc_tr_tagger1; - std::vector _tpc_tr_tagger2; - std::vector _tpc_tr_tagger3; - std::vector _tpc_diff_tagger_1; - std::vector _tpc_diff1_x; - std::vector _tpc_diff1_y; - std::vector _tpc_diff1_z; - std::vector _tpc_diff_tagger_2; - std::vector _tpc_diff2_x; - std::vector _tpc_diff2_y; - std::vector _tpc_diff2_z; std::vector _ptb_hlt_trigger; std::vector _ptb_hlt_timestamp; @@ -585,31 +564,10 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fTree->Branch("tpc_sp_ts0", "std::vector", &_tpc_sp_ts0); fTree->Branch("tpc_sp_ts1", "std::vector", &_tpc_sp_ts1); fTree->Branch("tpc_sp_score", "std::vector", &_tpc_sp_score); - fTree->Branch("tpc_sp_x", "std::vector", &_tpc_sp_x); - fTree->Branch("tpc_sp_y", "std::vector", &_tpc_sp_y); - fTree->Branch("tpc_sp_z", "std::vector", &_tpc_sp_z); - fTree->Branch("tpc_sp_tag", "std::vector", &_tpc_sp_tag); fTree->Branch("tpc_tr_matched", "std::vector", &_tpc_tr_matched); fTree->Branch("tpc_tr_ts0", "std::vector", &_tpc_tr_ts0); fTree->Branch("tpc_tr_ts1", "std::vector", &_tpc_tr_ts1); fTree->Branch("tpc_tr_score", "std::vector", &_tpc_tr_score); - fTree->Branch("tpc_tr_start_x", "std::vector", &_tpc_tr_start_x); - fTree->Branch("tpc_tr_start_y", "std::vector", &_tpc_tr_start_y); - fTree->Branch("tpc_tr_start_z", "std::vector", &_tpc_tr_start_z); - fTree->Branch("tpc_tr_dir_x", "std::vector", &_tpc_tr_dir_x); - fTree->Branch("tpc_tr_dir_y", "std::vector", &_tpc_tr_dir_y); - fTree->Branch("tpc_tr_dir_z", "std::vector", &_tpc_tr_dir_z); - fTree->Branch("tpc_tr_tagger1", "std::vector", &_tpc_tr_tagger1); - fTree->Branch("tpc_tr_tagger2", "std::vector", &_tpc_tr_tagger2); - fTree->Branch("tpc_tr_tagger3", "std::vector", &_tpc_tr_tagger3); - fTree->Branch("tpc_diff_tagger_1", "std::vector", &_tpc_diff_tagger_1); - fTree->Branch("tpc_diff1_x", "std::vector", &_tpc_diff1_x); - fTree->Branch("tpc_diff1_y", "std::vector", &_tpc_diff1_y); - fTree->Branch("tpc_diff1_z", "std::vector", &_tpc_diff1_z); - fTree->Branch("tpc_diff_tagger_2", "std::vector", &_tpc_diff_tagger_2); - fTree->Branch("tpc_diff2_x", "std::vector", &_tpc_diff2_x); - fTree->Branch("tpc_diff2_y", "std::vector", &_tpc_diff2_y); - fTree->Branch("tpc_diff2_z", "std::vector", &_tpc_diff2_z); if(!fDataMode) { fTree->Branch("tpc_truth_trackid", "std::vector", &_tpc_truth_trackid); @@ -1553,33 +1511,12 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_sp_ts0.resize(nTracks); _tpc_sp_ts1.resize(nTracks); _tpc_sp_score.resize(nTracks); - _tpc_sp_x.resize(nTracks); - _tpc_sp_y.resize(nTracks); - _tpc_sp_z.resize(nTracks); - _tpc_sp_tag.resize(nTracks); _tpc_tr_matchable.resize(nTracks); _tpc_tr_matched.resize(nTracks); _tpc_tr_good_match.resize(nTracks); _tpc_tr_ts0.resize(nTracks); _tpc_tr_ts1.resize(nTracks); _tpc_tr_score.resize(nTracks); - _tpc_tr_start_x.resize(nTracks); - _tpc_tr_start_y.resize(nTracks); - _tpc_tr_start_z.resize(nTracks); - _tpc_tr_dir_x.resize(nTracks); - _tpc_tr_dir_y.resize(nTracks); - _tpc_tr_dir_z.resize(nTracks); - _tpc_tr_tagger1.resize(nTracks); - _tpc_tr_tagger2.resize(nTracks); - _tpc_tr_tagger3.resize(nTracks); - _tpc_diff_tagger_1.resize(nTracks); - _tpc_diff1_x.resize(nTracks); - _tpc_diff1_y.resize(nTracks); - _tpc_diff1_z.resize(nTracks); - _tpc_diff_tagger_2.resize(nTracks); - _tpc_diff2_x.resize(nTracks); - _tpc_diff2_y.resize(nTracks); - _tpc_diff2_z.resize(nTracks); art::FindOneP tracksToPFPs(TPCTrackHandle, e, fTPCTrackModuleLabel); art::FindOneP pfpsToMetadata(PFPHandle, e, fPFPModuleLabel); @@ -1629,16 +1566,11 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: if(spacepoint.isNonnull()) { const anab::T0 spMatch = tracksToSPMatches.data(track.key()).ref(); - const art::Ptr cluster = spsToClusters.at(spacepoint.key()); - std::cout << "SPACE POINT TAG = " << cluster->Tagger() << std::endl; + _tpc_sp_matched[nActualTracks] = true; _tpc_sp_ts0[nActualTracks] = spacepoint->Ts0(); _tpc_sp_ts1[nActualTracks] = spacepoint->Ts1(); _tpc_sp_score[nActualTracks] = spMatch.TriggerConfidence(); - _tpc_sp_x[nActualTracks] = spacepoint->X(); - _tpc_sp_y[nActualTracks] = spacepoint->Y(); - _tpc_sp_z[nActualTracks] = spacepoint->Z(); - _tpc_sp_tag[nActualTracks] = cluster->Tagger(); } else { @@ -1646,225 +1578,16 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_sp_ts0[nActualTracks] = -std::numeric_limits::max(); _tpc_sp_ts1[nActualTracks] = -std::numeric_limits::max(); _tpc_sp_score[nActualTracks] = -std::numeric_limits::max(); - _tpc_sp_x[nActualTracks] = -std::numeric_limits::max(); - _tpc_sp_y[nActualTracks] = -std::numeric_limits::max(); - _tpc_sp_z[nActualTracks] = -std::numeric_limits::max(); - _tpc_sp_tag[nActualTracks] = -1; } + if(crttrack.isNonnull()) { const anab::T0 trackMatch = tracksToTrackMatches.data(track.key()).ref(); - geo::Point_t entry_tag1, exit_tag1, entry_tag2, exit_tag2, entry_tag3, exit_tag3; - bool hasTag1 = false; - bool hasTag2 = false; - bool hasTag3 = false; - int taggerNum1 = -1; - int taggerNum2 = -1; - int taggerNum3 = -1; - unsigned tag_i = 0; - for(auto const &tagger : crttrack->Taggers()) - { - if(tag_i == 0) { - _tpc_tr_tagger1[i] = tagger; - - for(auto const &[name, tagger2] : fCRTGeoAlg.GetTaggers()) - { - if (CRTCommonUtils::GetTaggerEnum(name) ==tagger) { - std::cout << "name = " << name << " " << CRTCommonUtils::GetTaggerEnum(name) << std::endl; - taggerNum1 = CRTCommonUtils::GetTaggerEnum(name); - - const geo::Point_t min(tagger2.minX, tagger2.minY, tagger2.minZ); - const geo::Point_t max(tagger2.maxX, tagger2.maxY, tagger2.maxZ); - - std::vector> distances; - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kX, min.X()), kX); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kX, max.X()), kX); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kY, min.Y()), kY); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kY, max.Y()), kY); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kZ, min.Z()), kZ); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kZ, max.Z()), kZ); - std::vector chosen_distances; - for(auto const& [k, plane] : distances) - { - const geo::Point_t intersection = start + k * dir; - if(CRTCommonUtils::IsInsideRectangle(min, max, intersection, plane)) - chosen_distances.push_back(k); - } - if(chosen_distances.size() == 0) - std::cout << "Looks like the TPC track wouldn't intersect this CRT plane. (Starting tag for CRT track))" << std::endl; - else if(chosen_distances.size() == 2) - { - entry_tag1 = start + chosen_distances[0] * dir; - exit_tag1 = start + chosen_distances[1] * dir; - hasTag1 = true; - if(chosen_distances[1] < chosen_distances[0]) - std::swap(entry_tag1, exit_tag1); - } - } //tagger match - } //tagger loop - - std::cout << "tagger = " << tagger <<" TPC has tag 1 in this tagger = " << hasTag1 <> distances; - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kX, min.X()), kX); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kX, max.X()), kX); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kY, min.Y()), kY); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kY, max.Y()), kY); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kZ, min.Z()), kZ); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kZ, max.Z()), kZ); - std::vector chosen_distances; - for(auto const& [k, plane] : distances) - { - const geo::Point_t intersection = start + k * dir; - if(CRTCommonUtils::IsInsideRectangle(min, max, intersection, plane)) - chosen_distances.push_back(k); - } - if(chosen_distances.size() == 0) - std::cout << "Looks like the TPC track wouldn't intersect this CRT plane. (Starting tag for CRT track))" << std::endl; - else if(chosen_distances.size() == 2) - { - entry_tag2 = start + chosen_distances[0] * dir; - exit_tag2 = start + chosen_distances[1] * dir; - hasTag2 = true; - if(chosen_distances[1] < chosen_distances[0]) - std::swap(entry_tag2, exit_tag2); - } - } //tagger match - } //tagger loop - std::cout << "tagger = " << tagger <<" TPC has tag 2 in this tagger = " << hasTag2<> distances; - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kX, min.X()), kX); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kX, max.X()), kX); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kY, min.Y()), kY); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kY, max.Y()), kY); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kZ, min.Z()), kZ); - distances.emplace_back(CRTCommonUtils::LinePlaneIntersection(start, dir, kZ, max.Z()), kZ); - std::vector chosen_distances; - for(auto const& [k, plane] : distances) - { - const geo::Point_t intersection = start + k * dir; - if(CRTCommonUtils::IsInsideRectangle(min, max, intersection, plane)) - chosen_distances.push_back(k); - } - if(chosen_distances.size() == 0) - std::cout << "Looks like the TPC track wouldn't intersect this CRT plane. (Starting tag for CRT track))" << std::endl; - else if(chosen_distances.size() == 2) - { - entry_tag3 = start + chosen_distances[0] * dir; - exit_tag3 = start + chosen_distances[1] * dir; - hasTag3 = true; - if(chosen_distances[1] < chosen_distances[0]) - std::swap(entry_tag3, exit_tag3); - } - } //tagger match - } //tagger loop - std::cout << "tagger = " << tagger <<" TPC has tag 3 in this tagger = " << hasTag3 <Ts0(); _tpc_tr_ts1[nActualTracks] = crttrack->Ts1(); _tpc_tr_score[nActualTracks] = trackMatch.TriggerConfidence(); - const geo::Point_t startcrt = crttrack->Start(); - const geo::Point_t endcrt = crttrack->End(); - _tpc_tr_start_x[nActualTracks] = startcrt.X(); - _tpc_tr_start_y[nActualTracks] = startcrt.Y(); - _tpc_tr_start_z[nActualTracks] = startcrt.Z(); - const geo::Vector_t dircrt = crttrack->Direction(); - _tpc_tr_dir_x[nActualTracks] = dircrt.X(); - _tpc_tr_dir_y[nActualTracks] = dircrt.Y(); - _tpc_tr_dir_z[nActualTracks] = dircrt.Z(); - if (hasTag3 && !hasTag2) { - entry_tag2 = entry_tag3; - exit_tag2 = exit_tag3; - taggerNum2 = taggerNum3; - } - std::cout << "So the four points then are..." << std::endl; - std::cout << "TPC tags = " << entry_tag1 << " and " << entry_tag2 << std::endl; - std::cout << "CRT tags = " << startcrt << " and " << endcrt << std::endl; - std::cout << "TPC Rs = " << entry_tag1.R() << " and " << entry_tag2.R() << std::endl; - std::cout << "CRT Rs = " << startcrt.R() << " and " << endcrt.R() << std::endl; - if ((entry_tag1.R() < entry_tag2.R()) && (startcrt.R() > endcrt.R())) { - std::cout << "Swapping..." << std::endl; - std::swap(entry_tag1, entry_tag2); - std::swap(taggerNum1, taggerNum2); - } - else if ((entry_tag1.R() > entry_tag2.R()) && (startcrt.R() < endcrt.R())) { - std::cout << "Swapping 2..." << std::endl; - std::swap(entry_tag1, entry_tag2); - std::swap(taggerNum1, taggerNum2); - } - std::cout << " " << std::endl; - std::cout << "So the four points then are..." << std::endl; - std::cout << "TPC tags = " << entry_tag1 << " and " << entry_tag2 << std::endl; - std::cout << "CRT tags = " << startcrt << " and " << endcrt << std::endl; - std::cout << " " << std::endl; - double diff1X = 0.0; - double diff1Y = 0.0; - double diff1Z = 0.0; - double diff2X = 0.0; - double diff2Y = 0.0; - double diff2Z = 0.0; - if (hasTag1 && (hasTag2 || hasTag3)) { - diff1X = entry_tag1.X() - startcrt.X(); - diff1Y = entry_tag1.Y() - startcrt.Y(); - diff1Z = entry_tag1.Z() - startcrt.Z(); - diff2X = entry_tag2.X() - endcrt.X(); - diff2Y = entry_tag2.Y() - endcrt.Y(); - diff2Z = entry_tag2.Z() - endcrt.Z(); - _tpc_diff_tagger_1[nActualTracks] = taggerNum1; - _tpc_diff1_x[nActualTracks] = diff1X; - _tpc_diff1_y[nActualTracks] = diff1Y; - _tpc_diff1_z[nActualTracks] = diff1Z; - _tpc_diff_tagger_2[nActualTracks] = taggerNum2; - _tpc_diff2_x[nActualTracks] = diff2X; - _tpc_diff2_y[nActualTracks] = diff2Y; - _tpc_diff2_z[nActualTracks] = diff2Z; - } - else { - _tpc_diff_tagger_1[nActualTracks] = -1; - _tpc_diff1_x[nActualTracks] = -std::numeric_limits::max(); - _tpc_diff1_y[nActualTracks] = -std::numeric_limits::max(); - _tpc_diff1_z[nActualTracks] = -std::numeric_limits::max(); - _tpc_diff_tagger_2[nActualTracks] = -1; - _tpc_diff2_x[nActualTracks] = -std::numeric_limits::max(); - _tpc_diff2_y[nActualTracks] = -std::numeric_limits::max(); - _tpc_diff2_z[nActualTracks] = -std::numeric_limits::max(); - } - std::cout << "In the 1st tagger " << taggerNum1 << std::endl;; - std::cout << " X TCP - CRT = " << diff1X << std::endl; - std::cout << " Y TCP - CRT = " << diff1Y << std::endl; - std::cout << " Z TCP - CRT = " << diff1Z << std::endl;; - std::cout << "In the 2nd tagger " << taggerNum2 << std::endl; - std::cout << " X TCP - CRT = " << diff2X << std::endl; - std::cout << " Y TCP - CRT = " << diff2Y << std::endl; - std::cout << " Z TCP - CRT = " << diff2Z << std::endl; - std::cout << "----------------------------------------------------" << std::endl; } else { @@ -1872,23 +1595,6 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_tr_ts0[nActualTracks] = -std::numeric_limits::max(); _tpc_tr_ts1[nActualTracks] = -std::numeric_limits::max(); _tpc_tr_score[nActualTracks] = -std::numeric_limits::max(); - _tpc_tr_start_x[nActualTracks] = -std::numeric_limits::max(); - _tpc_tr_start_y[nActualTracks] = -std::numeric_limits::max(); - _tpc_tr_start_z[nActualTracks] = -std::numeric_limits::max(); - _tpc_tr_dir_x[nActualTracks] = -std::numeric_limits::max(); - _tpc_tr_dir_y[nActualTracks] = -std::numeric_limits::max(); - _tpc_tr_dir_z[nActualTracks] = -std::numeric_limits::max(); - _tpc_tr_tagger1[i] = -1; - _tpc_tr_tagger2[i] = -1; - _tpc_tr_tagger3[i] = -1; - _tpc_diff_tagger_1[nActualTracks] = -1; - _tpc_diff1_x[nActualTracks] = -std::numeric_limits::max(); - _tpc_diff1_y[nActualTracks] = -std::numeric_limits::max(); - _tpc_diff1_z[nActualTracks] = -std::numeric_limits::max(); - _tpc_diff_tagger_2[nActualTracks] = -1; - _tpc_diff2_x[nActualTracks] = -std::numeric_limits::max(); - _tpc_diff2_y[nActualTracks] = -std::numeric_limits::max(); - _tpc_diff2_z[nActualTracks] = -std::numeric_limits::max(); } if(!fDataMode) @@ -1965,33 +1671,12 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_sp_ts0.resize(nActualTracks); _tpc_sp_ts1.resize(nActualTracks); _tpc_sp_score.resize(nActualTracks); - _tpc_sp_x.resize(nActualTracks); - _tpc_sp_y.resize(nActualTracks); - _tpc_sp_z.resize(nActualTracks); - _tpc_sp_tag.resize(nActualTracks); _tpc_tr_matchable.resize(nActualTracks); _tpc_tr_matched.resize(nActualTracks); _tpc_tr_good_match.resize(nActualTracks); _tpc_tr_ts0.resize(nActualTracks); _tpc_tr_ts1.resize(nActualTracks); _tpc_tr_score.resize(nActualTracks); - _tpc_tr_start_x.resize(nActualTracks); - _tpc_tr_start_y.resize(nActualTracks); - _tpc_tr_start_z.resize(nActualTracks); - _tpc_tr_dir_x.resize(nActualTracks); - _tpc_tr_dir_y.resize(nActualTracks); - _tpc_tr_dir_z.resize(nActualTracks); - _tpc_tr_tagger1.resize(nActualTracks); - _tpc_tr_tagger2.resize(nActualTracks); - _tpc_tr_tagger3.resize(nActualTracks); - _tpc_diff_tagger_1.resize(nActualTracks); - _tpc_diff1_x.resize(nActualTracks); - _tpc_diff1_y.resize(nActualTracks); - _tpc_diff1_z.resize(nActualTracks); - _tpc_diff_tagger_2.resize(nActualTracks); - _tpc_diff2_x.resize(nActualTracks); - _tpc_diff2_y.resize(nActualTracks); - _tpc_diff2_z.resize(nActualTracks); } DEFINE_ART_MODULE(sbnd::crt::CRTAnalysis) diff --git a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc index bdece4f77..842d68c59 100644 --- a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc +++ b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc @@ -475,7 +475,6 @@ namespace sbnd::crt { const art::Ptr TPCTrack = CRTSPstoTPCTracks.at(spacepoint.key()); if(TPCTrack.isNonnull()) { const anab::T0 t0Match = CRTSPstoTPCTracks.data(spacepoint.key()).ref(); - std::cout << "t0 SP match (confidence) = " << t0Match.TriggerConfidence() << std::endl; double t0SPMatchConfidence = t0Match.TriggerConfidence(); const geo::Point_t startTPC = TPCTrack->Start(); const geo::Vector_t dirTPC = TPCTrack->StartDirection(); @@ -507,8 +506,9 @@ namespace sbnd::crt { pt->SetLineStyle(0); pt->SetTextAlign(12); pt->SetBorderSize(0); - pt->AddText(Form("t0SPMatchConfidence=%g", t0SPMatchConfidence)); pt->Draw(); + if(fPrint) + std::cout << "t0 matching confidence for this spacepoint is " << t0SPMatchConfidence << std::endl; } else { continue; @@ -543,7 +543,7 @@ namespace sbnd::crt { auto tracksHandle = event.getValidHandle>(fTrackLabel); std::vector> tracksVec; art::fill_ptr_vector(tracksVec, tracksHandle); - art::FindOneP CRTTrackstoTPCTracks(tracksHandle, event, "crttrackmatching"); //fTrackMatchLabel + art::FindOneP CRTTrackstoTPCTracks(tracksHandle, event, "crttrackmatching"); for(auto track : tracksVec) @@ -570,7 +570,6 @@ namespace sbnd::crt { const art::Ptr TPCTrack = CRTTrackstoTPCTracks.at(track.key()); if(TPCTrack.isNonnull()) { const anab::T0 t0Match = CRTTrackstoTPCTracks.data(track.key()).ref(); - std::cout << "t0 match (confidence) = " << t0Match.TriggerConfidence() << std::endl; double t0MatchConfidence = t0Match.TriggerConfidence(); const geo::Point_t startTPC = TPCTrack->Start(); @@ -603,8 +602,9 @@ namespace sbnd::crt { pt->SetLineStyle(0); pt->SetTextAlign(12); pt->SetBorderSize(0); - pt->AddText(Form("t0MatchConfidence=%g", t0MatchConfidence)); pt->Draw(); + if(fPrint) + std::cout << "t0 matching confidence for this track is " << t0MatchConfidence << std::endl; } else { continue; @@ -650,7 +650,6 @@ namespace sbnd::crt { } } - std::cout << Form("%s", saveName.Data()) << std::endl; if(fSaveRoot) c1->SaveAs(Form("%s.root", saveName.Data())); diff --git a/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl b/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl index 6d7c5fcd2..1144a0df9 100644 --- a/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl +++ b/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl @@ -3,7 +3,7 @@ BEGIN_PROLOG crtspacepointmatchalg_sbnd: { TrackDirectionFrac: 0.5 - DCALimit: 140. + DCALimit: 70. MinTPCTrackLength: 10. DirMethod: 1 DCAuseBox: false @@ -26,9 +26,9 @@ crtspacepointmatchproducer_sbnd: crttrackmatchalg_sbnd: { - MaxAngleDiff: 0.8 - MaxDCA: 140. - MaxScore: 300. + MaxAngleDiff: 0.4 + MaxDCA: 70. + MaxScore: 150. MinTPCTrackLength: 5. SelectionMetric: "score" TPCTrackLabel: "pandoraTrack" diff --git a/sbndcode/Calibration/CRT/CalibService/run_caloskim_data.fcl b/sbndcode/Calibration/CRT/CalibService/run_caloskim_data.fcl deleted file mode 100644 index a0f818f4b..000000000 --- a/sbndcode/Calibration/CRT/CalibService/run_caloskim_data.fcl +++ /dev/null @@ -1,10 +0,0 @@ -#include "run_calorimetry_skim_sbnd.fcl" - -services.ParticleInventoryService: @erase -services.BackTrackerService: @erase - -physics.analyzers.caloskim.G4producer: "" -physics.analyzers.caloskim.SimChannelproducer: "" - -#outputs.outpid.dataTier: "reconstructed" -#outputs.outpid.fastCloning: true From e333d76a0bbd4575162fb7706213b7442dec11ae Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Tue, 25 Feb 2025 10:19:17 -0600 Subject: [PATCH 03/18] Move changes to independent fcl to preserve existing behaviour --- sbndcode/CRT/CRTAna/CRTAnalysis_module.cc | 37 ++++++++++--------- sbndcode/CRT/CRTAna/run_crtana_data.fcl | 5 +-- .../CRTAna/run_crtana_data_tpc_matching.fcl | 6 +++ 3 files changed, 26 insertions(+), 22 deletions(-) create mode 100644 sbndcode/CRT/CRTAna/run_crtana_data_tpc_matching.fcl diff --git a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc index 793d69eef..d4b1e1816 100644 --- a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc +++ b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc @@ -1563,38 +1563,39 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: const art::Ptr spacepoint = tracksToSPMatches.at(track.key()); const art::Ptr crttrack = tracksToTrackMatches.at(track.key()); + if(spacepoint.isNonnull()) { - const anab::T0 spMatch = tracksToSPMatches.data(track.key()).ref(); + const anab::T0 spMatch = tracksToSPMatches.data(track.key()).ref(); - _tpc_sp_matched[nActualTracks] = true; - _tpc_sp_ts0[nActualTracks] = spacepoint->Ts0(); - _tpc_sp_ts1[nActualTracks] = spacepoint->Ts1(); - _tpc_sp_score[nActualTracks] = spMatch.TriggerConfidence(); + _tpc_sp_matched[nActualTracks] = true; + _tpc_sp_ts0[nActualTracks] = spacepoint->Ts0(); + _tpc_sp_ts1[nActualTracks] = spacepoint->Ts1(); + _tpc_sp_score[nActualTracks] = spMatch.TriggerConfidence(); } else { - _tpc_sp_matched[nActualTracks] = false; - _tpc_sp_ts0[nActualTracks] = -std::numeric_limits::max(); - _tpc_sp_ts1[nActualTracks] = -std::numeric_limits::max(); - _tpc_sp_score[nActualTracks] = -std::numeric_limits::max(); + _tpc_sp_matched[nActualTracks] = false; + _tpc_sp_ts0[nActualTracks] = -std::numeric_limits::max(); + _tpc_sp_ts1[nActualTracks] = -std::numeric_limits::max(); + _tpc_sp_score[nActualTracks] = -std::numeric_limits::max(); } if(crttrack.isNonnull()) { - const anab::T0 trackMatch = tracksToTrackMatches.data(track.key()).ref(); + const anab::T0 trackMatch = tracksToTrackMatches.data(track.key()).ref(); - _tpc_tr_matched[nActualTracks] = true; - _tpc_tr_ts0[nActualTracks] = crttrack->Ts0(); - _tpc_tr_ts1[nActualTracks] = crttrack->Ts1(); - _tpc_tr_score[nActualTracks] = trackMatch.TriggerConfidence(); + _tpc_tr_matched[nActualTracks] = true; + _tpc_tr_ts0[nActualTracks] = crttrack->Ts0(); + _tpc_tr_ts1[nActualTracks] = crttrack->Ts1(); + _tpc_tr_score[nActualTracks] = trackMatch.TriggerConfidence(); } else { - _tpc_tr_matched[nActualTracks] = false; - _tpc_tr_ts0[nActualTracks] = -std::numeric_limits::max(); - _tpc_tr_ts1[nActualTracks] = -std::numeric_limits::max(); - _tpc_tr_score[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_matched[nActualTracks] = false; + _tpc_tr_ts0[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_ts1[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_score[nActualTracks] = -std::numeric_limits::max(); } if(!fDataMode) diff --git a/sbndcode/CRT/CRTAna/run_crtana_data.fcl b/sbndcode/CRT/CRTAna/run_crtana_data.fcl index 77e03ee4e..b44faf920 100644 --- a/sbndcode/CRT/CRTAna/run_crtana_data.fcl +++ b/sbndcode/CRT/CRTAna/run_crtana_data.fcl @@ -10,9 +10,6 @@ services.CRTCalibService: @local::crt_calib_service physics.analyzers.crtana.FEBDataModuleLabel: "crtdecoder" physics.analyzers.crtana.DataMode: true -physics.analyzers.crtana.NoTPC: false +physics.analyzers.crtana.NoTPC: true physics.analyzers.crtana.HasPTB: true physics.analyzers.crtana.HasTDC: true -physics.analyzers.crtana.TPCTrackModuleLabel: "pandoraTrack" -physics.analyzers.crtana.CRTSpacePointMatchingModuleLabel: "crtspacepointmatching" -physics.analyzers.crtana.CRTTrackMatchingModuleLabel: "crttrackmatching" diff --git a/sbndcode/CRT/CRTAna/run_crtana_data_tpc_matching.fcl b/sbndcode/CRT/CRTAna/run_crtana_data_tpc_matching.fcl new file mode 100644 index 000000000..d1e2d09fe --- /dev/null +++ b/sbndcode/CRT/CRTAna/run_crtana_data_tpc_matching.fcl @@ -0,0 +1,6 @@ +#include "run_crtana_data.fcl" + +physics.analyzers.crtana.NoTPC: false +physics.analyzers.crtana.TPCTrackModuleLabel: "pandoraTrack" +physics.analyzers.crtana.CRTSpacePointMatchingModuleLabel: "crtspacepointmatching" +physics.analyzers.crtana.CRTTrackMatchingModuleLabel: "crttrackmatching" From df8af68e2f4bff0a4808d684dea2804eb04e1c5a Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Wed, 26 Feb 2025 06:40:53 -0600 Subject: [PATCH 04/18] Clean up event display for matching --- .../CRT/CRTEventDisplay/CRTEventDisplayAlg.cc | 171 ++++++++---------- .../CRT/CRTEventDisplay/CRTEventDisplayAlg.h | 18 +- .../crteventdisplayalg_sbnd.fcl | 18 +- 3 files changed, 100 insertions(+), 107 deletions(-) diff --git a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc index 842d68c59..07625fdd8 100644 --- a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc +++ b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc @@ -28,7 +28,8 @@ namespace sbnd::crt { fClusterLabel = config.ClusterLabel(); fSpacePointLabel = config.SpacePointLabel(); fTrackLabel = config.TrackLabel(); - fTrackMatchLabel = config.TrackMatchLabel(); + fTPCSpacePointMatchLabel = config.TPCSpacePointMatchLabel(); + fTPCTrackMatchLabel = config.TPCTrackMatchLabel(); fSaveRoot = config.SaveRoot(); fSaveViews = config.SaveViews(); @@ -45,6 +46,7 @@ namespace sbnd::crt { fDrawClusters = config.DrawClusters(); fDrawSpacePoints = config.DrawSpacePoints(); fDrawTracks = config.DrawTracks(); + fDrawTPCMatching = config.DrawTPCMatching(); fChoseTaggers = config.ChoseTaggers(); fChosenTaggers = config.ChosenTaggers(); @@ -64,6 +66,7 @@ namespace sbnd::crt { fClusterColourInterval = config.ClusterColourInterval(); fSpacePointColour = config.SpacePointColour(); fTrackColour = config.TrackColour(); + fTPCMatchColour = config.TPCMatchColour(); fUseTs0 = config.UseTs0(); fMinTime = config.MinTime(); @@ -469,54 +472,6 @@ namespace sbnd::crt { const art::Ptr spacepoint = spacePointVec[0]; const geo::Point_t pos = spacepoint->Pos(); const geo::Point_t err = spacepoint->Err(); - try { - auto spacePointsHandle = event.getValidHandle>(fSpacePointLabel); - art::FindOneP CRTSPstoTPCTracks(spacePointsHandle, event, "crtspacepointmatching"); - const art::Ptr TPCTrack = CRTSPstoTPCTracks.at(spacepoint.key()); - if(TPCTrack.isNonnull()) { - const anab::T0 t0Match = CRTSPstoTPCTracks.data(spacepoint.key()).ref(); - double t0SPMatchConfidence = t0Match.TriggerConfidence(); - const geo::Point_t startTPC = TPCTrack->Start(); - const geo::Vector_t dirTPC = TPCTrack->StartDirection(); - TPolyLine3D *lineTPC = new TPolyLine3D(2); - geo::Point_t aTPC {0,0,0}; - geo::Point_t bTPC {0,0,0}; - int i = 0; - do - { - aTPC = startTPC + i * dirTPC; - ++i; - } - while(IsPointInsideBox(crtLims, aTPC)); - i = 0; - do - { - bTPC = startTPC + i * dirTPC; - --i; - } - while(IsPointInsideBox(crtLims, bTPC)); - lineTPC->SetPoint(0, aTPC.X(), aTPC.Y(), aTPC.Z()); - lineTPC->SetPoint(1, bTPC.X(), bTPC.Y(), bTPC.Z()); - lineTPC->SetLineColor(3); - lineTPC->SetLineWidth(fLineWidth); - lineTPC->Draw(); - TPaveText *pt = new TPaveText(0.05,0.85,0.35,0.65,"NB"); - pt->SetTextSize(0.02); - pt->SetFillStyle(0); - pt->SetLineStyle(0); - pt->SetTextAlign(12); - pt->SetBorderSize(0); - pt->Draw(); - if(fPrint) - std::cout << "t0 matching confidence for this spacepoint is " << t0SPMatchConfidence << std::endl; - } - else { - continue; - } - } catch(...) { - continue; - } - double rmin[3] = {pos.X() - err.X(), pos.Y() - err.Y(), pos.Z() - err.Z()}; double rmax[3] = {pos.X() + err.X(), pos.Y() + err.Y(), pos.Z() + err.Z()}; @@ -530,8 +485,43 @@ namespace sbnd::crt { << spacepoint->Ts0() << " (" << spacepoint->Ts0() - G4RefTime << ") or t1 = " << spacepoint->Ts1() << " (" << spacepoint->Ts1() - G4RefTime << ")" << " with PE " << spacepoint->PE() - << std::endl; - } + << std::endl; + + if(fDrawTPCMatching) + { + auto spacePointsHandle = event.getValidHandle>(fSpacePointLabel); + art::FindOneP CRTSPstoTPCTracks(spacePointsHandle, event, fTPCSpacePointMatchLabel); + + const art::Ptr TPCTrack = CRTSPstoTPCTracks.at(spacepoint.key()); + + if(TPCTrack.isNonnull()) + { + const anab::T0 t0Match = CRTSPstoTPCTracks.data(spacepoint.key()).ref(); + + const double t0MatchConfidence = t0Match.TriggerConfidence(); + const geo::Point_t startTPC = TPCTrack->Start(); + const geo::Point_t endTPC = TPCTrack->End(); + + TPolyLine3D *lineTPC = new TPolyLine3D(2); + lineTPC->SetPoint(0, startTPC.X(), startTPC.Y(), startTPC.Z()); + lineTPC->SetPoint(1, endTPC.X(), endTPC.Y(), endTPC.Z()); + lineTPC->SetLineColor(fTPCMatchColour); + lineTPC->SetLineWidth(fLineWidth); + lineTPC->Draw(); + + TPaveText *pt = new TPaveText(0.05,0.85,0.35,0.65,"NB"); + pt->SetTextSize(0.02); + pt->SetFillStyle(0); + pt->SetLineStyle(0); + pt->SetTextAlign(12); + pt->SetBorderSize(0); + pt->Draw(); + + if(fPrint) + std::cout << "t0 matching confidence for this spacepoint is " << t0MatchConfidence << std::endl; + } + } + } else if(spacePointVec.size() != 0) std::cout << "What an earth is going on here then..." << std::endl; } @@ -543,8 +533,6 @@ namespace sbnd::crt { auto tracksHandle = event.getValidHandle>(fTrackLabel); std::vector> tracksVec; art::fill_ptr_vector(tracksVec, tracksHandle); - art::FindOneP CRTTrackstoTPCTracks(tracksHandle, event, "crttrackmatching"); - for(auto track : tracksVec) { @@ -567,49 +555,6 @@ namespace sbnd::crt { const geo::Point_t start = track->Start(); const geo::Vector_t dir = track->Direction(); - const art::Ptr TPCTrack = CRTTrackstoTPCTracks.at(track.key()); - if(TPCTrack.isNonnull()) { - const anab::T0 t0Match = CRTTrackstoTPCTracks.data(track.key()).ref(); - double t0MatchConfidence = t0Match.TriggerConfidence(); - - const geo::Point_t startTPC = TPCTrack->Start(); - const geo::Vector_t dirTPC = TPCTrack->StartDirection(); - TPolyLine3D *lineTPC = new TPolyLine3D(2); - geo::Point_t aTPC {0,0,0}; - geo::Point_t bTPC {0,0,0}; - int i = 0; - do - { - aTPC = startTPC + i * dirTPC; - ++i; - } - while(IsPointInsideBox(crtLims, aTPC)); - i = 0; - do - { - bTPC = startTPC + i * dirTPC; - --i; - } - while(IsPointInsideBox(crtLims, bTPC)); - lineTPC->SetPoint(0, aTPC.X(), aTPC.Y(), aTPC.Z()); - lineTPC->SetPoint(1, bTPC.X(), bTPC.Y(), bTPC.Z()); - lineTPC->SetLineColor(2); - lineTPC->SetLineWidth(fLineWidth); - lineTPC->Draw(); - TPaveText *pt = new TPaveText(0.05,0.8,0.35,0.6,"NB"); - pt->SetTextSize(0.02); - pt->SetFillStyle(0); - pt->SetLineStyle(0); - pt->SetTextAlign(12); - pt->SetBorderSize(0); - pt->Draw(); - if(fPrint) - std::cout << "t0 matching confidence for this track is " << t0MatchConfidence << std::endl; - } - else { - continue; - } - TPolyLine3D *line = new TPolyLine3D(2); geo::Point_t a {0,0,0}; @@ -647,7 +592,39 @@ namespace sbnd::crt { << "\tat ts1 " << track->Ts1() << " (" << track->Ts1() - G4RefTime << ")\n" << "\tfrom three hits? " << track->Triple() << std::endl; - } + if(fDrawTPCMatching) + { + art::FindOneP CRTTrackstoTPCTracks(tracksHandle, event, fTPCTrackMatchLabel); + const art::Ptr TPCTrack = CRTTrackstoTPCTracks.at(track.key()); + + if(TPCTrack.isNonnull()) + { + const anab::T0 t0Match = CRTTrackstoTPCTracks.data(track.key()).ref(); + + const double t0MatchConfidence = t0Match.TriggerConfidence(); + const geo::Point_t startTPC = TPCTrack->Start(); + const geo::Point_t endTPC = TPCTrack->End(); + + TPolyLine3D *lineTPC = new TPolyLine3D(2); + lineTPC->SetPoint(0, startTPC.X(), startTPC.Y(), startTPC.Z()); + lineTPC->SetPoint(1, endTPC.X(), endTPC.Y(), endTPC.Z()); + lineTPC->SetLineColor(fTPCMatchColour); + lineTPC->SetLineWidth(fLineWidth); + lineTPC->Draw(); + + TPaveText *pt = new TPaveText(0.05,0.8,0.35,0.6,"NB"); + pt->SetTextSize(0.02); + pt->SetFillStyle(0); + pt->SetLineStyle(0); + pt->SetTextAlign(12); + pt->SetBorderSize(0); + pt->Draw(); + + if(fPrint) + std::cout << "t0 matching confidence for this track is " << t0MatchConfidence << std::endl; + } + } + } } if(fSaveRoot) diff --git a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h index b1875fb6e..bdc0ed9ed 100644 --- a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h +++ b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h @@ -87,8 +87,11 @@ namespace sbnd::crt { fhicl::Atom TrackLabel { Name("TrackLabel") }; - fhicl::Atom TrackMatchLabel { - Name("TrackMatchLabel") + fhicl::Atom TPCSpacePointMatchLabel { + Name("TPCSpacePointMatchLabel") + }; + fhicl::Atom TPCTrackMatchLabel { + Name("TPCTrackMatchLabel") }; fhicl::Atom SaveRoot { @@ -134,6 +137,9 @@ namespace sbnd::crt { fhicl::Atom DrawTracks { Name("DrawTracks") }; + fhicl::Atom DrawTPCMatching { + Name("DrawTPCMatching") + }; fhicl::Atom ChoseTaggers { Name("ChoseTaggers") @@ -185,6 +191,9 @@ namespace sbnd::crt { fhicl::Atom TrackColour { Name("TrackColour") }; + fhicl::Atom TPCMatchColour { + Name("TPCMatchColour") + }; fhicl::Atom UseTs0 { Name ("UseTs0") @@ -250,7 +259,8 @@ namespace sbnd::crt { art::InputTag fClusterLabel; art::InputTag fSpacePointLabel; art::InputTag fTrackLabel; - art::InputTag fTrackMatchLabel; + art::InputTag fTPCSpacePointMatchLabel; + art::InputTag fTPCTrackMatchLabel; bool fSaveRoot; bool fSaveViews; @@ -267,6 +277,7 @@ namespace sbnd::crt { bool fDrawClusters; bool fDrawSpacePoints; bool fDrawTracks; + bool fDrawTPCMatching; bool fChoseTaggers; std::vector fChosenTaggers; @@ -286,6 +297,7 @@ namespace sbnd::crt { int fClusterColourInterval; int fSpacePointColour; int fTrackColour; + int fTPCMatchColour; bool fUseTs0; double fMinTime; diff --git a/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl b/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl index f73721397..ffb78aa37 100644 --- a/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl +++ b/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl @@ -11,13 +11,14 @@ crteventdisplayalg_sbnd: MC: true - SimLabel: "largeant" - SimDepositLabel: "genericcrt" - StripHitLabel: "crtstrips" - ClusterLabel: "crtclustering" - SpacePointLabel: "crtspacepoints" - TrackLabel: "crttracks" - TrackMatchLabel: "crttrackmatching" + SimLabel: "largeant" + SimDepositLabel: "genericcrt" + StripHitLabel: "crtstrips" + ClusterLabel: "crtclustering" + SpacePointLabel: "crtspacepoints" + TrackLabel: "crttracks" + TPCTrackMatchLabel: "crttrackmatching" + TPCSpacePointMatchLabel: "crtspacepointmatching" SaveRoot: true SaveViews: false @@ -34,6 +35,7 @@ crteventdisplayalg_sbnd: DrawClusters: true DrawSpacePoints: true DrawTracks: true + DrawTPCMatching: false ## If chose taggers is set to true then the vector ChosenTaggers ## is used to determine which taggers are drawn. The integer values @@ -64,6 +66,7 @@ crteventdisplayalg_sbnd: ClusterColourInterval: 20 SpacePointColour: 6 TrackColour: 9 + TPCMatchColour: 8 UseTs0: false MinTime: -10000000000 @@ -82,6 +85,7 @@ crteventdisplayalg_sbnd_data.DrawStripHits: false crteventdisplayalg_sbnd_data.DrawClusters: false crteventdisplayalg_sbnd_data.DrawSpacePoints: true crteventdisplayalg_sbnd_data.DrawTracks: true +crteventdisplayalg_sbnd_data.DrawTPCMatching: true crteventdisplayalg_sbnd_data.MC: false crteventdisplayalg_sbnd_data.DrawModules: true crteventdisplayalg_sbnd_data.DrawTrueTracks: false From b2750b0ece38c828a69d14b77219ccb7e33f2815 Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Wed, 26 Feb 2025 08:18:55 -0600 Subject: [PATCH 05/18] Port changes to new fcl to preserve MC behaviour --- .../crttpcmatchingproducers_sbnd.fcl | 16 ++++++++++++++-- .../CRTTPCMatching/run_crttpcmatching_data.fcl | 6 ++++++ 2 files changed, 20 insertions(+), 2 deletions(-) create mode 100644 sbndcode/CRT/CRTTPCMatching/run_crttpcmatching_data.fcl diff --git a/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl b/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl index 1144a0df9..d1a044e0e 100644 --- a/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl +++ b/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl @@ -12,9 +12,12 @@ crtspacepointmatchalg_sbnd: MaxUncert: 20. TPCTrackLabel: "pandoraTrack" CRTSpacePointLabel: "crtspacepoints" - UseTs0: true + UseTs0: false } +crtspacepointmatchalg_data_sbnd: @local::crtspacepointmatchalg_sbnd +crtspacepointmatchalg_data_sbnd.UseTs0: true + crtspacepointmatchproducer_sbnd: { CRTSpacePointModuleLabel: "crtspacepoints" @@ -24,6 +27,9 @@ crtspacepointmatchproducer_sbnd: module_type: "CRTSpacePointMatching" } +crtspacepointmatchproducer_data_sbnd: @local::crtspacepointmatchproducer_sbnd +crtspacepointmatchproducer_data_sbnd.MatchingAlg: @local::crtspacepointmatchalg_data_sbnd + crttrackmatchalg_sbnd: { MaxAngleDiff: 0.4 @@ -32,9 +38,12 @@ crttrackmatchalg_sbnd: MinTPCTrackLength: 5. SelectionMetric: "score" TPCTrackLabel: "pandoraTrack" - UseTs0: true + UseTs0: false } +crttrackmatchalg_data_sbnd: @local::crttrackmatchalg_sbnd +crttrackmatchalg_data_sbnd.UseTs0: true + crttrackmatchproducer_sbnd: { CRTTrackModuleLabel: "crttracks" @@ -44,4 +53,7 @@ crttrackmatchproducer_sbnd: module_type: "CRTTrackMatching" } +crttrackmatchproducer_data_sbnd: @local::crttrackmatchproducer_sbnd +crttrackmatchproducer_data_sbnd.MatchingAlg: @local::crttrackmatchalg_sbnd + END_PROLOG diff --git a/sbndcode/CRT/CRTTPCMatching/run_crttpcmatching_data.fcl b/sbndcode/CRT/CRTTPCMatching/run_crttpcmatching_data.fcl new file mode 100644 index 000000000..d1723f5f6 --- /dev/null +++ b/sbndcode/CRT/CRTTPCMatching/run_crttpcmatching_data.fcl @@ -0,0 +1,6 @@ +#include "run_crttpcmatching.fcl" + +process_name: CRTTPCMatching + +physics.producers.crtspacepointmatching: @local::crtspacepointmatchproducer_data_sbnd +physics.producers.crttrackmatching: @local::crttrackmatchproducer_data_sbnd From 838ec3518a1525571c3151b468e250933005e319 Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Thu, 27 Feb 2025 08:17:42 -0600 Subject: [PATCH 06/18] More updates for visualising matching --- .../CRT/CRTEventDisplay/CRTEventDisplayAlg.cc | 378 ++++++++++++------ .../CRT/CRTEventDisplay/CRTEventDisplayAlg.h | 17 + .../crteventdisplayalg_sbnd.fcl | 38 +- 3 files changed, 284 insertions(+), 149 deletions(-) diff --git a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc index 07625fdd8..0e08c638c 100644 --- a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc +++ b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc @@ -30,6 +30,7 @@ namespace sbnd::crt { fTrackLabel = config.TrackLabel(); fTPCSpacePointMatchLabel = config.TPCSpacePointMatchLabel(); fTPCTrackMatchLabel = config.TPCTrackMatchLabel(); + fTPCTrackLabel = config.TPCTrackLabel(); fSaveRoot = config.SaveRoot(); fSaveViews = config.SaveViews(); @@ -47,6 +48,8 @@ namespace sbnd::crt { fDrawSpacePoints = config.DrawSpacePoints(); fDrawTracks = config.DrawTracks(); fDrawTPCMatching = config.DrawTPCMatching(); + fOnlyDrawMatched = config.OnlyDrawMatched(); + fDisplayMatchScore = config.DisplayMatchScore(); fChoseTaggers = config.ChoseTaggers(); fChosenTaggers = config.ChosenTaggers(); @@ -109,6 +112,16 @@ namespace sbnd::crt { fDrawClusters = tf; } + void CRTEventDisplayAlg::SetMinTime(double time) + { + fMinTime = time; + } + + void CRTEventDisplayAlg::SetMaxTime(double time) + { + fMaxTime = time; + } + void CRTEventDisplayAlg::SetPrint(bool tf) { fPrint = tf; @@ -151,6 +164,9 @@ namespace sbnd::crt { void CRTEventDisplayAlg::Draw(detinfo::DetectorClocksData const& clockData, const art::Event& event, const TString& saveName) { + auto geometryService = lar::providerFrom(); + auto const detProp = art::ServiceHandle()->DataFor(event); + if(fMC) fCRTBackTrackerAlg.SetupMaps(event); @@ -161,6 +177,7 @@ namespace sbnd::crt { TCanvas *c1 = new TCanvas("c1","",700,700); std::vector crtLims = fCRTGeoAlg.CRTLimits(); + std::vector crtLims2 = fCRTGeoAlg.CRTLimits(); crtLims[0] -= 100; crtLims[1] -= 100; crtLims[2] -= 100; crtLims[3] += 100; crtLims[4] += 100; crtLims[5] += 100; @@ -172,12 +189,8 @@ namespace sbnd::crt { if(fChoseTaggers && std::find(fChosenTaggers.begin(), fChosenTaggers.end(), CRTCommonUtils::GetTaggerEnum(name)) == fChosenTaggers.end()) continue; - double rmin[3] = {tagger.minX, - tagger.minY, - tagger.minZ}; - double rmax[3] = {tagger.maxX, - tagger.maxY, - tagger.maxZ}; + double rmin[3] = {tagger.minX, tagger.minY, tagger.minZ}; + double rmax[3] = {tagger.maxX, tagger.maxY, tagger.maxZ}; DrawCube(c1, rmin, rmax, fTaggerColour); } @@ -191,12 +204,8 @@ namespace sbnd::crt { if(fChoseTaggers && std::find(fChosenTaggers.begin(), fChosenTaggers.end(), CRTCommonUtils::GetTaggerEnum(module.taggerName)) == fChosenTaggers.end()) continue; - double rmin[3] = {module.minX, - module.minY, - module.minZ}; - double rmax[3] = {module.maxX, - module.maxY, - module.maxZ}; + double rmin[3] = {module.minX, module.minY, module.minZ}; + double rmax[3] = {module.maxX, module.maxY, module.maxZ}; if(fHighlightModules && std::find(fHighlightedModules.begin(), fHighlightedModules.end(), module.adID) == fHighlightedModules.end()) { @@ -212,13 +221,8 @@ namespace sbnd::crt { { const std::array febPos = fCRTGeoAlg.FEBWorldPos(module); - double rminFEB[3] = {febPos[0], - febPos[2], - febPos[4]}; - - double rmaxFEB[3] = {febPos[1], - febPos[3], - febPos[5]}; + double rminFEB[3] = {febPos[0], febPos[2], febPos[4]}; + double rmaxFEB[3] = {febPos[1], febPos[3], febPos[5]}; DrawCube(c1, rminFEB, rmaxFEB, fFEBColour); @@ -226,13 +230,8 @@ namespace sbnd::crt { { const std::array febCh0Pos = fCRTGeoAlg.FEBChannel0WorldPos(module); - double rminCh0[3] = {febCh0Pos[0], - febCh0Pos[2], - febCh0Pos[4]}; - - double rmaxCh0[3] = {febCh0Pos[1], - febCh0Pos[3], - febCh0Pos[5]}; + double rminCh0[3] = {febCh0Pos[0], febCh0Pos[2], febCh0Pos[4]}; + double rmaxCh0[3] = {febCh0Pos[1], febCh0Pos[3], febCh0Pos[5]}; DrawCube(c1, rminCh0, rmaxCh0, fFEBEndColour); } @@ -248,12 +247,8 @@ namespace sbnd::crt { if(fChoseTaggers && std::find(fChosenTaggers.begin(), fChosenTaggers.end(), fCRTGeoAlg.ChannelToTaggerEnum(strip.channel0)) == fChosenTaggers.end()) continue; - double rmin[3] = {strip.minX, - strip.minY, - strip.minZ}; - double rmax[3] = {strip.maxX, - strip.maxY, - strip.maxZ}; + double rmin[3] = {strip.minX, strip.minY, strip.minZ}; + double rmax[3] = {strip.maxX, strip.maxY, strip.maxZ}; DrawCube(c1, rmin, rmax, fTaggerColour, 1); } @@ -262,23 +257,15 @@ namespace sbnd::crt { // Draw the TPC with central cathode if(fDrawTPC) { - double rmin[3] = {fTPCGeoAlg.MinX(), - fTPCGeoAlg.MinY(), - fTPCGeoAlg.MinZ()}; - double rmax[3] = {-fTPCGeoAlg.CpaWidth(), - fTPCGeoAlg.MaxY(), - fTPCGeoAlg.MaxZ()}; + double rmin[3] = {fTPCGeoAlg.MinX(), fTPCGeoAlg.MinY(), fTPCGeoAlg.MinZ()}; + double rmax[3] = {-fTPCGeoAlg.CpaWidth(), fTPCGeoAlg.MaxY(), fTPCGeoAlg.MaxZ()}; DrawCube(c1, rmin, rmax, fTPCColour); - double rmin2[3] = {fTPCGeoAlg.CpaWidth(), - fTPCGeoAlg.MinY(), - fTPCGeoAlg.MinZ()}; - double rmax2[3] = {fTPCGeoAlg.MaxX(), - fTPCGeoAlg.MaxY(), - fTPCGeoAlg.MaxZ()}; + double rmin2[3] = {fTPCGeoAlg.CpaWidth(), fTPCGeoAlg.MinY(), fTPCGeoAlg.MinZ()}; + double rmax2[3] = {fTPCGeoAlg.MaxX(), fTPCGeoAlg.MaxY(), fTPCGeoAlg.MaxZ()}; DrawCube(c1, rmin2, rmax2, fTPCColour); } - + // Draw true track trajectories for visible particles that cross the CRT if(fDrawTrueTracks) { @@ -470,58 +457,121 @@ namespace sbnd::crt { if(spacePointVec.size() == 1) { const art::Ptr spacepoint = spacePointVec[0]; + const double spacepointTime = fUseTs0 ? spacepoint->Ts0() - G4RefTime : spacepoint->Ts1() - G4RefTime; + + double matchScore = -999.f; + bool foundMatch = false; + + if(fDrawTPCMatching) + { + auto spacePointsHandle = event.getValidHandle>(fSpacePointLabel); + auto tpcTracksHandle = event.getValidHandle>(fTPCTrackLabel); + art::FindOneP CRTSPstoTPCTracks(spacePointsHandle, event, fTPCSpacePointMatchLabel); + art::FindManyP TPCTracksToHits(tpcTracksHandle, event, fTPCTrackLabel); + + const art::Ptr TPCTrack = CRTSPstoTPCTracks.at(spacepoint.key()); + + if(TPCTrack.isNonnull()) + { + const std::vector> hits = TPCTracksToHits.at(TPCTrack.key()); + const int driftDirection = TPCGeoUtil::DriftDirectionFromHits(geometryService, hits); + const double shift = driftDirection * spacepointTime * detProp.DriftVelocity() * 1e-3; + + foundMatch = true; + + const anab::T0 t0Match = CRTSPstoTPCTracks.data(spacepoint.key()).ref(); + + matchScore = t0Match.TriggerConfidence(); + geo::Point_t startTPC = TPCTrack->Start(); + geo::Point_t endTPC = TPCTrack->End(); + const geo::Vector_t startDirTPC = TPCTrack->StartDirection(); + const geo::Vector_t endDirTPC = TPCTrack->EndDirection(); + + startTPC.SetX(startTPC.X() + shift); + endTPC.SetX(endTPC.X() + shift); + + TPolyLine3D *lineTPC = new TPolyLine3D(2); + lineTPC->SetPoint(0, startTPC.X(), startTPC.Y(), startTPC.Z()); + lineTPC->SetPoint(1, endTPC.X(), endTPC.Y(), endTPC.Z()); + lineTPC->SetLineColor(fTPCMatchColour); + lineTPC->SetLineWidth(fLineWidth); + lineTPC->Draw(); + + TPolyLine3D *lineTPCDashStart = new TPolyLine3D(2); + geo::Point_t intersectionStart {0,0,0}; + + int i = 0; + do + { + intersectionStart = startTPC + i * startDirTPC; + --i; + } + while(IsPointInsideBox(crtLims2, intersectionStart)); + + lineTPCDashStart->SetPoint(0, startTPC.X(), startTPC.Y(), startTPC.Z()); + lineTPCDashStart->SetPoint(1, intersectionStart.X(), intersectionStart.Y(), intersectionStart.Z()); + lineTPCDashStart->SetLineStyle(2); + lineTPCDashStart->SetLineColor(fTPCMatchColour); + lineTPCDashStart->SetLineWidth(fLineWidth - 1); + lineTPCDashStart->Draw(); + + TPolyLine3D *lineTPCDashEnd = new TPolyLine3D(2); + geo::Point_t intersectionEnd {0,0,0}; + + i = 0; + do + { + intersectionEnd = endTPC + i * endDirTPC; + ++i; + } + while(IsPointInsideBox(crtLims2, intersectionEnd)); + + lineTPCDashEnd->SetPoint(0, endTPC.X(), endTPC.Y(), endTPC.Z()); + lineTPCDashEnd->SetPoint(1, intersectionEnd.X(), intersectionEnd.Y(), intersectionEnd.Z()); + lineTPCDashEnd->SetLineStyle(2); + lineTPCDashEnd->SetLineColor(fTPCMatchColour); + lineTPCDashEnd->SetLineWidth(fLineWidth - 1); + lineTPCDashEnd->Draw(); + + if(fDisplayMatchScore) + { + TPaveText *pt = new TPaveText(0.4,0.86,0.85,0.9, "NDC"); + pt->SetTextSize(0.02); + pt->SetFillStyle(0); + pt->SetLineStyle(0); + pt->SetTextAlign(12); + pt->SetBorderSize(0); + pt->AddText(Form("SP Matching Score = %g", matchScore)); + pt->AddText(Form("SP Shift = %g cm (t = %g ns)", shift, spacepointTime)); + pt->Draw(); + } + } + } + const geo::Point_t pos = spacepoint->Pos(); const geo::Point_t err = spacepoint->Err(); double rmin[3] = {pos.X() - err.X(), pos.Y() - err.Y(), pos.Z() - err.Z()}; double rmax[3] = {pos.X() + err.X(), pos.Y() + err.Y(), pos.Z() + err.Z()}; - DrawCube(c1, rmin, rmax, fSpacePointColour); + if((fOnlyDrawMatched && foundMatch) || !fOnlyDrawMatched) + DrawCube(c1, rmin, rmax, fSpacePointColour); if(fPrint) - std::cout << "Space Point: (" - << rmin[0] << ", " << rmin[1] << ", " << rmin[2] << ") --> (" - << rmax[0] << ", " << rmax[1] << ", " << rmax[2] << ") at t0 = " - << spacepoint->Ts0() << " (" << spacepoint->Ts0() - G4RefTime << ") or t1 = " - << spacepoint->Ts1() << " (" << spacepoint->Ts1() - G4RefTime << ")" - << " with PE " << spacepoint->PE() - << std::endl; - - if(fDrawTPCMatching) - { - auto spacePointsHandle = event.getValidHandle>(fSpacePointLabel); - art::FindOneP CRTSPstoTPCTracks(spacePointsHandle, event, fTPCSpacePointMatchLabel); - - const art::Ptr TPCTrack = CRTSPstoTPCTracks.at(spacepoint.key()); - - if(TPCTrack.isNonnull()) - { - const anab::T0 t0Match = CRTSPstoTPCTracks.data(spacepoint.key()).ref(); - - const double t0MatchConfidence = t0Match.TriggerConfidence(); - const geo::Point_t startTPC = TPCTrack->Start(); - const geo::Point_t endTPC = TPCTrack->End(); - - TPolyLine3D *lineTPC = new TPolyLine3D(2); - lineTPC->SetPoint(0, startTPC.X(), startTPC.Y(), startTPC.Z()); - lineTPC->SetPoint(1, endTPC.X(), endTPC.Y(), endTPC.Z()); - lineTPC->SetLineColor(fTPCMatchColour); - lineTPC->SetLineWidth(fLineWidth); - lineTPC->Draw(); - - TPaveText *pt = new TPaveText(0.05,0.85,0.35,0.65,"NB"); - pt->SetTextSize(0.02); - pt->SetFillStyle(0); - pt->SetLineStyle(0); - pt->SetTextAlign(12); - pt->SetBorderSize(0); - pt->Draw(); - - if(fPrint) - std::cout << "t0 matching confidence for this spacepoint is " << t0MatchConfidence << std::endl; - } - } - } + { + std::cout << "Space Point: (" + << rmin[0] << ", " << rmin[1] << ", " << rmin[2] << ") --> (" + << rmax[0] << ", " << rmax[1] << ", " << rmax[2] << ") at t0 = " + << spacepoint->Ts0() << " (" << spacepoint->Ts0() - G4RefTime << ") or t1 = " + << spacepoint->Ts1() << " (" << spacepoint->Ts1() - G4RefTime << ")" + << " with PE " << spacepoint->PE(); + + if(foundMatch) + std::cout << " and found TPC track match with score: " << matchScore << std::endl; + else + std::cout << std::endl; + } + } else if(spacePointVec.size() != 0) std::cout << "What an earth is going on here then..." << std::endl; } @@ -543,7 +593,7 @@ namespace sbnd::crt { std::set taggers = track->Taggers(); - bool none = true; + bool none = fChoseTaggers; for(auto const& tagger : taggers) { if(fChoseTaggers && std::find(fChosenTaggers.begin(), fChosenTaggers.end(), tagger) != fChosenTaggers.end()) @@ -553,6 +603,94 @@ namespace sbnd::crt { if(none) continue; + double matchScore = -999.f; + bool foundMatch = false; + + if(fDrawTPCMatching) + { + art::FindOneP CRTTrackstoTPCTracks(tracksHandle, event, fTPCTrackMatchLabel); + const art::Ptr TPCTrack = CRTTrackstoTPCTracks.at(track.key()); + + auto tpcTracksHandle = event.getValidHandle>(fTPCTrackLabel); + art::FindManyP TPCTracksToHits(tpcTracksHandle, event, fTPCTrackLabel); + + if(TPCTrack.isNonnull()) + { + const std::vector> hits = TPCTracksToHits.at(TPCTrack.key()); + const int driftDirection = TPCGeoUtil::DriftDirectionFromHits(geometryService, hits); + const double shift = driftDirection * trackTime * detProp.DriftVelocity() * 1e-3; + + foundMatch = true; + + const anab::T0 t0Match = CRTTrackstoTPCTracks.data(track.key()).ref(); + + matchScore = t0Match.TriggerConfidence(); + geo::Point_t startTPC = TPCTrack->Start(); + geo::Point_t endTPC = TPCTrack->End(); + const geo::Vector_t startDirTPC = TPCTrack->StartDirection(); + const geo::Vector_t endDirTPC = TPCTrack->EndDirection(); + + startTPC.SetX(startTPC.X() + shift); + endTPC.SetX(endTPC.X() + shift); + + TPolyLine3D *lineTPC = new TPolyLine3D(2); + lineTPC->SetPoint(0, startTPC.X(), startTPC.Y(), startTPC.Z()); + lineTPC->SetPoint(1, endTPC.X(), endTPC.Y(), endTPC.Z()); + lineTPC->SetLineColor(fTPCMatchColour); + lineTPC->SetLineWidth(fLineWidth); + lineTPC->Draw(); + + TPolyLine3D *lineTPCDashStart = new TPolyLine3D(2); + geo::Point_t intersectionStart {0,0,0}; + + int i = 0; + do + { + intersectionStart = startTPC + i * startDirTPC; + --i; + } + while(IsPointInsideBox(crtLims2, intersectionStart)); + + lineTPCDashStart->SetPoint(0, startTPC.X(), startTPC.Y(), startTPC.Z()); + lineTPCDashStart->SetPoint(1, intersectionStart.X(), intersectionStart.Y(), intersectionStart.Z()); + lineTPCDashStart->SetLineStyle(2); + lineTPCDashStart->SetLineColor(fTPCMatchColour); + lineTPCDashStart->SetLineWidth(fLineWidth - 1); + lineTPCDashStart->Draw(); + + TPolyLine3D *lineTPCDashEnd = new TPolyLine3D(2); + geo::Point_t intersectionEnd {0,0,0}; + + i = 0; + do + { + intersectionEnd = endTPC + i * endDirTPC; + ++i; + } + while(IsPointInsideBox(crtLims2, intersectionEnd)); + + lineTPCDashEnd->SetPoint(0, endTPC.X(), endTPC.Y(), endTPC.Z()); + lineTPCDashEnd->SetPoint(1, intersectionEnd.X(), intersectionEnd.Y(), intersectionEnd.Z()); + lineTPCDashEnd->SetLineStyle(2); + lineTPCDashEnd->SetLineColor(fTPCMatchColour); + lineTPCDashEnd->SetLineWidth(fLineWidth - 1); + lineTPCDashEnd->Draw(); + + if(fDisplayMatchScore) + { + TPaveText *pt = new TPaveText(0.4,0.9,0.85,0.94, "NDC"); + pt->SetTextSize(0.02); + pt->SetFillStyle(0); + pt->SetLineStyle(0); + pt->SetTextAlign(12); + pt->SetBorderSize(0); + pt->AddText(Form("Track Matching Score = %g", matchScore)); + pt->AddText(Form("Track Shift = %g cm (t = %g ns)", shift, trackTime)); + pt->Draw(); + } + } + } + const geo::Point_t start = track->Start(); const geo::Vector_t dir = track->Direction(); @@ -581,50 +719,26 @@ namespace sbnd::crt { line->SetLineColor(fTrackColour); line->SetLineWidth(fLineWidth); - line->Draw(); + + if((fOnlyDrawMatched && foundMatch) || !fOnlyDrawMatched) + line->Draw(); if(fPrint) - std::cout << "Track at (" << start.X() << ", " << start.Y() << ", " << start.Z() << ")\n" - << "\twith direction (" << dir.X() << ", " << dir.Y() << ", " << dir.Z() << ")\n" - << "\tdrawn between (" << a.X() << ", " << a.Y() << ", " << a.Z() << ")\n" - << "\tand (" << b.X() << ", " << b.Y() << ", " << b.Z() << ")\n" - << "\tat ts0 " << track->Ts0() << " (" << track->Ts0() - G4RefTime << ")\n" - << "\tat ts1 " << track->Ts1() << " (" << track->Ts1() - G4RefTime << ")\n" - << "\tfrom three hits? " << track->Triple() << std::endl; - - if(fDrawTPCMatching) - { - art::FindOneP CRTTrackstoTPCTracks(tracksHandle, event, fTPCTrackMatchLabel); - const art::Ptr TPCTrack = CRTTrackstoTPCTracks.at(track.key()); - - if(TPCTrack.isNonnull()) - { - const anab::T0 t0Match = CRTTrackstoTPCTracks.data(track.key()).ref(); - - const double t0MatchConfidence = t0Match.TriggerConfidence(); - const geo::Point_t startTPC = TPCTrack->Start(); - const geo::Point_t endTPC = TPCTrack->End(); - - TPolyLine3D *lineTPC = new TPolyLine3D(2); - lineTPC->SetPoint(0, startTPC.X(), startTPC.Y(), startTPC.Z()); - lineTPC->SetPoint(1, endTPC.X(), endTPC.Y(), endTPC.Z()); - lineTPC->SetLineColor(fTPCMatchColour); - lineTPC->SetLineWidth(fLineWidth); - lineTPC->Draw(); - - TPaveText *pt = new TPaveText(0.05,0.8,0.35,0.6,"NB"); - pt->SetTextSize(0.02); - pt->SetFillStyle(0); - pt->SetLineStyle(0); - pt->SetTextAlign(12); - pt->SetBorderSize(0); - pt->Draw(); - - if(fPrint) - std::cout << "t0 matching confidence for this track is " << t0MatchConfidence << std::endl; - } - } - } + { + std::cout << "Track at (" << start.X() << ", " << start.Y() << ", " << start.Z() << ")\n" + << "\twith direction (" << dir.X() << ", " << dir.Y() << ", " << dir.Z() << ")\n" + << "\tdrawn between (" << a.X() << ", " << a.Y() << ", " << a.Z() << ")\n" + << "\tand (" << b.X() << ", " << b.Y() << ", " << b.Z() << ")\n" + << "\tat ts0 " << track->Ts0() << " (" << track->Ts0() - G4RefTime << ")\n" + << "\tat ts1 " << track->Ts1() << " (" << track->Ts1() - G4RefTime << ")\n" + << "\tfrom three hits? " << track->Triple(); + + if(foundMatch) + std::cout << " and found TPC track match with score: " << matchScore << std::endl; + else + std::cout << std::endl; + } + } } if(fSaveRoot) diff --git a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h index bdc0ed9ed..971582960 100644 --- a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h +++ b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h @@ -26,6 +26,7 @@ //larsoft #include "lardataalg/DetectorInfo/DetectorClocksData.h" +#include "lardata/DetectorInfoServices/DetectorPropertiesService.h" // sbnobj #include "sbnobj/SBND/CRT/CRTStripHit.hh" @@ -37,6 +38,7 @@ #include "sbndcode/Geometry/GeometryWrappers/TPCGeoAlg.h" #include "sbndcode/Geometry/GeometryWrappers/CRTGeoAlg.h" #include "sbndcode/CRT/CRTBackTracker/CRTBackTrackerAlg.h" +#include "sbndcode/CRT/CRTUtils/TPCGeoUtil.h" // ROOT #include "TPolyLine3D.h" @@ -93,6 +95,9 @@ namespace sbnd::crt { fhicl::Atom TPCTrackMatchLabel { Name("TPCTrackMatchLabel") }; + fhicl::Atom TPCTrackLabel { + Name("TPCTrackLabel") + }; fhicl::Atom SaveRoot { Name("SaveRoot"), @@ -140,6 +145,12 @@ namespace sbnd::crt { fhicl::Atom DrawTPCMatching { Name("DrawTPCMatching") }; + fhicl::Atom OnlyDrawMatched { + Name("OnlyDrawMatched") + }; + fhicl::Atom DisplayMatchScore { + Name("DisplayMatchScore") + }; fhicl::Atom ChoseTaggers { Name("ChoseTaggers") @@ -234,6 +245,9 @@ namespace sbnd::crt { void SetDrawStripHits(bool tf); void SetDrawClusters(bool tf); + void SetMinTime(double time); + void SetMaxTime(double time); + void SetPrint(bool tf); void SetHighlightedModules(std::vector hm); @@ -261,6 +275,7 @@ namespace sbnd::crt { art::InputTag fTrackLabel; art::InputTag fTPCSpacePointMatchLabel; art::InputTag fTPCTrackMatchLabel; + art::InputTag fTPCTrackLabel; bool fSaveRoot; bool fSaveViews; @@ -278,6 +293,8 @@ namespace sbnd::crt { bool fDrawSpacePoints; bool fDrawTracks; bool fDrawTPCMatching; + bool fOnlyDrawMatched; + bool fDisplayMatchScore; bool fChoseTaggers; std::vector fChosenTaggers; diff --git a/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl b/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl index ffb78aa37..90a2200a0 100644 --- a/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl +++ b/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl @@ -17,25 +17,28 @@ crteventdisplayalg_sbnd: ClusterLabel: "crtclustering" SpacePointLabel: "crtspacepoints" TrackLabel: "crttracks" - TPCTrackMatchLabel: "crttrackmatching" TPCSpacePointMatchLabel: "crtspacepointmatching" + TPCTrackMatchLabel: "crttrackmatching" + TPCTrackLabel: "pandoraTrack" SaveRoot: true SaveViews: false - DrawTaggers: true - DrawModules: false - DrawFEBs: false - DrawFEBEnds: false - DrawStrips: false - DrawTPC: true - DrawTrueTracks: true - DrawSimDeposits: false - DrawStripHits: true - DrawClusters: true - DrawSpacePoints: true - DrawTracks: true - DrawTPCMatching: false + DrawTaggers: true + DrawModules: false + DrawFEBs: false + DrawFEBEnds: false + DrawStrips: false + DrawTPC: true + DrawTrueTracks: true + DrawSimDeposits: false + DrawStripHits: true + DrawClusters: true + DrawSpacePoints: true + DrawTracks: true + DrawTPCMatching: false + OnlyDrawMatched: false + DisplayMatchScore: false ## If chose taggers is set to true then the vector ChosenTaggers ## is used to determine which taggers are drawn. The integer values @@ -86,11 +89,12 @@ crteventdisplayalg_sbnd_data.DrawClusters: false crteventdisplayalg_sbnd_data.DrawSpacePoints: true crteventdisplayalg_sbnd_data.DrawTracks: true crteventdisplayalg_sbnd_data.DrawTPCMatching: true +crteventdisplayalg_sbnd_data.OnlyDrawMatched: true crteventdisplayalg_sbnd_data.MC: false -crteventdisplayalg_sbnd_data.DrawModules: true crteventdisplayalg_sbnd_data.DrawTrueTracks: false -crteventdisplayalg_sbnd_data.ChoseTaggers: true -crteventdisplayalg_sbnd_data.ChosenTaggers: [ 0, 1, 2, 3, 4 ] +crteventdisplayalg_sbnd_data.ChoseTaggers: false crteventdisplayalg_sbnd_data.UseTs0: true +crteventdisplayalg_sbnd_data.MinTime: -1.5e6 +crteventdisplayalg_sbnd_data.MaxTime: 1.5e6 END_PROLOG From 6eab34146bc3ec00082199136da9dd660d67ac9b Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Thu, 27 Feb 2025 08:21:54 -0600 Subject: [PATCH 07/18] Add the possibility to save event displays whilst running matching --- sbndcode/CRT/CRTTPCMatching/CMakeLists.txt | 2 + .../CRTSpacePointMatching_module.cc | 28 +++++++- .../CRTTPCMatching/CRTTrackMatching_module.cc | 34 ++++++++-- .../crttpcmatchingproducers_sbnd.fcl | 68 +++++++++++-------- .../run_crttpcmatching_data.fcl | 11 +++ 5 files changed, 104 insertions(+), 39 deletions(-) diff --git a/sbndcode/CRT/CRTTPCMatching/CMakeLists.txt b/sbndcode/CRT/CRTTPCMatching/CMakeLists.txt index c1ea3afaa..c6a802419 100644 --- a/sbndcode/CRT/CRTTPCMatching/CMakeLists.txt +++ b/sbndcode/CRT/CRTTPCMatching/CMakeLists.txt @@ -11,6 +11,7 @@ simple_plugin( sbnobj::SBND_CRT sbndcode_GeoWrappers sbndcode_CRT_CRTTPCMatching + sbndcode_CRT_CRTEventDisplay ) simple_plugin( @@ -18,6 +19,7 @@ simple_plugin( sbnobj::SBND_CRT sbndcode_GeoWrappers sbndcode_CRT_CRTTPCMatching + sbndcode_CRT_CRTEventDisplay ) install_fhicl() diff --git a/sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatching_module.cc b/sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatching_module.cc index cde07b411..9866316aa 100644 --- a/sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatching_module.cc +++ b/sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatching_module.cc @@ -23,10 +23,12 @@ #include "lardata/Utilities/AssociationUtil.h" #include "lardata/DetectorInfoServices/DetectorPropertiesService.h" +#include "lardata/DetectorInfoServices/DetectorClocksService.h" #include "sbnobj/SBND/CRT/CRTSpacePoint.hh" #include "sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatchAlg.h" +#include "sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h" #include @@ -51,6 +53,10 @@ class sbnd::crt::CRTSpacePointMatching : public art::EDProducer { CRTSpacePointMatchAlg fMatchingAlg; art::InputTag fTPCTrackModuleLabel; art::InputTag fCRTSpacePointModuleLabel; + + CRTEventDisplayAlg fCRTEventDisplayAlg; + + bool fMakeEventDisplays; }; @@ -59,9 +65,12 @@ sbnd::crt::CRTSpacePointMatching::CRTSpacePointMatching(fhicl::ParameterSet cons , fMatchingAlg(p.get("MatchingAlg")) , fTPCTrackModuleLabel(p.get("TPCTrackModuleLabel")) , fCRTSpacePointModuleLabel(p.get("CRTSpacePointModuleLabel")) - { + , fCRTEventDisplayAlg(p.get("CRTEventDisplayAlg")) + , fMakeEventDisplays(p.get("MakeEventDisplays")) +{ + if(!fMakeEventDisplays) produces>(); - } +} void sbnd::crt::CRTSpacePointMatching::produce(art::Event& e) { @@ -104,6 +113,8 @@ void sbnd::crt::CRTSpacePointMatching::produce(art::Event& e) std::set used_crt_sps; + int nmatches = 0; + for(auto const &candidate : candidates) { if(used_crt_sps.count(candidate.thisSP.key()) == 0) @@ -111,11 +122,22 @@ void sbnd::crt::CRTSpacePointMatching::produce(art::Event& e) const anab::T0 t0(candidate.time * 1e3, 0, 0, 0, candidate.score); crtSPTPCTrackAssn->addSingle(candidate.thisSP, candidate.thisTrack, t0); + if(fMakeEventDisplays) + { + fCRTEventDisplayAlg.SetPrint(false); + fCRTEventDisplayAlg.SetMinTime(candidate.time * 1e3 - 1e3); + fCRTEventDisplayAlg.SetMaxTime(candidate.time * 1e3 + 1e3); + auto const clockData = art::ServiceHandle()->DataFor(e); + fCRTEventDisplayAlg.Draw(clockData, e, Form("crtEventDisplayEvent%iSPMatch%i", e.event(), nmatches)); + } + used_crt_sps.insert(candidate.thisSP.key()); + ++nmatches; } } - e.put(std::move(crtSPTPCTrackAssn)); + if(!fMakeEventDisplays) + e.put(std::move(crtSPTPCTrackAssn)); } DEFINE_ART_MODULE(sbnd::crt::CRTSpacePointMatching) diff --git a/sbndcode/CRT/CRTTPCMatching/CRTTrackMatching_module.cc b/sbndcode/CRT/CRTTPCMatching/CRTTrackMatching_module.cc index a5bc3e293..8cf234fc0 100644 --- a/sbndcode/CRT/CRTTPCMatching/CRTTrackMatching_module.cc +++ b/sbndcode/CRT/CRTTPCMatching/CRTTrackMatching_module.cc @@ -20,10 +20,12 @@ #include "lardata/Utilities/AssociationUtil.h" #include "lardata/DetectorInfoServices/DetectorPropertiesService.h" +#include "lardata/DetectorInfoServices/DetectorClocksService.h" #include "sbnobj/SBND/CRT/CRTTrack.hh" #include "sbndcode/CRT/CRTTPCMatching/CRTTrackMatchAlg.h" +#include "sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h" #include @@ -45,9 +47,13 @@ class sbnd::crt::CRTTrackMatching : public art::EDProducer { private: - CRTTrackMatchAlg fMatchingAlg; - art::InputTag fTPCTrackModuleLabel; - art::InputTag fCRTTrackModuleLabel; + CRTTrackMatchAlg fMatchingAlg; + art::InputTag fTPCTrackModuleLabel; + art::InputTag fCRTTrackModuleLabel; + + CRTEventDisplayAlg fCRTEventDisplayAlg; + + bool fMakeEventDisplays; }; @@ -56,9 +62,12 @@ sbnd::crt::CRTTrackMatching::CRTTrackMatching(fhicl::ParameterSet const& p) , fMatchingAlg(p.get("MatchingAlg")) , fTPCTrackModuleLabel(p.get("TPCTrackModuleLabel")) , fCRTTrackModuleLabel(p.get("CRTTrackModuleLabel")) - { + , fCRTEventDisplayAlg(p.get("CRTEventDisplayAlg")) + , fMakeEventDisplays(p.get("MakeEventDisplays")) +{ + if(!fMakeEventDisplays) produces>(); - } +} void sbnd::crt::CRTTrackMatching::produce(art::Event& e) { @@ -101,6 +110,8 @@ void sbnd::crt::CRTTrackMatching::produce(art::Event& e) std::set used_crt_tracks; + int nmatches = 0; + for(auto const &candidate : candidates) { if(used_crt_tracks.count(candidate.thisCRTTrack.key()) == 0) @@ -108,11 +119,22 @@ void sbnd::crt::CRTTrackMatching::produce(art::Event& e) const anab::T0 t0(candidate.time * 1e3, 0, 0, 0, candidate.score); crtTrackTPCTrackAssn->addSingle(candidate.thisCRTTrack, candidate.thisTPCTrack, t0); + if(fMakeEventDisplays) + { + fCRTEventDisplayAlg.SetPrint(false); + fCRTEventDisplayAlg.SetMinTime(candidate.time * 1e3 - 1e3); + fCRTEventDisplayAlg.SetMaxTime(candidate.time * 1e3 + 1e3); + auto const clockData = art::ServiceHandle()->DataFor(e); + fCRTEventDisplayAlg.Draw(clockData, e, Form("crtEventDisplayEvent%iTrackMatch%i", e.event(), nmatches)); + } + used_crt_tracks.insert(candidate.thisCRTTrack.key()); + ++nmatches; } } - e.put(std::move(crtTrackTPCTrackAssn)); + if(!fMakeEventDisplays) + e.put(std::move(crtTrackTPCTrackAssn)); } DEFINE_ART_MODULE(sbnd::crt::CRTTrackMatching) diff --git a/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl b/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl index d1a044e0e..4a5221635 100644 --- a/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl +++ b/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl @@ -1,18 +1,20 @@ +#include "crteventdisplayalg_sbnd.fcl" + BEGIN_PROLOG crtspacepointmatchalg_sbnd: { - TrackDirectionFrac: 0.5 - DCALimit: 70. - MinTPCTrackLength: 10. - DirMethod: 1 - DCAuseBox: false - DCAoverLength: false - PECut: 60. - MaxUncert: 20. - TPCTrackLabel: "pandoraTrack" - CRTSpacePointLabel: "crtspacepoints" - UseTs0: false + TrackDirectionFrac: 0.5 + DCALimit: 70. + MinTPCTrackLength: 10. + DirMethod: 1 + DCAuseBox: false + DCAoverLength: false + PECut: 60. + MaxUncert: 20. + TPCTrackLabel: "pandoraTrack" + CRTSpacePointLabel: "crtspacepoints" + UseTs0: false } crtspacepointmatchalg_data_sbnd: @local::crtspacepointmatchalg_sbnd @@ -20,25 +22,28 @@ crtspacepointmatchalg_data_sbnd.UseTs0: true crtspacepointmatchproducer_sbnd: { - CRTSpacePointModuleLabel: "crtspacepoints" - TPCTrackModuleLabel: @local::crtspacepointmatchalg_sbnd.TPCTrackLabel - PFPModuleLabel: "pandora" - MatchingAlg: @local::crtspacepointmatchalg_sbnd - module_type: "CRTSpacePointMatching" + CRTSpacePointModuleLabel: "crtspacepoints" + TPCTrackModuleLabel: @local::crtspacepointmatchalg_sbnd.TPCTrackLabel + PFPModuleLabel: "pandora" + MatchingAlg: @local::crtspacepointmatchalg_sbnd + CRTEventDisplayAlg: @local::crteventdisplayalg_sbnd + MakeEventDisplays: false + module_type: "CRTSpacePointMatching" } crtspacepointmatchproducer_data_sbnd: @local::crtspacepointmatchproducer_sbnd -crtspacepointmatchproducer_data_sbnd.MatchingAlg: @local::crtspacepointmatchalg_data_sbnd +crtspacepointmatchproducer_data_sbnd.MatchingAlg: @local::crtspacepointmatchalg_data_sbnd +crtspacepointmatchproducer_data_sbnd.CRTEventDisplayAlg: @local::crteventdisplayalg_sbnd_data crttrackmatchalg_sbnd: { - MaxAngleDiff: 0.4 - MaxDCA: 70. - MaxScore: 150. - MinTPCTrackLength: 5. - SelectionMetric: "score" - TPCTrackLabel: "pandoraTrack" - UseTs0: false + MaxAngleDiff: 0.4 + MaxDCA: 70. + MaxScore: 150. + MinTPCTrackLength: 5. + SelectionMetric: "score" + TPCTrackLabel: "pandoraTrack" + UseTs0: false } crttrackmatchalg_data_sbnd: @local::crttrackmatchalg_sbnd @@ -46,14 +51,17 @@ crttrackmatchalg_data_sbnd.UseTs0: true crttrackmatchproducer_sbnd: { - CRTTrackModuleLabel: "crttracks" - TPCTrackModuleLabel: @local::crttrackmatchalg_sbnd.TPCTrackLabel - PFPModuleLabel: "pandora" - MatchingAlg: @local::crttrackmatchalg_sbnd - module_type: "CRTTrackMatching" + CRTTrackModuleLabel: "crttracks" + TPCTrackModuleLabel: @local::crttrackmatchalg_sbnd.TPCTrackLabel + PFPModuleLabel: "pandora" + MatchingAlg: @local::crttrackmatchalg_sbnd + CRTEventDisplayAlg: @local::crteventdisplayalg_sbnd + MakeEventDisplays: false + module_type: "CRTTrackMatching" } crttrackmatchproducer_data_sbnd: @local::crttrackmatchproducer_sbnd -crttrackmatchproducer_data_sbnd.MatchingAlg: @local::crttrackmatchalg_sbnd +crttrackmatchproducer_data_sbnd.MatchingAlg: @local::crttrackmatchalg_data_sbnd +crttrackmatchproducer_data_sbnd.CRTEventDisplayAlg: @local::crteventdisplayalg_sbnd_data END_PROLOG diff --git a/sbndcode/CRT/CRTTPCMatching/run_crttpcmatching_data.fcl b/sbndcode/CRT/CRTTPCMatching/run_crttpcmatching_data.fcl index d1723f5f6..34ddd634b 100644 --- a/sbndcode/CRT/CRTTPCMatching/run_crttpcmatching_data.fcl +++ b/sbndcode/CRT/CRTTPCMatching/run_crttpcmatching_data.fcl @@ -4,3 +4,14 @@ process_name: CRTTPCMatching physics.producers.crtspacepointmatching: @local::crtspacepointmatchproducer_data_sbnd physics.producers.crttrackmatching: @local::crttrackmatchproducer_data_sbnd + +physics.producers.crtspacepointmatchingevd: @local::crtspacepointmatchproducer_data_sbnd +physics.producers.crtspacepointmatchingevd.MakeEventDisplays: true +physics.producers.crtspacepointmatchingevd.CRTEventDisplayAlg.DrawTracks: false +physics.producers.crtspacepointmatchingevd.CRTEventDisplayAlg.DisplayMatchScore: true + +physics.producers.crttrackmatchingevd: @local::crttrackmatchproducer_data_sbnd +physics.producers.crttrackmatchingevd.MakeEventDisplays: true +physics.producers.crttrackmatchingevd.CRTEventDisplayAlg.DisplayMatchScore: true + +physics.reco: [ crtspacepointmatching, crttrackmatching, crtspacepointmatchingevd, crttrackmatchingevd ] From 1f07826841681b53b810761cf3d5e7c9e24828fa Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Thu, 27 Feb 2025 09:08:18 -0600 Subject: [PATCH 08/18] Fix unphysical track matching --- .../CRT/CRTTPCMatching/CRTTrackMatchAlg.cc | 30 ++++++++++++++++++- .../CRT/CRTTPCMatching/CRTTrackMatchAlg.h | 3 ++ 2 files changed, 32 insertions(+), 1 deletion(-) diff --git a/sbndcode/CRT/CRTTPCMatching/CRTTrackMatchAlg.cc b/sbndcode/CRT/CRTTPCMatching/CRTTrackMatchAlg.cc index 2e61154b6..61bf6c524 100644 --- a/sbndcode/CRT/CRTTPCMatching/CRTTrackMatchAlg.cc +++ b/sbndcode/CRT/CRTTPCMatching/CRTTrackMatchAlg.cc @@ -106,7 +106,9 @@ namespace sbnd::crt { { std::vector> candidates; - const int driftDirection = TPCGeoUtil::DriftDirectionFromHits(fGeometryService, hits); + const int driftDirection = TPCGeoUtil::DriftDirectionFromHits(fGeometryService, hits); + const std::pair xLimits = TPCGeoUtil::XLimitsFromHits(fGeometryService, hits); + const std::pair t0MinMax = TrackT0Range(detProp, tpcTrack->Vertex().X(), tpcTrack->End().X(), driftDirection, xLimits); const geo::TPCID tpcID = hits[0]->WireID().asTPCID(); const geo::TPCGeo& tpcGeo = fGeometryService->GetElement(tpcID); @@ -119,6 +121,10 @@ namespace sbnd::crt { continue; const double crtTime = fUseTs0 ? crtTrack->Ts0() * 1e-3 : crtTrack->Ts1() * 1e-3; + + if(!(crtTime > t0MinMax.first - 10. && crtTime < t0MinMax.second + 10.)) + continue; + const double shift = driftDirection * crtTime * detProp.DriftVelocity(); geo::Point_t start = tpcTrack->Vertex(); @@ -321,4 +327,26 @@ namespace sbnd::crt { return aveDCA / usedPts; } + + std::pair CRTTrackMatchAlg::TrackT0Range(detinfo::DetectorPropertiesData const &detProp, const double startX, + const double endX, const int driftDirection, const std::pair xLimits) + { + if(driftDirection == 0) + return std::make_pair(0, 0); + + const double driftVel = driftDirection * detProp.DriftVelocity(); + + const double maxX = std::max(startX, endX); + const double maxLimit = std::max(xLimits.first, xLimits.second); + const double maxShift = maxLimit - maxX; + + const double minX = std::min(startX, endX); + const double minLimit = std::min(xLimits.first, xLimits.second); + const double minShift = minLimit - minX; + + const double t0max = maxShift / driftVel; + const double t0min = minShift / driftVel; + + return std::make_pair(std::min(t0min, t0max), std::max(t0min, t0max)); + } } diff --git a/sbndcode/CRT/CRTTPCMatching/CRTTrackMatchAlg.h b/sbndcode/CRT/CRTTPCMatching/CRTTrackMatchAlg.h index 35ae6eb04..bfd45c056 100644 --- a/sbndcode/CRT/CRTTPCMatching/CRTTrackMatchAlg.h +++ b/sbndcode/CRT/CRTTPCMatching/CRTTrackMatchAlg.h @@ -147,6 +147,9 @@ namespace sbnd::crt { double AveDCABetweenTracks(const art::Ptr &tpcTrack, const art::Ptr &crtTrack, const double shift); + std::pair TrackT0Range(detinfo::DetectorPropertiesData const &detProp, const double startX, + const double endX, const int driftDirection, const std::pair xLimits); + private: geo::GeometryCore const* fGeometryService; From e0f5f5d99b67efd34ec58b1857225377d3fafc55 Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Fri, 28 Feb 2025 13:01:51 -0600 Subject: [PATCH 09/18] Move general functions out of specific alg --- .../CRTTPCMatching/CRTSpacePointMatchAlg.cc | 58 +------------------ .../CRTTPCMatching/CRTSpacePointMatchAlg.h | 4 -- .../crttpcmatchingproducers_sbnd.fcl | 4 +- sbndcode/CRT/CRTUtils/CRTCommonUtils.cc | 46 +++++++++++++++ sbndcode/CRT/CRTUtils/CRTCommonUtils.h | 8 +++ 5 files changed, 58 insertions(+), 62 deletions(-) diff --git a/sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatchAlg.cc b/sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatchAlg.cc index 83c960458..ac4a04aa7 100644 --- a/sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatchAlg.cc +++ b/sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatchAlg.cc @@ -90,9 +90,9 @@ namespace sbnd::crt { std::pair startEndDir; if(fDirMethod==2) - startEndDir = AverageTrackDirections(track, fTrackDirectionFrac); + startEndDir = CRTCommonUtils::AverageTrackDirections(track, fTrackDirectionFrac); else - startEndDir = TrackDirections(detProp, track, fTrackDirectionFrac, crtTime, driftDirection); + startEndDir = CRTCommonUtils::TrackDirections(track); const geo::Vector_t startDir = startEndDir.first; const geo::Vector_t endDir = startEndDir.second; @@ -181,58 +181,4 @@ namespace sbnd::crt { else return CRTCommonUtils::SimpleDCA(crtSP, trackStart, trackDir); } - - std::pair CRTSpacePointMatchAlg::AverageTrackDirections(const art::Ptr &track, const double frac) - { - const unsigned N = track->NumberTrajectoryPoints(); - const unsigned NValid = track->CountValidPoints(); - const unsigned NValidFrac = std::floor(NValid * frac); - - geo::Vector_t forwardDir(0, 0, 0), backwardDir(0, 0, 0); - - unsigned iValidForward = 0, iValidBackward = 0; - - for(unsigned i = 0; i < N; ++i) - { - if(track->HasValidPoint(i) && iValidForward < NValidFrac) - { - forwardDir += track->DirectionAtPoint(i); - ++iValidForward; - } - - if(track->HasValidPoint(N - (i + 1)) && iValidBackward < NValidFrac) - { - backwardDir += track->DirectionAtPoint(N - (i + 1)); - ++iValidBackward; - } - } - - forwardDir /= NValidFrac; - backwardDir /= NValidFrac; - - return std::make_pair(forwardDir, backwardDir); - - } - - std::pair CRTSpacePointMatchAlg::TrackDirections(detinfo::DetectorPropertiesData const &detProp, const art::Ptr &track, - const double frac, const double CRTtime, const int driftDirection) - { - const unsigned N = track->NumberTrajectoryPoints(); - const unsigned NMid = std::floor(N * frac); - - geo::Point_t start = track->Start(); - geo::Point_t end = track->End(); - geo::Point_t mid = track->LocationAtPoint(NMid); - - const double xshift = driftDirection * CRTtime * detProp.DriftVelocity(); - - start.SetX(start.X() + xshift); - end.SetX(end.X() + xshift); - mid.SetX(mid.X() + xshift); - - const geo::Vector_t startDir = (mid - start).Unit(); - const geo::Vector_t endDir = (mid - end).Unit(); - - return std::make_pair(startDir, endDir); - } } diff --git a/sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatchAlg.h b/sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatchAlg.h index e9cb845f5..19bc01181 100644 --- a/sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatchAlg.h +++ b/sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatchAlg.h @@ -167,10 +167,6 @@ namespace sbnd::crt { const geo::Vector_t &trackDir, const art::Ptr &crtSP, const int driftDirection, const double t0, const art::Event &e); - std::pair AverageTrackDirections(const art::Ptr &track, const double frac); - std::pair TrackDirections(detinfo::DetectorPropertiesData const &detProp, const art::Ptr &track, - const double frac, const double CRTtime, const int driftDirection); - private: geo::GeometryCore const* fGeometryService; diff --git a/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl b/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl index 4a5221635..a02373416 100644 --- a/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl +++ b/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl @@ -4,10 +4,10 @@ BEGIN_PROLOG crtspacepointmatchalg_sbnd: { - TrackDirectionFrac: 0.5 + TrackDirectionFrac: 0.25 DCALimit: 70. MinTPCTrackLength: 10. - DirMethod: 1 + DirMethod: 2 DCAuseBox: false DCAoverLength: false PECut: 60. diff --git a/sbndcode/CRT/CRTUtils/CRTCommonUtils.cc b/sbndcode/CRT/CRTUtils/CRTCommonUtils.cc index 87409a24b..8d781463f 100644 --- a/sbndcode/CRT/CRTUtils/CRTCommonUtils.cc +++ b/sbndcode/CRT/CRTUtils/CRTCommonUtils.cc @@ -260,4 +260,50 @@ namespace sbnd::crt { return dP.R(); } + + std::pair CRTCommonUtils::AverageTrackDirections(const art::Ptr &track, const double frac) + { + const unsigned N = track->NumberTrajectoryPoints(); + const unsigned NValid = track->CountValidPoints(); + const unsigned NValidFrac = std::floor(NValid * frac); + + geo::Vector_t forwardDir(0, 0, 0), backwardDir(0, 0, 0); + + unsigned iValidForward = 0, iValidBackward = 0; + + for(unsigned i = 0; i < N; ++i) + { + if(track->HasValidPoint(i) && iValidForward < NValidFrac) + { + forwardDir += track->DirectionAtPoint(i); + ++iValidForward; + } + + if(track->HasValidPoint(N - (i + 1)) && iValidBackward < NValidFrac) + { + backwardDir += track->DirectionAtPoint(N - (i + 1)); + ++iValidBackward; + } + } + + forwardDir /= NValidFrac; + backwardDir /= NValidFrac; + + return std::make_pair(forwardDir, backwardDir); + } + + std::pair CRTCommonUtils::TrackDirections(const art::Ptr &track) + { + const unsigned N = track->NumberTrajectoryPoints(); + const unsigned NMid = std::floor(N * 0.5); + + geo::Point_t start = track->Start(); + geo::Point_t end = track->End(); + geo::Point_t mid = track->LocationAtPoint(NMid); + + const geo::Vector_t startDir = (mid - start).Unit(); + const geo::Vector_t endDir = (mid - end).Unit(); + + return std::make_pair(startDir, endDir); + } } diff --git a/sbndcode/CRT/CRTUtils/CRTCommonUtils.h b/sbndcode/CRT/CRTUtils/CRTCommonUtils.h index 251d3e02d..180113308 100644 --- a/sbndcode/CRT/CRTUtils/CRTCommonUtils.h +++ b/sbndcode/CRT/CRTUtils/CRTCommonUtils.h @@ -11,6 +11,8 @@ #include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h" +#include "lardataobj/RecoBase/Track.h" + #include "sbnobj/SBND/CRT/CRTEnums.hh" #include "sbnobj/SBND/CRT/CRTSpacePoint.hh" @@ -56,6 +58,12 @@ namespace sbnd::crt { // Returns the distance between an infinite line and a segment double LineSegmentDistance(const geo::Point_t &start1, const geo::Point_t &end1, const geo::Point_t &start2, const geo::Point_t &end2); + + // Returns the average direction for the first and last frac(tion) of trajectory points in a track + std::pair AverageTrackDirections(const art::Ptr &track, const double frac); + + // Returns the directions based on the line from the middle to either end of the track + std::pair TrackDirections(const art::Ptr &track); } } From 7cbfa5a287f36d0034ec34b14cd1b0541e848a2b Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Fri, 28 Feb 2025 13:03:26 -0600 Subject: [PATCH 10/18] Change how track directions are visualised --- .../CRT/CRTEventDisplay/CRTEventDisplayAlg.cc | 24 +++++++++++-------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc index 0e08c638c..6feaa0136 100644 --- a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc +++ b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc @@ -481,15 +481,17 @@ namespace sbnd::crt { const anab::T0 t0Match = CRTSPstoTPCTracks.data(spacepoint.key()).ref(); - matchScore = t0Match.TriggerConfidence(); - geo::Point_t startTPC = TPCTrack->Start(); - geo::Point_t endTPC = TPCTrack->End(); - const geo::Vector_t startDirTPC = TPCTrack->StartDirection(); - const geo::Vector_t endDirTPC = TPCTrack->EndDirection(); + matchScore = t0Match.TriggerConfidence(); + geo::Point_t startTPC = TPCTrack->Start(); + geo::Point_t endTPC = TPCTrack->End(); startTPC.SetX(startTPC.X() + shift); endTPC.SetX(endTPC.X() + shift); + const std::pair directions = CRTCommonUtils::AverageTrackDirections(TPCTrack, 0.2); + const geo::Vector_t startDirTPC = directions.first; + const geo::Vector_t endDirTPC = directions.second; + TPolyLine3D *lineTPC = new TPolyLine3D(2); lineTPC->SetPoint(0, startTPC.X(), startTPC.Y(), startTPC.Z()); lineTPC->SetPoint(1, endTPC.X(), endTPC.Y(), endTPC.Z()); @@ -624,15 +626,17 @@ namespace sbnd::crt { const anab::T0 t0Match = CRTTrackstoTPCTracks.data(track.key()).ref(); - matchScore = t0Match.TriggerConfidence(); - geo::Point_t startTPC = TPCTrack->Start(); - geo::Point_t endTPC = TPCTrack->End(); - const geo::Vector_t startDirTPC = TPCTrack->StartDirection(); - const geo::Vector_t endDirTPC = TPCTrack->EndDirection(); + matchScore = t0Match.TriggerConfidence(); + geo::Point_t startTPC = TPCTrack->Start(); + geo::Point_t endTPC = TPCTrack->End(); startTPC.SetX(startTPC.X() + shift); endTPC.SetX(endTPC.X() + shift); + const std::pair directions = CRTCommonUtils::AverageTrackDirections(TPCTrack, 0.2); + const geo::Vector_t startDirTPC = directions.first; + const geo::Vector_t endDirTPC = directions.second; + TPolyLine3D *lineTPC = new TPolyLine3D(2); lineTPC->SetPoint(0, startTPC.X(), startTPC.Y(), startTPC.Z()); lineTPC->SetPoint(1, endTPC.X(), endTPC.Y(), endTPC.Z()); From dd002aad17b010139219ceda55d9a2b8b53648e3 Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Fri, 28 Feb 2025 14:17:15 -0600 Subject: [PATCH 11/18] Make quicker option for MC with no backtracking --- sbndcode/CRT/CRTAna/CRTAnalysis_module.cc | 97 +++++++++++++++++------ 1 file changed, 74 insertions(+), 23 deletions(-) diff --git a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc index d4b1e1816..6b764988b 100644 --- a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc +++ b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc @@ -97,7 +97,7 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer { fCRTClusterModuleLabel, fCRTSpacePointModuleLabel, fCRTTrackModuleLabel, fTPCTrackModuleLabel, fCRTSpacePointMatchingModuleLabel, fCRTTrackMatchingModuleLabel, fPFPModuleLabel, fPTBModuleLabel, fTDCModuleLabel, fTimingReferenceModuleLabel; - bool fDebug, fDataMode, fNoTPC, fHasPTB, fHasTDC; + bool fDebug, fDataMode, fNoTPC, fHasPTB, fHasTDC, fTruthMatch; //! Adding some of the reco parameters to save corrections double fPEAttenuation, fTimeWalkNorm, fTimeWalkScale, fPropDelay; @@ -348,13 +348,14 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fNoTPC = p.get("NoTPC", false); fHasPTB = p.get("HasPTB", false); fHasTDC = p.get("HasTDC", false); + fTruthMatch = p.get("TruthMatch", true); //! Adding some of the reco parameters to save corrections fPEAttenuation = p.get("PEAttenuation", 1.0); fTimeWalkNorm = p.get("TimeWalkNorm", 0.0); fTimeWalkScale = p.get("TimeWalkScale", 0.0); fPropDelay = p.get("PropDelay", 0.0); - if(!fDataMode) + if(!fDataMode && fTruthMatch) fCRTBackTrackerAlg = CRTBackTrackerAlg(p.get("CRTBackTrackerAlg", fhicl::ParameterSet())); art::ServiceHandle fs; @@ -632,22 +633,29 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) { - if(!fDataMode) + _run = e.id().run(); + _subrun = e.id().subRun(); + _event = e.id().event(); + + if(fDebug) + std::cout << "CRTAnalysis -- This is event " << _run << "-" << _subrun << "-" << _event << std::endl; + + if(!fDataMode && fTruthMatch) { + if(fDebug) + std::cout << "CRTAnalysis -- Setting up back track maps" << std::endl; + fCRTBackTrackerAlg.SetupMaps(e); fCRTBackTrackerAlg.RunSpacePointRecoStatusChecks(e); fCRTBackTrackerAlg.RunTrackRecoStatusChecks(e); } - _run = e.id().run(); - _subrun = e.id().subRun(); - _event = e.id().event(); - - if(fDebug) std::cout << "This is event " << _run << "-" << _subrun << "-" << _event << std::endl; - _crt_timing_reference_type = -1; _crt_timing_reference_channel = -1; + if(fDebug) + std::cout << "CRTAnalysis -- Processing Timing Reference Info" << std::endl; + art::Handle TimingReferenceHandle; e.getByLabel(fTimingReferenceModuleLabel, TimingReferenceHandle); if(TimingReferenceHandle.isValid()) @@ -658,6 +666,9 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) if(fHasPTB) { + if(fDebug) + std::cout << "CRTAnalysis -- Processing PTB" << std::endl; + // Get PTBs art::Handle> PTBHandle; e.getByLabel(fPTBModuleLabel, PTBHandle); @@ -674,6 +685,9 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) if(fHasTDC) { + if(fDebug) + std::cout << "CRTAnalysis -- Processing TDC" << std::endl; + // Get TDCs art::Handle> TDCHandle; e.getByLabel(fTDCModuleLabel, TDCHandle); @@ -688,6 +702,9 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) AnalyseTDCs(TDCVec); } + if(fDebug) + std::cout << "CRTAnalysis -- Processing FEBDatas" << std::endl; + // Get FEBDatas art::Handle> FEBDataHandle; e.getByLabel(fFEBDataModuleLabel, FEBDataHandle); @@ -701,6 +718,9 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) // Fill FEBData variables AnalyseFEBDatas(FEBDataVec); + if(fDebug) + std::cout << "CRTAnalysis -- Processing StripHits" << std::endl; + // Get CRTStripHits art::Handle> CRTStripHitHandle; e.getByLabel(fCRTStripHitModuleLabel, CRTStripHitHandle); @@ -714,6 +734,9 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) // Fill CRTStripHit variables AnalyseCRTStripHits(e, CRTStripHitVec); + if(fDebug) + std::cout << "CRTAnalysis -- Processing Clusters" << std::endl; + // Get CRTClusters art::Handle> CRTClusterHandle; e.getByLabel(fCRTClusterModuleLabel, CRTClusterHandle); @@ -733,6 +756,9 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) // Fill CRTCluster variables AnalyseCRTClusters(e, CRTClusterVec, clustersToSpacePoints, clustersToStripHits); + if(fDebug) + std::cout << "CRTAnalysis -- Processing Tracks" << std::endl; + // Get CRTTracks art::Handle> CRTTrackHandle; e.getByLabel(fCRTTrackModuleLabel, CRTTrackHandle); @@ -786,6 +812,9 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) if(fDataMode) { + if(fDebug) + std::cout << "CRTAnalysis -- Processing TPC" << std::endl; + // Fill TPC matching variables AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, PFPHandle); @@ -795,6 +824,9 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) return; } + if(fDebug) + std::cout << "CRTAnalysis -- Processing MCParticles" << std::endl; + // Get MCParticles art::Handle> MCParticleHandle; e.getByLabel(fMCParticleModuleLabel, MCParticleHandle); @@ -808,6 +840,9 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) // Fill MCParticle variables AnalyseMCParticles(MCParticleVec); + if(fDebug) + std::cout << "CRTAnalysis -- Processing Sim Deposits" << std::endl; + // Get SimDeposits art::Handle> SimDepositHandle; e.getByLabel(fSimDepositModuleLabel, SimDepositHandle); @@ -821,21 +856,37 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) // Fill SimDeposit variables AnalyseSimDeposits(SimDepositVec); - // Get Map of TrueDeposits per tagger from BackTracker - std::map spacePointRecoStatusMap = fCRTBackTrackerAlg.GetSpacePointRecoStatusMap(); + if(fDebug) + std::cout << "CRTAnalysis -- Processing CRT Back Tracking" << std::endl; + + if(fTruthMatch) + { + // Get Map of TrueDeposits per tagger from BackTracker + std::map spacePointRecoStatusMap = fCRTBackTrackerAlg.GetSpacePointRecoStatusMap(); + + // Fill TrueDeposit variables + AnalyseTrueDepositsPerTagger(spacePointRecoStatusMap); - // Fill TrueDeposit variables - AnalyseTrueDepositsPerTagger(spacePointRecoStatusMap); + // Get Map of TrueDeposits from BackTracker + std::map> trackRecoStatusMap = fCRTBackTrackerAlg.GetTrackRecoStatusMap(); - // Get Map of TrueDeposits from BackTracker - std::map> trackRecoStatusMap = fCRTBackTrackerAlg.GetTrackRecoStatusMap(); + // Fill TrueDeposit variables + AnalyseTrueDeposits(trackRecoStatusMap); - // Fill TrueDeposit variables - AnalyseTrueDeposits(trackRecoStatusMap); + if(fDebug) + std::cout << "CRTAnalysis -- Processing TPC" << std::endl; - // Fill TPC matching variables - AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, PFPHandle, spacePointRecoStatusMap, trackRecoStatusMap); + // Fill TPC matching variables + AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, PFPHandle, spacePointRecoStatusMap, trackRecoStatusMap); + } + else + { + if(fDebug) + std::cout << "CRTAnalysis -- Processing TPC" << std::endl; + // Fill TPC matching variables + AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, PFPHandle); + } // Fill the Tree fTree->Fill(); } @@ -1082,7 +1133,7 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTStripHits(const art::Event &e, const std: _sh_saturated1[i] = hit->Saturated1(); _sh_saturated2[i] = hit->Saturated2(); - if(!fDataMode) + if(!fDataMode && fTruthMatch) { const CRTBackTrackerAlg::TruthMatchMetrics truthMatch = fCRTBackTrackerAlg.TruthMatching(e, hit); const std::vector localpos = fCRTGeoAlg.StripWorldToLocalPos(hit->Channel(), truthMatch.deposit.x, truthMatch.deposit.y, truthMatch.deposit.z); @@ -1212,7 +1263,7 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std:: } - if(!fDataMode) + if(!fDataMode && fTruthMatch) { const CRTBackTrackerAlg::TruthMatchMetrics truthMatch = fCRTBackTrackerAlg.TruthMatching(e, cluster); _cl_truth_trackid[i] = truthMatch.trackid; @@ -1445,7 +1496,7 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTTracks(const art::Event &e, const std::ve } } - if(!fDataMode) + if(!fDataMode && fTruthMatch) { const CRTBackTrackerAlg::TruthMatchMetrics truthMatch = fCRTBackTrackerAlg.TruthMatching(e, track); _tr_truth_trackid[i] = truthMatch.trackid; @@ -1598,7 +1649,7 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_tr_score[nActualTracks] = -std::numeric_limits::max(); } - if(!fDataMode) + if(!fDataMode && fTruthMatch) { const std::vector> trackHits = tracksToHits.at(track.key()); const int trackid = fCRTBackTrackerAlg.RollUpID(TruthMatchUtils::TrueParticleIDFromTotalRecoHits(clockData,trackHits,true)); From f7c00bc3cb6d5cd20ca1e98d00ae8183d6301cfe Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Fri, 28 Feb 2025 14:27:21 -0600 Subject: [PATCH 12/18] Finish the job --- sbndcode/CRT/CRTAna/CRTAnalysis_module.cc | 10 +++++----- sbndcode/CRT/CRTAna/run_crtana.fcl | 12 +++++++++++- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc index 6b764988b..cd99d1f0b 100644 --- a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc +++ b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc @@ -423,7 +423,7 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fTree->Branch("sh_adc2", "std::vector", &_sh_adc2); fTree->Branch("sh_saturated1", "std::vector", &_sh_saturated1); fTree->Branch("sh_saturated2", "std::vector", &_sh_saturated2); - if(!fDataMode) + if(!fDataMode && fTruthMatch) { fTree->Branch("sh_truth_trackid", "std::vector", &_sh_truth_trackid); fTree->Branch("sh_truth_completeness", "std::vector", &_sh_truth_completeness); @@ -446,7 +446,7 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fTree->Branch("cl_sh_feb_mac5_set", "std::vector>", &_cl_sh_feb_mac5_set); fTree->Branch("cl_sh_time_walk_set", "std::vector>", &_cl_sh_time_walk_set); fTree->Branch("cl_sh_prop_delay_set", "std::vector>", &_cl_sh_prop_delay_set); - if(!fDataMode) + if(!fDataMode && fTruthMatch) { fTree->Branch("cl_truth_trackid", "std::vector", &_cl_truth_trackid); fTree->Branch("cl_truth_completeness", "std::vector", &_cl_truth_completeness); @@ -479,7 +479,7 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fTree->Branch("cl_sp_ets1", "std::vector", &_cl_sp_ets1); fTree->Branch("cl_sp_complete", "std::vector", &_cl_sp_complete); - if(!fDataMode) + if(!fDataMode && fTruthMatch) { fTree->Branch("td_tag_trackid", "std::vector", &_td_tag_trackid); fTree->Branch("td_tag_pdg", "std::vector", &_td_tag_pdg); @@ -524,7 +524,7 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fTree->Branch("tr_tagger3", "std::vector", &_tr_tagger3); fTree->Branch("tr_channel_set", "std::vector>", &_tr_channel_set); fTree->Branch("tr_adc_set", "std::vector>", &_tr_adc_set); - if(!fDataMode) + if(!fDataMode && fTruthMatch) { fTree->Branch("tr_truth_trackid", "std::vector", &_tr_truth_trackid); fTree->Branch("tr_truth_completeness", "std::vector", &_tr_truth_completeness); @@ -569,7 +569,7 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fTree->Branch("tpc_tr_ts0", "std::vector", &_tpc_tr_ts0); fTree->Branch("tpc_tr_ts1", "std::vector", &_tpc_tr_ts1); fTree->Branch("tpc_tr_score", "std::vector", &_tpc_tr_score); - if(!fDataMode) + if(!fDataMode && fTruthMatch) { fTree->Branch("tpc_truth_trackid", "std::vector", &_tpc_truth_trackid); fTree->Branch("tpc_truth_pdg", "std::vector", &_tpc_truth_pdg); diff --git a/sbndcode/CRT/CRTAna/run_crtana.fcl b/sbndcode/CRT/CRTAna/run_crtana.fcl index 1f38761a9..0dff6f7a2 100644 --- a/sbndcode/CRT/CRTAna/run_crtana.fcl +++ b/sbndcode/CRT/CRTAna/run_crtana.fcl @@ -21,7 +21,17 @@ physics: { analyzers: { - crtana: @local::crtana_sbnd + crtana: + { + module_type: "CRTAnalysis" + CRTGeoAlg: @local::crtgeoalg_sbnd + CRTBackTrackerAlg: @local::crtbacktrackeralg_sbnd + PEAttenuation: @local::sbnd_crtsim.DetSimParams.NpeScaleShift + PropDelay: @local::sbnd_crtsim.DetSimParams.PropDelay + TimeWalkNorm: @local::sbnd_crtsim.DetSimParams.TDelayNorm + TimeWalkScale: @local::sbnd_crtsim.DetSimParams.TDelayScale + TruthMatch: false + } } ana: [ crtana ] From 6ea3aa389a77de49d4e7a1a493959192d2fb92b7 Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Mon, 3 Mar 2025 08:38:38 -0600 Subject: [PATCH 13/18] Add CRT-TPC matching to data reco2 workflow --- .../JobConfigurations/standard/reco/reco2_data.fcl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl b/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl index 738d8f5a8..b160f68ea 100644 --- a/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl +++ b/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl @@ -14,14 +14,15 @@ services: physics.producers: { @table::physics.producers - crtclustering: @local::crtclusterproducer_data_sbnd - crtspacepoints: @local::crtspacepointproducer_data_sbnd - crttracks: @local::crttrackproducer_data_sbnd - opt0finder: @local::sbnd_opt0_finder_data + crtclustering: @local::crtclusterproducer_data_sbnd + crtspacepoints: @local::crtspacepointproducer_data_sbnd + crttracks: @local::crttrackproducer_data_sbnd + crtspacepointmatching: @local::crtspacepointmatchproducer_data_sbnd + crttrackmatching: @local::crttrackmatchproducer_data_sbnd } physics.reco2: [ pandora, pandoraTrack, pandoraShower, pandoraShowerSBN, pandoraCaloData, pandoraPidData, caloskimCalorimetry, - crtclustering, crtspacepoints, crttracks, cvn, opt0finder] + crtclustering, crtspacepoints, crttracks, crtspacepointmatching, crttrackmatching, cvn ] #The next 3 lines need to be commented out once data use pandoraSCE by default. physics.producers.cvn.SliceLabel: "pandora" From 9132413a4b3b36ea134c93ea0a8198c94b7a9348 Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Tue, 4 Mar 2025 04:46:37 -0600 Subject: [PATCH 14/18] Merge remote-tracking branch 'origin/feature/hlay_change_crt_calib_ntuples' into feature/hlay_amoor_crt_tpc_matching_data --- sbndcode/CRT/CRTAna/crtana_sbnd.fcl | 1 + sbndcode/CRT/CRTAna/run_crtana.fcl | 12 +----------- sbndcode/CRT/CRTAna/run_crtana_data.fcl | 6 +----- 3 files changed, 3 insertions(+), 16 deletions(-) diff --git a/sbndcode/CRT/CRTAna/crtana_sbnd.fcl b/sbndcode/CRT/CRTAna/crtana_sbnd.fcl index 1cd48a612..5d0107236 100644 --- a/sbndcode/CRT/CRTAna/crtana_sbnd.fcl +++ b/sbndcode/CRT/CRTAna/crtana_sbnd.fcl @@ -12,6 +12,7 @@ crtana_sbnd: PropDelay: @local::sbnd_crtsim.DetSimParams.PropDelay TimeWalkNorm: @local::sbnd_crtsim.DetSimParams.TDelayNorm TimeWalkScale: @local::sbnd_crtsim.DetSimParams.TDelayScale + TruthMatch: false module_type: "CRTAnalysis" } diff --git a/sbndcode/CRT/CRTAna/run_crtana.fcl b/sbndcode/CRT/CRTAna/run_crtana.fcl index 0dff6f7a2..1f38761a9 100644 --- a/sbndcode/CRT/CRTAna/run_crtana.fcl +++ b/sbndcode/CRT/CRTAna/run_crtana.fcl @@ -21,17 +21,7 @@ physics: { analyzers: { - crtana: - { - module_type: "CRTAnalysis" - CRTGeoAlg: @local::crtgeoalg_sbnd - CRTBackTrackerAlg: @local::crtbacktrackeralg_sbnd - PEAttenuation: @local::sbnd_crtsim.DetSimParams.NpeScaleShift - PropDelay: @local::sbnd_crtsim.DetSimParams.PropDelay - TimeWalkNorm: @local::sbnd_crtsim.DetSimParams.TDelayNorm - TimeWalkScale: @local::sbnd_crtsim.DetSimParams.TDelayScale - TruthMatch: false - } + crtana: @local::crtana_sbnd } ana: [ crtana ] diff --git a/sbndcode/CRT/CRTAna/run_crtana_data.fcl b/sbndcode/CRT/CRTAna/run_crtana_data.fcl index b44faf920..b2ade27a9 100644 --- a/sbndcode/CRT/CRTAna/run_crtana_data.fcl +++ b/sbndcode/CRT/CRTAna/run_crtana_data.fcl @@ -8,8 +8,4 @@ services.BackTrackerService: @erase services.CRTChannelMapService: @local::crt_channel_map_standard services.CRTCalibService: @local::crt_calib_service -physics.analyzers.crtana.FEBDataModuleLabel: "crtdecoder" -physics.analyzers.crtana.DataMode: true -physics.analyzers.crtana.NoTPC: true -physics.analyzers.crtana.HasPTB: true -physics.analyzers.crtana.HasTDC: true +physics.analyzers.crtana: @local::crtana_data_sbnd From 53379c411c91d6e3fb38f124e00c944c89b3d8a0 Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Fri, 28 Feb 2025 12:57:51 -0600 Subject: [PATCH 15/18] Update analysis module to record more matched track info --- sbndcode/CRT/CRTAna/CRTAnalysis_module.cc | 201 +++++++++++++++++----- sbndcode/CRT/CRTAna/run_crtana.fcl | 8 +- 2 files changed, 162 insertions(+), 47 deletions(-) diff --git a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc index cd99d1f0b..b03fdd282 100644 --- a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc +++ b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc @@ -26,6 +26,7 @@ #include "larsim/Utils/TruthMatchUtils.h" #include "lardata/DetectorInfoServices/DetectorClocksService.h" +#include "lardata/DetectorInfoServices/DetectorPropertiesService.h" #include "sbnobj/SBND/CRT/FEBData.hh" #include "sbnobj/SBND/CRT/CRTStripHit.hh" @@ -39,6 +40,7 @@ #include "sbndcode/ChannelMaps/CRT/CRTChannelMapService.h" #include "sbndcode/CRT/CRTBackTracker/CRTBackTrackerAlg.h" #include "sbndcode/CRT/CRTUtils/CRTCommonUtils.h" +#include "sbndcode/CRT/CRTUtils/TPCGeoUtil.h" #include "sbndcode/Decoders/PTB/sbndptb.h" #include "sbndcode/Timing/SBNDRawTimingObj.h" @@ -284,33 +286,49 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer { std::vector _tr_truth_theta; std::vector _tr_truth_phi; - std::vector _tpc_start_x; - std::vector _tpc_start_y; - std::vector _tpc_start_z; - std::vector _tpc_end_x; - std::vector _tpc_end_y; - std::vector _tpc_end_z; - std::vector _tpc_dir_x; - std::vector _tpc_dir_y; - std::vector _tpc_dir_z; - std::vector _tpc_length; - std::vector _tpc_track_score; - std::vector _tpc_truth_trackid; - std::vector _tpc_truth_pdg; - std::vector _tpc_truth_energy; - std::vector _tpc_truth_time; - std::vector _tpc_sp_matchable; - std::vector _tpc_sp_matched; - std::vector _tpc_sp_good_match; - std::vector _tpc_sp_ts0; - std::vector _tpc_sp_ts1; - std::vector _tpc_sp_score; - std::vector _tpc_tr_matchable; - std::vector _tpc_tr_matched; - std::vector _tpc_tr_good_match; - std::vector _tpc_tr_ts0; - std::vector _tpc_tr_ts1; - std::vector _tpc_tr_score; + std::vector _tpc_start_x; + std::vector _tpc_start_y; + std::vector _tpc_start_z; + std::vector _tpc_end_x; + std::vector _tpc_end_y; + std::vector _tpc_end_z; + std::vector _tpc_start_dir_x; + std::vector _tpc_start_dir_y; + std::vector _tpc_start_dir_z; + std::vector _tpc_end_dir_x; + std::vector _tpc_end_dir_y; + std::vector _tpc_end_dir_z; + std::vector _tpc_length; + std::vector _tpc_track_score; + std::vector _tpc_truth_trackid; + std::vector _tpc_truth_pdg; + std::vector _tpc_truth_energy; + std::vector _tpc_truth_time; + std::vector _tpc_sp_matchable; + std::vector _tpc_sp_matched; + std::vector _tpc_sp_good_match; + std::vector _tpc_sp_xshift; + std::vector _tpc_sp_ts0; + std::vector _tpc_sp_ts1; + std::vector _tpc_sp_tagger; + std::vector _tpc_sp_nhits; + std::vector _tpc_sp_x; + std::vector _tpc_sp_y; + std::vector _tpc_sp_z; + std::vector _tpc_sp_score; + std::vector _tpc_tr_matchable; + std::vector _tpc_tr_matched; + std::vector _tpc_tr_good_match; + std::vector _tpc_tr_ts0; + std::vector _tpc_tr_ts1; + std::vector> _tpc_tr_taggers; + std::vector _tpc_tr_start_x; + std::vector _tpc_tr_start_y; + std::vector _tpc_tr_start_z; + std::vector _tpc_tr_end_x; + std::vector _tpc_tr_end_y; + std::vector _tpc_tr_end_z; + std::vector _tpc_tr_score; std::vector _ptb_hlt_trigger; std::vector _ptb_hlt_timestamp; @@ -556,18 +574,34 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fTree->Branch("tpc_end_x", "std::vector", &_tpc_end_x); fTree->Branch("tpc_end_y", "std::vector", &_tpc_end_y); fTree->Branch("tpc_end_z", "std::vector", &_tpc_end_z); - fTree->Branch("tpc_dir_x", "std::vector", &_tpc_dir_x); - fTree->Branch("tpc_dir_y", "std::vector", &_tpc_dir_y); - fTree->Branch("tpc_dir_z", "std::vector", &_tpc_dir_z); + fTree->Branch("tpc_start_dir_x", "std::vector", &_tpc_start_dir_x); + fTree->Branch("tpc_start_dir_y", "std::vector", &_tpc_start_dir_y); + fTree->Branch("tpc_start_dir_z", "std::vector", &_tpc_start_dir_z); + fTree->Branch("tpc_end_dir_x", "std::vector", &_tpc_end_dir_x); + fTree->Branch("tpc_end_dir_y", "std::vector", &_tpc_end_dir_y); + fTree->Branch("tpc_end_dir_z", "std::vector", &_tpc_end_dir_z); fTree->Branch("tpc_length", "std::vector", &_tpc_length); fTree->Branch("tpc_track_score", "std::vector", &_tpc_track_score); fTree->Branch("tpc_sp_matched", "std::vector", &_tpc_sp_matched); + fTree->Branch("tpc_sp_xshift", "std::vector", &_tpc_sp_xshift); fTree->Branch("tpc_sp_ts0", "std::vector", &_tpc_sp_ts0); fTree->Branch("tpc_sp_ts1", "std::vector", &_tpc_sp_ts1); + fTree->Branch("tpc_sp_tagger", "std::vector", &_tpc_sp_tagger); + fTree->Branch("tpc_sp_nhits", "std::vector", &_tpc_sp_nhits); + fTree->Branch("tpc_sp_x", "std::vector", &_tpc_sp_x); + fTree->Branch("tpc_sp_y", "std::vector", &_tpc_sp_y); + fTree->Branch("tpc_sp_z", "std::vector", &_tpc_sp_z); fTree->Branch("tpc_sp_score", "std::vector", &_tpc_sp_score); fTree->Branch("tpc_tr_matched", "std::vector", &_tpc_tr_matched); fTree->Branch("tpc_tr_ts0", "std::vector", &_tpc_tr_ts0); fTree->Branch("tpc_tr_ts1", "std::vector", &_tpc_tr_ts1); + fTree->Branch("tpc_tr_taggers", "std::vector>", &_tpc_tr_taggers); + fTree->Branch("tpc_tr_start_x", "std::vector", &_tpc_tr_start_x); + fTree->Branch("tpc_tr_start_y", "std::vector", &_tpc_tr_start_y); + fTree->Branch("tpc_tr_start_z", "std::vector", &_tpc_tr_start_z); + fTree->Branch("tpc_tr_end_x", "std::vector", &_tpc_tr_end_x); + fTree->Branch("tpc_tr_end_y", "std::vector", &_tpc_tr_end_y); + fTree->Branch("tpc_tr_end_z", "std::vector", &_tpc_tr_end_z); fTree->Branch("tpc_tr_score", "std::vector", &_tpc_tr_score); if(!fDataMode && fTruthMatch) { @@ -1547,9 +1581,12 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_end_x.resize(nTracks); _tpc_end_y.resize(nTracks); _tpc_end_z.resize(nTracks); - _tpc_dir_x.resize(nTracks); - _tpc_dir_y.resize(nTracks); - _tpc_dir_z.resize(nTracks); + _tpc_start_dir_x.resize(nTracks); + _tpc_start_dir_y.resize(nTracks); + _tpc_start_dir_z.resize(nTracks); + _tpc_end_dir_x.resize(nTracks); + _tpc_end_dir_y.resize(nTracks); + _tpc_end_dir_z.resize(nTracks); _tpc_length.resize(nTracks); _tpc_track_score.resize(nTracks); _tpc_truth_trackid.resize(nTracks); @@ -1559,14 +1596,27 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_sp_matchable.resize(nTracks); _tpc_sp_matched.resize(nTracks); _tpc_sp_good_match.resize(nTracks); + _tpc_sp_xshift.resize(nTracks); _tpc_sp_ts0.resize(nTracks); _tpc_sp_ts1.resize(nTracks); + _tpc_sp_tagger.resize(nTracks); + _tpc_sp_nhits.resize(nTracks); + _tpc_sp_x.resize(nTracks); + _tpc_sp_y.resize(nTracks); + _tpc_sp_z.resize(nTracks); _tpc_sp_score.resize(nTracks); _tpc_tr_matchable.resize(nTracks); _tpc_tr_matched.resize(nTracks); _tpc_tr_good_match.resize(nTracks); _tpc_tr_ts0.resize(nTracks); _tpc_tr_ts1.resize(nTracks); + _tpc_tr_taggers.resize(nTracks); + _tpc_tr_start_x.resize(nTracks); + _tpc_tr_start_y.resize(nTracks); + _tpc_tr_start_z.resize(nTracks); + _tpc_tr_end_x.resize(nTracks); + _tpc_tr_end_y.resize(nTracks); + _tpc_tr_end_z.resize(nTracks); _tpc_tr_score.resize(nTracks); art::FindOneP tracksToPFPs(TPCTrackHandle, e, fTPCTrackModuleLabel); @@ -1576,7 +1626,10 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: art::FindOneP spsToClusters(CRTSpacePointHandle, e, fCRTSpacePointModuleLabel); art::FindOneP tracksToTrackMatches(TPCTrackHandle, e, fCRTTrackMatchingModuleLabel); - const detinfo::DetectorClocksData clockData = art::ServiceHandle()->DataFor(e); + const detinfo::DetectorClocksData clockData = art::ServiceHandle()->DataFor(e); + const detinfo::DetectorPropertiesData detProp = art::ServiceHandle()->DataFor(e); + geo::GeometryCore const* geometryService = lar::providerFrom(); + int nActualTracks = 0; for(unsigned i = 0; i < nTracks; ++i) @@ -1599,10 +1652,17 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_end_y[nActualTracks] = end.Y(); _tpc_end_z[nActualTracks] = end.Z(); - const geo::Vector_t dir = track->StartDirection(); - _tpc_dir_x[nActualTracks] = dir.X(); - _tpc_dir_y[nActualTracks] = dir.Y(); - _tpc_dir_z[nActualTracks] = dir.Z(); + const std::pair dirs = CRTCommonUtils::AverageTrackDirections(track, 0.2); + + const geo::Vector_t start_dir = dirs.first; + _tpc_start_dir_x[nActualTracks] = start_dir.X(); + _tpc_start_dir_y[nActualTracks] = start_dir.Y(); + _tpc_start_dir_z[nActualTracks] = start_dir.Z(); + + const geo::Vector_t end_dir = dirs.first; + _tpc_end_dir_x[nActualTracks] = end_dir.X(); + _tpc_end_dir_y[nActualTracks] = end_dir.Y(); + _tpc_end_dir_z[nActualTracks] = end_dir.Z(); _tpc_length[nActualTracks] = track->Length(); @@ -1615,13 +1675,25 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: const art::Ptr spacepoint = tracksToSPMatches.at(track.key()); const art::Ptr crttrack = tracksToTrackMatches.at(track.key()); + const std::vector> trackHits = tracksToHits.at(track.key()); + if(spacepoint.isNonnull()) { - const anab::T0 spMatch = tracksToSPMatches.data(track.key()).ref(); + const anab::T0 spMatch = tracksToSPMatches.data(track.key()).ref(); + const art::Ptr cluster = spsToClusters.at(spacepoint.key()); + + const int driftDirection = TPCGeoUtil::DriftDirectionFromHits(geometryService, trackHits); + const double crtShiftingTime = fDataMode ? spacepoint->Ts0() * 1e-3 : spacepoint->Ts1() * 1e-3; _tpc_sp_matched[nActualTracks] = true; + _tpc_sp_xshift[nActualTracks] = driftDirection * crtShiftingTime * detProp.DriftVelocity(); _tpc_sp_ts0[nActualTracks] = spacepoint->Ts0(); _tpc_sp_ts1[nActualTracks] = spacepoint->Ts1(); + _tpc_sp_tagger[nActualTracks] = cluster->Tagger(); + _tpc_sp_nhits[nActualTracks] = cluster->NHits(); + _tpc_sp_x[nActualTracks] = spacepoint->X(); + _tpc_sp_y[nActualTracks] = spacepoint->Y(); + _tpc_sp_z[nActualTracks] = spacepoint->Z(); _tpc_sp_score[nActualTracks] = spMatch.TriggerConfidence(); } else @@ -1629,6 +1701,10 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_sp_matched[nActualTracks] = false; _tpc_sp_ts0[nActualTracks] = -std::numeric_limits::max(); _tpc_sp_ts1[nActualTracks] = -std::numeric_limits::max(); + _tpc_sp_tagger[nActualTracks] = -std::numeric_limits::max(); + _tpc_sp_x[nActualTracks] = -std::numeric_limits::max(); + _tpc_sp_y[nActualTracks] = -std::numeric_limits::max(); + _tpc_sp_z[nActualTracks] = -std::numeric_limits::max(); _tpc_sp_score[nActualTracks] = -std::numeric_limits::max(); } @@ -1640,18 +1716,40 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_tr_ts0[nActualTracks] = crttrack->Ts0(); _tpc_tr_ts1[nActualTracks] = crttrack->Ts1(); _tpc_tr_score[nActualTracks] = trackMatch.TriggerConfidence(); + + const std::set taggers = crttrack->Taggers(); + const geo::Point_t start = crttrack->Start(); + const geo::Point_t end = crttrack->End(); + + _tpc_tr_taggers[nActualTracks] = std::vector(); + + for(auto const& tagger : taggers) + _tpc_tr_taggers[nActualTracks].push_back(tagger); + + _tpc_tr_start_x[nActualTracks] = start.X(); + _tpc_tr_start_y[nActualTracks] = start.Y(); + _tpc_tr_start_z[nActualTracks] = start.Z(); + _tpc_tr_end_x[nActualTracks] = end.X(); + _tpc_tr_end_y[nActualTracks] = end.Y(); + _tpc_tr_end_z[nActualTracks] = end.Z(); } - else + else { _tpc_tr_matched[nActualTracks] = false; _tpc_tr_ts0[nActualTracks] = -std::numeric_limits::max(); _tpc_tr_ts1[nActualTracks] = -std::numeric_limits::max(); _tpc_tr_score[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_taggers[nActualTracks] = std::vector(); + _tpc_tr_start_x[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_start_y[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_start_z[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_end_x[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_end_y[nActualTracks] = -std::numeric_limits::max(); + _tpc_tr_end_z[nActualTracks] = -std::numeric_limits::max(); } if(!fDataMode && fTruthMatch) { - const std::vector> trackHits = tracksToHits.at(track.key()); const int trackid = fCRTBackTrackerAlg.RollUpID(TruthMatchUtils::TrueParticleIDFromTotalRecoHits(clockData,trackHits,true)); _tpc_truth_trackid[nActualTracks] = trackid; @@ -1708,9 +1806,12 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_end_x.resize(nActualTracks); _tpc_end_y.resize(nActualTracks); _tpc_end_z.resize(nActualTracks); - _tpc_dir_x.resize(nActualTracks); - _tpc_dir_y.resize(nActualTracks); - _tpc_dir_z.resize(nActualTracks); + _tpc_start_dir_x.resize(nActualTracks); + _tpc_start_dir_y.resize(nActualTracks); + _tpc_start_dir_z.resize(nActualTracks); + _tpc_end_dir_x.resize(nActualTracks); + _tpc_end_dir_y.resize(nActualTracks); + _tpc_end_dir_z.resize(nActualTracks); _tpc_length.resize(nActualTracks); _tpc_track_score.resize(nActualTracks); _tpc_truth_trackid.resize(nActualTracks); @@ -1719,12 +1820,24 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_truth_time.resize(nActualTracks); _tpc_sp_matchable.resize(nActualTracks); _tpc_sp_matched.resize(nActualTracks); + _tpc_sp_tagger.resize(nActualTracks); + _tpc_sp_nhits.resize(nActualTracks); + _tpc_sp_x.resize(nActualTracks); + _tpc_sp_y.resize(nActualTracks); + _tpc_sp_z.resize(nActualTracks); _tpc_sp_good_match.resize(nActualTracks); _tpc_sp_ts0.resize(nActualTracks); _tpc_sp_ts1.resize(nActualTracks); _tpc_sp_score.resize(nActualTracks); _tpc_tr_matchable.resize(nActualTracks); _tpc_tr_matched.resize(nActualTracks); + _tpc_tr_taggers.resize(nActualTracks); + _tpc_tr_start_x.resize(nActualTracks); + _tpc_tr_start_y.resize(nActualTracks); + _tpc_tr_start_z.resize(nActualTracks); + _tpc_tr_end_x.resize(nActualTracks); + _tpc_tr_end_y.resize(nActualTracks); + _tpc_tr_end_z.resize(nActualTracks); _tpc_tr_good_match.resize(nActualTracks); _tpc_tr_ts0.resize(nActualTracks); _tpc_tr_ts1.resize(nActualTracks); diff --git a/sbndcode/CRT/CRTAna/run_crtana.fcl b/sbndcode/CRT/CRTAna/run_crtana.fcl index 1f38761a9..3a1299d13 100644 --- a/sbndcode/CRT/CRTAna/run_crtana.fcl +++ b/sbndcode/CRT/CRTAna/run_crtana.fcl @@ -7,9 +7,11 @@ services: { TFileService: { fileName: "crtana_sbnd.root" } @table::sbnd_basic_services - ParticleInventoryService: @local::sbnd_particleinventoryservice - BackTrackerService: @local::sbnd_backtrackerservice - DetectorClocksService: @local::sbnd_detectorclocks + ParticleInventoryService: @local::sbnd_particleinventoryservice + BackTrackerService: @local::sbnd_backtrackerservice + DetectorClocksService: @local::sbnd_detectorclocks + LArPropertiesService: @local::sbnd_properties + DetectorPropertiesService: @local::sbnd_detproperties } source: From 557a93c90d7b35a2a06aea2368de72a5c981ef06 Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Fri, 7 Mar 2025 11:54:44 -0600 Subject: [PATCH 16/18] Store channel numbers in matching --- sbndcode/CRT/CRTAna/CRTAnalysis_module.cc | 118 ++++++++++++---------- 1 file changed, 66 insertions(+), 52 deletions(-) diff --git a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc index b03fdd282..885cce232 100644 --- a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc +++ b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc @@ -86,8 +86,10 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer { const art::FindOneP &spacePointsToClusters, const art::FindManyP &clustersToStripHits); void AnalyseTPCMatching(const art::Event &e, const art::Handle> &TPCTrackHandle, - const art::Handle> &CRTSpacePointModuleLabel, const art::Handle> &PFPHandle, - const std::map &spacePointRecoStatusMap = {}, const std::map> &trackRecoStatusMap = {}); + const art::Handle> &CRTSpacePointHandle, const art::Handle> &CRTClusterHandle, + const art::Handle> &PFPHandle, + const std::map &spacePointRecoStatusMap = {}, + const std::map> &trackRecoStatusMap = {}); private: @@ -286,49 +288,50 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer { std::vector _tr_truth_theta; std::vector _tr_truth_phi; - std::vector _tpc_start_x; - std::vector _tpc_start_y; - std::vector _tpc_start_z; - std::vector _tpc_end_x; - std::vector _tpc_end_y; - std::vector _tpc_end_z; - std::vector _tpc_start_dir_x; - std::vector _tpc_start_dir_y; - std::vector _tpc_start_dir_z; - std::vector _tpc_end_dir_x; - std::vector _tpc_end_dir_y; - std::vector _tpc_end_dir_z; - std::vector _tpc_length; - std::vector _tpc_track_score; - std::vector _tpc_truth_trackid; - std::vector _tpc_truth_pdg; - std::vector _tpc_truth_energy; - std::vector _tpc_truth_time; - std::vector _tpc_sp_matchable; - std::vector _tpc_sp_matched; - std::vector _tpc_sp_good_match; - std::vector _tpc_sp_xshift; - std::vector _tpc_sp_ts0; - std::vector _tpc_sp_ts1; - std::vector _tpc_sp_tagger; - std::vector _tpc_sp_nhits; - std::vector _tpc_sp_x; - std::vector _tpc_sp_y; - std::vector _tpc_sp_z; - std::vector _tpc_sp_score; - std::vector _tpc_tr_matchable; - std::vector _tpc_tr_matched; - std::vector _tpc_tr_good_match; - std::vector _tpc_tr_ts0; - std::vector _tpc_tr_ts1; - std::vector> _tpc_tr_taggers; - std::vector _tpc_tr_start_x; - std::vector _tpc_tr_start_y; - std::vector _tpc_tr_start_z; - std::vector _tpc_tr_end_x; - std::vector _tpc_tr_end_y; - std::vector _tpc_tr_end_z; - std::vector _tpc_tr_score; + std::vector _tpc_start_x; + std::vector _tpc_start_y; + std::vector _tpc_start_z; + std::vector _tpc_end_x; + std::vector _tpc_end_y; + std::vector _tpc_end_z; + std::vector _tpc_start_dir_x; + std::vector _tpc_start_dir_y; + std::vector _tpc_start_dir_z; + std::vector _tpc_end_dir_x; + std::vector _tpc_end_dir_y; + std::vector _tpc_end_dir_z; + std::vector _tpc_length; + std::vector _tpc_track_score; + std::vector _tpc_truth_trackid; + std::vector _tpc_truth_pdg; + std::vector _tpc_truth_energy; + std::vector _tpc_truth_time; + std::vector _tpc_sp_matchable; + std::vector _tpc_sp_matched; + std::vector> _tpc_sp_channel_set; + std::vector _tpc_sp_good_match; + std::vector _tpc_sp_xshift; + std::vector _tpc_sp_ts0; + std::vector _tpc_sp_ts1; + std::vector _tpc_sp_tagger; + std::vector _tpc_sp_nhits; + std::vector _tpc_sp_x; + std::vector _tpc_sp_y; + std::vector _tpc_sp_z; + std::vector _tpc_sp_score; + std::vector _tpc_tr_matchable; + std::vector _tpc_tr_matched; + std::vector _tpc_tr_good_match; + std::vector _tpc_tr_ts0; + std::vector _tpc_tr_ts1; + std::vector> _tpc_tr_taggers; + std::vector _tpc_tr_start_x; + std::vector _tpc_tr_start_y; + std::vector _tpc_tr_start_z; + std::vector _tpc_tr_end_x; + std::vector _tpc_tr_end_y; + std::vector _tpc_tr_end_z; + std::vector _tpc_tr_score; std::vector _ptb_hlt_trigger; std::vector _ptb_hlt_timestamp; @@ -583,6 +586,7 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fTree->Branch("tpc_length", "std::vector", &_tpc_length); fTree->Branch("tpc_track_score", "std::vector", &_tpc_track_score); fTree->Branch("tpc_sp_matched", "std::vector", &_tpc_sp_matched); + fTree->Branch("tpc_sp_channel_set", "std::vector>", &_tpc_sp_channel_set); fTree->Branch("tpc_sp_xshift", "std::vector", &_tpc_sp_xshift); fTree->Branch("tpc_sp_ts0", "std::vector", &_tpc_sp_ts0); fTree->Branch("tpc_sp_ts1", "std::vector", &_tpc_sp_ts1); @@ -850,7 +854,7 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) std::cout << "CRTAnalysis -- Processing TPC" << std::endl; // Fill TPC matching variables - AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, PFPHandle); + AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, CRTClusterHandle, PFPHandle); // Fill the Tree fTree->Fill(); @@ -911,7 +915,7 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) std::cout << "CRTAnalysis -- Processing TPC" << std::endl; // Fill TPC matching variables - AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, PFPHandle, spacePointRecoStatusMap, trackRecoStatusMap); + AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, CRTClusterHandle, PFPHandle, spacePointRecoStatusMap, trackRecoStatusMap); } else { @@ -919,7 +923,7 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) std::cout << "CRTAnalysis -- Processing TPC" << std::endl; // Fill TPC matching variables - AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, PFPHandle); + AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, CRTClusterHandle, PFPHandle); } // Fill the Tree fTree->Fill(); @@ -1567,8 +1571,10 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTTracks(const art::Event &e, const std::ve } void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art::Handle> &TPCTrackHandle, - const art::Handle> &CRTSpacePointHandle, const art::Handle> &PFPHandle, - const std::map &spacePointRecoStatusMap, const std::map> &trackRecoStatusMap) + const art::Handle> &CRTSpacePointHandle, const art::Handle> &CRTClusterHandle, + const art::Handle> &PFPHandle, + const std::map &spacePointRecoStatusMap, + const std::map> &trackRecoStatusMap) { std::vector> TPCTrackVec; art::fill_ptr_vector(TPCTrackVec, TPCTrackHandle); @@ -1595,6 +1601,7 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_truth_time.resize(nTracks); _tpc_sp_matchable.resize(nTracks); _tpc_sp_matched.resize(nTracks); + _tpc_sp_channel_set.resize(nTracks); _tpc_sp_good_match.resize(nTracks); _tpc_sp_xshift.resize(nTracks); _tpc_sp_ts0.resize(nTracks); @@ -1624,6 +1631,7 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: art::FindManyP tracksToHits(TPCTrackHandle, e, fTPCTrackModuleLabel); art::FindOneP tracksToSPMatches(TPCTrackHandle, e, fCRTSpacePointMatchingModuleLabel); art::FindOneP spsToClusters(CRTSpacePointHandle, e, fCRTSpacePointModuleLabel); + art::FindManyP clustersToStripHits(CRTClusterHandle, e, fCRTClusterModuleLabel); art::FindOneP tracksToTrackMatches(TPCTrackHandle, e, fCRTTrackMatchingModuleLabel); const detinfo::DetectorClocksData clockData = art::ServiceHandle()->DataFor(e); @@ -1679,8 +1687,9 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: if(spacepoint.isNonnull()) { - const anab::T0 spMatch = tracksToSPMatches.data(track.key()).ref(); - const art::Ptr cluster = spsToClusters.at(spacepoint.key()); + const anab::T0 spMatch = tracksToSPMatches.data(track.key()).ref(); + const art::Ptr cluster = spsToClusters.at(spacepoint.key()); + const std::vector> striphits = clustersToStripHits.at(cluster.key()); const int driftDirection = TPCGeoUtil::DriftDirectionFromHits(geometryService, trackHits); const double crtShiftingTime = fDataMode ? spacepoint->Ts0() * 1e-3 : spacepoint->Ts1() * 1e-3; @@ -1695,6 +1704,11 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: _tpc_sp_y[nActualTracks] = spacepoint->Y(); _tpc_sp_z[nActualTracks] = spacepoint->Z(); _tpc_sp_score[nActualTracks] = spMatch.TriggerConfidence(); + + _tpc_sp_channel_set[nActualTracks] = std::vector(); + + for(auto const &striphit : striphits) + _tpc_sp_channel_set[nActualTracks].push_back(striphit->Channel()); } else { From 7509b1f9fa8d73ff9b67426ce9e8bb3911a64b81 Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Fri, 21 Mar 2025 03:11:22 -0500 Subject: [PATCH 17/18] Remove geometry changes and extra fcls for testing --- .../CRTTPCMatching/run_crttpcmatching_data.fcl | 17 ----------------- 1 file changed, 17 deletions(-) delete mode 100644 sbndcode/CRT/CRTTPCMatching/run_crttpcmatching_data.fcl diff --git a/sbndcode/CRT/CRTTPCMatching/run_crttpcmatching_data.fcl b/sbndcode/CRT/CRTTPCMatching/run_crttpcmatching_data.fcl deleted file mode 100644 index 34ddd634b..000000000 --- a/sbndcode/CRT/CRTTPCMatching/run_crttpcmatching_data.fcl +++ /dev/null @@ -1,17 +0,0 @@ -#include "run_crttpcmatching.fcl" - -process_name: CRTTPCMatching - -physics.producers.crtspacepointmatching: @local::crtspacepointmatchproducer_data_sbnd -physics.producers.crttrackmatching: @local::crttrackmatchproducer_data_sbnd - -physics.producers.crtspacepointmatchingevd: @local::crtspacepointmatchproducer_data_sbnd -physics.producers.crtspacepointmatchingevd.MakeEventDisplays: true -physics.producers.crtspacepointmatchingevd.CRTEventDisplayAlg.DrawTracks: false -physics.producers.crtspacepointmatchingevd.CRTEventDisplayAlg.DisplayMatchScore: true - -physics.producers.crttrackmatchingevd: @local::crttrackmatchproducer_data_sbnd -physics.producers.crttrackmatchingevd.MakeEventDisplays: true -physics.producers.crttrackmatchingevd.CRTEventDisplayAlg.DisplayMatchScore: true - -physics.reco: [ crtspacepointmatching, crttrackmatching, crtspacepointmatchingevd, crttrackmatchingevd ] From e8849d210e229c05df17f649220906b911106b10 Mon Sep 17 00:00:00 2001 From: Henry Lay Date: Fri, 21 Mar 2025 05:34:08 -0500 Subject: [PATCH 18/18] Cleanup --- sbndcode/Decoders/run_decoders_job_no_choppy.fcl | 7 +++++++ sbndcode/JobConfigurations/standard/reco/reco2_data.fcl | 3 ++- 2 files changed, 9 insertions(+), 1 deletion(-) create mode 100644 sbndcode/Decoders/run_decoders_job_no_choppy.fcl diff --git a/sbndcode/Decoders/run_decoders_job_no_choppy.fcl b/sbndcode/Decoders/run_decoders_job_no_choppy.fcl new file mode 100644 index 000000000..167d4f695 --- /dev/null +++ b/sbndcode/Decoders/run_decoders_job_no_choppy.fcl @@ -0,0 +1,7 @@ +#include "run_decoders_job.fcl" + +physics.end_paths: [ stream1 ] + +# Stops the output of the choppy decoded file, only outputs the filtered decoded file +outputs.out2: @erase +physics.stream2: @erase diff --git a/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl b/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl index b160f68ea..2298b988d 100644 --- a/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl +++ b/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl @@ -19,10 +19,11 @@ physics.producers: crttracks: @local::crttrackproducer_data_sbnd crtspacepointmatching: @local::crtspacepointmatchproducer_data_sbnd crttrackmatching: @local::crttrackmatchproducer_data_sbnd + opt0finder: @local::sbnd_opt0_finder_data } physics.reco2: [ pandora, pandoraTrack, pandoraShower, pandoraShowerSBN, pandoraCaloData, pandoraPidData, caloskimCalorimetry, - crtclustering, crtspacepoints, crttracks, crtspacepointmatching, crttrackmatching, cvn ] + crtclustering, crtspacepoints, crttracks, crtspacepointmatching, crttrackmatching, cvn, opt0finder] #The next 3 lines need to be commented out once data use pandoraSCE by default. physics.producers.cvn.SliceLabel: "pandora"