%==========================================================================
% Muon_DK_MC.m
%==========================================================================
% Monte Carlo simulation of muon decay: mu+ ==> e+ nu_e nu_mu_bar
%
% Program written by Prof. Steven Errede 10/04/2012 09:40 hr
%
% Last updated:                          10/08/2012 16:50 hr {SME}
%
%==========================================================================
% The PDF for spin-polarized muon decay is: f(x) = 1/2*(1 + alpha*x)
% where x = cos(theta), theta = polar angle of the e+ relative to 
% the muon spin, which is || to the muon flight direction in the lab frame.
% 
% Thus, the randge of the random variable x = cos(theta) is: 
% -1 < x = cos(theta) < +1, since the range of theta is 0 < theta < pi.
%
% The corresponding CDF(x < X) = int_{x=-1}^{x=X} [f(x < X)dx] 
%                              = 1/2*[(X + 1) + 1/2*alpha*(X^2 - 1)]
% 
% n.b. for the sake of simplicity {here}, the MC simulation of a mu+ decay
% experiment is assumed to be *ideal*, i.e. it has 100% efficiency for
% observing/detecting mu+ decays over the *full* range of x = cos(theta).
%
%==========================================================================
clear all;
close all;

% Define the total # of MC events to generate per muon decay "experiment":
Nevts = 100000;

Cevt  = zeros(1,Nevts);
Xevt  = zeros(1,Nevts);

fprintf('\n');
fprintf('\n # MC muon decay events to be generated = %i \n',Nevts);

alpha = 0.75; % Asymmetry parameter for mu+ => e+ nu_e nu_mu_bar decay

%==========================================================================
%   For each MC mu+ => e+ nu_e nu_mu_bar event, we need to obtain a
% random number from the Uniform distribution U(0,1) for *this* CDF.
%
% Since the CDF(x < X) for muon DK is known: 
%
%       CDF(x < X) = 1/2*[(X+1) + 1/2*alpha*(X^2 - 1)]
%
% we can solve the resulting quadratic equation for X - and, requiring 
% X > -1, this uniquely selects out the "+" physical solution {here}:
%
%       X = [-1 + sqrt(1 - 2*alpha*(1 - 1/2*alpha - 2*CDF))]/alpha
%
% This X is in fact the corresponding random # x for *this* muon DK event!
%
%==========================================================================

%fprintf('\n');
for j=1:Nevts;
    CDFr    = random('unif', 0.0,1.0); % Random CDF(x < X)
    Cevt(j) = CDFr;
    Xevt(j) = (-1.0 + sqrt(1.0 - 2.0*alpha*(1.0 - 0.5*alpha - 2.0*CDFr)))/alpha;
   %fprintf('\n j, CDFr, Xr, = %i %f %f',j,CDFr,Xevt(j));
end;

%=========================================================
% Set up a 100 bin histogram on the x-interval [0.0,1.0]
% Define the histogram bin width, bin centers:
%=========================================================
 Mbins =  100;
 Wctr  = zeros(1,Mbins);

%Mhist = zeros(1,Mbins);
 Qhist = zeros(1,Mbins);

 Mhist_sig   = zeros(1,Mbins);
 Mhist_sighi = zeros(1,Mbins);
 Mhist_siglo = zeros(1,Mbins);
 
 Qhist_sig   = zeros(1,Mbins);
 Qhist_sighi = zeros(1,Mbins);
 Qhist_siglo = zeros(1,Mbins);
 
 Wlo   =  0.0;
 Whi   =  1.0;

 delW  = Whi - Wlo;
 dW    = delW/Mbins;

%fprintf('\n');
for i = 1:Mbins;
    Wctr(i) = Wlo + delW*((i-1)/Mbins) + 0.5*dW;
    
    % Print out the Wctr array:
    % fprintf('\n i, Wctr(i) = %i %f',i,Wctr(i));
end;

% Now histogram the above MC generated muon DK data:
Mhist = hist(Cevt,Wctr);

% Print out the Mhist array:
%fprintf('\n');
%for i = 1:Mbins;
%    fprintf('\n i, Mhist(i) = %i %i',i,Mhist(i));
%end
    
for i = 1:Mbins;
    Mhist_sig(i)   = sqrt(Mhist(i));
    Mhist_sighi(i) = Mhist(i) + Mhist_sig(i);
    Mhist_siglo(i) = Mhist(i) - Mhist_sig(i);
end;

