function[stda,stdt,mean_snn,timeabs]=sampling_stat_mcp(ffy,ref,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,taus,T,R,fs,ts,apjit]=params();
%
% Clear
clear stda stdt; 
disp('math time step')
tstep
% Minimum delay
% t0max=floor(4*tau);
% new
t0max=floor(1*tau);
%
% Silicon;  telec=5e-9; thole=20e-9;
% MCP;      trise=2e-10; tfall=2*trise
%
% Charge distributions
% Non central dist from Matlab 1 MIP MPV
% hist(random('nct',10,10000,1,20000),50);
%
% i0=(2.7e-10)*random('nct',10,1000,1,nn)
% figure(30); hist(i0); return
%
% disp('t0 distribution');
%
% Check with fixed amplitude i0=1*ones(1,nn)
% Check with fixed t0; t0=10*ones(1,nn)
% Stat noisy measurements at sn, amplitude i0(i), shifted by t0(i)
x=1; y=0; t0(1)=0;
t0;
i0;
disp('Sampling period');  
s=ts*tstep; n=size(ffy,2)*ts; 
% ref(1,:)=resample(ffy(stat+1,:),20,1,1000);
% ref(2,:)=resample(ffy(stat+2,:),20,1,1000);
size(ref)
r=zeros(2,n); 
r(1,:)=resample(ref(1,:),20,1,1000);  
r(2,:)=resample(ref(2,:),20,1,1000);
clear ref;
ref=zeros(2,n); 
ref(1,:)=r(1,:); ref(2,:)=r(2,:);
% Start iterations up to maxiter
maxiter=100; z=zeros(stat,maxiter); amp=z; time=z; resamp=z;
restime=z; cumsn=0; sni=0; it=0; 
%
amax=1e-5; tmax=0.05;
%
for i=1:stat
       disp('Evt');
       i
       % [mns maxsig(i)]=sample(i,tstep,fy,tau,stat);
       maxsig(i)=max(abs(ffy(i,:)));
       mns=ffy(i,:);
       mns(1)=0;
       a=maxsig(i);
       b=0;
       if snn~=inf
           cumsn=cumsn+snn(i);
           sni=sni+1;
       end
     for ii=1:2;
       % Reference delayed by t0, scaled to i0
       [ys,yps]=ref_delayed_mcp(a,b,stat,ref);
       % ys=-ys; yps=-yps;
       scalex=(1:size(ys,2))*ts*tstep;
       % plot(scalex,ys); 
       size(ys)
       mns=mns(1:size(ys,2));
       % mns=mns(20*tau:size(ys,2));
       % hold; plot(scalex,mns,'r'); 
       % figure(1); plot(scalex,ys); figure(2); plot(scalex,yps); figure(3); plot(scalex,mns);
       % return
       [x y]=iter(ys,yps,mns)
       % iter gets as x the fraction mns/ys, y the delay in unit of samples
       resamp(i,ii)=a*x; 
       restime(i,ii)=y*ts;
       amp(i,ii)=resamp(i,ii); a=amp(i,ii);
       time(i,ii)=sum(restime(i,:));
       % Sum to get the first guess and add-up small further corrections
       % b=round(sum(restime(i,:)));
       % disp('sum restime')
       b=sum(restime(i,:));
       it=it+1;
       % pause;
     end
   ii=3;
   % while (abs(amp(i,ii-1)-amp(i,ii-2))>mean(maxsig)*1e-5 | abs(time(i,ii-1)-time(i,ii-2))>=.8);
   while (abs(amp(i,ii-1)-amp(i,ii-2))>mean(maxsig)*amax | abs(time(i,ii-1)/ts-time(i,ii-2)/ts)>=tmax);
       i
       % a
       % b
       [ys,yps]=ref_delayed_mcp(a,b,stat,ref);
       figure(1000); scalex=(1:size(ys,2))*ts*tstep;
       % title('ref_delayed_mcp');
       plot(scalex,ys,'r'); hold; 
       plot(scalex,mns,'b'); title('while');
       pause(0.25);
       % ii
       hold;
       [x y]=iter(ys,yps,mns)
       it=it+1;
       resamp(i,ii)=a*x; restime(i,ii)=y*ts;
       amp(i,ii)=resamp(i,ii);  a=amp(i,ii);
       time(i,ii)=sum(restime(i,:));
       % disp('sum restime')
       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
        abs(amp(i,ii-1)-amp(i,ii-2))
        mean(maxsig)*amax
        abs(time(i,ii-1)-time(i,ii-2))
        figure(4); plot(ys,'b'); hold; plot(mns,'r');
        title('No convergence'); return
      end
   end
   i
   ii
   aa(i)=amp(i,ii-1);
   tt(i)=-time(i,ii-1);
   tt(i)=tt(i)-1;
   % Mystery...
   % disp('timeabs end loop');
   % timeabs(i)=tt(i)-t0(i)+tt(1);
   timeabs(i)=tt(i);
   % timeabs in units of tsteps
   timeabs(1:i);
   ab(i)=(aa(i)-maxsig(i))/maxsig(i);
   % disp('convergence reached')
   time(i,1:maxiter);
   amp(i,:);
   % Plot convergence
   pause(1);
   % sz=size(ys,2); scalex=ts*(1:sz)*tstep; s=ts*tstep;
   % figure(9); plot(scalex,ys,'r'); hold;  plot(scalex,mns,'b'); hold;
   % xlabel('s'); ylabel('Amplitude'); 
   % title('Sampled measurement and reference');
end
disp('average iterations')
it/stat
% Checks 
% disp('i0 residuals')% ab;
% disp('t0 residuals')% timeabs;
figure(10); subplot(1,2,1); hist(ab(2:stat));
subplot(1,2,2); hist(tstep*(timeabs(2:stat)));
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
