src/TRestRawSignalShapingProcess.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 "TRestRawSignalShapingProcess.h"
 
using namespace std;
 
#include <TFile.h>
#include <TMath.h>
 
ClassImp(TRestRawSignalShapingProcess);
 
TRestRawSignalShapingProcess::TRestRawSignalShapingProcess() { Initialize(); }
 
TRestRawSignalShapingProcess::TRestRawSignalShapingProcess(const char* configFilename) {
    Initialize();
    LoadConfigFromFile(configFilename);
}
 
TRestRawSignalShapingProcess::~TRestRawSignalShapingProcess() { delete fOutputSignalEvent; }
 
void TRestRawSignalShapingProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fInputSignalEvent = nullptr;
    fOutputSignalEvent = new TRestRawSignalEvent();
}
 
void TRestRawSignalShapingProcess::LoadConfig(const string& configFilename, const string& name) {
    LoadConfigFromFile(configFilename, name);
}
 
void TRestRawSignalShapingProcess::InitProcess() {
    /*
     * NOT IMPLEMENTED. TODO To use a generic response from a
     * predefined TRestDetectorSignal
     *
     * For the moment we do only a gaussian shaping"
     * /
 
    responseSignal = new TRestRawSignal();
 
    if( fShapingType == "responseFile" )
    {
        TString fullPathName = (TString) getenv("REST_PATH") + "/data/signal/" +
    fResponseFilename; TFile *f = new TFile(fullPathName); responseSignal =
    (TRestRawSignal *) f->Get("signal Response"); f->Close();
    }
 
    if( GetVerboseLevel() >= REST_Debug )
    {
        CreateCanvas();
        fCanvas->Draw();
    }
 
    if( fShapingType == "gaus" )
    {
        responseSignal->InitGaussian( 100, 100, 30, 200 );
 
        if( GetVerboseLevel() >= REST_Debug )
        {
            responseSignal->GetGraph()->Draw();
            fCanvas->Update();
            GetChar();
        }
    }
    */
}
 
TRestEvent* TRestRawSignalShapingProcess::ProcessEvent(TRestEvent* inputEvent) {
    fInputSignalEvent = (TRestRawSignalEvent*)inputEvent;
 
    if (fInputSignalEvent->GetNumberOfSignals() <= 0) {
        return nullptr;
    }
 
    std::vector<double> response;
    Int_t Nr = 0;
 
    if (fShapingType == "gaus") {
        Int_t cBin = (Int_t)(fShapingTime * 3.5);
        Nr = 2 * cBin;
        Double_t sigma = fShapingTime;
 
        response.resize(Nr);
 
        for (int i = 0; i < Nr; i++) {
            response[i] = TMath::Exp(-0.5 * (i - cBin) * (i - cBin) / sigma / sigma);
            response[i] = response[i] / TMath::Sqrt(2 * M_PI) / sigma;
        }
    } else if (fShapingType == "exponential") {
        Nr = (Int_t)(5 * fShapingTime);
 
        response.resize(Nr);
 
        for (int i = 0; i < Nr; i++) {
            Double_t coeff = ((Double_t)i) / fShapingTime;
            response[i] = TMath::Exp(-coeff);
        }
    } else if (fShapingType == "shaper") {
        Nr = (Int_t)(5 * fShapingTime);
 
        response.resize(Nr);
 
        for (int i = 0; i < Nr; i++) {
            Double_t coeff = ((Double_t)i) / fShapingTime;
            response[i] = TMath::Exp(-3. * coeff) * coeff * coeff * coeff;
        }
    } else if (fShapingType == "shaperSin") {
        Nr = (Int_t)(5 * fShapingTime);
 
        response.resize(Nr);
 
        for (int i = 0; i < Nr; i++) {
            Double_t coeff = ((Double_t)i) / fShapingTime;
            response[i] = TMath::Exp(-3. * coeff) * coeff * coeff * coeff * sin(coeff);
        }
    } else {
        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Warning)
            cout << "REST WARNING. Shaping type : " << fShapingType << " is not defined!!" << endl;
        return nullptr;
    }
 
    // Making sure that response integral is 1, and applying the gain
    Double_t sum = 0;
    for (int n = 0; n < Nr; n++) sum += response[n];
    for (int n = 0; n < Nr; n++) response[n] = response[n] * fShapingGain / sum;
 
    for (int n = 0; n < fInputSignalEvent->GetNumberOfSignals(); n++) {
        TRestRawSignal shapingSignal = TRestRawSignal();
        TRestRawSignal inSignal = *fInputSignalEvent->GetSignal(n);
        Int_t nBins = inSignal.GetNumberOfPoints();
 
        vector<double> out(nBins);
        for (int m = 0; m < nBins; m++) out[m] = 0;
 
        for (int m = 0; m < nBins; m++) {
            if (inSignal.GetData(m) >= 0) {
                if (fShapingType == "gaus") {
                    for (int n = -Nr / 2; m + n < nBins && n < Nr / 2; n++)
                        if (m + n >= 0) out[m + n] += response[n + Nr / 2] * inSignal.GetData(m);
                } else
                    for (int n = 0; m + n < nBins && n < Nr; n++)
                        out[m + n] += response[n] * inSignal.GetData(m);
            }
        }
 
        for (int i = 0; i < nBins; i++) {
            shapingSignal.AddPoint((Short_t)round(out[i]));
        }
        shapingSignal.SetSignalID(inSignal.GetSignalID());
 
        fOutputSignalEvent->AddSignal(shapingSignal);
    }
 
    return fOutputSignalEvent;
}
 
void TRestRawSignalShapingProcess::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();
}