TRestRawVetoAnalysisProcess.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 "TRestRawVetoAnalysisProcess.h"
 
using namespace std;
 
ClassImp(TRestRawVetoAnalysisProcess);
 
TRestRawVetoAnalysisProcess::TRestRawVetoAnalysisProcess() { Initialize(); }
 
TRestRawVetoAnalysisProcess::TRestRawVetoAnalysisProcess(const char* configFilename) {
    Initialize();
 
    LoadConfig(configFilename);
}
 
TRestRawVetoAnalysisProcess::~TRestRawVetoAnalysisProcess() {}
 
void TRestRawVetoAnalysisProcess::LoadDefaultConfig() {
    SetName(this->ClassName());
    SetTitle("Default config");
}
 
void TRestRawVetoAnalysisProcess::LoadConfig(const string& configFilename, const string& name) {
    if (LoadConfigFromFile(configFilename, name)) LoadDefaultConfig();
}
 
void TRestRawVetoAnalysisProcess::InitProcess() {}
 
void TRestRawVetoAnalysisProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fSignalEvent = nullptr;
}
 
TRestEvent* TRestRawVetoAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
    fSignalEvent = (TRestRawSignalEvent*)inputEvent;
 
    const auto run = GetRunInfo();
    if (run != nullptr) {
        fSignalEvent->InitializeReferences(run);
    }
 
    map<int, Double_t> VetoMaxPeakAmplitude_map;
    map<int, Double_t> VetoPeakTime_map;
 
    Int_t VetoAboveThreshold = 0;
    Int_t NVetoAboveThreshold = 0;
    Int_t VetoInTimeWindow = 0;
    Int_t NVetoInTimeWindow = 0;
 
    fSignalEvent->SetRange(fRange);
 
    VetoMaxPeakAmplitude_map.clear();
    VetoPeakTime_map.clear();
 
    // **************************************************************
    // if list of veto Ids without groups is given ******************
    // **************************************************************
 
    if (fReadoutMetadata == nullptr) {
        fReadoutMetadata = fSignalEvent->GetReadoutMetadata();
    }
 
    if (fReadoutMetadata != nullptr && fVetoSignalId[0] == -1) {
        // iterate over signals and get type
        fVetoSignalId.clear();
        for (int i = 0; i < fSignalEvent->GetNumberOfSignals(); i++) {
            TRestRawSignal* signal = fSignalEvent->GetSignal(i);
            auto channelType = fReadoutMetadata->GetTypeForChannelDaqId(signal->GetID());
            // uppercase
            transform(channelType.begin(), channelType.end(), channelType.begin(), ::toupper);
            if (channelType == "VETO") {
                fVetoSignalId.push_back(signal->GetID());
            }
        }
    }
 
    if (fVetoSignalId[0] != -1) {
        // iterate over vetoes
        for (double i : fVetoSignalId) {
            // Checks if channel (fVetoSignalId) participated in the event. If not, it
            // is -1
            if (fSignalEvent->GetSignalIndex(i) != -1) {
                // We extract the parameters from the veto signal
                TRestRawSignal* signal = fSignalEvent->GetSignalById(i);
                // Deal with noise
                signal->CalculateBaseLine(fBaseLineRange.X(), fBaseLineRange.Y(), "ROBUST");
                signal->InitializePointsOverThreshold(TVector2(fPointThreshold, fSignalThreshold),
                                                      fPointsOverThreshold);
 
                // Save two maps with (veto panel ID, max amplitude) and (veto panel ID, peak time)
                if (signal->GetPointsOverThreshold().size() >= (unsigned int)fPointsOverThreshold) {
                    // signal is not noise
                    VetoMaxPeakAmplitude_map[i] = signal->GetMaxPeakValue();
                } else {
                    // signal is noise
                    VetoMaxPeakAmplitude_map[i] = 0;
                }
                VetoPeakTime_map[i] = signal->GetMaxPeakBin();
                // We remove the signal from the event
                fSignalEvent->RemoveSignalWithId(i);
 
                // check if signal is above threshold
                if (signal->GetMaxPeakValue() > fThreshold) {
                    VetoAboveThreshold = 1;
                    NVetoAboveThreshold += 1;
                }
                // check if signal is in time window
                if (signal->GetMaxPeakBin() > fTimeWindow[0] && signal->GetMaxPeakBin() < fTimeWindow[1]) {
                    VetoInTimeWindow = 1;
                    NVetoInTimeWindow += 1;
                }
            }
        }
 
        SetObservableValue("PeakTime", VetoPeakTime_map);
        SetObservableValue("MaxPeakAmplitude", VetoMaxPeakAmplitude_map);
        if (fThreshold != -1) {
            SetObservableValue("VetoAboveThreshold", VetoAboveThreshold);
            SetObservableValue("NvetoAboveThreshold", NVetoAboveThreshold);
        }
        if (fTimeWindow[0] != -1) {
            SetObservableValue("VetoInTimeWindow", VetoInTimeWindow);
            SetObservableValue("NVetoInTimeWindow", NVetoInTimeWindow);
        }
    }
 
    // ***************************************************************
    // if the veto ids are defined within the veto groups ************
    // ***************************************************************
 
    // create observable names for veto groups
    for (const auto& fVetoGroupName : fVetoGroupNames) {
        fPeakTime.push_back("PeakTime_" + fVetoGroupName);
        fPeakAmp.push_back("MaxPeakAmplitude_" + fVetoGroupName);
    }
 
    if (fVetoSignalId[0] == -1) {
        // iterate over veto groups
        for (unsigned int i = 0; i < fVetoGroupNames.size(); i++) {
            // iterate over vetoes in each group
            vector<double> groupIds = StringToElements(fVetoGroupIds[i], ",");
            for (double groupId : groupIds) {
                // Checks if channel (groupIds) participated in the event. If not, it is -1
                if (fSignalEvent->GetSignalIndex(groupId) != -1) {
                    // We extract the parameters from the veto signal
                    TRestRawSignal* signal = fSignalEvent->GetSignalById(groupId);
                    // Deal with noise
                    signal->CalculateBaseLine(fBaseLineRange.X(), fBaseLineRange.Y(), "ROBUST");
                    signal->InitializePointsOverThreshold(TVector2(fPointThreshold, fSignalThreshold),
                                                          fPointsOverThreshold);
                    // Save two maps with (veto panel ID, max amplitude) and (veto panel ID, peak time)
                    if (signal->GetPointsOverThreshold().size() >= (unsigned int)fPointsOverThreshold) {
                        // signal is not noise
                        VetoMaxPeakAmplitude_map[groupId] = signal->GetMaxPeakValue();
                    } else {
                        // signal is noise
                        VetoMaxPeakAmplitude_map[groupId] = 0;
                    }
                    VetoPeakTime_map[groupId] = signal->GetMaxPeakBin();
                    // We remove the signal from the event
                    fSignalEvent->RemoveSignalWithId(groupId);
 
                    // check if signal is above threshold
                    if (signal->GetMaxPeakValue() > fThreshold) {
                        VetoAboveThreshold = 1;
                        NVetoAboveThreshold += 1;
                    }
                    // check if signal is in time window
                    if (signal->GetMaxPeakBin() > fTimeWindow[0] &&
                        signal->GetMaxPeakBin() < fTimeWindow[1]) {
                        VetoInTimeWindow = 1;
                        NVetoInTimeWindow += 1;
                    }
                }
            }
            SetObservableValue(fPeakTime[i], VetoPeakTime_map);
            SetObservableValue(fPeakAmp[i], VetoMaxPeakAmplitude_map);
 
            VetoMaxPeakAmplitude_map.clear();
            VetoPeakTime_map.clear();
        }
 
        if (fThreshold != -1) {
            SetObservableValue("VetoAboveThreshold", VetoAboveThreshold);
            SetObservableValue("NvetoAboveThreshold", NVetoAboveThreshold);
        }
        if (fTimeWindow[0] != -1) {
            SetObservableValue("VetoInTimeWindow", VetoInTimeWindow);
            SetObservableValue("NVetoInTimeWindow", NVetoInTimeWindow);
        }
    }
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
        fSignalEvent->PrintEvent();
 
        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) GetChar();
    }
 
    return fSignalEvent;
}
 
