#include "TCanvas.h"
#include "TPad.h"
#include "TMultiGraph.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TAxis.h"
#include "TLegend.h"
#include "TMath.h"
#include "TLatex.h"
#include "TROOT.h"
#include "TPaveStats.h"
#include "TStyle.h"
#include <cstdlib> 

char restit[100];
char resfile[100];
char flare[100];
double x1=55440;
double x2=55650;
float offx=1.3;
float offy=1.3;
float margr=0.08;
float w=3;
float margl=0.12;
float line=2;
int paramIndex = 0;
const int nGraphs = 6;
const int nRanges = 2;
const int nParameters = 5;
TGraphErrors *graphs[nGraphs];
int markerColors[] = {kBlue, kGreen, kOrange, kMagenta, kYellow+3, kBlack};
int fitColors[] = {kGreen, kYellow, kRed};
int markerStyles[] = {22, 29, 20, 33, 21, 33};
const char *dataFiles[] = {"se1fe.txt", "se2fe.txt", "se3fe.txt", "se4fe.txt", "se5fec.txt", "se6dtec.txt"};
const char *labels[nGraphs+1] = { "S1", "S2", "S3", "S4", "S5", "S6", "Fit"};
TF1 *flareFits[nRanges]; 
double chi2[nRanges+1];
double ndf[nRanges+1];
double redchi2[nRanges+1];
double totalFitParams[nRanges * nParameters];
double totalFitErrors[nRanges * nParameters];
double initialParams[nRanges][nParameters] = { 
	{0, 1.8e-6, 55490, 40, 30}, 
	{0, 1.e-6, 55580, 10, 20} };
double paramLimits[nRanges][nParameters][2] = { 
	{{0., 0.1e-6}, {0.5e-6, 4.e-6}, {55480, 55500}, {10, 50}, {10, 60}}, 
	{{0., 0.1e-6}, {0.5e-6, 1.1e-6}, {55575, 55690}, {5, 12}, {5, 30}}
    };
const char *paramNames[nParameters] = {"F_c", "F_0", "t_0", "T_r", "T_d"};
double subranges[nRanges][2] = { 
	{55450, 55550}, 
	{55570, 55600} };
double totalFitMin = subranges[0][0];
double totalFitMax = subranges[nRanges-1][1];
int bindays=3;
int binWidth;
int nbins = static_cast<int>((x2 - x1) / binWidth);




