#include <TFile.h>
#include <TProfile.h>
#include <TCanvas.h>
#include <TFrame.h>
#include <TROOT.h>
#include <TSystem.h>
#include <TCint.h>
#include <TRandom3.h>
#include <TRandom.h>
#include <iostream>
#include <TH1.h>
#include <TH2.h>
#include <TMath.h>
#include "TCanvas.h"
#include <vector>
#include <TTree.h>
#include <TKey.h>
#include <TStopwatch.h>
#include <fstream>
#include "stdio.h"
#include "TCanvas.h"
#include <TSpectrum.h>
#include "TF1.h"
#include "TVirtualFitter.h"
int main(int argc, char *argv[]){

	if(argc != 4 ){
		std::cout << "Syntax: <" << argv[0] << ">  <inputlist.txt> <output.root> <Output Title> <minX>, <maxX>"<< std::endl;
		return EXIT_FAILURE;
	}else{

		ifstream file;
		std::string fileName = argv[1];
		TString outputname = argv[2];
		TString outputTitle = argv[3];
		file.open(fileName.c_str());
		if (!file){
			std::cout << "Error in opening file" << std::endl;
			return EXIT_FAILURE;
		}
		std::string lineread;
		TFile *output;
		TFile * inputFile; 
		Int_t comp = 0;

		Int_t nspec = 11;
		TH1D* Spectra[nspec];
		TH1D* Spectra2[nspec];
		TH1D* Spectra3[nspec];
		TH1D* Spectra4[nspec];
		TH1D* Spectra5[nspec];
		TH1D* Spectra6[nspec];
		TH1D* Spectra7[nspec];
		TCanvas *AllSpectra = new TCanvas("AllSpectra", "All spectra overlaid", 1);
		TCanvas *CorrectedSpectra = new TCanvas("Corrected", "Corrected Spectra overlaid", 1);
		TCanvas *BackGroundSpectra = new TCanvas("Background", "Backgrounds Overlaid", 1);
		TCanvas *PeakSearches = new TCanvas("PeakSearches", "Peaks overlaid", 1);
		TCanvas *Calibrated = new TCanvas("Calibrated according to energy", "Peaks overlaid", 1);
		output = new TFile(outputname, "RECREATE", "All spectra");
		output->SetCompressionLevel(comp);

		
		Double_t meanMax = 0;
		Int_t j = 0;
		while(std::getline(file, lineread)){
			TString inputName = lineread;
			TString histNewName = inputName.Remove(inputName.Sizeof() - 6, 6);
			std::cout << "inserting histogram " << histNewName << " into " << outputname << std::endl;
			inputFile = new TFile(lineread.c_str());
			inputFile->cd();
			Spectra[j] = (TH1D*)gFile->Get("RawSpectrum");
			Spectra2[j] = (TH1D*)Spectra[j]->Clone();
			Spectra2[j]->SetName(histNewName);
			Spectra2[j]->Rebin(8);
			Spectra2[j]->ShowBackground(50, "BackOrder8");
			Spectra3[j] = (TH1D*)gDirectory->Get(histNewName+"_background");
			Spectra3[j]->SetName(histNewName+"_backgnd");
			output->cd();
			Spectra2[j]->Write();
			Spectra3[j]->Write();
			Spectra4[j] = (TH1D*)Spectra2[j]->Clone();
			Spectra4[j]->SetName(histNewName+"_corrected");
			Spectra4[j]->Add(Spectra3[j], -1.0);
			Spectra4[j]->GetXaxis()->SetRangeUser(200,2000);
			Spectra4[j]->GetYaxis()->SetRangeUser(-0.1,3);
			Spectra4[j]->Write();
			AllSpectra->cd();
			Spectra2[j]->Draw("SAME");
			BackGroundSpectra->cd();
			Spectra3[j]->Draw("SAME");
			CorrectedSpectra->cd();
			Spectra4[j]->Draw("SAME");

			PeakSearches->cd();
			Float_t fPositionX[100];
			Float_t fPositionY[100];   
			Spectra5[j] = (TH1D*)Spectra2[j]->Clone();
			Spectra5[j]->SetName(histNewName+"_peaks");
			Spectra5[j]->Draw();
			
			TSpectrum *s = new TSpectrum(10);	//FIXME: Can't initialise this with maxpeaks
			Int_t nfound = s->Search(Spectra5[j],3,"new",0.005);
			printf("Found %d candidate peaks to fit\n",nfound);
			Float_t *xpeaks = s->GetPositionX();
			TF1 * fit[nfound];

			for (Int_t i = 0; i < nfound; i++) {
				Double_t xp=xpeaks[i];
				//Int_t bin = 1 + Int_t(xp + 0.5);
				Int_t bin = Spectra5[j]->GetXaxis()->FindBin(xp);
				fPositionX[i] = Spectra5[j]->GetBinCenter(bin);
				fPositionY[i] = Spectra5[j]->GetBinContent(bin);
				fit[i]= new TF1(histNewName+"_fitfn", "gaus", xp-11,xp+11);
				Spectra6[j] = (TH1D*)Spectra4[j]->Clone();
				Spectra6[j]->SetName(histNewName+"_fit");
				Spectra6[j]->GetXaxis()->SetRangeUser(xp-11,xp+11);
				Spectra6[j]->Fit(fit[i],"R+");
				Spectra6[j]->GetXaxis()->SetTitle("ChannelNumber");
				Double_t sigma=(Double_t)(fit[i]->GetParameter(2));
				Double_t mean=(fit[i]->GetParameter(1));

				//Spectrum6->GetXaxis()->SetRangeUser(mean-(3*sigma),mean+(3*sigma));
				Spectra6[j]->GetXaxis()->SetRangeUser(0,4096);
				Double_t prob=Spectra6[j]->Integral();
				std::cout << histNewName << " peak number " << i << " has mean: " << mean << " sigma: " << sigma << " integral " << prob << std::endl;
				Spectra6[j]->Draw("SAME");
				if(mean > meanMax){meanMax = mean;}
			
			}

				j++;
		}

		std::cout << "meanMax: " << meanMax << std::endl;
		Calibrated->cd();
		for(int k = 0; k<nspec; k++){
		Spectra7[k]=(TH1D*)Spectra4[k]->Clone();
		Spectra7[k]->SetName((TString)Spectra4[k]->GetName()+"_calibrated");
		std::cout << "Bins: " << Spectra4[k]->GetNbinsX() << std::endl; 
		Spectra7[k]->SetBins(Spectra4[k]->GetNbinsX(), 0, (Spectra4[k]->GetNbinsX()/(meanMax /8.0))*5.39);
		Spectra7[k]->GetXaxis()->SetTitle("Energy (MeV)");
		Spectra7[k]->Draw("SAME");
		
		}
		AllSpectra->Write();
		BackGroundSpectra->Write();
		CorrectedSpectra->Write();
		PeakSearches->Write();
		Calibrated->Write();
		output->Write();
	}
}