#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <sstream>
#include <cmath>
#include "TRandom.h"
#include <math.h>
#include <TPolyLine.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <ctype.h>
#include "TTree.h"
#include <TROOT.h>
#include <TApplication.h>
#include <TRint.h>
#include <TSystem.h>
#include <TH1.h>
#include <TH2.h>
#include <TAxis.h>
#include <TGaxis.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TGraphAsymmErrors.h>
#include <TMultiGraph.h>
#include <TStyle.h>
#include <TKey.h>
#include <TLegend.h>
#include <TColor.h>
#include <TPad.h>
#include <TText.h>
#include <TPaveText.h>
#include <TBox.h>
#include <TLine.h>
#include <TMarker.h>
#include <TLatex.h>
#include <TMath.h>
#include <TF1.h>
#include <TFile.h>
#include <TClass.h>

#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "Math/Functor.h"

#include <TVirtualFitter.h>

using namespace std;

void RebinProperly(TH1F* hin, TH1F* hout);


Char_t *inputfile="6Li_endf.root";
Char_t *histo="xshighbin";
Char_t *outputfile="countsLi.root";
Char_t *output="counts";
int k=1;
 
void run()
{


  // Read the flux

    TFile *fflux = new TFile("evaluated_flux_EAR2_DEC2022.root", "read");
    TH1F *eval = (TH1F*)fflux->Get("h_flux_ear2");
    
    
//Read the cross section

    TFile *fcross=new TFile("10B_endf.root", "read");;
TH1F *hcross =(TH1F*)fcross->Get("xshighbin");
    
    
    // create counr histogram with logarithmic binning for x axis, commonly used for neutron energy histograms
    
          float Ene[20001];
          for(Int_t u=0;u<=20000;u++)
          {
              Ene[u]=pow(10,(float(u)-6000)/2000);
          }
          
          TH1F *hcounts =new TH1F("","histo",20000,Ene);
          hcounts->GetXaxis()->Set(20000,Ene);

// function to rebin the cross section histogram to the same binning as the counts histogram
    
    RebinProperly(hcross,hcounts);
    
    

    
    
  Int_t auxbin1;
  Float_t scaler1;
    
    // multiply by the n_TOF neutron flux (binning is in units of ExdPhi/dE, so independent of binning)

  for(Int_t i=1; i<=hcounts->GetNbinsX(); i++)
    {
      auxbin1 = eval->FindBin(hcounts->GetBinCenter(i));
      scaler1 = TMath::Log(hcounts->GetBinLowEdge(i+1)/hcounts->GetBinLowEdge(i));
      if(hcounts->GetBinContent(i)>0 && eval->GetBinContent(auxbin1)>0)hcounts->SetBinContent(i, hcounts->GetBinContent(i)*eval->GetBinContent(auxbin1)*scaler1);

    }
    
    
    
    
    for(int i=1;i<=hcounts->GetNbinsX();i++){
        if(hcounts->GetBinContent(i)==0)hcounts->SetBinContent(i,1);
        hcounts->SetBinError(i,0);
    }
 
    hcounts->Scale(1/7000.);
    
    

    
    hcounts->SetTitle("");
    hcounts->GetXaxis()->SetTitle("Neutron Energy (eV)");
    hcounts->GetYaxis()->SetTitle("Counts (arbitrary)");
    
TFile *fnew=new TFile(outputfile,"recreate");
    
    hcounts->Write(output);
    fnew->Close();
    
    
    
    fflux->Close();
    fcross->Close();

}







void RebinProperly(TH1F* hin, TH1F* hout){
    
    for(int i=1;i<=hout->GetNbinsX();i++)
    {
        float content=0;
        float error=0;
        float errorsquare=0;
        int zahler=0;
        float specedgelow=hout->GetBinLowEdge(i);
        float specedgehigh=hout->GetBinLowEdge(i+1);
        
        int binlow=hin->FindBin(specedgelow);
        int binup=hin->FindBin(specedgehigh);
        
        if(binlow==binup){content=hin->GetBinContent(binup);error=hin->GetBinError(binup);} // (stat) error probably under-estimated in this case
        if(binlow!=binup){
            
            content=content+hin->GetBinContent(binlow)*(hin->GetBinLowEdge(binlow+1)-specedgelow); //add first bin
            content=content+hin->GetBinContent(binup)*(-hin->GetBinLowEdge(binup)+specedgehigh); //add last bin
            //cout<<hin->GetBinLowEdge(binlow+1)-specedgelow<<" "<<-hin->GetBinLowEdge(binup)+specedgehigh<<endl;
            for(int q=binlow+1;q<binup;q++){content=content+hin->GetBinContent(q)*hin->GetBinWidth(q);} //add intermediate bins
            content=content/(specedgehigh-specedgelow);
            
            
            // error calculation
            errorsquare=errorsquare+hin->GetBinError(binlow)*(hin->GetBinLowEdge(binlow+1)-specedgelow)*hin->GetBinError(binlow)*(hin->GetBinLowEdge(binlow+1)-specedgelow);
            errorsquare=errorsquare+hin->GetBinError(binup)*(-hin->GetBinLowEdge(binup)+specedgehigh)*hin->GetBinError(binup)*(-hin->GetBinLowEdge(binup)+specedgehigh); //add last bin
            for(int q=binlow+1;q<binup;q++){errorsquare=errorsquare+hin->GetBinError(q)*hin->GetBinWidth(q)*hin->GetBinError(q)*hin->GetBinWidth(q);} //add intermediate bins
            error=sqrt(errorsquare)/(specedgehigh-specedgelow);
            
        }
        hout->SetBinContent(i,content);
        hout->SetBinError(i,error);
    }
    
    
}