TRestDetectorSignalRecoveryProcess.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 "TRestDetectorSignalRecoveryProcess.h"
 
using namespace std;
 
ClassImp(TRestDetectorSignalRecoveryProcess);
 
TRestDetectorSignalRecoveryProcess::TRestDetectorSignalRecoveryProcess() { Initialize(); }
 
TRestDetectorSignalRecoveryProcess::TRestDetectorSignalRecoveryProcess(const char* configFilename) {
    Initialize();
 
    if (LoadConfigFromFile(configFilename) == -1) LoadDefaultConfig();
 
    PrintMetadata();
}
 
TRestDetectorSignalRecoveryProcess::~TRestDetectorSignalRecoveryProcess() { delete fOutputSignalEvent; }
 
void TRestDetectorSignalRecoveryProcess::LoadDefaultConfig() {
    SetName("removeChannels-Default");
    SetTitle("Default config");
}
 
void TRestDetectorSignalRecoveryProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fInputSignalEvent = nullptr;
    fOutputSignalEvent = new TRestDetectorSignalEvent();
}
 
void TRestDetectorSignalRecoveryProcess::LoadConfig(const string& configFilename, const string& name) {
    if (LoadConfigFromFile(configFilename, name) == -1) LoadDefaultConfig();
}
 
void TRestDetectorSignalRecoveryProcess::InitProcess() {
    fReadout = GetMetadata<TRestDetectorReadout>();
 
    if (fReadout == nullptr) {
        RESTError << "TRestDetectorSignalRecoveryProcess. Readout has not been initialized!" << RESTendl;
        exit(-1);
    }
}
 
