#include "TCanvas.h"
#include "TPad.h"
#include "TMultiGraph.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TGraphErrors.h"
#include "TAxis.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TROOT.h"
#include "TPaveStats.h"
#include "TStyle.h"
#include "TFitResult.h"
#include "TF1.h"
#include "TH1D.h"
#include "TMath.h"
#include <assert.h>
#include <cmath>
#include <cstdlib> 
#include <cmath>
#include <iostream>
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <fstream>

char Detector[50];
char header[100];
char type[100];

int x1=274;
int x2=291;

double netarea0=2.5;
double mean0=278;
double sigma0=4;
double tail0=0;

double netarea1=2;
double mean1=277;
double sigma1=3;
double tail1=0;

double netarea2=3;
double mean2=278;
double sigma2=6;
double tail2=1;

double Integral;
double errIntegral;


void scangraphfit() {
	gStyle->SetOptFit(0); //no print stat
	//sprintf(Detector,"GePD");
	sprintf(Detector,"GeDD");
	//sprintf(type,"MC");
	sprintf(type,"Experimental");
	ofstream results;
	results.open("FitResults.txt",  ios::app);
	sprintf(header,"***************** Fit results %s- %s ***********************",type,Detector);
	//gStyle->SetOptFit();
   const char *data01; 
   const char *dataout; 
    if (strcmp(Detector, "GePD") == 0){
    	if (strcmp(type, "MC") == 0){
    		data01= "ScanGePDMCnonorm.txt";
    		dataout= "ScanGePDMCfitnonorm.png";
    	}	
    	else if (strcmp(type, "Experimental") == 0){
    		data01= "ScanGePDExpnonorm.txt";
    		dataout= "ScanGePDExpfitnonorm.png";
    	}	
	}
	else if (strcmp(Detector, "GeDD") == 0){
		if (strcmp(type, "MC") == 0){
    		data01= "ScanGeDDMCnonorm.txt";
    		dataout= "ScanGeDDMCfitnonorm.png";
    	}
    	else if (strcmp(type, "Experimental") == 0){
    		data01= "ScanGeDDExpnonorm.txt";
    		dataout= "ScanGeDDExpfitnonorm.png";
    	}	
	}
   TCanvas *c01 = new TCanvas("c01","multigraph",1280,1024);
  // c01->SetLogy();
	gStyle->SetOptStat(111110);
   float offx=1.3;
   float offy=1.3;
   float margr=0.08;
   float w=3;
   float margl=0.12;
   float line=2;
   gPad->SetLeftMargin(margl);
   gPad->SetRightMargin(margr);
   TGraphErrors *graph1 = new TGraphErrors(data01,"%lg %lg %lg %lg");
	
   printf("%d \n",graph1->GetN());
   graph1->SetMarkerColor(kBlue);
   graph1->SetLineColor(kBlue);
   graph1->SetMarkerStyle(22);
   graph1->SetMarkerSize(1.5);
   graph1->SetLineWidth(line);
 //  graph1->SetLineStyle(8);
   
   graph1->SetTitle("");
   graph1->GetXaxis()->SetTitle("E_{p} (keV)");
   if (strcmp(type, "MC") == 0) graph1->GetYaxis()->SetTitle("Counts");
   else if (strcmp(type, "Experimental") == 0) graph1->GetYaxis()->SetTitle("c/mC");
   graph1->GetYaxis()->SetTitleOffset(offy);
   graph1->GetXaxis()->SetTitleOffset(offx);
   graph1->GetYaxis()->SetTitleSize(40);
   graph1->GetYaxis()->SetTitleFont(43);
   graph1->GetYaxis()->SetLabelFont(43);
   graph1->GetYaxis()->SetLabelSize(40);
   graph1->GetXaxis()->SetTitleSize(40);
   graph1->GetXaxis()->SetTitleFont(43);
   graph1->GetXaxis()->SetLabelFont(43);
   graph1->GetXaxis()->SetLabelSize(40);
   TF1 *g1    = new TF1("g1","[0]*exp(-0.5*((x-[1])/([2]+[3]*exp(-(x-[1])/[1])))^2)",x1,x2);
   g1->SetParameter( 0, netarea0);
   g1->SetParameter( 1, mean0);
   g1->SetParameter( 2, sigma0);
   g1->SetParameter( 3, tail0);
   g1->SetParLimits( 0, netarea1,netarea2);
   g1->SetParLimits( 1, mean1,mean2);
   g1->SetParLimits( 2, sigma1,sigma2);
   g1->SetParLimits( 3, tail1,tail2);
   
    auto fitResult = graph1->Fit(g1, "S");
    c01->Update();
	graph1->Draw();
	auto covMatrix = fitResult->GetCovarianceMatrix();
    Integral=g1->Integral(x1, x2);
    double errIntegral = g1->IntegralError(x1,x2, fitResult->GetParams() , covMatrix.GetMatrixArray());

   graph1->Draw("AP");
   
   //graph1->SetMinimum(1.e-6);
   TLegend* leg = new TLegend(0.75, 0.73, .85, .88);
   leg->SetHeader(" ");
   leg->SetNColumns(1);
   
   leg->AddEntry(graph1, "MC", "ap");
   leg->AddEntry(g1, "Fit", "l");
   leg->SetBorderSize(0);
   leg->Draw();
   
   c01->Modified();
   c01->Update();
   c01->Print(dataout);
   
   double netarea = g1->GetParameter(0);
	double mean = g1->GetParameter(1);
	double sigma = g1->GetParameter(2);
	double tail = g1->GetParameter(3);
	double errnetarea = g1->GetParError(0);
	double errmean = g1->GetParError(1);
	double errsigma = g1->GetParError(2);
	double errtail = g1->GetParError(3);
	double chi2 = g1->GetChisquare();; //
	double ndf = g1->GetNDF(); 	
	double chi2red=chi2/ndf;

   
   results << header << endl;
   results << "Fit range: x1= " << x1 << "\t x2= " << x2 << endl;
   results << "Fit function= [0]*exp(-0.5*((x-[1])/([2]+[3]*exp(-(x-[1])/[1])))^2)" << endl;
   results << "Fit initialization net area = " << netarea0<< endl;
   results << "Fit initialization mean = " << mean0<< endl;
   results << "Fit initialization sigma = " << sigma0<< endl;
   results << "Fit initialization tail = " << tail0<< endl;
   results << "Fit net area initialization limits= " << netarea1 << "; " << netarea2 << endl;
   results << "Fit mean initialization limits= " << mean1 << "; " << mean2 << endl;
   results << "Fit sigma initialization limits= " << sigma1 << "; " << sigma2 << endl;
   results << "Fit tail initialization limits= " << tail1 << "; " << tail2 << endl;
   results << "Fit range= " << x1 << "; " << x2 << endl;
   results << "Net Area= " << netarea << " +/- " << errnetarea << endl;
	results << "Mean= " << mean <<  " +/- " << errmean << endl;
	results << "Sigma= " << sigma  <<  " +/- " << errsigma << endl;
	results << "Tail= " << tail  <<  " +/- " << errtail << endl;
	results << "chi2=  " << chi2 << endl;
	results << "chi2red=  " << chi2red << endl;
	results << "NdF=  " << ndf << endl;
	results << "Integral Area=  " << Integral << " +/- " << errIntegral << endl;
	results << "Plot=  " << dataout << endl;
	results <<  endl;
   
}