% Normalize the MC Mhist data to turn it into a MC PDF:
for i = 1:Mbins;
    Qhist(i)       = Mhist(i)/Nevts;
    Qhist_sig(i)   = sqrt(Mhist(i))/Nevts;
    Qhist_sighi(i) = Qhist(i) + Qhist_sig(i);
    Qhist_siglo(i) = Qhist(i) - Qhist_sig(i);
end;

% Print out the Qhist array:
%fprintf('\n');
%for i = 1:Mbins;
%    fprintf('\n i, Qhist(i) = %i %f',i,Qhist(i));
%end
    
% Plot the histogrammed/binned Mevt vs X data as a bar graph:
figure(01);
bar(Wctr,Mhist,'b');
axis tight;
grid on;
xlabel('X');
ylabel('CDF M(X)');
title('Muon DK MC: Random Uniform CDF M(X) vs. X  (X = cos(theta))')

% Plot the histogrammed/binned Mevt vs X data as a normal graph:
figure(02);
plot(Wctr,Mhist,'b',Wctr,Mhist,'b.',Wctr,Mhist_sighi,'m:',Wctr,Mhist_siglo,'m:');
axis tight;
grid on;
xlabel('X');
ylabel('CDF M(X)');
title('Muon DK MC: Random Uniform CDF M(X) vs. X  (X = cos(theta))')

% Plot the histogrammed/binned Mevt vs X data with errorbars:
figure(03);
errorbar(Wctr,Mhist,Mhist_sig,Mhist_sig);
axis tight;
grid on;
xlabel('X');
ylabel('CDF M(X)');
title('Muon DK MC: Random Uniform CDF M(X) vs. X  (X = cos(theta))')

% Plot the Random Uniform CDF's PDF data w. 1-sigma error bands:
figure(04);
plot(Wctr,Qhist,'b',Wctr,Qhist,'b.',Wctr,Qhist_sighi,'m:',Wctr,Qhist_siglo,'m:');
axis tight;
grid on;
xlabel('X');
ylabel('CDF F(X)');
title('Muon DK MC: Random Uniform CDF F(X) vs. X  (X = cos(theta))')

%=========================================================
% Set up a 200 bin histogram on the x-interval [-1.0,1.0]
% Define the histogram bin width, bin centers:
%=========================================================
 Nbins =  200;
 Xctr  = zeros(1,Nbins);

%Nhist = zeros(1,Nbins);
 
 Phist = zeros(1,Nbins);
 Chist = zeros(1,Nbins);
 
 Nhist_sig   = zeros(1,Nbins);
 Nhist_sighi = zeros(1,Nbins);
 Nhist_siglo = zeros(1,Nbins);
 
 Phist_sig   = zeros(1,Nbins);
 Phist_sighi = zeros(1,Nbins);
 Phist_siglo = zeros(1,Nbins);
 
 Hist_var = zeros(1,Nbins);
 Hist_sig = zeros(1,Nbins);
 Hist_cov = zeros(1,Nbins);

 Xlo   = -1.0;
 Xhi   =  1.0;

 delX  = Xhi - Xlo;
 dX    = delX/Nbins;

%fprintf('\n');
for i = 1:Nbins;
    Xctr(i) = Xlo + delX*((i-1)/Nbins) + 0.5*dX;
    
    % Print out the Xctr array:
    % fprintf('\n i, Xctr(i) = %i %f',i,Xctr(i));
end;

% Now histogram the above MC generated muon DK data:
Nhist = hist(Xevt,Xctr);

% Print out the Nhist array:
%fprintf('\n');
%for i = 1:Nbins;
%    fprintf('\n i, Nhist(i) = %i %i',i,Nhist(i));
%end
    
for i = 1:Nbins;
    Nhist_sig(i)   = sqrt(Nhist(i));
    Nhist_sighi(i) = Nhist(i) + Nhist_sig(i);
    Nhist_siglo(i) = Nhist(i) - Nhist_sig(i);
end;

% Plot the histogrammed/binned Nevt vs x data as a bar graph:
figure(11);
bar(Xctr,Nhist,'b');
axis tight;
grid on;
xlabel('x');
ylabel('N(x)');
title('Muon DK MC: N(x) vs. x  (x = cos(theta))')

% Plot the histogrammed/binned Nevt vs x data as a normal graph:
figure(12);
plot(Xctr,Nhist,'b',Xctr,Nhist,'b.',Xctr,Nhist_sighi,'m:',Xctr,Nhist_siglo,'m:');
axis tight;
grid on;
xlabel('x');
ylabel('N(x)');
title('Muon DK MC: N(x) vs. x  (x = cos(theta))')

