diff --git a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc index 718ad3540..885cce232 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" @@ -84,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: @@ -97,7 +101,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; @@ -284,33 +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_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_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; @@ -348,13 +369,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; @@ -422,7 +444,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); @@ -445,7 +467,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); @@ -478,7 +500,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); @@ -523,7 +545,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); @@ -555,29 +577,46 @@ 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); - if(!fDataMode) + 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); + 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) { fTree->Branch("tpc_truth_trackid", "std::vector", &_tpc_truth_trackid); fTree->Branch("tpc_truth_pdg", "std::vector", &_tpc_truth_pdg); 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); } } @@ -632,22 +671,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 +704,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 +723,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 +740,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 +756,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 +772,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 +794,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,8 +850,11 @@ 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); + AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, CRTClusterHandle, PFPHandle); // Fill the Tree fTree->Fill(); @@ -795,6 +862,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 +878,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 +894,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; - // Fill TrueDeposit variables - AnalyseTrueDepositsPerTagger(spacePointRecoStatusMap); + if(fTruthMatch) + { + // Get Map of TrueDeposits per tagger from BackTracker + std::map spacePointRecoStatusMap = fCRTBackTrackerAlg.GetSpacePointRecoStatusMap(); - // Get Map of TrueDeposits from BackTracker - std::map> trackRecoStatusMap = fCRTBackTrackerAlg.GetTrackRecoStatusMap(); + // Fill TrueDeposit variables + AnalyseTrueDepositsPerTagger(spacePointRecoStatusMap); - // Fill TrueDeposit variables - AnalyseTrueDeposits(trackRecoStatusMap); + // Get Map of TrueDeposits from BackTracker + std::map> trackRecoStatusMap = fCRTBackTrackerAlg.GetTrackRecoStatusMap(); - // Fill TPC matching variables - AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, PFPHandle, spacePointRecoStatusMap, trackRecoStatusMap); + // Fill TrueDeposit variables + AnalyseTrueDeposits(trackRecoStatusMap); + if(fDebug) + std::cout << "CRTAnalysis -- Processing TPC" << std::endl; + + // Fill TPC matching variables + AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, CRTClusterHandle, PFPHandle, spacePointRecoStatusMap, trackRecoStatusMap); + } + else + { + if(fDebug) + std::cout << "CRTAnalysis -- Processing TPC" << std::endl; + + // Fill TPC matching variables + AnalyseTPCMatching(e, TPCTrackHandle, CRTSpacePointHandle, CRTClusterHandle, PFPHandle); + } // Fill the Tree fTree->Fill(); } @@ -1082,7 +1171,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 +1301,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 +1534,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; @@ -1482,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); @@ -1496,9 +1587,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); @@ -1507,15 +1601,29 @@ 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); _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); @@ -1523,9 +1631,13 @@ 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); + 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) @@ -1548,10 +1660,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(); @@ -1561,9 +1680,90 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art:: else _tpc_track_score[nActualTracks] = -std::numeric_limits::max(); - if(!fDataMode) + 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 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; + + _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(); + + _tpc_sp_channel_set[nActualTracks] = std::vector(); + + for(auto const &striphit : striphits) + _tpc_sp_channel_set[nActualTracks].push_back(striphit->Channel()); + } + else + { + _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(); + } + + if(crttrack.isNonnull()) + { + 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(); + + 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 + { + _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; @@ -1587,49 +1787,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(); } } @@ -1642,9 +1820,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); @@ -1653,12 +1834,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/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 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: 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" diff --git a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.cc index 4ef671723..6feaa0136 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,9 @@ namespace sbnd::crt { fClusterLabel = config.ClusterLabel(); fSpacePointLabel = config.SpacePointLabel(); fTrackLabel = config.TrackLabel(); + fTPCSpacePointMatchLabel = config.TPCSpacePointMatchLabel(); + fTPCTrackMatchLabel = config.TPCTrackMatchLabel(); + fTPCTrackLabel = config.TPCTrackLabel(); fSaveRoot = config.SaveRoot(); fSaveViews = config.SaveViews(); @@ -39,6 +47,9 @@ namespace sbnd::crt { fDrawClusters = config.DrawClusters(); fDrawSpacePoints = config.DrawSpacePoints(); fDrawTracks = config.DrawTracks(); + fDrawTPCMatching = config.DrawTPCMatching(); + fOnlyDrawMatched = config.OnlyDrawMatched(); + fDisplayMatchScore = config.DisplayMatchScore(); fChoseTaggers = config.ChoseTaggers(); fChosenTaggers = config.ChosenTaggers(); @@ -58,6 +69,7 @@ namespace sbnd::crt { fClusterColourInterval = config.ClusterColourInterval(); fSpacePointColour = config.SpacePointColour(); fTrackColour = config.TrackColour(); + fTPCMatchColour = config.TPCMatchColour(); fUseTs0 = config.UseTs0(); fMinTime = config.MinTime(); @@ -100,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; @@ -142,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); @@ -152,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; @@ -163,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); } @@ -182,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()) { @@ -203,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); @@ -217,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); } @@ -239,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); } @@ -253,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) { @@ -461,22 +457,122 @@ 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(); + + 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()); + 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; + { + 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; @@ -499,7 +595,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()) @@ -509,6 +605,96 @@ 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(); + + 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()); + 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(); @@ -537,17 +723,25 @@ namespace sbnd::crt { line->SetLineColor(fTrackColour); line->SetLineWidth(fLineWidth); - 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((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(); + + if(foundMatch) + std::cout << " and found TPC track match with score: " << matchScore << std::endl; + else + std::cout << std::endl; + } } } diff --git a/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h b/sbndcode/CRT/CRTEventDisplay/CRTEventDisplayAlg.h index 9429bd033..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" @@ -87,6 +89,15 @@ namespace sbnd::crt { fhicl::Atom TrackLabel { Name("TrackLabel") }; + fhicl::Atom TPCSpacePointMatchLabel { + Name("TPCSpacePointMatchLabel") + }; + fhicl::Atom TPCTrackMatchLabel { + Name("TPCTrackMatchLabel") + }; + fhicl::Atom TPCTrackLabel { + Name("TPCTrackLabel") + }; fhicl::Atom SaveRoot { Name("SaveRoot"), @@ -131,6 +142,15 @@ namespace sbnd::crt { fhicl::Atom DrawTracks { Name("DrawTracks") }; + fhicl::Atom DrawTPCMatching { + Name("DrawTPCMatching") + }; + fhicl::Atom OnlyDrawMatched { + Name("OnlyDrawMatched") + }; + fhicl::Atom DisplayMatchScore { + Name("DisplayMatchScore") + }; fhicl::Atom ChoseTaggers { Name("ChoseTaggers") @@ -182,6 +202,9 @@ namespace sbnd::crt { fhicl::Atom TrackColour { Name("TrackColour") }; + fhicl::Atom TPCMatchColour { + Name("TPCMatchColour") + }; fhicl::Atom UseTs0 { Name ("UseTs0") @@ -222,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); @@ -247,6 +273,9 @@ namespace sbnd::crt { art::InputTag fClusterLabel; art::InputTag fSpacePointLabel; art::InputTag fTrackLabel; + art::InputTag fTPCSpacePointMatchLabel; + art::InputTag fTPCTrackMatchLabel; + art::InputTag fTPCTrackLabel; bool fSaveRoot; bool fSaveViews; @@ -263,6 +292,9 @@ namespace sbnd::crt { bool fDrawClusters; bool fDrawSpacePoints; bool fDrawTracks; + bool fDrawTPCMatching; + bool fOnlyDrawMatched; + bool fDisplayMatchScore; bool fChoseTaggers; std::vector fChosenTaggers; @@ -282,6 +314,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 5eb31515e..90a2200a0 100644 --- a/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl +++ b/sbndcode/CRT/CRTEventDisplay/crteventdisplayalg_sbnd.fcl @@ -11,28 +11,34 @@ crteventdisplayalg_sbnd: MC: true - SimLabel: "largeant" - SimDepositLabel: "genericcrt" - StripHitLabel: "crtstrips" - ClusterLabel: "crtclustering" - SpacePointLabel: "crtspacepoints" - TrackLabel: "crttracks" + SimLabel: "largeant" + SimDepositLabel: "genericcrt" + StripHitLabel: "crtstrips" + ClusterLabel: "crtclustering" + SpacePointLabel: "crtspacepoints" + TrackLabel: "crttracks" + 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 + 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 @@ -63,6 +69,7 @@ crteventdisplayalg_sbnd: ClusterColourInterval: 20 SpacePointColour: 6 TrackColour: 9 + TPCMatchColour: 8 UseTs0: false MinTime: -10000000000 @@ -81,11 +88,13 @@ 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.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 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/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/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/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; 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 77f35b4bc..a02373416 100644 --- a/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl +++ b/sbndcode/CRT/CRTTPCMatching/crttpcmatchingproducers_sbnd.fcl @@ -1,47 +1,67 @@ +#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.25 + DCALimit: 70. + MinTPCTrackLength: 10. + DirMethod: 2 + DCAuseBox: false + DCAoverLength: false + PECut: 60. + MaxUncert: 20. + TPCTrackLabel: "pandoraTrack" + CRTSpacePointLabel: "crtspacepoints" + UseTs0: false } +crtspacepointmatchalg_data_sbnd: @local::crtspacepointmatchalg_sbnd +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.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 +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_data_sbnd +crttrackmatchproducer_data_sbnd.CRTEventDisplayAlg: @local::crteventdisplayalg_sbnd_data + END_PROLOG 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); } } 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 11ba2052e..14803d26f 100644 --- a/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl +++ b/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl @@ -15,15 +15,18 @@ services: physics.producers: { @table::physics.producers - crtclustering: @local::crtclusterproducer_data_sbnd - crtspacepoints: @local::crtspacepointproducer_data_sbnd - crttracks: @local::crttrackproducer_data_sbnd + 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 + opt0finder: @local::sbnd_opt0_finder_data tpcpmtbarycentermatching: @local::TPCPMTBarycenterMatchProducer - opt0finder: @local::sbnd_opt0_finder_data + opt0finder: @local::sbnd_opt0_finder_data } physics.reco2: [ pandora, pandoraTrack, pandoraShower, pandoraShowerSBN, pandoraCaloData, pandoraPidData, caloskimCalorimetry, - crtclustering, crtspacepoints, crttracks, cvn, opt0finder, tpcpmtbarycentermatching] + crtclustering, crtspacepoints, crttracks, crtspacepointmatching, crttrackmatching, cvn, opt0finder, tpcpmtbarycentermatching ] #The next 3 lines need to be commented out once data use pandoraSCE by default. physics.producers.cvn.SliceLabel: "pandora"