#include<iostream>
#include<cmath>
#include<fstream>
#include<string>
#include<TCanvas.h>
#include <TFile.h>
#include<TApplication.h>
#include<TGraphErrors.h>
#include<TH1D.h>
#include<TF1.h>
#include<TStyle.h>
#include<TAxis.h>
#include<TLatex.h>
#include<TMath.h>
#include<TTree.h>
#include<cstdlib>
#include<stdlib.h>
#include<TRandom.h>
#include "TH2F.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TSpectrum.h"
#include<cstring>
#include<libgen.h>   //serve per basename
#include<sstream>  // serve per stringstream
#include<iomanip>
#include <stdio.h>

using namespace std;

// GLOBAL values /////////////////////////////////////////////////
string home = string(getenv("HOME"));
string path="Disang/";
string FileZiegler;
double Zgl[10][9];
const int Nelem = 2;
int Zt[Nelem] = {1, 40};
int At[Nelem] = {2, 91};
double w[Nelem] = {2, 1};
double Zp = 1;
double Ap = 1;
double Dpercm2 = 2.26e18;
//double Dpercm2 = 2.05e18;
double mp = 938.272034*1e3;
double md = 1875.612845*1e3;
double mu = mp*md/(mp+md);
double alpha = pow(137, -1);
double pi = acos(-1);
double tau;
// ==>>  mu/Tcm == mp/T
//double Exponent = 2*pi*alpha*Zp*Zt[0]*sqrt(mu/2);
double Exponent = 2*pi*alpha*Zp*Zt[0]*sqrt(mp/2);
//double Exponent = 25.6464;
//const int Npar = 3;
string DirDat = "Disang/";

double mass(int Z, int A, int flag) {
	int z, a;
  	double amu, massat, massnucl, exc;
  	double answ;
	string nucleus;
  	int check = 0;
  	ifstream in;
  	in.open("tabmass2021.dat");
  	while (!in.eof()) {
    		in>>nucleus>>z>>a>>exc>>amu>>massat>>massnucl;
    		if(z==Z && a==A) {
      			if(flag == 0) answ = amu;
      			if(flag == 1) answ = massat;
      			if(flag == 2) answ = massnucl;
      			check = 1;
      			break;
    		}
    		if(check == 1) break;
  	}
  	return answ;
}

///////////////////////////////////////////////////////////////////
// function stop
// input values: Z, A of proj,  Z, A of target, Tp of proj (keV)
// returns stopping cross section in units  eV/(10^15 atoms/cm^2
///////////////////////////////////////////////////////////////////
double stop(int Zp, int Ap, int Zt, int At, int J, double Tp) {
  	//warning: Tp must be in keV
  	double ziegler[9];
  	for (int j=0; j<9; j++) ziegler[j] = Zgl[J][j];
  	double epsil, Kp, L, bet2, Slow, Shigh;
  	if(Zp == 1) {                          // Hydrogen-like projectile..
    		Kp=Tp/Ap;
    		L=log(Kp);
    		bet2 = 1-1/pow((1+Kp/931494.062),2);
    		if(Kp<10) {
    	 		epsil = ziegler[0]*sqrt(Kp);
    		} else if(10 <= Kp < 1e3) {
      			Slow = ziegler[1]*pow(Kp,0.45);
      			Shigh = ziegler[2]/Kp*log(1+ziegler[3]/Kp+ziegler[4]*Kp);
      			epsil = Slow*Shigh/(Slow+Shigh);
		} else {
      			epsil = ziegler[5]/bet2*(log(ziegler[5]*bet2/(1-bet2))-bet2);
    		}
  	}
  	if(Zp>1) {                             // Helium-like projectile....
  		// Kp must be given in MeV...
  		// Kp is scaled with masses (compared to alpha)
  		// epsil is scaled with z^2 ratios (compared to alpha)
     		Kp = Tp*mass(2,4,1)/Ap*1e-03;
		
     		if(Kp < 1e1) {
       			Slow = ziegler[0]*pow((Kp*1e03),ziegler[1]);
       			Shigh = ziegler[2]/Kp*log(1+ziegler[3]/Kp+ziegler[4]*Kp);
       			epsil = Slow*Shigh/(Slow+Shigh)*pow((Zp/2),2); 
		} else {
        		L = log(1/Kp);
       			epsil = exp(ziegler[5]+ziegler[6]*L+ziegler[7]*pow(L,2)+ziegler[8]*pow(L,3))*pow((Zp/2),2);
     		}
  	}
  	// epsil is given in  eV/(10^15 atoms/cm^2)
  	return(epsil); 
}

double stopeff(double Tp) {
	double eps = (2*stop(Zp,Ap,Zt[0],At[0],0,Tp) + stop(Zp,Ap,Zt[1],At[1],1,Tp))/3;
	eps = eps*1e-15*1e-3;              // keV*cm2
	return eps;
}

