function[stda,stdt,mean_snn]=sampling_stat(fy,snn,i0,t0,stat);
clear figure; warning off;
% tstep=math time step, 
% tau = peaking time (number of samples), sn=signal/noise ratio, 
% stat=statistics size
%
% Get params
[tstep,trise,tfall,gain,NPE,Ishot,G,nadc,bw,tau,T,R,fs,ts,apjit]=params();
%
% Clear
clear stda stdt; 
disp('math time step')
tstep
t0max=floor(1*tau);
x=1; y=0; t0(1)=0;
t0;
i0;
disp('Sampling period');  
s=ts*tstep
%
% Start iterations up to maxiter
maxiter=500; z=zeros(stat,maxiter); amp=z; time=z; resamp=z; restime=z; cumsn=0; sni=0; it=0;
for i=1:stat
       disp('Evt');
       i
       [mns maxsig(i)]=sample(i,tstep,fy,tau,stat);
       mns(1)=0;
       a=maxsig(i)
       b=0;
       if snn~=inf
           cumsn=cumsn+snn(i);
           sni=sni+1;
       end
     for ii=1:2;
       ii
       % Reference delayed by t0, scaled to i0
       [ys,yps]=ref_delayed(a,b,fy);
       mns=mns(1:size(ys,2));
       % figure(11);
       scalex=(1:size(ys,2))*ts*tstep;
       % plot(scalex,ys,'b'); hold;  plot(scalex,mns,'r');
       % figure(12); plot(scalex,yps,'g');
       % pause;
       [x y]=iter(ys,yps,mns);
       % iter gets as x the fraction mns/ys, y the delay in unit of samples
       restime(i,ii)=y;
       a=a*x;
       time(i,ii)=sum(restime(i,:));
       % Sum to get the first guess and add-up small further corrections
       b=sum(restime(i,:));
       it=it+1;
     end
   ii=3;
   % pause
   while (abs(amp(i,ii-1)-amp(i,ii-2))>mean(maxsig)*1e-5 | abs(time(i,ii-1)-time(i,ii-2))>=1);
       a
       b
       [ys,yps]=ref_delayed(a,b,fy);
       % figure(1000); scalex=(1:size(ys,2))*ts*tstep;
       % title('ref_delayed');
       % plot(scalex,ys,'b'); hold; bar(scalex,mns,'r'); title('while');
       % return
       [x y]=iter(ys,yps,mns);
       it=it+1;
       resamp(i,ii)=a*x; restime(i,ii)=y;
       amp(i,ii)=resamp(i,ii);  a=amp(i,ii);
       time(i,ii)=sum(restime(i,:));
       b=sum(restime(i,:)); 
       if time(i,ii)> 1e5
           disp('time = NaN'); time(i,ii);
           % figure(4); bar(ys,'b'); hold; bar(mns,'r');
       end
       ii=ii+1; 
       if ii>=maxiter
        disp('No convergence');
        i
        ii
        figure(4); plot(ys,'b'); hold; plot(mns,'r');
        title('No convergence'); return
       end
      ii
   end
   if i==stat
       sz=size(ys,2); scalex=ts*(1:sz)*tstep; s=ts*tstep;
       figure(9); plot(scalex,ys,'r'); hold;  bar(scalex,mns,'b');
       xlabel('s'); ylabel('Amplitude'); 
       title('Sampled measurement and reference');
       % pause;
   end
   aa(i)=amp(i,ii-1);
   tt(i)=-time(i,ii-1);
   tt(i)=tt(i)-1;
   timeabs(i)=tt(i)-t0(i)-tt(1)
   timeabs(1:i);
   ab(i)=(aa(i)-maxsig(i))/maxsig(i);
   % disp('convergence reached')
   time; 
   amp;
end
disp('average iterations')
it/stat
figure(10); subplot(1,2,1); hist(ab(2:stat));
subplot(1,2,2); hist(tstep*(timeabs));
title('Error distributions')
% 
% Remove NaN's
%
p=1;
for i=1:stat-p+1
  if isnan(tt(i));
    tt(i:stat-1)=tt(i+1:stat); tt=tt(1:stat-p);
    p=p+1;
  end
end
stdt0=std(t0(2:stat))*tstep
stdt=std(tstep*timeabs(2:stat))
stda0=std(i0(2:stat))/mean(maxsig)
stda=std(ab(2:stat))
disp('mean_snn');
mean_snn=cumsn/sni
return
