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"
 
using namespace std;
 
#include <algorithm>
#include <iostream>
#include <vector>
 
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 = (TRestRawSignalEvent*)inputEvent;
 
    if (fInputEvent->GetNumberOfSignals() < fMinSignalsRequired) {
        for (int sgnl = 0; sgnl < fInputEvent->GetNumberOfSignals(); sgnl++) {
            fOutputEvent->AddSignal(*fInputEvent->GetSignal(sgnl));
        }
        return fOutputEvent;
    }
 
    // Event base line 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();
        baseLineMean += baseline;
    }
    Double_t Baseline = baseLineMean / fInputEvent->GetNumberOfSignals();
 
    if (fBlocks == 0) {
        Int_t N = fInputEvent->GetNumberOfSignals();
 
        // if (GetVerboseLevel() >= REST_Debug) N = 1;
        for (int sgnl = 0; sgnl < N; sgnl++) {
            fOutputEvent->AddSignal(*fInputEvent->GetSignal(sgnl));
        }
 
        Int_t nBins = fInputEvent->GetSignal(0)->GetNumberOfPoints();
        vector<Double_t> sgnlValues(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);
            }
 
            std::sort(sgnlValues.begin(), sgnlValues.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 += sgnlValues[i];
 
            binCorrection = binCorrection / norm;
 
            // Correction applied.
            for (Int_t sgnl = 0; sgnl < N; sgnl++)
                fOutputEvent->GetSignal(sgnl)->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 sgnl = 0; sgnl < N; sgnl++) {
                sigID = firstInBlock + sgnl;
                fInputEvent->GetSignalById(sigID)->CalculateBaseLine(20, 500);
                if (fInputEvent->GetSignalById(sigID)->GetBaseLineSigma() >= 3.3) {
                    // debug << "Baseline1: " <<
                    // fInputEvent->GetSignalById(sigID)->GetBaseLineSigma() <<
                    // endl;
                    fOutputEvent->AddSignal(*fInputEvent->GetSignalById(sigID));
                    nSign++;
                }
            }
 
            Int_t nBins = fInputEvent->GetSignal(0)->GetNumberOfPoints();
            vector<Double_t> sgnlValues(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) {
                        // debug << "Baseline2: " <<
                        // fInputEvent->GetSignalById(sigID)->GetBaseLineSigma() <<
                        // endl;
                        // debug << fOutputEvent->GetSignalById(sigID)->GetRawData(bin) <<
                        // endl;
                        sgnlValues[i] = fOutputEvent->GetSignalById(sigID)->GetRawData(bin);
                        i++;
                    }
                }
 
                std::sort(sgnlValues.begin(), sgnlValues.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 += sgnlValues[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)
                            ->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));
                }
            }
        }
        return fOutputEvent;
    }
    return nullptr;
}
 
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();
}