Skip to content

Commit

Permalink
update TRestRawCommonNoiseReductionProcess
Browse files Browse the repository at this point in the history
  • Loading branch information
@lobis
lobis committed May 7, 2024
1 parent 005cfb8 commit 11cba14
Show file tree
Hide file tree
Showing 2 changed files with 106 additions and 57 deletions.
20 changes: 14 additions & 6 deletions inc/TRestRawCommonNoiseReductionProcess.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,9 +23,9 @@
#ifndef RestCore_TRestRawCommonNoiseReductionProcess
#define RestCore_TRestRawCommonNoiseReductionProcess

#include <TRestEventProcess.h>
#include <TRestRawSignalEvent.h>

#include "TRestEventProcess.h"
#include "TRestRawSignal.h"

//! A process to subtract the common channels noise from RawSignal type
Expand All @@ -50,13 +50,13 @@ class TRestRawCommonNoiseReductionProcess : public TRestEventProcess {
/// Minimum number of signals required to apply the process.
Int_t fMinSignalsRequired = 200;

std::string fChannelType;
TRestRawReadoutMetadata* fReadoutMetadata = nullptr; //!

void Initialize() override;

void LoadDefaultConfig();

protected:
// add here the members of your event process

public:
RESTValue GetInputEvent() const override { return fInputEvent; }
RESTValue GetOutputEvent() const override { return fOutputEvent; }
Expand All @@ -72,9 +72,17 @@ class TRestRawCommonNoiseReductionProcess : public TRestEventProcess {
void PrintMetadata() override {
BeginPrintProcess();

if (!fChannelType.empty()) {
RESTMetadata << "channelType : " << fChannelType << RESTendl;
}

RESTMetadata << " mode : [" << fMode << "]";
if (fMode == 0) RESTMetadata << " --> Mode 0 activated." << RESTendl;
if (fMode == 1) RESTMetadata << " --> Mode 1 activated." << RESTendl;
if (fMode == 0) {
RESTMetadata << " --> Mode 0 activated." << RESTendl;
}
if (fMode == 1) {
RESTMetadata << " --> Mode 1 activated." << RESTendl;
}
RESTMetadata << " centerWidth : " << fCenterWidth << RESTendl;
RESTMetadata << "blocks : [" << fBlocks << "]" << RESTendl;
RESTMetadata << " Minimum number of signals : " << fMinSignalsRequired << RESTendl;
Expand Down
143 changes: 92 additions & 51 deletions src/TRestRawCommonNoiseReductionProcess.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,12 +85,12 @@
///
#include "TRestRawCommonNoiseReductionProcess.h"

using namespace std;

#include <algorithm>
#include <iostream>
#include <vector>

using namespace std;

ClassImp(TRestRawCommonNoiseReductionProcess);

///////////////////////////////////////////////
Expand All @@ -113,7 +113,9 @@ TRestRawCommonNoiseReductionProcess::TRestRawCommonNoiseReductionProcess() { Ini
TRestRawCommonNoiseReductionProcess::TRestRawCommonNoiseReductionProcess(const char* configFilename) {
Initialize();

if (LoadConfigFromFile(configFilename)) LoadDefaultConfig();
if (LoadConfigFromFile(configFilename)) {
LoadDefaultConfig();
}
}

///////////////////////////////////////////////
Expand Down Expand Up @@ -154,7 +156,9 @@ void TRestRawCommonNoiseReductionProcess::Initialize() {
/// corresponding TRestGeant4AnalysisProcess section inside the RML.
///
void TRestRawCommonNoiseReductionProcess::LoadConfig(const string& configFilename, const string& name) {
if (LoadConfigFromFile(configFilename, name)) LoadDefaultConfig();
if (LoadConfigFromFile(configFilename, name)) {
LoadDefaultConfig();
}
}

///////////////////////////////////////////////
Expand All @@ -169,38 +173,80 @@ TRestEvent* TRestRawCommonNoiseReductionProcess::ProcessEvent(TRestEvent* inputE
fInputEvent = (TRestRawSignalEvent*)inputEvent;

if (fInputEvent->GetNumberOfSignals() < fMinSignalsRequired) {
for (int sgnl = 0; sgnl < fInputEvent->GetNumberOfSignals(); sgnl++) {
fOutputEvent->AddSignal(*fInputEvent->GetSignal(sgnl));
for (int signal = 0; signal < fInputEvent->GetNumberOfSignals(); signal++) {
fOutputEvent->AddSignal(*fInputEvent->GetSignal(signal));
}
return fOutputEvent;
}

// Event base line determination.
if (fReadoutMetadata == nullptr) {
fReadoutMetadata = fInputEvent->GetReadoutMetadata();
}

if (fReadoutMetadata == nullptr && !fChannelType.empty()) {
cerr << "TRestRawCommonNoiseReductionProcess::ProcessEvent: readout metadata is null, cannot filter "
"the process by signal type"
<< endl;
exit(1);
}

vector<TRestRawSignal*> signalsToProcess;
vector<TRestRawSignal*> signalsToIgnore;

if (fChannelType.empty()) {
for (int signal = 0; signal < fInputEvent->GetNumberOfSignals(); signal++) {
signalsToProcess.push_back(fInputEvent->GetSignal(signal));
}
} else {
for (int signal = 0; signal < fInputEvent->GetNumberOfSignals(); signal++) {
TRestRawSignal* signalPtr = fInputEvent->GetSignal(signal);
const UShort_t signalId = signalPtr->GetSignalID();
const string channelType = fReadoutMetadata->GetTypeForChannelDaqId(signalId);

if (fChannelType == channelType) {
signalsToProcess.push_back(signalPtr);
} else {
signalsToIgnore.push_back(signalPtr);
}
}
}

// add signals to ignore to the event
for (auto signal : signalsToIgnore) {
fOutputEvent->AddSignal(*signal);
}

auto eventToProcess = TRestRawSignalEvent();
for (auto signal : signalsToProcess) {
eventToProcess.AddSignal(*signal);
}

// Event baseline determination.
Double_t baseLineMean = 0;
for (int sgnl = 0; sgnl < fInputEvent->GetNumberOfSignals(); sgnl++) {
fInputEvent->GetSignal(sgnl)->CalculateBaseLine(20, 150);
Double_t baseline = fInputEvent->GetSignal(sgnl)->GetBaseLine();
for (int signal = 0; signal < eventToProcess.GetNumberOfSignals(); signal++) {
eventToProcess.GetSignal(signal)->CalculateBaseLine(20, 150);
Double_t baseline = eventToProcess.GetSignal(signal)->GetBaseLine();
baseLineMean += baseline;
}
Double_t Baseline = baseLineMean / fInputEvent->GetNumberOfSignals();
Double_t Baseline = baseLineMean / eventToProcess.GetNumberOfSignals();

if (fBlocks == 0) {
Int_t N = fInputEvent->GetNumberOfSignals();
Int_t N = eventToProcess.GetNumberOfSignals();

// if (GetVerboseLevel() >= REST_Debug) N = 1;
for (int sgnl = 0; sgnl < N; sgnl++) {
fOutputEvent->AddSignal(*fInputEvent->GetSignal(sgnl));
for (int signal = 0; signal < N; signal++) {
fOutputEvent->AddSignal(*eventToProcess.GetSignal(signal));
}

Int_t nBins = fInputEvent->GetSignal(0)->GetNumberOfPoints();
vector<Double_t> sgnlValues(N, 0.0);
Int_t nBins = eventToProcess.GetSignal(0)->GetNumberOfPoints();
vector<Double_t> signalValues(N, 0.0);

for (Int_t bin = 0; bin < nBins; bin++) {
for (Int_t sgnl = 0; sgnl < N; sgnl++) {
sgnlValues[sgnl] = fOutputEvent->GetSignal(sgnl)->GetRawData(bin);
for (Int_t signal = 0; signal < N; signal++) {
signalValues[signal] = fOutputEvent->GetSignal(signal)->GetRawData(bin);
}

std::sort(sgnlValues.begin(), sgnlValues.end());
std::sort(signalValues.begin(), signalValues.end());

// Sorting the different methods
Int_t begin = 0, middle = 0, end = 0;
Expand All @@ -221,13 +267,13 @@ TRestEvent* TRestRawCommonNoiseReductionProcess::ProcessEvent(TRestEvent* inputE

// Calculation of the correction to be made to each TRestRawSignal
Double_t binCorrection = 0.0;
for (Int_t i = begin; i <= end; i++) binCorrection += sgnlValues[i];
for (Int_t i = begin; i <= end; i++) binCorrection += signalValues[i];

binCorrection = binCorrection / norm;

// Correction applied.
for (Int_t sgnl = 0; sgnl < N; sgnl++)
fOutputEvent->GetSignal(sgnl)->IncreaseBinBy(bin, Baseline - binCorrection);
for (Int_t signal = 0; signal < N; signal++)
fOutputEvent->GetSignal(signal)->IncreaseBinBy(bin, Baseline - binCorrection);
}

return fOutputEvent;
Expand All @@ -246,39 +292,39 @@ TRestEvent* TRestRawCommonNoiseReductionProcess::ProcessEvent(TRestEvent* inputE
nSign = 0;
// if (GetVerboseLevel() >= REST_Debug) N = 1;

for (Int_t sgnl = 0; sgnl < N; sgnl++) {
sigID = firstInBlock + sgnl;
fInputEvent->GetSignalById(sigID)->CalculateBaseLine(20, 500);
if (fInputEvent->GetSignalById(sigID)->GetBaseLineSigma() >= 3.3) {
for (Int_t signal = 0; signal < N; signal++) {
sigID = firstInBlock + signal;
eventToProcess.GetSignalById(sigID)->CalculateBaseLine(20, 500);
if (eventToProcess.GetSignalById(sigID)->GetBaseLineSigma() >= 3.3) {
// debug << "Baseline1: " <<
// fInputEvent->GetSignalById(sigID)->GetBaseLineSigma() <<
// eventToProcess.GetSignalById(sigID)->GetBaseLineSigma() <<
// endl;
fOutputEvent->AddSignal(*fInputEvent->GetSignalById(sigID));
fOutputEvent->AddSignal(*eventToProcess.GetSignalById(sigID));
nSign++;
}
}

Int_t nBins = fInputEvent->GetSignal(0)->GetNumberOfPoints();
vector<Double_t> sgnlValues(nSign, 0.0);
Int_t nBins = eventToProcess.GetSignal(0)->GetNumberOfPoints();
vector<Double_t> signalValues(nSign, 0.0);

// debug << "nSign: " << nSign << endl;

for (Int_t bin = 0; bin < nBins; bin++) {
int i = 0;
for (Int_t sgnl = 0; sgnl < N; sgnl++) {
sigID = firstInBlock + sgnl;
if (fInputEvent->GetSignalById(sigID)->GetBaseLineSigma() >= 3.3) {
for (Int_t signal = 0; signal < N; signal++) {
sigID = firstInBlock + signal;
if (eventToProcess.GetSignalById(sigID)->GetBaseLineSigma() >= 3.3) {
// debug << "Baseline2: " <<
// fInputEvent->GetSignalById(sigID)->GetBaseLineSigma() <<
// eventToProcess.GetSignalById(sigID)->GetBaseLineSigma() <<
// endl;
// debug << fOutputEvent->GetSignalById(sigID)->GetRawData(bin) <<
// endl;
sgnlValues[i] = fOutputEvent->GetSignalById(sigID)->GetRawData(bin);
signalValues[i] = fOutputEvent->GetSignalById(sigID)->GetRawData(bin);
i++;
}
}

std::sort(sgnlValues.begin(), sgnlValues.end());
std::sort(signalValues.begin(), signalValues.end());

// Sorting the different methods
Int_t begin = 0, middle = 0, end = 0;
Expand All @@ -299,21 +345,23 @@ TRestEvent* TRestRawCommonNoiseReductionProcess::ProcessEvent(TRestEvent* inputE

// Calculation of the correction to be made to each TRestRawSignal
Double_t binCorrection = 0.0;
for (Int_t i = begin; i <= end; i++) binCorrection += sgnlValues[i];
for (Int_t i = begin; i <= end; i++) {
binCorrection += signalValues[i];
}

binCorrection = binCorrection / norm;

// Correction applied.
for (Int_t sgnl = 0; sgnl < N; sgnl++) {
if (fInputEvent->GetSignalById(firstInBlock + sgnl)->GetBaseLineSigma() >= 3.3) {
fOutputEvent->GetSignalById(firstInBlock + sgnl)
for (Int_t signal = 0; signal < N; signal++) {
if (eventToProcess.GetSignalById(firstInBlock + signal)->GetBaseLineSigma() >= 3.3) {
fOutputEvent->GetSignalById(firstInBlock + signal)
->IncreaseBinBy(bin, Baseline - binCorrection);
}
}
}
for (int sgnl = 0; sgnl < N; sgnl++) {
if (fInputEvent->GetSignalById(firstInBlock + sgnl)->GetBaseLineSigma() < 3.3) {
fOutputEvent->AddSignal(*fInputEvent->GetSignalById(firstInBlock + sgnl));
for (int signal = 0; signal < N; signal++) {
if (eventToProcess.GetSignalById(firstInBlock + signal)->GetBaseLineSigma() < 3.3) {
fOutputEvent->AddSignal(*eventToProcess.GetSignalById(firstInBlock + signal));
}
}
}
Expand All @@ -326,11 +374,4 @@ TRestEvent* TRestRawCommonNoiseReductionProcess::ProcessEvent(TRestEvent* inputE
/// \brief Function to include required actions after all events have been
/// processed. This method will write the channels histogram.
///
void TRestRawCommonNoiseReductionProcess::EndProcess() {
// Function to be executed once at the end of the process
// (after all events have been processed)

// Start by calling the EndProcess function of the abstract class.
// Comment this if you don't want it.
// TRestEventProcess::EndProcess();
}
void TRestRawCommonNoiseReductionProcess::EndProcess() {}

0 comments on commit 11cba14

Please sign in to comment.