%==========================================================================
% Muon_DK_Asym_MLM_Fit.m
%==========================================================================
%
% Program written by Prof. Steven Errede 10/04/2012 10:40 hr
%
% Last updated:                          10/08/2012 16:45 hr {SME}
%
%==========================================================================
%
% Reads in MC muon decay positron x = cos(theta) data output from the
% Muon_DK_MC.m program, then carries out MLM determination of estimator 
% alpha* and +_1 sigma uncertainties on the muon weak decay asymmetry 
% parameter alpha, in fully-polarized mu+ => e+ nu_e nu_mu_bar decays.
%
%==========================================================================
clear all;
close all;

%====================================================
% Open/read in MC muon DK data formatted *.txt file:
%====================================================
%Nevts =  1000; % n.b. the total # events to read in *must* be apriori known...
%fid   = fopen('Muon_DK_Data_Alpha_eq_1p0_1K_Evts.txt','rt'); 

%Nevts = 10000; % n.b. the total # events to read in *must* be apriori known...
%fid   = fopen('Muon_DK_Data_Alpha_eq_1p0_10K_Evts.txt','rt'); 

 Nevts = 10000; % n.b. the total # events to read in *must* be apriori known...
 fid   = fopen('Muon_DK_Data_Alpha_eq_0p75_10K_Evts.txt','rt'); 

 Xevt  = zeros(1,Nevts);
 
 fprintf('\n ')
 fprintf('\n Opening & reading in MC muon DK *.txt data file... ')

%fprintf('\n ')
for j=1:Nevts; 
   line    =  fgets(fid);
   Xevt(j) = sscanf(line,'%f');
  %fprintf('\n %i %12.8f',j,Xevt(j));
end

 fclose(fid);

 fprintf('\n Closed the MC muon DK *.txt data file... ')

%====================================================

%=====================================================================
% Now make histogram of these events - do the histo plots for read-in
% data look similar/identical to those during MC event generation?
%=====================================================================

%=========================================================
% 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(01);
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(02);
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(03);
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(04);
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(05);
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(06);
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(07);
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(08);
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})');

