#include <filesystem>
#include <fstream>
#include <iostream>
#include <iterator>
#include <sstream>
#include <string>

using namespace std;

#include "ROOT/RDataFrame.hxx"
#include "TCanvas.h"
#include "TDirectory.h"
#include "TF1.h"
#include "TFile.h"
#include "TGaxis.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
#include "TKey.h"
#include "TMath.h"
#include "TPad.h"
#include "TPaveStats.h"
#include "TPaveText.h"
#include "TProfile.h"
#include "TROOT.h"
#include "TString.h"
#include "TStyle.h"
#include "TSystem.h"
#include "TTree.h"

TString filename_base;
//-------------------------------Parametwers--------------------------------
Double_t sample_inte_ns = 0.4;
Int_t cor_num = 10;
Int_t peak_bin = 50;
Int_t corr_amp_start[4] = {500, 500, 500, 500};
Int_t corr_amp_stop[4] = {3000, 2500, 3000, 3000};
Int_t time_bin = 50;
Double_t corr_time_sta_1st = -500;
Double_t corr_time_end_1st = 4000;
Int_t corr_time_sta = -2000, corr_time_end = 2000;

void draw_amp_and_time(Int_t n_corr, vector<Double_t> &deltaT,
                       vector<vector<Double_t>> &peak, Int_t pbin,
                       Int_t *amp_start, Int_t *amp_stop, Int_t tbin,
                       Int_t time_sta, Int_t time_end, TString label,
                       TString label2, TString titlepart, Int_t *targetID) {
  Int_t valid_evt = deltaT.size();
  TCanvas *c_time = new TCanvas("c_time", "Time Resolution", 800, 600);
  TH1D *hist_T = new TH1D("h_time", titlepart + " " + label2 + " slewing", 100,
                          time_sta, time_end);
  for (Int_t i = 0; i < valid_evt; i++) {
    hist_T->Fill(deltaT[i]);
  }
  hist_T->Fit("gaus", "Q");
  hist_T->GetXaxis()->SetTitle("Time[ps]");
  hist_T->GetXaxis()->CenterTitle();
  hist_T->GetYaxis()->SetTitle("Count");
  hist_T->GetYaxis()->CenterTitle();
  gPad->Update();
  // c_time->Print(label + "_resolution_" + label2 + ".pdf");
  // Double_t sigma = ((TF1 *)hist_T->FindObject("gaus"))->GetParameter(2);
  // ofstream out(filename_base + "_time_resolution.txt", ios::app);
  // out << label << "\t" << sigma << endl;
  // out.close();
  delete c_time;
  delete hist_T;

  TCanvas *c_slewing;
  switch (n_corr) {
  case 2:
    c_slewing = new TCanvas("c_slewing", "slewing", 800, 400);
    c_slewing->Divide(2, 1);
    break;
  case 3:
    c_slewing = new TCanvas("c_slewing", "slewing", 1200, 400);
    c_slewing->Divide(3, 1);
    break;
  case 4:
    c_slewing = new TCanvas("c_slewing", "slewing", 800, 800);
    c_slewing->Divide(2, 2);
    break;
  default:
    break;
  }

  TH2D *hist2D[4]{};
  for (Int_t ch = 0; ch < n_corr; ch++) {
    TString name = "Sig" + TString::Itoa(targetID[ch], 10);
    hist2D[ch] = new TH2D(name, name, pbin, amp_start[ch], amp_stop[ch], tbin,
                          time_sta, time_end);
    for (Int_t i = 0; i < valid_evt; i++) {
      hist2D[ch]->Fill(peak[ch][i], deltaT[i]);
    }
    c_slewing->cd(ch + 1);
    hist2D[ch]->Draw("colz");
    hist2D[ch]->GetXaxis()->SetTitle("Amplitude[mV]");
    hist2D[ch]->GetXaxis()->CenterTitle();
    hist2D[ch]->GetYaxis()->SetTitle("Time[ps]");
    hist2D[ch]->GetYaxis()->CenterTitle();
    gPad->Update();
  }
  // c_slewing->Print(label + "_" + label2 + ".pdf");
  delete c_slewing;
  for (Int_t ch = 0; ch < n_corr; ch++) {
    delete hist2D[ch];
  }
}