% Plot the histogrammed/binned Nevt vs x data with errorbars:
figure(13);
errorbar(Xctr,Nhist,Nhist_sig,Nhist_sig);
axis tight;
grid on;
xlabel('x');
ylabel('N(x)');
title('Muon DK MC: N(x) vs. x  (x = cos(theta))')

% Normalize the MC Nhist data to turn it into a MC PDF:
for i = 1:Nbins;
    Phist(i)       = Nhist(i)/Nevts;
    Phist_sig(i)   = sqrt(Nhist(i))/Nevts;
    Phist_sighi(i) = Phist(i) + Phist_sig(i);
    Phist_siglo(i) = Phist(i) - Phist_sig(i);
end;

% Print out the Phist array:
%fprintf('\n');
%for i = 1:Nbins;
%    fprintf('\n i, Phist(i) = %i %f',i,Phist(i));
%end
    
% Plot the PDF data w. 1-sigma error bands:
figure(14);
plot(Xctr,Phist,'b',Xctr,Phist,'b.',Xctr,Phist_sighi,'m:',Xctr,Phist_siglo,'m:');
axis tight;
grid on;
xlabel('x');
ylabel('f(x)');
title('Muon DK MC: PDF f(x) vs. x  (x = cos(theta))')

% Calculate CDF from Phist PDF data:
PDF_sum = 0.0;
for i = 1:Nbins;
    PDF_sum  = PDF_sum + Phist(i);
    Chist(i) = PDF_sum;
end;

% Plot the CDF data:
figure(15);
plot(Xctr,Chist,'b');
axis tight;
grid on;
xlabel('x');
ylabel('CDF(x < X)');
title('Muon DK MC: CDF(x < X) vs. X  (X = cos(theta))')

%==========================================================================
% Calculate the variance and covariance associated w/ fluctuations in the
% Nbins histogram bins of the Nevts(x) vs. x histogram {due to multinomial 
% probability distribution associated w/ the Nbins histogram bins}:
%==========================================================================
for i = 1:Nbins;
    Hist_var(i) =   Nevts*Phist(i)*Phist(i);
    Hist_sig(i) = sqrt(Hist_var(i));
    for j = 1:Nbins;
        Hist_cov(i,j) = -Nevts*Phist(i)*Phist(j); % n.b. (i,j) symmetric!
    end;
    % n.b. the diagonal elements of cov(x_i,x_j) don't exist!
    Hist_cov(i,i) = 0.0;
end;

% Plot the Nevt histogram bin variance data:
figure(16);
plot(Xctr,Hist_var,'b');
axis tight;
grid on;
xlabel('x');
ylabel('Nevt Histogram Bin var(x)');
title('Muon DK MC: Nevt Histogram of Bin Variance(x) vs. x  (x = cos(theta))')

% Plot the Nevt histogram bin sigma data:
figure(17);
plot(Xctr,Hist_sig,'b');
axis tight;
grid on;
xlabel('x');
ylabel('Nevt Histogram Bin sigma(x)');
title('Muon DK MC: Nevt Histogram Bin Sigma(x) vs. x  (x = cos(theta))')

% Plot the Nevt histogram bin covariance data as a 3-D surface:
figure(18);
surf(Xctr,Xctr,Hist_cov);
axis tight;
% Specify the "best" azimuthal and elevation angles (deg.) 
% for viewing this plot... the default angles are [-37.5,30.0]
view([115.0,35.0]);
shading interp;
xlabel('x_{I}');
ylabel('x_{J}');
zlabel('Nevt Histogram Bin cov(x_{I},x_{J})');
title ('Muon DK MC: 3-D Plot of Nevt Histogram Bin Covariance(x_{I},x_{J}) vs. x_{I} vs. x_{J})');

%==========================================================================
% Optionally write out the MC muon DK data to formatted *.txt file,
% e.g. if needed/wanted for "offline" analysis...
%==========================================================================

 Iwrite = 1; % write out MC muon DK events to *.txt file
%Iwrite = 0; % don't write out any events....
if (Iwrite == 1);
    fprintf('\n');
    fprintf('\n Writing out the muon DK data to formatted *.txt file ');
    fid = fopen('Muon_DK_data.txt','wt');

   %for j=1:Nevts; % write out all x_i = cos(theta_i) data for muon DK events...
    for j=1:10000; % only write out a sub-set of the events...
        fprintf(fid,'%12.8f \n',Xevt(j));
       %fprintf(fid,'%12.8f %12.8f \n',Cevt(j),Xevt(j));
    end

    fclose(fid);
end;

%==========================================================================
beep();
fprintf('\n');
fprintf('\n Muon_DK_MC calculations completed...');
fprintf('\n');
%==========================================================================