Add files via upload

muon.m

@@ -0,0 +1,214 @@
+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')
+
+
+