TRestDataSetCalibration.cxx

/*************************************************************************
 * This file is part of the REST software framework.                     *
 *                                                                       *
 * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza)                *
 * For more information see https://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 https://www.gnu.org/licenses/.                            *
 * For the list of contributors see $REST_PATH/CREDITS.                  *
 *************************************************************************/
 
 
#include "TRestDataSetCalibration.h"
 
#include "TRestDataSet.h"
 
ClassImp(TRestDataSetCalibration);
 
TRestDataSetCalibration::TRestDataSetCalibration() { Initialize(); }
 
TRestDataSetCalibration::TRestDataSetCalibration(const char* configFilename, std::string name)
    : TRestMetadata(configFilename) {
    LoadConfigFromFile(fConfigFileName, name);
    Initialize();
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata();
}
 
TRestDataSetCalibration::~TRestDataSetCalibration() {}
 
void TRestDataSetCalibration::Initialize() { SetSectionName(this->ClassName()); }
 
void TRestDataSetCalibration::InitFromConfigFile() {
    Initialize();
    TRestMetadata::InitFromConfigFile();
 
    TiXmlElement* peakDefinition = GetElement("peak");
    while (peakDefinition != nullptr) {
        std::string energy = GetFieldValue("energy", peakDefinition);
        if (energy.empty() || energy == "Not defined") {
            RESTError << "< peak variable key does not contain energy!" << RESTendl;
            exit(1);
        } else {
            fEnergyPeaks.push_back(StringToDouble(energy));
        }
 
        peakDefinition = GetNextElement(peakDefinition);
    }
 
    if (fEnergyPeaks.empty()) {
        RESTError << "Energy peaks not provided, exiting..." << RESTendl;
        exit(1);
    }
 
    if (fOutputFileName == "") fOutputFileName = GetParameter("outputFileName", "");
}
 
void TRestDataSetCalibration::Calibrate() {
    PrintMetadata();
 
    TRestDataSet dataSet;
    dataSet.Import(fDataSetName);
 
    dataSet.PrintMetadata();
 
    std::unique_ptr<TH1F> spectrum;
    std::unique_ptr<TGraph> gr;
    std::unique_ptr<TF1> linearFit;
 
    if (fCalibFile.empty()) {
        auto histo = dataSet.GetDataFrame().Histo1D(
            {"spectrum", "spectrum", fNBins, fCalibRange.X(), fCalibRange.Y()}, fCalObservable);
        spectrum = std::unique_ptr<TH1F>(static_cast<TH1F*>(histo->DrawClone()));
 
        // Get position of the maximum
        const double max = spectrum->GetBinCenter(spectrum->GetMaximumBin());
        // Get ratio between maximum and energy peak
        const double ratio = max / fEnergyPeaks.front();
 
        RESTDebug << "Max peak position " << max << RESTendl;
 
        std::vector<TF1*> gauss;
        gr = std::unique_ptr<TGraph>(new TGraph());
        gr->SetName("grFit");
 
        int c = 0;
        for (const auto& energy : fEnergyPeaks) {
            std::string fName = "g" + std::to_string(c);
            double pos = energy * ratio;
            double start = pos * 0.85;
            double end = pos * 1.15;
            TF1* g = new TF1(fName.c_str(), "gaus", start, end);
            RESTDebug << "Fitting gaussian " << pos << " " << start << " " << end << RESTendl;
            spectrum->Fit(g, "R+");
            gauss.emplace_back(g);
            gr->SetPoint(c, g->GetParameter(1), energy);
            fCalibPeaks.push_back(g->GetParameter(1));
            fCalibFWHM.push_back(g->GetParameter(2) * 235. / g->GetParameter(1));
            c++;
        }
 
        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) GetChar();
 
        gr->SetPoint(c, 0, 0);
 
        gr->Draw("AP");
        linearFit = std::unique_ptr<TF1>(new TF1("linearFit", "pol1"));
        gr->Fit("linearFit", "S");
    } else {
        TFile* f = TFile::Open(fCalibFile.c_str());
        if (f == nullptr) {
            RESTError << "Cannot open calibration file " << fCalibFile << RESTendl;
            exit(1);
        }
        RESTInfo << "Opening " << fCalibFile << RESTendl;
        gr = std::unique_ptr<TGraph>(f->Get<TGraph>("grFit"));
        linearFit = std::unique_ptr<TF1>(f->Get<TF1>("linearFit"));
    }
 
    fSlope = linearFit->GetParameter(1);
    fIntercept = linearFit->GetParameter(0);
 
    RESTDebug << "Slope " << fSlope << " Intercept " << fIntercept << RESTendl;
 
    auto df = dataSet.GetDataFrame();
 
    auto calibrate = [&linearFit](double val) {
        return linearFit->GetParameter(1) * val + linearFit->GetParameter(0);
    };
    df = df.Define("calib_Energy", calibrate, {fCalObservable});
 
    dataSet.SetDataFrame(df);
 
    if (!fOutputFileName.empty()) {
        if (TRestTools::GetFileNameExtension(fOutputFileName) == "root") {
            dataSet.Export(fOutputFileName);
            TFile* f = TFile::Open(fOutputFileName.c_str(), "UPDATE");
            this->Write();
            if (gr) gr->Write();
            if (linearFit) linearFit->Write();
            // if(lFit)lFit->Write();
            // spectrumFit->Write();
            if (spectrum) spectrum->Write();
            f->Close();
        }
    }
}
 
void TRestDataSetCalibration::PrintMetadata() {
    TRestMetadata::PrintMetadata();
 
    RESTMetadata << " Energy peaks to calibrate: ";
    for (const auto& peak : fEnergyPeaks) RESTMetadata << peak << " keV; ";
    RESTMetadata << RESTendl;
    RESTMetadata << " Calibrated peak position: ";
    for (const auto& peak : fCalibPeaks) RESTMetadata << peak << " ADC; ";
    RESTMetadata << RESTendl;
    RESTMetadata << " Calibrated FWHM: ";
    for (const auto& fwhm : fCalibFWHM) RESTMetadata << fwhm << " %; ";
    RESTMetadata << RESTendl;
    RESTMetadata << "Observable to calibrate: " << fCalObservable << RESTendl;
    RESTMetadata << "Calibration range: (" << fCalibRange.X() << ", " << fCalibRange.Y() << ")" << RESTendl;
    RESTMetadata << "Number of bins: " << fNBins << RESTendl;
    RESTMetadata << "Slope: " << fSlope << " keV/ADC" << RESTendl;
    RESTMetadata << "Intercept: " << fIntercept << " keV" << RESTendl;
    RESTMetadata << "----" << RESTendl;
}