diff --git a/include/MHit.h b/include/MHit.h index a44ac6b6..8398a9a9 100644 --- a/include/MHit.h +++ b/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 m_StripHits; + + //! List of nearest neighbor strip hits associated with this hit + vector m_NearestNeighborStripHits; //! Position of the hit MVector m_Position; diff --git a/include/MModuleStripPairingMultiRoundChiSquare.h b/include/MModuleStripPairingMultiRoundChiSquare.h index ccdb8e72..24b3a3d8 100644 --- a/include/MModuleStripPairingMultiRoundChiSquare.h +++ b/include/MModuleStripPairingMultiRoundChiSquare.h @@ -87,7 +87,13 @@ class MModuleStripPairingMultiRoundChiSquare : public MModule float ChargeTrappingCorrection(unsigned int d, const vector>& StripHits); //! Divide an event's strip hits by detector and LV/HV side - vector>> CollectStripHits(MReadOutAssembly* Event); + tuple>>, bool> CollectStripHits(MReadOutAssembly* Event); + + //! Divide an event's nearest neighbor strip hits by detector and LV/HV side + vector>> 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>>& StripHits); @@ -97,6 +103,7 @@ class MModuleStripPairingMultiRoundChiSquare : public MModule //! Evaluate the reduced chi square for all possible strip pairings tuple>, vector>, double> EvaluateAllCombinations(unsigned int d, const vector>>>>& Combinations, const vector>>& StripHits); + //! Create hits bool CreateHits(unsigned int d, MReadOutAssembly* Event, const vector>>& StripHits, const vector>& BestLVSideCombo, const vector>& BestHVSideCombo); //! Return the order of indices resulting from sorting a vector diff --git a/include/MStripHit.h b/include/MStripHit.h index 5e35ad44..36952849 100644 --- a/include/MStripHit.h +++ b/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 diff --git a/src/MHit.cxx b/src/MHit.cxx index 260352c4..bcf97ae3 100644 --- a/src/MHit.cxx +++ b/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 "<= c_Error) cout<<"Error in MHit::AddNearestNeighborStripHit: Strip hit is nullptr"<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 "<>> MModuleStripPairingMultiRoundChiSquare::CollectStripHits(MReadOutAssembly* Event) +tuple>>, bool> MModuleStripPairingMultiRoundChiSquare::CollectStripHits(MReadOutAssembly* Event) { // Split hits by detector ID vector DetectorIDs; // List of detector IDs vector>> 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> List; // list of sides, list of hits + List.push_back(vector()); // LV + List.push_back(vector()); // 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> List; // list of sides, list of hits - List.push_back(vector()); // LV - List.push_back(vector()); // 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>> MModuleStripPairingMultiRoundChiSquare::CollectNearestNeighborStripHits(MReadOutAssembly* Event) +{ + + // Split hits by detector ID + vector DetectorIDs; // List of detector IDs + vector>> 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> List; // list of sides, list of hits + List.push_back(vector()); // LV + List.push_back(vector()); // HV + List[Side].push_back(SH); + NNStripHits.push_back(List); + DetectorIDs.push_back(SH->GetDetectorID()); + } } } - return StripHits; + return NNStripHits; } //////////////////////////////////////////////////////////////////////////////// @@ -276,7 +329,7 @@ vector>> MModuleStripPairingMultiRoundChiSquare::Colle bool MModuleStripPairingMultiRoundChiSquare::EventSelection(MReadOutAssembly* Event, const vector>>& 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>> NNStripHits = CollectNearestNeighborStripHits(Event); // List of detectors, list of sides, list of strip hits + + vector 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> StripIDs; // list of sides, list of strips + StripIDs.push_back(vector()); // LV + StripIDs.push_back(vector()); // 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>> 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; diff --git a/src/MStripHit.cxx b/src/MStripHit.cxx index 9d165168..cec82b01 100644 --- a/src/MStripHit.cxx +++ b/src/MStripHit.cxx @@ -87,6 +87,7 @@ void MStripHit::Clear() m_IsGuardRing = false; m_IsNearestNeighbor = false; + m_IsAmbiguousNearestNeighbor = false; m_HasFastTiming = false; m_HasCalibratedTiming = false; @@ -109,8 +110,8 @@ bool MStripHit::Parse(const MString& Line, int Version) if (g_Verbosity >= c_Error) cout<<"Error in MStripHit::Parse: line too short with length "<= c_Error) cout<<"Error in MStripHit::Parse: line starts with '"<= c_Error) cout<<"Error in MStripHit::Parse: line starts with '"<& Origins) bool MStripHit::StreamDat(ostream& S, int Version) { - // Stream the strip hit in Nuclearizer's DAT format - - S<<"SH " - <GetDetectorID()<<" " - <<((m_ReadOutElement->IsLowVoltageStrip() == true) ? "l" : "h")<<" " - <GetStripID()<<" " - <GetDetectorID()<<" " + <<((m_ReadOutElement->IsLowVoltageStrip() == true) ? "l" : "h")<<" " + <GetStripID()<<" " + <GetDetectorID()<<" " + <<((m_ReadOutElement->IsLowVoltageStrip() == true) ? "l" : "h")<<" " + <GetStripID()<<" " + <