TRestRawCommonNoiseReductionProcess.cxx

/*************************************************************************
 * This file is part of the REST software framework.                     *
 *                                                                       *
 * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza)                *
 * For more information see http://gifna.unizar.es/trex                  *
 *                                                                       *
 * REST is free software: you can redistribute it and/or modify          *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation, either version 3 of the License, or     *
 * (at your option) any later version.                                   *
 *                                                                       *
 * REST is distributed in the hope that it will be useful,               *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
 * GNU General Public License for more details.                          *
 *                                                                       *
 * You should have a copy of the GNU General Public License along with   *
 * REST in $REST_PATH/LICENSE.                                           *
 * If not, see http://www.gnu.org/licenses/.                             *
 * For the list of contributors see $REST_PATH/CREDITS.                  *
 *************************************************************************/
 
#include "TRestRawCommonNoiseReductionProcess.h"
 
#include <algorithm>
#include <iostream>
#include <vector>
 
using namespace std;
 
ClassImp(TRestRawCommonNoiseReductionProcess);
 
TRestRawCommonNoiseReductionProcess::TRestRawCommonNoiseReductionProcess() { Initialize(); }
 
TRestRawCommonNoiseReductionProcess::TRestRawCommonNoiseReductionProcess(const char* configFilename) {
    Initialize();
 
    if (LoadConfigFromFile(configFilename)) {
        LoadDefaultConfig();
    }
}
 
TRestRawCommonNoiseReductionProcess::~TRestRawCommonNoiseReductionProcess() {
    delete fInputEvent;
    delete fOutputEvent;
}
 
void TRestRawCommonNoiseReductionProcess::LoadDefaultConfig() { SetTitle("Default config"); }
 
void TRestRawCommonNoiseReductionProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fInputEvent = nullptr;
    fOutputEvent = new TRestRawSignalEvent();
}
 
void TRestRawCommonNoiseReductionProcess::LoadConfig(const string& configFilename, const string& name) {
    if (LoadConfigFromFile(configFilename, name)) {
        LoadDefaultConfig();
    }
}
 
void TRestRawCommonNoiseReductionProcess::InitProcess() {}
 
TRestEvent* TRestRawCommonNoiseReductionProcess::ProcessEvent(TRestEvent* inputEvent) {
    fInputEvent = dynamic_cast<TRestRawSignalEvent*>(inputEvent);
 
    const auto run = GetRunInfo();
    if (run != nullptr) {
        fInputEvent->InitializeReferences(run);
    }
 
    if (fInputEvent->GetNumberOfSignals() < fMinSignalsRequired) {
        for (int signal = 0; signal < fInputEvent->GetNumberOfSignals(); signal++) {
            fOutputEvent->AddSignal(*fInputEvent->GetSignal(signal));
        }
        return fOutputEvent;
    }
 
    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 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 / eventToProcess.GetNumberOfSignals();
 
    if (fBlocks == 0) {
        Int_t N = eventToProcess.GetNumberOfSignals();
 
        // if (GetVerboseLevel() >= REST_Debug) N = 1;
        for (int signal = 0; signal < N; signal++) {
            fOutputEvent->AddSignal(*eventToProcess.GetSignal(signal));
        }
 
        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 signal = 0; signal < N; signal++) {
                signalValues[signal] = fOutputEvent->GetSignal(signal)->GetRawData(bin);
            }
 
            std::sort(signalValues.begin(), signalValues.end());
 
            // Sorting the different methods
            Int_t begin = 0, middle = 0, end = 0;
            middle = (Int_t)N / 2;
            Double_t norm = 1.0;
 
            if (fMode == 0) {
                // We take only the middle one
                begin = (Int_t)((double_t)N / 2.0);
                end = begin;
                norm = 1.;
            } else if (fMode == 1) {
                // We take the average of the TRestDetectorSignals at the Center
                begin = middle - (Int_t)(N * fCenterWidth * 0.01);
                end = middle + (Int_t)(N * fCenterWidth * 0.01);
                norm = (Double_t)end - begin;
            }
 
            // Calculation of the correction to be made to each TRestRawSignal
            Double_t binCorrection = 0.0;
            for (Int_t i = begin; i <= end; i++) binCorrection += signalValues[i];
 
            binCorrection = binCorrection / norm;
 
            // Correction applied.
            for (Int_t signal = 0; signal < N; signal++)
                fOutputEvent->GetSignal(signal)->IncreaseBinBy(bin, Baseline - binCorrection);
        }
 
        return fOutputEvent;
    } else if (fBlocks == 1) {
        Int_t N = 68;
        Int_t nBlocks = 8;
        Int_t firstID = 578;
        Int_t gap = 4;
 
        Int_t firstInBlock;
        Int_t nSign;
        Int_t sigID;
 
        for (int block = 0; block < nBlocks; block++) {
            firstInBlock = firstID + block * (N + gap);
            nSign = 0;
            // if (GetVerboseLevel() >= REST_Debug) N = 1;
 
            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: " <<
                    // eventToProcess.GetSignalById(sigID)->GetBaseLineSigma() <<
                    // endl;
                    fOutputEvent->AddSignal(*eventToProcess.GetSignalById(sigID));
                    nSign++;
                }
            }
 
            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 signal = 0; signal < N; signal++) {
                    sigID = firstInBlock + signal;
                    if (eventToProcess.GetSignalById(sigID)->GetBaseLineSigma() >= 3.3) {
                        // debug << "Baseline2: " <<
                        // eventToProcess.GetSignalById(sigID)->GetBaseLineSigma() <<
                        // endl;
                        // debug << fOutputEvent->GetSignalById(sigID)->GetRawData(bin) <<
                        // endl;
                        signalValues[i] = fOutputEvent->GetSignalById(sigID)->GetRawData(bin);
                        i++;
                    }
                }
 
                std::sort(signalValues.begin(), signalValues.end());
 
                // Sorting the different methods
                Int_t begin = 0, middle = 0, end = 0;
                middle = (Int_t)nSign / 2;
                Double_t norm = 1.0;
 
                if (fMode == 0) {
                    // We take only the middle one
                    begin = (Int_t)((double_t)nSign / 2.0);
                    end = begin;
                    norm = 1.;
                } else if (fMode == 1) {
                    // We take the average of the TRestDetectorSignals at the Center
                    begin = middle - (Int_t)(nSign * fCenterWidth * 0.01);
                    end = middle + (Int_t)(nSign * fCenterWidth * 0.01);
                    norm = (Double_t)end - begin;
                }
 
                // Calculation of the correction to be made to each TRestRawSignal
                Double_t binCorrection = 0.0;
                for (Int_t i = begin; i <= end; i++) {
                    binCorrection += signalValues[i];
                }
 
                binCorrection = binCorrection / norm;
 
                // Correction applied.
                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 signal = 0; signal < N; signal++) {
                if (eventToProcess.GetSignalById(firstInBlock + signal)->GetBaseLineSigma() < 3.3) {
                    fOutputEvent->AddSignal(*eventToProcess.GetSignalById(firstInBlock + signal));
                }
            }
        }
        return fOutputEvent;
    }
    return nullptr;
}
 
void TRestRawCommonNoiseReductionProcess::EndProcess() {}