TRestDataSetPlot.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 "TRestDataSetPlot.h"
 
#include "TCanvas.h"
#include "TDirectory.h"
#include "TStyle.h"
 
ClassImp(TRestDataSetPlot);
 
TRestDataSetPlot::TRestDataSetPlot() { Initialize(); }
 
TRestDataSetPlot::TRestDataSetPlot(const char* configFilename, std::string name)
    : TRestMetadata(configFilename) {
    Initialize();
    LoadConfigFromFile(fConfigFileName, name);
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata();
}
 
TRestDataSetPlot::~TRestDataSetPlot() {}
 
void TRestDataSetPlot::Initialize() { SetSectionName(this->ClassName()); }
 
void TRestDataSetPlot::InitFromConfigFile() {
    TRestMetadata::InitFromConfigFile();
 
    if (fDataSetName == "") fDataSetName = GetParameter("inputFileName", "");
    if (fOutputFileName == "") fOutputFileName = GetParameter("outputFileName", "");
 
    fCut = ReadCut(fCut);
 
    ReadPlotInfo();
    ReadPanelInfo();
}
 
TRestCut* TRestDataSetPlot::ReadCut(TRestCut* cut, TiXmlElement* ele) {
    TiXmlElement* cutele = GetElement("addCut", ele);
    while (cutele != nullptr) {
        std::string cutName = GetParameter("name", cutele, "");
        if (!cutName.empty()) {
            if (cut == nullptr) {
                cut = (TRestCut*)InstantiateChildMetadata("TRestCut", cutName);
            } else {
                cut->AddCut((TRestCut*)InstantiateChildMetadata("TRestCut", cutName));
            }
        }
        cutele = GetNextElement(cutele);
    }
 
    return cut;
}
 
void TRestDataSetPlot::ReadPanelInfo() {
    if (!fPanels.empty()) {
        RESTWarning << "Plot metadata already initialized" << RESTendl;
    }
 
    TiXmlElement* panelele = GetElement("panel");
    while (panelele != nullptr) {
        std::string active = GetParameter("value", panelele, "ON");
        if (ToUpper(active) != "ON") continue;
 
        PanelInfo panel;
        panel.font_size = StringToDouble(GetParameter("font_size", panelele, "0.1"));
        panel.precision = StringToInteger(GetParameter("precision", panelele, "2"));
        panel.delimiter = GetParameter("delimiter", panelele, ": ");
 
        panel.panelCut = ReadCut(panel.panelCut, panelele);
 
        TiXmlElement* labelele = GetElement("variable", panelele);
        while (labelele != nullptr) {
            std::array<std::string, 3> label;
            label[0] = GetParameter("value", labelele, "");
            label[1] = GetParameter("label", labelele, "");
            label[2] = GetParameter("units", labelele, "");
            double posX = StringToDouble(GetParameter("x", labelele, "0.1"));
            double posY = StringToDouble(GetParameter("y", labelele, "0.1"));
 
            panel.variablePos.push_back(std::make_pair(label, TVector2(posX, posY)));
 
            labelele = GetNextElement(labelele);
        }
        TiXmlElement* metadata = GetElement("metadata", panelele);
        while (metadata != nullptr) {
            std::array<std::string, 3> label;
            label[0] = GetParameter("value", metadata, "");
            label[1] = GetParameter("label", metadata, "");
            label[2] = GetParameter("units", metadata, "");
            double posX = StringToDouble(GetParameter("x", metadata, "0.1"));
            double posY = StringToDouble(GetParameter("y", metadata, "0.1"));
 
            panel.metadataPos.push_back(std::make_pair(label, TVector2(posX, posY)));
 
            metadata = GetNextElement(metadata);
        }
        TiXmlElement* observable = GetElement("observable", panelele);
        while (observable != nullptr) {
            std::array<std::string, 3> label;
            label[0] = GetParameter("value", observable, "");
            label[1] = GetParameter("label", observable, "");
            label[2] = GetParameter("units", observable, "");
            double posX = StringToDouble(GetParameter("x", observable, "0.1"));
            double posY = StringToDouble(GetParameter("y", observable, "0.1"));
 
            panel.obsPos.push_back(std::make_pair(label, TVector2(posX, posY)));
 
            observable = GetNextElement(observable);
        }
        TiXmlElement* expression = GetElement("expression", panelele);
        while (expression != nullptr) {
            std::array<std::string, 4> label;
            label[0] = GetParameter("value", expression, "");
            label[1] = GetParameter("label", expression, "");
            label[2] = GetParameter("units", expression, "");
            label[3] = GetParameter("precision", expression, IntegerToString(panel.precision));
            double posX = StringToDouble(GetParameter("x", expression, "0.1"));
            double posY = StringToDouble(GetParameter("y", expression, "0.1"));
 
            panel.expPos.push_back(std::make_pair(label, TVector2(posX, posY)));
 
            expression = GetNextElement(expression);
        }
 
        fPanels.push_back(panel);
        panelele = GetNextElement(panelele);
    }
}
 
