TRestDataSetOdds.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 "TRestDataSetOdds.h"
 
#include "TRestDataSet.h"
 
ClassImp(TRestDataSetOdds);
 
TRestDataSetOdds::TRestDataSetOdds() { Initialize(); }
 
TRestDataSetOdds::TRestDataSetOdds(const char* configFilename, std::string name)
    : TRestMetadata(configFilename) {
    LoadConfigFromFile(fConfigFileName, name);
    Initialize();
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata();
}
 
TRestDataSetOdds::~TRestDataSetOdds() {}
 
void TRestDataSetOdds::Initialize() { SetSectionName(this->ClassName()); }
 
void TRestDataSetOdds::InitFromConfigFile() {
    Initialize();
    TRestMetadata::InitFromConfigFile();
 
    TiXmlElement* obsDefinition = GetElement("observable");
    while (obsDefinition != nullptr) {
        std::string obsName = GetFieldValue("name", obsDefinition);
        if (obsName.empty() || obsName == "Not defined") {
            RESTError << "< observable variable key does not contain a name!" << RESTendl;
            exit(1);
        } else {
            fObsName.push_back(obsName);
        }
 
        std::string range = GetFieldValue("range", obsDefinition);
        if (range.empty() || range == "Not defined") {
            RESTError << "< observable key does not contain a range value!" << RESTendl;
            exit(1);
        } else {
            TVector2 roi = StringTo2DVector(range);
            fObsRange.push_back(roi);
        }
 
        std::string nBins = GetFieldValue("nBins", obsDefinition);
        if (nBins.empty() || nBins == "Not defined") {
            RESTError << "< observable key does not contain a nBins value!" << RESTendl;
            exit(1);
        } else {
            fObsNbins.push_back(StringToInteger(nBins));
        }
 
        obsDefinition = GetNextElement(obsDefinition);
    }
 
    if (fObsName.empty() || fObsRange.empty()) {
        RESTError << "No observables provided, exiting..." << RESTendl;
        exit(1);
    }
 
    if (fOutputFileName == "") fOutputFileName = GetParameter("outputFileName", "");
 
    fCut = (TRestCut*)InstantiateChildMetadata("TRestCut");
}
 
