#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;

	Double_t Emin = 0; 
	Double_t Emax = 10;
	Double_t bin_size = 0.25;
	Int_t nb = (Emax-Emin)/bin_size;
	
	Double_t Bn = 8.020;
	Double_t E_thres = 0.52; 
	
	Double_t B = 1;
	Double_t Bcal = 0.0;
	
	
	TH1D *LDT = new TH1D("LDT","Total Level density",nb,Emin,Emax);
	TH1D *LDJ = new TH1D("LDJ","0<J<2 Level density",nb,Emin,Emax);
	TH1D *RSF = new TH1D("RSF","Partial gamma width via E gamma",nb,Emin,Emax);
	TH1D *rsf = new TH1D("rsf","Gamma strength function",nb,Emin,Emax);
	TH1D *gRSF = new TH1D("gRSF","Total Gamma",nb,Emin,Emax);
	TH1D *mRSF = new TH1D("mRSF","Partial gamma width for intermediate",nb,Emin,Emax);
	
	TH1D *SD = new TH1D("SD","spin distribution",nb,Emin,Emax);
	TH1D *Igg = new TH1D("Igg","Cascade intensity",nb,Emin,Emax);
	TH1D *IggC = new TH1D("IggC","Calculated Cascade intensity",nb,Emin,Emax);
	
	TH1D *LD_RIPL3 = new TH1D("LD_RIPL3","Counting LD", nb, Emin, Emax);
	TH1D *LD_RIPL3_LS = new TH1D("LD_RIPL3_LS","Counting LD J<2", nb, Emin, Emax);
	
	TH1D *tmpld = new TH1D("tmpld","tmp ld",nb,Emin,Emax);
	TH1D *tmprsf = new TH1D("tmprsf","tmp rsf",nb,Emin,Emax);
	
	TH1D *addLD = new TH1D("addLD","adding ld logarithm",nb,Emin,Emax);
	TH1D *addrsf = new TH1D("addrsf","adding rsf logarithm",nb,Emin,Emax);
	
	TH1D *simLD = new TH1D("simLD","LD simulation", nb, Emin, Emax);
	TH1D *simrsf = new TH1D("simrsf","RSF simulation", nb, Emin, Emax);
	
	
void read_Iexp()
{
ifstream if1; 
	double t[9];
	
	if1.open("unresolve.in");
	while (if1 >> t[0] >> t[1])
	{
		Igg->SetBinContent(Igg->GetXaxis()->FindBin(t[0]/1000),Igg->GetBinContent(Igg->GetXaxis()->FindBin(t[0]/1000))+t[1]);
	}
	if1.close();
	
	for (Int_t i = 1; i<=nb; i++)
	{
		Igg->SetBinError(i,0.5*Igg->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])
	{
		Igg->SetBinContent(Igg->GetXaxis()->FindBin(t[0]/1000),Igg->GetBinContent(Igg->GetXaxis()->FindBin(t[0]/1000))+t[7]);
		Igg->SetBinError(Igg->GetXaxis()->FindBin(t[0]/1000),sqrt(pow(Igg->GetBinError(Igg->GetXaxis()->FindBin(t[0]/1000)),2.)+pow(t[8],2.)));
	}
	if1.close();
	
	//Normalize to 1 relative
	Igg->Scale(1/Igg->Integral(3,30));
}


void cal_SD()
	{
	// 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;
	// Level density from RIPL3
	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);
	}
	
	//Spin distribution
	for (int i = 1; i<=nb; i++)
	{
		Double_t E = SD->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);
			}			
			SD->SetBinContent(i,sum2/sum);
		}
		else SD->SetBinContent(i,LD_RIPL3_LS->GetBinContent(i)/LD_RIPL3->GetBinContent(i));
	}
	for (int i = 1; i<=nb; i++)
	{
		if (SD->GetBinContent(i)==0) SD->SetBinContent(i,1.);
	}
	}
	

void initial_LDRSF()
{
	for (int i=1; i<=nb; i++)
	{
		LDJ->SetBinContent(i,1.0e2);
		rsf->SetBinContent(i,1.0e-8);
	}
}

