#include <iostream>
#include <fstream>

gROOT->Reset();

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

void plotMk_3D_ProjZ_bin_D0KPi() {

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

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

  int min = zminBin;
  int max = zmaxBin;
  //requireSetMax = true;
  //int min = -100;
  //int max = 4000;

  int nTotBins = nbinsx * nbinsy;
  enum {nbinsTot = 50};

  //################################################
  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);
	}
      }
    }
  }

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


  //Define histograms and fill them
  cout << "defining histograms to plot" << endl;
  TH1D *OnPeak[nbinsTot]; 
  TH1D *OffPeak[nbinsTot]; 
  TH1D *Dlnu[nbinsTot]; 
  TH1D *Dstar[nbinsTot]; 
  TH1D *Ddstar [nbinsTot];
  TH1D *DPi[nbinsTot]; 
  TH1D *DiffB[nbinsTot]; 
  TH1D *CascL[nbinsTot]; 
  TH1D *LMisID[nbinsTot]; 
  TH1D *Comb[nbinsTot];
  TH1D *RealD[nbinsTot]; 
  TH1D *Cont[nbinsTot]; 
  for (int l=0; l<nbinsTot; l++) {
    OnPeak[l] = new TH1D(TString("OnPeak")+l,histTitle,nbinsz,zbins);
    OffPeak[l] = new TH1D(TString("OffPeak")+l,histTitle,nbinsz,zbins);
    Dlnu[l] = new TH1D(TString("Dlnu")+l,histTitle,nbinsz,zbins);
    Dstar[l] = new TH1D(TString("Dstar")+l,histTitle,nbinsz,zbins);
    Ddstar[l] = new TH1D(TString("Ddstar")+l,histTitle,nbinsz,zbins);
    DPi[l] = new TH1D(TString("DPi")+l,histTitle,nbinsz,zbins);
    DiffB[l] = new TH1D(TString("DiffB")+l,histTitle,nbinsz,zbins);
    CascL[l] = new TH1D(TString("CascL")+l,histTitle,nbinsz,zbins);
    LMisID[l] = new TH1D(TString("LMisID")+l,histTitle,nbinsz,zbins);
    Comb[l] = new TH1D(TString("Comb")+l,histTitle,nbinsz,zbins);
    RealD[l] = new TH1D(TString("RealD")+l,histTitle,nbinsz,zbins);
    Cont[l] = new TH1D(TString("Cont")+l,histTitle,nbinsz,zbins);
  }

  cout << "filling output histograms" << endl;

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

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

  THStack hAll[nbinsTot];

  for (int l=0; l<nTotBins; l++) {
    //cout << "l = " << l << endl;

    cx->cd(l+1);

    //cout << "setting marker and color\n";
    OnPeak[l]->SetMarkerStyle(29);  
    OffPeak[l]->SetMarkerStyle(29);  
    Dlnu[l]->SetFillColor(2);
    Dstar[l]->SetFillColor(3);
    Ddstar[l]->SetFillColor(4);
    DPi[l]->SetFillColor(5);
    DiffB[l]->SetFillColor(6);
    CascL[l]->SetFillColor(7) ;
    LMisID[l]->SetFillColor(8) ;
    Comb[l]->SetFillColor(21);
    RealD[l]->SetFillColor(9);
    Cont[l]->SetFillColor(12);
    
    //cout << "stacking\n";
    hAll[l] = new THStack("hAll"+l,histTitle);
    hAll[l].Add(Cont[l]);
    hAll[l].Add(RealD[l]);
    hAll[l].Add(Comb[l]);
    hAll[l].Add(LMisID[l]);
    hAll[l].Add(CascL[l]);
    hAll[l].Add(DiffB[l]);
    hAll[l].Add(DPi[l]);
    hAll[l].Add(Ddstar[l]);
    hAll[l].Add(Dstar[l]);
    hAll[l].Add(Dlnu[l]);

    hAll[l].SetMinimum(min);
    if (requireSetMax)  hAll[l].SetMaximum(max);
    else hAll[l].SetMaximum(OnPeak[l]->GetMaximum()*1.2);
    hAll[l].Draw("hist");
    OnPeak[l]->Draw("e1psames");
    
  }

  cx->Print(outPlotName);




}