Implementing NN into Strip Pairing

include/MHit.h

@@ -52,9 +52,13 @@ class MHit
   //! Return strip hit i or nullptr if i is out of bounds
   //! Ownership stays elsewhere
   MStripHit* GetStripHit(unsigned int i);
+  //! Return nearest neighbor strip hit i or nullptr if i is out of bounds
+  MStripHit* GetNearestNeighborStripHit(unsigned int i);
   //! Add a strip hit
   //! Ownership stays elsewhere
   void AddStripHit(MStripHit* StripHit);
+  //! Add a nearest neighbor strip hit
+  void AddNearestNeighborStripHit(MStripHit* StripHit);
   //! Remove strip hit i without deleting it
   void RemoveStripHit(unsigned int i);
   //! Remove a strip hit without deleting it
@@ -172,6 +176,9 @@ class MHit
   //! List of strip hits contributing to this hit
   //! Ownership stays elsewhere
   vector<MStripHit*> m_StripHits;
+
+  //! List of nearest neighbor strip hits associated with this hit
+  vector<MStripHit*> m_NearestNeighborStripHits;
 
   //! Position of the hit
   MVector m_Position;

include/MModuleStripPairingMultiRoundChiSquare.h

@@ -87,7 +87,13 @@ class MModuleStripPairingMultiRoundChiSquare : public MModule
     float ChargeTrappingCorrection(unsigned int d, const vector<vector<MStripHit*>>& StripHits);
 
   //! Divide an event's strip hits by detector and LV/HV side
-    vector<vector<vector<MStripHit*>>> CollectStripHits(MReadOutAssembly* Event);
+    tuple<vector<vector<vector<MStripHit*>>>, bool> CollectStripHits(MReadOutAssembly* Event);
+
+  //! Divide an event's nearest neighbor strip hits by detector and LV/HV side
+    vector<vector<vector<MStripHit*>>> CollectNearestNeighborStripHits(MReadOutAssembly* Event);
+
+  //! Assign nearest neighbor strip hits to their associated hits
+    void AssignNearestNeighbors(MReadOutAssembly* Event);
 
   //! Read in strip hits on each side for each detector and perform quality selections
     bool EventSelection(MReadOutAssembly* Event, const vector<vector<vector<MStripHit*>>>& StripHits);
@@ -97,6 +103,7 @@ class MModuleStripPairingMultiRoundChiSquare : public MModule
 
   //! Evaluate the reduced chi square for all possible strip pairings
     tuple<vector<vector<unsigned int>>, vector<vector<unsigned int>>, double> EvaluateAllCombinations(unsigned int d, const vector<vector<vector<vector<vector<unsigned int>>>>>& Combinations, const vector<vector<vector<MStripHit*>>>& StripHits);
+
   //! Create hits
     bool CreateHits(unsigned int d, MReadOutAssembly* Event, const vector<vector<vector<MStripHit*>>>& StripHits, const vector<vector<unsigned int>>& BestLVSideCombo, const vector<vector<unsigned int>>& BestHVSideCombo);
     //! Return the order of indices resulting from sorting a vector

include/MStripHit.h

@@ -128,8 +128,13 @@ class MStripHit
   void IsNearestNeighbor(bool NearestNeighbor) { m_IsNearestNeighbor = NearestNeighbor; }
   //! Return whether the strip is a nearest-neighbor hit
   bool IsNearestNeighbor() const { return m_IsNearestNeighbor; }
-
-  //! Set whether the strip has passed the fast threshold
+
+  //! Set if this is an ambiguous neighbor (ie. associated with multiple hits)
+  void IsAmbiguousNearestNeighbor(bool AmbiguousNearestNeighbor) { m_IsAmbiguousNearestNeighbor = AmbiguousNearestNeighbor; }
+  //! Return boolean indicating whether strip is an ambiguous nearest neighbor (default is false for triggered strips)
+  bool IsAmbiguousNearestNeighbor() const { return m_IsAmbiguousNearestNeighbor; }
+
+  //! Set the Fast Timing flag
   void HasFastTiming(bool FastTiming) { m_HasFastTiming = FastTiming; }
   //! Return whether the strip has passed the fast threshold
   bool HasFastTiming() const { return m_HasFastTiming; }
@@ -207,6 +212,9 @@ class MStripHit
   bool m_IsNearestNeighbor;
   //! True if the strip has triggered
   bool m_HasTriggered;
+  //! True if the nearest neighbor strip hit is associated with multiple strip paired hits
+  bool m_IsAmbiguousNearestNeighbor;
+
   //! True if the hit has fast timing
   bool m_HasFastTiming;
   //! True if the hit has calibrated timing