Int_t TRestRawVetoAnalysisProcess::GetGroupIndex(const string& groupName) {
    auto it = find(fVetoGroupNames.begin(), fVetoGroupNames.end(), groupName);
    if (it != fVetoGroupNames.end()) {
        return it - fVetoGroupNames.begin();
    }
    return -1;
}
 
string TRestRawVetoAnalysisProcess::GetGroupIds(const string& groupName) {
    Int_t index = GetGroupIndex(groupName);
    if (index != -1) {
        return fVetoGroupIds[index];
    }
    return "-1";
}
 
void TRestRawVetoAnalysisProcess::InitFromConfigFile() {
    fBaseLineRange = StringTo2DVector(GetParameter("baseLineRange", "(5,55)"));
    fRange = StringTo2DVector(GetParameter("range", "(5,507)"));
    fThreshold = StringToInteger(GetParameter("threshold", "-1"));
    fTimeWindow = StringToElements(GetParameter("timeWindow", "-1,-1"), ",");
    if (fTimeWindow.size() != 2) {
        cout << "Error: timeWindow has to consist of two comma-separated values." << endl;
        GetChar();
    }
    std::vector<double> potpars = StringToElements(GetParameter("PointsOverThresholdPars", "1.5,1.5,4"), ",");
    fPointThreshold = potpars[0];
    fSignalThreshold = potpars[1];
    fPointsOverThreshold = (Int_t)potpars[2];
 
    // **************************************************************
    // ***** Vetoes are defined as a single list ********************
    // **************************************************************
 
    fVetoSignalId = StringToElements(GetParameter("vetoSignalId", "-1"), ",");
 
    // **************************************************************
    // ***** Vetoes are defined in groups ***************************
    // **************************************************************
 
    // Read all the info from the veto group definitions
 
    TiXmlElement* vetoDefinition = GetElement("vetoGroup");
 
    while (vetoDefinition != nullptr) {
        fVetoGroupNames.push_back(GetFieldValue("name", vetoDefinition));
        fVetoGroupIds.push_back(GetFieldValue("signalIDs", vetoDefinition));
        vetoDefinition = GetNextElement(vetoDefinition);
    }
 
    // Stop, in case signalIDs and groups are defined separately
    if (fVetoSignalId[0] != -1 && !fVetoGroupNames.empty()) {
        cout << "Error: veto groups and veto IDs defined separately!" << endl;
        GetChar();
    }
}
 
