close all; clear all; clc;
which fit -all; % curvefit is need!
% First we set the parameters of the program:
N_exp = 10000;          % Number of "experiments" (i.e. drawing from random distributions)
save_plots = false;
rng(42)
s = rng;

mu1   =  3.5;
sig1  =  0.4;
mu2   =  0.8;
sig2  =  0.2;
rho12 =  0.5;           %Correlation parameter!

if ~(abs(rho12) <= 1.0)
    error('Correlation factor not in interval [-1,1]');
end

% Define the parameters needed for the transformation:
theta = 0.5 * atan( 2.0 * rho12 * sig1 * sig2 / ( sig1^2 - sig2^2 ) );
sigu = sqrt( abs( (((sig1*cos(theta))^2) - (sig2*sin(theta))^2 ) / ( (cos(theta))^2 - sin(theta)^2) ) );
sigv = sqrt( abs( (((sig2*cos(theta))^2) - (sig1*sin(theta))^2 ) / ( (cos(theta))^2 - sin(theta)^2) ) );

% Produce random numbers with the (possible) correlation:
u = normrnd(0.0, sigu, [N_exp,1]);
v = normrnd(0.0, sigv, [N_exp,1]);
x1_all = mu1 + cos(theta)*u - sin(theta)*v;
x2_all = mu2 + sin(theta)*u + cos(theta)*v;
x12_all = [x1_all, x2_all];

y_all = x1_all - 2*x2_all;         % Silly formula - you have to put this in yourself!

plot([0,6],[0,0], '--'); hold on;
Xedges = [0:6.0/120:6.0]; Yedges = [-1.0:4/80:3.0];
histogram2(x1_all, x2_all,Xedges, Yedges, 'DisplayStyle','tile');
xlabel('x'); ylabel('y'); colorbar;
d.Entries = ['Entries: ' num2str(length(x12_all))];
d.x_mean = ['Mean x: ' num2str(mean(x1_all))];
d.y_mean = ['Mean y: ' num2str(mean(x2_all))];
d.x_std= ['Std x: ' num2str(std_with_ddof(x1_all, 1))];
d.y_std = ['Std y: ' num2str(std_with_ddof(x2_all, 1))];
txt = struct2cell(d);
text(0.5,2.3,txt);

xmin = 0.0; xmax = 4.0; bins=100.;
edges = [xmin:(xmax-xmin)/bins:xmax];
figure; h = histogram(y_all,edges);
bin_centers = (edges(1:length(edges)-1)+edges(2:length(edges)))/2;
counts = h.Values; s_counts = sqrt(counts);

txt = ['    Mean y_all = ' num2str(mean(y_all)), '    std y_all = ' num2str(std_with_ddof(y_all, 1))];
disp(txt);
gaussEqn = 'N*1.0/(sigma*sqrt(2*pi))*exp(-0.5*((x-mu)/sigma)^2)';
startPoints = [5000 1.0 2.0];
f = fit(bin_centers(counts>0)', counts(counts>0)', gaussEqn,'Start', startPoints, 'Weights', s_counts(counts>0));
hold on; plot(f,bin_centers, counts)
xlabel('x'); ylabel('#');
your_result = [];
chi2=0;
for i=1:numel(counts)
    if counts(i) == 0 
            continue
    end
    obs = gaussian(bin_centers(i), f.N, f.mu, f.sigma);
    chi2=chi2+(obs-counts(i))^2/s_counts(i)^2;
end

uncertainties = confint(f, 0.95); % 

clearvars d
d.Entries = ['Entries: ' num2str(length(y_all))];
d.mean = ['Mean: ' num2str(mean(y_all))];
d.std = ['Std: ' num2str(std_with_ddof(y_all, 1))];
d.chi2 = ['chi2: ' num2str(chi2)];
d.prob = ['Prob: ' 'calculate and input'];
d.NDF = ['NDF: ' num2str(length(counts(counts>0)))];
d.N = ['N: ' num2str(f.N) ' +/- ' num2str(round(uncertainties(2)-uncertainties(1),3))];
d.mu = ['\mu: ' num2str(f.mu) ' +/- ' num2str(round(uncertainties(4)-uncertainties(3),3))];
d.sigma = ['\sigma: ' num2str(f.sigma) ' +/- ' num2str(round(uncertainties(6)-uncertainties(5),3))];

txt = struct2cell(d);
text(min(bin_centers)+std(bin_centers)/40,max(counts)-std(counts),txt);

% Standard Deviation with degrees of freedom
function std = std_with_ddof(array, ddof)
    array_mu = mean(array);
    arr_sum = sum((array-array_mu).^2);
    std = sqrt(arr_sum/(length(array)-ddof));
end

% Gaussian function
function gauss_array = gaussian(x, N, mu, sigma)
    gauss_array = N * 1.0 / (sigma*sqrt(2*pi)) * exp(-0.5* (x-mu)^2/sigma^2);
end
















