#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <map>

#include "TFile.h"

#include "IUAmpTools/AmpToolsInterface.h"
#include "IUAmpTools/FitResults.h"
#include "IUAmpTools/Histogram1D.h"
#include "IUAmpTools/Histogram2D.h"

#include "AMPTOOLS_DATAIO/ROOTDataReader.h"
#include "AMPTOOLS_AMPS/Uniform.h"

#include "HCA_Amplitudes/HCAmplitude_Chain.h"
#include "HCA_Amplitudes/BreitWigner.h"
#include "HCA_Amplitudes/TChannel_Basic.h"

#include "PlotGenerator_d2pi/PlotGenerator_d2pi.h"

using namespace std;
using namespace CLHEP;

int main(int argc, char* argv[])
{
	if(argc < 4)
	{
		cout << "Usage:" << endl << endl;
		cout << "\tplotResults <fit results name> <output file name> <Bin Mean Sqrt(S)>" << endl << endl;
		return 0;
	}

	// parse the command line parameters
	string locResultsName(argv[1]);
	string outname(argv[2]);
	istringstream locSqrtSStream(argv[3]);
	double locBinMeanSqrtS = 0.0;
	locSqrtSStream >> locBinMeanSqrtS;
	cout << "Fit results file name		= " << locResultsName << endl;
	cout << "Output file name		= " << outname << endl;
	cout << "Bin sqrt S = " << locBinMeanSqrtS << endl << endl;

	// load the results and display the configuration info
	FitResults locFitResults(locResultsName);
	if(!locFitResults.valid())
	{
		cout << "Invalid fit results in file:  " << locResultsName << endl;
		return 1;
	}
	locFitResults.configInfo()->display();
	string locReactionName = locFitResults.reactionList()[0];

	// set up an output Root file
	TFile* locOutputROOTFile = new TFile(outname.c_str(), "RECREATE");
	TH1::AddDirectory(kFALSE);

	// set up a PlotGenerator and make plots
	AmpToolsInterface::registerDataReader(ROOTDataReader());
	AmpToolsInterface::registerAmplitude(BreitWigner());
	AmpToolsInterface::registerAmplitude(HCAmplitude_Chain());
	AmpToolsInterface::registerAmplitude(TChannel_Basic());
	AmpToolsInterface::registerAmplitude(Uniform());
	
	PlotGenerator_d2pi locPlotGenerator(locFitResults, locBinMeanSqrtS);
	locPlotGenerator.enableReaction(locReactionName);
	vector<string> amps = locPlotGenerator.uniqueAmplitudes();

	// turn off all amplitudes
	for(unsigned int i = 0; i < amps.size(); ++i)
		locPlotGenerator.disableAmp(i);

	// loop over amplitude configurations
	for(unsigned int iamp = 0; iamp <= amps.size(); ++iamp)
	{
		if(iamp == amps.size())
		{
			// turn on all amplitudes
			for(unsigned int i = 0; i < amps.size(); ++i)
				locPlotGenerator.enableAmp(i);
		}
		else
		{
			// turn on current amplitude
			locPlotGenerator.enableAmp(iamp);
			if(iamp != 0) //turn off previous one
				locPlotGenerator.disableAmp(iamp - 1);
		}

		// loop over different histograms
		for(unsigned int ivar = 0; ivar	< locPlotGenerator.dNumHists; ivar++)
		{
			// loop over data, accMC, and genMC
			for(unsigned int iplot = 0; iplot < PlotGenerator::kNumTypes; iplot++)
			{
				if((iamp != amps.size()) && (iplot == PlotGenerator::kData))
					continue;

				ostringstream locHistNameStream;
				locHistNameStream << "h_" << ivar << "_";

				if(iplot == PlotGenerator::kData)
					locHistNameStream << "data";
				else if(iplot == PlotGenerator::kAccMC)
					locHistNameStream << "accmc";
				else if(iplot == PlotGenerator::kGenMC)
					locHistNameStream << "genmc";

				if (iamp < amps.size())
					locHistNameStream << "_Amp" << iamp;

				Histogram* locHist = locPlotGenerator.projection(ivar, locReactionName, iplot);
				Histogram1D* loc1DHist = dynamic_cast<Histogram1D*>(locHist);
				Histogram2D* loc2DHist = dynamic_cast<Histogram2D*>(locHist);
				if(loc1DHist != NULL)
				{
					TH1F* thist = static_cast<TH1F*>(loc1DHist->toRoot());
					thist->SetName(locHistNameStream.str().c_str());
					locOutputROOTFile->cd();
					thist->Write();
				}
				else if(loc2DHist != NULL)
				{
					TH2F* thist = static_cast<TH2F*>(loc2DHist->toRoot());
					thist->SetName(locHistNameStream.str().c_str());
					locOutputROOTFile->cd();
					thist->Write();
				}
			}
		}
	}

	locOutputROOTFile->Close();

	// print results to the screen
	cout << "TOTAL EVENTS = " << locFitResults.intensity().first << " +- " << locFitResults.intensity().second << endl;
	vector<string> fullamps = locPlotGenerator.fullAmplitudes();
	for(unsigned int i = 0; i < fullamps.size(); ++i)
	{
		vector<string> useamp;
		useamp.push_back(fullamps[i]);
		double locFitFraction = locFitResults.intensity(useamp).first / locFitResults.intensity().first;
		double locFitFractionUncertainty = locFitResults.intensity(useamp).second / locFitResults.intensity().first;
		cout << "FIT FRACTION " << i+1 << " = " << locFitFraction <<	" +- " << locFitFractionUncertainty << endl;
	}

	return 0;
}