#include <iostream>
#include <cmath>
#include <iomanip>
#include <fstream>
#include <string>
#include <sstream>
#include <TH1.h>
#include <TFile.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TRoot.h>
#include <TGraphErrors.h>
#include <TRandom.h>

using namespace std;

	//gROOT->Reset();
	Double_t Emin = 0; 
	Double_t Emax = 10;
	Double_t bin_size = 0.25;
	Int_t nb = (Emax-Emin)/bin_size;
	
	Double_t rhoBn; //For normalization LD
	
	Double_t Eg_threshold_1 = 1.0;   //Only when extracting
	Double_t Eg_threshold = 1.0;      //For getting data
	Double_t Bn = 8.020;
	
	TH1D *I_exp = new TH1D("I_exp","Igg experiment",nb,Emin,Emax);
	TH1D *I_th  = new TH1D("I_th" ,"Igg theory"    ,nb,Emin,Emax);
	
	TH1D *LD   = new TH1D("LD"  ,"LD extract"   ,nb,Emin,Emax);    //experimental LD (0<=J<=2)
	TH1D *LD1   = new TH1D("LD1"  ,"LD1 extract"   ,nb,Emin,Emax); //Total LD
	TH1D *RSF   = new TH1D("RSF"  ,"RSF extract"   ,nb,Emin,Emax); //RSF * Eg^3
	TH1D *rsf   = new TH1D("RSF"  ,"rsf extract"   ,nb,Emin,Emax); //RSF 
	
	TH1D *LD_RIPL3 = new TH1D("LD_RIPL3","Discrete Level Density from RIPL3", nb, Emin, Emax);
	TH1D *LD_RIPL3_LS = new TH1D("LD_RIPL3_LS","Discrete Level Density from RIPL3 for spin from 0 to 2", nb, Emin, Emax);
	
	TH1D *LD_BSFG  = new TH1D("LD_BSFG","BSFG Level Density", nb, Emin, Emax);
	
	TH1D *his_gJ = new TH1D("his_gJ","gJ histogram",nb,Emin,Emax);
	
	
void read_Iexp()
{
	ifstream if1; 
	double t[9];
	
	if1.open("unresolve.in");
	while (if1 >> t[0] >> t[1])
	{
		I_exp->SetBinContent(I_exp->GetXaxis()->FindBin(t[0]/1000),I_exp->GetBinContent(I_exp->GetXaxis()->FindBin(t[0]/1000))+t[1]);
	}
	if1.close();
	
	for (Int_t i = 1; i<nb; i++)
	{
		I_exp->SetBinError(i,sqrt(I_exp->GetBinContent(i)));
	}
	
	if1.open("resolve.in");
	
	while (if1 >> t[0] >> t[1] >> t[2] >> t[3] >> t[4] >> t[5] >> t[6] >> t[7] >> t[8])
	{
		I_exp->SetBinContent(I_exp->GetXaxis()->FindBin(t[0]/1000),I_exp->GetBinContent(I_exp->GetXaxis()->FindBin(t[0]/1000))+t[7]);
		I_exp->SetBinError(I_exp->GetXaxis()->FindBin(t[0]/1000),I_exp->GetBinError(I_exp->GetXaxis()->FindBin(t[0]/1000))+t[8]);
	}
	if1.close();
	
	I_exp->Scale(1/1e6); //Normalized to 1 captures
	
	Double_t sum = I_exp->Integral();
	for (int i = 1; i<=nb; i++) {I_exp->SetBinContent(i,I_exp->GetBinContent(i)/sum);}
	
	TFile *f = new TFile("Igg.root","RECREATE");
	I_exp->Write();
	f->Close();
	
//	I_exp->Draw("");
//	I_exp->Draw("HIST SAME");
	
}


Double_t fLD(Double_t x, Double_t *par)
{
   Double_t tmp = par[LD->GetXaxis()->FindBin(x)-1];
   return tmp;
}

Double_t fRSF(Double_t x, Double_t *par)
{
   Double_t tmp = par[RSF->GetXaxis()->FindBin(x)-1];
   return tmp;
}

