CRT Blobs

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<art::Ptr<CRTTrack>> &CRTTrackVec, const art::FindManyP<CRTSpacePoint> &tracksToSpacePoints,
                         const art::FindOneP<CRTCluster> &spacePointsToClusters, const art::FindManyP<CRTStripHit> &clustersToStripHits);
 
+  void AnalyseCRTBlobs(std::vector<art::Ptr<CRTBlob>> &CRTBlobVec);
+
   void AnalyseTPCMatching(const art::Event &e, const art::Handle<std::vector<recob::Track>> &TPCTrackHandle,
                           const art::Handle<std::vector<CRTSpacePoint>> &CRTSpacePointHandle, const art::Handle<std::vector<CRTCluster>> &CRTClusterHandle,
                           const art::Handle<std::vector<recob::PFParticle>> &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<double>                _tr_truth_theta;
   std::vector<double>                _tr_truth_phi;
 
+  // crt blob information
+  std::vector<double>           _bl_ts0;
+  std::vector<double>           _bl_ets0;
+  std::vector<double>           _bl_ts1;
+  std::vector<double>           _bl_ets1;
+  std::vector<double>           _bl_pe;
+  std::vector<int>              _bl_nsps;
+  std::vector<std::vector<int>> _bl_nsps_per_tagger;
+
+  // tpc track information (including crt matching)
   std::vector<double>                _tpc_start_x;
   std::vector<double>                _tpc_start_y;
   std::vector<double>                _tpc_start_z;
@@ -341,12 +354,14 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   std::vector<double>                _tpc_tr_end_z;
   std::vector<double>                _tpc_tr_score;
 
+  // ptb information (trigger board)
   std::vector<uint64_t> _ptb_hlt_trigger;
   std::vector<uint64_t> _ptb_hlt_timestamp;
 
   std::vector<uint64_t> _ptb_llt_trigger;
   std::vector<uint64_t> _ptb_llt_timestamp;
 
+  // spec tdc information (timing board)
   std::vector<uint32_t>    _tdc_channel;
   std::vector<uint64_t>    _tdc_timestamp;
   std::vector<uint64_t>    _tdc_offset;
@@ -364,6 +379,7 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
   fCRTClusterModuleLabel            = p.get<std::string>("CRTClusterModuleLabel", "crtclustering");
   fCRTSpacePointModuleLabel         = p.get<std::string>("CRTSpacePointModuleLabel", "crtspacepoints");
   fCRTTrackModuleLabel              = p.get<std::string>("CRTTrackModuleLabel", "crttracks");
+  fCRTBlobModuleLabel               = p.get<std::string>("CRTBlobModuleLabel", "crtblobs");
   fTPCTrackModuleLabel              = p.get<std::string>("TPCTrackModuleLabel", "pandoraSCETrack");
   fCRTSpacePointMatchingModuleLabel = p.get<std::string>("CRTSpacePointMatchingModuleLabel", "crtspacepointmatchingSCE");
   fCRTTrackMatchingModuleLabel      = p.get<std::string>("CRTTrackMatchingModuleLabel", "crttrackmatchingSCE");
@@ -376,8 +392,8 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
   fNoTPC                            = p.get<bool>("NoTPC", false);
   fHasPTB                           = p.get<bool>("HasPTB", false);
   fHasTDC                           = p.get<bool>("HasTDC", false);
+  fHasBlobs                         = p.get<bool>("HasBlobs", false);
   fTruthMatch                       = p.get<bool>("TruthMatch", true);
-  //! Adding some of the reco parameters to save corrections
   fPEAttenuation                    = p.get<double>("PEAttenuation", 1.0);
   fTimeWalkNorm                     = p.get<double>("TimeWalkNorm",  0.0);
   fTimeWalkScale                    = p.get<double>("TimeWalkScale", 0.0);
@@ -578,6 +594,17 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
       fTree->Branch("tr_truth_phi", "std::vector<double>", &_tr_truth_phi);
     }
 
