#include <fstream>
#include <iostream>
#include <string>
#include <math.h>
#include "TFile.h"
#include "TMath.h"
#include "TTree.h"
#include "TRandom.h"
#include <stdlib.h>
#include <sstream>
#include <stdio.h>
#include <stdlib.h>
#include <cmath>


using namespace std;

void histogram1()

{
   //gROOT->Reset();
    double_t MACS2_sum=0.;
    double_t MACS1_sum=0.;
    double_t MACS_sum=0.;
    auto c1 = new TCanvas ("c1","c1",0,0,700,250);
        c1->SetGrid();
        c1->SetLogy();
        c1->SetLogx();
   

        auto mg  = new TMultiGraph();

        auto Xaxis = mg->GetXaxis();
        Xaxis->SetTickLength(0.02);
        Xaxis->SetTitle("Neutron Energy [eV]");
        Xaxis->SetLimits(0.0001,10000000000);
  
     
    TGraph *gr1 = new TGraph("test.txt");

     
        mg->Add(gr1);

        mg->Draw("AL");
    
    
   // TSpline3 *s = new TSpline3("Test",gr1,"b1e1",0,0);
    TSpline5 *s = new TSpline5("Test",gr1);

    for(int k=10000;k<=100000;k=k+5000){
        
    
    for (double_t i=0; i<100000; i=i+1)
        {
            double E=i*209./210.;
            double CS=gr1->Eval(E);
        
            
            //double T1=CS;
            // double T2=E*exp(-E/50000)*10;
           // double T_total= (E)*CS*exp(-E/50000)*10;
            double_t MACS=(E)*CS*exp(-E/k)*1;
            MACS_sum=MACS+MACS_sum;
           // cout<<MACS<<endl;
            
            double MACS1=(E/2.1)*CS*exp(-E/2.1/k)*1;
           double MACS2=(E/2.1)*exp(-E/2.1/k)*1;
            MACS1_sum=MACS1+MACS1_sum;
            MACS2_sum=MACS2+MACS2_sum;
        

        
       // cout<<"E= "<<E<<"  cs= "<<CS<<" MACS=  "<<MACS<<"  MACS_sum=   "<<MACS_sum<<endl;
       // cout<<" MACS=  "<<MACS<<endl;
    }
   
        MACS_sum=(2/sqrt(3.14))*MACS_sum/k/k;
    //double_t CF= MACS_sum/0.003069 ;
    // double_t    MACS_EXP=CF*.00183;
       double_t  CF= .00183/0.003069;
        double_t MACS_EXP = MACS_sum*CF;
        std::cout<<"k= "<<k/1000<<"   MACS_sum=  "<<MACS_sum*1000<<"    CF= "<<CF<<"  MACS_EXP=  "<<MACS_EXP*1000<<endl;

    
   // std::cout<<" MACS1_sum=  "<<MACS1_sum<<endl;
   // std::cout<<" MACS2_sum=  "<<MACS2_sum<<endl;
        
       

    
       double  MACS_sum_INT=(2/sqrt(3.14))*MACS1_sum/MACS2_sum;
    double_t CF1= MACS_sum_INT/0.003069 ;
     double_t    MACS_EXP_INT=CF1*.00183;
        
   //std::cout<<"k= "<<k/1000<<"   MACS_sum_INT=  "<<MACS_sum_INT*1000<<"    CF1=" <<CF1<<"  MACS_EXP=         "<<MACS_EXP_INT*1000<<endl;
// std::cout<<"  "<<endl;
        // std::cout<<"  "<<endl;
    
    }


 
    s->SetNpx(5 * gr1->GetN() ); 
    
    s->SetLineColor(kRed);
    
    s->Draw("same");


    // add ROOT interpolator

    std::vector<double> xp(gr1->GetX(), gr1->GetX() + gr1->GetN() ); 
    std::vector<double> yp(gr1->GetY(), gr1->GetY() + gr1->GetN() ); 
    ROOT::Math::Interpolator itp(xp.size(), ROOT::Math::Interpolation::kCSPLINE);
    itp.SetData(xp,yp);

    // compute points to interpolate
    std::vector<double> x2(xp.size()-1);
    std::vector<double> y2(xp.size()-1);
    for (size_t i = 0; i < xp.size()-1 ; ++i) {
       x2[i] = xp[i] + 0.5*(xp[i+1]-xp[i]);
       y2[i] = itp.Eval(x2[i] );
    }


    TGraph * gr2 = new TGraph(x2.size(), x2.data(), y2.data() );
    gr2->SetLineColor(kGreen);
    gr2->Draw("L SAME"); 
}