Double_t fIgg_ext(Double_t *x, Double_t *par)
{
	Double_t E1 = x[0];
	Double_t sum = 0.;
	for (int i = LD->GetXaxis()->FindBin(Eg_threshold_1); i<=LD->GetXaxis()->FindBin(Bn-Eg_threshold_1); i++)
	{
		Double_t E = LD->GetXaxis()->GetBinCenter(i);
		sum = sum + fLD(Bn-E,par)*fRSF(E,&par[nb]);
	}
	if ((E1>Eg_threshold_1)&&(E1<(Bn-Eg_threshold_1))) return fLD(Bn-E1,par)*fRSF(E1,&par[nb])/sum;
	else return 0;
}

void extract()
{
	Double_t ini_par = 10.0; //10.0
	TF1 *fitfunc = new TF1("fitfunc",fIgg_ext,Emin,Emax,nb*2);
	for (int i = 0; i<nb*2; i++) fitfunc->SetParameter(i,ini_par);
	//for (int i = 0; i<nb*2; i++) fitfunc->SetParLimits(i,0,1e10);
	I_exp->Fit("fitfunc","");
	Double_t *par; par = new Double_t[nb*2];
	fitfunc->GetParameters(&par[0]);
	for (int i=0; i<nb; i++)
	{
		LD->SetBinContent(i+1,par[i]);
	}
	for (int i=40; i<nb*2; i++)
	{
		RSF->SetBinContent(i-nb+1,par[i]);
	}
	//Modifiy rho + LD from bin to point
	for (int i=1; i<=nb; i++)
	{
		LD->SetBinContent(i, LD->GetBinContent(i)*bin_size);
		RSF->SetBinContent(i, RSF->GetBinContent(i)*bin_size);
	}
}

void cal_Igg()
{
	Double_t sum = 0.;
	for (int i = LD->GetXaxis()->FindBin(Eg_threshold_1); i<=LD->GetXaxis()->FindBin(Bn-Eg_threshold_1); i++)
	{
		Double_t E = LD->GetXaxis()->GetBinCenter(i);
		Int_t k_Ei = LD->GetXaxis()->FindBin(Bn-E);
		sum = sum + RSF->GetBinContent(i)*LD->GetBinContent(k_Ei);
	}
	for (int i = 1; i<=nb; i++)
	{
		Double_t E = LD->GetXaxis()->GetBinCenter(i);
		Int_t k_Ei = LD->GetXaxis()->FindBin(Bn-E);
		if ((E>Eg_threshold_1)&&(E<(Bn-Eg_threshold_1))) I_th->SetBinContent(i,LD->GetBinContent(k_Ei)*RSF->GetBinContent(i)/sum); 
		else I_th->SetBinContent(i,0);
	}
}

