TRestRawSignalFittingProcess.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 "TRestRawSignalFittingProcess.h"
 
using namespace std;
 
ClassImp(TRestRawSignalFittingProcess);
 
TRestRawSignalFittingProcess::TRestRawSignalFittingProcess() { Initialize(); }
 
TRestRawSignalFittingProcess::TRestRawSignalFittingProcess(const char* configFilename) {
    Initialize();
 
    if (LoadConfigFromFile(configFilename)) LoadDefaultConfig();
}
 
TRestRawSignalFittingProcess::~TRestRawSignalFittingProcess() {}
 
void TRestRawSignalFittingProcess::LoadDefaultConfig() { SetTitle("Default config"); }
 
void TRestRawSignalFittingProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fRawSignalEvent = nullptr;
}
 
void TRestRawSignalFittingProcess::LoadConfig(const string& configFilename, const string& name) {
    if (LoadConfigFromFile(configFilename, name)) LoadDefaultConfig();
}
 
void TRestRawSignalFittingProcess::InitProcess() {
    // fSignalFittingObservables = TRestEventProcess::ReadObservables();
}
 
TRestEvent* TRestRawSignalFittingProcess::ProcessEvent(TRestEvent* inputEvent) {
    // no need for verbose copy now
    fRawSignalEvent = (TRestRawSignalEvent*)inputEvent;
 
    RESTDebug << "TRestRawSignalFittingProcess::ProcessEvent. Event ID : " << fRawSignalEvent->GetID()
              << RESTendl;
 
    Double_t SigmaMean = 0;
    vector<Double_t> Sigma(fRawSignalEvent->GetNumberOfSignals());
    Double_t RatioSigmaMaxPeakMean = 0;
    vector<Double_t> RatioSigmaMaxPeak(fRawSignalEvent->GetNumberOfSignals());
    Double_t ChiSquareMean = 0;
    vector<Double_t> ChiSquare(fRawSignalEvent->GetNumberOfSignals());
 
    map<int, Double_t> baselineFit;
    map<int, Double_t> amplitudeFit;
    map<int, Double_t> shapingTimeFit;
    map<int, Double_t> peakPositionFit;
 
    baselineFit.clear();
    amplitudeFit.clear();
    shapingTimeFit.clear();
    peakPositionFit.clear();
 
    for (int s = 0; s < fRawSignalEvent->GetNumberOfSignals(); s++) {
        TRestRawSignal* singleSignal = fRawSignalEvent->GetSignal(s);
 
        int MaxPeakBin = singleSignal->GetMaxPeakBin();
 
        // ShaperSin function (AGET theoretic curve) times logistic function
        TF1* f = new TF1("fit",
                         "[0]+[1]*TMath::Exp(-3. * (x-[3])/[2]) * "
                         "(x-[3])/[2] * (x-[3])/[2] * (x-[3])/[2] * "
                         "sin((x-[3])/[2])/(1+TMath::Exp(-10000*(x-[3])))",
                         0, 511);
        f->SetParameters(0, 2000, 70, 80);
        // f->SetParameters(0, 0);  // Initial values adjusted from Desmos
        // f->SetParLimits(0, 150, 350);
        // f->SetParameters(1, 2000);
        // f->SetParLimits(1, 30, 90000);
        // f->SetParameters(2, 70);
        // f->SetParLimits(2, 10, 80);
        // f->SetParameters(3, 80);
        // f->SetParLimits(3, 150, 250);
        f->SetParNames("Baseline", "Amplitude", "ShapingTime", "PeakPosition");
 
        // Create histogram from signal
        Int_t nBins = singleSignal->GetNumberOfPoints();
        TH1D* h = new TH1D("histo", "Signal to histo", nBins, 0, nBins);
 
        for (int i = 0; i < nBins; i++) {
            h->Fill(i, singleSignal->GetRawData(i));
        }
 
        // Fit histogram with ShaperSin
        h->Fit(f, "NWW", "", 0, 511);  // Options: R->fit in range, N->No draw, Q->Quiet
 
        Double_t sigma = 0;
        for (int j = MaxPeakBin - 145; j < MaxPeakBin + 165; j++) {
            sigma += (singleSignal->GetRawData(j) - f->Eval(j)) * (singleSignal->GetRawData(j) - f->Eval(j));
        }
        Sigma[s] = TMath::Sqrt(sigma / (145 + 165));
        RatioSigmaMaxPeak[s] = Sigma[s] / singleSignal->GetRawData(MaxPeakBin);
        RatioSigmaMaxPeakMean += RatioSigmaMaxPeak[s];
        SigmaMean += Sigma[s];
        ChiSquare[s] = f->GetChisquare();
        ChiSquareMean += ChiSquare[s];
 
        baselineFit[singleSignal->GetID()] = f->GetParameter(0);
        amplitudeFit[singleSignal->GetID()] = f->GetParameter(1);
        shapingTimeFit[singleSignal->GetID()] = f->GetParameter(2);
        peakPositionFit[singleSignal->GetID()] = f->GetParameter(3);
 
        fShaping = f->GetParameter(2);
        fStartPosition = f->GetParameter(3);
        fBaseline = f->GetParameter(0);
        fAmplitude = f->GetParameter(1);
 
        h->Delete();
    }
 
    SetObservableValue("FitBaseline_map", baselineFit);
    SetObservableValue("FitAmplitude_map", amplitudeFit);
    SetObservableValue("FitShapingTime_map", shapingTimeFit);
    SetObservableValue("FitPeakPosition_map", peakPositionFit);
 
    SigmaMean = SigmaMean / (fRawSignalEvent->GetNumberOfSignals());
    SetObservableValue("FitSigmaMean", SigmaMean);
 
    Double_t sigmaMeanStdDev = 0;
    for (int k = 0; k < fRawSignalEvent->GetNumberOfSignals(); k++) {
        sigmaMeanStdDev += (Sigma[k] - SigmaMean) * (Sigma[k] - SigmaMean);
    }
    Double_t SigmaMeanStdDev = TMath::Sqrt(sigmaMeanStdDev / fRawSignalEvent->GetNumberOfSignals());
    SetObservableValue("FitSigmaStdDev", SigmaMeanStdDev);
 
    ChiSquareMean = ChiSquareMean / fRawSignalEvent->GetNumberOfSignals();
    SetObservableValue("FitChiSquareMean", ChiSquareMean);
 
    RatioSigmaMaxPeakMean = RatioSigmaMaxPeakMean / fRawSignalEvent->GetNumberOfSignals();
    SetObservableValue("FitRatioSigmaMaxPeakMean", RatioSigmaMaxPeakMean);
 
    RESTDebug << "SigmaMean: " << SigmaMean << RESTendl;
    RESTDebug << "SigmaMeanStdDev: " << SigmaMeanStdDev << RESTendl;
    RESTDebug << "ChiSquareMean: " << ChiSquareMean << RESTendl;
    RESTDebug << "RatioSigmaMaxPeakMean: " << RatioSigmaMaxPeakMean << RESTendl;
    for (int k = 0; k < fRawSignalEvent->GetNumberOfSignals(); k++) {
        RESTDebug << "Standard deviation of signal number " << k << ": " << Sigma[k] << RESTendl;
        RESTDebug << "Chi square of fit signal number " << k << ": " << ChiSquare[k] << RESTendl;
        RESTDebug << "Sandard deviation divided by amplitude of signal number " << k << ": "
                  << RatioSigmaMaxPeak[k] << RESTendl;
    }
 
    // This will affect the calculation of observables, but not the stored
    // TRestRawSignal data.
    // fRawSignalEvent->SetBaseLineRange(fBaseLineRange);
    // fRawSignalEvent->SetRange(fIntegralRange);
 
    /* Perhaps we want to identify the points over threshold where to apply the
  fitting?
     * Then, we might need to initialize points over threshold
     *
  for (int s = 0; s < fRawSignalEvent->GetNumberOfSignals(); s++) {
    TRestRawSignal* sgnl = fRawSignalEvent->GetSignal(s);
 
  TRestRawSignal
    sgnl->InitializePointsOverThreshold(TVector2(fPointThreshold,
  fSignalThreshold),
                                        fNPointsOverThreshold);
 
  }
    */
 
    // If cut condition matches the event will be not registered.
    if (ApplyCut()) return nullptr;
 
    return fRawSignalEvent;
}
 
void TRestRawSignalFittingProcess::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();
}