Merge pull request #51 from DUNE/50-missing-hits-and-starts-and-ends-…

…of-tracks

[Config]

    Changed MergeTracks default to True.

[General]

    Hits trailing the end of a Hough Track can now be added to the track they 'should' correspond to. Currently this defaults to on, will be switchable later.

config/TMS_Default_Config.toml

@@ -44,10 +44,12 @@
     # Do we run DBSCAN clustering before the Hough transform and then run Hough on each cluster?
     FirstCluster = false
     # Do we merge adjacent tracks afterwards?
-    MergeTracks = false
+    MergeTracks = true
     # Do we run AStar afterwards on start and end points
     # To undo the greediness in the Hough line finding?
     RunAStarCleanup = true
+    ExtrapolateDist = 4
+    ExtrapolateLimit = 20
 
   [Recon.AStar]
     #CostMetric = "Manhattan"
@@ -82,4 +84,4 @@
 # Draw PDF of "event display". Slows down reco considerably
 [Applications]
   DrawPDF = false
-  MaximumNEvents = -1
+  MaximumNEvents = 100

scripts/Reco/draw_spill.py

@@ -10,6 +10,7 @@
 ROOT.TH1.AddDirectory(False)
 #ROOT.TImage.AddDirectory(False)
 ROOT.TH2.AddDirectory(False)
+import numpy
 
 
 ROOT.gStyle.SetPalette(52) # greyscale
@@ -174,6 +175,16 @@ def inside_tms(x, y, z):
                 elif reco_hit_energy < 0:
                     print(f"Found unexpected reco_hit_energy < 0: {reco_hit_energy}")
 
+                # Extract all hits in clusters and tracks that are not -999 for easier access later on
+                TotalHitsInClusters = sum(numpy.array([event.nHitsInCluster[i] for i in range(min(25, event.nClusters))]))
+                FilteredClusterHitPos = numpy.empty(TotalHitsInClusters*2)
+                j = 0
+                for i in range(len(event.ClusterHitPos)):
+                  if event.ClusterHitPos[i] != -999:
+                    FilteredClusterHitPos[j] = event.ClusterHitPos[i]
+                    j += 1
+
+                usedClusters = 0
                 for cluster in range(min(25, event.nClusters)):
                     cluster_energy = event.ClusterEnergy[cluster]
                     cluster_time = event.ClusterTime[cluster]
@@ -185,7 +196,7 @@ def inside_tms(x, y, z):
                     cluster_total_std_dev = max(cluster_pos_z_std_dev, cluster_pos_x_std_dev)
                     e = int(min(255, 255 * cluster_energy / 10.0))
                     t = int(min(255, 255 * cluster_time / 10000.0))
-                    color = ROOT.kBlack # ROOT.TColor.GetColor(e, t, 128)
+                    color = ROOT.kAzure #kBlack # ROOT.TColor.GetColor(e, t, 128)
                     x = cluster_pos_z
                     y = cluster_pos_x
                     #print("cluster:", cluster, e, t, 255 * cluster_energy / 10.0, 255 * cluster_time / 10000.0, color, x, y, cluster_total_std_dev)
@@ -194,7 +205,28 @@ def inside_tms(x, y, z):
                     marker.SetMarkerColor(color)
                     marker.SetMarkerSize(20 * cluster_total_std_dev)
                     markers.append(marker)
-
+                    marker = 0
+
+                    for ClusterHit in range(event.nHitsInCluster[cluster]):
+                      hit_z = FilteredClusterHitPos[usedClusters*2 + ClusterHit*2 + 0] / 1000.0
+                      hit_x = FilteredClusterHitPos[usedClusters*2 + ClusterHit*2 + 1] / 1000.0
+
+                      marker = ROOT.TMarker(hit_z, hit_x, 21)
+                      marker.SetMarkerColor(ROOT.kAzure-8)
+                      marker.SetMarkerSize(0.5)
+                      markers.append(marker)
+                      marker = 0
+                    usedClusters += event.nHitsInCluster[cluster]
+
+                TotalHitsInTracks = sum(numpy.array([event.nHitsInTrack[i] for i in range(event.nLines)]))
+                FilteredTrackHitPos = numpy.empty(TotalHitsInTracks*2)
+                j = 0
+                for i in range(len(event.TrackHitPos)):
+                  if event.TrackHitPos[i] != -999:
+                    FilteredTrackHitPos[j] = event.TrackHitPos[i]
+                    j += 1
+
+                usedTracks = 0    
                 for line in range(event.nLines):
                     #track_z = event.TrackHitPos[line*2 + 0] / 1000.0
                     #track_x = event.TrackHitPos[line*2 + 1] / 1000.0
