Add Hit Reconstruction from Axial and Stereo Tracking Layers to rLDMX Tracking (#1705)

* Modified LinearSeedFinder.cxx and associated header file with code for new 3D hit reconstruction. Currently still in debug stage

* Updated tracking to use the output of the 3D Hit reconstruction algorithm. Also updated uncertainty in y-positions of the 3D recoil hits to reflect the calculation we do to derive them. The performance of the 3D hit reconstruction is currently sub-optimal and will need to be improved.

* Small modification to reconstruction algorithm and general code cleanup

* Noticed error in LinearSeedFinder. The final x position of the reconstructed hit should be u_axial - dx, but we were not projecting properly.

* Modified LinearSeedFinder.cxx and associated header file with code for new 3D hit reconstruction. Currently still in debug stage

* Updated tracking to use the output of the 3D Hit reconstruction algorithm. Also updated uncertainty in y-positions of the 3D recoil hits to reflect the calculation we do to derive them. The performance of the 3D hit reconstruction is currently sub-optimal and will need to be improved.

* Small modification to reconstruction algorithm and general code cleanup

* Noticed error in LinearSeedFinder. The final x position of the reconstructed hit should be u_axial - dx, but we were not projecting properly.

* Apply clang-format

* Fixed GOLD config.py todefine the histograms in LinearSeedFinder.cxx

* Removed LinearSeedFinder histogram declarations

---------

Co-authored-by: github-actions[bot] <github-actions[bot]@users.noreply.github.com>

Author: fdelzanno

.github/validation_samples/reduced_ldmx/config.py

@@ -111,6 +111,7 @@
 rSeedTracking.layer12_midpoint = layer12_mid
 rSeedTracking.layer23_midpoint = layer23_mid
 rSeedTracking.layer34_midpoint = layer34_mid
+rSeedTracking.recoil_uncertainty = [0.006, 0.085]
 rSeedTracking.ecal_distance_threshold = 15.0
 
 rTracking = reducedTracking.LinearTrackFinder("LinearTrackFinder")

EventDisplay/python/rLDMX_EventDisplay.py

@@ -27,9 +27,10 @@ def addBranch(tree, ldmx_class, branch_name):
 
     return branch
 
-def trackPlotter(tree, event_number, tag, save):
+def trackPlotter(tree, event_number, tag, save, save_tag):
 
     recoilSimHits = addBranch(tree, 'SimTrackerHit', 'RecoilSimHits_{}'.format(tag))
+    targetSP = addBranch(tree, 'SimTrackerHit', 'TargetScoringPlaneHits_{}'.format(tag))
     ecalRecHit = addBranch(tree, 'EcalHit', 'EcalRecHits_{}'.format(tag))
     digiRecoil = addBranch(tree, 'Measurement', 'DigiRecoilSimHits_{}'.format(tag))
     recoilTruth = addBranch(tree, 'StraightTrack', 'LinearRecoilTruthTracks_{}'.format(tag))
@@ -61,6 +62,10 @@ def trackPlotter(tree, event_number, tag, save):
     ecal_end_NOISE_y = []
     ecal_end_NOISE_z = []
 
+    targetSP_x = []
+    targetSP_y = []
+    targetSP_z = []
+    
     trackParams = []
     truthTrackParams = []
 
@@ -69,6 +74,12 @@ def trackPlotter(tree, event_number, tag, save):
         recoilSim_x.append(particle.getPosition()[0])
         recoilSim_y.append(particle.getPosition()[1])
 
+    for sp_particle in targetSP:
+        if (sp_particle.getPosition()[2] > 0):
+            targetSP_z.append(sp_particle.getPosition()[2])
+            targetSP_x.append(sp_particle.getPosition()[0])
+            targetSP_y.append(sp_particle.getPosition()[1])
+        
     for x_digi in digiRecoil:
         zpos_digi_tot.append(x_digi.getGlobalPosition()[0])
         xpos_digi_tot.append(x_digi.getGlobalPosition()[1])
@@ -102,14 +113,18 @@ def trackPlotter(tree, event_number, tag, save):
     ecalRecHits = np.column_stack((ecal_end_z, ecal_end_x, ecal_end_y))
     ecalRecHits_noise = np.column_stack((ecal_end_NOISE_z, ecal_end_NOISE_x, ecal_end_NOISE_y))
     recoilSim = np.column_stack((recoilSim_z, recoilSim_x, recoilSim_y))
+    targetSP_hits = np.column_stack((targetSP_z, targetSP_x, targetSP_y))
     first_sensor_pos = np.column_stack((first_sensor_z, first_sensor_x, first_sensor_y))
     second_sensor_pos = np.column_stack((second_sensor_z, second_sensor_x, second_sensor_y))
 
     fig = plt.figure(figsize=(8, 8))
     ax = fig.add_subplot(111, projection='3d')
 
+    #These values come from uncertainty propogation in the function that reconstructs
+    #the 3D sensor point from the axial-stereo combo. See Tracking/LinearSeedFinder.cxx
+    #for the function definitions
     x_sensor_uncertainty = 0.006 #mm
-    y_sensor_uncertainty = 20./np.sqrt(12) # mm
+    y_sensor_uncertainty = 0.085 # mm
 
     for zi, xi, yi in zip(first_sensor_pos[:,0], first_sensor_pos[:,1], first_sensor_pos[:,2]):
         ax.plot([zi, zi], [xi - x_sensor_uncertainty, xi + x_sensor_uncertainty], [yi, yi], color='black',linewidth=2)  # x error
@@ -122,6 +137,7 @@ def trackPlotter(tree, event_number, tag, save):
     ax.scatter(digiPoints[:,0], digiPoints[:,1], digiPoints[:,2], c='b', marker='o', label='DigiRecoil SimHit', s=25)
     ax.scatter(ecalRecHits[:, 0], ecalRecHits[:, 1], ecalRecHits[:, 2], c='purple', label='ECalRecHit', s=50, alpha=0.5)
     ax.scatter(recoilSim[:, 0], recoilSim[:, 1], recoilSim[:, 2], marker='o', c='gray', label='RecoilSimHits', s=100, alpha=0.25)
+    ax.scatter(targetSP_hits[:, 0], targetSP_hits[:, 1], targetSP_hits[:, 2], marker='o', c='yellow', label='Target SP Hits', s=50, alpha=0.5)
     ax.scatter(first_sensor_pos[:,0], first_sensor_pos[:,1], first_sensor_pos[:,2], marker='o', c='red', label='First Sensor Point', s=30, alpha=0.75)
     ax.scatter(second_sensor_pos[:,0], second_sensor_pos[:,1], second_sensor_pos[:,2], marker='o', c='orange', label='Second Sensor Point', s=30, alpha=0.75)
 
@@ -143,9 +159,9 @@ def trackPlotter(tree, event_number, tag, save):
     z_range = np.linspace(min(ecalRecHits[:,2]) - 5, max(ecalRecHits[:,2]) + 5, 50)
     Y, Z = np.meshgrid(y_range, z_range)
     X = np.full(Y.shape, ecalRecHits[0][0])
-
+        
     plane = ax.plot_surface(X, Y, Z, color='red', alpha=0.1, edgecolor='none')
-    
+        
     for recHit in ecalRecHits:
         center_z, center_x, center_y = recHit[0], recHit[1], recHit[2]
         radius = 3.87
@@ -169,16 +185,16 @@ def trackPlotter(tree, event_number, tag, save):
 
     ax.legend(fontsize=10)
 
-    plt.savefig(f'{save}/{tag}_eventDisplay_eventID_{event_number}.png')
+    plt.savefig(f'{save}/{tag}_eventDisplay_eventID_{event_number}_{save_tag}.png')
 
 #    UNCOMMENT THESE LINES IF YOU WANT A ZOOMED IN PLOT (i.e. just recoil points)
 #    ax.set_xlim(0, 50)
 #
-#    x_low_lim = -10
-#    x_up_lim = -5
-#    y_low_lim = 14
-#    y_up_lim = 17
-#    
+#    x_low_lim = -6
+#    x_up_lim = -4
+#    y_low_lim = -18.5
+#    y_up_lim = -17
+#
 #    x_target = 0.0
 #    y_target_range = np.linspace(x_low_lim, x_up_lim, 50)
 #    z_target_range = np.linspace(y_low_lim, y_up_lim, 50)
@@ -194,16 +210,19 @@ def trackPlotter(tree, event_number, tag, save):
 
 def main():
     tree = r.TChain("LDMX_Events")
-    tree.Add('events_5000_rLDMXV1.root')
+    tree.Add('events_5000_rLDMX_v1_yZERO_new3Dpoints.root')
 
     nentries = tree.GetEntries()
     print("nentries = ", nentries)
 
     tag = 'rLDMX_v1'
-    save_loc = 'event_displays/rLDMX_v1'
-    event_number = 10
+    save_loc = 'plots/event_displays/'
+    event_number = 26
+    save_tag = 'yZERO_new3Dpoints'
 
-    trackPlotter(tree, event_number, tag, save_loc)
+#    random_numbers = random.sample(range(0, 101), 10)
+#    for number in random_numbers:
+    trackPlotter(tree, event_number, tag, save_loc, save_tag)
 
 if __name__ == "__main__":
     main()

Tracking/include/Tracking/Event/Measurement.h

@@ -111,7 +111,7 @@ class Measurement {
   /// Add a trackId to the internal vector
   void addTrackId(int trkId) { trackIds_.push_back(trkId); };
   /// @return the sim particle IDs that compose the measurement
-  std::vector<unsigned int> getTrackIds() { return trackIds_; };
+  std::vector<unsigned int> getTrackIds() const { return trackIds_; };
 
   ///  @return The energy deposited in the sensor where the measurement took
   ///  place.

Tracking/include/Tracking/Reco/LinearSeedFinder.h

@@ -69,19 +69,6 @@ class LinearSeedFinder : public TrackingGeometryUser {
           recoil_two,
       const std::array<double, 3> ecal_one);
 
-  // Function to combine Recoil layer points into "real" sensor points
-  std::pair<std::vector<std::tuple<std::array<double, 3>, ldmx::Measurement,
-                                   std::optional<ldmx::Measurement>>>,
-            std::vector<std::tuple<std::array<double, 3>, ldmx::Measurement,
-                                   std::optional<ldmx::Measurement>>>>
-  combineMultiGlobalHits(const std::vector<ldmx::Measurement> &hit_collection);
-
-  // Function to do weighted averaging when combining Recoil layer points
-  std::vector<std::tuple<std::array<double, 3>, ldmx::Measurement,
-                         std::optional<ldmx::Measurement>>>
-  midPointCalculation(const std::vector<ldmx::Measurement> &layer1,
-                      const std::vector<ldmx::Measurement> &layer2);
-
   // Fitting function: fit a straight line in 3D using 3 points (1 degree of
   // freedom)
   std::tuple<double, double, double, double, std::vector<double>> fit3DLine(
@@ -93,6 +80,26 @@ class LinearSeedFinder : public TrackingGeometryUser {
   double calculateDistance(const std::array<double, 3> &point1,
                            const std::array<double, 3> &point2);
 
+  // Do 3D hit reconstruction using an axial measurement and a stereo
+  // measurement according to geometric projections
+  Acts::Vector3 simple3DHitV2(const ldmx::Measurement &axial,
+                              const Acts::Surface &axial_surface,
+                              const ldmx::Measurement &stereo,
+                              const Acts::Surface &stereo_surface,
+                              const ldmx::SimTrackerHit &hitOnTarget,
+                              std::vector<ldmx::SimTrackerHit> pair_sim_hits);
+
+  // Makes all combinations of sensor measurements to use in the seeding
+  std::vector<std::tuple<
+      std::array<double, 3>,
+      std::tuple<ldmx::Measurement, ldmx::SimTrackerHit, ldmx::SimTrackerHit>,
+      std::optional<std::tuple<ldmx::Measurement, ldmx::SimTrackerHit,
+                               ldmx::SimTrackerHit>>>>
+  processMeasurements(
+      const std::vector<std::tuple<ldmx::Measurement, ldmx::SimTrackerHit,
+                                   ldmx::SimTrackerHit>> &measurements,
+      const geo::TrackersTrackingGeometry &tg);
+
   // Calculate chi2 of the fit
   double globalChiSquare(const std::array<double, 3> &first_sensor,
                          const std::array<double, 3> &second_sensor,
@@ -103,6 +110,15 @@ class LinearSeedFinder : public TrackingGeometryUser {
   // enough points to fit)
   int uniqueLayersHit(const std::vector<ldmx::Measurement> &digi_points);
 
+  // Helper functions:
+  std::array<double, 3> convertToLdmxStdArray(const Acts::Vector3 &vec);
+  std::tuple<Acts::Vector3, Acts::Vector3, Acts::Vector3> getSurfaceVectors(
+      const Acts::Surface &surface);
+  double dotProduct(const Acts::Vector3 &v1, const Acts::Vector3 &v2);
+  std::array<double, 3> getPointAtZ(std::array<double, 3> target,
+                                    std::array<double, 3> measurement,
+                                    double z_target);
+
   double processing_time_{0.};
   long n_events_{0};
   unsigned int n_seeds_{0};
@@ -125,7 +141,7 @@ class LinearSeedFinder : public TrackingGeometryUser {
   double layer34_midpoint_{27.5};
   double ecal_first_layer_z_threshold_{250.0};
 
-  std::vector<double> recoil_uncertainty_{0.006, 5.7735};
+  std::vector<double> recoil_uncertainty_{0.006, 0.085};
 
   // Check failures
   long n_missing_{0};

Tracking/python/reducedTracking.py

@@ -37,7 +37,7 @@ def __init__(self, instance_name="LinearSeedFinder"):
         self.input_hits_collection = 'DigiRecoilSimHits'
         self.input_rec_hits_collection = 'EcalRecHits'
         self.out_seed_collection = 'LinearRecoilSeedTracks'
-        self.recoil_uncertainty = [0.006, 20./(12**(0.5))]
+        self.recoil_uncertainty = [0.006, 0.085]
         self.ecal_uncertainty = 3.87
         self.ecal_first_layer_z_threshold = 250.0
         self.ecal_distance_threshold = 15.0