#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <sstream>

#include <TStyle.h>
#include <TLegend.h>
#include <TFile.h>
#include <TTree.h>
#include <TBranch.h>
#include <TH1D.h>
#include <TCanvas.h>
#include <TLatex.h>
#include <TString.h>

//Variables definition
Float_t MEM_LL_sum_sg;
Float_t MEM_LL_sum_bg;

Float_t MEM_LL_sum_sg2;
Float_t MEM_LL_sum_bg2;

Float_t MEMdf;
Float_t MEMd2f;

Double_t MEMd;
Double_t MEMd2;

Double_t nomWeight_weight_mc;
Double_t nomWeight_weight_pu;
Double_t nomWeight_weight_elec;
Double_t nomWeight_weight_elec_trigger;
Double_t nomWeight_weight_muon;
Double_t nomWeight_weight_muon_trigger;
Double_t nomWeight_weight_btag;
Double_t nomWeight_weight_norm;
Double_t nomWeight_weight_jvt;

Double_t nomWeight_weight_mc2;
Double_t nomWeight_weight_pu2;
Double_t nomWeight_weight_elec2;
Double_t nomWeight_weight_elec_trigger2;
Double_t nomWeight_weight_muon2;
Double_t nomWeight_weight_muon_trigger2;
Double_t nomWeight_weight_btag2;
Double_t nomWeight_weight_norm2;
Double_t nomWeight_weight_jvt2;

