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
% timeabs in absolute time units
% Get params
[tstep,trise,tfall,gain,NPE,Ishot,G,nadc,bw,tau,taus,T,R,fs,ts,apjit]=params();
%
ts
% 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
r(1,:)=resample(ref(1,:),round((1/fs)/tstep),1,10);  
r(2,:)=resample(ref(2,:),round((1/fs)/tstep),1,10);
clear ref;
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; amax=1e-4; tmax=1;
restime=z; cumsn=0; sni=0; it=0; 
%
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);
       mns=mns(1:size(ys,2));
       % mns=mns(20*tau:size(ys,2));
       if ii==1
       % figure(15);
       % plot(scalex,mns,'b');
       % hold;  
       % plot(scalex,ys,'r'); % figure(2); plot(scalex,yps); 
       % hold;
       % plot(scalex,mns);
       % pause
       % hold
       end
       [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;
       %
       % Modif 3/27/09
       %
       % restime=y;
       amp(i,ii)=resamp(i,ii); a=amp(i,ii);
       time(i,ii)=sum(restime(i,:));
       b=sum(restime(i,:));
       it=it+1;
     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,'b'); hold; 
       % plot(scalex,mns,'r'); title('while');
       ii
       [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
   aa(i)=amp(i,ii-1);
   tt(i)=-time(i,ii-1);
   disp('tt');
   tt(i)=tt(i)-1;
   tt(i)=tt(i)*tstep
   clear figure(15)
   % pause
   % -1=mystery...
   % disp('timeabs end loop');
   % timeabs(i)=tt(i)-t0(i)+tt(1);
   timeabs(i)=tt(i);
   % timeabs in absolute units ....  not tsampling !
   timeabs(1:i);
   ab(i)=(aa(i)-maxsig(i))/maxsig(i);
   % disp('convergence reached')
   time(i,1:maxiter);
   amp(i,:);
   % Plot convergence
   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(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(timeabs(2:stat))
stda0=std(i0(2:stat))/mean(maxsig)
stda=std(ab(2:stat))
disp('mean_snn');
mean_snn=cumsn/sni
return
