The rootfiles countsB.root and countsLi.root contain histograms of count spectra that would be expected at n_TOF for a 10B and a 6Li target, respectively. The units in y are arbitrary. This can be used to check, if we can reproduce the expected trend with our Li and B measurements. So the histograms called "counts" can be directly compared to the histograms called "energy" in Sili_deed.c (they should have the same binning). You need to scale the histogram to get a decent overlap. You can use this also to estimate the neutron energy calibration - the thermal bump at low energy will give you a good idea of the flight path length. Structure at high energy will give you a better idea on the offset.
The files 6Li_endf.root and 10B_endf.root are the original cross section files, and countrate_calc.C is the file used to produce the count spectra.
|
#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);
}
}
|