%------------------------------------------------------------
% 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 = '0deg_60s.Spe';
draw = 1;
NSIG = 2; % Number of standard deviations to include in the count
          % from the fitted Gaussian

data = load(filename);
D = length(data);
sig = sqrt(data);

%------------------------------------------------------------
% Optional plotting, to see what the program is doing
%------------------------------------------------------------
if (draw) 
    hold off;
    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(0:(D-1),data,sig,'o');
    xlim([left right]);
    xlabel('MCA channel number');
    ylabel('Number detected');
    headline = sprintf('Data for %s',filename);
    title(headline);
end

%------------------------------------------------------------
% 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)) + ...
                             a(4)*(exp(-((x-a(5))/a(6)).^2));
a0 = [maximum maximum_index maximum_index-HWHM_index...
      maximum*.14 maximum_index-70 maximum_index-HWHM_index];
sgn=0; dYda={@(x,a) 1};
[a,aerr,chisq,yfit,corr] = levmar(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));