% Reads Ascii gamma ray spectral file from Hadr03 for thermal capture
% data is three columns of ascii numbers representing a histogram of gamma ray numbers vs energy from a run of Hadr03
% example Hadr03 macro is Ni58_nCapture.mac and the ascii output lines in it are:
% /analysis/h1/set 2  10000  0. 10000 keV	#gamma
% /analysis/h1/setAscii 2
% Plot histogram of data
% J.E.McFee 2020-06-15, 2021-12-29 (modified from ReadAscii_Scint_LaCl3Ce_Spectrum_AddBroadening_HistPlot.m)

format long g
format compact

%main program

% read the data
% 3 columns separated by spaces: bin number; energy (keV); counts in bin
% for Gd155, Gd157 capture I have generated 10^7 0.025eV neutrons with 10000 bins of width 1keV

%input_file = './ascii_files/Gd155_nCapture_107_ed.ascii';
input_file = './ascii_files/Gd157_nCapture_107_ed.ascii';
%input_file = './ascii_files/Gd157_nCapture_107_default.ascii';

%lineCountCommand = ['wc -l', ' ', input_file]; % use Linux system command to get number of lines in file
%[s,w] = system(lineCountCommand);
%[NLineString Fname] = strtok( w, ' ');
%NLines = str2num(NLineString);
%number_events = NLines;

input_file1 = strcat(input_file, ".gz");
file_list = gunzip(input_file1);
input_file = file_list{1};

if(exist(input_file)!=2)
  fprintf("File does not exist! %s\n",input_file);
end

p1= load("-ascii",input_file);
nrows = size(p1,1);
ncols = size(p1,2);
fprintf("File name=%s rows %d cols %d\n", input_file, nrows, ncols);  

file_list = gzip(input_file);
if(exist(input_file)==2)
  delete(input_file);
  %fprintf("File deleted %s\n",input_file);
end

if (ncols == 3)
  y1 = p1(:,1);     % bin number
  y2 = p1(:,2);     % energy (keV)
  y3 = p1(:,3);     % counts in bin
else
  printf("Columns ~= 3: rows %d cols %d\n", nrows, ncols);
endif
  


%figure(1); plot(y2,y3,".")  % Plot histogram
figure(1); bar(y2,y3)  % Plot histogram
xlabel('Gamma Ray Energy (keV)');       %  add axis labels and plot title
ylabel('Frequency of occurrence');
title('Gamma Ray Energy Spectrum');

%figure(2); plot(y2,y3,".")  % Plot log histogram
figure(2); bar(y2,y3)  % Plot log histogram
warning ("off", "Octave:negative-data-log-axis");
set (gca, 'Yscale', 'log');
xlabel('Gamma Ray Energy (keV)');       %  add axis labels and plot title
ylabel('Frequency of occurrence');
title('Gamma Ray Energy Spectrum');

Eg = y2(y3~=0); % get rid of zero entries
Ig = y3(y3~=0); % get rid of zero entries
Eg1 = Eg.*0.001;
Ig1 = Ig./Ig(length(Ig)).*100.0;

%figure(3); plot(Eg,Ig1,".")  % Plot histogram
figure(3); bar(Eg,Ig1)  % Plot histogram
xlabel('Gamma Ray Energy (keV)');       %  add axis labels and plot title
ylabel('Frequency of occurrence');
title('Gamma Ray Energy Spectrum zeros removed');

% output file of nonzero energy, intensity pairs
%fid = fopen('Ni58_spectrum_sdn_Compacted.txt', 'w');
%fprintf(fid, '%u %u \n', [Eg1 Ig1]');
%fclose(fid);