void MEMd1 () {

  //Styling
  gStyle->SetOptStat(0);
  gStyle->SetLineWidth(1);
  gStyle->SetLegendBorderSize(0);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetPadGridX(0);
  gStyle->SetPadGridY(0);
  gStyle->SetPadTickX(1);
  gStyle->SetPadTickY(1);
  Double_t norm = 1;
  Bool_t isNDC = kTRUE;

  TFile *S_sgmodel = new TFile("sg_MEM-000000-tthbb-200-6ji4bi_0.root");
  TTree *S_sgtree = (TTree*)S_sgmodel->Get("tthbb_sg");
  S_sgtree->SetBranchAddress("MEM_LL_sum_sg",&MEM_LL_sum_sg);

  TFile *S_bgmodel = new TFile("bg_MEM-000000-tthbb-200-6ji4bi_0.root");
  TTree *S_bgtree = (TTree*)S_bgmodel->Get("tthbb_bg");
  S_bgtree->SetBranchAddress("MEM_LL_sum_bg",&MEM_LL_sum_bg);

  TFile *S_originalrootfile = new TFile("tthbb_200_TREE.root");
  TTree *treeS = (TTree*)S_originalrootfile->Get("tree");
  treeS->SetBranchAddress("nomWeight_weight_mc",&nomWeight_weight_mc);
	treeS->SetBranchAddress("nomWeight_weight_pu",&nomWeight_weight_pu);
	treeS->SetBranchAddress("nomWeight_weight_elec",&nomWeight_weight_elec);
	treeS->SetBranchAddress("nomWeight_weight_muon",&nomWeight_weight_muon);
	treeS->SetBranchAddress("nomWeight_weight_elec_trigger",&nomWeight_weight_elec_trigger);
	treeS->SetBranchAddress("nomWeight_weight_muon_trigger",&nomWeight_weight_muon_trigger);
	treeS->SetBranchAddress("nomWeight_weight_btag",&nomWeight_weight_btag);
	treeS->SetBranchAddress("nomWeight_weight_norm",&nomWeight_weight_norm);
  treeS->SetBranchAddress("nomWeight_weight_jvt",&nomWeight_weight_jvt);

  TFile *B_sgmodel = new TFile("sg_MEM-000000-ttbar-bg-6ji4bi_0.root");
  TTree *B_sgtree = (TTree*)B_sgmodel->Get("ttbar_sg");
  B_sgtree->SetBranchAddress("MEM_LL_sum_sg",&MEM_LL_sum_sg2);

  TFile *B_bgmodel = new TFile("bg_MEM-000000-ttbar-bg-6ji4bi_0.root");
  TTree *B_bgtree = (TTree*)B_bgmodel->Get("ttbar_bg");
  B_bgtree->SetBranchAddress("MEM_LL_sum_bg",&MEM_LL_sum_bg2);

  TFile *B_originalrootfile = new TFile("ttbar_TREE.root");
  TTree *treeB = (TTree*)B_originalrootfile->Get("tree");
  treeB->SetBranchAddress("nomWeight_weight_mc",&nomWeight_weight_mc2);
  treeB->SetBranchAddress("nomWeight_weight_pu",&nomWeight_weight_pu2);
  treeB->SetBranchAddress("nomWeight_weight_elec",&nomWeight_weight_elec2);
  treeB->SetBranchAddress("nomWeight_weight_muon",&nomWeight_weight_muon2);
  treeB->SetBranchAddress("nomWeight_weight_elec_trigger",&nomWeight_weight_elec_trigger2);
  treeB->SetBranchAddress("nomWeight_weight_muon_trigger",&nomWeight_weight_muon_trigger2);
  treeB->SetBranchAddress("nomWeight_weight_btag",&nomWeight_weight_btag2);
  treeB->SetBranchAddress("nomWeight_weight_norm",&nomWeight_weight_norm2);
  treeB->SetBranchAddress("nomWeight_weight_jvt",&nomWeight_weight_jvt2);

  TH1D *hs  = new TH1D("hs","Matrix Element Discriminant;MEM_{ D1};Events (norm. to unity)",50,0,1);
  hs->SetTitleOffset(1.5,"Y");
  hs->SetFillColor(kRed);
  hs->SetLineWidth(2);
  hs->SetFillStyle(3354);
  hs->SetLineColor(kRed);

  TH1D *hb  = new TH1D("hb","Matrix Element Discriminant;MEM_{ D1};Events (norm. to unity)",50,0,1);
  hb->SetTitleOffset(1.5,"Y");
  hb->SetFillColor(kAzure-8);
  hb->SetLineWidth(2);
  hb->SetFillStyle(1001);
  hb->SetLineColor(kBlue);

  unsigned int NEN = 1000;

  for (unsigned int k = 0; k < NEN; k++) {
    S_sgtree->GetEntry(k);
    S_bgtree->GetEntry(k);
    treeS->GetEntry(k);
    //MEMdf = 1/(1 + exp(- (MEM_LL_sum_sg - MEM_LL_sum_bg)));
    MEMdf = pow(MEM_LL_sum_sg - MEM_LL_sum_bg,10);
    MEMd = (Double_t)MEMdf;
    hs->Fill(MEMd,nomWeight_weight_btag*nomWeight_weight_elec*nomWeight_weight_elec_trigger*nomWeight_weight_jvt*nomWeight_weight_mc*nomWeight_weight_muon*nomWeight_weight_muon_trigger*nomWeight_weight_norm*nomWeight_weight_pu);
  }

  for (unsigned int k = 0; k < NEN; k++) {
    B_sgtree->GetEntry(k);
    B_bgtree->GetEntry(k);
    treeB->GetEntry(k);
    //MEMd2f = 1/(1 + exp(- (MEM_LL_sum_sg2 - MEM_LL_sum_bg2)));
    MEMd2f = pow(MEM_LL_sum_sg2 - MEM_LL_sum_bg2,10);
    MEMd2 = (Double_t)MEMd2f;
    hb->Fill(MEMd2,nomWeight_weight_btag2*nomWeight_weight_elec2*nomWeight_weight_elec_trigger2*nomWeight_weight_jvt2*nomWeight_weight_mc2*nomWeight_weight_muon2*nomWeight_weight_muon_trigger2*nomWeight_weight_norm2*nomWeight_weight_pu2);
  }

  TCanvas *c3;
  c3 = NULL;
  c3 = new TCanvas("c3","c3");

  hs->Scale(1/hs->GetSumOfWeights());
  hb->Scale(1/hb->GetSumOfWeights());

  double ymax;
  if (hb->GetMaximum() > hs->GetMaximum()) ymax = 1.4*hb->GetMaximum();
  else ymax = 1.4*hs->GetMaximum();
  hb->SetMaximum(ymax);

  hb->Draw("");
  hs->Draw("same");
  TLegend *lj = new TLegend(.52, .70,.87,.85);
  lj->AddEntry("hb","Total background","f");
  lj->AddEntry("hs","ttH(#rightarrow bb) (m_{H} = 200 GeV)","f");
  lj->SetEntrySeparation(0);
  lj->Draw();
  TLatex t(.14, .76, "#splitline{#it{ATLAS} #bf{Internal}}{#scale[0.8]{#bf{6j-in. 4b-in.}}}");
  t.SetNDC();
  t.DrawClone();
  c3->Print("tthbb_200.pdf");
  delete c3;
  c3 = NULL;
}