//////////////////////////////////////////////////////////////////////////////////////////////////////
///
/// Anthony Hillairet, October 2006
///
//////////////////////////////////////////////////////////////////////////////////////////////////////

// Include de C++
using namespace std;

#include <iostream>
#include <stdlib.h>
#include <string.h>
// #include <iomanip>
#include <fstream>
#include <sstream>

// ROOT Includes
#include "TROOT.h"
#include "TApplication.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TF1.h"
#include "TMath.h"
#include "TGraphErrors.h"
#include "TLatex.h"
#include "TPaveText.h"
#include "TPaveLabel.h"

// GSL Includes
#include "gsl/gsl_cdf.h"

#define P_SFC			29.79
// #define SIGMASCALE		8.0		// To have a reasonable chisquare
#define SIGMASCALE		1.0		// Disabled
#define STARTB1			879.0	// Arbitraly chosen start for B1NMR at 29.79MeV/C
#define POW_BKGD		2.5		// Pow_bckg
#define P_PRECISION		0.00001	// Maximal difference between P_SFC and P_edge

Double_t Background (Double_t *x, Double_t *par)
{
	return ( pow( x[0]/ par[0], par[1] ) * par[2] );
}

Double_t Signal (Double_t *x, Double_t *par)
{
	Double_t Part1 = pow( x[0]/ par[4], par[0] );
	Double_t Part2 = 1 - gsl_cdf_ugaussian_P((x[0] - par[2] ) / par[3]);
	return ( Part1 * Part2 * par[1] );
}

Double_t FitFunction(Double_t *x, Double_t *par)
{
	return ( Signal(x, par) + Background(x, &par[4]) );
} 