void TRestDataSetPlot::ReadPlotInfo() {
    if (!fPlots.empty()) {
        RESTWarning << "Plot metadata already initialized" << RESTendl;
        return;
    }
 
    TiXmlElement* plotele = GetElement("plot");
    while (plotele != nullptr) {
        std::string active = GetParameter("value", plotele, "ON");
        if (ToUpper(active) != "ON") continue;
        int N = fPlots.size();
        PlotInfo plot;
        plot.name = RemoveWhiteSpaces(GetParameter("name", plotele, "plot_" + ToString(N)));
        plot.title = GetParameter("title", plotele, plot.name);
        plot.logX = StringToBool(GetParameter("logX", plotele, "false"));
        plot.logY = StringToBool(GetParameter("logY", plotele, "false"));
        plot.logZ = StringToBool(GetParameter("logZ", plotele, "false"));
        plot.gridY = StringToBool(GetParameter("gridY", plotele, "false"));
        plot.gridX = StringToBool(GetParameter("gridX", plotele, "false"));
        plot.normalize = StringToDouble(GetParameter("norm", plotele, ""));
        plot.scale = GetParameter("scale", plotele, "");
        plot.labelX = GetParameter("xlabel", plotele, "");
        plot.labelY = GetParameter("ylabel", plotele, "");
        plot.marginBottom = StringToDouble(GetParameter("marginBottom", plotele, "0.15"));
        plot.marginTop = StringToDouble(GetParameter("marginTop", plotele, "0.07"));
        plot.marginLeft = StringToDouble(GetParameter("marginLeft", plotele, "0.25"));
        plot.marginRight = StringToDouble(GetParameter("marginRight", plotele, "0.1"));
        plot.legendOn = StringToBool(GetParameter("legend", plotele, "OFF"));
        plot.stackDrawOption = GetParameter("stackOption", plotele, "nostack");
        // plot.annotationOn = StringToBool(GetParameter("annotation", plotele, "OFF"));
        plot.xOffset = StringToDouble(GetParameter("xOffset", plotele, "0"));
        plot.yOffset = StringToDouble(GetParameter("yOffset", plotele, "0"));
        plot.timeDisplay = StringToBool(GetParameter("timeDisplay", plotele, "OFF"));
        plot.save = RemoveWhiteSpaces(GetParameter("save", plotele, ""));
 
        TiXmlElement* histele = GetElement("histo", plotele);
        if (histele == nullptr) {
            histele = plotele;
        }
        while (histele != nullptr) {
            HistoInfo hist;
            hist.name = RemoveWhiteSpaces(GetParameter("name", histele, plot.name));
            hist.drawOption = GetParameter("option", histele, "colz");
            TiXmlElement* varele = GetElement("variable", histele);
            while (varele != nullptr) {
                hist.variable.push_back(GetParameter("name", varele));
                std::string rangeStr = GetParameter("range", varele);
                hist.range.push_back(StringTo2DVector(rangeStr));
                hist.nBins.push_back(StringToInteger(GetParameter("nbins", varele)));
                varele = GetNextElement(varele);
            }
 
            hist.lineColor = GetIDFromMapString(ColorIdMap, GetParameter("lineColor", histele, "602"));
            hist.lineWidth = StringToInteger(GetParameter("lineWidth", histele, "1"));
            hist.lineStyle = GetIDFromMapString(LineStyleMap, GetParameter("lineStyle", histele, "1"));
            hist.fillStyle = GetIDFromMapString(FillStyleMap, GetParameter("fillStyle", histele, "1001"));
            hist.fillColor = GetIDFromMapString(ColorIdMap, GetParameter("fillColor", histele, "0"));
            hist.statistics = StringToBool(GetParameter("stats", histele, "OFF"));
            // hist.weight = GetParameter("weight", histele, "");
            hist.histoCut = ReadCut(hist.histoCut, histele);
            plot.histos.push_back(hist);
 
            if (histele == plotele) {
                break;
            }
            histele = GetNextElement(histele);
        }
 
        fPlots.push_back(plot);
        plotele = GetNextElement(plotele);
    }
}
 