void multiflares() {
	gStyle->SetOptFit(0);
	ofstream results;
   sprintf(flare,"S1");
   sprintf(restit,"***************** Fit Results %s***********************",flare);
   sprintf(resfile,"FitResults%s.txt",flare);
   TString outfile = TString::Format("Flare%s.png",flare).Data(); 
   results.open(resfile,  ios::app);

   TCanvas *c01 = new TCanvas("c01","graph",1280,1024);
   //c01->SetGrid();
   gPad->SetLeftMargin(margl);
   gPad->SetRightMargin(margr);
   auto mg = new TMultiGraph();
  
  for (int i = 0; i < nGraphs; i++) {
        graphs[i] = new TGraphErrors(dataFiles[i], "%lg %lg %lg %lg");
        printf("Graph %d has %d points\n", i, graphs[i]->GetN());
        graphs[i]->SetMarkerColor(markerColors[i]);
        graphs[i]->SetLineColor(markerColors[i]);
        graphs[i]->SetMarkerStyle(markerStyles[i]);
        graphs[i]->SetMarkerSize(1);
        graphs[i]->SetLineWidth(2);
        mg->Add(graphs[i]);
    }
    
   mg->GetYaxis()->SetTitleOffset(offy);
   mg->GetXaxis()->SetTitleOffset(offx);
   mg->GetYaxis()->SetTitleSize(40);
   mg->GetYaxis()->SetTitleFont(43);
   mg->GetYaxis()->SetLabelFont(43);
   mg->GetYaxis()->SetLabelSize(40);
   mg->GetXaxis()->SetTitleSize(40);
   mg->GetXaxis()->SetTitleFont(43);
   mg->GetXaxis()->SetLabelFont(43);
   mg->GetXaxis()->SetLabelSize(40);
   mg->SetTitle("");
   mg->GetXaxis()->SetTitle("Time [MJD]");
   mg->GetYaxis()->SetTitle("Flux [10^{-6} ph cm^{-2}s^{-1}]");
   mg->GetXaxis()->SetLimits(x1,x2);
   mg->SetMinimum(0);
   mg->SetMaximum(3.e-6);
   //mg->GetXaxis()->SetNdivisions(nbins); 
   mg->Draw("ap");  
   
 for (int i = 0; i < nRanges; i++) {
        flareFits[i] = new TF1(Form("flareFit%d", i), "[0] + [1] * pow((exp(([2] - x)/[3]) + exp((x - [2])/ [4])),-1)", subranges[i][0], subranges[i][1]);
        flareFits[i]->SetParameters(initialParams[i]);
        for (int j = 0; j < nParameters; j++) {
            flareFits[i]->SetParName(j, Form("%s%d", paramNames[j], i));
             flareFits[i]->SetParLimits(j, paramLimits[i][j][0], paramLimits[i][j][1]);
        }
    }
   
   for (int j = 0; j < nRanges; j++) {
            flareFits[j]->SetRange(subranges[j][0], subranges[j][1]);
            mg->Fit(flareFits[j], "R+");
            chi2[j] = flareFits[j]->GetChisquare();
            ndf[j] = flareFits[j]->GetNDF();
            redchi2[j] = chi2[j]/ndf[j];
            flareFits[j]->SetLineColor(fitColors[j]);
            flareFits[j]->Draw("same");
        }
        	TF1 *totalFit = new TF1("totalFit", [](double *x, double *par) {
        double sum = 0;
        for (int i = 0; i < nRanges; ++i) {
            sum += par[i * nParameters] + par[i * nParameters + 1] * pow((exp((par[i * nParameters + 2] - x[0]) / par[i * nParameters + 3]) + exp((x[0] - par[i * nParameters + 2]) / par[i * nParameters + 4])),-1);
        }
        return sum;
    }, totalFitMin, totalFitMax, nRanges * nParameters);
   
	for (int i = 0; i < nRanges; i++) {
        for (int j = 0; j < flareFits[i]->GetNpar(); j++) {
            totalFit->SetParameter(i * nParameters + j, flareFits[i]->GetParameter(j));
        }
    }
    for (int i = 0; i < nRanges; i++) {
 		   for (int j = 0; j < nParameters; j++) {
        	paramIndex = i * nParameters + j;
        	totalFit->SetParName(paramIndex, Form("%s%d", paramNames[j], i+1));
    	}	
	}

    mg->Fit(totalFit, "R+");
     chi2[nRanges] = totalFit->GetChisquare();
            ndf[nRanges] = totalFit->GetNDF();
            redchi2[nRanges] = chi2[nRanges]/ndf[nRanges];
    totalFit->SetLineColor(fitColors[2]);
    totalFit->Draw("same");
    c01->Modified();
    c01->Update();
     
  
  for (int i = 0; i < nRanges; i++) {
    for (int j = 0; j < flareFits[i]->GetNpar(); ++j) {
        totalFitParams[i * nParameters + j] = totalFit->GetParameter(i * 5 + j);
        totalFitErrors[i * nParameters + j] = totalFit->GetParError(i * 5 + j);
    }
 }
   totalFit->GetYaxis()->SetTitleOffset(offy);
   totalFit->GetXaxis()->SetTitleOffset(offx);
   
   TLegend* leg = new TLegend(0.75, 0.75, .85, .85);
   leg->SetHeader("Legend");
   leg->SetNColumns(1);
   for (int i = 0; i < nGraphs; ++i) {
   		 leg->AddEntry(graphs[i], labels[i], "ap");
	}
   leg->Draw();
   gPad->Update();
    c01->Modified();
   c01->Update();
   c01->Print(outfile);
	results << "***********GENERAL INFORMATION**************"  << endl;
	results << "Plot=  " << outfile << endl;
	for (int i = 0; i < nRanges; ++i) {
		results << "Fit range " << i+1 << "[" << subranges[i][0] << ", " << subranges[i][1] << "]" << endl;
		
	}
    results << "Total fit range: [" << totalFitMin << ", " << totalFitMax << "]" << endl;
    results << endl;
    results << "***********INIZIALIZATION PARAMETERS**************" <<  endl;
    for (int i = 0; i < nRanges; ++i) {
    	results << endl;
    	results << "Inizialization parameters for subrange " << i+1 << ":" << endl;
    	for (int j = 0; j < nParameters; j++) {
        	results << paramNames[j] << " = " << initialParams[i][j] << " limited in : [" << paramLimits[i][j][0] << ", " << paramLimits[i][j][1] << "]" << endl;
	    }
	}
    results << endl;
    results << "***********PARTIAL FIT RESULTS **************" <<  endl;
    for (int i = 0; i < nRanges; ++i) {
	results << endl;
        results << "Fit parameters for subrange " << i+1 << ":" << endl;
        for (int j = 0; j < nParameters; ++j) {
            int paramIndex = i * nParameters + j;
            results << paramNames[j] << i << " = " << flareFits[i]->GetParameter(j) << " +/- " << flareFits[i]->GetParError(j) << endl;
        }
        results << "Chisquare = " << chi2[i] << " "  << "Ndf " << ndf[i] << " "<<  "Reduced Chisquare " << redchi2[i] << endl;
    }
    results << endl;
    results << "***********TOTAL FIT RESULTS**************" <<  endl;
   for (int i = 0; i < nRanges; ++i) {
   	results << endl;
    results << "Fit parameters for subrange " << i+1 << " results:" << std::endl;
    for (int j = 0; j < nParameters; ++j) {
       paramIndex = i * nParameters + j;
        results << paramNames[j] << i << " = " << totalFit->GetParameter(paramIndex) << " +/- " << totalFit->GetParError(paramIndex) << std::endl;
    }
}
	results << "Chisquare = " << chi2[nRanges] << " "  << "Ndf " << ndf[nRanges] << " "<<  "Reduced Chisquare " << redchi2[nRanges] << endl;
	results <<  endl;  
    delete c01;

}