void gJ()
{
	// Double_t P = 0.69 + 0.68; // First Oslo pub
	Double_t P = 1.82998;     //RIPL2 = 0.417  //RIPL3 (GC model) 1.82998   0.417
	Double_t a = 20.31525;      //RIPL2 = 17.588 //RIPL3 (GC model) 20.31525  17.588
	Double_t A = 172;         //Yb 172
	Double_t sig2_172yb = 0.0888*sqrt(a*(Bn-P))*pow(A,2.0/3.0); //First Oslo
	Double_t D0 = 5.8e-6; //RIPL3 Level spacing at Bn (MeV)
	Double_t It = .5; //Spin of target nucleus Yb171;
	
	//Calculate rho(Bn) for 172Yb
	rhoBn = 2.*sig2_172yb/D0/((It+1.)*exp(-1.*pow(It+1,2)/(2.*sig2_172yb))+It*exp(-1.*It*It/2/sig2_172yb));
	cout<<"rhoBn = "<<rhoBn<<endl;
	
	//RIPL3 LD counting
	ifstream if1;
	if1.open("DLS.ripl3");
	Double_t x1, x2; 
	while (if1>>x1>>x2)
	{
		LD_RIPL3->Fill(x1); //Total
		if (abs(x2)<=2.0) LD_RIPL3_LS->Fill(x1); //Limit Spin
	}
	if1.close();
	for (int i = 1; i<=nb; i++)
	{
		LD_RIPL3->SetBinContent(i,LD_RIPL3->GetBinContent(i)/bin_size);
		LD_RIPL3_LS->SetBinContent(i,LD_RIPL3_LS->GetBinContent(i)/bin_size);
	
	}
	//End RIPL3
	
	
	//Calculation the probability of J from 0 to 2;
	for (int i = 1; i<=nb; i++)
	{
		Double_t E = his_gJ->GetXaxis()->GetBinCenter(i);
		Double_t U = E - P; 
		if (U>0)
		{		
			Double_t sig2 = 0.0888*sqrt(a*U)*pow(A,2.0/3.0);
			Double_t sum  = 0.0;
			Double_t sum2 = 0.0;
			Double_t J = 0.0;
			for (int k= 0; k<=50; k++)
			{
				J = J + k*1.0;
				sum = sum + ((2.0*J+1)/(2.0*sig2))*exp(-1.0*pow(J+.5,2)/2.0/sig2);
				if ((J==0.0)||(J==1.0)||(J==2.0)) sum2 = sum2+((2.0*J+1)/(2.0*sig2))*exp(-1.0*pow(J+.5,2)/2.0/sig2);
			}			
			his_gJ->SetBinContent(i,sum2/sum);
		}
		else his_gJ->SetBinContent(i,LD_RIPL3_LS->GetBinContent(i)/LD_RIPL3->GetBinContent(i));
	}
	
	
	//BSFG LD model
	for (int i = 1; i<=nb; i++)
	{
		Double_t sol, sig2;
		Double_t E = LD_BSFG->GetXaxis()->GetBinCenter(i);
		Double_t U = E - P; 
		
		if (U<=0) LD_BSFG->SetBinContent(i,0);
		else
		{
			sig2 = 0.0888*sqrt(a*U)*pow(A,2.0/3.0);
			sol = sqrt(3.14)*exp(2*sqrt(a*U))/12.0/pow(a,1.0/4.0)/pow(U,5.0/4.0)/sqrt(2.0*3.14)/sqrt(sig2);
			LD_BSFG->SetBinContent(i,sol);
		}	
		
	}	
	}
	
Double_t md_E1(Double_t x, Double_t *par)
{
	// 0 sigE1 ; GE1 1; EE1 2
	return 8.6e-8*0.7*par[0]*pow(par[1],2.)*(pow(x,2.)+pow(2*3.14*par[3],2.))/par[2]/(pow(pow(x,2.)-pow(par[2],2.),2.));
}

Double_t md_M1(Double_t x, Double_t *par)
{
	// 0 sigE1 ; GE1 1; EE1 2
	return 8.6e-8*par[0]*x*pow(par[1],2.)/(pow((pow(x,2.)-pow(par[2],2.)),2.)+pow(x*par[1],2));
}

Double_t md_E2(Double_t x, Double_t *par)
{
	// 0 sigE1 ; GE1 1; EE1 2
	return (5./3.)*8.6e-8*par[0]*x*pow(par[1],2.)/(pow((pow(x,2.)-pow(par[2],2.)),2.)+pow(x*par[1],2))/pow(x,2);
}

Double_t md_E1M1E2(Double_t *x, Double_t *par)
{
	return md_E1(x[0],par)+md_M1(x[0],&par[4])+pow(x[0],2.)*md_E2(x[0],&par[7]);
}

Double_t frsf(Double_t x, Double_t *par)
{
	Double_t tmp = par[rsf->GetXaxis()->FindBin(x)-1];
	//Double_t tmp = md_E1M1E2(x,par);
   return tmp;
}

