diff --git a/source/framework/analysis/inc/TRestDataSetCalibration.h b/source/framework/analysis/inc/TRestDataSetCalibration.h new file mode 100644 index 000000000..475cc562e --- /dev/null +++ b/source/framework/analysis/inc/TRestDataSetCalibration.h @@ -0,0 +1,89 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef REST_TRestDataSetCalibration +#define REST_TRestDataSetCalibration + +#include "TRestMetadata.h" + +/// This class is meant to perform the calibration of different runs +class TRestDataSetCalibration : public TRestMetadata { + private: + /// Name of the output file + std::string fOutputFileName = ""; + + /// Name of the dataSet inside the config file + std::string fDataSetName = ""; + + /// Vector containing expected energy peaks in keV must be sorted + std::vector fEnergyPeaks; + + /// Vector containing calibrated peaks in ADCs + std::vector fCalibPeaks; + + /// Vector containing calibrated sigma in ADCs + std::vector fCalibFWHM; + + /// Name of the calibration file to be used + std::string fCalibFile = ""; + + /// Calibration variable to be used + std::string fCalObservable = ""; + + /// Range to be calibrated + TVector2 fCalibRange; + + /// Number of bins used in the calibration + Int_t fNBins; + + /// Slope from the calibration fit + Double_t fSlope = 0; + + /// Intercept of the calibration fit + Double_t fIntercept = 0; + + void Initialize() override; + void InitFromConfigFile() override; + + public: + void PrintMetadata() override; + + void Calibrate(); + + inline void SetDataSetName(const std::string& dSName) { fDataSetName = dSName; } + inline void SetOutputFileName(const std::string& outName) { fOutputFileName = outName; } + inline void SetCalibrationFile(const std::string& calibFile) { fCalibFile = calibFile; } + + inline auto GetCalibPeaks() const { return fCalibPeaks; } + inline auto GetCalibFWHM() const { return fCalibFWHM; } + + inline double GetSlope() const { return fSlope; } + inline double GetIntercept() const { return fIntercept; } + inline std::string GetCalObservable() const { return fCalObservable; } + + TRestDataSetCalibration(); + TRestDataSetCalibration(const char* configFilename, std::string name = ""); + ~TRestDataSetCalibration(); + + ClassDefOverride(TRestDataSetCalibration, 1); +}; +#endif diff --git a/source/framework/analysis/inc/TRestDataSetOdds.h b/source/framework/analysis/inc/TRestDataSetOdds.h new file mode 100644 index 000000000..61210a3b9 --- /dev/null +++ b/source/framework/analysis/inc/TRestDataSetOdds.h @@ -0,0 +1,75 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef REST_TRestDataSetOdds +#define REST_TRestDataSetOdds + +#include "TH1F.h" +#include "TRestCut.h" +#include "TRestMetadata.h" + +/// This class is meant to compute the log odds for different datasets +class TRestDataSetOdds : public TRestMetadata { + private: + /// Name of the output file + std::string fOutputFileName = ""; + + /// Name of the dataSet inside the config file + std::string fDataSetName = ""; + + /// Vector containing different obserbable names + std::vector fObsName; + + /// Vector containing different obserbable ranges + std::vector fObsRange; + + /// Vector containing number of bins for the different observables + std::vector fObsNbins; + + /// Name of the odds file to be used to get the PDF + std::string fOddsFile = ""; + + /// Cuts over the dataset for PDF selection + TRestCut* fCut = nullptr; + + /// Map containing the PDF of the different observables + std::map fHistos; //! + + void Initialize() override; + void InitFromConfigFile() override; + + public: + void PrintMetadata() override; + + void ComputeLogOdds(); + + inline void SetDataSetName(const std::string& dSName) { fDataSetName = dSName; } + inline void SetOutputFileName(const std::string& outName) { fOutputFileName = outName; } + inline void SetOddsFile(const std::string& oddsFile) { fOddsFile = oddsFile; } + + TRestDataSetOdds(); + TRestDataSetOdds(const char* configFilename, std::string name = ""); + ~TRestDataSetOdds(); + + ClassDefOverride(TRestDataSetOdds, 1); +}; +#endif diff --git a/source/framework/analysis/src/TRestDataSetCalibration.cxx b/source/framework/analysis/src/TRestDataSetCalibration.cxx new file mode 100644 index 000000000..4f441c73f --- /dev/null +++ b/source/framework/analysis/src/TRestDataSetCalibration.cxx @@ -0,0 +1,293 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +///////////////////////////////////////////////////////////////////////// +/// TRestDataSetCalibration performs the calibration of a TRestDataSet, the +/// calibration is performed over different peaks provided in the config +/// file. The first peak provided in the config file should correspond +/// to the maximum in the spectra. The expected position of the rest +/// of the peaks are estimated with respect to the maximum. +/// +/// A summary of the basic parameters is described below: +/// * **calObservable**: Name of the observable to be calibrated +/// * **nBins**: Number of bins for the calibration spectra +/// * **calibRange**: Calibration range inside calObservable that +/// will be used for the spectrum +/// * **dataSetName**: Name of the dataSet to be calibrated +/// * **calibFile**: Name of the calibration file to be used for +/// the calibration parameters, if empty it performs the calibration +/// over the dataSet provided. +/// +/// The energy peaks to calibrate are given using peak metadata: +/// * **energy**: Energy of the peak(s) to be calibrated in keV, +/// note that the first peak should correspond to the maximum +/// +/// The following metadata members are used as output after the fit: +/// * **calibPeaks**: Position of the mean value of the peak after +/// the fit using the same order as the energy peaks +/// * **fCalibFWHM**: FWHM in percentage for the different peaks, +/// using the same order as the energy peaks. +/// +/// If an output file name is provided a new dataset will be generated +/// with the corresponding metadata and calibration constants, while +/// adding a new observable, calib_energy. In addition, the spectrum +/// and the fit results will be stored in the output file. +/// +/// +/// ### Examples +/// RML example to calibrate a dataset with a peak at 5.89 keV +/// \code +/// +/// +/// +/// +/// +/// +/// +/// \endcode +/// +/// Example to perform calibration over a calibration dataSet using restRoot: +/// \code +/// [0] TRestDataSetCalibration cal ("calibration.rml"); +/// [1] cal.SetDataSetName("myDataSet.root"); +/// [2] cal.SetOutputFileName("myCalibratedDataSet.root"); +/// [3] cal.Calibrate(); +/// \endcode +/// +/// Example to perform calibration over a background dataSet using restRoot: +/// \code +/// [0] TRestDataSetCalibration cal ("calibration.rml"); +/// [1] cal.SetDataSetName("myBackgroundDataSet.root"); +/// [2] cal.SetOutputFileName("myBackgroundCalibratedDataSet.root"); +/// [3] cal.SetCalibrationFile("myCalibratedDataSet.root"); +/// [4] cal.Calibrate(); +/// \endcode +/// +///---------------------------------------------------------------------- +/// +/// REST-for-Physics - Software for Rare Event Searches Toolkit +/// +/// History of developments: +/// +/// 2023-03: First implementation of TRestDataSetCalibration +/// JuanAn Garcia +/// +/// \class TRestDataSetCalibration +/// \author: JuanAn Garcia e-mail: juanangp@unizar.es +/// +///
+/// + +#include "TRestDataSetCalibration.h" + +#include "TRestDataSet.h" + +ClassImp(TRestDataSetCalibration); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestDataSetCalibration::TRestDataSetCalibration() { Initialize(); } + +///////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param configFilename A const char* that defines the RML filename. +/// \param name The name of the metadata section. It will be used to find the +/// corresponding TRestDataSetCalibration section inside the RML. +/// +TRestDataSetCalibration::TRestDataSetCalibration(const char* configFilename, std::string name) + : TRestMetadata(configFilename) { + LoadConfigFromFile(fConfigFileName, name); + Initialize(); + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata(); +} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestDataSetCalibration::~TRestDataSetCalibration() {} + +/////////////////////////////////////////////// +/// \brief Function to initialize input/output event members and define +/// the section name +/// +void TRestDataSetCalibration::Initialize() { SetSectionName(this->ClassName()); } + +/////////////////////////////////////////////// +/// \brief Function to initialize some variables from +/// configfile, in case that the variable is not found +/// in the rml it checks the input file pattern. +/// +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", ""); +} + +///////////////////////////////////////////// +/// \brief Performs the calibration of a given dataSet. If +/// no calibration file is provided, it performs the gaussian +/// fit to the different parameters provided in the config file. +/// Otherwise, it takes the calibration constants from the +/// provided calibration file. +/// +void TRestDataSetCalibration::Calibrate() { + PrintMetadata(); + + TRestDataSet dataSet; + dataSet.Import(fDataSetName); + + dataSet.PrintMetadata(); + + std::unique_ptr spectrum; + std::unique_ptr gr; + std::unique_ptr linearFit; + + if (fCalibFile.empty()) { + auto histo = dataSet.GetDataFrame().Histo1D( + {"spectrum", "spectrum", fNBins, fCalibRange.X(), fCalibRange.X()}, fCalObservable); + spectrum = std::unique_ptr(static_cast(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 gauss; + gr = std::unique_ptr(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(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(f->Get("grFit")); + linearFit = std::unique_ptr(f->Get("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.SetDataSet(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(); + } + } +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestDataSetCalibration +/// +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; +} diff --git a/source/framework/analysis/src/TRestDataSetOdds.cxx b/source/framework/analysis/src/TRestDataSetOdds.cxx new file mode 100644 index 000000000..2aa4d18c8 --- /dev/null +++ b/source/framework/analysis/src/TRestDataSetOdds.cxx @@ -0,0 +1,281 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +///////////////////////////////////////////////////////////////////////// +/// TRestDataSetOdds performs the log odds of the different observables given +/// in the config file and for a particular dataSet. To perform the log odds +/// first the probability density funcion (PDF) is obtained for a set of +/// observables in the desired range. Later on, the log odds is computed as +/// log(1. - odds) - log(odds) obtaining a number which is proportional to +/// how likely is an event with respect the desired distribution; lower the number, +/// more likely is the event to the input distribution. New observables are created in +/// the output dataSet odds_obserbable and the addition of all of them in odds_total. +/// If an input odds file is provided, the different PDFs are retrieved from the input +/// file. +/// +/// A summary of the basic parameters is described below: +/// * **dataSetName**: Name of the dataSet to be computed +/// * **oddsFile**: Name of the input odds file to be used for the definition of the +/// different PDFs. +/// +/// The different observables, range and nBins are added as follow: +/// \code +/// +/// \endcode +/// * **name**: Name of the observable PDF to be computed +/// * **range**: Range of the observable PDF to be used +/// * **nBins**: Number of bins for the observable PDF +/// +/// In addition a TRestCut is used as input for the generation of PDFs, check TRestCut +/// class for more info. +/// +/// ### Examples +/// Example of RML config file: +/// \code +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// \endcode +/// +/// Example to compute the the odds over a dataSet using restRoot: +/// \code +/// [0] TRestDataSetOdds odds ("odds.rml"); +/// [1] odds.SetDataSetName("myDataSet.root"); +/// [2] odds.SetOutputFileName("myComputedOddsDataSet.root"); +/// [3] odds.ComputeLogOdds(); +/// \endcode +/// +/// Example to compute the the odds over a dataSet with input odds file using restRoot: +/// \code +/// [0] TRestDataSetOdds odds ("odds.rml"); +/// [1] odds.SetDataSetName("myDataSet.root"); +/// [2] odds.SetOutputFileName("myComputedOddsDataSet.root"); +/// [3] odds.odds.SetOddsFile("myOddsDataSet.root"); +/// [4] odds.ComputeLogOdds(); +/// \endcode +/// +///---------------------------------------------------------------------- +/// +/// REST-for-Physics - Software for Rare Event Searches Toolkit +/// +/// History of developments: +/// +/// 2023-03: First implementation of TRestDataSetOdds +/// JuanAn Garcia +/// +/// \class TRestDataSetOdds +/// \author: JuanAn Garcia e-mail: juanangp@unizar.es +/// +///
+/// + +#include "TRestDataSetOdds.h" + +#include "TRestDataSet.h" + +ClassImp(TRestDataSetOdds); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestDataSetOdds::TRestDataSetOdds() { Initialize(); } + +///////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param configFilename A const char* that defines the RML filename. +/// \param name The name of the metadata section. It will be used to find the +/// corresponding TRestDataSetOdds section inside the RML. +/// +TRestDataSetOdds::TRestDataSetOdds(const char* configFilename, std::string name) + : TRestMetadata(configFilename) { + LoadConfigFromFile(fConfigFileName, name); + Initialize(); + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata(); +} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestDataSetOdds::~TRestDataSetOdds() {} + +/////////////////////////////////////////////// +/// \brief Function to initialize input/output event members and define +/// the section name +/// +void TRestDataSetOdds::Initialize() { SetSectionName(this->ClassName()); } + +/////////////////////////////////////////////// +/// \brief Function to initialize some variables from +/// configfile +/// +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"); +} + +///////////////////////////////////////////// +/// \brief This function computes the log odds for +/// a given dataSet. If no calibration odds file is +/// provided it computes the PDF for the given +/// observables. Otherwise, it takes the PDF from the +/// input file. This function generate different observables +/// odds_obsName and the addition of all of them for a further +/// processing, which is stored in odds_total observable. +/// +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(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](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.SetDataSet(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(); + } + } +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members of TRestDataSetOdds +/// +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; +}