#include  "source.h"

using namespace std;

TF1* f_T_particle;
TF1* f_LorentzFactor_e ;
TF1* f_Beta_e;
TF1* f_STP_e;
TF1* f_STP_e_over_rho;

double MEP_medium = MEP_Si; // MEP_Calculator(Z_Si); //mean potential energy
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 RME_e = 0.51099895069;//RME_calculator(m_electron); // rest mass energy


void fittingcode3() {
    

    TCanvas *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
    
    string dataname = "Si_estar.txt";
    TGraph *gr1 = new TGraph(Form("%s", dataname.c_str()) , "%lg%*lg%*lg%lg%*lg%*lg%*lg"); //total STP (MeV. cm^2 /g)
 
    
    if(gr1){
        gr1->SetName("gr1");
        gr1->SetTitle("estar, pstar, astar");
        gr1->SetMarkerSize(0.6);
        gr1->SetMarkerStyle(20);
        gr1->SetMarkerColor(kRed);
        //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");

    }
    else {cout << "check gr1!" << endl; }
    

    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, 100.);
    // Print initial parameter values for debugging
    f_STP_e_over_rho->Print();
    
    
    f_STP_e_over_rho->SetNpx(5000);
     
     gr1->Fit(f_STP_e_over_rho, "MR");
     gr1->GetXaxis()->SetRangeUser(1.e-2, 1.);
     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);
     
    f_STP_e_over_rho->SetLineColor(kBlack);

    gr1->Draw("AP");

    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 result2= sqrt( (pow(x[0]/RME_e+1. ,2.) -1. ) /(pow(x[0]/RME_e +1. ,2.)) ) ;
    double T=x[0];
    double gamma=(T/RME_e) +1.;
    double result2 = sqrt(1 -1/gamma/gamma);
    return result2;
}

double STP_e(double *x,double *par){ //Unit is MeV/m

    
    double gamma=(x[0]/RME_e) +1.;
    double beta = sqrt(1.-1./gamma/gamma);
      
    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);
    return par[0] * arg  ; //MeV * cm^2 /g
   
}