/* vim:set noexpandtab tabstop=4 wrap */
#include "TROOT.h"
#include "TSystem.h"
#include "TSystemFile.h"
#include "TSystemDirectory.h"
#include "TApplication.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "THStack.h"
#include "TNtuple.h"
#include "TParticlePDG.h"
#include "TDatabasePDG.h"
#include "TVector3.h"
#include "TMath.h"
#include "TCanvas.h"
#include <string>
#include <iostream>

#ifndef VERBOSE
#define VERBOSE 1
#endif

void stacktest(){
	const char* outdir="/tmp";
	
	// histograms to compare distributions
	TH1D truenuKEhist("truenuKEhist","True Neutrino KE",50,0,2000);
	truenuKEhist.FillRandom("gaus",1000);
	TH1D truemuKEhist("truemuKEhist","True Muon KE",50,0,2000);
	truemuKEhist.FillRandom("landau",1000);
	
	TH1D reconuKEhist("reconuKEhist","Reco Neutrino KE",50,0,2000);
	reconuKEhist.FillRandom("gaus",1000);
	TH1D recomuKEhist("recomuKEhist","Reco Muon KE",50,0,2000);
	recomuKEhist.FillRandom("landau",1000);
	
	// Save PNGs of the true and reco histograms overlaid
	TCanvas c1;
	THStack* thestack = new THStack("thestack","Stacked Histos");
	thestack->Add(&truenuKEhist);
	thestack->Add(&reconuKEhist);
	thestack->Draw("nostack");
	thestack->GetXaxis()->SetTitle("Neutrino Energy [MeV]");
	thestack->GetYaxis()->SetTitle("Num Events");
	c1.SaveAs(TString::Format("%s/reconstructed_neutrino_energy.png",outdir));
	delete thestack;
	
#ifdef VERBOSE
	cout<<"making 2"<<endl;
#endif
	THStack* thestack2 = new THStack("thestack","Stacked Histos");
	thestack2->Add(&truemuKEhist);
	thestack2->Add(&recomuKEhist);
// XXX uncomment either of the below lines, and the program sefaults....
//	thestack2->GetXaxis()->SetTitle("Muon Energy [MeV]");
//	thestack2->GetYaxis()->SetTitle("Num Events");
	thestack2->Draw("nostack");
	c1.SaveAs(TString::Format("%s/reconstructed_muon_energy.png",outdir));
	delete thestack2;
}