void Pcal()
{
	for (int i=1; i<=nb; i++)
	{
		Double_t E = rsf->GetXaxis()->GetBinCenter(i);
		RSF->SetBinContent(i,rsf->GetBinContent(i)*E*E*E);
	}
	for (int i=1; i<=nb; i++)
	{
		Double_t E = rsf->GetXaxis()->GetBinCenter(i);
		gRSF->SetBinContent(i,RSF->Integral(1,i));
	}
	Double_t Ef[6]={0., .078, 1.043, 1.117, 1.156, 1.197};
	for (int i=1; i<=nb; i++)
	{
		Double_t E = rsf->GetXaxis()->GetBinCenter(i);
		mRSF->SetBinContent(i,0.0);
		for (int j = 0; j<6; j++)
		{
			if (Bn-E-Ef[j]>0)
			{
				Int_t kE2 = mRSF->GetXaxis()->FindBin(Bn-E-Ef[j]);
				mRSF->SetBinContent(i,mRSF->GetBinContent(i)+RSF->GetBinContent(kE2));
			}				
		}
	}
}
void cal_Igg()
{
	//Pcal();
	Int_t kBn = Igg->GetXaxis()->FindBin(Bn);
	
	for (int i = 1; i<=nb; i++)
	{
		if ((i>=3)&&(i<=30))
		{
		Double_t E = RSF->GetXaxis()->GetBinCenter(i);
		Int_t kEi = RSF->GetXaxis()->FindBin(Bn-E);
		Double_t gim,gi,gmf,gm,rho;
		gim = RSF->GetBinContent(i);
		gi = gRSF->GetBinContent(kBn);
		gmf = mRSF->GetBinContent(i);
		gm = gRSF->GetBinContent(kEi);
		rho = LDJ->GetBinContent(kEi);
		IggC->SetBinContent(i,0.25*gim*gmf*rho/gi/gm);
		}
		else
		IggC->SetBinContent(i,0.0);
	}
}

void cal_B()
{
	Double_t tmp = 0.0;
	Int_t kBn = RSF->GetXaxis()->FindBin(Bn);
	for (int i = 1; i<=kBn; i++)
	{
		Double_t E = RSF->GetXaxis()->GetBinCenter(i);
		Int_t kEi = RSF->GetXaxis()->FindBin(Bn-E);
		tmp=tmp+RSF->GetBinContent(i)*LDJ->GetBinContent(kEi)*0.25;
	}
	Bcal=tmp;
}

void sim_igg()
{
	for (int i = 1; i<=nb; i++)
	{
		Double_t E = simLD->GetXaxis()->GetBinCenter(i);
		simrsf->SetBinContent(i, pow(10.,-9)*2.5*E);
		simLD->SetBinContent(i,exp(E*1.2));
		if (E<=2) simLD->SetBinContent(i,10.);
	}
	for (int i = 1; i<=nb; i++)
	{
		LDJ->SetBinContent(i,simLD->GetBinContent(i));
		rsf->SetBinContent(i,simrsf->GetBinContent(i));
	}
	Pcal();
	cal_Igg();
	for (int i = 1; i<=nb; i++)
	{
		Igg->SetBinContent(i,IggC->GetBinContent(i));
	}
}

Double_t Chi2_I()
{
	Double_t tmp=0;
	for (int i=3; i<=30; i++)
	{
		tmp=tmp+pow(Igg->GetBinContent(i)-IggC->GetBinContent(i),2.);
	}
	return tmp;
}

Double_t Chi2_b()
{
	return (Bcal-B)*(Bcal-B);
}

