#include <iostream>
#include <fstream>

gROOT->Reset();

// Setup file including parameters
#include "plotMk_3D_setup.h"

void plotMk_3D_ProjZ_D0KPi() {

  // Root macro to produce a projection plot of 3D histogram
  //###########################################################
  // Change these values to meet your condition
  //###########################################################

  TString outPlotName("plot_"+oNameBase+"_"+varz+"_D0KPi.eps");
  TString histTitle(varz+"_D0KPi");

  // Set max value of plot
  int min = zmin;
  int max = zmax;
  //setMax = true;
  //int min = -100;
  //int max = 40000;

  //################################################
  cout << "start";
  // Open input files and get 3D histos
  TFile *f[nFiles_D0KPi];
  TH3D * inHist[nFiles_D0KPi];
  for (int l=0; l<nFiles_D0KPi; l++) {
    f[l] = new TFile(inFile_D0KPi[l]);
    //f[l]->ls();
    inHist[l] = (TH3D *)f[l]->Get("sub;1");
  }

  // Readout 3D histos into arrays
  cout << "reading histos into arrays" << endl;
  double num_3D[nFiles_D0KPi][nbinsx][nbinsy][nbinsz];
  double err_3D[nFiles_D0KPi][nbinsx][nbinsy][nbinsz];
  for (int l=0; l<nFiles_D0KPi; l++) {
    for (int i=0; i<nbinsx; i++) {
      for (int j=0; j<nbinsy; j++) {
	for (int k=0; k<nbinsz; k++) {
	  num_3D[l][i][j][k] = inHist[l]->GetBinContent(i+1,j+1,k+1);
	  err_3D[l][i][j][k] = inHist[l]->GetBinError(i+1,j+1,k+1);
	}
      }
    }
  }

  // 1D arrays to be projected
  cout << "define summed 1D arrays" << endl;
  double num_1D[nFiles_D0KPi][nbinsz];
  double var_1D[nFiles_D0KPi][nbinsz];
  double err_1D[nFiles_D0KPi][nbinsz];
  for (int l=0; l<nFiles_D0KPi; l++) {
    for (int j=0; j<nbinsz; j++) {
      num_1D[l][j] = 0;
      var_1D[l][j] = 0;
      err_1D[l][j] = 0;
    }
  }
  
  for (int l=0; l<nFiles_D0KPi; l++) {
    for (int k=0; k<nbinsz; k++) {
      for (int i=0; i<nbinsx; i++) {
	for (int j=0; j<nbinsy; j++) {
	  num_1D[l][k] += c_D0KPi[l] * num_3D[l][i][j][k];
	  var_1D[l][k] += c_D0KPi[l]*c_D0KPi[l]*err_3D[l][i][j][k]*err_3D[l][i][j][k];
	}
      }
      err_1D[l][k] = sqrt(var_1D[l][k]);
    }
  }

  // Sum over componemts
  cout << "getting combined array" << endl;
  double num_OnPeak[nbinsz];
  double num_OffPeak[nbinsz];
  double num_Dlnu[nbinsz];
  double num_Dstar[nbinsz];
  double num_Ddstar[nbinsz];
  double num_DPi[nbinsz];
  double num_DiffB[nbinsz];
  double num_CascL[nbinsz];
  double num_LMisID[nbinsz];
  double num_Comb[nbinsz];
  double num_RealD[nbinsz];
  double num_Cont[nbinsz];
  double var_OnPeak[nbinsz];
  double var_OffPeak[nbinsz];
  double var_Dlnu[nbinsz];
  double var_Dstar[nbinsz];
  double var_Ddstar[nbinsz];
  double var_DPi[nbinsz];
  double var_DiffB[nbinsz];
  double var_CascL[nbinsz];
  double var_LMisID[nbinsz];
  double var_Comb[nbinsz];
  double var_RealD[nbinsz];
  double var_Cont[nbinsz];
  double err_OnPeak[nbinsz];
  double err_OffPeak[nbinsz];
  double err_Dlnu[nbinsz];
  double err_Dstar[nbinsz];
  double err_Ddstar[nbinsz];
  double err_DPi[nbinsz];
  double err_DiffB[nbinsz];
  double err_CascL[nbinsz];
  double err_LMisID[nbinsz];
  double err_Comb[nbinsz];
  double err_RealD[nbinsz];
  double err_Cont[nbinsz];
  for (int j=0; j<nbinsz; j++) {
    num_OnPeak[j] = num_1D[0][j];
    var_OnPeak[j] = var_1D[0][j];
    err_OnPeak[j] = err_1D[0][j];
    if (useOneHist) {
      num_OffPeak[j] = num_1D[1][j];
      var_OffPeak[j] = var_1D[1][j];
    } else {
    num_Dlnu[j] 
      = (num_1D[29][j] - num_1D[30][j]*rho2 + num_1D[31][j]*rho2*rho2)*RN;
    var_Dlnu[j] 
      = (var_1D[29][j] + var_1D[30][j]*rho2*rho2
	 + var_1D[31][j]*rho2*rho2*rho2*rho2)*RN;
    }
    err_Dlnu[j] = sqrt(var_Dlnu[j]);
    num_Dlnu[j] = num_1D[2][j];
    var_Dlnu[j] = var_1D[2][j];
    err_Dlnu[j] = err_1D[2][j];
    num_Dstar[j] = num_1D[27][j] + num_1D[28][j];
    var_Dstar[j] = var_1D[27][j] + var_1D[28][j];
    err_Dstar[j] = sqrt(var_Dstar[j]);
    num_Ddstar[j] 
      = num_1D[3][j] + num_1D[4][j] + num_1D[5][j] + num_1D[6][j]
      + num_1D[7][j] + num_1D[8][j] + num_1D[9][j] + num_1D[10][j];
    var_Ddstar[j] 
      = var_1D[3][j] + var_1D[4][j] + var_1D[5][j] + var_1D[6][j]
      + var_1D[7][j] + var_1D[8][j] + var_1D[9][j] + var_1D[10][j];
    err_Ddstar[j] = sqrt(var_Ddstar[j]);
    num_DPi[j] 
      = num_1D[11][j] + num_1D[12][j] + num_1D[13][j] + num_1D[14][j];
    var_DPi[j] 
      = var_1D[11][j] + var_1D[12][j] + var_1D[13][j] + var_1D[14][j];
    err_DPi[j] = sqrt(var_DPi[j]);
    num_DiffB[j] = num_1D[15][j] + num_1D[16][j];
    var_DiffB[j] = var_1D[15][j] + var_1D[16][j];
    err_DiffB[j] = sqrt(var_DiffB[j]);
    num_CascL[j] = num_1D[17][j] + num_1D[18][j];
    var_CascL[j] = var_1D[17][j] + var_1D[18][j];
    err_CascL[j] = sqrt(var_CascL[j]);
    num_LMisID[j] = num_1D[19][j] + num_1D[20][j];
    var_LMisID[j] = var_1D[19][j] + var_1D[20][j];
    err_LMisID[j] = sqrt(var_LMisID[j]);
    num_Comb[j] = num_1D[21][j] + num_1D[22][j];
    var_Comb[j] = var_1D[21][j] + var_1D[22][j];
    err_Comb[j] = sqrt(var_Comb[j]);
    num_RealD[j] = num_1D[23][j];
    var_RealD[j] = var_1D[23][j];
    err_RealD[j] = err_1D[23][j];
    num_Cont[j] = num_1D[24][j] + num_1D[25][j] + num_1D[26][j];
    var_Cont[j] = var_1D[24][j] + var_1D[25][j] + var_1D[26][j];
    err_Cont[j] = sqrt(var_Cont[j]);
  }

  //Define 1D histograms and fill them
  TH1D *OnPeak = new TH1D("OnPeak",histTitle,nbinsz,zbins);
  TH1D *OffPeak = new TH1D("OffPeak",histTitle,nbinsz,zbins);
  TH1D *Dlnu = new TH1D("Dlnu",histTitle,nbinsz,zbins);
  TH1D *Dstar = new TH1D("Dstar",histTitle,nbinsz,zbins);
  TH1D *Ddstar = new TH1D("Ddstar",histTitle,nbinsz,zbins);
  TH1D *DPi = new TH1D("DPi",histTitle,nbinsz,zbins);
  TH1D *DiffB = new TH1D("DiffB",histTitle,nbinsz,zbins);
  TH1D *CascL = new TH1D("CascL",histTitle,nbinsz,zbins);
  TH1D *LMisID = new TH1D("LMisID",histTitle,nbinsz,zbins);
  TH1D *Comb = new TH1D("Comb",histTitle,nbinsz,zbins);
  TH1D *RealD = new TH1D("RealD",histTitle,nbinsz,zbins);
  TH1D *Cont = new TH1D("Cont",histTitle,nbinsz,zbins);

  for (int j=0; j<nbinsz; j++) {
    OnPeak->SetBinContent(j+1,num_OnPeak[j]);
    OnPeak->SetBinError(j+1,err_OnPeak[j]);
    OffPeak->SetBinContent(j+1,num_OffPeak[j]);
    OffPeak->SetBinError(j+1,err_OffPeak[j]);
    Dlnu->SetBinContent(j+1,num_Dlnu[j]);
    Dlnu->SetBinError(j+1,err_Dlnu[j]);
    Dstar->SetBinContent(j+1,num_Dstar[j]);
    Dstar->SetBinError(j+1,err_Dstar[j]);
    Ddstar->SetBinContent(j+1,num_Ddstar[j]);
    Ddstar->SetBinError(j+1,err_Ddstar[j]);
    DPi->SetBinContent(j+1,num_DPi[j]);
    DPi->SetBinError(j+1,err_DPi[j]);
    DiffB->SetBinContent(j+1,num_DiffB[j]);
    DiffB->SetBinError(j+1,err_DiffB[j]);
    CascL->SetBinContent(j+1,num_CascL[j]);
    CascL->SetBinError(j+1,err_CascL[j]);
    LMisID->SetBinContent(j+1,num_LMisID[j]);
    LMisID->SetBinError(j+1,err_LMisID[j]);
    Comb->SetBinContent(j+1,num_Comb[j]);
    Comb->SetBinError(j+1,err_Comb[j]);
    RealD->SetBinContent(j+1,num_RealD[j]);
    RealD->SetBinError(j+1,err_RealD[j]);
    Cont->SetBinContent(j+1,num_Cont[j]);
    Cont->SetBinError(j+1,err_Cont[j]);
  }

  //#############################################################
  cout << "strt plotting\n";
  cout << "setting marker and color\n";
  OnPeak->SetMarkerStyle(29);  
  OffPeak->SetFillColor(12);  
  Dlnu->SetFillColor(2);
  Dstar->SetFillColor(3);
  Ddstar->SetFillColor(4);
  DPi->SetFillColor(5);
  DiffB->SetFillColor(6);
  CascL->SetFillColor(7) ;
  LMisID->SetFillColor(8) ;
  Comb->SetFillColor(21);
  RealD->SetFillColor(9);
  Cont->SetFillColor(12);

  cout << "stacking\n";
  THStack hAll("hAll",histTitle);
  hAll.Add(Cont);
  hAll.Add(RealD);
  hAll.Add(Comb);
  hAll.Add(LMisID);
  hAll.Add(CascL);
  hAll.Add(DiffB);
  hAll.Add(DPi);
  hAll.Add(Ddstar);
  hAll.Add(Dstar);
  hAll.Add(Dlnu);

  THStack hOffPeak("hOffPeak",histTitle);
  hOffPeak.Add(OffPeak);
  hOffPeak.Add(Comb);
  hOffPeak.Add(LMisID);
  hOffPeak.Add(CascL);
  hOffPeak.Add(DiffB);
  hOffPeak.Add(DPi);
  hOffPeak.Add(Ddstar);
  hOffPeak.Add(Dstar);
  hOffPeak.Add(Dlnu);

  // Create a canvas and plot
  cout << "plotting\n";
  TCanvas * cx = new TCanvas("cx",histTitle,0,0,800,400);
  cx->Divide(2,1);
  gStyle->SetOptStat(000000);

  cx->cd(1);

  hAll.SetMinimum(min);
  if (setMax) {
    hAll.SetMaximum(max);
  } else {
    hAll.SetMaximum(OnPeak->GetMaximum()*1.2);
  }
  hAll.Draw("hist");
  OnPeak->Draw("e1psames");

  TLegend *legendx1 = new TLegend(0.70,0.50,0.90,0.92);
  legendx1->SetTextSize(.02);
  legendx1->SetFillColor(0);
  TLegendEntry *plex1;
  plex1 = legendx1->AddEntry(OnPeak, "OnPeak data","P");
  plex1 = legendx1->AddEntry(Dlnu,"B->Dlnu","F");
  plex1 = legendx1->AddEntry(Dstar,"B->D*lnu","F");
  plex1 = legendx1->AddEntry(Ddstar,"B->D**lnu","F");
  plex1 = legendx1->AddEntry(DPi,"B->D/D* Pilnu","F");
  plex1 = legendx1->AddEntry(DiffB,"D l from different B","F");
  plex1 = legendx1->AddEntry(CascL,"B->D->l","F");
  plex1 = legendx1->AddEntry(LMisID,"lepton Miss ID","F");
  plex1 = legendx1->AddEntry(Comb,"combinatrial BB","F");
  plex1 = legendx1->AddEntry(RealD,"ccbar with real D","F");
  plex1 = legendx1->AddEntry(Cont,"continuum","F");

  legendx1->Draw();

  cx->cd(2);

  hOffPeak.SetMinimum(min);
  if (setMax) {
    hOffPeak.SetMaximum(max);
  } else {
    hOffPeak.SetMaximum(OnPeak->GetMaximum()*1.2);
  }
  hOffPeak.Draw("hist");
  OnPeak->Draw("e1psames");

  TLegend *legendx2 = new TLegend(0.70,0.50,0.90,0.92);
  legendx2->SetTextSize(.02);
  legendx2->SetFillColor(0);
  TLegendEntry *plex2;
  plex2 = legendx2->AddEntry(OnPeak, "OnPeak data","P");
  plex2 = legendx2->AddEntry(Dlnu,"B->Dlnu","F");
  plex2 = legendx2->AddEntry(Dstar,"B->D*lnu","F");
  plex2 = legendx2->AddEntry(Ddstar,"B->D**lnu","F");
  plex2 = legendx2->AddEntry(DPi,"B->D/D* Pilnu","F");
  plex2 = legendx2->AddEntry(DiffB,"D l from different B","F");
  plex2 = legendx2->AddEntry(CascL,"B->D->l","F");
  plex2 = legendx2->AddEntry(LMisID,"lepton Miss ID","F");
  plex2 = legendx2->AddEntry(Comb,"combinatrial BB","F");
  plex2 = legendx2->AddEntry(OffPeak,"OffPeak data","F");

  legendx2->Draw();

  cx->Print(outPlotName);

}