%========================================================================== % 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'); %==========================================================================