src/MHit.cxx

@@ -81,6 +81,7 @@ void MHit::Clear()
   m_EnergyResolution = g_DoubleNotDefined;
 
   m_StripHits.clear();
+  m_NearestNeighborStripHits.clear();
   m_Origins.clear();
 
   m_CrossTalk = false;
@@ -114,6 +115,24 @@ MStripHit* MHit::GetStripHit(unsigned int i)
 ////////////////////////////////////////////////////////////////////////////////
 
 
+MStripHit* MHit::GetNearestNeighborStripHit(unsigned int i)
+{
+  // Return strip hit i
+
+  if (i < m_NearestNeighborStripHits.size()) {
+    return m_StripHits[i];
+  }
+
+  if (g_Verbosity >= c_Error) cout<<"Error in MHit::GetNearestNeighborStripHit: Strip hit index "<<i<<" is out of bounds: "<<m_NearestNeighborStripHits.size()<<" strip hits available"<<endl;
+
+  return nullptr;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+
+
+
 void MHit::AddStripHit(MStripHit* StripHit)
 {
   // Add a strip hit
@@ -129,6 +148,21 @@ void MHit::AddStripHit(MStripHit* StripHit)
 ////////////////////////////////////////////////////////////////////////////////
 
 
+void MHit::AddNearestNeighborStripHit(MStripHit* StripHit)
+{
+  // Add a strip hit
+
+  if (StripHit != nullptr) {
+    m_NearestNeighborStripHits.push_back(StripHit);
+  } else {
+    if (g_Verbosity >= c_Error) cout<<"Error in MHit::AddNearestNeighborStripHit: Strip hit is nullptr"<<endl;
+  }
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+
+
 void MHit::RemoveStripHit(unsigned int i)
 {
   // Remove a strip hit without deleting it
@@ -177,6 +211,9 @@ bool MHit::StreamDat(ostream& S, int Version)
     for (auto SH : m_StripHits) {
       SH->StreamDat(S, 0);
     }
+    for (auto SH : m_NearestNeighborStripHits) {
+      SH->StreamDat(S, 0);
+    }
   } else {
     if (g_Verbosity >= c_Error) cout<<"Error in MHit::StreamDat: Stream version "<<Version<<" not handled"<<endl;
     return false;

src/MModuleStripPairingMultiRoundChiSquare.cxx

@@ -234,40 +234,93 @@ float MModuleStripPairingMultiRoundChiSquare::ChargeTrappingCorrection(unsigned
 ////////////////////////////////////////////////////////////////////////////////
 
 //! Divide an event's strip hits by detector and LV/HV side
-vector<vector<vector<MStripHit*>>> MModuleStripPairingMultiRoundChiSquare::CollectStripHits(MReadOutAssembly* Event)
+tuple<vector<vector<vector<MStripHit*>>>, bool> MModuleStripPairingMultiRoundChiSquare::CollectStripHits(MReadOutAssembly* Event)
 {
 
   // Split hits by detector ID
   vector<unsigned int> DetectorIDs; // List of detector IDs
   vector<vector<vector<MStripHit*>>> StripHits; // list of detector IDs, list of sides (LV and HV), list of strip hits
+  bool IncludingNearestNeighbors = false;
 
   for (unsigned int sh = 0; sh < Event->GetNStripHits(); ++sh) { // Populate StripHits with this event's strip hits
     MStripHit* SH = Event->GetStripHit(sh);
-    unsigned int Side = (SH->IsLowVoltageStrip() == true) ? 0 : 1;
-
-    // Check if detector is on list
-    bool DetectorFound = false;
-    unsigned int DetectorPos = 0;
-    for (unsigned int d = 0; d < DetectorIDs.size(); ++d) {
-      if (DetectorIDs[d] == SH->GetDetectorID()) {
-        DetectorFound = true;
-        DetectorPos = d;
+
+    // Separate out the triggered and NN strip hits
+    if (SH->IsNearestNeighbor() == false) {
+
+      unsigned int Side = (SH->IsLowVoltageStrip() == true) ? 0 : 1;
+
+      // Check if detector is on list
+      bool DetectorFound = false;
+      unsigned int DetectorPos = 0;
+      for (unsigned int d = 0; d < DetectorIDs.size(); ++d) {
+        if (DetectorIDs[d] == SH->GetDetectorID()) {
+          DetectorFound = true;
+          DetectorPos = d;
+        }
+      }
+
+      // Once the correct detector is found, add strip hit to StripHits
+      if (DetectorFound == true) {
+        StripHits[DetectorPos][Side].push_back(SH);
+      } else { // If encountering a new detector, initialize list of sides/hits corresponding to that detector
+        vector<vector<MStripHit*>> List; // list of sides, list of hits
+        List.push_back(vector<MStripHit*>()); // LV
+        List.push_back(vector<MStripHit*>()); // HV
+        List[Side].push_back(SH);
+        StripHits.push_back(List);
+        DetectorIDs.push_back(SH->GetDetectorID());
       }
     }
+    else {
+      IncludingNearestNeighbors = true;
+    }
+  }
+  return {StripHits, IncludingNearestNeighbors};
+}
+
+////////////////////////////////////////////////////////////////////////////////
 
-    // Once the correct detector is found, add strip hit to StripHits
-    if (DetectorFound == true) {
-      StripHits[DetectorPos][Side].push_back(SH);
-    } else { // If encountering a new detector, initialize list of sides/hits corresponding to that detector
-      vector<vector<MStripHit*>> List; // list of sides, list of hits
-      List.push_back(vector<MStripHit*>()); // LV
-      List.push_back(vector<MStripHit*>()); // HV
-      List[Side].push_back(SH);
-      StripHits.push_back(List);
-      DetectorIDs.push_back(SH->GetDetectorID());
+//! Divide an event's nearest neighbor strip hits by detector and LV/HV side
+vector<vector<vector<MStripHit*>>> MModuleStripPairingMultiRoundChiSquare::CollectNearestNeighborStripHits(MReadOutAssembly* Event)
+{
+
+  // Split hits by detector ID
+  vector<unsigned int> DetectorIDs; // List of detector IDs
+  vector<vector<vector<MStripHit*>>> NNStripHits; // list of detector IDs, list of sides (LV and HV), list of strip hits
+
+  for (unsigned int sh = 0; sh < Event->GetNStripHits(); ++sh) { // Populate StripHits with this event's NN strip hits
+    MStripHit* SH = Event->GetStripHit(sh);
+
+    // Separate out the triggered and NN strip hits
+    if (SH->IsNearestNeighbor() == true) {
+
+      unsigned int Side = (SH->IsLowVoltageStrip() == true) ? 0 : 1;
+
+      // Check if detector is on list
+      bool DetectorFound = false;
+      unsigned int DetectorPos = 0;
+      for (unsigned int d = 0; d < DetectorIDs.size(); ++d) {
+        if (DetectorIDs[d] == SH->GetDetectorID()) {
+          DetectorFound = true;
+          DetectorPos = d;
+        }
+      }
+
+      // Once the correct detector is found, add strip hit to NNStripHits
+      if (DetectorFound == true) {
+        NNStripHits[DetectorPos][Side].push_back(SH);
+      } else { // If encountering a new detector, initialize list of sides/hits corresponding to that detector
+        vector<vector<MStripHit*>> List; // list of sides, list of hits
+        List.push_back(vector<MStripHit*>()); // LV
+        List.push_back(vector<MStripHit*>()); // HV
+        List[Side].push_back(SH);
+        NNStripHits.push_back(List);
+        DetectorIDs.push_back(SH->GetDetectorID());
+      }
     }
   }
-  return StripHits;
+  return NNStripHits;
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -276,7 +329,7 @@ vector<vector<vector<MStripHit*>>> MModuleStripPairingMultiRoundChiSquare::Colle
 bool MModuleStripPairingMultiRoundChiSquare::EventSelection(MReadOutAssembly* Event, const vector<vector<vector<MStripHit*>>>& StripHits)
 {
 
-  // Limit the number of strip hits on each side
+  // Limit the number of (triggered) strip hits on each side
   for (unsigned int d = 0; d < StripHits.size(); ++d) { // Detector loop
     for (unsigned int side = 0; side <= 1; ++side) { // Side loop
       if (StripHits[d][side].size() > m_MaximumStrips) {
@@ -711,6 +764,74 @@ bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOut
 
 ////////////////////////////////////////////////////////////////////////////////
 
+//! Assign nearest neighbor strip hits to their appropriate hit
+void MModuleStripPairingMultiRoundChiSquare::AssignNearestNeighbors(MReadOutAssembly* Event) {
+
+  vector<vector<vector<MStripHit*>>> NNStripHits = CollectNearestNeighborStripHits(Event); // List of detectors, list of sides, list of strip hits
+
+  vector<MStripHit*> AssignedNeighbors; // List of all the assigned NN strip hits, in order to check if NNs are double counted
+
+  for (unsigned int h = 0; h < Event->GetNHits(); h++) {
+    vector<vector<int>> StripIDs; // list of sides, list of strips
+    StripIDs.push_back(vector<int>()); // LV
+    StripIDs.push_back(vector<int>()); // HV
+    bool AssignedDetector = false;
+    int DetectorID; // Define detector ID where hit took place
+    for (unsigned int sh = 0; sh < Event->GetHit(h)->GetNStripHits(); sh++) {
+      // Collect all the strip hits in a hit and split them by side
+      MStripHit* SH = Event->GetHit(h)->GetStripHit(sh);
+      unsigned int Side = (SH->IsLowVoltageStrip() == true) ? 0 : 1;
+      StripIDs[Side].push_back(SH->GetStripID());
+
+      if (AssignedDetector == false) {
+        DetectorID = SH->GetDetectorID();
+        AssignedDetector = true;
+      }
+    }
+    // For each side, find the edge strip hit. i.e if there's charge sharing between strips 4, 5, and 6, the edges will be 4 and 6
+    int LeftEdgeLV = *min_element(StripIDs[0].begin(), StripIDs[0].end());
+    int RightEdgeLV = *max_element(StripIDs[0].begin(), StripIDs[0].end());
+    int LeftEdgeHV = *min_element(StripIDs[1].begin(), StripIDs[1].end());
+    int RightEdgeHV = *max_element(StripIDs[1].begin(), StripIDs[1].end());
+
+    // If there are two hits that are one strip hit apart, then the NN strip hit will be added to both hits
+
+    // Define the LV neighbors
+    for (unsigned int sh = 0; sh < NNStripHits[DetectorID][0].size(); sh++) {
+      MStripHit* NNSH = NNStripHits[DetectorID][0][sh];
+      if ((NNSH->GetStripID() == LeftEdgeLV - 1) or (NNSH->GetStripID() == RightEdgeLV + 1)) {
+        Event->GetHit(h)->AddNearestNeighborStripHit(NNSH);
+        // If NN strip hit is not yet assigned to a hit, then add it to the list of assigned neighbors
+        if (find(AssignedNeighbors.begin(), AssignedNeighbors.end(), NNSH) == AssignedNeighbors.end()) {
+          AssignedNeighbors.push_back(NNSH);
+        }
+        // If it has already been assigned to a hit, then flag that strip hit as an ambiguous nearest neighbor
+        else {
+          NNSH->IsAmbiguousNearestNeighbor(true);
+        }
+      }
+    }
+
+    // Define the HV neighbors
+    for (unsigned int sh = 0; sh < NNStripHits[DetectorID][1].size(); sh++) {
+      MStripHit* NNSH = NNStripHits[DetectorID][1][sh];
+      if ((NNSH->GetStripID() == LeftEdgeHV - 1) or (NNSH->GetStripID() == RightEdgeHV + 1)) {
+        Event->GetHit(h)->AddNearestNeighborStripHit(NNSH);
+      }
+      // If NN strip hit is not yet assigned to a hit, then add it to the list of assigned neighbors
+      if (find(AssignedNeighbors.begin(), AssignedNeighbors.end(), NNSH) == AssignedNeighbors.end()) {
+        AssignedNeighbors.push_back(NNSH);
+      }
+      // If it has already been assigned to a hit, then flag that strip hit as an ambiguous nearest neighbor
+      else {
+        NNSH->IsAmbiguousNearestNeighbor(true);
+      }
+    }
+  }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
 //! Main data analysis routine, which updates the event to a new level
 bool MModuleStripPairingMultiRoundChiSquare::AnalyzeEvent(MReadOutAssembly* Event)
 {
@@ -722,7 +843,7 @@ bool MModuleStripPairingMultiRoundChiSquare::AnalyzeEvent(MReadOutAssembly* Even
   }
 
   // Collect strip hits from input event
-  vector<vector<vector<MStripHit*>>> StripHits = CollectStripHits(Event); // List of detectors, list of sides, list of strip hits
+  auto [StripHits, IncludingNearestNeighbors] = CollectStripHits(Event); // List of detectors, list of sides, list of strip hits (and bool saying if running with nearest neighbors or not
 
   // Perform some event selections
   bool CheckStripHits = EventSelection(Event, StripHits);
@@ -840,6 +961,11 @@ bool MModuleStripPairingMultiRoundChiSquare::AnalyzeEvent(MReadOutAssembly* Even
 
   } // End Detector loop
 
+  // If there are NN strips, assign them to their appropriate hits
+  if (IncludingNearestNeighbors == true) {
+    AssignNearestNeighbors(Event);
+  }
+
   Event->SetAnalysisProgress(MAssembly::c_StripPairing);
 
   return true;