+  if(fHasBlobs)
+    {
+      fTree->Branch("bl_ts0", "std::vector<double>", &_bl_ts0);
+      fTree->Branch("bl_ets0", "std::vector<double>", &_bl_ets0);
+      fTree->Branch("bl_ts1", "std::vector<double>", &_bl_ts1);
+      fTree->Branch("bl_ets1", "std::vector<double>", &_bl_ets1);
+      fTree->Branch("bl_pe", "std::vector<double>", &_bl_pe);
+      fTree->Branch("bl_nsps", "std::vector<int>", &_bl_nsps);
+      fTree->Branch("bl_nsps_per_tagger", "std::vector<std::vector<int>>", &_bl_nsps_per_tagger);
+    }
+
   if(!fNoTPC)
     {
       fTree->Branch("tpc_start_x", "std::vector<double>", &_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<std::vector<CRTBlob>> CRTBlobHandle;
+      e.getByLabel(fCRTBlobModuleLabel, CRTBlobHandle);
+      if(!CRTBlobHandle.isValid()){
+        std::cout << "CRTBlob product " << fCRTBlobModuleLabel << " not found..." << std::endl;
+        throw std::exception();
+      }
+
+      std::vector<art::Ptr<CRTBlob>> 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<art::Ptr<CRTBlob>> &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<int>(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<std::vector<recob::Track>> &TPCTrackHandle,
                                                 const art::Handle<std::vector<CRTSpacePoint>> &CRTSpacePointHandle, const art::Handle<std::vector<CRTCluster>> &CRTClusterHandle,
                                                 const art::Handle<std::vector<recob::PFParticle>> &PFPHandle,

sbndcode/CRT/CRTReco/CMakeLists.txt

@@ -37,4 +37,10 @@ simple_plugin(
               Eigen3::Eigen
 )
 
+simple_plugin(
+              CRTBlobProducer module
+              lardata::Utilities
+              sbnobj::SBND_CRT
+)
+
 install_fhicl()

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 <numeric>
+
+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<art::Ptr<CRTSpacePoint>> &spacePointVec, const art::FindOneP<CRTCluster> &spacePointsToCluster);
+
+  std::vector<std::pair<CRTBlob, std::set<unsigned>>> CreateBlobCandidates(const std::vector<art::Ptr<CRTSpacePoint>> &spacePointVec,
+                                                                           const art::FindOneP<CRTCluster> &spacePointsToCluster);
+
+  void TimeErrorCalculator(const std::vector<double> &times, double &mean, double &err);
+
+  void OrderBlobCandidates(std::vector<std::pair<CRTBlob, std::set<unsigned>>> &blobCandidates);
+
+  std::vector<std::pair<sbnd::crt::CRTBlob, std::set<unsigned>>> ChoseBlobs(std::vector<std::pair<CRTBlob, std::set<unsigned>>> &blobCandidates);
+
+private:
+
+  std::string      fCRTSpacePointModuleLabel;
+  double           fCoincidenceTimeRequirement;
+  bool             fUseTs0;
+};
+
+
+sbnd::crt::CRTBlobProducer::CRTBlobProducer(fhicl::ParameterSet const& p)
+  : EDProducer{p}
+  , fCRTSpacePointModuleLabel(p.get<std::string>("CRTSpacePointModuleLabel"))
+  , fCoincidenceTimeRequirement(p.get<double>("CoincidenceTimeRequirement"))
+  , fUseTs0(p.get<bool>("UseTs0"))
+{
+  produces<std::vector<CRTBlob>>();
+  produces<art::Assns<CRTSpacePoint, CRTBlob>>();
+}
+
+void sbnd::crt::CRTBlobProducer::produce(art::Event& e)
+{
+  auto blobVec            = std::make_unique<std::vector<CRTBlob>>();
+  auto blobSpacePointAssn = std::make_unique<art::Assns<CRTSpacePoint, CRTBlob>>();
+
+  art::Handle<std::vector<CRTSpacePoint>> CRTSpacePointHandle;
+  e.getByLabel(fCRTSpacePointModuleLabel, CRTSpacePointHandle);
+
+  std::vector<art::Ptr<CRTSpacePoint>> CRTSpacePointVec;
+  art::fill_ptr_vector(CRTSpacePointVec, CRTSpacePointHandle);
+
+  art::FindOneP<CRTCluster> spacePointsToCluster(CRTSpacePointHandle, e, fCRTSpacePointModuleLabel);
+
+  OrderSpacePoints(CRTSpacePointVec, spacePointsToCluster);
+
+  std::vector<std::pair<CRTBlob, std::set<unsigned>>> blobCandidates = CreateBlobCandidates(CRTSpacePointVec, spacePointsToCluster);
+
+  OrderBlobCandidates(blobCandidates);
+
+  std::vector<std::pair<CRTBlob, std::set<unsigned>>> 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<art::Ptr<CRTSpacePoint>> &spacePointVec, const art::FindOneP<CRTCluster> &spacePointsToCluster)
+{
+  std::sort(spacePointVec.begin(), spacePointVec.end(), 
+            [&](const art::Ptr<CRTSpacePoint> &a, const art::Ptr<CRTSpacePoint> &b) -> bool {
+              if(fUseTs0)
+                return a->Ts0() < b->Ts0();
+              else
+                return a->Ts1() < b->Ts1();
+            });
+}
+
+std::vector<std::pair<sbnd::crt::CRTBlob, std::set<unsigned>>> sbnd::crt::CRTBlobProducer::CreateBlobCandidates(const std::vector<art::Ptr<CRTSpacePoint>> &spacePointVec,
+                                                                                                                const art::FindOneP<CRTCluster> &spacePointsToCluster)
+{
+  std::vector<std::pair<CRTBlob, std::set<unsigned>>> candidates;
+
+  for(unsigned i = 0; i < spacePointVec.size(); ++i)
+    {
+      const art::Ptr<CRTSpacePoint> primarySpacePoint = spacePointVec[i];
+      const art::Ptr<CRTCluster> primaryCluster       = spacePointsToCluster.at(primarySpacePoint.key());
+
+      std::set<unsigned> used_spacepoints = { i };
+      std::vector<double> t0s             = { primarySpacePoint->Ts0() };
+      std::vector<double> t1s             = { primarySpacePoint->Ts1() };
+      std::vector<double> pes             = { primarySpacePoint->PE() };
+
+      std::map<CRTTagger, int> 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<CRTSpacePoint> secondarySpacePoint = spacePointVec[ii];
+          const art::Ptr<CRTCluster> 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> &times, 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<std::pair<CRTBlob, std::set<unsigned>>> &blobCandidates)
+{
+  std::sort(blobCandidates.begin(), blobCandidates.end(),
+            [&](const std::pair<CRTBlob, std::set<unsigned>> &a, const std::pair<CRTBlob, std::set<unsigned>> &b) -> bool {
+              return fUseTs0 ? a.first.Ts0Err() < b.first.Ts0Err() : a.first.Ts1Err() < b.first.Ts1Err();
+            });
+}
+
+std::vector<std::pair<sbnd::crt::CRTBlob, std::set<unsigned>>> sbnd::crt::CRTBlobProducer::ChoseBlobs(std::vector<std::pair<CRTBlob, std::set<unsigned>>> &blobCandidates)
+{
+  std::vector<std::pair<sbnd::crt::CRTBlob, std::set<unsigned>>> chosenBlobs;
+
+  std::set<unsigned> 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)