// Include standard stuff
#include <cstdlib>
#include <cwchar>
#include <fstream>
#include <iostream>
#include <sstream>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <sys/stat.h>

// Include ROOT stuff
#include "Math/IFunction.h"
#include "TCanvas.h"
#include "TClonesArray.h"
#include "TCutG.h"
#include "TDirectory.h"
#include "TF1.h"
#include "TF2.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
#include "TLatex.h"
#include "TObject.h"
#include "TPDF.h"
#include "TProfile.h"
#include "TRatioPlot.h"
#include "TTree.h"
#include <Fit/BinData.h>
#include <Fit/Chi2FCN.h>
#include <Fit/Fitter.h>
#include <HFitInterface.h>
#include <Math/WrappedMultiTF1.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TStyle.h>

using namespace std;
TPDF pdff("result.pdf");

TFile *outfile = new TFile("data.root", "READONLY");

TGraph *simul1;
Double_t InterpolatedGraph(Double_t *x, Double_t *par) {
  return par[0] + par[1] * (simul1->Eval(x[0]));
};

TGraph *simul2;
Double_t InterpolatedGraph2(Double_t *x, Double_t *par) {
  return par[0] + par[1] * (simul2->Eval(x[0]));
};

int ipar[2] = {
    0, // intercept common parameter
    1  // slope common parameter
};

// Create the GlobalCHi2 structure

struct GlobalChi2 {
  GlobalChi2(ROOT::Math::IMultiGenFunction &f1,
             ROOT::Math::IMultiGenFunction &f2)
      : fChi2_1(&f1), fChi2_2(&f2) {}

  double operator()(const double *par) const {
    double p[2];
    for (int i = 0; i < 2; ++i)
      p[i] = par[ipar[i]];

    return (*fChi2_1)(p) + (*fChi2_2)(p);
  }

  const ROOT::Math::IMultiGenFunction *fChi2_1;
  const ROOT::Math::IMultiGenFunction *fChi2_2;
};

char filename[4];

int combinedFitTH2() {

  TH2F *htemp1 = (TH2F *)outfile->Get("htemp1");
  TH2F *htemp2 = (TH2F *)outfile->Get("htemp2");

  TCanvas *canvas = new TCanvas("canvas");
  canvas->Print("result.pdf[");

  canvas->SetLogz();

  simul1 = new TGraph("f1.txt", "%lg %lg");
  TF1 *f1 = new TF1("f1", InterpolatedGraph, 0, 200, 2);
  f1->SetParameters(0, 1);
  f1->SetRange(0, 600, 70, 1400);
  f1->SetParName(0, "intercept");
  f1->SetParName(1, "slope");

  simul2 = new TGraph("f2.txt", "%lg %lg");
  TF1 *f2 = new TF1("f2", InterpolatedGraph2, 0, 20, 2);
  f2->SetParameters(0, 1);
  f2->SetRange(0, 0, 20, 200);
  f2->SetParName(0, "intercept");
  f2->SetParName(1, "slope");

  ROOT::Math::WrappedMultiTF1 wf1(*f1, 0);
  ROOT::Math::WrappedMultiTF1 wf2(*f2, 0);

  ROOT::Fit::DataOptions opt;

  // set the data range
  ROOT::Fit::DataRange range1;
  range1.SetRange(0, 70, 600, 1400);
  ROOT::Fit::BinData data1(opt, range1);
  ROOT::Fit::FillData(data1, htemp1);

  ROOT::Fit::DataRange range2;
  range2.SetRange(0, 15, 0, 200);
  ROOT::Fit::BinData data2(opt, range2);
  ROOT::Fit::FillData(data2, htemp2);

  ROOT::Fit::Chi2Function chi2_f1(data1, wf1);
  ROOT::Fit::Chi2Function chi2_f2(data2, wf2);

  GlobalChi2 globalChi2(chi2_f1, chi2_f2);

  ROOT::Fit::Fitter fitter;

  const int Npar = 2;
  double par0[Npar] = {0, 1};

  // create before the parameter settings in order to fix or set range on them
  fitter.Config().SetParamsSettings(Npar, par0);
  // fitter.Config().ParSettings(0).SetLimits(0,200);
  // fitter.Config().ParSettings(1).SetLimits(0,100);

  fitter.Config().MinimizerOptions().SetPrintLevel(0);
  fitter.Config().SetMinimizer("Minuit2", "Migrad");

  // fit FCN function directly
  // (specify optionally data size and flag to indicate that is a chi2 fit)
  fitter.FitFCN(Npar, globalChi2, par0, data1.Size() + data2.Size(), true);
  ROOT::Fit::FitResult result = fitter.Result();
  result.Print(std::cout);

  std::cout << "Combined fit Chi2 = " << result.Chi2() << std::endl;

  // f1->SetFitResult(result, ipar);
  f1->SetLineColor(kBlue);
  f2->SetLineColor(kRed);
  // hSi->GetListOfFunctions()->Add(f1);
  // hSi->Draw();
  // f1->Draw("same");

  htemp1->Draw("colz");
  htemp2->Draw("samecolz");
  f1->Draw("same");
  f2->Draw("same");

  TLatex *tex3 = new TLatex(0, 713.64,
                            Form("Blue: a*f(x) + b, a = %.2f b = %.2f",
                                 f1->GetParameter(1), f1->GetParameter(0)));
  tex3->SetTextSize(0.03);
  tex3->SetLineWidth(1);
  tex3->Draw();

  TLatex *tex1 = new TLatex(0, 753.64,
                            Form("Red: a*f(x) + b, a = %.2f b = %.2f",
                                 f2->GetParameter(1), f2->GetParameter(0)));
  tex1->SetTextSize(0.03);
  tex1->SetLineWidth(1);
  tex1->Draw();

  canvas->Update();
  canvas->Print("result.pdf", "");

  canvas->Print("result.pdf]");
  return 0;
}