#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"




void plot(){
    // replace this with name of your singles file
    TFile *f=new TFile("Single_Si29_ded.root","read");
    TString histo, norm, protons;
    
    ///file with a list of the histogram names for one detector
    // format Zn_T_u3_C6D6_runxxxxxx Zn_A_u3_SILI_runxxxxxx Zn_h_info_C6D6_runxxxxxx
    ifstream inni("Si29_C6D64_HistList.dat");
    TCanvas *c=new TCanvas();

    /// replace this with the tof around a resonance, i.e. binlow is the start of the resonance, binhigh is the end of the resonance
    // you can choose any resonance, best is a large one which is usally at big tof values
    /// this will be different for Au and for Zn.

float binlow=2250000;
float binhigh=2400000;

    TLegend *legend = new TLegend(0.65, 0.38, 0.89, 0.89);
    legend->SetBorderSize(1);
    legend->SetFillColor(0);
    legend->SetTextSize(0.03);
    int count2=0;
    int count_protons=0;
    
    
    //// replace this with the number of lines in your file "list"
int lines=59;

            float value[lines];
            float value2[lines];
            float valuesili[lines];
            float valuepkup[lines];
    float staterror[lines];
    while(1){
        
        inni>>histo>>norm>>protons;
        if ( ! inni ) break;

        
        char *s;
        if(count2==0) s="hist";
        if(count2>0)s="hist same";
        count2++;
        cout<<histo<<" "<<norm<<endl;
        
        
        TH1F *h=(TH1F*)f->Get(histo);
        //TH1F *hnorm=(TH1F*)f->Get(norm);
        TH1F *hpro=(TH1F*)f->Get(protons);
        
        
        int bin1=h->FindBin(binlow); //Ge70Bins
        int bin2=h->FindBin(binhigh);
        
        // this gives the total number of count in the resonance
        
        int counter2=h->Integral(bin1,bin2);

      
    
    value[count2-1]=float(counter2)/hpro->GetBinContent(4)*7E12;
        staterror[count2-1]=sqrt(float(counter2))/hpro->GetBinContent(4)*7E12;
             value2[count2-1]=float(counter2);
             valuesili[count2-1]=float(hpro->GetBinContent(4)/hpro->GetBinContent(8));
            valuepkup[count2-1]=float(hpro->GetBinContent(4)/hpro->GetBinContent(7));
   
        h->Scale(1/hpro->GetBinContent(4)*7E12);
        h->Rebin(50);
        h->SetLineColor(count2);
        h->DrawCopy(s);
        legend->AddEntry(h,histo, "l");
    }
    float sum=0;
    float weight=0;
     float weightsum=0;   
    for(int i=0;i<count2;i++){

        weight=value2[i]/(value[i]*value[i]);

       sum=sum+value[i]*weight;


        weightsum=weightsum+weight;

        cout<<value[i]<<" "<<weight<<" "<<i<<endl;
    }


cout<< sum/weightsum <<endl;
   
    float deviation=0;
double stddev=0;

    for(int i=0;i<count2;i++){
        deviation=(sum/weightsum-value[i])/sum*weightsum*100;

stddev=(sum/weightsum-value[i])*(sum/weightsum-value[i])+stddev;

        cout<<" "<<deviation<<" % "<<sqrt(value2[i])/value2[i]*100<<" % counting "<<deviation/(sqrt(value2[i])/value2[i]*100)<< "   "<<i<<endl;
    }
    deviation=sqrt(1/count2*deviation);
    cout<<count2<<" "<<sum/count2<<endl;

double stddev2=sqrt(1/(float(count2-1))*stddev);

//STOP PROTON SCALING

cout<<"std dev is "<<stddev2/(sum/weightsum)*100<<" %"<<endl;
    legend->Draw();
cout<<"                                "<<endl;
    for(int i=0;i<count2;i++){

        cout<<" "<<valuesili[i]<<"  Prot/SILI    "<<valuepkup[i]<<"  Prot/PKUP   "<<valuesili[i]/valuepkup[i]<<" SILI/PKUP"<<endl;
    }


TH1F *hprosi=new TH1F("","",50,0,50);
TH1F *hprook=new TH1F("","",50,0,50);
TH1F *hsipk=new TH1F("","",50,0,50);
    TH1F *hcounts=new TH1F("","",50,0,50);

for(int i=0;i<count2;i++){
hprosi->SetBinContent(i+1,valuesili[i]);
hprook->SetBinContent(i+1,valuepkup[i]);
hsipk->SetBinContent(i+1,valuesili[i]/valuepkup[i]);
    hcounts->SetBinContent(i+1,value[i]);
    hcounts->SetBinError(i+1,staterror[i]);
    }

TCanvas *c3=new TCanvas;
c3->Divide(2,2);
c3->cd(1);
hprosi->Draw();
c3->cd(2);
hprook->Draw();
c3->cd(3);
hsipk->Draw();
    c3->cd(4);
    hcounts->Draw();
    
    
    TCanvas *c5=new TCanvas;
    hcounts->Draw();
}










