CRT-TPC Matching with Data

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"
 }
 

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:

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"

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<art::InputTag> TrackLabel {
         Name("TrackLabel")
       };
+      fhicl::Atom<art::InputTag> TPCSpacePointMatchLabel {
+        Name("TPCSpacePointMatchLabel")
+      };
+      fhicl::Atom<art::InputTag> TPCTrackMatchLabel {
+        Name("TPCTrackMatchLabel")
+      };
+      fhicl::Atom<art::InputTag> TPCTrackLabel {
+        Name("TPCTrackLabel")
+      };
 
       fhicl::Atom<bool> SaveRoot {
         Name("SaveRoot"),
@@ -131,6 +142,15 @@ namespace sbnd::crt {
       fhicl::Atom<bool> DrawTracks {
         Name("DrawTracks")
       };
+      fhicl::Atom<bool> DrawTPCMatching {
+        Name("DrawTPCMatching")
+      };
+      fhicl::Atom<bool> OnlyDrawMatched {
+        Name("OnlyDrawMatched")
+      };
+      fhicl::Atom<bool> DisplayMatchScore {
+        Name("DisplayMatchScore")
+      };
 
       fhicl::Atom<bool> ChoseTaggers {
         Name("ChoseTaggers")
@@ -182,6 +202,9 @@ namespace sbnd::crt {
       fhicl::Atom<int> TrackColour {
         Name("TrackColour")
       };
+      fhicl::Atom<int> TPCMatchColour {
+        Name("TPCMatchColour")
+      };
 
       fhicl::Atom<bool> 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<int> 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<int> fChosenTaggers;
@@ -282,6 +314,7 @@ namespace sbnd::crt {
     int fClusterColourInterval;
     int fSpacePointColour;
     int fTrackColour;
+    int fTPCMatchColour;
 
     bool   fUseTs0;
     double fMinTime;

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

sbndcode/CRT/CRTTPCMatching/CMakeLists.txt

@@ -11,13 +11,15 @@ simple_plugin(
               sbnobj::SBND_CRT
               sbndcode_GeoWrappers
               sbndcode_CRT_CRTTPCMatching
+              sbndcode_CRT_CRTEventDisplay
 )
 
 simple_plugin(
               CRTSpacePointMatching module
               sbnobj::SBND_CRT
               sbndcode_GeoWrappers
               sbndcode_CRT_CRTTPCMatching
+              sbndcode_CRT_CRTEventDisplay
 )
 
 install_fhicl()

sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatchAlg.cc

@@ -90,9 +90,9 @@ namespace sbnd::crt {
       std::pair<geo::Vector_t, geo::Vector_t> 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<geo::Vector_t, geo::Vector_t> CRTSpacePointMatchAlg::AverageTrackDirections(const art::Ptr<recob::Track> &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<geo::Vector_t, geo::Vector_t> CRTSpacePointMatchAlg::TrackDirections(detinfo::DetectorPropertiesData const &detProp, const art::Ptr<recob::Track> &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);
-  }
 }

sbndcode/CRT/CRTTPCMatching/CRTSpacePointMatchAlg.h

@@ -167,10 +167,6 @@ namespace sbnd::crt {
                                  const geo::Vector_t &trackDir, const art::Ptr<CRTSpacePoint> &crtSP,
                                  const int driftDirection, const double t0, const art::Event &e);
 
-    std::pair<geo::Vector_t, geo::Vector_t> AverageTrackDirections(const art::Ptr<recob::Track> &track, const double frac);
-    std::pair<geo::Vector_t, geo::Vector_t> TrackDirections(detinfo::DetectorPropertiesData const &detProp, const art::Ptr<recob::Track> &track,
-                                                            const double frac, const double CRTtime, const int driftDirection);
-
   private:
 
     geo::GeometryCore const* fGeometryService;

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 <memory>
 
@@ -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<fhicl::ParameterSet>("MatchingAlg"))
   , fTPCTrackModuleLabel(p.get<art::InputTag>("TPCTrackModuleLabel"))
   , fCRTSpacePointModuleLabel(p.get<art::InputTag>("CRTSpacePointModuleLabel"))
-  {
+  , fCRTEventDisplayAlg(p.get<fhicl::ParameterSet>("CRTEventDisplayAlg"))
+  , fMakeEventDisplays(p.get<bool>("MakeEventDisplays"))
+{
+  if(!fMakeEventDisplays)
     produces<art::Assns<CRTSpacePoint, recob::Track, anab::T0>>();
-  }
+}
 
 void sbnd::crt::CRTSpacePointMatching::produce(art::Event& e)
 {
@@ -104,18 +113,31 @@ void sbnd::crt::CRTSpacePointMatching::produce(art::Event& e)
 
   std::set<int> used_crt_sps;
 
+  int nmatches = 0;
+
   for(auto const &candidate : candidates)
     {
       if(used_crt_sps.count(candidate.thisSP.key()) == 0)
         {
           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<detinfo::DetectorClocksService const>()->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)