function [vpp,sig,mv,bound]=histfitlandau(val,passo,inizio,fine) % function to compute and display a Landau Distribution % for INFO visit www.cern.ch/meroli % special thanks to Daniele %%% input % val --> acquired data % passo --> bin size % inizio --> start point for the fit % fine --> last point for the fit %%%output % vpp --> most probable value for the Landau Distribution % sig --> sigma for the Landau Distribution % mv --> mean value for the Landau Distribution % bound --> bounds for vpp and sig %% Example % x(1:1000)=randn(1,1000); % x(1001:2000)=randn(1,1000).*4+2; %[vpp,sig,mv,bound]=histfitlandau(x,0.3,-3,15); clf; if(nargin==1) [y1,x1]=hist(val); passo=(max(x1)-min(x1))/calcnbins(val); end; if nargin>3 [y,x]=hist_(val,passo,min(val),fine); else [y,x]=hist_(val,passo); end; x=x+passo/2; if nargin>2; inz=find(x<inizio); if isempty(inz) inz=1; end; y(1:inz(end))=0; end; a0=max(y); mpv0=min(x(find(y==a0))); sigma0=std(val); ftype=fittype('a*landau(x,mpv,sigma,0)'); fitres=fit(x,y,ftype,'StartPoint',[a0 mpv0 sigma0]); bound=confint(fitres); xfit=x(1):(passo/10):x(end); yfit=fitres(xfit); line(xfit,yfit,'LineWidth',2,'LineStyle','-','Color','r'); picco=max(yfit); mpv=xfit(find(yfit==picco)); mpv=mpv(end); mv=mean(val); errMPV=mean(abs(bound(:,2)-mpv)); errSig=mean(abs(bound(:,3)-fitres.sigma)); line([mpv mpv],[0 picco],'LineWidth',2,'LineStyle','--','Color','r'); %text(mpv,0,[' \leftarrow' num2str(mpv)],'rotation',-90,'FontSize',12'); testo={['MPV: ', num2str(mpv), '\pm', num2str(errMPV), ];['\sigma = ', num2str(fitres.sigma),'\pm', num2str(errSig),];['Entries: ' num2str(length(val))]}; %testo={['MPV: ', num2str(mpv) ];['\sigma = ', num2str(fitres.sigma)];['Entries: ' num2str(length(val))]}; annotation('textbox',[.61 .8 .1 .1],'FitHeightToText','on','Interpreter','tex','String',testo,'FontName','Century Gothic','Fontsize',18); if nargout vpp=mpv; sig=fitres.sigma; mv=mean(val); bound=confint(fitres); bound(:,1)=[]; end; end function nbins = calcnbins(x, method, minimum, maximum) error(nargchk(1, 4, nargin)); if ~isnumeric(x) && ~islogical(x) error('calcnbins:invalidX', 'The X argument must be numeric or logical.') end if ~isreal(x) x = real(x); warning('calcnbins:complexX', 'Imaginary parts of X will be ignored.'); end % Ignore dimensions of x. if ~isvector(x) x = x(:); warning('calcnbins:nonvectorX', 'X will be coerced to a vector.'); end nanx = isnan(x); if any(nanx) x = x(~nanx); warning('calcnbins:nanX', 'Values of X equal to NaN will be ignored.'); end if nargin < 2 || isempty(method) method = 'middle'; end if ~ischar(method) error('calcnbins:invalidMethod', 'The method argument must be a char array.'); end validmethods = {'fd'; 'scott'; 'sturges'; 'all'; 'middle'}; methodmatches = strmatch(lower(method), validmethods); nmatches = length(methodmatches); if nmatches~=1 error('calnbins:unknownMethod', 'The method specified is unknown or ambiguous.'); end method = validmethods{methodmatches}; if nargin < 3 || isempty(minimum) minimum = 1; end if nargin < 4 || isempty(maximum) maximum = Inf; end % Perform the calculation switch(method) case 'fd' nbins = calcfd(x); case 'scott' nbins = calcscott(x); case 'sturges' nbins = calcsturges(x); case 'all' nbins.fd = calcfd(x); nbins.scott = calcscott(x); nbins.sturges = calcsturges(x); case 'middle' nbins = median([calcfd(x) calcscott(x) calcsturges(x)]); end % Calculation details function nbins = calcfd(x) h = diff(prctile0(x, [25; 75])); %inter-quartile range if h == 0 h = 2*median(abs(x-median(x))); %twice median absolute deviation end if h > 0 nbins = ceil((max(x)-min(x))/(2*h*length(x)^(-1/3))); else nbins = 1; end nbins = confine2range(nbins, minimum, maximum); end function nbins = calcscott(x) h = 3.5*std(x)*length(x)^(-1/3); if h > 0 nbins = ceil((max(x)-min(x))/h); else nbins = 1; end nbins = confine2range(nbins, minimum, maximum); end function nbins = calcsturges(x) nbins = ceil(log2(length(x)) + 1); nbins = confine2range(nbins, minimum, maximum); end function y = confine2range(x, lower, upper) y = ceil(max(x, lower)); y = floor(min(y, upper)); end function y = prctile0(x, prc) % Simple version of prctile that only operates on vectors, and skips % the input checking (In particluar, NaN values are now assumed to % have been removed.) lenx = length(x); if lenx == 0 y = []; return end if lenx == 1 y = x; return end function foo = makecolumnvector(foo) if size(foo, 2) > 1 foo = foo'; end end sortx = makecolumnvector(sort(x)); posn = prc.*lenx/100 + 0.5; posn = makecolumnvector(posn); posn = confine2range(posn, 1, lenx); y = interp1q((1:lenx)', sortx, posn); end end function [h,intervalli]=hist_(x,passo,mi,mx,n,noplot,shift) set(gca,'Fontsize',16) if(nargin==1) passo=(max(x)-min(x))/calcnbins(x); mi=min(x); mx=max(x)+passo; n='n'; noplot=0; shift=0; end; if(nargin==2) if ischar(passo) n=passo; passo=(max(x)-min(x))/calcnbins; mi=min(x); mx=max(x)+passo; noplot=0; shift=0; else mi=min(x); mx=max(x)+passo; n='n'; noplot=0; shift=0; end; end; if(nargin==3) if ischar(mi) n=mi; mi=min(x); mx=max(x)+passo; noplot=0; shift=0; else mx=max(x)+passo; n='n'; noplot=0; shift=0; end; end; if(nargin==4) if ischar(mx) n=mx; mx=max(x)+passo; noplot=0; shift=0; else n='n'; noplot=0; shift=0; end; end; if(nargin==5) noplot=0; shift=0; end; if(nargin==6) shift=0; end; if isempty(mi) mi=min(x); end if isempty(mx) mx=max(x)+passo; end intervalli=mi:passo:mx; h=histc(x,intervalli); if (any(strcmp(n,{'n' 'N'}))) % h=h/(sum(h)*passo); % h=h/(sum(h)); end; if shift intervalli=(intervalli-passo/2)'; end % bar(intervalli,h); if (noplot==0) stairs(intervalli,h,'k','LineWidth',2); end; if nargout h=h(:); intervalli=intervalli(:); end; end function y=landau(x,mpv,sigma,norm) if (nargin<4) norm=0; end; mpshift = -0.22278298; mpv = mpv - mpshift * sigma; [r,c]=size(x); y=zeros(r,c); for(j=1:c) for(i=1:r) y(i,j)=landau_(x(i,j),mpv,sigma,norm); end; end; end function y=landau_(x,mpv,sigma,norm) p1 = [0.4259894875,-0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253]; q1 = [1.0 ,-0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063]; p2 = [0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 0.0001283617211]; q2 = [1.0 , 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714]; p3 = [0.1788544503, 0.09359161662,0.006325387654, 0.00006611667319,-0.000002031049101]; q3 = [1.0 , 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675]; p4 = [0.9874054407, 118.6723273, 849.2794360, -743.7792444, 427.0262186]; q4 = [1.0 , 106.8615961, 337.6496214, 2016.712389, 1597.063511]; p5 = [1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.94910]; q5 = [1.0 , 156.9424537, 3745.310488, 9834.698876, 66924.28357]; p6 = [1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109]; q6 = [1.0 , 651.4101098, 56974.73333, 165917.4725, -2815759.939]; a1 = [0.04166666667,-0.01996527778, 0.02709538966]; a2 = [-1.845568670,-4.284640743]; if (sigma <= 0) y=0; return; end; v = (x-mpv)/sigma; if (v < -5.5) u = exp(v+1.0); if (u < 1e-10) y=0.0; return; end; ue = exp(-1/u); us = sqrt(u); den = 0.3989422803*(ue/us)*(1+(a1(1)+(a1(2)+a1(3)*u)*u)*u); elseif(v < -1) u = exp(-v-1); den = exp(-u)*sqrt(u)*(p1(1)+(p1(2)+(p1(3)+(p1(4)+p1(5)*v)*v)*v)*v)/(q1(1)+(q1(2)+(q1(3)+(q1(4)+q1(5)*v)*v)*v)*v); elseif(v < 1) den = (p2(1)+(p2(2)+(p2(3)+(p2(4)+p2(5)*v)*v)*v)*v)/(q2(1)+(q2(2)+(q2(3)+(q2(4)+q2(5)*v)*v)*v)*v); elseif(v < 5) den = (p3(1)+(p3(2)+(p3(3)+(p3(4)+p3(5)*v)*v)*v)*v)/(q3(1)+(q3(2)+(q3(3)+(q3(4)+q3(5)*v)*v)*v)*v); elseif(v < 12) u = 1/v; den = u*u*(p4(1)+(p4(2)+(p4(3)+(p4(4)+p4(5)*u)*u)*u)*u)/(q4(1)+(q4(2)+(q4(3)+(q4(4)+q4(5)*u)*u)*u)*u); elseif(v < 50) u = 1/v; den = u*u*(p5(1)+(p5(2)+(p5(3)+(p5(4)+p5(5)*u)*u)*u)*u)/(q5(1)+(q5(2)+(q5(3)+(q5(4)+q5(5)*u)*u)*u)*u); elseif(v < 300) u = 1/v; den = u*u*(p6(1)+(p6(2)+(p6(3)+(p6(4)+p6(5)*u)*u)*u)*u)/(q6(1)+(q6(2)+(q6(3)+(q6(4)+q6(5)*u)*u)*u)*u); else u = 1/(v-v*log(v)/(v+1)); den = u*u*(1+(a2(1)+a2(2)*u)*u); end; if (not(norm)) y=den; return; end; y=den/sigma; end % end -->