void TRestDataSetPlot::GenerateDataSetFromFilePattern(TRestDataSet& dataSet) {
    std::vector<std::string> obsList;
 
    // Add obserbables from global cut
    if (fCut != nullptr) {
        const auto paramCut = fCut->GetParamCut();
        for (const auto& [param, condition] : paramCut) {
            obsList.push_back(param);
        }
    }
 
    // Add obserbables from plot info, both variables and cuts
    for (const auto& plots : fPlots) {
        for (const auto& hist : plots.histos) {
            for (const auto& var : hist.variable) {
                obsList.push_back(var);
            }
            if (hist.histoCut == nullptr) continue;
            const auto paramCut = hist.histoCut->GetParamCut();
            for (const auto& [param, condition] : paramCut) {
                obsList.push_back(param);
            }
        }
    }
 
    std::map<std::string, TRestDataSet::RelevantQuantity> quantity;
 
    for (auto& panel : fPanels) {
        // Add metadata and observables from expression key from panel info
        for (auto& [key, posLabel] : panel.expPos) {
            // look for metadata which are surrounded by [] but not [[]] (variables)
            auto&& [exp, label, units, precision] = key;
            std::string text = exp;
            while (text.find_last_of('[') != std::string::npos) {
                int squareBracketCorrector = 0;
                size_t posOpen = text.find_last_of('[');
                size_t posClose = text.find_first_of(']', posOpen);
                if (posOpen != 0) {
                    if (text[posOpen - 1] == '[') {
                        squareBracketCorrector = 1;
                    }
                }
                std::string varOrMeta = text.substr(posOpen - squareBracketCorrector,
                                                    posClose + 1 - posOpen + 2 * squareBracketCorrector);
                if (squareBracketCorrector == 0) {
                    // Here varOrMeta is a metadata
                    // Add metadata to relevant quantity from the panel info
                    TRestDataSet::RelevantQuantity quant;
                    quant.metadata = varOrMeta;
                    quant.strategy = "unique";
                    quantity[label] = quant;
                }
                text = Replace(text, varOrMeta, "1");
            }
 
            // look for observables (characterized by having a _ in the name)
            while (text.find("_") != std::string::npos) {
                size_t pos_ = text.find("_");
                size_t beginning = text.find_last_of(" -+*/)(^%", pos_) + 1;
                size_t end = text.find_first_of(" -+*/)(^%", pos_);
                std::string obs = text.substr(beginning, end - beginning);
                text = Replace(text, obs, "1");
                obsList.push_back(obs);
            }
 
            if (!(isAExpression(text) || isANumber(text)))
                RESTWarning << "The expression " << exp
                            << " has not been correctly parsed into variables, metadata and observables"
                            << RESTendl;
        }
 
        // Add observables from observable key of panel info
        for (auto& [key, posLabel] : panel.obsPos) {
            auto&& [obs, label, units] = key;
            obsList.push_back(obs);
        }
        // Add relevant quantity to metadata from the panel info
        for (auto& [key, posLabel] : panel.metadataPos) {
            auto&& [metadata, label, units] = key;
            TRestDataSet::RelevantQuantity quant;
            quant.metadata = metadata;
            quant.strategy = "unique";
            quantity[label] = quant;
        }
    }
 
    // Remove duplicated observables if any
    std::sort(obsList.begin(), obsList.end());
    obsList.erase(std::unique(obsList.begin(), obsList.end()), obsList.end());
    dataSet.SetFilePattern(fDataSetName);
    dataSet.SetObservablesList(obsList);
    dataSet.SetQuantity(quantity);
    dataSet.GenerateDataSet();
}
 
