peakwidth_pf(s, startsd, nsubdiv)

function imxs = peakwidth_pf(s, startsd, nsubdiv)
% 
% This function estimated the left half-width of the peaks accross TOF
% record
%
% OUTPUT:
%
% imxs - point of the maximum, and estimated width (two columns);
%
% INPUT:
% s - signal equally sampled in time;
% startsd - starting point for subdivision;
% nsubdiv - number of subdivisions for peak width estimate;
% 
% USAGE:
% imxs = peakwidth_pf(s, startsd, nsubdiv);
% Dependency: none
%

drange_thr = 20;
Nr=size(s,1);
Nc=size(s,2);

psdiv=floor((Nr-startsd)/nsubdiv);
mxs=[0];sgmmxs=[0];
for k=1:nsubdiv-1
%startsd+1+(k-1)*psdiv
%startsd+k*psdiv
ss=s(startsd+1+(k-1)*psdiv:startsd+k*psdiv);
[mx,imx]=max(ss);
sgm=gpdf_hw(ss, imx); % find sigma of a peak
if sgm>0
  mxs=[mxs;imx+startsd+1+(k-1)*psdiv];
sgmmxs=[sgmmxs;sgm];
end;
end;
%whos mxs
%whos sgmmxs

mxs=mxs(2:end); sgmmxs=sgmmxs(2:end);

imxs=[mxs,sgmmxs];

[mmx,immx]=max(s(imxs(:,1)));
%% look only at big signals

r1=imxs(:,1)-imxs(:,2); r2=imxs(:,1)+imxs(:,2);
maxsp=r1*0;
for k1=1:length(r1)
maxsp(k1)=max(s(r1(k1):r2(k1)));    
end;

kimx=find((s(imxs(:,1)).*imxs(:,2)>mmx*imxs(immx,2)/drange_thr) & (s(imxs(:,1))> 0.5*maxsp(:)));
mxs=mxs(kimx); sgmmxs=sgmmxs(kimx);
imxs=imxs(kimx,:);

% sigmafit=1; sigmafit=polyfit(mxs,sgmmxs,2);

t=1:Nr; t=t';
% plot(t, s*max(sgmmxs)/max(s), t(mxs), sgmmxs, 'rs');

return;