///////////////////////////////////////////////////////////////////
// function dedx
// Tp in keV
// summing stopping over compund elements
// converting  eV/(10^15 atoms/cm^2) into keV*cm2/g
///////////////////////////////////////////////////////////////////
 double dedx(int Zp, int Ap, int Nelem, double Tp) {
	double SP = 0;
  	double SP_parz;
  	double PM = 0;
  	for (int i=0; i<Nelem; i++) {
    		SP_parz = stop(Zp, Ap, Zt[i], At[i], i, Tp);
    		SP = SP + SP_parz*w[i]; 
    		PM = PM + At[i]*w[i];
  	}
  	SP = SP*6.023e5/PM;
  	return(SP);
}

///////////////////////////////////////////////////////////////////
// function range
// computing range [cm2/g] in compound
// Tp in keV
// dedx in keV*cm2/g
///////////////////////////////////////////////////////////////////
double Range(int Zp, int Ap, int Nelem, double Tp) {
	double Tmin = 0.01;
  	double dT = Tmin;
  	double T = Tp;
  	double R = 0;
  	while (T>=Tmin) {
  		R = R +  1/dedx(Zp, Ap, Nelem, T);
		T = T - dT;
  	}
  	R = R*dT;
  	return R;
}

///////////////////////////////////////////////////////////////////
// function out
// computing resdual energy in compound
// Tp in keV
// dedx in keV*cm2/g
///////////////////////////////////////////////////////////////////
double Tfin(int Zp, int Ap, int Nelem, double Tp) {
	double dx = tau*1e-3;
	double x = dx;
	double T = Tp;
	while (x<tau) {
		  T = T - dedx(Zp, Ap, Nelem, T)*dx;
		  x = x + dx;		  //cout<<tau<<"   "<<x<<"   "<<T<<endl;
	}
	double Tout = T;
	return Tout;
}

double F(double T) {
//  usiamo sempre Tlab
	double XX = exp(-Exponent/sqrt(T))/(T*stopeff(T));
	return XX;
}

double Legendre (int n, double x) {
	const int nmax = 10;
	double p[nmax];
	int k;
	p[0] = 1;
	if(n==0) return p[0];
	p[1] = x;
	if(n==1) return p[1];
	k = 2;
	while (k<nmax) {
		p[k] = (2.0*k-1.0)/k*x*p[k-1]-(k-1.0)/k*p[k-2];
		if(n == k) break;
		k++;	
	}
	return p[k];
}

double fitfun(double *x, double *par, int Npar) {
  	double xx = x[0];
	double Sum = 0;
	for (int k=1; k<Npar; k++) {
		Sum = Sum + par[k]*Legendre(k, xx);
	}
  	return par[0]*(1 + Sum);
}

/*double disang(double x, double a1, double a2) {
	double W = 1 + a1*Legendre(1,x) + a2*Legendre(2,x);

	return W;
}*/

double disang(double x, int nPar, double *a) {
  double W = 1.0;
  for (int j = 1; j <= nPar; j++) {
    W += a[j] * Legendre(j, x);
  }
  return W;
}

