#include <iostream>
#include <fstream>

gROOT->Reset();

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

void histMk_3DMC_DchKPiPi() {

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

  TString outHistName("hist_AllMC_"+oNameBase+"_DchKPiPi.root");
  TString histTitle("All MC (DchKPiPi)");

  double r_lumi = 1; //lumi_Run3_OnPeak/lumi_Run2_OnPeak;


  //################################################
  cout << "start";
  // Open input files and get 3D histos
  TFile *f[nFiles_DchKPiPi];
  TH3D * inHist[nFiles_DchKPiPi];
  for (int l=2; l<nFiles_DchKPiPi; l++) {
    f[l] = new TFile(inFile_DchKPiPi[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_DchKPiPi][nbinsx][nbinsy][nbinsz];
  double err_3D[nFiles_DchKPiPi][nbinsx][nbinsy][nbinsz];
  for (int l=2; 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[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 all MC
  double num_3DMC[nbinsx][nbinsy][nbinsz];
  double var_3DMC[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_3DMC[i][j][k] = 0;
	var_3DMC[i][j][k] = 0;
      }
    }
  }
  for (int i=0; i<nbinsx; i++) {
    for (int j=0; j<nbinsy; j++) {
      for (int k=0; k<nbinsz; k++) {
	for (int l=2; l<=27; l++) {
	  num_3DMC[i][j][k] += r_lumi*num_3D[l][i][j][k];
	  var_3DMC[i][j][k] 
	    += r_lumi*r_lumi*err_3D[l][i][j][k]*err_3D[l][i][j][k];
	}
      }
    }
  }

  // Define new hist and fill it
  TH3D * sub 
    = new TH3D("sub", histTitle, nbinsx, xbins, nbinsy, ybins, nbinsz, zbins);
  for (int i=0; i<nbinsx; i++) {
    for (int j=0; j<nbinsy; j++) {
      for (int k=0; k<nbinsz; k++) {
	double err_3DMC = sqrt(var_3DMC[i][j][k]);
	sub->SetBinContent(i+1,j+1,k+1,num_3DMC[i][j][k]);
	sub->SetBinError(i+1,j+1,k+1,err_3DMC);
      }
    }
  }

  TFile *f_out = new TFile(outHistName, "Recreate");

  sub->Write();
 
  f_out->Write();
  f_out->Close();

  cout << "histogram " << outHistName << " is saved." << endl;

}