Double_t FG(Double_t *x, Double_t *par)
{
	Double_t E1 = x[0];
	Double_t Ef[6]={0, 0.078, 1.043, 1.117,1.156, 1.197};
	if ((E1>Eg_threshold)&&(E1<Bn-Eg_threshold))
	{
		Double_t Gim = frsf(E1,par)*bin_size;
	
		Double_t Gi=0.; 
	for (int i=1; i<=rsf->GetXaxis()->FindBin(Bn); i++)
	{
		Double_t E = rsf->GetXaxis()->GetBinCenter(i);
		Gi = Gi + frsf(E,par)*bin_size;
	}
	
	Double_t Ei = Bn-E1;
	Double_t Gm = 0.;
	for (int i=1; i<=rsf->GetXaxis()->FindBin(Ei); i++)
	{
		Double_t E = rsf->GetXaxis()->GetBinCenter(i);
		Gm = Gm+frsf(E,par)*bin_size;
	}
	
	Double_t Gmf = 0.;
	for (int i = 0; i<6; i++)
	{
		if ((Bn-E1-Ef[i])>Eg_threshold)
		Gmf = Gmf + frsf(Bn-E1-Ef[i],par)*bin_size;
	}
	
	return (Gim*Gmf/(Gi*Gm))*par[40]; //*exp(par[41]*E1);  //1e3; 1.05
	}
	else return 0.;
	
}


Double_t norm_LD(Double_t *x, Double_t *par)
{
	Int_t i = LD->GetXaxis()->FindBin(x[0]); 
	if (x[0]>=Eg_threshold) return par[0]*LD->GetBinContent(i)*exp(par[1]*x[0])/his_gJ->GetBinContent(i); 
	else 
	return 0;
}

Double_t extra_LD1(Double_t *x, Double_t *par)
{
	Double_t sig2 = 0.0888*sqrt(par[0]*(x[0]-1.37))*pow(172.0,2.0/3.0);
	return sqrt(3.14)*exp(2*sqrt(par[0]*(x[0]-1.37)))/12.0/pow(par[0],1.0/4.0)/pow((x[0]-1.37),5.0/4.0)/sqrt(2.0*3.14)/sqrt(sig2);
}

Double_t extra_LD2(Double_t *x, Double_t *par)  //Constant temperature
{
	return (1/par[0])*exp((x[0]-par[1])/par[0]);
}

Double_t extra_RSF(Double_t *x, Double_t *par)
{
	return par[0]*pow(10,par[1]*x[0]+par[2])*pow(x[0],3);
}

Double_t E2_RSF(Double_t *x, Double_t *par)
{
	Int_t i = RSF->GetXaxis()->FindBin(x[0]);
	return RSF->GetBinContent(i)/2/3.14/(par[0]*pow(x[0],3.)+par[1]*pow(x[0],5.));
}

void cal_err()
{
	//gROOT->Reset();
	read_Iexp();
	TH1D *h_tmp1 = new TH1D("h_tmp1","histogram1",nb,Emin,Emax);
	TH1D *h_tmp2_LD = new TH1D("h_tmp2_LD","histogram2",nb,Emin,Emax);
	TH1D *h_tmp3_RSF = new TH1D("h_tmp3_RSF","histogram3",nb,Emin,Emax);
	TH1D *h_tmp4_LD = new TH1D("h_tmp4_LD","Relative Error of LD",nb,Emin,Emax);
	TH1D *h_tmp5_RSF = new TH1D("h_tmp5_RSF","Relative Error of RSF",nb,Emin,Emax);
	h_tmp1 = (TH1D*)I_exp->Clone();
	extract();
	h_tmp2_LD = (TH1D*)LD->Clone();
	h_tmp3_RSF = (TH1D*)RSF->Clone();
	

	//Modify 
	Int_t nloop=100;
	for (int n = 1; n<=nloop; n++)
	{
		cout<<"Loop number "<<n<<endl;
		for (int i = 1; i<=nb; i++)
		{
			Double_t rand_value = gRandom->Gaus(0,h_tmp1->GetBinError(i));			
			I_exp->SetBinContent(i,h_tmp1->GetBinContent(i)+rand_value);
		}
		extract();
		for (int i = 1; i<=nb; i++)
		{
			Double_t ld_dif2 = pow(LD->GetBinContent(i)-h_tmp2_LD->GetBinContent(i),2.);
			Double_t rsf_dif2 = pow(RSF->GetBinContent(i)-h_tmp3_RSF->GetBinContent(i),2.);
			//cout<<sqrt(rsf_dif2)<<"\t"<<h_tmp3_RSF->GetBinContent(i)<<endl<<"\t"<<RSF->GetBinContent(i)<<endl;
			h_tmp4_LD->SetBinContent(i,h_tmp4_LD->GetBinContent(i)+ld_dif2);
			h_tmp5_RSF->SetBinContent(i,h_tmp5_RSF->GetBinContent(i)+rsf_dif2);
		}
	}
	
	for (int i = 1; i<=nb; i++)
		{
			h_tmp4_LD->SetBinContent(i,sqrt(h_tmp4_LD->GetBinContent(i))/sqrt(nloop*1.0)/abs(LD->GetBinContent(i)));
			h_tmp5_RSF->SetBinContent(i,sqrt(h_tmp5_RSF->GetBinContent(i))/sqrt(nloop*1.0)/abs(RSF->GetBinContent(i)));
		}			

		//h_tmp4_LD->GetXaxis()->SetRangeUser(0.5,7.5);
		//h_tmp4_LD->Draw();
		
	TFile *f = new TFile("Exp_result.root","RECREATE");
	h_tmp4_LD->Write();
	h_tmp5_RSF->Write();
	f->Close();		
}


