#include "Riostream.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TH1F.h"
#include "TH2F"
#include "TF1.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TMinuit.h"
#include "TString.h"
#include <fstream>
#include "TAxis.h"
#include "TArrayD.h"
#include "TFractionFitter.h"
#include "TVirtualFitter.h"
#include "Fit/FitConfig.h"
#include "Fit/Fitter.h"
#include "TFitResult.h"

using namespace std;

TH1F *hist_data;
TH1F *hist_605;
TH1F *hist_613;
TH1F *hist_692;
TH1F *hist_712;
TH1F *hist_ground;

void TFraction()
{
	TFile *f = new TFile("histograms.root","READ");

	f->ls();

	hist_data   = (TH1F *)f->Get("hist1");
	hist_605    = (TH1F *)f->Get("hist_605");
	hist_613    = (TH1F *)f->Get("hist_613");
	hist_692    = (TH1F *)f->Get("hist_692");
	hist_712    = (TH1F *)f->Get("hist_712");
	hist_ground = (TH1F *)f->Get("hist_ground");
	
//////////////////////////////////////////////////////////////////////////////////////
///////////////////////Fits MC fractions to data histogram.//////////////////////////
	const int nh = 5;
	const int CM = 100;
	

	double par1[CM], par2[CM], par3[CM], par4[CM], par5[CM];
	double br[CM][nh],brer[CM][nh];
	double chisq[CM];
	double br_add[CM],br_sub[CM];
	double before_add[CM];	
	
	for(int j=0; j<CM;j++){	
	
		//par1[j] = 0.05 + j/800.0;
		//par2[j] = 0.06 + j/800.0;
		//par3[j] = 0.031+ j/800.0;
		par4[j] = 0.020 + j/800.0;
		//par5[j] = 0.50 + j/600.0;
	
		TVirtualFitter::SetDefaultFitter("Minuit");

		TObjArray *mc = new TObjArray(nh);        // MC histograms are put in this array
		mc->Add(hist_605);
		mc->Add(hist_613);
		mc->Add(hist_692);
		mc->Add(hist_712);
		mc->Add(hist_ground);

		TFractionFitter* fit = new TFractionFitter(hist_data, mc); // initialise
		fit -> SetRangeX(hist_data->GetXaxis()->FindBin(1.8),hist_data->GetXaxis()->FindBin(16.0));			
		cout <<" " << endl;
		TVirtualFitter* vFit = fit->GetFitter();
		
//////////////////////////////////////////////////////////////////////////////////////////////////////
		
		//double vstart[nh] = {par1[j],  0.134139, 0.092894, 0.080129, 0.578215};
		//double vstart[nh] = {0.114623, par2[j],  0.092894, 0.080129, 0.578215};
		//double vstart[nh] = {0.114623, 0.134139, par3[j],   0.080129, 0.578215};	
		double vstart[nh] = {0.114623, 0.134139, 0.092894, par4[j],   0.578215 + 0.080129- par4[j]};
		//double vstart[nh] = {0.114623, 0.134139, 0.092894, 0.080129, par5[j]};
		
		before_add[j] = 0.114623 + 0.134139 + 0.092894 + par4[j] + 0.578215 + 0.080129- par4[j];
		//double vstart[nh] = {0.114623, 0.134139, 0.092894, 0.080129, 0.578215};
    		double step[nh]    = {0.000001, 0.000001, 0.000001, 0.000001, 0.000001};
    	
		fit->GetFitter()->SetParameter(0, "6.05 branch", vstart[0], step[0], 0, 1);
		fit->GetFitter()->SetParameter(1, "6.13 branch", vstart[1], step[1], 0, 1);
		fit->GetFitter()->SetParameter(2, "6.92 branch", vstart[2], step[2], 0, 1);
		fit->GetFitter()->SetParameter(3, "7.12 branch", vstart[3], step[3], 0, 1);
		fit->GetFitter()->SetParameter(4, "g.st branch", vstart[4], step[4], 0, 1);
	
	
	
		//fit->GetFitter()->FixParameter(0);
		//fit->GetFitter()->FixParameter(1);
		//fit->GetFitter()->FixParameter(2);
		fit->GetFitter()->FixParameter(3);
		//fit->GetFitter()->FixParameter(4);
	
//////////////////////////////////////////////////////////////////////////////////////////////////

		int status = fit->Fit();  // perform the fit
		cout << endl;
		cout << "Fit status: " << status << endl;
		
		fit->ErrorAnalysis(0.5);	// minos error calculation
	
		chisq[j] = fit->GetChisquare();
	
		for (int iz=0; iz<nh; iz++){
		
            	fit -> GetResult(iz,br[j][iz],brer[j][iz]);
        	}
	}//for loop
	
///////////////////////////////////////////////////////////////////////////////////////////////////
	
	ofstream outfile;
	ofstream outfile1;
	ofstream outfile2;
	ofstream outfile_add;
		
	outfile.open("br_712_vs_chi.dat",ios::binary);
	//outfile1.open("br_712_add.dat",ios::binary);
	//outfile2.open("br_712_sub.dat",ios::binary);
	outfile_add.open("all_712_add.dat",ios::binary);
	
	//outfile.open("br_692_vs_chi.dat",ios::binary);
	//outfile1.open("br_692_add.dat",ios::binary);
	//outfile2.open("br_692_sub.dat",ios::binary);	
	//outfile_add.open("all_692_add.dat",ios::binary);
	
	double all_br_add[CM];
	for(int ii=0; ii<CM; ii++){
	
		outfile << fixed << setprecision(6) << br[ii][0] <<"\t"<< br[ii][1] <<"\t"<< br[ii][2] <<"\t"<< br[ii][3] <<"\t"<< br[ii][4] <<"\t"<< chisq[ii] << endl;
		cout    << fixed << setprecision(6) << br[ii][0] <<"\t"<< br[ii][1] <<"\t"<< br[ii][2] <<"\t"<< br[ii][3] <<"\t"<< br[ii][4] <<"\t"<< chisq[ii] << endl;	
		
		//br_add[ii] = br[ii][2] + br[ii][3];
		//br_sub[ii] = br[ii][2] - br[ii][3];
		all_br_add[ii] = br[ii][0] + br[ii][1] + br[ii][2] + br[ii][3] + br[ii][4];
		
		outfile1 << fixed << setprecision(6) << br_add[ii] << "\t" << chisq[ii] << endl;
		outfile2 << fixed << setprecision(6) << br_sub[ii] << "\t" << chisq[ii] << endl;
		outfile_add << fixed << setprecision(6) <<before_add[ii] <<"\t"<< all_br_add[ii] << "\t" << chisq[ii] <<endl;
		
			
	} 
	outfile.close();
	outfile1.close();
	outfile2.close();
	outfile_add.close();			
	
}