%% Make a data set and use gaussian fitting tool to get simulated counts
%file path to atom probe lab source folder
addpath();
%set seed
rng(1)
%newPos = 'PeakWidth20_20200504.pos'
% make masses %totalAtomCount, totalNoiseCount, ionInput(Name   at%  wth  +   ++   +++ ++++), peakShape, minimumAbundance
[m,p2,ionNames, eleTable] = massSpecSimulator(1E6, 0, {'Zr',83,20,0,1;'ZrH',16.72,20,0,1;'Nb',0.28,20,0,1;}, 'norm');

x = (20+20).*rand(1000000,1)-20;
y = (20+20).*rand(1000000,1)-20;
z = (70-0).*rand(1000000,1)+0;


%fileName = newPos
%savepos2(x',y',z',m',fileName);

% "overlap matrix"? non-overlapped peaks to overlapped peak counts
M = [1,0,0,0,0,0,0,0,0,0,0;0,1,1,0,0,0,0,0,0,0,0;0,0,0,1,1,0,0,0,0,0,0;0,0,0,0,0,1,1,0,0,0,0;0,0,0,0,0,0,0,1,0,0,0;0,0,0,0,0,0,0,0,1,0,0;0,0,0,0,0,0,0,0,0,1,0;0,0,0,0,0,0,0,0,0,0,1];

p2 = sortrows(p2,1);

simulatedCounts = M*p2(:,2);



% Gaussian fitting:
%% Fit: 'untitled fit 1'.
[y,x]=histcounts(m,10000);
y(end+1)=y(end);

[xData, yData] = prepareCurveData( x, y );

% Set up fittype and options.
ft = fittype( 'gauss8' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Lower = [-Inf -Inf 0 -Inf -Inf 0 -Inf -Inf 0 -Inf -Inf 0 -Inf -Inf 0 -Inf -Inf 0 -Inf -Inf 0 -Inf -Inf 0];
opts.StartPoint = [1213 44.9452127 0.0540620037016038 525 45.44988606 0.054238810107707 512.389509753997 46 0.0737770546323975 479 46.5 0.0519399104695197 413 47 0.0495630284413687 296.253364768199 47.5 0.109762847680175 239.043834338999 48 0.0709510409331983 212.485024036864 48.5 0.066676937565712];

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
% plot fit model:
figure
plot(xData,yData,xData,fitresult(xData))
legend('Data','Fitted')
xlabel('m/z (Da)')
ylabel('Ion Count')
f = reshape(coeffvalues(fitresult)',3,[])';

% sort peaks in f
f = sortrows(f,2);

quantifiedCounts = f(:,1);

[simulatedCounts./sum(simulatedCounts) quantifiedCounts./sum(quantifiedCounts) OP.rangeCounts(:,1)./sum(OP.rangeCounts(:,1))]

figure
mass2charge=(45:0.5:48.5);
bar(mass2charge,[simulatedCounts./sum(simulatedCounts) quantifiedCounts./sum(quantifiedCounts) OP.rangeCounts(:,1)./sum(OP.rangeCounts(:,1))])
legend('Simulated','Fitted','AtomProbeLab')
ylabel('Fraction of counts')
xlabel('m/z (Da)')

figure
plot(xData,yData)
xlabel('m/z (Da)')
ylabel('Ion Count')

figure
plot(xData, log(yData))
xlabel('m/z (Da)')
ylabel('log(Ion Count)')