void conduct_slewing(Int_t n_corr, Int_t n_fit, vector<Double_t> &deltaT,
                     vector<vector<Double_t>> &peak, Int_t cor_num,
                     TString label, TString titlepart, Int_t *targetID) {
  Int_t valid_evt = deltaT.size();

  draw_amp_and_time(n_corr, deltaT, peak, peak_bin, corr_amp_start,
                    corr_amp_stop, time_bin, corr_time_sta_1st,
                    corr_time_end_1st, label, "before", titlepart, targetID);

  //
  TCanvas c_temp("c_temp", "c_temp", 800, 600);
  // c_temp.Print(label + "_resolution_after.pdf[");
  // c_temp.SetWindowSize(n_corr * 400, 400);
  // c_temp.Print(label + "_after.pdf[");

  TCanvas *c_work = new TCanvas("work", "work", n_corr * 400, 800);
  c_work->Divide(n_corr, 2);
  // c_work->Print(label + "_process.pdf[");
  Double_t min_amp = *min_element(corr_amp_start, corr_amp_start + n_corr);
  Double_t max_amp = *max_element(corr_amp_stop, corr_amp_stop + n_corr);
  TString formula = "[0]+[1]/sqrt(x)+[2]/x";
  for (Int_t i = 0; i < n_fit; i++) {
    formula += TString::Format("+[%d]*x**%d", i + 3, i + 1);
  }
  TF1 *myfit = new TF1("myfit", formula, min_amp, max_amp);

  for (Int_t ir = 1; ir <= cor_num; ir++) {
    TH2D *before_corr[4]{};
    TH2D *after_corr[4]{};
    for (Int_t ch = 0; ch < n_corr; ch++) {
      TString label_process = "round_" + TString::Itoa(ir, 10) + "_ch_" +
                              TString::Itoa(targetID[ch], 10);
      if (ir == 1 && ch == 0) {
        before_corr[ch] =
            new TH2D(label_process + "_before", label_process + "_before",
                     peak_bin, corr_amp_start[ch], corr_amp_stop[ch], time_bin,
                     corr_time_sta_1st, corr_time_end_1st);
      } else {
        before_corr[ch] =
            new TH2D(label_process + "_before", label_process + "_before",
                     peak_bin, corr_amp_start[ch], corr_amp_stop[ch], time_bin,
                     corr_time_sta, corr_time_end);
      }

      for (Int_t i = 0; i < valid_evt; i++) {
        before_corr[ch]->Fill(peak[ch][i], deltaT[i]);
      }
      c_work->cd(ch + 1);
      before_corr[ch]->Draw("colz");
      before_corr[ch]->GetXaxis()->SetTitle("Amplitude[mV]");
      before_corr[ch]->GetXaxis()->CenterTitle();
      before_corr[ch]->GetYaxis()->SetTitle("Time[ps]");
      before_corr[ch]->GetYaxis()->CenterTitle();
      TProfile *htemp = before_corr[ch]->ProfileX();
      htemp->SetLineColor(kBlack);
      htemp->SetMarkerColor(kBlack);
      htemp->SetMarkerStyle(21);
      htemp->SetMarkerSize(0.6);
      c_work->Update();
      htemp->Draw("same");
      htemp->Fit("myfit", "Q", "QRsame", corr_amp_start[ch], corr_amp_stop[ch]);
      Double_t par[6]{};
      myfit->GetParameters(par);
      for (Int_t i = 0; i < valid_evt; i++) {
        deltaT[i] = deltaT[i] - myfit->Eval(peak[ch][i]);
      };

      after_corr[ch] =
          new TH2D(label_process + "_after", label_process + "_after", time_bin,
                   corr_amp_start[ch], corr_amp_stop[ch], peak_bin,
                   corr_time_sta, corr_time_end);

      for (Int_t i = 0; i < valid_evt; i++) {
        after_corr[ch]->Fill(peak[ch][i], deltaT[i]);
      }
      c_work->cd(ch + n_corr + 1);
      after_corr[ch]->Draw("colz");
      after_corr[ch]->GetXaxis()->SetTitle("Amplitude[mV]");
      after_corr[ch]->GetXaxis()->CenterTitle();
      after_corr[ch]->GetYaxis()->SetTitle("Time[ps]");
      after_corr[ch]->GetYaxis()->CenterTitle();
      c_work->Update();
    }
    // c_work->Print(label + "_process.pdf");
    for (Int_t ch = 0; ch < n_corr; ch++) {
      delete before_corr[ch];
      delete after_corr[ch];
    }
    draw_amp_and_time(n_corr, deltaT, peak, peak_bin, corr_amp_start,
                      corr_amp_stop, time_bin, corr_time_sta, corr_time_end,
                      label, "after", titlepart, targetID);
  }
  // c_work->Print(label + "_process.pdf]");

  // c_temp.Print(label + "_resolution_after.pdf]");
  // c_temp.Print(label + "_after.pdf]");

  delete c_work;
}

