#include "RooMomentMorph.h"
#include "RooWorkspace.h"
#include "RooRealVar.h"
#include "TROOT.h"
#include "TH1.h"
#include "TH1D.h"
#include "TFile.h"

#include <algorithm>
#include <fstream>
#include <string>

auto getHistnames(std::string const& path) {

    std::map<std::string, double> histnames;

    std::ifstream openFile(path);

    if(!openFile)
    {
        std::cout << "Error opening file!" << std::endl;
    } 

    std::string key;
    std::string valueString;

    while(true) {
        if (!std::getline(openFile, key, ':')) break;

        std::getline(openFile, valueString, '\n');
        histnames[key] = std::stof(valueString);
    }

    return histnames;
}

void script_fast() {

    using namespace RooFit;

    gROOT->ProcessLine("gErrorIgnoreLevel = 1001;"); // Kill some CHI2 warnings
    RooMsgService::instance().setGlobalKillBelow(WARNING); // Kill some long RooFit Info log

    // Load the morphing

    TFile morph_file("morphing.root");

    auto workspace = morph_file.Get<RooWorkspace>("MorphWorkspace");
    auto delayVar  = workspace->var("delayVar");
    auto bxidVar   = workspace->var("bxidVar");
    auto morphing  = static_cast<RooMomentMorph*>(workspace->pdf("morph"));
    morphing->useHorizontalMorphing(true); // False works fine, True produces memory leak

    morph_file.Close();

    // Load the histograms

    TFile hist_file("hists.root");

    std::string TH1_dir = "DQMData/Run 1/Ph2TkBXHist/Run summary/Hist1D";

    auto histnames = getHistnames("hists.yml");

    // Bunch crossing range
    std::vector<double> bxRange;
    for (int i = 0; i < 500; ++i) {
        bxRange.push_back(0.1 * i);
    }

    // Interpolate the true hists
    std::vector<TH1*> trueHists;
    for (int i = 0; i < bxRange.size(); ++i) {
        delayVar->setVal(bxRange[i]);
        auto histname = std::string("trueHist") + std::to_string(bxRange[i]);
        trueHists.push_back(morphing->createHistogram(histname.c_str(),*bxidVar));
        std::cout << "Created morphed histogram " << i << std::endl;
    }

    // Function to estimate the reco delay
    auto getRecoDelay = [&](TH1 const& recoHist) {

        std::vector<double> chi2;

        for (int i = 0; i < bxRange.size(); ++i) {
            chi2.push_back(recoHist.Chi2Test(trueHists[i],"UU CHI2"));
        }
        auto argmin = std::min_element( chi2.begin(), chi2.end() ) - chi2.begin();
        // print for debgging
        //std::cout << chi2[argmin] << std::endl;
        return bxRange[argmin];
    };


    // Run the "algo"

    for (auto const& item : histnames) {
        auto const& name = item.first;
        auto const& true_delay = item.second;

        auto h = hist_file.Get<TH1>((TH1_dir+'/'+name).c_str());
        auto reco_delay = getRecoDelay(*h);
        std::cout << "True = " << true_delay << ", reco = " << reco_delay << std::endl;
    }

    hist_file.Close();
}