void main()
{
	gROOT->Reset();
	//Read experimental data and extract LD and RSF
	read_Iexp();
	extract();
	TCanvas *c1 = new TCanvas("c1","Experimental data and first extract order");
	c1->Divide(2,2);
	c1->cd(1); I_exp->Draw("HIST");
	c1->cd(2); c1->cd(2)->SetLogy();
	LD->Draw();
	c1->cd(3); c1->cd(3)->SetLogy();
	RSF->Draw();
	cal_Igg(); //Reconstruct Igg from extracted LD and RSF
	c1->cd(4);
	I_th->Draw();
	// //End First step
	
	// //Normalize LD
	gJ(); //Prepare Data for Normalize;
	TCanvas *c2 = new TCanvas("c2","BSFG, RIPl3 and Spin distribution");
	c2->Divide(1,2);
	c2->cd(1); c2->cd(1)->SetLogy();
	LD_BSFG->Draw();
	LD_BSFG->SetLineColor(kGreen);
	LD_RIPL3->Draw("same");
	c2->cd(2);
	his_gJ->Draw();
	
	//Normalization process
	Double_t A, alpha, sv_alpha=0;
	TCanvas *c3 = new TCanvas("c3","Normalization");
	c3->cd(); c3->Divide(2,2); c3->cd(1); c3->cd(1)->SetLogy();
	for (int j = 1; j<=2; j++)	
	{
		// First fit with BSFG
	TF1 *fitf = new TF1("fitf",norm_LD,3.0,7.0,2);//5
	fitf->SetParameters(1,0);
	LD_BSFG->Fit("fitf","R0");
	A = fitf->GetParameter(0);
	alpha = fitf->GetParameter(1);
	sv_alpha=sv_alpha+alpha;
	for (int i = 1; i<=nb; i++)
	{
		Double_t E = LD->GetXaxis()->GetBinCenter(i);
		if (E<Eg_threshold) 
		{
			LD->SetBinContent(i,LD_RIPL3_LS->GetBinContent(i));
			LD1->SetBinContent(i,LD_RIPL3->GetBinContent(i));
			RSF->SetBinContent(i,0);
		};
		if ((E>=Eg_threshold)&&(E<=Bn-Eg_threshold))
		{
			LD->SetBinContent(i,A*LD->GetBinContent(i)*exp(alpha*E));
			LD1->SetBinContent(i,LD->GetBinContent(i)/his_gJ->GetBinContent(i));
			RSF->SetBinContent(i,RSF->GetBinContent(i)*exp(alpha*E));
		};
		if (E>Bn-Eg_threshold)
		{
			LD->SetBinContent(i,LD_BSFG->GetBinContent(i)*his_gJ->GetBinContent(i));
			LD1->SetBinContent(i,LD->GetBinContent(i)/his_gJ->GetBinContent(i));
			RSF->SetBinContent(i, 0);
		};
	}
	}

	// Fit with discrete data 
	for (int j = 1; j<=10; j++)
	{
	Double_t ELDthres = 2.;
	TF1 *fitf1 = new TF1("fitf1",norm_LD,1,ELDthres,2);
	fitf1->SetParameters(1,0);
	LD_RIPL3->Fit("fitf1","R0");
	A = fitf1->GetParameter(0);
	alpha = fitf1->GetParameter(1);
	sv_alpha=sv_alpha+alpha;
	for (int i = RSF->GetXaxis()->FindBin(Eg_threshold); i<=LD->GetXaxis()->FindBin(Bn); i++)
	{
		Double_t E = LD->GetXaxis()->GetBinCenter(i);
		LD->SetBinContent(i,A*LD->GetBinContent(i)*exp(alpha*E));
		LD1->SetBinContent(i,LD->GetBinContent(i)/his_gJ->GetBinContent(i));
		RSF->SetBinContent(i, RSF->GetBinContent(i)*exp(alpha*E));
	}
	}

		
	// Fit with point
	Double_t naE[5], naR[5];
	for (int i=0; i<4; i++) 
	{
		naE[i]=LD_RIPL3->GetXaxis()->GetBinCenter(i+5);
		naR[i]=LD_RIPL3->GetBinContent(i+5);
	}
	naE[4]=Bn; naR[4]=rhoBn;
	TGraph *grNorm = new TGraph(5,naE,naR); 

	for (int j=1; j<=100; j++)	
	{
	
	TF1 *fitf4 = new TF1("fitf4",norm_LD,1,8.0,2);
	fitf4->SetParameters(1,0);
	grNorm->Fit("fitf4");
	A = fitf4->GetParameter(0);
	alpha = fitf4->GetParameter(1);
	sv_alpha = sv_alpha;
	for (int i = RSF->GetXaxis()->FindBin(Eg_threshold); i<=LD->GetXaxis()->FindBin(Bn); i++)
	{
		Double_t E = LD->GetXaxis()->GetBinCenter(i);
		LD->SetBinContent(i,A*LD->GetBinContent(i)*exp(alpha*E));
		LD1->SetBinContent(i,LD->GetBinContent(i)/his_gJ->GetBinContent(i));
		RSF->SetBinContent(i, RSF->GetBinContent(i)*exp(alpha*E));
	}
	}

	LD_BSFG->Draw();
	LD1->Draw("same");
	LD_RIPL3->Draw("same");
	
	//Check LD vs Oslo
	//Read Data Oslo
	Double_t x[100], y[100], ex[100], ey[100];
	Double_t a1, a2, a3, a4;
	Int_t n;
	ifstream if1;
	if1.open("rho.oslo");
	n = 0; 
	while (if1 >> a1 >> a2 >> a3 >> a4)
		{
			x[n] = a2; y[n]=a3; ex[n]=0; ey[n]=a4; n++;
		}
	if1.close();
	TGraphErrors *LD_oslo = new TGraphErrors(n,x,y,ex,ey);
	
	if1.open("rsf.oslo");
	n = 0; 
	while (if1 >> a1 >> a2 >> a3 >> a4)
		{
			x[n] = a2; y[n]=a3; ex[n]=0; ey[n]=a4; n++;
		}
	if1.close();
	TGraphErrors *RSF_oslo = new TGraphErrors(n,x,y,ex,ey);
	//End Read Data Oslo
	
	//Compare this work with Oslo
	c3->cd(2);c3->cd(2)->SetLogy();
	for (int i=1; i<nb; i++) {LD1->SetBinError(i,LD1->GetBinContent(i)*.45);}
	LD1->Draw();
	LD_oslo->Draw("same");
	//End Compare
cal_Igg();	
	//ReCreate RSF from model strength function
	// TF1 *FGFit = new TF1("FGfit",FG,1,7,11);
	// Double_t ipar[11] = {302,4.8,15.5,0.35,1.76,4,7.5,6.8,4.05,11.33,10};
	// FGFit->SetParameters(ipar);
	// c3->cd(3);//c3->cd(3)->SetLogy();
	// FGFit->Draw();
	

	//Fitting for search rsf
	
	c3->cd(3);c3->cd(3)->SetLogy();
	
	TF1 *fmd_E1M1E2 = new TF1("fmd_E1M1E2",md_E1M1E2,0,10,10);
	fmd_E1M1E2->SetParameters(302,4.8,15.5,0.35,1.76,4,7.5,6.8,4.05,11.33);
	TF1 *FGfit = new TF1("FGfit",FG,1,7,nb);
	
	for (int i = 1; i<=nb; i++) FGfit->SetParameter(i-1,fmd_E1M1E2->Eval(rsf->GetXaxis()->GetBinCenter(i))*pow(rsf->GetXaxis()->GetBinCenter(i),3.));
	FGfit->SetParameter(40,1e3);
	
	// for (int k = 1; k<10; k++)
	// {
		// for (int i = 1; i<=nb; i++) FGfit->SetParLimits(i-1,0,1e7);
		// //FGfit->SetParLimits(40,100,2000);
		// //FGfit->SetParLimits(41,2.3777,2.4);
		// RSF->Fit("FGfit");
	// }
	
	//Double_t fgpar[40];
	//FGfit->GetParameters(fgpar);
	//for (int i=1;i<=nb;i++) {rsf->SetBinContent(i,fgpar[i-1]/pow(0.125+(i-1)*0.25,3.));}
	
	//rsf->Draw();
	FGfit->Draw();
	cout<<"alpha "<<sv_alpha<<endl;
	//for (int i = 1; i<=nb; i++) FGfit->SetParameter(i-1,1.e10*fmd_E1M1E2->Eval(rsf->GetXaxis()->GetBinCenter(i)));
	
	//for (int i = 1; i<=nb; i++) FGfit->SetParameter(i-1,1.e10*fmd_E1M1E2->Eval(rsf->GetXaxis()->GetBinCenter(i)));
	//for (int i = 1; i<=nb; i++) FGfit->SetParameter(i-1,1);
	//for (int i = 1; i<=nb; i++) FGfit->SetParLimits(i-1,1e-11,1e-5);
	//for (int i = 1; i<=nb; i++) cout<<i<<"\t"<<FGfit->Eval(rsf->GetXaxis()->GetBinCenter(i))<<endl;
	//Double_t fgpar[41];
	//for (int k =1; k<40; k++)
	//{
	//I_exp->Fit("FGfit","R");
	//FGfit->GetParameters(fgpar);
	//FGfit->SetParameter(40,fgpar[40]);
	//cout<<"fit number "<<k<<" ============= "<<FGfit->GetChisquare()<<endl;
	//}
	//for (int i = 1; i<=nb; i++) rsf->SetBinContent(i,fgpar[i-1]);
	
	
	//for (int i = 1; i<=nb; i++) rsf->SetBinContent(i,rsf->GetBinContent(i)/pow(0.125+(i-1)*0.25,3.));

	//rsf->Draw();
	//FGfit->Draw("same");
	
	//c1->cd(3);
	//FGfit->Draw("same");
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	//Double_t fgpar[40];
	//FGfit->GetParameters(fgpar);
	//for (int i = 1; i<nb; i++) rsf->SetBinContent(i,fgpar[i-1]);
	//FGfit->Draw();

// //Extrapolation for RSF
	// TF1 *exfit3 = new TF1("exfit3",extra_RSF,Eg_threshold,Bn-Eg_threshold,2);
	// exfit3->SetParameters(1,1);
	// Double_t tmp_par[2];
	// RSF->Fit("exfit3","R0");
	// exfit3->GetParameters(tmp_par);
	// TF1 *exfit3p = new TF1("exfit3p",extra_RSF,Emin,Emax,2);
	// exfit3p->SetParameters(tmp_par);
	// for (int i = 1; i<=nb; i++)
	// {
		// Double_t E = RSF->GetXaxis()->GetBinCenter(i);
		// if ((E<Eg_threshold)||(E>(Bn-Eg_threshold))) RSF->SetBinContent(i,exfit3->Eval(E));
	// }
	// c3->cd(2); c3->cd(2)->SetLogy();
	// RSF->Draw();
	// exfit3p->Draw("same");
	
	
	// // //Determine B factor;
	// Double_t D0 = 5.8e-6; //RIPL3
	// Double_t AG = 0;
	// for (int i=RSF->GetXaxis()->FindBin(0.01); i<=RSF->GetXaxis()->FindBin(Bn); i++)
	// {AG = AG + RSF->GetBinContent(i)*LD1->GetBinContent(LD1->GetXaxis()->FindBin(Bn-LD1->GetXaxis()->GetBinCenter(i)))*bin_size*his_gJ->GetBinContent(i);}
	// Double_t GgA = 75e-9;     //Oslo 2004: 75 meV (10)
	// B = GgA*4*3.14/D0/AG;
	// cout<<"B = "<<B<<endl;
	
	// for (int i = 1; i<=nb; i++)
	// {
		// Double_t E = LD1->GetXaxis()->GetBinCenter(i);
		// RSF->SetBinContent(i,B*RSF->GetBinContent(i));
		// rsf->SetBinContent(i,RSF->GetBinContent(i)/2/3.14/pow(E,3.));
	// }
	
	
	
	// ofstream of2; 
// of2.open("Test.txt");
// for (int i = 1; i<=nb; i++) {of2<<(double(i-1)*0.25+0.125)<<"\t"<<rsf->GetBinContent(i)<<endl;}

// of2<<endl; of2<<endl;
// of2<<8.020<<"\t"<<rhoBn<<endl;
// of2.close();

	
	
	// c3->cd(3); c3->cd(3)->SetLogy();
	// for (int i=1; i<=nb; i++) {rsf->SetBinError(i,rsf->GetBinContent(i)*.45);}
	// rsf->Draw();
	// RSF_oslo->Draw("same");
	
	// c3->cd(4);c3->cd(4)->SetLogy();
	// for (int i=1; i<nb; i++) {LD1->SetBinError(i,LD1->GetBinContent(i)*.45);}
	// LD1->Draw();
	// LD_oslo->Draw("same");
	
	// TF1 *E2fit = new TF1("E2fit",E2_RSF,1,7,2);
	// E2fit->SetParameters(1,0);
	// RSF_oslo->Fit("E2fit");
		
	// of1.close();
	
	
	// //Fit RSF with model
	// TCanvas *c4 = new TCanvas("c4","Fit RSF with model");
	// c4->cd()->SetLogy();
	
	// rsf->Draw();
	
	// TF1 *fmd_E1 = new TF1("fmd_E1",md_E1,1,7,4);
	// fmd_E1->SetParameters(239,2.6,12.25,0.35);
	// //fmd_E1->Draw("same");
	
	// TF1 *fmd_M1 = new TF1("fmd_M1",md_M1,1,7,3);
	// fmd_M1->SetParameters(1.76,4,7.5);
	// //fmd_M1->Draw("same");
	
	// TF1 *fmd_E2 = new TF1("fmd_E2",md_E2,1,7,3);
	// fmd_E2->SetParameters(6.8,4.05,11.33);
	// //fmd_E2->Draw("same");
	
	// TF1 *fmd_E1M1E2 = new TF1("fmd_E1M1E2",md_E1M1E2,1,7,10);
	// Double_t ipar[10] = {302,4.8,15.5,0.35,1.76,4,7.5,6.8,4.05,11.33};
	// fmd_E1M1E2->SetParameters(302,4.8,15.5,0.35,1.76,4,7.5,6.8,4.05,11.33);
	// fmd_E1M1E2->Draw("same");
	// //for (int i = 0; i<7; i++) {fmd_E1M1E2->FixParameter(i,ipar[i]);}
	// //for (int i = 7; i<10; i++) {fmd_E1M1E2->SetParLimits(i,0,10000);}
	// // for (int i = 0; i<10; i++) 
	// // {
		// // if ((i!=0)&&(i!=4)&&(i!=7))
		// // fmd_E1M1E2->FixParameter(i,ipar[i]);
		// // else 
			// // fmd_E1M1E2->SetParLimits(i,0,1000);
	// // }
	
	
	// //rsf->Fit("fmd_E1M1E2","RB");
	
	
	
}