Merge pull request #44 from DUNE/kleykamp_areal_density_fix

Fix areal density by using geometry to tag y and v bars, and fixing CalculateTrackLength

src/TMS_Bar.cpp

@@ -45,6 +45,15 @@ bool TMS_Bar::FindModules(double xval, double yval, double zval) {
     // We've found the plane number
     if (NodeName.find(TMS_Const::TMS_ModuleLayerName) != std::string::npos) {
       PlaneNumber = geom->GetCurrentNode()->GetNumber();
+      // There are two rotations of bars, and their names are literally "modulelayervol1" and "modulelayervol2"
+      if (NodeName.find(TMS_Const::TMS_ModuleLayerName1) != std::string::npos) BarOrient = kYBar;
+      else if (NodeName.find(TMS_Const::TMS_ModuleLayerName2) != std::string::npos) BarOrient = kVBar;
+      else if (NodeName.find(TMS_Const::TMS_ModuleLayerName3) != std::string::npos) BarOrient = kXBar;
+      else if (NodeName.find(TMS_Const::TMS_ModuleLayerName4) != std::string::npos) BarOrient = kUBar;
+      else {
+        std::cout<<"Error: Do not understand TMS_ModuleLayerName '"<<NodeName<<"'. This bar will have type kError."<<std::endl;
+        BarOrient = kError;
+      }
     }
 
     // This is the furthest down hit we have: scintillator bar
@@ -101,16 +110,15 @@ bool TMS_Bar::FindModules(double xval, double yval, double zval) {
   if (PlaneNumber % 2 == 0) BarOrient = kXBar;
   else BarOrient = kYBar;
   */
-  BarOrient = kYBar;
 
   // If this is a y-bar, remove the y coordinate
-  if (BarOrient == kXBar) {
+  if (BarOrient == kXBar || BarOrient == kUBar) {
     x = -99999000;
     // Flip the widths
     double tempyw = yw;
     yw = xw;
     xw = tempyw; 
-  } else if (BarOrient == kYBar) {
+  } else if (BarOrient == kYBar || BarOrient == kVBar) {
     y = -99999000;
     // Don't need to flip the widths because they're already correct (yw = large, xw = 4cm)
   } else {

src/TMS_Constants.h

@@ -90,6 +90,10 @@ namespace TMS_Const {
   const std::string TMS_EDepSim_VolumeName = "volTMS";
   // To find in z
   const std::string TMS_ModuleLayerName = "modulelayervol";
+  const std::string TMS_ModuleLayerName1 = "modulelayervol1"; // y orientation
+  const std::string TMS_ModuleLayerName2 = "modulelayervol2"; // v orientation
+  const std::string TMS_ModuleLayerName3 = "modulelayervol3"; // x orientation
+  const std::string TMS_ModuleLayerName4 = "modulelayervol4"; // u orientation
   // To find scintillator "box"
   const std::string TMS_ModuleName = "ModuleBoxvol";
   // To find scintillator "box"

src/TMS_Geom.h

@@ -304,12 +304,6 @@ class TMS_Geom {
         std::cout << "  thickness*density = " << density*thickness << std::endl;
         */
 
-        // Skip if density or thickness is small
-        if (density*thickness < 0.1) {
-          //std::cout << "  Skipping material, to little path length to bother" << std::endl;
-          continue;
-        }
-
         TotalPathLength += density*thickness;
         TotalLength += thickness;
 

src/TMS_Reco.cpp

@@ -523,24 +523,39 @@ void TMS_TrackFinder::CalculateTrackEnergy() {
 void TMS_TrackFinder::CalculateTrackLength() {
   // Look at the reconstructed tracks
   if (HoughCandidates.size() == 0) return;
-
   // Loop over each Hough Candidate and find the track length
   for (auto it = HoughCandidates.begin(); it != HoughCandidates.end(); ++it) {
-    double total = 0;
-
-    // Sort by increasing z
-    std::sort((*it).begin(), (*it).end(), TMS_Hit::SortByZInc);
-
-    for (auto hit = (*it).begin(); hit != (*it).end(); ++hit) {
-      auto nexthit = *(hit+1);
-      // Use the geometry to calculate the track length between hits
-      TVector3 point1((*hit).GetNotZ(), -200, (*hit).GetZ());
-      TVector3 point2(nexthit.GetNotZ(), -200, nexthit.GetZ());
-      if ((point2-point1).Mag() > 100) continue;
-      double tracklength = TMS_Geom::GetInstance().GetTrackLength(point1, point2);
-      total += tracklength;
+    double final_total = 0;
+    int max_n_nodes_used = 0;
+
+    // Loop through all bar types so we're only calculating along the same plane rotation type
+    // Otherwise it would overestimate the track length by zig-zagging between y and v
+    TMS_Bar::BarType bartypes[] = {TMS_Bar::kXBar, TMS_Bar::kYBar, TMS_Bar::kUBar, TMS_Bar::kVBar, TMS_Bar::kError};
+    for (auto itbartype = std::begin(bartypes); itbartype != std::end(bartypes); ++itbartype) {
+      // Get only the nodes with the current bar type
+      auto track_hits_only_matching_bar = ProjectHits((*it), (*itbartype));
+      double total = 0;
+      int n_nodes = 0;
+      for (auto hit = track_hits_only_matching_bar.begin(); hit != track_hits_only_matching_bar.end() \
+          && (hit+1) != track_hits_only_matching_bar.end(); ++hit) {
+        auto nexthit = *(hit+1);
+        // Use the geometry to calculate the track length between hits
+        TVector3 point1((*hit).GetNotZ(), -200, (*hit).GetZ());
+        TVector3 point2(nexthit.GetNotZ(), -200, nexthit.GetZ());
+        double tracklength = TMS_Geom::GetInstance().GetTrackLength(point1, point2);
+        total += tracklength;
+        n_nodes += 1;
+      }
+      // Currently tracks are made with ProjectHits so tracks will only have one (y) bar exclusively
+      // So taking the only nonzero node is a good proxy
+      // In the future we might want to take an average or take the max length, etc, or really do it in 3d
+      // But this captures possible future cases where tracks might include both y and v views
+      if (n_nodes > max_n_nodes_used) {
+        final_total = total;
+        max_n_nodes_used = n_nodes;
+      }
     }
-    TrackLength.push_back(total);
+    TrackLength.push_back(final_total);
   }
 }
 
@@ -873,7 +888,7 @@ std::vector<std::vector<TMS_Hit> > TMS_TrackFinder::FindClusters(const std::vect
 std::vector<TMS_Hit> TMS_TrackFinder::RunHough(const std::vector<TMS_Hit> &TMS_Hits) {
 
   // Check if we're in XZ view
-  bool IsXZ = ((TMS_Hits[0].GetBar()).GetBarType() == TMS_Bar::kYBar);
+  bool IsXZ = ((TMS_Hits[0].GetBar()).GetBarType() == TMS_Bar::kYBar || (TMS_Hits[0].GetBar()).GetBarType() == TMS_Bar::kVBar);
 
   // Recalculate Hough parameters event by event... not fully tested!
   bool VariableHough = false;
@@ -1395,7 +1410,7 @@ std::vector<TMS_Hit> TMS_TrackFinder::RunAstar(const std::vector<TMS_Hit> &TMS_x
 
   // Remember which orientation these hits are
   // needed when we potentially skip the air gap in xz (but not in yz!)
-  bool IsXZ = ((TMS_xz[0].GetBar()).GetBarType() == TMS_Bar::kYBar);
+  bool IsXZ = ((TMS_xz[0].GetBar()).GetBarType() == TMS_Bar::kYBar || (TMS_xz[0].GetBar()).GetBarType() == TMS_Bar::kVBar);
   // Reset remembering where gaps are in xz
   if (IsXZ) PlanesNearGap.clear();