%------------------------------------------------------------
% count_mca
%   Yields the number of hits from MCA histogram in a 
%   systematic way. Internal knob includes how many standard
%   deviations from the main hump to take.
%------------------------------------------------------------
% Usage: N = count_mca('20deg_10min.Spe',1);
%------------------------------------------------------------
filename = '10deg_60_50.Spe';
draw = 1;
NSIG = 2; % Number of standard deviations to include in the count
          % from the fitted Gaussian

loadeddata = load(filename);

%------------------------------------------------------------
% Different binning
%------------------------------------------------------------
L = length(loadeddata);
k = 5;
D = floor(L/k);
data = zeros(1,D);
sig = zeros(1,D);
for i = 1:D
    data(i) = sum(loadeddata((k*(i-1)+1):(k*i)));
    if (data(i) == 0)
        sig(i) = 1;
    else
        sig(i) = sqrt(data(i));
    end
end
bins = 0:k:(k*(D-1));



%------------------------------------------------------------
% Optional plotting, to see what the program is doing
%------------------------------------------------------------
if (draw) 
    figure;
%     for i = 1:1:D
%         if(data(i) ~= 0)
%             break;
%         end
%     end
%     left = i;
%     for i = D:-1:1
%         if(data(i) ~= 0)
%             break;
%         end
%     end
%     right = i;
    errorbar(bins,data,sig,'o');
%     xlim([left right]);
    xlim([400 1500]);
    fs = 15;
    xlabel('MCA channel number','fontsize',fs);
    ylabel('Frequency','fontsize',fs);
    title('Histogram of \alpha-particle scattering at 10^\circ over 60.50s','fontsize',fs);
end

function_gaussian   = @(x,a) a(1)*(exp(-((x-a(2))/a(3)).^2));
a0 = [20 1127 85];
sgn=0; dYda={@(x,a) 1};
l = 180;
r = 260;
[a,aerr,chisq,yfit,corr] = levmarold(bins(l:r)',data(l:r)',sig(l:r)',...
                                    function_gaussian,a0,dYda,sgn);
hold on;
x = linspace(400,1500);
y = function_gaussian(x,a);
plot(x,y,'r','linewidth',2);
ylim([0 35]);
text(480,28,'Gaussian fit over channels 900:1300','fontsize',fs);
text(480,26,'yields \chi^2/dof = 0.93 and','fontsize',fs);
text(550,22,'\mu = 1123 \pm 3','fontsize',fs);
text(550,20,'\sigma = 76 \pm 3','fontsize',fs);


% 
% %------------------------------------------------------------
% % Fitting: Get approximate parameters
% %------------------------------------------------------------
% for i = 1:D
%     if(sig(i)==0)
%         sig(i) = 1;
%     end
% end % Arbitrarily set std of 0 to be 1
% mca = (0:(D-1))';
% [maximum maximum_index] = max(data);
% for i = 1:D
%     if (data(i) > maximum/2)
%         break;
%     end
% end % Get approximate half maximum index
% HWHM_index = i;
% function_gaussian   = @(x,a) a(1)*(exp(-((x-a(2))/a(3)).^2));
% a0 = [maximum maximum_index maximum_index-HWHM_index];
% sgn=0; dYda={@(x,a) 1};
% [a,aerr,chisq,yfit,corr] = levmarold(mca,data,sig,...
%                                     function_gaussian,a0,dYda,sgn);
% low_index  = round(a(2)-NSIG*a(3));
% high_index = round(a(2)+NSIG*a(3));
% if (draw)
%     hold on;
%     plot(mca,yfit,'r','linewidth',3);
% %     plot([low_index low_index],[0 maximum],'g','linewidth',3');
% %     plot([high_index high_index],[0 maximum],'g','linewidth',3');
%     right = max([high_index right]);
%     left  = min([low_index left]);
%     xlim([left right]);
% end
% 
% % Now we count and report
% N = sum(data(low_index:high_index));