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

using namespace std;

double C = 1e-2;

double LD(double E)
{
	
	//E = E/1000;
	double a = 17.0;
	double Umin = 9/(4*a);
	double U = E - Umin;
	if (E-Umin<=0) return 2;
	if (E-Umin>0) return C*pow(U,-3/2)*exp(2*sqrt(a*U));
}

double F(double E)
{
	//E = E/1000;
	return C*pow(E,4.2);
}

void main()
{
	//Creat simulation G matrix, called HG
	const Int_t nb1 = 100, nb2 = 100;
	double Ex_min = 3; 
	double Eg_min = 1;
	double err = 0.05;
	TH2D *HG = new TH2D("HG","G matrix",nb1,0,10,nb2,0,10);
	for (int i = HG->GetYaxis()->FindBin(3);i<=HG->GetYaxis()->FindBin(8); i++)
		for (int j = HG->GetXaxis()->FindBin(1); j<=HG->GetXaxis()->FindBin(8); j++)
		{
			double Ex = HG->GetYaxis()->GetBinCenter(i);
			double Eg = HG->GetXaxis()->GetBinCenter(j);
			if ((i==81)&&(j==81))
			{
				double G = LD(Ex-Eg)*F(Eg);
				double sum = 0;
				for (int k=HG->GetXaxis()->FindBin(Eg_min); k<=i; k++)
				{
					double Eg1 = HG->GetXaxis()->GetBinCenter(k);
					sum = sum + LD(Ex-Eg1)*F(Eg1);
				}
				//HG->SetBinContent(j,i,G/sum);
				//HG->SetBinError(j,i,0.05*G/sum);
				//if (HG->GetBinContent(j,i)==0) cout<<i<<"\t"<<j<<endl;
				//double tmp1 = LD(Ex-Eg);
				//double tmp2 = F(Eg);
				cout<<sum<<endl;
				cout<<LD(Ex-Eg)<<endl;
				cout<<F(Eg)<<endl;
				cout<<LD(Ex-Eg)*F(Eg)<<endl;
			}
		}
	HG->Draw("COLZ");
 f = new TFile("test.root","RECREATE");
 HG->Write();
}