diff --git a/include/MHit.h b/include/MHit.h index a44ac6b6..d0646d5e 100644 --- a/include/MHit.h +++ b/include/MHit.h @@ -47,14 +47,19 @@ class MHit // Strip hits: - //! Return the number of strip hits + //! Return the number of triggered strip hits unsigned int GetNStripHits() const { return m_StripHits.size(); } //! Return strip hit i or nullptr if i is out of bounds //! Ownership stays elsewhere + //! This includes only triggered strip hits MStripHit* GetStripHit(unsigned int i); - //! Add a strip hit + //! Return nearest neighbor strip hit i or nullptr if i is out of bounds + MStripHit* GetNearestNeighborStripHit(unsigned int i); + //! Add a triggered 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 +177,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..a28ab712 100644 --- a/include/MModuleStripPairingMultiRoundChiSquare.h +++ b/include/MModuleStripPairingMultiRoundChiSquare.h @@ -86,8 +86,14 @@ class MModuleStripPairingMultiRoundChiSquare : public MModule //! Function to apply charge trapping correction float ChargeTrappingCorrection(unsigned int d, const vector>& StripHits); - //! Divide an event's strip hits by detector and LV/HV side - vector>> CollectStripHits(MReadOutAssembly* Event); + //! Divide an event's triggered strip hits by detector and LV/HV side + vector>> CollectTriggeredStripHits(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 2e1aa9d9..0719ed56 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; } @@ -200,6 +205,9 @@ class MStripHit bool m_IsGuardRing; //! True if the hit is a nearest neighbor hit bool m_IsNearestNeighbor; + //! 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..e1c09a21 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_NearestNeighborStripHits[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::FindNewCombinations(const vector>>& OldOnes, const vector& StripHits, bool RoundTwo) +vector>> MModuleStripPairingMultiRoundChiSquare::FindNewCombinations(const vector>>& OldOnes, const vector& TriggeredStripHits, bool RoundTwo) { // Define new vector of ints NewOnes vector>> NewOnes; // of of @@ -167,7 +168,7 @@ vector>> MModuleStripPairingMultiRoundChiSquare::Fin // Reserve once since this temporary vector mirrors NewCombinedStrips in size. NewCombinedAsIDs.reserve(NewCombinedStrips.size()); for (unsigned int s = 0; s < NewCombinedStrips.size(); ++s) { - NewCombinedAsIDs.push_back(StripHits[NewCombinedStrips[s]]->GetStripID()); // Translates the hit number to the actual strip ID + NewCombinedAsIDs.push_back(TriggeredStripHits[NewCombinedStrips[s]]->GetStripID()); // Translates the hit number to the actual strip ID } sort(NewCombinedAsIDs.begin(), NewCombinedAsIDs.end()); @@ -234,59 +235,108 @@ float MModuleStripPairingMultiRoundChiSquare::ChargeTrappingCorrection(unsigned //////////////////////////////////////////////////////////////////////////////// //! Divide an event's strip hits by detector and LV/HV side -vector>> MModuleStripPairingMultiRoundChiSquare::CollectStripHits(MReadOutAssembly* Event) +vector>> MModuleStripPairingMultiRoundChiSquare::CollectTriggeredStripHits(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 + vector>> TriggeredStripHits; // list of detector IDs, list of sides (LV and HV), list of triggered strip hits - for (unsigned int sh = 0; sh < Event->GetNStripHits(); ++sh) { // Populate StripHits with this event's strip hits + for (unsigned int sh = 0; sh < Event->GetNStripHits(); ++sh) { // Populate TriggeredStripHits with this event's triggered 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 TriggeredStripHits + if (DetectorFound == true) { + TriggeredStripHits[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); + TriggeredStripHits.push_back(List); + DetectorIDs.push_back(SH->GetDetectorID()); } } + } + return TriggeredStripHits; +} + +//////////////////////////////////////////////////////////////////////////////// + +//! Divide an event's nearest neighbor strip hits by detector and LV/HV side +vector>> MModuleStripPairingMultiRoundChiSquare::CollectNearestNeighborStripHits(MReadOutAssembly* Event) +{ - // 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()); + // 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 nearest neighbor strip hits + + for (unsigned int sh = 0; sh < Event->GetNStripHits(); ++sh) { // Populate StripHits with this event's nearest neighbor strip hits + MStripHit* SH = Event->GetStripHit(sh); + + // Separate out the triggered and nearest neighbor 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; } //////////////////////////////////////////////////////////////////////////////// //! Read in strip hits on each side for each detector and perform quality selections -bool MModuleStripPairingMultiRoundChiSquare::EventSelection(MReadOutAssembly* Event, const vector>>& StripHits) +bool MModuleStripPairingMultiRoundChiSquare::EventSelection(MReadOutAssembly* Event, const vector>>& TriggeredStripHits) { - // Limit the number of strip hits on each side - for (unsigned int d = 0; d < StripHits.size(); ++d) { // Detector loop + // Limit the number of (triggered) strip hits on each side + for (unsigned int d = 0; d < TriggeredStripHits.size(); ++d) { // Detector loop for (unsigned int side = 0; side <= 1; ++side) { // Side loop - if (StripHits[d][side].size() > m_MaximumStrips) { - Event->SetStripPairingError("More than maximum number of strip hits allowed on one side (" + to_string(StripHits[d][side].size()) + ")"); + if (TriggeredStripHits[d][side].size() > m_MaximumStrips) { + Event->SetStripPairingError("More than maximum number of strip hits allowed on one side (" + to_string(TriggeredStripHits[d][side].size()) + ")"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } // Check if one side of the detector has no strip hits - if (StripHits[d][side].size() == 0) { + if (TriggeredStripHits[d][side].size() == 0) { Event->SetStripPairingError("One detector side has no strip hits"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; @@ -300,7 +350,7 @@ bool MModuleStripPairingMultiRoundChiSquare::EventSelection(MReadOutAssembly* Ev //////////////////////////////////////////////////////////////////////////////// //! Find all strip combinations for each detector on LV and HV sides given seed combinations -void MModuleStripPairingMultiRoundChiSquare::FindAllCombinations(unsigned int d, vector>>>>& Combinations, const vector>>& StripHits, bool RoundTwo) +void MModuleStripPairingMultiRoundChiSquare::FindAllCombinations(unsigned int d, vector>>>>& Combinations, const vector>>& TriggeredStripHits, bool RoundTwo) { for (unsigned int side = 0; side <= 1; ++side) { // Side loop (LV and HV) @@ -311,7 +361,7 @@ void MModuleStripPairingMultiRoundChiSquare::FindAllCombinations(unsigned int d, while (CombinationsAdded == true) { CombinationsAdded = false; - NewCombinations = FindNewCombinations(Combinations[d][side], StripHits[d][side], RoundTwo); + NewCombinations = FindNewCombinations(Combinations[d][side], TriggeredStripHits[d][side], RoundTwo); //cout<<"Size: "<>, vector>, double> MModuleStripPairingMultiRoundChiSquare::EvaluateAllCombinations(unsigned int d, const vector>>>>& Combinations, const vector>>& StripHits) +tuple>, vector>, double> MModuleStripPairingMultiRoundChiSquare::EvaluateAllCombinations(unsigned int d, const vector>>>>& Combinations, const vector>>& TriggeredStripHits) { double BestChiSquare = numeric_limits::max(); @@ -387,11 +437,11 @@ tuple>, vector>, double> MModul // Add up LV energy and energy resolution for grouping of strips for (unsigned int entry = 0; entry < LVSideCombo[en].size(); ++entry) { // Entry is on the strip level - LVEnergy += StripHits[d][0][LVSideCombo[en][entry]]->GetEnergy(); - LVResolution += pow(StripHits[d][0][LVSideCombo[en][entry]]->GetEnergyResolution(), 2); + LVEnergy += TriggeredStripHits[d][0][LVSideCombo[en][entry]]->GetEnergy(); + LVResolution += pow(TriggeredStripHits[d][0][LVSideCombo[en][entry]]->GetEnergyResolution(), 2); // Add strip to current hit pairing - CurrentHitPairing[0].push_back(StripHits[d][0][LVSideCombo[en][entry]]); + CurrentHitPairing[0].push_back(TriggeredStripHits[d][0][LVSideCombo[en][entry]]); } // Repeats for HV side @@ -399,11 +449,11 @@ tuple>, vector>, double> MModul double HVResolution = 0; for (unsigned int entry = 0; entry < HVSideCombo[ep].size(); ++entry) { - HVEnergy += StripHits[d][1][HVSideCombo[ep][entry]]->GetEnergy(); - HVResolution += pow(StripHits[d][1][HVSideCombo[ep][entry]]->GetEnergyResolution(), 2); + HVEnergy += TriggeredStripHits[d][1][HVSideCombo[ep][entry]]->GetEnergy(); + HVResolution += pow(TriggeredStripHits[d][1][HVSideCombo[ep][entry]]->GetEnergyResolution(), 2); // Add strip to current hit pairing - CurrentHitPairing[1].push_back(StripHits[d][1][HVSideCombo[ep][entry]]); + CurrentHitPairing[1].push_back(TriggeredStripHits[d][1][HVSideCombo[ep][entry]]); } // Apply charge trapping correction for each LV/HV pairing @@ -439,7 +489,7 @@ tuple>, vector>, double> MModul //////////////////////////////////////////////////////////////////////////////// //! Create hits -bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOutAssembly* Event, const vector>>& StripHits, const vector>& BestLVSideCombo, const vector>& BestHVSideCombo) +bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOutAssembly* Event, const vector>>& TriggeredStripHits, const vector>& BestLVSideCombo, const vector>& BestHVSideCombo) { @@ -481,7 +531,7 @@ bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOut bool AllAdjacent = true; for (unsigned int sh = 0; sh < BestLVSideCombo[h].size() - 1; ++sh) { - if (StripHits[d][0][BestLVSideCombo[h][sh]]->GetStripID() + 1 != StripHits[d][0][BestLVSideCombo[h][sh + 1]]->GetStripID()) { + if (TriggeredStripHits[d][0][BestLVSideCombo[h][sh]]->GetStripID() + 1 != TriggeredStripHits[d][0][BestLVSideCombo[h][sh + 1]]->GetStripID()) { AllAdjacentLV = false; AllAdjacent = false; break; @@ -489,7 +539,7 @@ bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOut } for (unsigned int sh = 0; sh < BestHVSideCombo[h].size() - 1; ++sh) { - if (StripHits[d][1][BestHVSideCombo[h][sh]]->GetStripID() + 1 != StripHits[d][1][BestHVSideCombo[h][sh + 1]]->GetStripID()) { + if (TriggeredStripHits[d][1][BestHVSideCombo[h][sh]]->GetStripID() + 1 != TriggeredStripHits[d][1][BestHVSideCombo[h][sh + 1]]->GetStripID()) { AllAdjacentHV = false; AllAdjacent = false; break; @@ -502,11 +552,11 @@ bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOut // Add up energy and energy resolution for each grouping of strips for (unsigned int sh = 0; sh < BestLVSideCombo[h].size(); ++sh) { //cout<<"x-pos: "<GetNonStripPosition()<GetEnergy(); - LVEnergyRes += StripHits[d][0][BestLVSideCombo[h][sh]]->GetEnergyResolution() * StripHits[d][0][BestLVSideCombo[h][sh]]->GetEnergyResolution(); + LVEnergy += TriggeredStripHits[d][0][BestLVSideCombo[h][sh]]->GetEnergy(); + LVEnergyRes += TriggeredStripHits[d][0][BestLVSideCombo[h][sh]]->GetEnergyResolution() * TriggeredStripHits[d][0][BestLVSideCombo[h][sh]]->GetEnergyResolution(); // Add strip to current hit pairing - CurrentHitPairing[0].push_back(StripHits[d][0][BestLVSideCombo[h][sh]]); + CurrentHitPairing[0].push_back(TriggeredStripHits[d][0][BestLVSideCombo[h][sh]]); } LVEnergyResTotal += LVEnergyRes; @@ -515,11 +565,11 @@ bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOut LVEnergyRes = sqrt(LVEnergyRes); for (unsigned int sh = 0; sh < BestHVSideCombo[h].size(); ++sh) { - HVEnergy += StripHits[d][1][BestHVSideCombo[h][sh]]->GetEnergy(); - HVEnergyRes += StripHits[d][1][BestHVSideCombo[h][sh]]->GetEnergyResolution() * StripHits[d][1][BestHVSideCombo[h][sh]]->GetEnergyResolution(); + HVEnergy += TriggeredStripHits[d][1][BestHVSideCombo[h][sh]]->GetEnergy(); + HVEnergyRes += TriggeredStripHits[d][1][BestHVSideCombo[h][sh]]->GetEnergyResolution() * TriggeredStripHits[d][1][BestHVSideCombo[h][sh]]->GetEnergyResolution(); // Add strip to current hit pairing - CurrentHitPairing[1].push_back(StripHits[d][1][BestHVSideCombo[h][sh]]); + CurrentHitPairing[1].push_back(TriggeredStripHits[d][1][BestHVSideCombo[h][sh]]); } // Apply charge trapping correction for each LV/HV pairing @@ -560,10 +610,10 @@ bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOut // Populate hit with strip hits for (unsigned int sh = 0; sh < BestLVSideCombo[h].size(); ++sh) { - Hit->AddStripHit(StripHits[d][0][BestLVSideCombo[h][sh]]); + Hit->AddStripHit(TriggeredStripHits[d][0][BestLVSideCombo[h][sh]]); } for (unsigned int sh = 0; sh < BestHVSideCombo[h].size(); ++sh) { - Hit->AddStripHit(StripHits[d][1][BestHVSideCombo[h][sh]]); + Hit->AddStripHit(TriggeredStripHits[d][1][BestHVSideCombo[h][sh]]); } // Check for charge sharing on either side @@ -587,9 +637,9 @@ bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOut // Assign hit energy based on energy measured on HV side for (unsigned int sh = 0; sh < BestHVSideCombo[h].size(); ++sh) { - Energy = StripHits[d][1][BestHVSideCombo[h][sh]]->GetEnergy(); + Energy = TriggeredStripHits[d][1][BestHVSideCombo[h][sh]]->GetEnergy(); EnergyTotal += Energy; - EnergyResolution = StripHits[d][1][BestHVSideCombo[h][sh]]->GetEnergyResolution(); + EnergyResolution = TriggeredStripHits[d][1][BestHVSideCombo[h][sh]]->GetEnergyResolution(); HVEnergy = Energy; LVEnergy = Energy; @@ -615,9 +665,9 @@ bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOut HVEnergies.push_back(HVEnergy); // Populate hit with strip hits - Hit->AddStripHit(StripHits[d][1][BestHVSideCombo[h][sh]]); + Hit->AddStripHit(TriggeredStripHits[d][1][BestHVSideCombo[h][sh]]); for (unsigned int sh = 0; sh < BestLVSideCombo[h].size(); ++sh) { - Hit->AddStripHit(StripHits[d][0][BestLVSideCombo[h][sh]]); + Hit->AddStripHit(TriggeredStripHits[d][0][BestLVSideCombo[h][sh]]); } // Check if there's charge sharing on LV side (currently no way to tell if there's charge sharing on side that doesn't have multiple hits on single strip) if (BestLVSideCombo[h].size() > 1) { @@ -634,9 +684,9 @@ bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOut // Assign hit energy based on energy measured on LV side for (unsigned int sh = 0; sh < BestLVSideCombo[h].size(); ++sh) { - Energy = StripHits[d][0][BestLVSideCombo[h][sh]]->GetEnergy(); + Energy = TriggeredStripHits[d][0][BestLVSideCombo[h][sh]]->GetEnergy(); EnergyTotal += Energy; - EnergyResolution = StripHits[d][0][BestLVSideCombo[h][sh]]->GetEnergyResolution(); + EnergyResolution = TriggeredStripHits[d][0][BestLVSideCombo[h][sh]]->GetEnergyResolution(); HVEnergy = Energy; LVEnergy = Energy; @@ -660,9 +710,9 @@ bool MModuleStripPairingMultiRoundChiSquare::CreateHits(unsigned int d, MReadOut HVEnergies.push_back(HVEnergy); // Populate hit with strip hits - Hit->AddStripHit(StripHits[d][0][BestLVSideCombo[h][sh]]); + Hit->AddStripHit(TriggeredStripHits[d][0][BestLVSideCombo[h][sh]]); for (unsigned int sh = 0; sh < BestHVSideCombo[h].size(); ++sh) { - Hit->AddStripHit(StripHits[d][1][BestHVSideCombo[h][sh]]); + Hit->AddStripHit(TriggeredStripHits[d][1][BestHVSideCombo[h][sh]]); } // Check if there's charge sharing on HV side (currently no way to tell if there's charge sharing on side that doesn't have multiple hits on single strip) if (BestHVSideCombo[h].size() > 1) { @@ -711,23 +761,102 @@ 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 NN 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 + // Loop over the triggered strip hits assigned to a hit + for (unsigned int sh = 0; sh < Event->GetHit(h)->GetNStripHits(); sh++) { + // Collect all the triggered 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) { - // Check if there are actually any strip hits + // Check if there are actually any strip hits (triggered or un-triggered) if (Event->GetNStripHits() == 0) { Event->SetStripPairingError("No strip hits"); Event->SetAnalysisProgress(MAssembly::c_StripPairing); return false; } - - // Collect strip hits from input event - vector>> StripHits = CollectStripHits(Event); // List of detectors, list of sides, list of strip hits + + // Flag if the event includes any NN strip hits + bool IncludingNearestNeighbors = false; + for (unsigned int sh = 0; sh < Event->GetNStripHits(); ++sh) { // loop over the strip hits of an event (both triggered and NN) + MStripHit* SH = Event->GetStripHit(sh); + + if (SH->IsNearestNeighbor() == true) { + IncludingNearestNeighbors = true; + break; + } + } + // Collect triggered strip hits from input event + vector>> TriggeredStripHits = CollectTriggeredStripHits(Event); // List of detectors, list of sides, list of triggered strip hits // Perform some event selections - bool CheckStripHits = EventSelection(Event, StripHits); + bool CheckTriggeredStripHits = EventSelection(Event, TriggeredStripHits); - if (CheckStripHits == false) { + if (CheckTriggeredStripHits == false) { return false; } @@ -736,7 +865,7 @@ bool MModuleStripPairingMultiRoundChiSquare::AnalyzeEvent(MReadOutAssembly* Even // Define Combinations as list of: detector IDs, sides, strip combinations, set of grouped strips (ex. charge sharing on adjacent strip), strip vector>>>> Combinations; - for (unsigned int d = 0; d < StripHits.size(); ++d) { // Detector loop + for (unsigned int d = 0; d < TriggeredStripHits.size(); ++d) { // Detector loop Combinations.push_back(vector>>>()); // Create 4D vectors for each detector Combinations[d].push_back(vector>>()); // Create 3D vectors within each detector vector for each side Combinations[d].push_back(vector>>()); @@ -744,7 +873,7 @@ bool MModuleStripPairingMultiRoundChiSquare::AnalyzeEvent(MReadOutAssembly* Even // Create the seed combinations that will be added (and expanded upon) for each side of each detector for (unsigned int s = 0; s <= 1; ++s) { vector> Combination; - for (unsigned int sh = 0; sh < StripHits[d][s].size(); ++sh) { + for (unsigned int sh = 0; sh < TriggeredStripHits[d][s].size(); ++sh) { vector CombinedStrips = { sh }; // Use grouping of single strip hit as seed for subsequent groupings // sort(CombinedStrips.begin(), CombinedStrips.end()); Combination.push_back(CombinedStrips); @@ -753,23 +882,23 @@ bool MModuleStripPairingMultiRoundChiSquare::AnalyzeEvent(MReadOutAssembly* Even } } - for (unsigned int d = 0; d < StripHits.size(); ++d) { // Detector loop + for (unsigned int d = 0; d < TriggeredStripHits.size(); ++d) { // Detector loop // Begin the first round of strip pairing bool RoundTwo = false; // Find all possible combinations based on the above seed combination - FindAllCombinations(d, Combinations, StripHits, RoundTwo); + FindAllCombinations(d, Combinations, TriggeredStripHits, RoundTwo); // Evaluate reduced chi square for all combinations and select best LV/HV combinations - auto [BestLVSideCombo, BestHVSideCombo, BestChiSquare] = EvaluateAllCombinations(d, Combinations, StripHits); + auto [BestLVSideCombo, BestHVSideCombo, BestChiSquare] = EvaluateAllCombinations(d, Combinations, TriggeredStripHits); if (BestChiSquare > ChiSquareThreshold) { RoundTwo = true; // Repeat strip pairing, now allowing groupings of non adjacent strips - FindAllCombinations(d, Combinations, StripHits, RoundTwo); - auto [BestLVSideComboRoundTwo, BestHVSideComboRoundTwo, BestChiSquareRoundTwo] = EvaluateAllCombinations(d, Combinations, StripHits); + FindAllCombinations(d, Combinations, TriggeredStripHits, RoundTwo); + auto [BestLVSideComboRoundTwo, BestHVSideComboRoundTwo, BestChiSquareRoundTwo] = EvaluateAllCombinations(d, Combinations, TriggeredStripHits); // Update best LV/HV combos if a better pairing is found if (BestChiSquareRoundTwo < BestChiSquare) { @@ -803,7 +932,7 @@ bool MModuleStripPairingMultiRoundChiSquare::AnalyzeEvent(MReadOutAssembly* Even int index = BestHVSideCombo.size(); for (unsigned int h = index; h < BestLVSideCombo.size(); h++) { for (unsigned int sh = 0; sh < BestLVSideCombo[h].size(); ++sh) { - UnpairedEnergy += StripHits[d][0][BestLVSideCombo[h][sh]]->GetEnergy(); // Add all the energy on the unpaired strips + UnpairedEnergy += TriggeredStripHits[d][0][BestLVSideCombo[h][sh]]->GetEnergy(); // Add all the energy on the unpaired strips } } Event->SetStripPairing_QualityFlag("Best strip pairing leaves at least one grouping of strips unpaired (unpaired energy: " + to_string(UnpairedEnergy) + ")"); @@ -813,14 +942,14 @@ bool MModuleStripPairingMultiRoundChiSquare::AnalyzeEvent(MReadOutAssembly* Even int index = BestLVSideCombo.size(); for (unsigned int h = index; h < BestHVSideCombo.size(); h++) { for (unsigned int sh = 0; sh < BestHVSideCombo[h].size(); ++sh) { - UnpairedEnergy += StripHits[d][1][BestHVSideCombo[h][sh]]->GetEnergy(); // Add all the energy on the unpaired strips + UnpairedEnergy += TriggeredStripHits[d][1][BestHVSideCombo[h][sh]]->GetEnergy(); // Add all the energy on the unpaired strips } } Event->SetStripPairing_QualityFlag("Best strip pairing leaves at least one grouping of strips unpaired (unpaired energy: " + to_string(UnpairedEnergy) + ")"); } // Populate hits with best strip paired combination - bool PopulateHits = CreateHits(d, Event, StripHits, BestLVSideCombo, BestHVSideCombo); + bool PopulateHits = CreateHits(d, Event, TriggeredStripHits, BestLVSideCombo, BestHVSideCombo); if (PopulateHits == false) { return false; @@ -840,6 +969,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 2efdc4f6..9d55187e 100644 --- a/src/MStripHit.cxx +++ b/src/MStripHit.cxx @@ -86,6 +86,7 @@ void MStripHit::Clear() m_IsGuardRing = false; m_IsNearestNeighbor = false; + m_IsAmbiguousNearestNeighbor = false; m_HasFastTiming = false; m_HasCalibratedTiming = false; @@ -99,7 +100,7 @@ void MStripHit::Clear() bool MStripHit::Parse(const MString& Line, int Version) { - // Parse the hit from a string starting with SH + // Parse the hit from a string starting with SH or NN Clear(); @@ -108,8 +109,9 @@ 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()<<" " - <<(m_IsNearestNeighbor == false)<<" " - <GetDetectorID()<<" " + <<((m_ReadOutElement->IsLowVoltageStrip() == true) ? "l" : "h")<<" " + <GetStripID()<<" " + <<(m_IsNearestNeighbor == false)<<" " + <GetDetectorID()<<" " + <<((m_ReadOutElement->IsLowVoltageStrip() == true) ? "l" : "h")<<" " + <GetStripID()<<" " + <<(m_IsNearestNeighbor == false)<<" " + <