void TRestDataSetPlot::PlotCombinedCanvas() {
    TRestDataSet dataSet;
    dataSet.EnableMultiThreading(true);
 
    // Import dataSet
    dataSet.Import(fDataSetName);
 
    // If dataSet is not valid, try to generate it, but note that this is deprecated
    if (dataSet.GetTree() == nullptr) {
        RESTWarning << "Cannot import dataSet, trying to generate it with pattern " << fDataSetName
                    << RESTendl;
        RESTWarning << "Note that the generation of a dataSet inside TRestDataSetPlot is deplecated. Check "
                       "TRestDataSet documentation to generate a dataSet"
                    << RESTendl;
        GenerateDataSetFromFilePattern(dataSet);
        if (dataSet.GetTree() == nullptr) {
            RESTError << "Cannot generate dataSet " << RESTendl;
            exit(1);
        }
    }
 
    // Perform the global cut over the dataSet
    dataSet.SetDataFrame(dataSet.MakeCut(fCut));
 
    TCanvas combinedCanvas(this->GetName(), this->GetName(), 0, 0, fCanvasSize.X(), fCanvasSize.Y());
    combinedCanvas.Divide((Int_t)fCanvasDivisions.X(), (Int_t)fCanvasDivisions.Y(),
                          fCanvasDivisionMargins.X(), fCanvasDivisionMargins.Y());
 
    gStyle->SetPalette(fPaletteStyle);
 
    const double duration = dataSet.GetTotalTimeInSeconds();
    const double startTime = dataSet.GetStartTime();
    const double endTime = dataSet.GetEndTime();
 
    // Fill general variables for the panel
    std::map<std::string, std::string> paramMap;
    paramMap["[[startTime]]"] = ToDateTimeString(startTime);
    paramMap["[[endTime]]"] = ToDateTimeString(endTime);
 
    // DataSet quantity is used to replace metadata parameters
    const auto quantity = dataSet.GetQuantity();
 
    int canvasIndex = 1;
 
    for (auto& panel : fPanels) {
        combinedCanvas.cd(canvasIndex);
        // Gets a dataFrame with the panel cut
        auto dataFrame = dataSet.MakeCut(panel.panelCut);
        const int entries = *dataFrame.Count();
        const double meanRate = entries / duration;
        const double runLength = duration / 3600.;
        paramMap["[[runLength]]"] = DoubleToString(runLength, "%.9e");
        paramMap["[[entries]]"] = IntegerToString(entries);
        paramMap["[[meanRate]]"] = DoubleToString(meanRate, "%.9e");
 
        paramMap["[[cutNames]]"] = "";
        paramMap["[[cuts]]"] = "";
        if (fCut) {
            for (const auto& cut : fCut->GetCuts()) {
                if (paramMap["[[cutNames]]"].empty())
                    paramMap["[[cutNames]]"] += cut.GetName();
                else
                    paramMap["[[cutNames]]"] += "," + (std::string)cut.GetName();
                if (paramMap["[[cuts]]"].empty())
                    paramMap["[[cuts]]"] += cut.GetTitle();
                else
                    paramMap["[[cuts]]"] += " && " + (std::string)cut.GetTitle();
            }
        }
 
        paramMap["[[panelCutNames]]"] = "";
        paramMap["[[panelCuts]]"] = "";
        if (panel.panelCut) {
            for (const auto& cut : panel.panelCut->GetCuts()) {
                if (paramMap["[[panelCutNames]]"].empty())
                    paramMap["[[panelCutNames]]"] += cut.GetName();
                else
                    paramMap["[[panelCutNames]]"] += "," + (std::string)cut.GetName();
                if (paramMap["[[panelCuts]]"].empty())
                    paramMap["[[panelCuts]]"] += cut.GetTitle();
                else
                    paramMap["[[panelCuts]]"] += " && " + (std::string)cut.GetTitle();
            }
        }
 
        RESTInfo << "Global cuts: " << paramMap["[[cuts]]"] << RESTendl;
        if (!paramMap["[[panelCuts]]"].empty())
            RESTInfo << "Additional panel cuts: " << paramMap["[[panelCuts]]"] << RESTendl;
 
        // Replace panel variables and generate a TLatex label
        for (const auto& [key, posLabel] : panel.variablePos) {
            auto&& [variable, label, units] = key;
            bool found = false;
            std::string var = variable;
            if (!var.empty()) {
                for (const auto& [param, val] : paramMap) {
                    if (var == param) {
                        size_t pos = 0;
                        var = Replace(var, param, val, pos);
                        found = true;
                        break;
                    }
                }
                if (!found) RESTWarning << "Variable " << variable << " not found" << RESTendl;
            }
            if (isANumber(var)) {  // if it is a number, convert it to number to apply precision
                if (var.find('.') == std::string::npos) {
                    int dblVar = StringToInteger(var);
                    var = StringWithPrecision(dblVar, panel.precision);
                } else {
                    double dblVar = StringToDouble(var);
                    var = StringWithPrecision(dblVar, panel.precision);
                }
            }
            std::string lab =
                label + panel.delimiter.Data() + StringWithPrecision(var, panel.precision) + " " + units;
            panel.text.emplace_back(new TLatex(posLabel.X(), posLabel.Y(), lab.c_str()));
        }
 
        // Replace metadata variables and generate a TLatex label
        for (const auto& [key, posLabel] : panel.metadataPos) {
            auto&& [metadata, label, units] = key;
            std::string value = "";
 
            for (const auto& [name, quant] : quantity) {
                if (quant.metadata == metadata) value = quant.value;
            }
 
            if (value.empty()) {
                RESTWarning << "Metadata quantity " << metadata << " not found in dataSet" << RESTendl;
                continue;
            }
 
            std::string lab =
                label + panel.delimiter.Data() + StringWithPrecision(value, panel.precision) + " " + units;
            panel.text.emplace_back(new TLatex(posLabel.X(), posLabel.Y(), lab.c_str()));
        }
 
        // Replace observable variables and generate a TLatex label
        for (const auto& [key, posLabel] : panel.obsPos) {
            auto&& [obs, label, units] = key;
            auto value = *dataFrame.Mean(obs);
 
            std::string lab =
                label + panel.delimiter.Data() + StringWithPrecision(value, panel.precision) + " " + units;
            panel.text.emplace_back(new TLatex(posLabel.X(), posLabel.Y(), lab.c_str()));
        }
 
        // Replace any expression and generate TLatex label
        for (const auto& [key, posLabel] : panel.expPos) {
            auto&& [text, label, units, precision] = key;
            std::string var = text;
 
            // find and split the text into subtexts which are surrounded by {}
            std::vector<std::string> subtexts;
            while (var.find_last_of('{') != std::string::npos) {
                size_t posOpen = var.find_last_of('{');
                size_t posClose = var.find_first_of('}', posOpen);
                if (posClose == std::string::npos) {
                    RESTWarning << "Unmatched { in expression: " << var << RESTendl;
                    break;
                }
                std::string subtext = var.substr(posOpen + 1, posClose - posOpen - 1);
                subtexts.push_back(subtext);
                // get rid of "{"+subtext+"}" from the var
                var.erase(posOpen, posClose - posOpen + 1);
            }
            var = text;  // reset var to the original text
 
            // get precision formatting
            std::vector<std::string> precisionParts;
            for (const auto& part : Split(precision, ",", false, true)) {
                precisionParts.push_back(part);
            }
            if (precisionParts.size() != subtexts.size()) {
                RESTDebug << "Not enough precision values provided for the expression: `" << text
                          << "`. Using " << precisionParts.back() << " for last subtexts." << RESTendl;
                precisionParts.resize(subtexts.size(), precisionParts.back());  // fill with last value
            }
 
            size_t precisionIndex = precisionParts.size() - 1;
            for (const auto& subtext : subtexts) {
                std::string subVar = subtext;
                // replace variables
                for (const auto& [param, val] : paramMap) {
                    subVar = Replace(subVar, param, val);
                }
                // replace metadata
                for (const auto& [name, quant] : quantity) {
                    subVar = Replace(subVar, "[" + name + "]", quant.value);  // metadata are surrounded by []
                    subVar = Replace(subVar, name, quant.value);              // just in case ?
                }
                // replace observables
                for (const auto& obs : dataFrame.GetColumnNames()) {
                    if (var.find(obs) == std::string::npos) continue;
                    // here there should be a checking that the mean(obs) can be calculated
                    // (checking obs data type?)
                    double value = *dataFrame.Mean(obs);
                    subVar = Replace(subVar, obs, DoubleToString(value));
                }
                subVar = Replace(subVar, "[", "(");
                subVar = Replace(subVar, "]", ")");
                subVar = EvaluateExpression(subVar);
                if (isANumber(subVar)) {
                    double value = StringToDouble(subVar);
                    std::string prec = precisionParts[precisionIndex--];
                    if (prec.find('%') != std::string::npos) {  // it is a format e.g. %.2f
                        subVar = DoubleToString(value, prec);
                    } else if (isANumber(prec)) {  // it is a number
                        subVar = StringWithPrecision(value, StringToInteger(prec));
                    } else {
                        RESTWarning << "Unknown precision format: " << prec
                                    << ". Using panels precision: " << panel.precision << RESTendl;
                        subVar = StringWithPrecision(value, panel.precision);
                    }
                }
                var = Replace(var, "{" + subtext + "}", subVar);
            }
 
            std::string lab = label + panel.delimiter.Data() + var + " " + units;
            panel.text.emplace_back(new TLatex(posLabel.X(), posLabel.Y(), lab.c_str()));
        }
 
        // Draw the labels inside the pad
        for (const auto& text : panel.text) {
            text->SetTextColor(1);
            text->SetTextSize(panel.font_size);
            text->Draw("same");
        }
        canvasIndex++;
    }
 
    for (auto& plots : fPlots) {
        // Histograms are added to a THStack and will be ploted later on
        combinedCanvas.cd(canvasIndex);
        plots.hs = new THStack(plots.name.c_str(), plots.title.c_str());
        if (plots.legendOn) plots.legend = new TLegend(fLegendX1, fLegendY1, fLegendX2, fLegendY2);
        for (auto& hist : plots.histos) {
            auto dataFrame = dataSet.MakeCut(hist.histoCut);
            if (hist.variable.front() == "timeStamp") {
                hist.range.front().SetX(startTime);
                hist.range.front().SetY(endTime);
            }
            // 1-D Histograms
            if (hist.variable.size() == 1) {
                auto histo = dataFrame.Histo1D({hist.name.c_str(), hist.name.c_str(), hist.nBins.front(),
                                                hist.range.front().X(), hist.range.front().Y()},
                                               hist.variable.front());
                hist.histo = static_cast<TH1*>(histo->DrawClone());
                // 2-D Histograms
            } else if (hist.variable.size() == 2) {
                auto histo = dataFrame.Histo2D(
                    {hist.name.c_str(), hist.name.c_str(), hist.nBins.front(), hist.range.front().X(),
                     hist.range.front().Y(), hist.nBins.back(), hist.range.back().X(), hist.range.back().Y()},
                    hist.variable.front(), hist.variable.back());
                hist.histo = static_cast<TH1*>(histo->DrawClone());
            } else {
                RESTError << "Only 1D or 2D histograms are supported " << RESTendl;
                continue;
            }
            hist.histo->SetLineColor(hist.lineColor);
            hist.histo->SetLineWidth(hist.lineWidth);
            hist.histo->SetLineStyle(hist.lineStyle);
            hist.histo->SetFillColor(hist.fillColor);
            hist.histo->SetFillStyle(hist.fillStyle);
            // If stats are on histos must we drawn
            if (hist.statistics) {
                hist.histo->SetStats(true);
                hist.histo->Draw();
                combinedCanvas.Update();
            } else {
                hist.histo->SetStats(false);
            }
            // Normalize histos
            if (plots.normalize > 0) {
                const double integral = hist.histo->Integral();
                if (integral > 0) hist.histo->Scale(plots.normalize / integral);
            }
            // Scale histos
            if (plots.scale != "") {
                std::string inputScale = plots.scale;
                double binSize = hist.histo->GetXaxis()->GetBinWidth(1);
                double entries = hist.histo->GetEntries();
                double runLength = dataSet.GetTotalTimeInSeconds() / 3600.;  // in hours
                double integral = hist.histo->Integral("width");
 
                inputScale = Replace(inputScale, "binSize", DoubleToString(binSize));
                inputScale = Replace(inputScale, "entries", DoubleToString(entries));
                inputScale = Replace(inputScale, "runLength", DoubleToString(runLength));
                inputScale = Replace(inputScale, "integral", DoubleToString(integral));
 
                std::string scale = "1./(" + inputScale + ")";
                hist.histo->Scale(StringToDouble(EvaluateExpression(scale)));  // -1 if 'scale' isn't valid
            }
 
            // Add histos to the THStack
            plots.hs->Add(hist.histo, hist.drawOption.c_str());
            // Add histos to the legend
            if (plots.legend != nullptr) plots.legend->AddEntry(hist.histo, hist.histo->GetName(), "lf");
        }
    }
 
    // This function do the actual drawing of the THStack with the different options
    for (auto& plots : fPlots) {
        if (plots.hs == nullptr) continue;
        // TPad parameters
        TPad* targetPad = (TPad*)combinedCanvas.cd(canvasIndex);
        targetPad->SetLogx(plots.logX);
        targetPad->SetLogy(plots.logY);
        targetPad->SetLogz(plots.logZ);
        targetPad->SetGridx(plots.gridX);
        targetPad->SetGridy(plots.gridY);
        targetPad->SetLeftMargin(plots.marginLeft);
        targetPad->SetRightMargin(plots.marginRight);
        targetPad->SetBottomMargin(plots.marginBottom);
        targetPad->SetTopMargin(plots.marginTop);
 
        // HStack draw parameters
        plots.hs->Draw(plots.stackDrawOption.c_str());
        plots.hs->GetXaxis()->SetTitle(plots.labelX.c_str());
        plots.hs->GetYaxis()->SetTitle(plots.labelY.c_str());
        plots.hs->GetXaxis()->SetLabelSize(1.1 * plots.hs->GetXaxis()->GetLabelSize());
        plots.hs->GetYaxis()->SetLabelSize(1.1 * plots.hs->GetYaxis()->GetLabelSize());
        plots.hs->GetXaxis()->SetTitleSize(1.1 * plots.hs->GetXaxis()->GetTitleSize());
        plots.hs->GetYaxis()->SetTitleSize(1.1 * plots.hs->GetYaxis()->GetTitleSize());
 
        if (plots.timeDisplay) plots.hs->GetXaxis()->SetTimeDisplay(1);
        if (plots.legend != nullptr) plots.legend->Draw();
 
        targetPad->Update();
        combinedCanvas.Update();
        canvasIndex++;
    }
 
    // Preview plot. User can make some changed before saving
    if (!REST_Display_CompatibilityMode && fPreviewPlot) {
        combinedCanvas.Resize();
        GetChar();
    }
 
    // Save single pads if save is marked
    for (auto& plots : fPlots) {
        if (plots.save.empty()) continue;
        std::unique_ptr<TCanvas> canvas(new TCanvas());
        canvas->SetLogx(plots.logX);
        canvas->SetLogy(plots.logY);
        canvas->SetLogz(plots.logZ);
        canvas->SetGridx(plots.gridX);
        canvas->SetGridy(plots.gridY);
        canvas->SetLeftMargin(plots.marginLeft);
        canvas->SetRightMargin(plots.marginRight);
        canvas->SetBottomMargin(plots.marginBottom);
        canvas->SetTopMargin(plots.marginTop);
        plots.hs->Draw(plots.stackDrawOption.c_str());
        canvas->Print(plots.save.c_str());
    }
 
    // Save combined canvas
    if (!fOutputFileName.empty()) {
        for (const auto& [name, quant] : quantity) {
            size_t pos = 0;
            fOutputFileName = Replace(fOutputFileName, quant.metadata, quant.value, pos);
        }
        combinedCanvas.Print(fOutputFileName.c_str());
        // In case of root file save also the histograms
        if (TRestTools::GetFileNameExtension(fOutputFileName) == "root") {
            std::unique_ptr<TFile> f(TFile::Open(fOutputFileName.c_str(), "UPDATE"));
            for (auto& plots : fPlots) {
                for (auto& hist : plots.histos) {
                    hist.histo->Write();
                }
            }
            this->Write();
        }
    }
 
    CleanUp();
}
 
