//CHOOSE  1keV binning for consistency

#include <stdio.h>
#include<locale.h>
#include <iostream>
#include "TCanvas.h"
#include "TFitResult.h"
#include "TDirectory.h"
#include "TAxis.h"
#include "TH1D.h"
#include "TStyle.h"  //for gStyle decleration
#include <fstream>
#include "TF1.h"
#include "TFormula.h"
#include "TMath.h"
#include "Math/DistFuncMathCore.h"
#include "TROOT.h"
#include "Math/Minimizer.h"
#include "Minuit2/MnUserParameterState.h"
#include "TMatrixD.h"

#include <TROOT.h>
#include <TFile.h>
#include <TGraphErrors.h>
#include <TF1.h>

using namespace std;


const double PI = 3.141592653589793;
const double Avogadro   =  6.022E23; //atoms/mole
const double c_square   =  931.5 ;  //MeV/u
const double FourPiRsquare= 4.0* PI * (2.818*1.E-15) *  (2.818*1.E-15); //m^2
const double m_electron =  0.00054858 ; //in akb, u;  9.11E-31;  //kg

double A_Si       = 28.0  ;
double rho_Si     = 2.33*1.E3 ;
double N_Si       = (rho_Si*Avogadro)/ (A_Si*1.0E-3);
double Z_Si       = 14.0 ;
double MEP_Si     = 172.*1.0E-6 ;

TF1* f_T_particle;
TF1* f_LorentzFactor_e ;
TF1* f_Beta_e;
TF1* f_STP_e;
TF1* f_STP_e_over_rho;
double RME_calculator(double m_particle);
double T_particle(double *x,double *par);
double LorentzFactor_e(double *x,double *par);
double Beta_e(double *x, double *par);
double STP_e(double *x, double *par);
double STP_e_over_rho(double *x, double *par);

double MEP_medium = MEP_Si; // MEP_Calculator(Z_Si); //mean potential energy
double RME_e = RME_calculator(m_electron); // rest mass energy


void STP_Fit() {
    
    
    
    auto *c1 = new TCanvas("c1", "c1",294,78,1343,770);
    c1->cd();
    //c1->SetLogy();
    //c1->SetLogx();
    ////////////  Inner Frame
    //gPad->SetLogy(0); TO GET BACK TO NORMAL MODE
    
    //Fit works between 0.01 to 100 MeV
    string dataname ;//example = "Al_estar.txt";
    cout<<"Please put the name of datafile as Si_estar.txt"<<endl;
    cin >> dataname;
    auto *gr1 = new TGraph(Form("/Users/spyhunter0066/Library/CloudStorage/OneDrive-harran.edu.tr/MacbookCplusplus/RangeCalculators/Fitting/%s", dataname.c_str()) , "%lg%*lg%*lg%lg%*lg%*lg%*lg"); //total STP (MeV. cm^2 /g)

    
    //======================== For Electron @Si =============================================
    //*************************************************************************************
    //  TF1::TF1(const char* name, void* fcn, Double_t xmin, Double_t xmax, Int_t npar); !!!
    //  We can get the address of a function by just writing the function’s name without parentheses. !!!
 
     gr1->GetXaxis()->SetRangeUser(1.e-2, 20.);
     gr1->GetYaxis()->SetRangeUser(1.e-2, 50.);
     gr1->GetXaxis()->SetTitle("T_particle(MeV)");
     gr1->GetYaxis()->SetTitle(Form("dE/dx(MeV.cm^{%d}/g )",2));
     gr1->GetXaxis()->SetLabelSize(0.03);
     gr1->GetXaxis()->SetLabelOffset(0.001);
     gr1->GetYaxis()->SetLabelSize(0.03);
     gr1->GetYaxis()->SetLabelOffset(0.001);
     gr1->GetXaxis()->SetTitleOffset(1.5);
     gr1->GetXaxis()->SetTitleSize(0.025);
     gr1->Draw("AP");

    
    
    double min_Limit= 1.e-2;
    double max_Limit= 20. ;
    f_T_particle = new TF1("f_T_particle", T_particle, min_Limit,max_Limit,0);
    f_LorentzFactor_e = new TF1("f_LorentzFactor_e", LorentzFactor_e, min_Limit,max_Limit,0 ); //between 1 and 10
    f_Beta_e = new TF1("f_Beta_e",Beta_e, min_Limit,max_Limit,0); //between 0 and 1
    f_STP_e = new TF1("f_STP_e", STP_e, min_Limit, max_Limit,0);
    f_STP_e_over_rho= new TF1("f_STP_e_over_rho", STP_e_over_rho, min_Limit, max_Limit,1);
    f_STP_e_over_rho->SetParameter(0, 0.);

    gr1->Fit(f_STP_e_over_rho, "MR");

    // f_STP_e_over_rho->SetLineColor(kBlack);
    // f_STP_e_over_rho->Draw("same");
    
    // Get the number of parameters of the fitted function
    int nParams = f_STP_e_over_rho->GetNpar();
    double ParStore[nParams];
    // Print the parameters
    for (int i = 0; i < nParams; ++i) {
        ParStore[i] = f_STP_e_over_rho->GetParameter(i);
       cout << "Parameter " << i << ": " << ParStore[i]  << endl;
    }
    
}

double T_particle(double *x,double  *par){
    return x[0];
}

double RME_calculator(double m_particle){
    double RME;
    RME= m_particle * c_square; //u* (MeV/u) = MeV
    return RME; //MeV
    
}

double LorentzFactor_e(double *x,double  *par){
    //cout <<"LorentzFactor_e= " <<((x[0] / RME_e) +1.0) << endl;
    //cout <<"RME_electron= " << RME_e << endl;
    double result1= ((x[0] / RME_e) +1.0);
    return result1;
}

double Beta_e(double *x, double *par){
    //  cout << "Beta_e= " << ... << endl;
    double T = x[0];
    double gamma=(T/RME_e) +1.;
    double beta = sqrt( 1.- (1./ pow(gamma, 2.)) );
    return  beta;
}

double STP_e(double *x,double *par){ //Unit is MeV/m
//cout << "f_Beta_e= " <<f_Beta_e->Eval(x[0]) << "  f_LorentzFactor_e= " << f_LorentzFactor_e->Eval(x[0]) << endl;
    double T = x[0];
    double gamma=(T/RME_e) +1.;
    double beta =   sqrt( 1.- (1./ pow(gamma, 2.)) );
    double result3 = FourPiRsquare* RME_e* N_Si * Z_Si * (1. /pow(beta,2)) * ( log( beta*gamma*sqrt(gamma-1.)*RME_e/MEP_medium) + ((1./(2.*pow(gamma,2))) * ( pow(gamma-1,2.)/8. +1. -log(2.)*(pow(gamma,2) + (2.*gamma) -1.))) );
    
    return result3  ;
}

double STP_e_over_rho(double *x,double *par){ // (MeV .m^2/kg) * 10 = (MeV.cm^2 /g)
    double arg = 10. * (f_STP_e->Eval(x[0]) / rho_Si);
    printf("par[0] =%f\n",par[0]);
    return  par[0] * arg  ; //MeV * cm^2 /g
}