int main(int argc, char **argv)
{
    gStyle	->SetPadBorderMode	( 0			);
    gStyle	->SetOptTitle		( kTRUE		);	// Display histogram titles.
    gStyle	->SetCanvasColor	( 10		);	// Background color (white)
    gStyle	->SetPadColor		( 10		);	// Background color (white)
	gStyle	->SetFillColor		( 10		);	// Background color (white)
    gStyle	->SetTitleFillColor	( 10		);	// Background color (white)
	gStyle	->SetPaperSize		( 20., 24.	);
	gROOT	->SetBatch();
	
	string InFile	= argv[1];
	string Title	= argv[2];
	string TmpOutFile	= argv[3];
	string OutFile	= TmpOutFile + ".ps";
	string RootFile	= TmpOutFile + ".root";

	TApplication TheApp("App",&argc,argv);
	TCanvas	*Canv = new TCanvas("Canv","Canv");

	int pts	= 0;
	Double_t 	P	[100];
	Double_t 	Run	[100];
	Double_t 	B1	[100];
	Double_t 	B2	[100];
	Double_t	M	[100];
	Double_t	T1	[100];
	Double_t 	MoT	[100];
	Double_t 	yErr[100];

	// Open the ascii file
	ifstream Input(InFile.c_str());
	if ( !Input )
	{
		cerr<<" *** ERROR *** ASCII data file can't be opened"<<endl;
		return 1;
	}	
	// read the tab in arrays
	string line;
	while ( getline( Input, line ) )
	{
		if ( line[0] == '#' )
		{
			continue;
		}
		istringstream iss(line);
		iss>>P[pts]>>Run[pts]>>B1[pts]>>B2[pts]>>M[pts]>>T1[pts]>>MoT[pts];
		//  assumes no uncertainty in t1ion (it is _not_ statistically Poisson)
		if ( MoT[pts] == -1 )
		{
			cout<<" *** WARNING *** M/T1ion == -1  => Point skipped"<<endl;
			continue;
		}
		yErr[pts]	= SIGMASCALE * MoT[pts] * sqrt( 1/(M[pts]) );
		pts++;
	}
	Input.close();
	
	// define the function
	TF1 *FitFunc	= new TF1("FitFunc", FitFunction,28.,31.,7);

	// set the parameters for the first fit
	FitFunc		->SetParameter(0, 1.0		);		// Pow_mu
	FitFunc		->SetParameter(1, 1.0		);		// Rate
	FitFunc		->SetParameter(2, P_SFC		);		// P_edge
	FitFunc		->SetParameter(3, 1.0		);		// W_edge
	FitFunc		->FixParameter(4, P_SFC		);		// P_SFC
	FitFunc		->FixParameter(5, POW_BKGD	);		// Pow_bkgd
	FitFunc		->SetParameter(6, 1.0		);		// bkgd

	Double_t NewB1	= STARTB1;
	bool Done		= false;
	Double_t P_edge;

	// The goal is to match P_edge and P_SFC
	// The loop is performed until the precision is satisfying
	while ( ! Done )
	{
		for ( int i = 0; i < pts; i++ )
		{
			P[i]= P_SFC * B1[i] / NewB1;
		}
		TGraphErrors *G0	= new TGraphErrors(pts, P, MoT, 0, yErr);
		G0	->Fit ("FitFunc","0Q");
		P_edge	= FitFunc	->GetParameter(2);
		if ( TMath::Abs( P_SFC - P_edge ) > P_PRECISION )
		{
			// B1 variation is function of the difference P_SFC - P_edge to go quickly to the good value
			// It seems weird to mix momentum and magnetic field but it works very well
			NewB1 = NewB1 - ((P_SFC - P_edge));
		}
		else
			Done = true;
	}
	Double_t B1NMR = NewB1 * 29.6/P_SFC;
	Double_t MomBite = 100. * 2.355 * FitFunc->GetParameter(3) / P_SFC;
	
	// TGraphErrors used for the display
	TFile *RootOut = new TFile(RootFile.c_str(), "RECREATE");
	TGraphErrors *G	= new TGraphErrors(pts, P, MoT, 0, yErr);
	G	->Fit ("FitFunc","Q");
	
	Canv->	cd();
	Canv->	cd();

	// Fitted equation printed
	TPaveText *Tpt	= new TPaveText(0.6,0.7,0.97,0.9,"r");
	Tpt	->	SetTextAlign(12);
	Tpt	->	SetTextSize(0.025);
	Tpt	->	SetBorderSize(1);
	Tpt	->	AddText(0.05, 0.73, "y = #left(#frac{x}{P_{sfc}}#right)^{U} #times Rate #times #left[1 - Gauss_{cdf}#left(#frac{x-P_{edge}}{W_{edge}}#right)#right]");
	Tpt	->	AddText(0.05, 0.23, "                                           + #left(#frac{x}{P_{sfc}}#right)^{B} #times Bkgd");

	// Parameters printed
	char* Text	= new char[100];
	TPaveText *Tpt2	= new TPaveText(0.7,0.22,0.97,0.7,"r");
	Tpt2->	SetTextAlign(12);
	Tpt2->	SetTextSize(0.025);
	Tpt2->	SetBorderSize(1);
	sprintf (Text, " #frac{#chi^{2}}{dof}    = %1.3f", FitFunc->GetChisquare()/FitFunc->GetNDF());
	Tpt2->	AddText(0.05, 0.93, Text);
	sprintf (Text, " P_{sfc}    = %2.2f (Fixed)", FitFunc->GetParameter(4));
	Tpt2->	AddText(0.05, 0.79, Text);
	sprintf (Text, " U       = %1.3f #pm %1.3f", FitFunc->GetParameter(0), FitFunc->GetParError(0));
	Tpt2->	AddText(0.05, 0.69, Text);
	sprintf (Text, " Rate  = %1.4f #pm %1.4f", FitFunc->GetParameter(1), FitFunc->GetParError(1));
	Tpt2->	AddText(0.05, 0.59, Text);
	sprintf (Text, " P_{edge}  = %2.4f #pm %1.4f", FitFunc->GetParameter(2), FitFunc->GetParError(2));
	Tpt2->	AddText(0.05, 0.49, Text);
	sprintf (Text, " W_{edge} = %1.4f #pm %1.4f", FitFunc->GetParameter(3), FitFunc->GetParError(3));
	Tpt2->	AddText(0.05, 0.39, Text);
	sprintf (Text, " B         = %1.1f (Fixed)", FitFunc->GetParameter(5));
	Tpt2->	AddText(0.05, 0.29, Text);
	sprintf (Text, " Bkgd  = %1.4e #pm %1.2e", FitFunc->GetParameter(6), FitFunc->GetParError(6));
	Tpt2->	AddText(0.05, 0.19, Text);
	sprintf (Text, " #frac{dp}{p}       = %1.4f %%", MomBite);
	Tpt2->	AddText(0.05, 0.08, Text);

	// New B1NMR value printed
	TPaveText *Tpt3	= new TPaveText(0.15,0.2,0.4,0.4,"r");
	Tpt3->	SetTextAlign(12);
	Tpt3->	SetTextSize(0.04);
	Tpt3->	SetBorderSize(0);
	Tpt3->	AddText(0.05, 0.8, "New B1NMR value");
	Tpt3->	AddText(0.05, 0.55, "for 29.6MeV/c :");
	sprintf (Text, "B1  = %3.2f G", B1NMR);
	Tpt3->	AddText(0.05, 0.2, Text);
	cout<<endl<<" **********************************"<<endl;
	cout<<"   New B1NMR value for 29.6MeV/c"<<endl<<"   ==> "<<Text<<endl;
	cout<<" **********************************"<<endl<<endl;

	Tpt	->	Draw();
	Tpt2->	Draw();
	Tpt3->	Draw();
	G	->	Draw("AP");
	G	->	SetTitle(Title.c_str());
	Tpt	->	Draw();
	Tpt2->	Draw();
	Tpt3->	Draw();

	G->Write();
	RootOut->Close();

	Canv->	Print(OutFile.c_str());
}