void Dis(int Npar) {
	string targ_state = "solid";
  	if (Zp == 1) {
  		FileZiegler = "01-H.dat";
  	} else {
  		if (targ_state == "gas") {
			FileZiegler = "02-He-gas.dat";
		} else {
			FileZiegler = "02-He.dat";
		}
  	}
  	FileZiegler = path + FileZiegler;
  	cout<<endl;
  	cout<<"Useds Ziegler file: "<<FileZiegler<<endl;
  	cout<<endl;
  	int z;
  	double temp[12];
  	ifstream in0(FileZiegler.c_str());
  	int flag = 0;
  	while (!in0.eof()) {
  		in0 >> z;
		for (int k=0; k<12; k++) in0 >> temp[k];
		for (int m=0; m<Nelem; m++) {
			if (z == Zt[m]) {
				for (int i=0; i<9; i++) Zgl[m][i] = temp[i];
			}
		}
  	}
  	in0.close();
  	cout<<"Zt    ------------------------------- Ziegler parameters ----------------------------------"<<endl;
  	for (int j=0; j<Nelem; j++) {
  		cout<<left<<setw(6)<<Zt[j];
  		for (int m=0; m<9; m++) cout<<left<<setw(10)<<Zgl[j][m];;
		cout<<endl;
  	}
  	cout << endl;

	string FileOut = "Disang_Sfactor_Dresda.txt";
	ofstream out(FileOut.c_str());
	const int NE = 4;
	double E[NE] ={500, 608, 700, 800};
	string E_s[NE];	
	string zz[NE];
	ostringstream ss[NE], nn[NE];
	for (int j=0; j<NE; j++) {
		nn[j] << j;
		ss[j] << E[j];
		E_s[j] = ss[j].str();
		zz[j] = nn[j].str();
		
	}
	string FileIn = "Yield_Dresda.txt";
	FileIn = DirDat + FileIn;
	double PM = 0;
	for (int j=0; j<Nelem; j++) {
		PM = PM + At[j]*w[j];
	}
	double Avo = 6.022e23;
	double DperCompound = 2;
	double Compoundpercm2 = Dpercm2/DperCompound;
	tau = PM*Compoundpercm2/Avo;          // tau target [g/cm2]
	const int Ndet = 7;

	double teta_MB2     = 36;
	double teta_EB17    = 39;
	double teta_EB18AF  = 65;
	double teta_EB18BGE = 90;
	double teta_MB1     = 118;
	double teta_EB18CD  = 115;
	double teta_Ron100  = 143.5;
	double tetadeg[Ndet] = {teta_MB2, teta_EB17, teta_EB18AF, teta_EB18BGE, teta_MB1, teta_EB18CD, teta_Ron100};
	double teta[Ndet];
	//double Dteta[Ndet];
	double X[Ndet], DX[Ndet];
	double pi = acos(-1);
	for (int j=0; j<Ndet; j++) {
		teta[j] = tetadeg[j]*pi/180;
		X[j] = cos(teta[j]);
		DX[j] = fabs(0.01*X[j]);
	}

	string Detector[Ndet] = {"teta_MB2, teta_EB17", "teta_EB18AF", "teta_EB18BGE", "teta_MB1", "teta_EB18CD", "teta_Ron100"};
	double Tp[NE], Q[NE], C[NE][Ndet], DC[NE][Ndet];
	double y[NE][Ndet], Dy[NE][Ndet];
	string dummy_s;
	double dummy;
	int LineToJump = 1;
	ifstream in(FileIn.c_str());
	for (int j=0; j<LineToJump; j++) in.ignore(256, '\n');
	for (int j=0; j<Ndet; j++) {
		in>>dummy_s>>dummy;
		for (int k=0; k<NE; k++) {
			in>>y[k][j]>>Dy[k][j];	
		}		
	}
	in.close();

	TCanvas *ZZ[NE];
	int N;
	string GrName;
	double Y[Ndet], DY[Ndet];
	TGraphErrors *G[NE];
	TF1 *fit[NE];
	double Ebeam;
	double a[NE][Npar], Da[NE][Npar];
	double A0[NE], dA0[NE];
	double S[NE][Ndet], dS[NE][Ndet]; 
	double Eout[NE], Eeff[NE], Ecmeff[NE];
	double W[NE][Ndet];
	double Sum;

	for (int n=0; n<NE; n++) {
		ZZ[n] = new TCanvas(zz[n].c_str(), zz[n].c_str(), 800, 600);
		ZZ[n]->cd();
		Ebeam = E[n];
		GrName = E_s[n]+"keV";
		for (int k=0; k<Ndet; k++) {
			Y[k] = y[n][k];
			DY[k] = Dy[n][k];
		}
		G[n] = new  TGraphErrors(Ndet, X, Y, DX, DY);
		G[n]->SetTitle(GrName.c_str());
		G[n]->Draw("AP");
		string ParName[Npar];
		char *id = new char[10];
		for (int j=0; j<Npar; j++) {
			sprintf(id, "a%d", j);
			ParName[j] = id;
		}
		fit[n] = new TF1("fit", fitfun(Npar), -1, 1, 3);
		for (int j=0; j<Npar; j++) {
  			fit[n]->SetParName(j,ParName[j].c_str());
		}
		G[n]->Fit(fit[n], "R");
		fit[n]->Draw("same");
		for (int j=0; j<Npar; j++) {
			a[n][j] = fit[n]->GetParameter(j);
			Da[n][j] = fit[n]->GetParError(j);
		}
		A0[n] = a[n][0];
		dA0[n] = Da[n][0];
		Eout[n] = Tfin(Zp, Ap, Nelem, E[n]);
		double dT = (E[n]-Eout[n])/1000;
		Eeff[n] = (E[n]+Eout[n])/2;
		Ecmeff[n] = Eeff[n]*md/(mp+md);
		//out<<endl;
		out<<"Ebeam = "<<E[n]<<"        Eeff = "<<Eeff[n]<<"       Ecmeff = "<<Ecmeff[n]<<endl;
		out<<left<<setw(15)<<"teta"<<left<<setw(15)<<"S [eV*b]"<<left<<setw(15)<<"dS [eV*b]"<<endl;
		double T = E[n];
		Sum = 0;
		while (T>Eout[n]) {
			Sum = Sum + F(T);
			T = T - dT;
		}
		Sum = Sum*dT;
		for (int k=0; k<Ndet; k++) {		
			W[n][k] = disang(X[k], Npar, a[n]);
			//W[n][k] = disang(X[k], a[n][1], a[n][2]);
			S[n][k] = 1e24*1e3*Y[k]/W[n][k]/Sum;
			dS[n][k] = 1e24*1e3*DY[k]/W[n][k]/Sum;
			out<<left<<setw(15)<<tetadeg[k]<<left<<setw(15)\
			<<S[n][k]<<left<<setw(15)<<dS[n][k]<<endl;
		}
		out<<endl;
	}
	out.close();
	cout<<"Results on file "<<FileOut<<endl;
//	app.Run(true);		
}