//C++ libraries
#include <iostream>
#include <locale>
#include <string>
#include <vector>
#include <iomanip>
//ROOT libraries
#include "TFormula.h"
#include "TMath.h"
#include "TROOT.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TLatex.h"
#include "TAxis.h"
#include "TPaveText.h"
#include "TText.h"

#define PI 3.14159265

//using namespace std;

int roottest(){
     
    std::vector<float>T_energyLIST;
    float T_energy_max = 20.;
    float T_energy_min = 0.;
    float increment = 0.1;
    float total_number = (T_energy_max- T_energy_min) / increment ;
    
    for (unsigned int i = 0; i <total_number ; i++ ){
        T_energyLIST.push_back(T_energy_max- (i*increment));
    }
    
     std::cout << "T_energyLIST.size():  " << T_energyLIST.size() << std::endl;
     
     /*
     for (unsigned int i = 0; i < T_energyLIST.size(); i++ ){
     std::cout << "T_energy values: " << T_energyLIST[i] << std::endl;
     }
     */
    float z_particle = 2.0 ;
    float M_particle = 4.002603 ; // in unit of u, amu
    float rho_material= 2330. ; // kg/m^3
    float A_material = 0.028 ; // kg/ mole
    float Z_material = 14.0 ;
    float I_material = 0.000172 ; //  MeV
    float area= TMath::Power(10.,-28.); //m^2, 4.pi.r^2
    float mc2= 0.511; //MeV
    float NA= 6.022*TMath::Power(10.,23.); //atoms / mole
    float c2= 931.5 ; // MeV/u
    float N_material =  (rho_material * NA )/ A_material ;
    //Arrays
    std::vector<float> gamma;
    std::vector<float> Beta;
    std::vector<float> V_particle;
    std::vector<float> AverageEnergyLossPerUnitDistance_O;
    std::vector<float> AverageEnergyLossPerUnitDistance_I;
    std::vector<float> result;
    
    for (unsigned int i = 0; i < T_energyLIST.size(); i++ ){
        gamma.push_back(T_energyLIST[i] / (M_particle * c2) + 1. ) ;
        Beta.push_back( TMath::Sqrt(1.- (1./(gamma[i]* gamma[i]) ) ) ) ;
        V_particle.push_back( Beta[i] * 3. * TMath::Power(10.,8) ) ;
        AverageEnergyLossPerUnitDistance_O.push_back( area * TMath::Power(z_particle, 2.)* mc2 * N_material * Z_material * (1./ TMath::Power(Beta[i], 2.)) );
        AverageEnergyLossPerUnitDistance_I.push_back( TMath::Log(2.* mc2* TMath::Power(Beta[i], 2.)* TMath::Power(gamma[i], 2.)* (1./ I_material) )- TMath::Power(Beta[i], 2.) );
        result.push_back(AverageEnergyLossPerUnitDistance_O[i] * AverageEnergyLossPerUnitDistance_I[i] );
    }
    
        
        std::cout << std::setw(20) <<  std::setprecision(6)  << "Gamma" << "\t" <<
        std::setw(20) <<  std::setprecision(6)  << "Beta" <<"\t" <<
        std::setw(20) <<  std::setprecision(6)  << "V_particle (m/s)" << "\t" <<
        std::setw(20) <<  std::setprecision(6)  << "dE/dx 1. part"  << "\t" <<
        std::setw(20) <<  std::setprecision(6)  << "dE/dx 2. part" << "\t" <<
        std::setw(20) <<  std::setprecision(6)  << "dE/dx (MeV/m)" << "\t" <<
        std::setw(20) <<  std::setprecision(6)  << "(dE/dx)*(1/rho), (MeV*m^{2}/kg)" << std::endl;
    
    
    for (unsigned int i = 0; i < T_energyLIST.size(); i++ ){
    std::cout << std::setw(20) <<  std::setprecision(6)  << gamma[i] << "\t" <<
                 std::setw(20) <<  std::setprecision(6)  << Beta[i] <<"\t"   <<
                 std::setw(20) <<  std::setprecision(6)  << V_particle[i] << "\t"  <<
                 std::setw(20) <<  std::setprecision(6)  << AverageEnergyLossPerUnitDistance_O[i]  << "\t" <<
                 std::setw(20) <<  std::setprecision(6)  << AverageEnergyLossPerUnitDistance_I[i] << "\t"  <<
                 std::setw(20) <<  std::setprecision(6)  <<  result[i] << "\t"  <<
                 std::setw(20) <<  std::setprecision(6)  <<  result[i]*(1/rho_material) << std::endl;

       }
    
    // create graph
      int n= T_energyLIST.size();
      double x[n], y[n];
      for (Int_t i=0;i<n;i++) {
         x[i] = T_energyLIST[i];
         y[i] = result[i];
      }

    TCanvas *c1 = new TCanvas("c1", "Graph Draw Options",200,185,600,400);
    c1->Range(-5.459234,-3905860,25.05125,1148714);
    c1->SetFillColor(0);
    c1->SetBorderMode(0);
    c1->SetBorderSize(2);
    c1->SetLeftMargin(0.1789298);
    c1->SetFrameBorderMode(0);
    c1->SetFrameBorderMode(0);
    
    
    TGraph *gr  = new TGraph(n,x,y);
    gr->GetXaxis()->SetTitle("T_particle");
    gr->GetXaxis()->CenterTitle(true);
    gr->GetXaxis()->SetLabelFont(42);
    gr->GetXaxis()->SetTitleOffset(1);
    gr->GetXaxis()->SetTitleFont(42);
    gr->GetYaxis()->SetTitle("dE/dx");
    gr->GetYaxis()->CenterTitle(true);
    gr->GetYaxis()->SetLabelFont(42);
    gr->GetYaxis()->SetTitleFont(42);
    gr->SetLineWidth(1);
    gr->SetLineColor(kRed);
    gr->SetMarkerStyle(20);
    gr->SetMarkerSize(0.5);
    gr->SetMarkerColor(kBlack);
    gr->Draw("ACP");
    
    TPaveText *pt = new TPaveText(0.3453512,0.9342,0.6546488,0.995,"blNDC");
    pt->SetName("title");
    pt->SetBorderSize(0);
    pt->SetFillColor(0);
    pt->SetFillStyle(0);
    pt->SetTextFont(42);
    TText *pt_LaTex = pt->AddText("dE/dx vs T_particle");
    pt->Draw();
    c1->Modified();
    c1->cd();
    c1->SetSelected(c1);
    
    /*
    // superimpose the second graph by leaving out the axis option "A"
      gr2->SetLineWidth(3);
      gr2->SetMarkerStyle(21);
      gr2->SetLineColor(2);
      gr2->Draw("CP");
    */
    std::cout<<"done"<<std::endl;
    
    return 0;
    
}