Here are all the files needed for calculation of the yield. There is some information that still needs to be added to the energyandyield.c file, for example the atomic mass of your target in amu, the neutron separation energy for the compound nucleus in MeV (this is needed due to using the Weighting technique). Also, you need to add some lines  to read the tof histogram, which will then be converted in a neutron energy hisogram and then divided by the flux.

For the conversion to neutron energy we assume an approximate flight path for now. We will determine this more accurately once we have the yields and fit resonances. Emma will need to use a different flux than Annie, because we have a different spallation target. Emma's flux is preliminary, so we also should look at the SILI detectors at some point.

There are some more instances where names for rootfiles and ascii files need to be added. In the first instance, please try and add the missing information and try running the code. I haven't tested it so there may be still problem.  To run it you will need to download the 3 rootfiles into the same folder (or move somewhere else and adjust the path name in the .c file.

PLEASE set variable THRESHSCALER to 1 for now. We dont know yet what it is :-)

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


}