void plotau(){
    // replace this with name of your singles file
    TFile *f=new TFile("Single_Au_ded.root","read");
    TString histo, norm, protons;
    
    ///file with a list of the histogram names for one detector
    // format Zn_T_u3_C6D6_runxxxxxx Zn_A_u3_SILI_runxxxxxx Zn_h_info_C6D6_runxxxxxx
    ifstream inni("Au_C6D64_HistList.dat");
    TCanvas *c=new TCanvas();

    /// replace this with the tof around a resonance, i.e. binlow is the start of the resonance, binhigh is the end of the resonance
    // you can choose any resonance, best is a large one which is usally at big tof values
    /// this will be different for Au and for Zn.

float binlow=12000000;
float binhigh=14000000;

    TLegend *legend = new TLegend(0.65, 0.38, 0.89, 0.89);
    legend->SetBorderSize(1);
    legend->SetFillColor(0);
    legend->SetTextSize(0.03);
    int count2=0;
    int count_protons=0;
    
    
    //// replace this with the number of lines in your file "list"
int lines=59;

            float value[lines];
            float value2[lines];
            float valuesili[lines];
            float valuepkup[lines];
    float staterror[lines];
    while(1){
        
        inni>>histo>>norm>>protons;
        if ( ! inni ) break;

        
        char *s;
        if(count2==0) s="hist";
        if(count2>0)s="hist same";
        count2++;
        cout<<histo<<" "<<norm<<endl;
        
        
        TH1F *h=(TH1F*)f->Get(histo);
        //TH1F *hnorm=(TH1F*)f->Get(norm);
        TH1F *hpro=(TH1F*)f->Get(protons);
        
        
        int bin1=h->FindBin(binlow); //Ge70Bins
        int bin2=h->FindBin(binhigh);
        
        // this gives the total number of count in the resonance
        
        int counter2=h->Integral(bin1,bin2);

      
    
    value[count2-1]=float(counter2)/hpro->GetBinContent(4)*7E12;
        staterror[count2-1]=sqrt(float(counter2))/hpro->GetBinContent(4)*7E12;
             value2[count2-1]=float(counter2);
             valuesili[count2-1]=float(hpro->GetBinContent(4)/hpro->GetBinContent(8));
            valuepkup[count2-1]=float(hpro->GetBinContent(4)/hpro->GetBinContent(7));
   
        h->Scale(1/hpro->GetBinContent(4)*7E12);
        h->Rebin(50);
        h->SetLineColor(count2);
        h->DrawCopy(s);
        legend->AddEntry(h,histo, "l");
    }
    float sum=0;
    float weight=0;
     float weightsum=0;
    for(int i=0;i<count2;i++){

        weight=value2[i]/(value[i]*value[i]);

       sum=sum+value[i]*weight;


        weightsum=weightsum+weight;

        cout<<value[i]<<" "<<weight<<" "<<i<<endl;
    }


cout<< sum/weightsum <<endl;
   
    float deviation=0;
double stddev=0;

    for(int i=0;i<count2;i++){
        deviation=(sum/weightsum-value[i])/sum*weightsum*100;

stddev=(sum/weightsum-value[i])*(sum/weightsum-value[i])+stddev;

        cout<<" "<<deviation<<" % "<<sqrt(value2[i])/value2[i]*100<<" % counting "<<deviation/(sqrt(value2[i])/value2[i]*100)<< "   "<<i<<endl;
    }
    deviation=sqrt(1/count2*deviation);
    cout<<count2<<" "<<sum/count2<<endl;

double stddev2=sqrt(1/(float(count2-1))*stddev);

//STOP PROTON SCALING

cout<<"std dev is "<<stddev2/(sum/weightsum)*100<<" %"<<endl;
    legend->Draw();
cout<<"                                "<<endl;
    for(int i=0;i<count2;i++){

// g++ -Wall -O3 -g List_Formatter.cpp -o List_Formatter && ./List_Formatter

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>

using namespace std;

int main(){
	
	string sample_name = "Au"; //A prefix for the TTOFSort output file histogram names, e.g. Zn1 is the first block of zinc runs
	
    int runlist[]={118026, 118027, 118028, -1};
		
	for(int detn = 1; detn <= 4; detn++){ //Loop over 4 C6D6 detectors
		
		string outname = sample_name; //Create output file name in the format Sample_C6D6#_HistList.dat
		outname += "_C6D6";
		stringstream formatter;
		formatter << detn;
		outname += formatter.str(); //This string stream now contains detector name integer converted to string
		outname += "_HistList.dat";
		
		ofstream ofile;
		ofile.open(outname.c_str());
		if(!ofile.is_open()){cout << "ERROR: Cannot open output file " << outname << endl; return 0;}
		
		for(int run_index = 0; run_index < (int)(sizeof(runlist) / sizeof(int)); run_index++){ //Loop over each of the runs
			
			if(runlist[run_index] != -1){
				
				stringstream run_str;
				run_str << runlist[run_index]; //Save run number as string
				
				string histname1 = sample_name;
				histname1 += "_T_u";
				histname1 += formatter.str();
				histname1 += "_C6D6_run";
				histname1 += run_str.str();
				
				string histname2 = sample_name;
				histname2 += "_A_u";
				histname2 += formatter.str();
				histname2 += "_SILI_run";
				histname2 += run_str.str();
				
				string histname3 = sample_name;
				histname3 += "_h_info_C6D6_run";
				histname3 += run_str.str();
				
				run_str.str(""); //Clear the stringstream
				
				ofile << histname1 << " " << histname2 << " " << histname3 << endl;
			}
		}
		
		formatter.str(""); //Clear the stringstream
		
		ofile.close();
	}
	
	return 0;
}

#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;


}