void TRestDataSetPlot::CleanUp() {
    for (auto& plots : fPlots) {
        for (auto& hist : plots.histos) {
            delete hist.histo;
        }
        delete plots.hs;
        delete plots.legend;
    }
 
    for (auto& panel : fPanels) {
        for (auto& text : panel.text) {
            delete text;
        }
        panel.text.clear();
    }
}
 
Int_t TRestDataSetPlot::GetIDFromMapString(const std::map<std::string, int>& mapStr, const std::string& in) {
    if (in.find_first_not_of("0123456789") == std::string::npos) {
        return StringToInteger(in);
    }
    auto it = mapStr.find(in);
    if (it != mapStr.end()) {
        return it->second;
    } else {
        RESTWarning << "cannot find ID with name \"" << in << "\"" << RESTendl;
    }
    return -1;
}
 
void TRestDataSetPlot::PrintMetadata() {
    TRestMetadata::PrintMetadata();
 
    RESTMetadata << "DataSet name: " << fDataSetName << RESTendl;
    RESTMetadata << "PaletteStyle: " << fPaletteStyle << RESTendl;
    if (fPreviewPlot) RESTMetadata << "Preview plot is ACTIVE" << RESTendl;
    RESTMetadata << "Canvas size: (" << fCanvasSize.X() << " ," << fCanvasSize.Y() << ")" << RESTendl;
    RESTMetadata << "Canvas divisions: (" << fCanvasDivisions.X() << " ," << fCanvasDivisions.Y() << ")"
                 << RESTendl;
    RESTMetadata << "-------------------" << RESTendl;
    for (const auto& plot : fPlots) {
        RESTMetadata << "-------------------" << RESTendl;
        RESTMetadata << "Plot name/title: " << plot.name << " " << plot.title << RESTendl;
        RESTMetadata << "Save string: " << plot.save << RESTendl;
        RESTMetadata << "Set log X,Y,Z: " << plot.logX << ", " << plot.logY << ", " << plot.logZ << RESTendl;
        RESTMetadata << "Stack draw Option: " << plot.stackDrawOption << RESTendl;
        if (plot.legend) RESTMetadata << "Legend is ON" << RESTendl;
        if (plot.timeDisplay) RESTMetadata << "Time display is ON" << RESTendl;
        RESTMetadata << "Labels X,Y: " << plot.labelX << ", " << plot.labelY << RESTendl;
        for (const auto& hist : plot.histos) {
            RESTMetadata << "****************" << RESTendl;
            RESTMetadata << "Histo name: " << hist.name << RESTendl;
            RESTMetadata << "Draw Option: " << hist.drawOption << RESTendl;
            RESTMetadata << "Histogram size: " << hist.variable.size() << " with parameters:" << RESTendl;
            for (size_t i = 0; i < hist.variable.size(); i++) {
                RESTMetadata << "\t" << i << " " << hist.variable[i] << ", " << hist.nBins[i] << ", "
                             << hist.range[i].X() << ", " << hist.range[i].Y() << RESTendl;
            }
            RESTMetadata << "****************" << RESTendl;
        }
    }
    RESTMetadata << "-------------------" << RESTendl;
    for (auto& panel : fPanels) {
        RESTMetadata << "-------------------" << RESTendl;
        RESTMetadata << "Panel font size/precision " << panel.font_size << ", " << panel.precision
                     << RESTendl;
        RESTMetadata << "****************" << RESTendl;
        for (auto& [key, posLabel] : panel.variablePos) {
            auto&& [obs, label, units] = key;
            RESTMetadata << "Label variable " << obs << ", label " << label << ", units " << units << " Pos ("
                         << posLabel.X() << ", " << posLabel.Y() << ")" << RESTendl;
        }
        RESTMetadata << "****************" << RESTendl;
        for (auto& [key, posLabel] : panel.metadataPos) {
            auto&& [obs, label, units] = key;
            RESTMetadata << "Label metadata " << obs << ", label " << label << ", units " << units << " Pos ("
                         << posLabel.X() << ", " << posLabel.Y() << ")" << RESTendl;
        }
        RESTMetadata << "****************" << RESTendl;
        for (auto& [key, posLabel] : panel.obsPos) {
            auto&& [obs, label, units] = key;
            RESTMetadata << "Label Observable " << obs << ", label " << label << ", units " << units
                         << " Pos (" << posLabel.X() << ", " << posLabel.Y() << ")" << RESTendl;
        }
        RESTMetadata << "****************" << RESTendl;
        for (auto& [key, posLabel] : panel.expPos) {
            auto&& [obs, label, units, precision] = key;
            RESTMetadata << "Label Expression " << obs << ", label " << label << ", units " << units
                         << ", precision " << precision << " Pos (" << posLabel.X() << ", " << posLabel.Y()
                         << ")" << RESTendl;
        }
        RESTMetadata << "****************" << RESTendl;
    }
    RESTMetadata << "-------------------" << RESTendl;
 
    RESTMetadata << RESTendl;
}