#include <TFile.h>
#include <TH1F.h>
#include <TPad.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string>
using namespace std;

////////////////////
// Set parameters //
////////////////////
double Qgate1 = 8.20; // Qgate min, MeV
double Qgate2 = 9.20; // Qgate max, MeV
double threshold = 120; // Ethreshold, keV
int resonance[5] = {1, 2, 3, 4, 5};  // resonances (by number) to plot, watch for 0 or 1 numbering
////////////////////

double resbinsx, resbinsy, resbinsz;
vector<double> respos;
vector<double> resbinmin;
vector<double> resbinmax;
int numresonances, Qgatebin1, Qgatebin2;
char rootfilename[70], histoname[60];
TH1D *M[10], *Mraw[10];
TH1F *Eg[10][10], *Esum[10][10];


int viewspectra( void ) {

	ifstream binedgein("justbinedges.txt", ifstream::in);				// read in resonance binning data
	if(!binedgein.is_open()) {
		cout << " Bin edges file (binedges.txt) not found." << endl;
		exit(0);
		}
	cout << "--> Reading bin edges : " << endl;
	while(!binedgein.eof()){
		binedgein  >> resbinsx >> resbinsy >> resbinsz;
		if(binedgein.eof()) break;
		respos.push_back (resbinsx);
		resbinmin.push_back (resbinsy);
		resbinmax.push_back (resbinsz);
		}
	binedgein.close();
	numresonances = respos.size();
	printf("Check Number of Resonances: %d \n", numresonances);

 TCanvas *canvas = new TCanvas("canvas","multipads",900,700);		// create canvas for viewing
 gStyle->SetOptStat(0);
 canvas->Divide(3,3,0,0);


 sprintf(rootfilename,"Mo97_2010_10ns_%3.0f_%3.0f_%3.0fkeV_resonances_BGm.root",Qgate1*100,Qgate2*100, threshold);
 cout << "Root File name: " << rootfilename << endl;			// read from resonance root file
 TFile *f1 = TFile::Open( rootfilename , "READ");													
	for (int j=0;j<5;j++){		// loop over resonances
		for (int k=2;k<=4;k++){		// loop over multiplicities
			int i = resonance[j];	
			sprintf(histoname,"resonance%0#7.2feV/Egamma_%0#7.2feV_Mcl%i",respos[i],respos[i],k);
//			cout << histoname << endl;
			Eg[j][k] = (TH1F*) f1->Get(histoname);
			Eg[j][k]->SetDirectory(0);
			Qgatebin1 = Eg[j][k]->GetXaxis->FindBin(Qgate1);
			Qgatebin2 = Eg[j][k]->GetXaxis->FindBin(Qgate2);
			Eg[j]k]->Scale(1/Eg[j][k]->Integral(Qgatebin1,Qgatebin2));
			canvas->cd(k-1);											
			if (j==0) {Eg[j][k]->Draw();}
			else {Eg[j][k]->Draw("same");}

			sprintf(histoname,"resonance%0#7.2feV/Esum_%0#7.2feV_Mcl_%i",respos[i],respos[i],k);
			cout << histoname << endl;
			Esum[j][k] = (TH1F*) f1->Get(histoname);
			Esum[j][k]->SetDirectory(0);
			Esum[j][k]->Scale(1/Esum[j][k]->Integral(Qgatebin1,Qgatebin2));
			canvas->cd(k+2);											
			if (j==0) {Esum[j][k]->Draw();}
			else {Esum[j][k]->Draw("same");}
			}

		sprintf(histoname,"resonance%0#7.2feV/m_final_%0#7.2feV",respos[i],respos[i]);
//		cout << histoname << endl;
		M[j] = (TH1D*) f1->Get(histoname);
		M[j]->SetDirectory(0);
//		cout << M[j]->Integral() << endl;
		canvas->cd(7);											
		if (j==0) {M[j]->Draw();}
		else {M[j]->Draw("same");}

		sprintf(histoname,"resonance%0#7.2feV/m_raw_%0#7.2feV",respos[i],respos[i]);
//		cout << histoname << endl;
		Mraw[j] = (TH1D*) f1->Get(histoname);
		Mraw[j]->SetDirectory(0);
//		cout << Mraw[j]->Integral() << endl;
		canvas->cd(7);											
		if (j==0) {Mraw[j]->Draw();}
		else {Mraw[j]->Draw("same");}
		}
	f1->Close();
return 0;
}