%==========================================================================
% Muon_DK_Asym.m
%==========================================================================
%
% Analytic program written by Prof. Steven Errede 10/04/2012 08:40 hr
%
% Last updated:                                   10/08/2012 16:40 hr {SME}
%
%==========================================================================
%
% Plots the PDF(x) and the CDF(x<X) for muon decay mu+ ==> e+ nu_e nu_mu_bar
%
% 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)]
% 
%==========================================================================
clear all;
close all;

alpha_m1 = -1.00;
alpha_oo =  0.00;
alpha_mc =  0.75;
alpha_p1 = +1.00;

Npts     = 20000;

Xmu      = zeros(1,Npts);

PDF_m1   = zeros(1,Npts);
PDF_oo   = zeros(1,Npts);
PDF_mc   = zeros(1,Npts);
PDF_p1   = zeros(1,Npts);

CDF_m1   = zeros(1,Npts);
CDF_oo   = zeros(1,Npts);
CDF_mc   = zeros(1,Npts);
CDF_p1   = zeros(1,Npts);


Xlo  = -1.0;
Xhi  = +1.0;

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

X    = Xlo;
for j = 1:Npts;
    Xmu(j) = X;
    
    PDF_m1(j) = 0.5*(1.0 + alpha_m1*X);
    PDF_oo(j) = 0.5*(1.0 + alpha_oo*X);
    PDF_mc(j) = 0.5*(1.0 + alpha_mc*X);
    PDF_p1(j) = 0.5*(1.0 + alpha_p1*X);
    
    CDF_m1(j) = 0.5*((X + 1.0) + 0.5*alpha_m1*(X^2 - 1.0));
    CDF_oo(j) = 0.5*((X + 1.0) + 0.5*alpha_oo*(X^2 - 1.0));
    CDF_mc(j) = 0.5*((X + 1.0) + 0.5*alpha_mc*(X^2 - 1.0));
    CDF_p1(j) = 0.5*((X + 1.0) + 0.5*alpha_p1*(X^2 - 1.0));
    
    X = X + dX;
end

figure(01);
plot(Xmu,PDF_m1,'m',Xmu,PDF_oo,'g',Xmu,PDF_mc,'k',Xmu,PDF_p1,'b');
axis tight;
grid on;
xlabel('x = cos(theta)');
ylabel('f(x)');
title('Muon Decay Asymmetry: f(x) vs. x');
legend('{\alpha} = -1','{\alpha} =  0','{\alpha} = 0.75','{\alpha} = +1');

figure(02);
semilogy(Xmu,PDF_m1,'m',Xmu,PDF_oo,'g',Xmu,PDF_mc,'k',Xmu,PDF_p1,'b');
axis tight;
grid on;
xlabel('x = cos(theta)');
ylabel('f(x)');
title('Muon Decay Asymmetry: f(x) vs. x');
legend('{\alpha} = -1','{\alpha} =  0','{\alpha} = 0.75','{\alpha} = +1');



figure(03);
plot(Xmu,CDF_m1,'m',Xmu,CDF_oo,'g',Xmu,CDF_mc,'k',Xmu,CDF_p1,'b');
axis tight;
grid on;
xlabel('X');
ylabel('CDF(x < X)');
title('Muon Decay Asymmetry: CDF(x < X) vs. X');
legend('{\alpha} = -1','{\alpha} =  0','{\alpha} = 0.75','{\alpha} = +1');

figure(04);
semilogy(Xmu,CDF_m1,'m',Xmu,CDF_oo,'g',Xmu,CDF_mc,'k',Xmu,CDF_p1,'b');
axis tight;
grid on;
xlabel('X');
ylabel('CDF(x < X)');
title('Muon Decay Asymmetry: CDF(x < X) vs. X');
legend('{\alpha} = -1','{\alpha} =  0','{\alpha} = 0.75','{\alpha} = +1');

%==========================================================================
% Make plots specifically for MC muon DK expt with alpha = 0.75:
%==========================================================================
Npts     = 2000;
alpha_mc = 0.75;
Xmc      = zeros(1,Npts);
PDF_mc   = zeros(1,Npts);
CDF_mc   = zeros(1,Npts);

Xlo  = -1.0;
Xhi  = +1.0;

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

X    = Xlo;
for j = 1:Npts;
    Xmc(j) = X;
    
    PDF_mc(j) = 0.5*(1.0 + alpha_mc*X);

    CDF_mc(j) = 0.5*((X + 1.0) + 0.5*alpha_mc*(X^2 - 1.0));
    
    X = X + dX;
end

figure(11);
plot(Xmc,PDF_mc,'b');
axis tight;
axis([-1.0 1.0 0.0 1.0]);
grid on;
xlabel('x = cos(theta)');
ylabel('f(x)');
title('Muon Decay Asymmetry: f(x) vs. x [{\alpha} = 0.75]');

% A "randomly" chosen point on CDF(x < X) vs. X graph:
xg0(1) =  0.6950;
yg0(1) =  0.7506;

xg1(1) =  0.6950;
xg1(2) =  0.6950;

yg1(1) = -1.0000;
yg1(2) =  0.7506;

xg2(1) = -1.0000;
xg2(2) =  0.6950;

yg2(1) =  0.7506;
yg2(2) =  0.7506;

figure(13);
plot(Xmc,CDF_mc,'b',xg0,yg0,'bo',xg0,yg0,'m*',xg1,yg1,'m',xg2,yg2,'m');
axis tight;
axis([-1.0 1.0 0.0 1.0]);
grid on;
xlabel('X');
ylabel('CDF(x < X)');
title('Muon Decay Asymmetry: CDF(x < X) vs. X [{\alpha} = 0.75]');

%==========================================================================
beep();
fprintf('\n');
fprintf('\n Muon_DK_Asym calculations finished \n');
%==========================================================================