#include <iostream>
#include <fstream>

gROOT->Reset();

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

void plotMk_3DMC_ProjY_DchKPiPi() {

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

  TString outPlotName("plot_"+oNameBase+"_"+vary+"_DchKPiPi.eps");
  TString histTitle(vary+"_DchKPiPi");

  // Set max value of plot
  int min = ymin;
  int max = ymax;
  //setMax = true;
  //int min = -100;
  //int max = 40000;

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

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

  // 1D arrays to be projected
  cout << "define summed 1D arrays" << endl;
  double num_1D_V1[nFiles_DchKPiPi][nbinsy];
  double var_1D_V1[nFiles_DchKPiPi][nbinsy];
  double err_1D_V1[nFiles_DchKPiPi][nbinsy];
  for (int l=0; l<nFiles_DchKPiPi; l++) {
    for (int j=0; j<nbinsy; j++) {
      num_1D_V1[l][j] = 0;
      var_1D_V1[l][j] = 0;
      err_1D_V1[l][j] = 0;
    }
  }
  
  for (int l=0; l<nFiles_DchKPiPi; l++) {
    for (int j=0; j<nbinsy; j++) {
      for (int i=0; i<nbinsx; i++) {
	for (int k=0; k<nbinsz; k++) {
	  num_1D_V1[l][j] += c_DchKPiPi[l] * num_3D_V1[l][i][j][k];
	  var_1D_V1[l][j] 
	    += c_DchKPiPi[l]*c_DchKPiPi[l]*err_3D_V1[l][i][j][k]*err_3D_V1[l][i][j][k];
	}
      }
      err_1D_V1[l][j] = sqrt(var_1D_V1[l][j]);
    }
  }


  //Define 1D histograms and fill them
  TH1D *outHist_V1[nFiles_DchKPiPi];
  for (int l=0; l<nFiles_DchKPiPi; l++) {
    outHist_V1[l] = new TH1D("outHist_V1"+l,histTitle,nbinsy,ybins);
  }

  for (int l=0; l<nFiles_DchKPiPi; l++) {
    for (int j=0; j<nbinsy; j++) {
      outHist_V1[l]->SetBinContent(j+1,num_1D_V1[l][j]);
      outHist_V1[l]->SetBinError(j+1,err_1D_V1[l][j]);
    }
  }

  // Open input files and get 3D histos
  TFile *f_V2[nFiles_DchKPiPi];
  TH3D * inHist_V2[nFiles_DchKPiPi];
  for (int l=0; l<nFiles_DchKPiPi; l++) {
    f_V2[l] = new TFile(inFile_V2_DchKPiPi[l]);
    f_V2[l]->ls();
    inHist_V2[l] = (TH3D *)f_V2[l]->Get("sub;1");
  }

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

  // 1D arrays to be projected
  cout << "define summed 1D arrays" << endl;
  double num_1D_V2[nFiles_DchKPiPi][nbinsy];
  double var_1D_V2[nFiles_DchKPiPi][nbinsy];
  double err_1D_V2[nFiles_DchKPiPi][nbinsy];
  for (int l=0; l<nFiles_DchKPiPi; l++) {
    for (int j=0; j<nbinsy; j++) {
      num_1D_V2[l][j] = 0;
      var_1D_V2[l][j] = 0;
      err_1D_V2[l][j] = 0;
    }
  }
  
  for (int l=0; l<nFiles_DchKPiPi; l++) {
    for (int j=0; j<nbinsy; j++) {
      for (int i=0; i<nbinsx; i++) {
	for (int k=0; k<nbinsz; k++) {
	  num_1D_V2[l][j] += c_DchKPiPi[l] * num_3D_V2[l][i][j][k];
	  var_1D_V2[l][j] 
	    += c_DchKPiPi[l]*c_DchKPiPi[l]*err_3D_V2[l][i][j][k]*err_3D_V2[l][i][j][k];
	}
      }
      err_1D_V2[l][j] = sqrt(var_1D_V2[l][j]);
    }
  }


  //Define 1D histograms and fill them
  TH1D *outHist_V2[nFiles_DchKPiPi];
  for (int l=0; l<nFiles_DchKPiPi; l++) {
    outHist_V2[l] = new TH1D("outHist_V2"+l,histTitle,nbinsy,ybins);
  }

  for (int l=0; l<nFiles_DchKPiPi; l++) {
    for (int j=0; j<nbinsy; j++) {
      outHist_V2[l]->SetBinContent(j+1,num_1D_V2[l][j]);
      outHist_V2[l]->SetBinError(j+1,err_1D_V2[l][j]);
    }
  }

  //#############################################################
  cout << "plotting\n";
  // Create a canvas and plot
  TCanvas * cx = new TCanvas("cx",histTitle,0,0,700,500);
  cx->Divide(7,5);
  gStyle->SetOptStat(000000);

  for (int l=0; l<nFiles_DchKPiPi; l++) {
    cx->cd(l+1);

    outHist_V1[l]->SetLineColor(2);  
    outHist_V2[l]->SetLineColor(3);
    //outHist_V1[l]->SetMaximum(7000);
    outHist_V1[l]->Draw("hist");
    outHist_V2[l]->Draw("histsames");
  }
    
  cx->Print(outPlotName);

}