From 58c9dcf868d228c4a38057bbd1efc715b2631bfd Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Thu, 2 Mar 2023 12:02:37 +0100 Subject: [PATCH 1/6] GetAdjacentSignalIds now supports two dead channel identification --- src/TRestDetectorSignalRecoveryProcess.cxx | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/TRestDetectorSignalRecoveryProcess.cxx b/src/TRestDetectorSignalRecoveryProcess.cxx index 051f6f65..584c2423 100644 --- a/src/TRestDetectorSignalRecoveryProcess.cxx +++ b/src/TRestDetectorSignalRecoveryProcess.cxx @@ -207,6 +207,10 @@ TRestEvent* TRestDetectorSignalRecoveryProcess::ProcessEvent(TRestEvent* evInput return fOutputSignalEvent; } +/////////////////////////////////////////////// +/// \brief It returns the channel daq id of the adjacent readout channels. It will properly +/// identify that we got two dead channels, but no more than two. +/// void TRestDetectorSignalRecoveryProcess::GetAdjacentSignalIds(Int_t signalId, Int_t& idLeft, Int_t& idRight) { idLeft = -1; idRight = -1; @@ -225,6 +229,14 @@ void TRestDetectorSignalRecoveryProcess::GetAdjacentSignalIds(Int_t signalId, In idLeft = mod->GetChannel(readoutChannelID - 1)->GetDaqID(); idRight = mod->GetChannel(readoutChannelID + 1)->GetDaqID(); + // If idLeft is a dead channel we take the previous channel + if (std::find(fChannelIds.begin(), fChannelIds.end(), idLeft) != fChannelIds.end()) + idLeft = mod->GetChannel(readoutChannelID - 2)->GetDaqID(); + + // If idRight is a dead channel we take the next channel + if (std::find(fChannelIds.begin(), fChannelIds.end(), idLeft) != fChannelIds.end()) + idRight = mod->GetChannel(readoutChannelID + 2)->GetDaqID(); + return; } } From f026925836ab9d385bc9d2fa7ce53e981706d1bc Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Thu, 2 Mar 2023 12:17:21 +0100 Subject: [PATCH 2/6] TRestDetectorSignalRecoveryProcess::GetAdjacentSignalIds now returns a status integer --- inc/TRestDetectorSignalRecoveryProcess.h | 2 +- src/TRestDetectorSignalRecoveryProcess.cxx | 26 ++++++++++++++++++---- 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/inc/TRestDetectorSignalRecoveryProcess.h b/inc/TRestDetectorSignalRecoveryProcess.h index 472719ac..f34b4824 100644 --- a/inc/TRestDetectorSignalRecoveryProcess.h +++ b/inc/TRestDetectorSignalRecoveryProcess.h @@ -46,7 +46,7 @@ class TRestDetectorSignalRecoveryProcess : public TRestEventProcess { void LoadDefaultConfig(); - void GetAdjacentSignalIds(Int_t signalId, Int_t& idLeft, Int_t& idRight); + int GetAdjacentSignalIds(Int_t signalId, Int_t& idLeft, Int_t& idRight); public: any GetInputEvent() const override { return fInputSignalEvent; } diff --git a/src/TRestDetectorSignalRecoveryProcess.cxx b/src/TRestDetectorSignalRecoveryProcess.cxx index 584c2423..f87621c2 100644 --- a/src/TRestDetectorSignalRecoveryProcess.cxx +++ b/src/TRestDetectorSignalRecoveryProcess.cxx @@ -58,10 +58,16 @@ /// History of developments: /// /// 2017-November: First implementation of TRestDetectorSignalRecoveryProcess. +/// Originally it was TRestRawSignalRecoveryProcess. /// Javier Galan /// +/// 2023-February: Update of the process so that it preserves energy and takes into +/// account two consecutive dead channels. +/// Ana Quintana/Javier Galan +/// /// \class TRestDetectorSignalRecoveryProcess /// \author Javier Galan +/// \author Ana Quintana /// ///
/// @@ -211,7 +217,14 @@ TRestEvent* TRestDetectorSignalRecoveryProcess::ProcessEvent(TRestEvent* evInput /// \brief It returns the channel daq id of the adjacent readout channels. It will properly /// identify that we got two dead channels, but no more than two. /// -void TRestDetectorSignalRecoveryProcess::GetAdjacentSignalIds(Int_t signalId, Int_t& idLeft, Int_t& idRight) { +/// It returns a value that determines the case we are considering. +/// +/// 0: Something went wrong (channels not found?) +/// 1: There is only one dead channel +/// 2: There is an additional dead channel to the left +/// 3: There is an additional dead channel to the right +/// +int TRestDetectorSignalRecoveryProcess::GetAdjacentSignalIds(Int_t signalId, Int_t& idLeft, Int_t& idRight) { idLeft = -1; idRight = -1; @@ -230,15 +243,20 @@ void TRestDetectorSignalRecoveryProcess::GetAdjacentSignalIds(Int_t signalId, In idRight = mod->GetChannel(readoutChannelID + 1)->GetDaqID(); // If idLeft is a dead channel we take the previous channel - if (std::find(fChannelIds.begin(), fChannelIds.end(), idLeft) != fChannelIds.end()) + if (std::find(fChannelIds.begin(), fChannelIds.end(), idLeft) != fChannelIds.end()) { idLeft = mod->GetChannel(readoutChannelID - 2)->GetDaqID(); + return 2; + } // If idRight is a dead channel we take the next channel - if (std::find(fChannelIds.begin(), fChannelIds.end(), idLeft) != fChannelIds.end()) + if (std::find(fChannelIds.begin(), fChannelIds.end(), idLeft) != fChannelIds.end()) { idRight = mod->GetChannel(readoutChannelID + 2)->GetDaqID(); + return 3; + } - return; + return 1; } } } + return 0; } From d660e71aee2ee3c655677b664a722325fe4181d9 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Thu, 2 Mar 2023 12:23:55 +0100 Subject: [PATCH 3/6] TRestDetectorSignalRecoveryProcess. Adding documentation with what the method will do now --- src/TRestDetectorSignalRecoveryProcess.cxx | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/TRestDetectorSignalRecoveryProcess.cxx b/src/TRestDetectorSignalRecoveryProcess.cxx index f87621c2..c8473b0c 100644 --- a/src/TRestDetectorSignalRecoveryProcess.cxx +++ b/src/TRestDetectorSignalRecoveryProcess.cxx @@ -43,6 +43,22 @@ /// using the charge of the adjacent readout channels, /// \f$s_i = 0.5 \times (s_{i-1} + s_{i+1})\f$ /// +/// Now we preserve the energy of the event, we apply the following +/// correction to the dead channel (i) and the adjacent channels. +/// +/// \f$s_{i-1} = 0.5 * s_{i-1}\f$ +/// \f$s_i = 0.5 \times (s_{i-1} + s_{i+1})\f$ +/// \f$s_{i+1} = 0.5 * s_{i+1}\f$ +// +/// The algorithm will also identify in case we have up to two +/// consecutive dead channels. In that case it will apply the following +/// correction, where `i` and `i+1` are the dead channels. +/// +/// \f$s_{i-1} = 0.5 * s_{i-1}\f$ +/// \f$s_i = 1/6 \times (2*s_{i-1} + s_{i+2})\f$ +/// \f$s_{i+1} = 1/6 \times (s_{i-1} + 2*s_{i+2})\f$ +/// \f$s_{i+2} = 0.5 * s_{i+2}\f$ +/// /// This process will access the information of the decoding stored in /// the TRestDetectorReadout definition to assure that the righ signal ids, /// corresponding to the adjacent channels, are used in the calculation. From 0d7d8b1242cad36dab6d6c72ec5c28efd612ec7a Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Thu, 2 Mar 2023 12:39:28 +0100 Subject: [PATCH 4/6] TRestDetectorSignalRecoveryProcess. Adding logic to preserve energy and consider two dead channels --- src/TRestDetectorSignalRecoveryProcess.cxx | 44 +++++++++++++++++++--- 1 file changed, 39 insertions(+), 5 deletions(-) diff --git a/src/TRestDetectorSignalRecoveryProcess.cxx b/src/TRestDetectorSignalRecoveryProcess.cxx index c8473b0c..c9afae6c 100644 --- a/src/TRestDetectorSignalRecoveryProcess.cxx +++ b/src/TRestDetectorSignalRecoveryProcess.cxx @@ -187,7 +187,7 @@ TRestEvent* TRestDetectorSignalRecoveryProcess::ProcessEvent(TRestEvent* evInput Int_t idL; Int_t idR; for (unsigned int x = 0; x < fChannelIds.size(); x++) { - GetAdjacentSignalIds(fChannelIds[x], idL, idR); + Int_t type = GetAdjacentSignalIds(fChannelIds[x], idL, idR); RESTDebug << "Channel id : " << fChannelIds[x] << " Left : " << idL << " Right : " << idR << RESTendl; if (fOutputSignalEvent->GetSignalIndex(fChannelIds[x]) > 0) @@ -198,19 +198,53 @@ TRestEvent* TRestDetectorSignalRecoveryProcess::ProcessEvent(TRestEvent* evInput TRestDetectorSignal* leftSgnl = fInputSignalEvent->GetSignalById(idL); TRestDetectorSignal* rightSgnl = fInputSignalEvent->GetSignalById(idR); - if (leftSgnl == nullptr && rightSgnl == nullptr) continue; + /// If the dead channel has no charge on left and right then we do not + /// correct. We could think about correction here? It means it is + /// the first or last channel. + /// + if (leftSgnl == nullptr || rightSgnl == nullptr) continue; TRestDetectorSignal* recoveredSignal = new TRestDetectorSignal(); recoveredSignal->SetID(fChannelIds[x]); - if (leftSgnl != nullptr) { + if (type == 1) // Only one dead channel + { for (int n = 0; n < leftSgnl->GetNumberOfPoints(); n++) recoveredSignal->IncreaseAmplitude(leftSgnl->GetPoint(n) / 2.); - } - if (rightSgnl != nullptr) { for (int n = 0; n < rightSgnl->GetNumberOfPoints(); n++) recoveredSignal->IncreaseAmplitude(rightSgnl->GetPoint(n) / 2.); + + /// We removed the charge that we place at the dead channel + /// This could be optional using a metadata parameter + leftSgnl->IncreaseAmplitude(-leftSgnl->GetPoint(n) / 2); + rightSgnl->IncreaseAmplitude(-rightSgnl->GetPoint(n) / 2); + } + + if (type == 2 || type == 3) // We got two dead-channels + { + if (type == 2) // The dead channel is the one at the right + { + for (int n = 0; n < leftSgnl->GetNumberOfPoints(); n++) + recoveredSignal->IncreaseAmplitude(leftSgnl->GetPoint(n) / 6.); + + for (int n = 0; n < rightSgnl->GetNumberOfPoints(); n++) + recoveredSignal->IncreaseAmplitude(2 * rightSgnl->GetPoint(n) / 6.); + } + + if (type == 3) // The dead channel is the one at the left + { + for (int n = 0; n < leftSgnl->GetNumberOfPoints(); n++) + recoveredSignal->IncreaseAmplitude(2 * leftSgnl->GetPoint(n) / 6.); + + for (int n = 0; n < rightSgnl->GetNumberOfPoints(); n++) + recoveredSignal->IncreaseAmplitude(rightSgnl->GetPoint(n) / 6.); + } + + /// We removed the charge that we place at the dead channel + /// In this case we remove a 25% because we will enter twice in this loop + leftSgnl->IncreaseAmplitude(-leftSgnl->GetPoint(n) / 4); + rightSgnl->IncreaseAmplitude(-rightSgnl->GetPoint(n) / 4); } fOutputSignalEvent->AddSignal(*recoveredSignal); From b0d85f89282020580fbb1694ce89cee1fcbc6ed5 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Thu, 2 Mar 2023 12:49:24 +0100 Subject: [PATCH 5/6] TRestDetectorSignalRecoveryProcess. Fixing compilation errors --- src/TRestDetectorSignalRecoveryProcess.cxx | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/TRestDetectorSignalRecoveryProcess.cxx b/src/TRestDetectorSignalRecoveryProcess.cxx index c9afae6c..85d6aed6 100644 --- a/src/TRestDetectorSignalRecoveryProcess.cxx +++ b/src/TRestDetectorSignalRecoveryProcess.cxx @@ -209,16 +209,17 @@ TRestEvent* TRestDetectorSignalRecoveryProcess::ProcessEvent(TRestEvent* evInput if (type == 1) // Only one dead channel { - for (int n = 0; n < leftSgnl->GetNumberOfPoints(); n++) + for (int n = 0; n < leftSgnl->GetNumberOfPoints(); n++) { recoveredSignal->IncreaseAmplitude(leftSgnl->GetPoint(n) / 2.); + /// Energy preserved. This could be optional using a new metadata member + leftSgnl->IncreaseAmplitude(-leftSgnl->GetPoint(n) / 2); + } - for (int n = 0; n < rightSgnl->GetNumberOfPoints(); n++) + for (int n = 0; n < rightSgnl->GetNumberOfPoints(); n++) { recoveredSignal->IncreaseAmplitude(rightSgnl->GetPoint(n) / 2.); - - /// We removed the charge that we place at the dead channel - /// This could be optional using a metadata parameter - leftSgnl->IncreaseAmplitude(-leftSgnl->GetPoint(n) / 2); - rightSgnl->IncreaseAmplitude(-rightSgnl->GetPoint(n) / 2); + /// Energy preserved. This could be optional using a new metadata member + rightSgnl->IncreaseAmplitude(-rightSgnl->GetPoint(n) / 2); + } } if (type == 2 || type == 3) // We got two dead-channels @@ -243,8 +244,10 @@ TRestEvent* TRestDetectorSignalRecoveryProcess::ProcessEvent(TRestEvent* evInput /// We removed the charge that we place at the dead channel /// In this case we remove a 25% because we will enter twice in this loop - leftSgnl->IncreaseAmplitude(-leftSgnl->GetPoint(n) / 4); - rightSgnl->IncreaseAmplitude(-rightSgnl->GetPoint(n) / 4); + for (int n = 0; n < leftSgnl->GetNumberOfPoints(); n++) + leftSgnl->IncreaseAmplitude(-leftSgnl->GetPoint(n) / 4); + for (int n = 0; n < rightSgnl->GetNumberOfPoints(); n++) + rightSgnl->IncreaseAmplitude(-rightSgnl->GetPoint(n) / 4); } fOutputSignalEvent->AddSignal(*recoveredSignal); From 75dc330e8149ffbe224cc70679e046914afc0153 Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Thu, 2 Mar 2023 13:05:17 +0100 Subject: [PATCH 6/6] TRestDetectorSignalRecoveryProcess. Fixing a compilation problem --- src/TRestDetectorSignalRecoveryProcess.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/TRestDetectorSignalRecoveryProcess.cxx b/src/TRestDetectorSignalRecoveryProcess.cxx index 85d6aed6..e2ecd582 100644 --- a/src/TRestDetectorSignalRecoveryProcess.cxx +++ b/src/TRestDetectorSignalRecoveryProcess.cxx @@ -212,13 +212,13 @@ TRestEvent* TRestDetectorSignalRecoveryProcess::ProcessEvent(TRestEvent* evInput for (int n = 0; n < leftSgnl->GetNumberOfPoints(); n++) { recoveredSignal->IncreaseAmplitude(leftSgnl->GetPoint(n) / 2.); /// Energy preserved. This could be optional using a new metadata member - leftSgnl->IncreaseAmplitude(-leftSgnl->GetPoint(n) / 2); + leftSgnl->IncreaseAmplitude(-1. * leftSgnl->GetPoint(n) / 2); } for (int n = 0; n < rightSgnl->GetNumberOfPoints(); n++) { recoveredSignal->IncreaseAmplitude(rightSgnl->GetPoint(n) / 2.); /// Energy preserved. This could be optional using a new metadata member - rightSgnl->IncreaseAmplitude(-rightSgnl->GetPoint(n) / 2); + rightSgnl->IncreaseAmplitude(-1. * rightSgnl->GetPoint(n) / 2); } } @@ -245,9 +245,9 @@ TRestEvent* TRestDetectorSignalRecoveryProcess::ProcessEvent(TRestEvent* evInput /// We removed the charge that we place at the dead channel /// In this case we remove a 25% because we will enter twice in this loop for (int n = 0; n < leftSgnl->GetNumberOfPoints(); n++) - leftSgnl->IncreaseAmplitude(-leftSgnl->GetPoint(n) / 4); + leftSgnl->IncreaseAmplitude(-1. * leftSgnl->GetPoint(n) / 4); for (int n = 0; n < rightSgnl->GetNumberOfPoints(); n++) - rightSgnl->IncreaseAmplitude(-rightSgnl->GetPoint(n) / 4); + rightSgnl->IncreaseAmplitude(-1. * rightSgnl->GetPoint(n) / 4); } fOutputSignalEvent->AddSignal(*recoveredSignal);