// To run:
// root -l -q -b plotTiming.C
#include <fstream>
#include <string>
#include <sstream>
#include <iostream>
#include <ostream>
#include <algorithm>
#include <iomanip>
#include <TLorentzVector.h>
#include <TLorentzRotation.h>
#include <TVector3.h>
#include <TRotation.h>
#include <TMath.h>
#include <math.h>
#include <TFile.h>
#include <TChain.h>
#include <TTree.h>
#include <vector>
#include <set>
#include <TH1F.h>
#include <TH2F.h>
#include <TRandom3.h>
#include <exception>
#include <cmath>
#include <TLorentzVector.h>
#include <TStyle.h>
#include "TROOT.h"
#include "TH1F.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TMVA/MethodBase.h"
#include "TMVA/Tools.h"
using namespace std;


//#include "TCanvas.h"
//#include "TStyle.h"
//#include "TH1.h"
//#include <TChain.h>
//#include "TGaxis.h"
//#include "TRandom.h"
//#include "TTree.h"
//#include "TBranch.h"
//#include "TFile.h"
//#include "TF1.h"
//#include <sstream>

void overlay_hist(){

  //TChain for ttH
  std::ifstream file ("list_ttH.txt");
  std::string str;
  std::string file_contents;
  std::string treename;
  
  TChain *ch = new TChain ("nominal");
  while (std::getline(file,str))
    {
      file_contents = str;
      const char* filename = file_contents.c_str();
      std::cout << "name: " << filename << "\n";
      TFile* f = TFile::Open(filename);
      TTree* tree = (TTree*)f->Get("nominal");
      treename = file_contents + "/nominal";
      const char* name = treename.c_str();
      ch->AddFile(name);
    }
  
 //TChain for ttbar
  std::ifstream file1 ("list_ttbar.txt");
  std::string str1;
  std::string file_contents1;
  std::string treename1;

  TChain *ch2 = new TChain ("nominal");
  while (std::getline(file1,str1))
    {
      file_contents1 = str1;
      const char* filename1 = file_contents1.c_str();
      std::cout << "name: " << filename1 << "\n";
      TFile* f1 = TFile::Open(filename1);
      TTree* tree1 = (TTree*)f1->Get("nominal");
      treename1 = file_contents1 + "/nominal";
      const char* name1 = treename1.c_str();
      ch2->AddFile(name1);
    }

    
  TFile *newfile = new TFile("/gs/project/atlas-mcgill/ssaha/AnalysisTop-2.4.17_NonPrompt/merged_ttbar_R_new.root");
  TTree *ttHTree = new TTree("tth_final", "ttH_weighted_cs");
  TTree *ttbarTree = new TTree("ttbar_final", "ttbar_weighted_cs");

  
  ttHTree = (TTree*)ch->CloneTree(0);
  ttbarTree = (TTree*)ch2->CloneTree(0);

  
  Double_t cross_section;
  Int_t mcChannelNumber;
  mcChannelNumber = 0;
  cross_section = 0;

 ch->SetBranchAddress("mcChannelNumber", &mcChannelNumber);
 ch2->SetBranchAddress("mcChannelNumber", &mcChannelNumber);

  ttHTree->Branch("cross_section", &cross_section);
  ttbarTree->Branch("cross_section", &cross_section);


  //tth loop
  for(Int_t ievt =0; ievt < ch->GetEntries(); ievt++){
    ch->GetEntry(ievt);
    if(mcChannelNumber == 343365) cross_section = 0.05343;
    if(mcChannelNumber == 343366) cross_section = 0.22276;
    if(mcChannelNumber == 343367) cross_section = 0.23082;
    ttHTree->Fill();
  }
  
 //ttbar loop
  for(Int_t ievt =0; ievt < ch2->GetEntries(); ievt++){
    ch2->GetEntry(ievt);
    cross_section = 377.9932 * 1.1949;
    ttbarTree->Fill();
  }



  // TTree *t1 = (TTree *)f1->Get("nominal");
  TH1F *h1 = new TH1F ("ttbar","log_likelihood", 100, -90, -25);
  ttbarTree->Project("ttbar","klfitter_logLikelihood","weight_mc");
  h1->Scale(1/h1->Integral());

  
  
  // TTree  *t2 = (TTree *)f2->Get("nominal"); //loads the Tree named "nominal"
  TH1F *h2 = new TH1F ("ttH","log_likelihood", 100, -90, -25);
  ttHTree->Project("ttH","klfitter_logLikelihood","weight_mc");  //Draws the branch named "kl.." to a histogram "h22"
  h2->Scale(1/h2->Integral());//normalizes the histogram.
  
  //Overlap area calculation
  
  TH1F *h31 = new TH1F("h31","subtracted_histo",100,-90,-25);
  h31->Add(h1,h2,1.,-1.);
  //TCanvas *c1 = new TCanvas("h31","h31");
  //h31->Draw();
  
  //loop manually over the bins to calculate absolute integral
  
  TH1F *h3 = new TH1F(*h1);
  h3->SetNameTitle("h3", "min(h1, h2)");
  h3->Reset("M");
  for (int i = 1; i <= h3->GetNbinsX(); i++)
    h3->Fill(h3->GetBinCenter(i),
	     ((h1->GetBinContent(i) < h2->GetBinContent(i)) ? h2->GetBinContent(i)- h1->GetBinContent(i):h1->GetBinContent(i)- h2->GetBinContent(i)));
  std::cout << h3->Integral() << std::endl;
  std::cout <<"width_integral"<<" "<< h3->Integral("width") << std::endl;//By default, the integral is computed as the sum of bin contents in the range. if option "width" is specified, the integral is the sum of the bin contents multiplied by the bin width in x.
  //TCanvas *c2 = new TCanvas("h3","h3");
  //h3->GetYaxis()->SetRangeUser(-0.3,0.3);
  //h3->Draw();
  
TLatex *l =
  new TLatex(0.5, 0.5,
             TString::Format("#int%s = %g", h3->GetName(), h3->Integral()));
 l->SetNDC(kTRUE);
 l->Draw();
 
 
 Double_t integral = 0;
 integral= TMVA::gTools().GetSeparation(h1,h2);
 std::cout<<"TMVA"<<integral<<std::endl;
 
 
 //Draw the plots and save them. 
 
 TCanvas *c0 = new TCanvas("loglikelihood","loglikelihood");
 h1->Draw();
 h1->SetLineColor(kRed);
 gPad->Update();
 TPaveStats* st1 = (TPaveStats*) h1->FindObject("stats");
 
 h2->Draw();
 gPad->Update();
 TPaveStats* st2 = (TPaveStats*) h2->FindObject("stats");
 st1->SetX1NDC(0.1);
 st1->SetX2NDC(0.3);
 st1->SetY1NDC(0.7);
 st1->SetY2NDC(0.9);
 st2->SetX1NDC(0.3);
 st2->SetX2NDC(0.5);
 st2->SetY1NDC(0.7);
 st2->SetY2NDC(0.9);
 st1->SetTextColor(kRed);
 st2->SetTextColor(kBlue);
 
 h2->SetLineColor(kBlue);
 h1->Draw("");
 h2->Draw("SAME");
 
 c0->Update();
 c0->SaveAs("loglikelihood_ttbar_ttH.gif");
 
}