#include "TCanvas.h"
#include "TBrowser.h"
#include "TH2F.h"
#include "TH1F.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TMath.h"
#include <fstream>
#include "TFrame.h"
#include "TSystem.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TFile.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TBox.h"
#include "TRandom.h"
#include "TObject.h"
#include "TObjString.h"
#include <iostream>
#include <cstdio>
#include <string>
#include <sstream>
#include "TLine.h"
#include "TTree.h"
#include "TBrowser.h"
#include "TF1.h"
#include <TStyle.h>
#include <THStack.h>
#include <TPad.h>
#include "TRandom.h"


using namespace std;

TH1F *toftoene_fixedbin(TH1F *htof,Double_t L,Double_t offset,Int_t binoffset);
void WriteHistogramFile(TString filename,TH1F *h);

float Sn= 0  ; // neutron separation energy in MeV
float Amass= 0; // mass of your nucleus in atomic mass units

/// L flight path in cm, assumed
/// offset ... any offset in time of flight from daq (N/A here, hence 0)
/// binoffset this is just to exclude very low tof bins which are at too high energy to be of interest 10000 should be ok

void run(){

    /// LOAD TIME OF FLIGHT SPECTRUM
    
    TH1F *hene=toftoene_fixedbin(NAMEOFTOFSPECTRUM,18500,0,10000);
    

    hene->Draw();  /// draw tof histogram concerted to energy
    
    TH1F *h_me=(TH1F*)hene->Clone();
    
    //// NOW moving on to division by flux
     /// the neutron flux is given per bunch assuming a standard proton intensity of 7E12 protons / bunch. Hence the division is by no of protons and then multiplied by 7E12 to get to the number of standard bunches
    
    
   h_me->Scale(1/(h_mon->GetBinContent(4)/7E12));

    
    
    /// file that contains the neutron flux. We are looking at the flux in isolethargic units
    TFile *fflux = new TFile("nTOF-Ph2_fluence_2009-2011_6Dec2011.root", "read"); //for EMMA: evalFlux_prelim.root
    
    TH1F *sim = (TH1F*)fflux->Get("hNFluenceEVALUATED2011"); // FOR EMMA: hEval_Abs
    
    //Neutron flux that is independent of the binning of the histogram. The number of neutrons between E1 and E1 can be found by multiplying the value in the histogram by the logarithmic bin width, i.e. log(E2)-log(E1)
      
        TH1F *yield1 = (TH1F*)h_me->Clone();
        
        cout<<"DIVIDE BY FLUX"<<endl;
        for(Int_t i=1; i<=h_me->GetNbinsX(); i++)
        {
            
            float auxbin1 = sim->FindBin(h_me->GetBinCenter(i));  // find the bin in the flux histogram
            float scaler1 = TMath::Log(h_me->GetBinLowEdge(i+1)/h_me->GetBinLowEdge(i)); // determine the binwidth of count histogram

            yield1->SetBinContent(i,h_me->GetBinContent(i)/(sim->GetBinContent(auxbin1)*scaler1));
            // this is the division by the flux. The number of coutns gets divided by the number of neutrons in that particular bin
            
            if(sim->GetBinContent(auxbin1)<=0)yield1->SetBinContent(i,0); // this is in case the flux histgram has a content <=0, so we do not divide by a negative or 0
           
             if(sim->GetBinContent(auxbin1)>0){
                 
            float binerr_rel=sqrt(pow(h_me->GetBinError(i)/h_me->GetBinContent(i),2) + pow(sim->GetBinError(auxbin1)/sim->GetBinContent(auxbin1),2));
                 
                 // this calculates the relative uncertainty taking into consideration the sample uncertainties and the uncertaities in the flux histogram
                 
            float binerr=binerr_rel*yield1->GetBinContent(i);
                  // calculation of absolute uncertainty
                 
            if(h_me->GetBinContent(i)==0)binerr=h_me->GetBinError(i)/(sim->GetBinContent(auxbin1)*scaler1);
                 /// this is again to avoid divition by  0

            yield1->SetBinError(i,binerr);
            }



        }


        for(int g=1;g<=h_me->GetNbinsX();g++)
        {
    yield1->SetBinContent(g,yield1->GetBinContent(g)*threshscaler/(Sn+  Amass/(1.0086649+Amass)*yield1->GetBinCenter(g)/1e6));
    yield1->SetBinError(g,yield1->GetBinError(g)*threshscaler/(Sn+  Amass/(1.0086649+Amass)*yield1->GetBinCenter(g)/1e6));
    }
     
    
   // this is an additional correction since the neutron beam size varies slightly with neutrom energy. Emma, you can skip that for now as it will be different for you, and we dont have that information yet.
    
        cout << "BIF CORRECTION" <<endl;
        
        TFile *bif=new TFile("BIF_2011_norm.root","read");
        TH1F *hbif=(TH1F*)bif->Get("histo");
        
        
        
        for(int f=1;f<=yield1->GetNbinsX();f++){
            float center=yield1->GetBinCenter(f);
            int binni=hbif->FindBin(center);
            
            yield1->SetBinContent(f,yield1->GetBinContent(f)/hbif->GetBinContent(binni));

            float binerr_rel=sqrt(pow(yield1->GetBinError(f)/yield1->GetBinContent(f),2) + pow(hbif->GetBinError(binni)/hbif->GetBinContent(binni),2));
            float binerr=binerr_rel*yield1->GetBinContent(f);
            if(yield1->GetBinContent(f)==0)binerr=yield1->GetBinError(f)/hbif->GetBinContent(binni);

            yield1->SetBinError(f,binerr);
        }

    cout <<"BIF CORRECTION DONE"<<endl;
    
    
    
        
    yield1->GetXaxis()->SetTitle("energy (eV)");
    yield1->GetYaxis()->SetTitle("yield");

    cout <<"WRITE FILE"<<endl;
    
    // we are writing the yield as root file and as a data file, since that is the input for resonance fitting.

    TFile *gefinalyield= new TFile(NAME, "update");
    yield1->Write(histosam);
    gefinalyield->Write();
    gefinalyield->Close();

    cout <<"ROOT FILE WRITTEN"<<endl;

    WriteHistogramFile(FILENAME,yield1);

    cout <<"TEXT FILE WRITTEN"<<endl;
    
    
}





