#include <iostream>
#include <fstream>
#include <cstring>
#include <vector>
#include "TCanvas.h"
#include "TFile.h"
#include "TSystem.h"
#include "TStyle.h"
#include "TApplication.h"
#include "TMultiGraph.h"
#include "TRint.h"
#include "TH1.h"
#include "TF1.h"
#include "TH2.h"
#include "TH3.h"
#include "TKey.h"
#include "TMath.h"
#include "TPDF.h"
#include <TGaxis.h>
#include <TLatex.h>
#include <TLegend.h>
#include <TLegendEntry.h>
#include <TGraphErrors.h>
#include <Rtypes.h>
#include "TH2D.h"
#include "THStack.h"
#include "TPaveStats.h"
#include "TPaletteAxis.h"
#include "TRandom.h"
#include <math.h>
#include "TList.h"
#include "Math/WrappedMultiTF1.h"
#include "HFitInterface.h"
#include "Fit/Fitter.h"
#include "Fit/BinData.h"
#include "Fit/Chi2FCN.h"
using namespace std;

///// Here I find the right parameters for tau, tau is parameter 0, and now both of them are the same.
int iparB[3]  = {0,1,2};
int iparSB[3] = {0,3,4};

// Create the GlobalCHi2 structure
struct GlobalChi2 {
        GlobalChi2(  ROOT::Math::IMultiGenFunction & f1,
                     ROOT::Math::IMultiGenFunction & f2) :
                fChi2_1(&f1), fChi2_2(&f2) {
        }
        // parameter vector is first background (in common 1 and 2)
        // and then is signal (only in 2)
        double operator() (const double *par) const {
                double p1[3];
                for (int i = 0; i < 3; ++i) p1[i] = par[iparB[i] ];
                double p2[3];
                for (int i = 0; i < 3; ++i) p2[i] = par[iparSB[i] ];
                return (*fChi2_1)(p1) + (*fChi2_2)(p2);
        }
        const ROOT::Math::IMultiGenFunction * fChi2_1;
        const ROOT::Math::IMultiGenFunction * fChi2_2;
};

std::vector < pair < int, int >>
runinfo;
unsigned int
StartRun (12061),
EndRun (19698);
float firstday = TDatime::GetGlobalDayFromDate (20150401); //"2015-04-01 00:00:00"); // The day when UAr fully filled
void run_days ()
{
        std::ifstream in ("/home/iftikhar/Templates/RunTimeStamp.txt",
                          std::ios::in);
        if (in.fail ())
        {
                cout << "File: RunTimeStamp.txt not exist." << endl;
                return;
        }
        while (true)
        {
                unsigned int RunId;
                unsigned int year,month,day,hour,minute, second;
                double sec;
                in >> RunId >> year >> month >> day >> hour >> minute >> sec;
                if (in.eof ())
                        break;
                second = (unsigned int) sec;
                if (RunId < StartRun)
                        continue;
                if (RunId > EndRun)
                        continue;
                TDatime
                date (year, month, day, hour, minute, second);
                int
                        dday = TDatime::GetGlobalDayFromDate (date.GetDate ()) - firstday;
                std::pair < int, int >
                run = { RunId, dday };
                runinfo.push_back (run);
                //cout<<"Run: "<<RunId<<" day: "<<dday<<endl;
        }
}
Int_t GetDays (Int_t run_id)
{
        auto it = std::find_if (runinfo.begin (), runinfo.end (),
                                [run_id] (const std::pair < int, int >&run)
                                { return run.first == run_id;}
                                );
        return it->second; // return days fropm 2015-04-01
}

