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);
        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" + obsName;
            const int nBins = fObsNbins[i];
            auto histo =
                DF.Histo1D({histName.c_str(), histName.c_str(), nBins, range.X(), range.Y()}, obsName);
            TH1F* h = static_cast<TH1F*>(histo->DrawClone());
            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);
        }
        std::cout << "Opening " << fOddsFile << std::endl;
        for (size_t i = 0; i < fObsName.size(); i++) {
            const std::string obsName = fObsName[i];
            const std::string histName = "h" + obsName;
            TH1F* h = (TH1F*)f->Get(histName.c_str());
            fHistos[obsName] = h;
        }
    }
 
    auto df = dataSet.GetDataFrame();
    std::string totName = "";
    for (const auto& [obsName, histo] : fHistos) {
        const std::string oddsName = "odds_" + 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);
        };
        df = df.Define(oddsName, GetLogOdds, {obsName});
        auto h = df.Histo1D(oddsName);
 
        if (!totName.empty()) totName += "+";
        totName += oddsName;
    }
 
    df = df.Define("odds_total", totName);
 
    dataSet.SetDataFrame(df);
 
    if (!fOutputFileName.empty()) {
        if (TRestTools::GetFileNameExtension(fOutputFileName) == "root") {
            dataSet.Export(fOutputFileName);
            TFile* f = TFile::Open(fOutputFileName.c_str(), "UPDATE");
            this->Write();
            for (const auto& [obsName, histo] : fHistos) histo->Write();
            f->Close();
        }
    }
}
 
void TRestDataSetOdds::PrintMetadata() {
    TRestMetadata::PrintMetadata();
 
    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;
}