void TRestDataSetOdds::ComputeLogOdds() {
    PrintMetadata();
 
    TRestDataSet dataSet;
    dataSet.Import(fDataSetName);
 
    if (fOddsFile.empty()) {
        auto DF = dataSet.MakeCut(fCut);
        RESTInfo << "Generating PDFs for dataset: " << fDataSetName << RESTendl;
        for (size_t i = 0; i < fObsName.size(); i++) {
            const std::string obsName = fObsName[i];
            const TVector2 range = fObsRange[i];
            const std::string histName = "h" + std::string(GetName()) + "_" + obsName;
            const int nBins = fObsNbins[i];
            RESTDebug << "\tGenerating PDF for " << obsName << " with range: (" << range.X() << ", "
                      << range.Y() << ") and nBins: " << nBins << RESTendl;
            auto histo =
                DF.Histo1D({histName.c_str(), histName.c_str(), nBins, range.X(), range.Y()}, obsName);
            TH1F* h = static_cast<TH1F*>(histo->DrawClone());
            RESTDebug << "\tNormalizing by integral = " << h->Integral() << RESTendl;
            h->Scale(1. / h->Integral());
            fHistos[obsName] = h;
        }
    } else {
        TFile* f = TFile::Open(fOddsFile.c_str());
        if (f == nullptr) {
            RESTError << "Cannot open calibration odds file " << fOddsFile << RESTendl;
            exit(1);
        }
        RESTInfo << "Opening " << fOddsFile << " as oddsFile." << RESTendl;
        for (size_t i = 0; i < fObsName.size(); i++) {
            const std::string obsName = fObsName[i];
            const std::string histName = "h" + std::string(GetName()) + "_" + obsName;
            TH1F* h = (TH1F*)f->Get(histName.c_str());
            fHistos[obsName] = h;
        }
    }
 
    auto df = dataSet.GetDataFrame();
    std::string totName = "";
    RESTDebug << "Computing log odds from " << fDataSetName << RESTendl;
    for (const auto& [obsName, histo] : fHistos) {
        const std::string oddsName = std::string(GetName()) + "_" + obsName;
        auto GetLogOdds = [&histo = histo](double val) {
            double odds = histo->GetBinContent(histo->GetXaxis()->FindBin(val));
            if (odds == 0) return 1000.;
            return log(1. - odds) - log(odds);
        };
 
        if (df.GetColumnType(obsName) != "Double_t") {
            RESTWarning << "Column " << obsName << " is not of type 'double'. It will be converted."
                        << RESTendl;
            df = df.Redefine(obsName, "static_cast<double>(" + obsName + ")");
        }
        df = df.Define(oddsName, GetLogOdds, {obsName});
        auto h = df.Histo1D(oddsName);
 
        if (!totName.empty()) totName += "+";
        totName += oddsName;
    }
 
    RESTDebug << "Computing total log odds" << RESTendl;
    RESTDebug << "\tTotal log odds = " << totName << RESTendl;
    df = df.Define(std::string(GetName()) + "_total", totName);
 
    dataSet.SetDataFrame(df);
 
    if (!fOutputFileName.empty()) {
        if (TRestTools::GetFileNameExtension(fOutputFileName) == "root") {
            RESTDebug << "Exporting dataset to " << fOutputFileName << RESTendl;
            dataSet.Export(fOutputFileName);
            TFile* f = TFile::Open(fOutputFileName.c_str(), "UPDATE");
            this->Write();
            RESTDebug << "Writing histograms to " << fOutputFileName << RESTendl;
            for (const auto& [obsName, histo] : fHistos) histo->Write();
            f->Close();
        }
    }
}
 
std::vector<std::tuple<std::string, TVector2, int>> TRestDataSetOdds::GetOddsObservables() {
    std::vector<std::tuple<std::string, TVector2, int>> obs;
    for (size_t i = 0; i < fObsName.size(); i++) {
        if (i >= fObsName.size() || i >= fObsRange.size() || i >= fObsNbins.size()) {
            RESTError << "Sizes for observables names, ranges and bins do not match!" << RESTendl;
            break;
        }
        obs.push_back(std::make_tuple(fObsName[i], fObsRange[i], fObsNbins[i]));
    }
    return obs;
}
 
void TRestDataSetOdds::AddOddsObservable(const std::string& name, const TVector2& range, int nbins) {
    fObsName.push_back(name);
    fObsRange.push_back(range);
    fObsNbins.push_back(nbins);
}
 
void TRestDataSetOdds::SetOddsObservables(const std::vector<std::tuple<std::string, TVector2, int>>& obs) {
    fObsName.clear();
    fObsRange.clear();
    fObsNbins.clear();
    for (const auto& [name, range, nbins] : obs) AddOddsObservable(name, range, nbins);
}
 
void TRestDataSetOdds::WriteHistograms(TFile* f) const {
    if (!f) return;
    f->cd();
    for (const auto& [name, histo] : fHistos) {
        if (histo) histo->Write();
    }
}
 
void TRestDataSetOdds::PrintMetadata() {
    TRestMetadata::PrintMetadata();
 
    // if (fCut) fCut->PrintMetadata();
    if (!fOddsFile.empty()) RESTMetadata << " Odds file: " << fOddsFile << RESTendl;
    RESTMetadata << " DataSet file: " << fDataSetName << RESTendl;
 
    RESTMetadata << " Observables to compute: " << RESTendl;
    for (size_t i = 0; i < fObsName.size(); i++) {
        RESTMetadata << fObsName[i] << "; Range: (" << fObsRange[i].X() << ", " << fObsRange[i].Y()
                     << "); nBins: " << fObsNbins[i] << RESTendl;
    }
    RESTMetadata << "----" << RESTendl;
}