@@ -223,6 +255,7 @@ def inside_tms(x, y, z):
                         lines.SetLineWidth(4)
                         lines.SetLineColorAlpha(ROOT.kBlue, 0.15)
                         markers.append(lines)
+                        lines = 0
 
                     # Add markers for front and end of track
                     offset = 0.075
@@ -234,6 +267,17 @@ def inside_tms(x, y, z):
                     marker_end.SetMarkerColor(ROOT.kRed)
                     marker_end.SetMarkerSize(2)
                     markers.append(marker_end)
+
+                    for TrackHit in range(event.nHitsInTrack[line]):
+                      hit_z = FilteredTrackHitPos[usedTracks*2 + TrackHit*2 + 0] / 1000.0
+                      hit_x = FilteredTrackHitPos[usedTracks*2 + TrackHit*2 + 1] / 1000.0
+
+                      marker = ROOT.TMarker(hit_z, hit_x, 21)
+                      marker.SetMarkerColor(ROOT.kPink-3)
+                      marker.SetMarkerSize(0.5)
+                      markers.append(marker)
+                      marker = 0
+                    usedTracks += event.nHitsInTrack[line]
 
 
         minimum_energy_to_print = 0

src/TMS_Manager.cpp

@@ -41,6 +41,8 @@ TMS_Manager::TMS_Manager() {
   _RECO_HOUGH_HitMult = toml::find<double>(data, "Recon", "Hough", "HitMult");
   _RECO_HOUGH_MergeTracks = toml::find<bool>(data, "Recon", "Hough", "MergeTracks");
   _RECO_HOUGH_RunAStar = toml::find<bool>(data, "Recon", "Hough", "RunAStarCleanup");
+  _RECO_HOUGH_ExtrapolateDist = toml::find<int>(data, "Recon", "Hough", "ExtrapolateDist");
+  _RECO_HOUGH_ExtrapolateLimit = toml::find<int>(data, "Recon", "Hough", "ExtrapolateLimit");
   _RECO_HOUGH_FirstCluster = toml::find<bool>(data, "Recon", "Hough", "FirstCluster");
   _RECO_HOUGH_MinDist = toml::find<double>(data, "Recon", "Hough", "MinDist");
 

src/TMS_Manager.h

@@ -35,6 +35,8 @@ class TMS_Manager {
     double Get_Reco_HOUGH_HitMult() { return _RECO_HOUGH_HitMult; };
     double Get_Reco_HOUGH_MergeTracks() { return _RECO_HOUGH_MergeTracks; };
     bool Get_Reco_HOUGH_RunAStar() { return _RECO_HOUGH_RunAStar; };
+    int Get_Reco_HOUGH_ExtrapolateDist() { return _RECO_HOUGH_ExtrapolateDist; };
+    int Get_Reco_HOUGH_ExtrapolateLimit() { return _RECO_HOUGH_ExtrapolateLimit; };
     bool Get_Reco_HOUGH_FirstCluster() { return _RECO_HOUGH_FirstCluster; };
     double Get_Reco_HOUGH_MinDist() { return _RECO_HOUGH_MinDist; };
 
@@ -86,6 +88,8 @@ class TMS_Manager {
     double _RECO_HOUGH_HitMult;
     bool _RECO_HOUGH_MergeTracks;
     bool _RECO_HOUGH_RunAStar;
+    int _RECO_HOUGH_ExtrapolateDist;
+    int _RECO_HOUGH_ExtrapolateLimit;
     bool _RECO_HOUGH_FirstCluster;
     double _RECO_HOUGH_MinDist;