TRestDetectorSignalToRawSignalProcess.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 "TRestDetectorSignalToRawSignalProcess.h"
 
using namespace std;
 
ClassImp(TRestDetectorSignalToRawSignalProcess);
 
TRestDetectorSignalToRawSignalProcess::TRestDetectorSignalToRawSignalProcess() { Initialize(); }
 
TRestDetectorSignalToRawSignalProcess::TRestDetectorSignalToRawSignalProcess(const char* configFilename) {
    Initialize();
    LoadConfig(configFilename);
}
 
TRestDetectorSignalToRawSignalProcess::~TRestDetectorSignalToRawSignalProcess() {
    delete fOutputRawSignalEvent;
}
 
void TRestDetectorSignalToRawSignalProcess::LoadConfig(const string& configFilename, const string& name) {
    LoadConfigFromFile(configFilename, name);
}
 
void TRestDetectorSignalToRawSignalProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fInputSignalEvent = nullptr;
    fOutputRawSignalEvent = new TRestRawSignalEvent();
}
 
void TRestDetectorSignalToRawSignalProcess::InitProcess() {
    if (fTriggerMode == "fixed") {
        fTimeStart = fTriggerFixedStartTime - fTriggerDelay * fSampling;
        fTimeEnd = fTimeStart + fNPoints * fSampling;
    }
 
    const set<string> validTriggerModes = {"firstDeposit", "integralThreshold", "fixed"};
    if (validTriggerModes.count(fTriggerMode.Data()) == 0) {
        RESTError << "Trigger mode set to: '" << fTriggerMode
                  << "' which is not a valid trigger mode. Please use one of the following trigger modes: ";
        for (const auto& triggerMode : validTriggerModes) {
            RESTError << triggerMode << " ";
        }
        RESTError << RESTendl;
        exit(1);
    }
 
    if (IsLinearCalibration()) {
        const auto range = numeric_limits<Short_t>::max() - numeric_limits<Short_t>::min();
        fCalibrationGain = range * (fCalibrationRange.Y() - fCalibrationRange.X()) /
                           (fCalibrationEnergy.Y() - fCalibrationEnergy.X());
        fCalibrationOffset = range * (fCalibrationRange.X() - fCalibrationGain * fCalibrationEnergy.X()) +
                             numeric_limits<Short_t>::min();
    }
}
 