TRestEvent* TRestDetectorSignalRecoveryProcess::ProcessEvent(TRestEvent* evInput) {
    fInputSignalEvent = (TRestDetectorSignalEvent*)evInput;
 
    for (int n = 0; n < fInputSignalEvent->GetNumberOfSignals(); n++)
        fOutputSignalEvent->AddSignal(*fInputSignalEvent->GetSignal(n));
 
    Int_t idL;
    Int_t idR;
    for (unsigned int x = 0; x < fChannelIds.size(); x++) {
        Int_t type = GetAdjacentSignalIds(fChannelIds[x], idL, idR);
        RESTDebug << "Channel id : " << fChannelIds[x] << " Left : " << idL << " Right : " << idR << RESTendl;
 
        if (idL == -1 || idR == -1) continue;
 
        TRestDetectorSignal* leftSgnl = fInputSignalEvent->GetSignalById(idL);
        TRestDetectorSignal* rightSgnl = fInputSignalEvent->GetSignalById(idR);
 
 
        if (leftSgnl == nullptr || rightSgnl == nullptr) continue;
 
        TRestDetectorSignal recoveredSignal;
        recoveredSignal.SetID(fChannelIds[x]);
 
        if (type == 1)  // Only one dead channel
        {
            for (int n = 0; n < leftSgnl->GetNumberOfPoints(); n++) {
                recoveredSignal.IncreaseAmplitude(leftSgnl->GetTime(n), leftSgnl->GetData(n) / 2.);
                leftSgnl->IncreaseAmplitude(leftSgnl->GetTime(n), -1. * leftSgnl->GetData(n) / 2.);
            }
 
            for (int n = 0; n < rightSgnl->GetNumberOfPoints(); n++) {
                recoveredSignal.IncreaseAmplitude(rightSgnl->GetTime(n), rightSgnl->GetData(n) / 2.);
                rightSgnl->IncreaseAmplitude(rightSgnl->GetTime(n), -1. * rightSgnl->GetData(n) / 2.);
            }
        } else if (type == 2 || type == 3) {  // We got two dead-channels
            if (type == 2)                    // The other dead channel is the one at the left
            {
                for (int n = 0; n < leftSgnl->GetNumberOfPoints(); n++)
                    recoveredSignal.IncreaseAmplitude(leftSgnl->GetTime(n), leftSgnl->GetData(n) / 6.);
 
                for (int n = 0; n < rightSgnl->GetNumberOfPoints(); n++)
                    recoveredSignal.IncreaseAmplitude(rightSgnl->GetTime(n), 2 * rightSgnl->GetData(n) / 6.);
            }
 
            if (type == 3)  // The other dead channel is the one at the right
            {
                for (int n = 0; n < leftSgnl->GetNumberOfPoints(); n++)
                    recoveredSignal.IncreaseAmplitude(leftSgnl->GetTime(n), 2 * leftSgnl->GetData(n) / 6.);
 
                for (int n = 0; n < rightSgnl->GetNumberOfPoints(); n++)
                    recoveredSignal.IncreaseAmplitude(rightSgnl->GetTime(n), rightSgnl->GetData(n) / 6.);
            }
 
            for (int n = 0; n < leftSgnl->GetNumberOfPoints(); n++)
                leftSgnl->IncreaseAmplitude(leftSgnl->GetTime(n), -1. * leftSgnl->GetData(n) / 4.);
            for (int n = 0; n < rightSgnl->GetNumberOfPoints(); n++)
                rightSgnl->IncreaseAmplitude(rightSgnl->GetTime(n), -1. * rightSgnl->GetData(n) / 4.);
        } else {
            RESTWarning << "Adjacent channels for signal " << fChannelIds[x] << " not found " << RESTendl;
            continue;
        }
 
        if (fOutputSignalEvent->GetSignalIndex(fChannelIds[x]) > 0)
            fOutputSignalEvent->RemoveSignalWithId(fChannelIds[x]);
 
        fOutputSignalEvent->AddSignal(recoveredSignal);
 
        /*cout << "Channel recovered!! " << endl;
        if( leftSgnl != nullptr && rightSgnl != nullptr )
            for( int n = 0; n < nPoints; n++ )
                cout << "Sample " << n << " : " << leftSgnl->GetData(n) << " + " << rightSgnl->GetData(n) << "
        = " << recoveredSignal.GetData(n) << endl;*/
 
        RESTDebug << "Channel recovered!! " << RESTendl;
        if (leftSgnl != nullptr && rightSgnl != nullptr)
            for (int n = 0; n < leftSgnl->GetNumberOfPoints(); n++)
                RESTDebug << "Sample " << n << " : " << leftSgnl->GetData(n) << " + " << rightSgnl->GetData(n)
                          << " = " << recoveredSignal.GetData(n) << RESTendl;
    }
 
    RESTDebug << "Channels after : " << fOutputSignalEvent->GetNumberOfSignals() << RESTendl;
 
    return fOutputSignalEvent;
}
 
int TRestDetectorSignalRecoveryProcess::GetAdjacentSignalIds(Int_t signalId, Int_t& idLeft, Int_t& idRight) {
    idLeft = -1;
    idRight = -1;
 
    for (int p = 0; p < fReadout->GetNumberOfReadoutPlanes(); p++) {
        TRestDetectorReadoutPlane* plane = fReadout->GetReadoutPlane(p);
        for (size_t m = 0; m < plane->GetNumberOfModules(); m++) {
            TRestDetectorReadoutModule* mod = plane->GetModule(m);
            // We iterate over all readout modules searching for the one that contains
            // our signal id
            if (mod->IsDaqIDInside(signalId)) {
                // If we find it we use the readoutModule id, and the signalId
                // corresponding to the physical readout channel
                Int_t readoutChannelID = mod->DaqToReadoutChannel(signalId);
 
                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();
                    return 2;
                }
 
                // If idRight is a dead channel we take the next channel
                if (std::find(fChannelIds.begin(), fChannelIds.end(), idRight) != fChannelIds.end()) {
                    idRight = mod->GetChannel(readoutChannelID + 2)->GetDaqID();
                    return 3;
                }
                return 1;
            }
        }
    }
    return 0;
}