void TRestRawVetoAnalysisProcess::PrintMetadata() {
    BeginPrintProcess();
 
    if ((fVetoSignalId.empty() || fVetoSignalId[0] != -1) && fVetoGroupNames.empty()) {
        RESTMetadata
            << "Veto analysis process will apply to veto signals (to be determined during processing)"
            << RESTendl;
        return;
    }
 
    for (unsigned int i = 0; i < fVetoGroupNames.size(); i++) {
        RESTMetadata << "Veto group " << fVetoGroupNames[i] << " signal IDs: " << fVetoGroupIds[i]
                     << RESTendl;
    }
 
    if (fVetoSignalId[0] != -1) {
        for (double i : fVetoSignalId) {
            RESTMetadata << "Veto signal ID: " << i << RESTendl;
        }
    } else {
        RESTMetadata << " " << RESTendl;
        RESTMetadata << "All veto signal IDs: ";
        for (unsigned int i = 0; i < fVetoGroupIds.size() - 1; i++) {
            RESTMetadata << fVetoGroupIds[i] << ",";
        }
        RESTMetadata << fVetoGroupIds[fVetoGroupIds.size() - 1] << RESTendl;
    }
    if (fThreshold != -1) {
        RESTMetadata << "Veto threshold: " << fThreshold << RESTendl;
    }
    if (fTimeWindow[0] != -1) {
        RESTMetadata << "Peak time window: (" << fTimeWindow[0] << ", " << fTimeWindow[1] << ")" << RESTendl;
    }
    RESTMetadata << "Noise reduction: Points over Threshold parameters = (" << fPointThreshold << ", "
                 << fSignalThreshold << ", " << fPointsOverThreshold << ")" << RESTendl;
 
    EndPrintProcess();
}