TRestEvent* TRestDetectorSignalToRawSignalProcess::ProcessEvent(TRestEvent* inputEvent) {
    fInputSignalEvent = (TRestDetectorSignalEvent*)inputEvent;
 
    if (fInputSignalEvent->GetNumberOfSignals() <= 0) {
        return nullptr;
    }
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
        fOutputRawSignalEvent->PrintEvent();
    }
 
    fOutputRawSignalEvent->SetID(fInputSignalEvent->GetID());
    fOutputRawSignalEvent->SetSubID(fInputSignalEvent->GetSubID());
    fOutputRawSignalEvent->SetTimeStamp(fInputSignalEvent->GetTimeStamp());
    fOutputRawSignalEvent->SetSubEventTag(fInputSignalEvent->GetSubEventTag());
 
    if (fTriggerMode == "firstDeposit") {
        fTimeStart = fInputSignalEvent->GetMinTime() - fTriggerDelay * fSampling;
        fTimeEnd = fInputSignalEvent->GetMinTime() + (fNPoints - fTriggerDelay) * fSampling;
    } else if (fTriggerMode == "integralThreshold") {
        bool thresholdReached = false;
        for (Double_t t = fInputSignalEvent->GetMinTime() - fNPoints * fSampling;
             t <= fInputSignalEvent->GetMaxTime() + fNPoints * fSampling; t = t + 0.5) {
            Double_t energy = fInputSignalEvent->GetIntegralWithTime(t, t + (fSampling * fNPoints) / 2.);
 
            if (energy > fIntegralThreshold) {
                fTimeStart = t - fTriggerDelay * fSampling;
                fTimeEnd = t + (fNPoints - fTriggerDelay) * fSampling;
                thresholdReached = true;
            }
        }
        if (!thresholdReached) {
            RESTWarning << "Integral threshold for trigger not reached" << RESTendl;
            fTimeStart = 0;
            fTimeEnd = fNPoints * fSampling;
        }
    }
 
    RESTDebug << "fTimeStart : " << fTimeStart << " us " << RESTendl;
    RESTDebug << "fTimeEnd : " << fTimeEnd << " us " << RESTendl;
 
    for (int n = 0; n < fInputSignalEvent->GetNumberOfSignals(); n++) {
        vector<Double_t> data(fNPoints, fCalibrationOffset);
        TRestDetectorSignal* signal = fInputSignalEvent->GetSignal(n);
        Int_t signalID = signal->GetSignalID();
 
        for (int m = 0; m < signal->GetNumberOfPoints(); m++) {
            Double_t t = signal->GetTime(m);
            Double_t d = signal->GetData(m);
 
            if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug && n < 3 && m < 5) {
                cout << "Signal : " << n << " Sample : " << m << " T : " << t << " Data : " << d << endl;
            }
 
            if (t > fTimeStart && t < fTimeEnd) {
                // convert physical time (in us) to timeBin
                Int_t timeBin = (Int_t)round((t - fTimeStart) / fSampling);
 
                if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Warning) {
                    if (timeBin < 0 || timeBin > fNPoints) {
                        cout << "Time bin out of range!!! bin value : " << timeBin << endl;
                        timeBin = 0;
                    }
                }
 
                RESTDebug << "Adding data : " << signal->GetData(m) << " to Time Bin : " << timeBin
                          << RESTendl;
                data[timeBin] += fCalibrationGain * signal->GetData(m);
            }
        }
 
        if (IsShapingEnabled()) {
            const auto shapingFunction = [](Double_t t) -> Double_t {
                if (t <= 0) {
                    return 0;
                }
                // function is normalized such that its absolute maximum is 1.0
                // max is at x = 1.1664004483744728
                return TMath::Exp(-3.0 * t) * TMath::Power(t, 3.0) * TMath::Sin(t) * 22.68112123672292;
            };
 
            vector<Double_t> dataAfterShaping(fNPoints, fCalibrationOffset);
            for (int i = 0; i < fNPoints; i++) {
                const Double_t value = data[i] - fCalibrationOffset;
                if (value <= 0) {
                    // Only positive values are possible, 0 means no signal in this bin
                    continue;
                }
                for (int j = 0; j < fNPoints; j++) {
                    dataAfterShaping[j] += value * shapingFunction(((j - i) * fSampling) / fShapingTime);
                }
            }
            data = dataAfterShaping;
        }
 
        TRestRawSignal rawSignal;
        rawSignal.SetSignalID(signalID);
        for (int x = 0; x < fNPoints; x++) {
            double value = round(data[x]);
            if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Warning) {
                if (value < numeric_limits<Short_t>::min() || value > numeric_limits<Short_t>::max()) {
                    RESTWarning << "value (" << value << ") is outside short range ("
                                << numeric_limits<Short_t>::min() << ", " << numeric_limits<Short_t>::max()
                                << ")" << RESTendl;
                }
            }
 
            if (value < numeric_limits<Short_t>::min()) {
                value = numeric_limits<Short_t>::min();
                fOutputRawSignalEvent->SetOK(false);
            } else if (value > numeric_limits<Short_t>::max()) {
                value = numeric_limits<Short_t>::max();
                fOutputRawSignalEvent->SetOK(false);
            }
            rawSignal.AddPoint((Short_t)value);
        }
 
        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
            rawSignal.Print();
        }
        RESTDebug << "Adding signal to raw signal event" << RESTendl;
 
        fOutputRawSignalEvent->AddSignal(rawSignal);
    }
 
    RESTDebug << "TRestDetectorSignalToRawSignalProcess. Returning event with N signals "
              << fOutputRawSignalEvent->GetNumberOfSignals() << RESTendl;
 
    return fOutputRawSignalEvent;
}