#include "TFile.h"
#include "TCanvas.h"
#include "TH1.h"

#define SBSCALE_DPL 0.5
#define SBSCALE_DST 0.1
#define REBIN

TFile *dpl;
TFile *dst;

TCanvas *dpl_canvas;
TCanvas *dst_canvas;

TCanvas *dpl_sub_canvas;
TCanvas *dst_sub_canvas;

TCanvas *dpl_final_canvas;
TCanvas *dst_final_canvas;

TH1F *dpl_srs,*dpl_sbrs,*dpl_sws,*dpl_sbws;
TH1F *dst_srs,*dst_sbrs,*dst_sws,*dst_sbws;

TH1F *dpl_srs_sws,*dpl_sbrs_sbws,dpl_srs_sbrs,dpl_sws_sbws,*dpl_final,*dpl_final_scaled;
TH1F *dst_srs_sws,*dst_sbrs_sbws,dst_srs_sbrs,dst_sws_sbws,*dst_final,*dst_final_scaled;

void subtract()
{
  dpl=new TFile("bdSS_dkpp.root","read");
  dst=new TFile("bdSS_dkps.root","read");

  dpl_canvas=new TCanvas("dpl_canvas","D^{+} histograms");
  dst_canvas=new TCanvas("dst_canvas","D^{*} histograms");

  dpl_canvas->Divide(2,2);
  dst_canvas->Divide(2,2);

  dpl->cd();
  dpl_srs =(TH1F*)dpl->FindObjectAny("hDD0MassDPi1RS_3")->Clone("dpl_srs");
  dpl_sbrs=(TH1F*)dpl->FindObjectAny("hDD0MassDPi1WS_3")->Clone("dpl_sbrs");
  dpl_sws =(TH1F*)dpl->FindObjectAny("hDD0MassDPi2RS_3")->Clone("dpl_sws");
  dpl_sbws=(TH1F*)dpl->FindObjectAny("hDD0MassDPi2WS_3")->Clone("dpl_sbws");
  // dpl_sbws=(TH1F*)dpl->FindObjectAny("hDD0MassDPi2RS_3")->Clone("dpl_sbws");


  dst->cd();
  dst_srs =(TH1F*)dst->FindObjectAny("hDD0MassDPi1RS_8")->Clone("dst_srs");
  dst_sbrs=(TH1F*)dst->FindObjectAny("hDD0MassDPi1WS_8")->Clone("dst_sbrs");
  dst_sws =(TH1F*)dst->FindObjectAny("hDD0MassDPi2RS_8")->Clone("dst_sws");
  dst_sbws=(TH1F*)dst->FindObjectAny("hDD0MassDPi2WS_8")->Clone("dst_sbws");
  // dst_sbws=(TH1F*)dst->FindObjectAny("hDD0MassDPi2RS_8")->Clone("dst_sbws");

  dpl_srs ->Sumw2();
  dpl_sbrs->Sumw2();
  dpl_sws ->Sumw2();
  dpl_sbws->Sumw2();

  dst_srs ->Sumw2();
  dst_sbrs->Sumw2();
  dst_sws ->Sumw2();
  dst_sbws->Sumw2();

  dpl_srs ->SetTitle("D^+ Signal Right Sign  ");
  dpl_sbrs->SetTitle("D^+ Sideband Right Sign");
  dpl_sws ->SetTitle("D^+ Signal Wrong Sign  ");
  dpl_sbws->SetTitle("D^+ Sideband Wrong Sign");

  dst_srs ->SetTitle("D^* Signal Right Sign  ");
  dst_sbrs->SetTitle("D^* Sideband Right Sign");
  dst_sws ->SetTitle("D^* Signal Wrong Sign  ");
  dst_sbws->SetTitle("D^* Sideband Wrong Sign");

#ifdef REBIN
  dpl_srs ->Rebin(4);
  dpl_sbrs->Rebin(4);
  dpl_sws ->Rebin(4);
  dpl_sbws->Rebin(4);

  dst_srs ->Rebin(4);
  dst_sbrs->Rebin(4);
  dst_sws ->Rebin(4);
  dst_sbws->Rebin(4);
#endif

  dpl_canvas->cd(1); dpl_srs ->DrawCopy();
  dpl_canvas->cd(2); dpl_sws ->DrawCopy();
  dpl_canvas->cd(3); dpl_sbrs->DrawCopy();
  dpl_canvas->cd(4); dpl_sbws->DrawCopy();

  dst_canvas->cd(1); dst_srs ->DrawCopy();
  dst_canvas->cd(2); dst_sws ->DrawCopy();
  dst_canvas->cd(3); dst_sbrs->DrawCopy();
  dst_canvas->cd(4); dst_sbws->DrawCopy();


  dpl_srs_sws=(TH1F*)dpl_srs->Clone("dpl_srs_sws");
  dpl_srs_sws->Add(dpl_sws,-1.0);
  dpl_sbrs_sbws=(TH1F*)dpl_sbrs->Clone("dpl_sbrs_sbws");
  dpl_sbrs_sbws->Add(dpl_sbws,-1.0);
  dpl_final=(TH1F*)dpl_srs_sws->Clone("dpl_final");
  dpl_final->Add(dpl_sbrs_sbws,-1.0*SBSCALE_DPL);
  dpl_srs_sbrs=(TH1F*) dpl_srs->Clone("dpl_srs_sbrs");
  dpl_srs_sbrs->Add(dpl_sbrs,-1.0*SBSCALE_DPL);
  dpl_sws_sbws=(TH1F*)dpl_sws->Clone("dpl_sws_sbws");
  dpl_sws_sbws->Add(dpl_sbws,-1.0*SBSCALE_DPL);

  dpl_srs_sws  ->SetTitle("D^+ Signal RS-WS");
  dpl_sbrs_sbws->SetTitle("D^+ Sideband RS-WS");
  dpl_srs_sbrs->SetTitle("D^+ Signal RS - a*Sideband RS");
  dpl_sws_sbws->SetTitle("D^+ Signal WS - a*Sideband WS");

  dst_srs_sws=(TH1F*)dst_srs->Clone("dst_srs_sws");
  dst_srs_sws->Add(dst_sws,-1.0);
  dst_sbrs_sbws=(TH1F*)dst_sbrs->Clone("dst_sbrs_sbws");
  dst_sbrs_sbws->Add(dst_sbws,-1.0);
  dst_final=(TH1F*)dst_srs_sws->Clone("dst_final");
  dst_final->Add(dst_sbrs_sbws,-1.0*SBSCALE_DST);
  dst_srs_sbrs=(TH1F*)dst_srs->Clone("dst_srs_sbrs");
  dst_srs_sbrs->Add(dst_sbrs,-1.0*SBSCALE_DST);
  dst_sws_sbws=(TH1F*)dst_sws->Clone("dst_sws_sbws");
  dst_sws_sbws->Add(dst_sbws,-1.0*SBSCALE_DST);

  dst_srs_sws  ->SetTitle("D^* Signal RS-WS");
  dst_sbrs_sbws->SetTitle("D^* Sideband RS-WS");
  dst_srs_sbrs->SetTitle("D^* Signal RS - a*Sideband RS");
  dst_sws_sbws->SetTitle("D^* Signal WS - a*Sideband WS");

  dpl_sub_canvas=new TCanvas("dpl_sub_canvas","D^{+} histograms");
  dst_sub_canvas=new TCanvas("dst_sub_canvas","D^{*} histograms");

  dpl_sub_canvas->Divide(2,2);
  dpl_sub_canvas->cd(1); dpl_srs_sbrs ->Draw();
  dpl_sub_canvas->cd(2); dpl_sws_sbws ->Draw();
  dpl_sub_canvas->cd(3); dpl_final    ->Draw();

  dst_sub_canvas->Divide(2,2);
  dst_sub_canvas->cd(1); dst_srs_sbrs ->Draw();
  dst_sub_canvas->cd(2); dst_sws_sbws ->Draw();
  dst_sub_canvas->cd(3); dst_final    ->Draw();

  dpl_final_canvas=new TCanvas("dpl_final_canvas","D^{+} scaled histograms");
  dst_final_canvas=new TCanvas("dst_final_canvas","D^{*} scaled histograms");

  dst_final_scaled=(TH1F*)dst_final->Clone("dst_final_scaled");
  dst_final_scaled->Scale(76.75/dst_final->Integral());

  dpl_final_scaled=(TH1F*)dpl_final->Clone("dpl_final_scaled");
  dpl_final_scaled->Scale(345.04/dpl_final->Integral());

  dpl_final->SetTitle("D^+ Final Subtraction");
  dst_final->SetTitle("D^* Final Subtraction");

  dpl_final_scaled->SetTitle("D^+ Final Subtraction, Scaled to our integral");
  dst_final_scaled->SetTitle("D^* Final Subtraction, Scaled to our integral");


  dpl_final_canvas->Divide(2,1);
  dpl_final_canvas->cd(1); dpl_final->Draw();
  dpl_final_scaled->SetMinimum(-10.0); dpl_final_scaled->SetMaximum(35.0);
  dpl_final_canvas->cd(2); dpl_final_scaled->Draw();

  dst_final_canvas->Divide(2,1);
  dst_final_canvas->cd(1); dst_final->Draw();
  dst_final_scaled->SetMinimum(-5.0); dst_final_scaled->SetMaximum(22.0);
  dst_final_canvas->cd(2); dst_final_scaled->Draw();
}