muon.m

num=xlsread('muon.xlsx');
x=(num(:,1));
x=x./1000;
step=0.4;
X=0:step:max(x);
figure(1)
h=histogram(x,X);
xlabel('Decay Time (us)')
ylabel('Observed Counts')
hold on
y=h.Values;
XX=(step/2):step:max(x);
XX=XX(1:end-1);
errorbar(XX,y,sqrt(y),'o')


%linear and nonlinear fit before background subtraction
%plot(XX,log(y),'o')
fun1=@(v) v(1)*exp(-v(2)*XX)-y;
x0=[4000,1];
[v,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(fun1,x0);


Jacobian = full(jacobian); 
varv=resnorm*inv(Jacobian'*Jacobian)/length(XX);
stdv=sqrt(diag(varv)); 


Y1=v(1)*exp(-v(2)*XX);
plot(XX,Y1,'g')
legend('Observed Decays','Uncertainty','Best Fit Without Background Subtraction')
title('Nonlinear Regression with no Consideration of Bacground Effects')
ylabel('Bin Count')
xlabel('Decay Time (microseconds)')
%hold on
%plot(XX,y)



%%
%linear and nonlinear fit after background subtraction
clearvars v resnorm residual stdv jacobian

figure(2)
count=sum(y(XX>12));
avg=count/(length(y(XX>12)));
Y=y-avg;
Y(Y<0)=0;
%plot(XX,log(Y),'o')
h=histogram(x,X);
xlabel('Decay Time (us)')
ylabel('Observed Counts')
hold on
y=h.Values;
errorbar(XX,y,sqrt(y),'o')




%fun2=@(v) v(1)*exp(-v(2)*XX)+v(3);%-y;
%v=[50,5,1];
%fun2=@(v,x) (v(1)*exp(-v(2)*x)+v(3));
%nlinfit(XX,y,fun2,X0)

fun2=@(v) v(1)*exp(-v(2)*XX)+v(3)-y;
X0=[4000,5,10];

[v,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(fun2,X0);
ci=nlparci(v,residual,'Jacobian',jacobian);
v=v';


Jacobian = full(jacobian); 
varv=resnorm*inv(Jacobian'*Jacobian)/length(XX);
stdv=sqrt(diag(varv)); 





Y=v(1)*exp(-v(2)*XX)+v(3);
YY=log(v(1))-v(2)*XX+log(v(3));
plot(XX,Y,'g')
legend('Observed Decays','Uncertainty','Best Fit With Background Subtraction')
title('Nonlinear Regression with Consideration of Bacground Effects')
ylabel('Bin Count')
xlabel('Decay Time (microseconds)')
%hold on 
%plot(XX,y)





%% Maximum likelihood
%without background subtraction
lambda=0.275:0.0001:0.282;
%lambda=0.1:0.01:0.6;
P=zeros(length(lambda),length(x));
mP=zeros(1,length(lambda));
for i=1:length(lambda)
    for j=1:length(x)
        P(i,j)=lambda(i)*exp(-lambda(i)*x(j));
    end
end
P=log(P);

for i=1:length(lambda)
    mP(i)=sum(P(i,:));
end

figure(5)
plot(lambda,mP)
hold on
yy=(max(mP(:))-0.5)*ones(1,length(lambda));
plot(lambda,yy)
vs=(mP((mP-yy>-0.04 & mP-yy<0.04)));
Vs=zeros(1,2);
for i=1:length(vs)
    Vs(i)=lambda(mP==vs(i));
end
plot(Vs,yy(1:2),'o')


lambda1=lambda(mP==max(mP));
siga=abs(lambda1-Vs(1));
sigb=abs(lambda1-Vs(2));
yy1=get(gca,'ylim');
plot([lambda1,lambda1],yy1);
title('Maximum Likelihood Approach with no Consideration of Background')
xlabel('Lambda (s^-1)')
ylabel('Sum of log of Probabilities')
legend('Sum of log of Probabilities','Maximum-0.5')



%% with background subtraction
lambda=0.48:0.0002:0.50;
b=0.0095:0.00002:0.0101;

P=zeros(length(lambda),length(b),length(x));
mP=zeros(length(lambda),length(b));
for i=1:length(lambda)
    for j=1:length(b)
        for k=1:length(x)
            P(i,j,k)=lambda(i)*(1-b(j)*20)*exp(-lambda(i)*x(k))+b(j);
        end
    end
end
P=log(P);

for i=1:length(lambda)
    for j=1:length(b)
        mP(i,j)=sum(P(i,j,:));
    end
end

figure(6)
[XXX,YYY]=meshgrid(b,lambda);
surf(XXX,YYY,mP)
xlabel('Background')
ylabel('Lambda (s^-1)')
zlabel('Sum of log of Probabilities')
hold on

M=max(mP(:));
[row,col]=find(mP==M);


lambda2=lambda(row);
background=b(col);


%% Uncertainty
L=M-0.5;
[row,col]=find(mP>L-0.03 & mP<L+0.03);
xp=b(col);
yp=lambda(row);
zp=L*ones(length(row),length(row));
xp=xp'; yp=yp';
k=boundary(xp,yp);
LL=length(k);
zp2=L*ones(LL,1);
plot3(xp(k),yp(k),zp2,'r','Linewidth',2)
title('Maximum Likelihood Approach with Consideration of Background')
hold off


figure(7)
contourf(YYY,XXX,mP)
hold on
scatter(yp,xp,'b')
xlabel('Lambda (s^-1)')
ylabel('Background')
zlabel('Sum of log of Probabilities')


yy1=get(gca,'ylim');
h1=plot([lambda2,lambda2],yy1);
sig1=abs(lambda2-min(yp));
sig2=abs(lambda2-max(yp));

ul=xp(yp==max(yp));
ll=xp(yp==min(yp));

h2=plot([lambda2-sig1,lambda2],[ll(1),ll(1)]);

h3=plot([lambda2+sig2,lambda2],[ul(3),ul(3)]);
%legend('','''Value of Lambda','Sigma 1','Sigma 2')
legend([h1 h2 h3],{'Value of Lambda Determined','Sigma (-)','Sigma (+)'})
title('Uncertainties in the Multidimensional Maximum Likelihood Approach')