%==========================================================================
% Determine *estimator* of alpha parameter in muon DK via the MLM
% For muon DK: PDF f(x) = 1/2*(1 + alpha*x), -1 < alpha < +1
% First, calculate ln(Like(alpha) vs. alpha over full/allowed alpha range:
%==========================================================================
 fprintf('\n ')
 fprintf('\n Calculating the Ln{Likelihood} vs. Alpha parameter... ')

Nsteps   = 30000;

alfa_lo  = -1.0;
alfa_hi  = +1.0;
del_alfa = alfa_hi - alfa_lo;
dalfa    = del_alfa/Nsteps;

Alpha  = zeros(1,Nsteps);
LnLike = zeros(1,Nsteps);

alfa  = alfa_lo;
for n = 1:Nsteps;
    log_sum = 0.0;
    for j = 1:Nevts;
        f_of_x  = 0.5*(1.0 + alfa*Xevt(j));
        
        if (f_of_x <= 0.0); % *must* NOT be negative !!!
            fprintf('\n ');
            fprintf('\n Warning! f(x) <= 0 !!! j, alpha, x, f(x) = %i %f %f %f',j,alfa,Xevt(j),f_of_x);
        end
        log_sum = log_sum + log(f_of_x);
    end;
    
     Alpha(n) = alfa;
    LnLike(n) = log_sum;
    
    alfa = alfa + dalfa;
end;

% Plot LnLike vs. Alpha:
figure(11);
plot(Alpha,LnLike,'b');
axis tight;
grid on;
xlabel('Alpha');
ylabel('Ln(Likelihood(Alpha))');
title('Muon DK MC: Ln(Likelihood(Alpha)) vs. Alpha')

%==========================================================================
% From above plot, window in on narrow range of ln(Like(alpha) vs. alpha:
%==========================================================================
Nsteps   = 10000;

alfa_lo  = 0.70;
alfa_hi  = 0.80;
del_alfa = alfa_hi - alfa_lo;
dalfa    = del_alfa/Nsteps;

Alpha  = zeros(1,Nsteps);
LnLike = zeros(1,Nsteps);
Parab  = zeros(1,Nsteps);
LnLike_max_expt = -1.0e20;

alfa  = alfa_lo;
for n = 1:Nsteps;
    log_sum = 0.0;
    for j = 1:Nevts;
        f_of_x  = 0.5*(1.0 + alfa*Xevt(j));
        
        if (f_of_x <= 0.0); % *must* NOT be negative !!!
            fprintf('\n ');
            fprintf('\n Warning! f(x) <= 0 !!! j, alpha, x, f(x) = %i %f %f %f',j,alfa,Xevt(j),f_of_x);
        end
        log_sum = log_sum + log(f_of_x);
    end;
    
    % Simultaneously, determine the global maximum of LnLike:
    % n.b. LnLike < 0 !!!
    if (log_sum > LnLike_max_expt);
        LnLike_max_expt = log_sum;
         Alpha_max_expt =    alfa;
    end
    
     Alpha(n) = alfa;
    LnLike(n) = log_sum;
    
    alfa = alfa + dalfa;
end;

fprintf('\n ');
fprintf('\n Global MAXIMUM LogLike(alpha_max),alpha_max =  %f %f',LnLike_max_expt,Alpha_max_expt);

%==========================================================================
% Determine +_ 1 sigma hi/lo points on Alpha_max by determining the hi/lo
% alpha points where LnLike(alpha) is reduced by 0.5 from LnLike(alpha_max)
%==========================================================================
  LnLike_lo = LnLike_max_expt - 0.5;
  LnLike_hi = LnLike_max_expt - 0.5;

  fprintf('\n ');
  fprintf('\n -1 Sigma LO LnLike = %f',LnLike_lo);
  fprintf('\n +1 Sigma HI LnLike = %f',LnLike_hi);
 
  Del_LnLike_lo_old = 1.0e20;
  Del_LnLike_hi_old = 1.0e20;
 
 %fprintf('\n ');
for n = 1:Nsteps;
    Del_LnLike_lo = abs(LnLike(n) - LnLike_lo);
    if ((Del_LnLike_lo < Del_LnLike_lo_old) && (Alpha(n) < Alpha_max_expt));
         Del_LnLike_lo_old = Del_LnLike_lo;
         
        %fprintf('\n n, Del_LnLike_lo, Alpha(n),LnLike(n) = %i %f %f %f',n,Del_LnLike_lo,Alpha(n),LnLike(n));
         
         LnLike_lo_expt = LnLike(n);
          Alpha_lo_expt =  Alpha(n);
    end
    
    Del_LnLike_hi = abs(LnLike(n) - LnLike_hi);
    if ((Del_LnLike_hi < Del_LnLike_hi_old) && (Alpha(n) > Alpha_max_expt));
         Del_LnLike_hi_old = Del_LnLike_hi;
         
        %fprintf('\n n, Del_LnLike_hi, Alpha(n),LnLike(n) = %i %f %f %f',n,Del_LnLike_hi,Alpha(n),LnLike(n));

        LnLike_hi_expt = LnLike(n);
         Alpha_hi_expt =  Alpha(n);
    end    
end;

fprintf('\n ');
fprintf('\n -1 Sigma  LOW   LogLike(alpha_lo), alpha_lo  = %f %f',LnLike_lo_expt,Alpha_lo_expt);
fprintf('\n +1 Sigma  HIGH  LogLike(alpha_hi), alpha_hi  = %f %f',LnLike_hi_expt,Alpha_hi_expt);

Sigma_Alpha_lo_expt  = Alpha_max_expt -  Alpha_lo_expt;
Sigma_Alpha_hi_expt  = Alpha_hi_expt  - Alpha_max_expt;

Sigma_Alpha_avg_expt = 0.5*(Sigma_Alpha_hi_expt + Sigma_Alpha_lo_expt);

fprintf('\n ');
fprintf('\n -1 Sigma_alpha_lo_expt = %f',Sigma_Alpha_lo_expt);
fprintf('\n +1 Sigma_alpha_hi_expt = %f',Sigma_Alpha_hi_expt);
fprintf('\n   <Sigma_alpha_expt>   = %f',Sigma_Alpha_avg_expt);

xg(1) = alfa_lo;
xg(2) = alfa_hi;

yg(1) = LnLike_lo;
yg(2) = LnLike_hi;

% Plot the LnLike(Alpha) vs. Alpha data in the narrowed-down Alpha range:
figure(12);
plot(Alpha,LnLike,'b',Alpha_max_expt,LnLike_max_expt,'bo',Alpha_max_expt,LnLike_max_expt,'m*',Alpha_lo_expt,LnLike_lo_expt,'bo',Alpha_lo_expt,LnLike_lo_expt,'m*',Alpha_hi_expt,LnLike_hi_expt,'bo',Alpha_hi_expt,LnLike_hi_expt,'m*',xg,yg,'k');
axis tight;
grid on;
xlabel('Alpha');
ylabel('Ln(Likelihood(Alpha))');
title('Muon DK MC: Ln(Likelihood(Alpha)) vs. Alpha')

%=================================
% Now compute the parabola curve:
%=================================
 A  = -1.0/(2.0*Sigma_Alpha_avg_expt^2);
 B  = -LnLike_max_expt;

alfa  = alfa_lo;
for n = 1:Nsteps;
     Alpha(n) = alfa;
     Parab(n) = A*(alfa - Alpha_max_expt)^2 - B;
    
     alfa = alfa + dalfa;
end

% Plot the LnLike(Alpha) vs. Alpha with parabola overlay in the narrowed-down Alpha range:
figure(13);
plot(Alpha,LnLike,'b',Alpha,Parab,'m',Alpha_max_expt,LnLike_max_expt,'bo',Alpha_max_expt,LnLike_max_expt,'m*',Alpha_lo_expt,LnLike_lo_expt,'bo',Alpha_lo_expt,LnLike_lo_expt,'m*',Alpha_hi_expt,LnLike_hi_expt,'bo',Alpha_hi_expt,LnLike_hi_expt,'m*',xg,yg,'k');
axis tight;
grid on;
xlabel('Alpha');
ylabel('Ln(Likelihood(Alpha))');
title('Muon DK MC: Ln(Likelihood(Alpha)) vs. Alpha')
legend('ln(Like)','Parabola');

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