void normalized()
{
        gStyle->SetOptStat(0);
        gStyle->SetOptFit(11111);
        run_days();
        TFile *file = new TFile("/home/iftikhar/Templates/UAr_70d_mt2_SLAD_v3_5_0_Hist.root");
        TFile *Background_data =new TFile("/home/iftikhar/Templates/UAr_500d_SLAD_v3_5_0_merged_v4_Hist.root");

        TH2F *Bkg_hist;
        TH1F *Bkg_live;
        TH1D *hLivetime;
        TH2F *Ar_hist;

        Ar_hist       = (TH2F *) file->Get   ("hNePerRun");
        Bkg_live      = (TH1F *) Background_data->Get("hLivetimePerRun");
        Bkg_hist      = (TH2F *) Background_data->Get("hNePerRun_bkg"  );
        TString hname = "hNePerRun";
        TH2D *htmp    = (TH2D*) file->Get(hname.Data());
        hLivetime     = (TH1D*) file->Get("hLivetimePerRun");

        TH1D *h_runs_peak1    =new TH1D ("Runs_peak 1",  "Counts in each run for peak 1", 1121,21,123);// 12061 - 0.5,13181 + 0.5);
        TH1D *h_runs_peak2    =new TH1D ("Runs_peak 2",  "Counts in each run for peak 2", 1121,21,123);// 12061 - 0.5,13181 + 0.5);
        TH1D *h_days_peak1    =new TH1D ("Days_peak 1",  "Counts in each Day for Peak 1", 1121, 21, 123);
        TH1D *h_days_peak2    =new TH1D ("Days_peak 2",  "Counts in each Day for peak 2", 1121, 21, 123);
        TH1D *h_rate_peak1    =new TH1D ("Rate peak 1",  "Rate for Peak 1", 1121, 21, 123);
        TH1D *h_rate_peak2    =new TH1D ("Rate peak 2",  "Rate for Peak 2", 1121, 21, 123);
        TH1D *hlive_days      =new TH1D ("livetimedays", "Livetime vs days",1121, 21, 123);
        TH1D *h_runs1         =new TH1D ("Runs1", "Counts in each run for peak 1", 1121, 21,123);//12061 - 0.5,13181 + 0.5);

        Double_t lowbin=  Bkg_live->FindBin(13204.);
        Double_t highbin=  Bkg_live->FindBin (19698.);
        TH1D *bkg_slice=Bkg_hist->ProjectionY("slice_bkg",lowbin,highbin);
        double counts_integral_five_hundred  =Bkg_live->Integral (lowbin,highbin);
        Double_t lowb=  hLivetime->FindBin(12061.);
        Double_t highb=  hLivetime->FindBin (13181.);

        //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
        Int_t dt_days = 1;
        Int_t day_min = (GetDays(htmp->GetXaxis()->GetBinCenter(1))/dt_days)*dt_days; // convert lower bin edge in unit of dt_days
        Int_t day_max= ((GetDays(htmp->GetXaxis()->GetBinCenter(htmp->GetNbinsX()))/dt_days)+1)*dt_days;;// convert higher bin edge in unit of dt_days
        Int_t nbins = (day_max-day_min)/dt_days+1;
        Int_t nbiny = htmp->GetNbinsY();
        Int_t Ne_min = htmp->GetYaxis()->GetBinLowEdge(1);
        Int_t Ne_max = htmp->GetYaxis()->GetBinUpEdge(nbiny);

        TString title = Form("; days since 2015-04-01; %s", htmp->GetYaxis()->GetTitle());
        TH2D *runs_days = new TH2D((hname + "_merged").Data(), title.Data(), nbins, day_min-0.5, day_max+0.5, nbiny, Ne_min,Ne_max);
        TH2D *RunId_livetime = new TH2D("hLTmarged", title.Data(), nbins, day_min-0.5, day_max+0.5, nbiny,Ne_min,Ne_max);
        TH2D *h1 = new TH2D("h1", title.Data(), nbins, day_min-0.5, day_max+0.5, nbiny,Ne_min,Ne_max);

        for (int ibin = 1; ibin <= htmp->GetNbinsX(); ++ibin) { // loop for run
                Double_t day =  GetDays(htmp->GetXaxis()->GetBinCenter(ibin));
                Double_t lvt = hLivetime->GetBinContent(ibin);
                for (int iybin = 1; iybin <= htmp->GetNbinsY(); ++iybin) { // loop for time
                        Double_t entry = htmp->GetBinContent(ibin, iybin);
                        Double_t Ne = htmp->GetYaxis()->GetBinCenter(iybin);
                        runs_days->Fill(day, Ne, entry);
                        RunId_livetime->Fill(day, Ne, lvt);
                }
        }
        TCanvas *c8 = new TCanvas ("c8", "", 800, 600);
        c8->Divide(2,2);
        c8->cd(1);
        runs_days->SetTitle("Events in each day");
        runs_days->Draw("colz");
        c8->cd(2);
        RunId_livetime->Draw("colz");
        RunId_livetime->SetTitle("Livetime");
        c8->cd(3);
        h1 = (TH2D*)runs_days->Clone();
        h1->SetTitle("Rate for each day");
        h1->Divide(RunId_livetime);// The problem is here in my opinion.
        h1->Draw("colz");


        TCanvas *c_l = new TCanvas ("c_l", "background subtraction with scaling", 1800, 1200); c_l->Divide(2,2);
        for (int64_t ibin = 1; ibin <=99; ++ibin) {
                //        TCanvas *c = new TCanvas(TString::Format("Days_%ld", (ibin+21) + 1),TString::Format("Days_%ld", (ibin+21) + 1), 780, 600);
                TH1D *five_hundred_days = (TH1D*)bkg_slice->Clone("five_hundred_days");
                TH1D *px_slices = h1->ProjectionY (Form ("slices_%zu", ibin), ibin, ibin, "");

                TH2D *no_subtraction_slices = (TH2D*)px_slices->Clone("no_subtraction_slices");
                five_hundred_days->Scale(1/counts_integral_five_hundred);

                c_l->cd(ibin)->SetLogy();
                no_subtraction_slices->Draw("HIST");
                no_subtraction_slices->GetYaxis ()->SetTitle ("Events");
                no_subtraction_slices->SetTitle(Form("Day %ld", ibin+22));
                px_slices->Add(five_hundred_days,-1);
                px_slices->SetLineColor(kGreen);
                px_slices->SetLineWidth(2);
                px_slices->Draw("HIST same");
                five_hundred_days->SetLineColor(kRed);
                five_hundred_days->SetLineWidth(2);
                five_hundred_days->Draw("HIST same");


                Int_t day =runs_days->GetXaxis ()->GetBinCenter (ibin);

                //Legends
                auto legend = new TLegend(0.7,0.7,0.9,0.9);
                legend->SetHeader("Ar-37 and Background","c");
                legend->AddEntry(no_subtraction_slices,"No subtraction");
                legend->AddEntry(five_hundred_days,"Scaled background");
                legend->AddEntry(px_slices," Ar^{37} data - Scaled background");
                legend->Draw();

                Int_t xmin_peak1 = px_slices->GetXaxis()->FindBin(6.);
                Int_t xmax_peak1 = px_slices->GetXaxis()->FindBin(19.);
                Int_t xmin_peak2 = px_slices->GetXaxis()->FindBin(25.);
                Int_t xmax_peak2 = px_slices->GetXaxis()->FindBin(70.);

                double counts_integral_peak1 =px_slices->Integral(xmin_peak1,xmax_peak1);
                if (counts_integral_peak1<0) continue;
                double counts_integral_peak2 =px_slices->Integral(xmin_peak2,xmax_peak2);
                if (counts_integral_peak2<0) continue;
                //cout<<counts_integral_peak1<<endl;

                h_runs_peak1->Fill(day,counts_integral_peak1);
                for(int i=1; i<=h_runs_peak1->GetNbinsX(); ++i) {
                        Double_t entry1=h_runs_peak1->GetBinContent(i);
                        Double_t error1=TMath::Sqrt(entry1);
                        h_runs_peak1->SetBinContent(i,entry1);
                        h_runs_peak1->SetBinError(i, error1);
                }
                h_runs_peak2->Fill(day,counts_integral_peak2);
                for(int i=1; i<=h_runs_peak2->GetNbinsX(); ++i) {
                        Double_t entry2=h_runs_peak2->GetBinContent(i);
                        Double_t error2=TMath::Sqrt(entry2);
                        h_runs_peak2->SetBinContent(i,entry2);
                        h_runs_peak2->SetBinError(i, error2);
                }
        }
        for (int64_t ibin = 1; ibin <=hLivetime->GetXaxis()->GetNbins(); ++ibin) {
                Int_t run_id = hLivetime->GetXaxis ()->GetBinCenter (ibin);
                double days = GetDays(run_id);
                Double_t live = hLivetime->GetBinContent (ibin);
                h_runs1->Fill(days,live);
                for(int i=1; i<=h_runs1->GetNbinsX(); ++i) {
                        Double_t entry0=h_runs1->GetBinContent(i);
                        Double_t error0=TMath::Sqrt(entry0);
                        h_runs1->SetBinContent(i,entry0);
                        h_runs1->SetBinError(i, error0);
                }
        }

        h_rate_peak1 = (TH1D *) h_runs_peak1->Clone ();
        h_rate_peak1->Sumw2();

        double xmin_days_peak1 = 21., xmax_days_peak1 = 123.;
        //TF1 *fexp1 = new TF1("fexp1", "[N_{1}]/[#tau]*exp(-x/[#tau])+[Constant]", xmin_days_peak1, xmax_days_peak1);
        TF1 *fexp1 = new TF1("fexp1", "[2]/[0]*exp(-x/[0])+[1]", xmin_days_peak1, xmax_days_peak1);
        fexp1->SetParameter(0, 43.);
        fexp1->SetParameter(1, 0.0002);
        fexp1->SetParameter(2, 0.023);
        fexp1->SetParNames("#tau","Constant","N_{1}");
        h_rate_peak2->Sumw2();
        h_rate_peak2 = (TH1D *) h_runs_peak2->Clone ();
        double xmin_days_peak2 = 21., xmax_days_peak2 = 123.;
        TF1 *fexp2 = new TF1("fexp2", "[2]/[0]*exp(-x/[0])+[1]", xmin_days_peak2,xmax_days_peak2);
        fexp2->SetParameter(0, 43.);
        fexp2->SetParameter(1, 0.002);
        fexp2->SetParameter(2, 0.18);
        fexp2->SetParNames("#tau","Constant","N_{1}");
        ROOT::Math::WrappedMultiTF1 wfB(*fexp1,1);
        ROOT::Math::WrappedMultiTF1 wfSB(*fexp2,1);
        ROOT::Fit::DataOptions opt;
        ROOT::Fit::DataRange rangeB;

        // set the data range
        rangeB.SetRange(21,123);
        ROOT::Fit::BinData dataB(opt,rangeB);
        ROOT::Fit::FillData(dataB, h_rate_peak1);

        ROOT::Fit::DataRange rangeSB;
        rangeSB.SetRange(21,123);
        ROOT::Fit::BinData dataSB(opt,rangeSB);
        ROOT::Fit::FillData(dataSB, h_rate_peak2);

        ROOT::Fit::Chi2Function chi2_B(dataB, wfB);
        ROOT::Fit::Chi2Function chi2_SB(dataSB, wfSB);

        GlobalChi2 globalChi2(chi2_B, chi2_SB);

        ROOT::Fit::Fitter fitter;
        //  double par0[Npar] = {25,19,0.5,80,18};//  double par0[Npar] = {49,-39,0.9,55,18};
        const int Npar = 5;

        double par0[Npar] = {43,0.02,0.0001,0.2,0.001};//18}; ///parameters {0,1,2}{0,3,4},//Here tau is fixed only
        // create before the parameter settings in order to fix or set range on them
        fitter.Config().SetParamsSettings(5,par0);
        fitter.Config().ParSettings(0).SetLimits(35,55);
        fitter.FitFCN(5,globalChi2,0,dataB.Size()+dataSB.Size(),true);
        ROOT::Fit::FitResult result = fitter.Result();
        result.Print(std::cout);

        TCanvas * c3 = new TCanvas("c3","Simultaneous fit of two histograms, Same constant,5 pars",10,10,800,700);
        c3->Divide(1,2);
        c3->SetGrid();

        c3->cd(1);
        gStyle->SetOptFit(11111);
        h_rate_peak1->GetYaxis()->SetRangeUser(0,0.00067);
        h_rate_peak1->SetNameTitle ("h_rate_peak1", "Rate for Peak 1");
        h_rate_peak1->GetXaxis ()->SetTitle ("Days_{From 2015-04-01}");
        h_rate_peak1->GetYaxis ()->SetTitle ("Rate (Hz)");
        h_rate_peak1->SetMarkerColor (kRed);
        h_rate_peak1->SetMarkerStyle (20);
        fexp1->SetFitResult( result, iparB);
        fexp1->SetRange(rangeB().first, rangeB().second);
        fexp1->SetLineColor(kBlue);
        h_rate_peak1->GetListOfFunctions()->Add(fexp1);

        h_rate_peak1->Draw("E1");

        c3->cd(2);
        h_rate_peak2->GetYaxis()->SetRangeUser(0,0.0045);
        h_rate_peak2->SetMarkerColor (kMagenta);
        h_rate_peak2->SetMarkerStyle (20);
        h_rate_peak2->SetNameTitle ("h_rate_peak2", "Rate for Peak 2");
        h_rate_peak2->GetXaxis ()->SetTitle ("Days_{From 2015-04-01}");
        h_rate_peak2->GetYaxis ()->SetTitle ("Rate (Hz)");
        fexp2->SetFitResult( result, iparSB);
        fexp2->SetRange(rangeSB().first, rangeSB().second);
        fexp2->SetLineColor(kRed);
        h_rate_peak2->GetListOfFunctions()->Add(fexp2);
        h_rate_peak2->Draw("E1");
}