TH1F *toftoene_fixedbin(TH1F *htof,Double_t L,Double_t offset,Int_t binoffset){
    
    
int number=htof->GetNbinsX();
float lowedge;
float eneedge[number-binoffset];

const double Mass=939565560.81;            //Neutron Mass
const double c0=29.972458;                 // speed of light

cout<<number<<endl;

for(int i=binoffset;i<=number;i++){
    lowedge=htof->GetBinLowEdge(i)+L/c0+offset;
    eneedge[number-i]=(1/pow(1-((L/(lowedge)/c0)*(L/(lowedge)/c0)),0.5)-1)*Mass;
  
    if(i==30000)cout<<lowedge<<" "<<eneedge[number-i]<<" "<<number-i<<endl;
}


TH1F *hene =new TH1F("Counts","Counts",number-binoffset-1,eneedge);
hene->GetXaxis()->Set(number-1-binoffset,eneedge);


for(int i=1;i<=hene->GetNbinsX();i++){
  
    hene->SetBinContent(i,htof->GetBinContent(number-i));
    hene->SetBinError(i,htof->GetBinError(number-i));

}
    
    
    return hene;

}




void WriteHistogramFile(TString filename,TH1F *h){

  
  FILE * fout;
  
  fout = fopen(filename,"w");
  
  for(Int_t i=1;i<=h->GetNbinsX();i++)
    {
double bincenternew=h->GetBinCenter(i);
double binContent=h->GetBinContent(i);
//if(binContent<0)binContent=0;
double binError=h->GetBinError(i);
if(binError<=0)binError=1e-8;
  if(bincenternew>0.02571)fprintf(fout,"%20.10f%20.10f%20.10f\n",bincenternew,binContent,binError);
      
    }
  fclose(fout);
  
  cout << "*************************************"<< endl;
  cout << " File " << filename  <<" written, with format: E_center_bin, Content, error "<<endl;
  cout << "*************************************"<< endl;


}