Double_t gf(Double_t C, Double_t U, Double_t sig, Double_t E)
{
	return C*exp(-1*(E-U)*(E-U)/sig/sig);
}
void addhis()
{
	Double_t C1 = -0.2+gRandom->Uniform(1)*0.4;
	Double_t C2 = -0.2+gRandom->Uniform(1)*0.4;
	Double_t U1 = gRandom->Uniform(1)*Bn;
	Double_t U2 = gRandom->Uniform(1)*Bn;
	Double_t S1 = gRandom->Uniform(1)*(0.3*Bn-bin_size)+bin_size;
	Double_t S2 = gRandom->Uniform(1)*(0.3*Bn-bin_size)+bin_size;
	
	for (int i = 1; i<=nb; i++)
	{
		Double_t E = addLD->GetXaxis()->GetBinCenter(i);
		addLD->SetBinContent(i,gf(C1,U1,S1,E));
		addrsf->SetBinContent(i,gf(C2,U2,S2,E));
	}
}

void extract()
{
	initial_LDRSF();
	Pcal();
	cal_Igg();
	//cal_B();
	Double_t X2I = Chi2_I();
	Double_t X2b = Chi2_b();
	cout<<X2I<<"\t"<<X2b<<endl;
	
	for (int id =1; id<=100; id++)
	{
		for (int i=1; i<=nb; i++)
		{
			tmpld->SetBinContent(i,LDJ->GetBinContent(i));
			tmprsf->SetBinContent(i,rsf->GetBinContent(i));
		}
		
		label1: 
		if (X2I<=1.e-5) goto kt; //Diff < 0.3%
		
		addhis();
		for (int i=1; i<=nb; i++)
		{
			LDJ->SetBinContent(i,pow(10.,log10(tmpld->GetBinContent(i))+addLD->GetBinContent(i)));
			rsf->SetBinContent(i,pow(10.,log10(tmprsf->GetBinContent(i))+addrsf->GetBinContent(i)));
		}
		
		Pcal();
		cal_Igg();
		//cal_B();
		
		Double_t X2It = Chi2_I(); Double_t X2bt = Chi2_b();
		
		if (X2It<X2I)
		{X2I = X2It; cout<<"RNUmber = "<<id<<"\t"<<X2I<<endl;}
		else goto label1;
	}
	kt: ;
	//Save extracted histogram to file
	
	TFile f("kq.root","recreate");
	f.cd();
	LDJ->Write("LDJ1");
	rsf->Write("rsf1");
	f.Close();
}

Double_t FNorm(Double_t *x, Double_t *par)
{
	Int_t k = LDJ->GetXaxis()->FindBin(x[0]);
	Double_t tmp = LDJ->GetBinContent(k);
	return par[0]*exp(par[1]*x[0])*tmp;
}


void main()
{
	//read_Iexp();
	sim_igg();
	TCanvas *c1 = new TCanvas("c1","Experimental Igg");
	c1->Divide(2,2); 
	c1->cd(1); c1->cd(1)->SetLogy(); simLD->Draw();
	c1->cd(2); c1->cd(2)->SetLogy(); simrsf->Draw();
	c1->cd(3);  Igg->Draw();
	c1->cd(4);  IggC->Draw();
	//extract();
	
	TFile f("kq.root");	
	LDJ = (TH1D*)f.Get("LDJ");
    LDJ->SetDirectory(nullptr);
	LDJ->Draw();
	f.Close();
	
	
	// TCanvas *c2 = new TCanvas("c2","after extracting");
	// c2->Divide(2,2);
	// c2->cd(1); c2->cd(1)->SetLogy(); LDJ->Draw();
	// c2->cd(2); c2->cd(2)->SetLogy(); rsf->Draw();
	// c2->cd(3); Igg->Draw(); IggC->SetLineColor(kRed); IggC->Draw("same");
	
	// TF1 *ffit = new TF1("ffit",FNorm,0.5, 7.5, 2);
	// ffit->SetParameters(10,2);
	// Double_t A, alp; 
	// for (int j =1; j<=10; j++)
	// {
		// simLD->Fit("ffit","RW");
		// A = ffit->GetParameter(0);
		// alp = ffit->GetParameter(1);
		// for (int i=1; i<=nb; i++)
		// {
			// Double_t E = LDJ->GetXaxis()->GetBinCenter(i);
			// LDJ->SetBinContent(i,LDJ->GetBinContent(i)*A*exp(E*alp));
			// rsf->SetBinContent(i,rsf->GetBinContent(i)*exp(E*alp));
		// }
	// }

	
	
	
}