TRestRawSignalAnalysisProcess.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 "TRestRawSignalAnalysisProcess.h"
 
using namespace std;
 
ClassImp(TRestRawSignalAnalysisProcess);
 
TRestRawSignalAnalysisProcess::TRestRawSignalAnalysisProcess() { Initialize(); }
 
TRestRawSignalAnalysisProcess::~TRestRawSignalAnalysisProcess() {}
 
void TRestRawSignalAnalysisProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fInputEvent = nullptr;
}
 
void TRestRawSignalAnalysisProcess::InitProcess() {
    if (fSignalsRange.X() != -1 && fSignalsRange.Y() != -1) {
        fRangeEnabled = true;
    }
}
 
void TRestRawSignalAnalysisProcess::InitFromConfigFile() {
    TRestEventProcess::InitFromConfigFile();
    const auto filterType = GetParameter("channelType", "");
    if (!filterType.empty()) {
        fChannelTypes.insert(filterType);
    }
}
 
TRestEvent* TRestRawSignalAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
    fInputEvent = dynamic_cast<TRestRawSignalEvent*>(inputEvent);
 
    const auto run = GetRunInfo();
    if (run != nullptr) {
        fInputEvent->InitializeReferences(run);
    }
 
    if (fReadoutMetadata == nullptr) {
        fReadoutMetadata = fInputEvent->GetReadoutMetadata();
    }
 
    if (fReadoutMetadata == nullptr && !fChannelTypes.empty()) {
        cerr << "TRestRawSignalAnalysisProcess::ProcessEvent: readout metadata is null, cannot filter "
                "the process by signal type"
             << endl;
        exit(1);
    }
 
    auto event = fInputEvent->GetSignalEventForTypes(fChannelTypes, fReadoutMetadata);
 
    // we save some complex typed analysis result
    map<int, Double_t> baseline;
    map<int, Double_t> baselinesigma;
    map<int, Double_t> ampsgn_maxmethod;
    map<int, Double_t> ampsgn_intmethod;
    map<int, int> risetime;
    map<int, int> peak_time;
    map<int, int> npointsot;
    vector<int> saturatedchnId;
 
    baseline.clear();
    baselinesigma.clear();
    ampsgn_maxmethod.clear();
    ampsgn_intmethod.clear();
    risetime.clear();
    npointsot.clear();
 
    Int_t nGoodSignals = 0;
 
    // This will affect the calculation of observables, but not the stored
    // TRestRawSignal data.
    event.SetBaseLineRange(fBaseLineRange, fBaseLineOption);
    event.SetRange(fIntegralRange);
 
    for (int s = 0; s < event.GetNumberOfSignals(); s++) {
        TRestRawSignal* sgnl = event.GetSignal(s);
 
        sgnl->InitializePointsOverThreshold(TVector2(fPointThreshold, fSignalThreshold),
                                            fPointsOverThreshold);
 
        if (fRangeEnabled && (sgnl->GetID() < fSignalsRange.X() || sgnl->GetID() > fSignalsRange.Y())) {
            continue;
        }
 
        // We do not want that signals that are not identified as such contribute to
        // define our observables
        // nkx: we still need to store all the signals in baseline/rise time maps in
        // case for noise analysis
        // if (sgnl->GetPointsOverThreshold().size() < 2) continue;
        if (sgnl->GetPointsOverThreshold().size() >= 2) nGoodSignals++;
 
        // Now TRestRawSignal returns directly baseline subtracted values
        baseline[sgnl->GetID()] = sgnl->GetBaseLine();
        baselinesigma[sgnl->GetID()] = sgnl->GetBaseLineSigma();
        ampsgn_intmethod[sgnl->GetID()] = sgnl->GetThresholdIntegral();
        ampsgn_maxmethod[sgnl->GetID()] = sgnl->GetMaxPeakValue();
        risetime[sgnl->GetID()] = sgnl->GetRiseTime();
        peak_time[sgnl->GetID()] = sgnl->GetMaxPeakBin();
        npointsot[sgnl->GetID()] = sgnl->GetPointsOverThreshold().size();
        if (sgnl->IsADCSaturation()) saturatedchnId.push_back(sgnl->GetID());
    }
 
    SetObservableValue("pointsoverthres_map", npointsot);
    SetObservableValue("risetime_map", risetime);
    SetObservableValue("peak_time_map", peak_time);
    SetObservableValue("baseline_map", baseline);
    SetObservableValue("baselinesigma_map", baselinesigma);
    SetObservableValue("max_amplitude_map", ampsgn_maxmethod);
    SetObservableValue("thr_integral_map", ampsgn_intmethod);
    SetObservableValue("SaturatedChannelID", saturatedchnId);
 
    Double_t baseLineMean = event.GetBaseLineAverage();
    SetObservableValue("BaseLineMean", baseLineMean);
 
    Double_t baseLineSigma = event.GetBaseLineSigmaAverage();
    SetObservableValue("BaseLineSigmaMean", baseLineSigma);
 
    Double_t timeDelay = event.GetMaxTime() - event.GetMinTime();
    SetObservableValue("TimeBinsLength", timeDelay);
 
    Int_t nSignals = event.GetNumberOfSignals();
    SetObservableValue("NumberOfSignals", nSignals);
    SetObservableValue("NumberOfGoodSignals", nGoodSignals);
 
    // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    // SubstractBaselines
    // After this the signal gets zero-ed, for the following analysis
    // Keep in mind, to add raw signal analysis, we must write code at before
    // This is where most of the problems occur
    // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    // Javier: I believe we should not substract baseline in the analysis process
    // then ...
    // ... of course we need to consider baseline substraction for each
    // observable. TRestRawSignal methods
    // should do that internally. I have updated that to be like that, but we need
    // to be with open eyes for
    // some period.
    // Baseline substraction will always happen when we transfer a TRestRawSignal
    // to TRestDetectorSignal
    //
    // We do not substract baselines then now, as it was done before
    //
    // fInputEvent->SubstractBaselines(fBaseLineRange.X(), fBaseLineRange.Y());
    //
    // Methods in TRestRawSignal have been updated to consider baseline.
    // TRestRawSignal now implements that internally. We need to define the
    // baseline range, and the range
    // where calculations take place. All we need is to call at some point to the
    // following methods.
    //
    // TRestRawSignalEvent::SetBaseLineRange and TRestRawSignalEvent::SetRange.
    //
    // Then, if any method accepts a different range it will be given in the
    // method name,
    // for example: GetIntegralInRange( Int_t startBin, Int_t endBin );
    //
 
    Double_t fullIntegral = event.GetIntegral();
    SetObservableValue("FullIntegral", fullIntegral);
 
    Double_t thrIntegral = event.GetThresholdIntegral();
    SetObservableValue("ThresholdIntegral", thrIntegral);
 
    Double_t riseSlope = event.GetRiseSlope();
    SetObservableValue("RiseSlopeAvg", riseSlope);
 
    Double_t slopeIntegral = event.GetSlopeIntegral();
    SetObservableValue("SlopeIntegral", slopeIntegral);
 
    Double_t rateOfChange = riseSlope / slopeIntegral;
    if (slopeIntegral == 0) rateOfChange = 0;
    SetObservableValue("RateOfChangeAvg", rateOfChange);
 
    Double_t riseTime = event.GetRiseTime();
    SetObservableValue("RiseTimeAvg", riseTime);
 
    Double_t tripleMaxIntegral = event.GetTripleMaxIntegral();
    SetObservableValue("TripleMaxIntegral", tripleMaxIntegral);
 
    Double_t integralRatio = (fullIntegral - thrIntegral) / (fullIntegral + thrIntegral);
    SetObservableValue("IntegralBalance", integralRatio);
 
    Double_t maxValue = 0;
    Double_t minValue = 1.e6;
    Double_t maxValueIntegral = 0;
    Double_t minDownValue = 1.e6;
 
    Double_t minPeakTime = 1000;  // TODO substitute this for something better
    Double_t maxPeakTime = 0;
    Double_t peakTimeAverage = 0;
 
    for (int s = 0; s < event.GetNumberOfSignals(); s++) {
        TRestRawSignal* sgnl = event.GetSignal(s);
 
        if (fRangeEnabled && (sgnl->GetID() < fSignalsRange.X() || sgnl->GetID() > fSignalsRange.Y()))
            continue;
 
        if (sgnl->GetPointsOverThreshold().size() > 1) {
            Double_t value = event.GetSignal(s)->GetMaxValue();
            maxValueIntegral += value;
 
            if (value > maxValue) maxValue = value;
            if (value < minValue) minValue = value;
 
            Double_t peakBin = sgnl->GetMaxPeakBin();
            peakTimeAverage += peakBin;
 
            if (minPeakTime > peakBin) minPeakTime = peakBin;
            if (maxPeakTime < peakBin) maxPeakTime = peakBin;
        }
        Double_t mindownvalue = event.GetSignal(s)->GetMinValue();
        if (mindownvalue < minDownValue) {
            minDownValue = mindownvalue;
        }
    }
 
    if (nGoodSignals > 0) peakTimeAverage /= nGoodSignals;
 
    Double_t ampIntRatio = thrIntegral / maxValueIntegral;
    if (maxValueIntegral == 0) ampIntRatio = 0;
 
    SetObservableValue("AmplitudeIntegralRatio", ampIntRatio);
    SetObservableValue("MinPeakAmplitude", minValue);
    SetObservableValue("MaxPeakAmplitude", maxValue);
    SetObservableValue("PeakAmplitudeIntegral", maxValueIntegral);
 
    SetObservableValue("MinEventValue", minDownValue);
 
    Double_t amplitudeRatio = maxValueIntegral / maxValue;
    if (maxValue == 0) amplitudeRatio = 0;
 
    SetObservableValue("AmplitudeRatio", amplitudeRatio);
    SetObservableValue("MaxPeakTime", maxPeakTime);
    SetObservableValue("MinPeakTime", minPeakTime);
 
    Double_t peakTimeDelay = maxPeakTime - minPeakTime;
 
    SetObservableValue("MaxPeakTimeDelay", peakTimeDelay);
    SetObservableValue("AveragePeakTime", peakTimeAverage);
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
        for (const auto& i : fObservablesDefined) {
            fAnalysisTree->PrintObservable(i.second);
        }
    }
 
    // If cut condition matches the event will be not registered.
    if (ApplyCut()) {
        return nullptr;
    }
 
    return fInputEvent;
}