diff --git a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc index 948d3fd4c..29f6614c9 100644 --- a/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc +++ b/sbndcode/CRT/CRTAna/CRTAnalysis_module.cc @@ -33,6 +33,7 @@ #include "sbnobj/SBND/CRT/CRTCluster.hh" #include "sbnobj/SBND/CRT/CRTSpacePoint.hh" #include "sbnobj/SBND/CRT/CRTTrack.hh" +#include "sbnobj/SBND/CRT/CRTBlob.hh" #include "sbnobj/SBND/Timing/DAQTimestamp.hh" #include "sbndcode/Geometry/GeometryWrappers/CRTGeoService.h" @@ -87,6 +88,8 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer { void AnalyseCRTTracks(const art::Event &e, const std::vector> &CRTTrackVec, const art::FindManyP &tracksToSpacePoints, const art::FindOneP &spacePointsToClusters, const art::FindManyP &clustersToStripHits); + void AnalyseCRTBlobs(std::vector> &CRTBlobVec); + void AnalyseTPCMatching(const art::Event &e, const art::Handle> &TPCTrackHandle, const art::Handle> &CRTSpacePointHandle, const art::Handle> &CRTClusterHandle, const art::Handle> &PFPHandle, @@ -104,10 +107,10 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer { CRTBackTrackerAlg fCRTBackTrackerAlg; std::string fMCParticleModuleLabel, fSimDepositModuleLabel, fFEBDataModuleLabel, fCRTStripHitModuleLabel, - fCRTClusterModuleLabel, fCRTSpacePointModuleLabel, fCRTTrackModuleLabel, fTPCTrackModuleLabel, + fCRTClusterModuleLabel, fCRTSpacePointModuleLabel, fCRTTrackModuleLabel, fCRTBlobModuleLabel, fTPCTrackModuleLabel, fCRTSpacePointMatchingModuleLabel, fCRTTrackMatchingModuleLabel, fPFPModuleLabel, fPTBModuleLabel, fTDCModuleLabel, fTimingReferenceModuleLabel; - bool fDebug, fDataMode, fNoTPC, fHasPTB, fHasTDC, fTruthMatch; + bool fDebug, fDataMode, fNoTPC, fHasPTB, fHasTDC, fHasBlobs, fTruthMatch; //! Adding some of the reco parameters to save corrections double fPEAttenuation, fTimeWalkNorm, fTimeWalkScale, fPropDelay; @@ -296,6 +299,16 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer { std::vector _tr_truth_theta; std::vector _tr_truth_phi; + // crt blob information + std::vector _bl_ts0; + std::vector _bl_ets0; + std::vector _bl_ts1; + std::vector _bl_ets1; + std::vector _bl_pe; + std::vector _bl_nsps; + std::vector> _bl_nsps_per_tagger; + + // tpc track information (including crt matching) std::vector _tpc_start_x; std::vector _tpc_start_y; std::vector _tpc_start_z; @@ -341,12 +354,14 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer { std::vector _tpc_tr_end_z; std::vector _tpc_tr_score; + // ptb information (trigger board) std::vector _ptb_hlt_trigger; std::vector _ptb_hlt_timestamp; std::vector _ptb_llt_trigger; std::vector _ptb_llt_timestamp; + // spec tdc information (timing board) std::vector _tdc_channel; std::vector _tdc_timestamp; std::vector _tdc_offset; @@ -364,6 +379,7 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fCRTClusterModuleLabel = p.get("CRTClusterModuleLabel", "crtclustering"); fCRTSpacePointModuleLabel = p.get("CRTSpacePointModuleLabel", "crtspacepoints"); fCRTTrackModuleLabel = p.get("CRTTrackModuleLabel", "crttracks"); + fCRTBlobModuleLabel = p.get("CRTBlobModuleLabel", "crtblobs"); fTPCTrackModuleLabel = p.get("TPCTrackModuleLabel", "pandoraSCETrack"); fCRTSpacePointMatchingModuleLabel = p.get("CRTSpacePointMatchingModuleLabel", "crtspacepointmatchingSCE"); fCRTTrackMatchingModuleLabel = p.get("CRTTrackMatchingModuleLabel", "crttrackmatchingSCE"); @@ -376,8 +392,8 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fNoTPC = p.get("NoTPC", false); fHasPTB = p.get("HasPTB", false); fHasTDC = p.get("HasTDC", false); + fHasBlobs = p.get("HasBlobs", 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); @@ -578,6 +594,17 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p) fTree->Branch("tr_truth_phi", "std::vector", &_tr_truth_phi); } + if(fHasBlobs) + { + fTree->Branch("bl_ts0", "std::vector", &_bl_ts0); + fTree->Branch("bl_ets0", "std::vector", &_bl_ets0); + fTree->Branch("bl_ts1", "std::vector", &_bl_ts1); + fTree->Branch("bl_ets1", "std::vector", &_bl_ets1); + fTree->Branch("bl_pe", "std::vector", &_bl_pe); + fTree->Branch("bl_nsps", "std::vector", &_bl_nsps); + fTree->Branch("bl_nsps_per_tagger", "std::vector>", &_bl_nsps_per_tagger); + } + if(!fNoTPC) { fTree->Branch("tpc_start_x", "std::vector", &_tpc_start_x); @@ -835,6 +862,23 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e) // Fill CRTTrack variables AnalyseCRTTracks(e, CRTTrackVec, tracksToSpacePoints, spacepointsToClusters, clustersToStripHits); + if(fHasBlobs) + { + // Get CRTBlobs + art::Handle> CRTBlobHandle; + e.getByLabel(fCRTBlobModuleLabel, CRTBlobHandle); + if(!CRTBlobHandle.isValid()){ + std::cout << "CRTBlob product " << fCRTBlobModuleLabel << " not found..." << std::endl; + throw std::exception(); + } + + std::vector> CRTBlobVec; + art::fill_ptr_vector(CRTBlobVec, CRTBlobHandle); + + // Fill CRTBlob variables + AnalyseCRTBlobs(CRTBlobVec); + } + if(fNoTPC) { // Fill the Tree @@ -1585,6 +1629,34 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTTracks(const art::Event &e, const std::ve } } +void sbnd::crt::CRTAnalysis::AnalyseCRTBlobs(std::vector> &CRTBlobVec) +{ + const unsigned nBlobs = CRTBlobVec.size(); + + _bl_ts0.resize(nBlobs); + _bl_ets0.resize(nBlobs); + _bl_ts1.resize(nBlobs); + _bl_ets1.resize(nBlobs); + _bl_pe.resize(nBlobs); + _bl_nsps.resize(nBlobs); + _bl_nsps_per_tagger.resize(nBlobs, std::vector(7)); + + for(unsigned i = 0; i < nBlobs; ++i) + { + const auto blob = CRTBlobVec[i]; + + _bl_ts0[i] = blob->Ts0(); + _bl_ets0[i] = blob->Ts0Err(); + _bl_ts1[i] = blob->Ts1(); + _bl_ets1[i] = blob->Ts1Err(); + _bl_pe[i] = blob->PE(); + _bl_nsps[i] = blob->TotalSpacePoints(); + + for(unsigned j = 0; j < 7; ++j) + _bl_nsps_per_tagger[i][j] = blob->SpacePointsInTagger((CRTTagger)j); + } +} + void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art::Handle> &TPCTrackHandle, const art::Handle> &CRTSpacePointHandle, const art::Handle> &CRTClusterHandle, const art::Handle> &PFPHandle, diff --git a/sbndcode/CRT/CRTReco/CMakeLists.txt b/sbndcode/CRT/CRTReco/CMakeLists.txt index 3ea86c28c..ad5cfc7cb 100644 --- a/sbndcode/CRT/CRTReco/CMakeLists.txt +++ b/sbndcode/CRT/CRTReco/CMakeLists.txt @@ -37,4 +37,10 @@ simple_plugin( Eigen3::Eigen ) +simple_plugin( + CRTBlobProducer module + lardata::Utilities + sbnobj::SBND_CRT +) + install_fhicl() diff --git a/sbndcode/CRT/CRTReco/CRTBlobProducer_module.cc b/sbndcode/CRT/CRTReco/CRTBlobProducer_module.cc new file mode 100644 index 000000000..0ea91e935 --- /dev/null +++ b/sbndcode/CRT/CRTReco/CRTBlobProducer_module.cc @@ -0,0 +1,219 @@ +//////////////////////////////////////////////////////////////////////// +// Class: CRTBlobProducer +// Plugin Type: producer +// File: CRTBlobProducer_module.cc +// +// Author: Henry Lay (h.lay@sheffield.ac.uk) +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "canvas/Persistency/Common/FindOneP.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "lardata/Utilities/AssociationUtil.h" + +#include "sbnobj/SBND/CRT/CRTCluster.hh" +#include "sbnobj/SBND/CRT/CRTSpacePoint.hh" +#include "sbnobj/SBND/CRT/CRTBlob.hh" + +#include + +namespace sbnd::crt { + class CRTBlobProducer; +} + +class sbnd::crt::CRTBlobProducer : public art::EDProducer { +public: + explicit CRTBlobProducer(fhicl::ParameterSet const& p); + + CRTBlobProducer(CRTBlobProducer const&) = delete; + CRTBlobProducer(CRTBlobProducer&&) = delete; + CRTBlobProducer& operator=(CRTBlobProducer const&) = delete; + CRTBlobProducer& operator=(CRTBlobProducer&&) = delete; + + void produce(art::Event& e) override; + + void OrderSpacePoints(std::vector> &spacePointVec, const art::FindOneP &spacePointsToCluster); + + std::vector>> CreateBlobCandidates(const std::vector> &spacePointVec, + const art::FindOneP &spacePointsToCluster); + + void TimeErrorCalculator(const std::vector ×, double &mean, double &err); + + void OrderBlobCandidates(std::vector>> &blobCandidates); + + std::vector>> ChoseBlobs(std::vector>> &blobCandidates); + +private: + + std::string fCRTSpacePointModuleLabel; + double fCoincidenceTimeRequirement; + bool fUseTs0; +}; + + +sbnd::crt::CRTBlobProducer::CRTBlobProducer(fhicl::ParameterSet const& p) + : EDProducer{p} + , fCRTSpacePointModuleLabel(p.get("CRTSpacePointModuleLabel")) + , fCoincidenceTimeRequirement(p.get("CoincidenceTimeRequirement")) + , fUseTs0(p.get("UseTs0")) +{ + produces>(); + produces>(); +} + +void sbnd::crt::CRTBlobProducer::produce(art::Event& e) +{ + auto blobVec = std::make_unique>(); + auto blobSpacePointAssn = std::make_unique>(); + + art::Handle> CRTSpacePointHandle; + e.getByLabel(fCRTSpacePointModuleLabel, CRTSpacePointHandle); + + std::vector> CRTSpacePointVec; + art::fill_ptr_vector(CRTSpacePointVec, CRTSpacePointHandle); + + art::FindOneP spacePointsToCluster(CRTSpacePointHandle, e, fCRTSpacePointModuleLabel); + + OrderSpacePoints(CRTSpacePointVec, spacePointsToCluster); + + std::vector>> blobCandidates = CreateBlobCandidates(CRTSpacePointVec, spacePointsToCluster); + + OrderBlobCandidates(blobCandidates); + + std::vector>> chosenBlobs = ChoseBlobs(blobCandidates); + + for(auto const& [blob, spIDs] : chosenBlobs) + { + blobVec->push_back(blob); + + for(auto const& spID : spIDs) + util::CreateAssn(*this, e, *blobVec, CRTSpacePointVec[spID], *blobSpacePointAssn); + } + + e.put(std::move(blobVec)); + e.put(std::move(blobSpacePointAssn)); +} + +void sbnd::crt::CRTBlobProducer::OrderSpacePoints(std::vector> &spacePointVec, const art::FindOneP &spacePointsToCluster) +{ + std::sort(spacePointVec.begin(), spacePointVec.end(), + [&](const art::Ptr &a, const art::Ptr &b) -> bool { + if(fUseTs0) + return a->Ts0() < b->Ts0(); + else + return a->Ts1() < b->Ts1(); + }); +} + +std::vector>> sbnd::crt::CRTBlobProducer::CreateBlobCandidates(const std::vector> &spacePointVec, + const art::FindOneP &spacePointsToCluster) +{ + std::vector>> candidates; + + for(unsigned i = 0; i < spacePointVec.size(); ++i) + { + const art::Ptr primarySpacePoint = spacePointVec[i]; + const art::Ptr primaryCluster = spacePointsToCluster.at(primarySpacePoint.key()); + + std::set used_spacepoints = { i }; + std::vector t0s = { primarySpacePoint->Ts0() }; + std::vector t1s = { primarySpacePoint->Ts1() }; + std::vector pes = { primarySpacePoint->PE() }; + + std::map taggers; + for(int i = 0; i < 7; ++i) + taggers[(CRTTagger)i] = 0; + ++taggers[primaryCluster->Tagger()]; + + for(unsigned ii = i+1; ii < spacePointVec.size(); ++ii) + { + const art::Ptr secondarySpacePoint = spacePointVec[ii]; + const art::Ptr secondaryCluster = spacePointsToCluster.at(secondarySpacePoint.key()); + + const double tdiff_prim_sec = fUseTs0 ? secondarySpacePoint->Ts0() - primarySpacePoint->Ts0() : secondarySpacePoint->Ts1() - primarySpacePoint->Ts1(); + + if(tdiff_prim_sec > fCoincidenceTimeRequirement) + break; + + used_spacepoints.insert(ii); + t0s.push_back(secondarySpacePoint->Ts0()); + t1s.push_back(secondarySpacePoint->Ts1()); + pes.push_back(secondarySpacePoint->PE()); + ++taggers[secondaryCluster->Tagger()]; + } + + double t0, et0; + TimeErrorCalculator(t0s, t0, et0); + + double t1, et1; + TimeErrorCalculator(t1s, t1, et1); + + const double pe = std::accumulate(pes.begin(), pes.end(), 0.); + + const CRTBlob blob(t0, et0, t1, et1, pe, taggers); + + candidates.emplace_back(blob, used_spacepoints); + } + + return candidates; +} + +void sbnd::crt::CRTBlobProducer::TimeErrorCalculator(const std::vector ×, double &mean, double &err) +{ + double sum = 0.; + for(auto const &time : times) + sum += time; + + mean = sum / times.size(); + + double summed_var = 0.; + for(auto const &time : times) + summed_var += std::pow((time - mean), 2); + + err = std::sqrt(summed_var / (times.size() - 1)); +} + +void sbnd::crt::CRTBlobProducer::OrderBlobCandidates(std::vector>> &blobCandidates) +{ + std::sort(blobCandidates.begin(), blobCandidates.end(), + [&](const std::pair> &a, const std::pair> &b) -> bool { + return fUseTs0 ? a.first.Ts0Err() < b.first.Ts0Err() : a.first.Ts1Err() < b.first.Ts1Err(); + }); +} + +std::vector>> sbnd::crt::CRTBlobProducer::ChoseBlobs(std::vector>> &blobCandidates) +{ + std::vector>> chosenBlobs; + + std::set used; + + for(auto const& [blob, spIDs] : blobCandidates) + { + bool keep = true; + for(auto const& spID : spIDs) + { + if(used.count(spID) != 0) + keep = false; + } + + if(keep) + { + chosenBlobs.emplace_back(blob, spIDs); + + for(auto const& spID : spIDs) + used.insert(spID); + } + } + + return chosenBlobs; +} + +DEFINE_ART_MODULE(sbnd::crt::CRTBlobProducer) diff --git a/sbndcode/CRT/CRTReco/crtrecoproducers_sbnd.fcl b/sbndcode/CRT/CRTReco/crtrecoproducers_sbnd.fcl index 9aa8e2a26..0b4ca1c3c 100644 --- a/sbndcode/CRT/CRTReco/crtrecoproducers_sbnd.fcl +++ b/sbndcode/CRT/CRTReco/crtrecoproducers_sbnd.fcl @@ -78,4 +78,15 @@ crttrackproducer_data_sbnd: @local::crttrackproducer_sbnd crttrackproducer_data_sbnd.UseTs0: true crttrackproducer_data_sbnd.MaskedTaggers: [ ] +crtblobproducer_sbnd: +{ + CRTSpacePointModuleLabel: "crtspacepoints" + CoincidenceTimeRequirement: 100. + UseTs0: false + module_type: "CRTBlobProducer" +} + +crtblobproducer_data_sbnd: @local::crtblobproducer_sbnd +crtblobproducer_data_sbnd.UseTs0: true + END_PROLOG diff --git a/sbndcode/CRT/CRTReco/run_crtblobreco.fcl b/sbndcode/CRT/CRTReco/run_crtblobreco.fcl new file mode 100644 index 000000000..7d46331d2 --- /dev/null +++ b/sbndcode/CRT/CRTReco/run_crtblobreco.fcl @@ -0,0 +1,38 @@ +#include "services_sbnd.fcl" +#include "rootoutput_sbnd.fcl" +#include "crtrecoproducers_sbnd.fcl" + +process_name: CRTBlobReco + +services: +{ + @table::sbnd_services +} + +source: +{ + module_type: RootInput +} + +outputs: +{ + out1: + { + @table::sbnd_rootoutput + dataTier: "reconstructed" + } +} + +physics: +{ + producers: + { + crtblobs: @local::crtblobproducer_sbnd + } + + reco: [ crtblobs ] + stream1: [ out1 ] + + trigger_paths: [ reco ] + end_paths: [ stream1 ] +} diff --git a/sbndcode/CRT/CRTReco/run_crtreco.fcl b/sbndcode/CRT/CRTReco/run_crtreco.fcl index 3a625c1b7..585238a0a 100644 --- a/sbndcode/CRT/CRTReco/run_crtreco.fcl +++ b/sbndcode/CRT/CRTReco/run_crtreco.fcl @@ -33,9 +33,10 @@ physics: crtclustering: @local::crtclusterproducer_sbnd crtspacepoints: @local::crtspacepointproducer_sbnd crttracks: @local::crttrackproducer_sbnd + crtblobs: @local::crtblobproducer_sbnd } - reco: [ crtstrips, crtclustering, crtspacepoints, crttracks ] + reco: [ crtstrips, crtclustering, crtspacepoints, crttracks, crtblobs ] stream1: [ out1 ] trigger_paths: [ reco ] diff --git a/sbndcode/CRT/CRTReco/run_crtreco_data.fcl b/sbndcode/CRT/CRTReco/run_crtreco_data.fcl index 2493c9b33..9e179b108 100644 --- a/sbndcode/CRT/CRTReco/run_crtreco_data.fcl +++ b/sbndcode/CRT/CRTReco/run_crtreco_data.fcl @@ -10,6 +10,7 @@ physics.producers.crtstrips: @local::crtstriphitproducer_data_sbnd physics.producers.crtclustering: @local::crtclusterproducer_data_sbnd physics.producers.crtspacepoints: @local::crtspacepointproducer_data_sbnd physics.producers.crttracks: @local::crttrackproducer_data_sbnd +physics.producers.crtblobs: @local::crtblobproducer_data_sbnd outputs.out1.outputCommands: [ "keep *_*_*_*", "drop *_daq_*_*", diff --git a/sbndcode/CRT/CRTReco/run_crtrecoana.fcl b/sbndcode/CRT/CRTReco/run_crtrecoana.fcl index 983502ca8..15a37b22d 100644 --- a/sbndcode/CRT/CRTReco/run_crtrecoana.fcl +++ b/sbndcode/CRT/CRTReco/run_crtrecoana.fcl @@ -37,6 +37,7 @@ physics: crtclustering: @local::crtclusterproducer_sbnd crtspacepoints: @local::crtspacepointproducer_sbnd crttracks: @local::crttrackproducer_sbnd + crtblobs: @local::crtblobproducer_sbnd } analyzers: @@ -44,10 +45,13 @@ physics: crtana: @local::crtana_sbnd } - reco: [ crtstrips, crtclustering, crtspacepoints, crttracks ] + reco: [ crtstrips, crtclustering, crtspacepoints, crttracks, crtblobs ] ana: [ crtana ] stream1: [ out1 ] trigger_paths: [ reco ] end_paths: [ ana, stream1 ] } + +physics.analyzers.crtana.HasBlobs: true +physics.analyzers.crtana.NoTPC: true diff --git a/sbndcode/CRT/CRTReco/run_crtrecoana_data.fcl b/sbndcode/CRT/CRTReco/run_crtrecoana_data.fcl index 11bbc53cf..c249a7d99 100644 --- a/sbndcode/CRT/CRTReco/run_crtrecoana_data.fcl +++ b/sbndcode/CRT/CRTReco/run_crtrecoana_data.fcl @@ -10,5 +10,8 @@ physics.producers.crtstrips: @local::crtstriphitproducer_data_sbnd physics.producers.crtclustering: @local::crtclusterproducer_data_sbnd physics.producers.crtspacepoints: @local::crtspacepointproducer_data_sbnd physics.producers.crttracks: @local::crttrackproducer_data_sbnd +physics.producers.crtblobs: @local::crtblobproducer_data_sbnd physics.analyzers.crtana: @local::crtana_data_sbnd +physics.analyzers.crtana.HasBlobs: true +physics.analyzers.crtana.NoTPC: true diff --git a/sbndcode/CRT/CRTReco/run_crtrecoana_notpc.fcl b/sbndcode/CRT/CRTReco/run_crtrecoana_notpc.fcl deleted file mode 100644 index 621b8d45a..000000000 --- a/sbndcode/CRT/CRTReco/run_crtrecoana_notpc.fcl +++ /dev/null @@ -1,3 +0,0 @@ -#include "run_crtrecoana.fcl" - -physics.analyzers.crtana.NoTPC: true diff --git a/sbndcode/CRT/CRTReco/run_crtrecoana_notpc_no_flat_tracks.fcl b/sbndcode/CRT/CRTReco/run_crtrecoana_notpc_no_flat_tracks.fcl deleted file mode 100644 index 74a8504d7..000000000 --- a/sbndcode/CRT/CRTReco/run_crtrecoana_notpc_no_flat_tracks.fcl +++ /dev/null @@ -1,3 +0,0 @@ -#include "run_crtrecoana_notpc.fcl" - -physics.producers.crttracks.MaskedTaggers: [ 0 ]