void slewing_t0_trg(TString label = "t0_RUN1", Double_t threshold = 0) {
  vector<vector<Double_t>> peak{2}; // mV
  vector<vector<Double_t>> time{2}; // ns->ps
  vector<Double_t> deltaT{};        // delta T ps

  TString filename = label + "_timing.root";
  TString treename = "time";
  TFile datafile(filename, "read");
  TTree *tree = datafile.Get<TTree>(treename);
  Double_t t0, trg1, trg2, amp1, amp2;
  tree->SetBranchAddress("t0", &t0);
  tree->SetBranchAddress("trg1", &trg1);
  tree->SetBranchAddress("trg2", &trg2);
  tree->SetBranchAddress("amp1", &amp1);
  tree->SetBranchAddress("amp2", &amp2);
  TF1 *f1 = nullptr;
  TF1 *f2 = nullptr;
  tree->SetBranchAddress("f1", &f1);
  tree->SetBranchAddress("f2", &f2);
  Long64_t entries = tree->GetEntries();
  for (Long64_t ev = 0; ev < entries; ev++) {
    tree->GetEntry(ev);
    Double_t t1 = 0;
    Double_t t2 = 0;
    if (threshold == 0) {
      t1 = trg1;
      t2 = trg2;
    } else if (amp1 >= threshold + 30 && amp2 >= threshold + 30) {
      t1 = f1->GetX(threshold * 4.096, 0, f1->GetMaximumX()) * sample_inte_ns;
      t2 = f2->GetX(threshold * 4.096, 0, f2->GetMaximumX()) * sample_inte_ns;
    } else {
      continue;
    }

    if (TMath::IsNaN(t1) || TMath::IsNaN(t2) || t1 == 0 || t2 == 0)
      continue;
    if (amp1 < 1500 || amp2 < 1500)
      continue;

    peak[0].push_back(amp1);
    time[0].push_back(t1 * 1000);
    peak[1].push_back(amp2);
    time[1].push_back(t2 * 1000);
  }
  Int_t valid_evt = time[0].size();
  for (Int_t i = 0; i < valid_evt; i++) {
    Double_t delta_t = (time[1][i] - time[0][i]) / 2;
    deltaT.push_back(delta_t);
  }

  label += "_time_" + TString::Itoa((Int_t)threshold, 10);
  filename_base = label;
  TString label_for_slewing = filename_base + TString::Format("_%d%d", 1, 2);
  TString title_for_slewing = TString::Format("(T%d-T%d)/2", 2, 1);
  Int_t tarID2[2]{1, 2};

  conduct_slewing(2, 3, deltaT, peak, cor_num, label_for_slewing,
                  title_for_slewing, tarID2);
  datafile.Close();
}

void slewing() {
  TH1::AddDirectory(false);
  gStyle->SetOptFit(1110);
  gStyle->SetOptStat(10);
  TGaxis::SetMaxDigits(3);

  sample_inte_ns = 0.4;
  cor_num = 5;
  peak_bin = 50;
  corr_amp_start[0] = 1500;
  corr_amp_start[1] = 1500;
  corr_amp_start[2] = 1500;
  corr_amp_start[3] = 1500;
  corr_amp_stop[0] = 3000;
  corr_amp_stop[1] = 2500;
  corr_amp_stop[2] = 2500;
  corr_amp_stop[3] = 2500;
  time_bin = 50;
  corr_time_sta_1st = 1000;
  corr_time_end_1st = 6000;
  corr_time_sta = -2000;
  corr_time_end = 2000;

  slewing_t0_trg("t0_RUN1", 0);
  cout << "Just before the end of the whole program" << endl;
}