function[stdt,mean_snn]=cfd_stat(fy,snn,maxsig,i0,t0,stat);
clear figure;
% tstep=math time step, n=number of math time steps
% tau = peaking time (number of samples), sn=signal/noise ratio, 
% stat=statistics size, nech = number of samples
%
% [G,nadc,tau,elec_noise,tstep]=elec_params();
% [trise,tfall,gain,Ishot,NPE,T]=detector_params(tau,tstep);
%
[tstep,trise,tfall,gain,NPE,Ishot,G,nadc,bw,tau,T,R,fs,ts,apjit]=params();
[thresh,frac,delay]=cfd_params(maxsig)
clear figure(1); clear figure(10);
disp('tstep')
tstep
%
% Check with fixed amplitude i0=1*ones(1,size(fy,2)-2))
% Check with fixed t0; t0=10*ones(1,size(fy,2)-2))
% Stat noisy measurements at sn, amplitude i0(i), shifted by t0(i)
%
x=1; y=0; t0(1)=0;
xx=(1:size(fy,2));
% Real CFD
% figure(8); plot(fy,'r'); hold;  plot(xx,th);   title('Noisy shaper');
% pause;
% Check noisy shaper 
% figure(15); plot(fy,'r'); title('Noisy shaper');    % return %
% pause
scalex=0:tstep:(size(fy,2)-1)*tstep;
delay=round(delay/tstep);
%
for i=1:stat
    disp('Evt');
    i
 % disp('t0'); t0(i)
 % disp('i0'); i0(i)
 %
 % CFD signal
 %
 sfy=size(fy,2); dfy=[zeros(1,delay) fy(i,:)];
 dfy=dfy(1:sfy); fracfy=-frac*fy(i,:); cfdsig=dfy+fracfy;
 % cfd signals
 % clear figure(10);
 % plot(scalex, fy(i,:),'blue'); hold; plot(scalex, dfy,'black'); 
 % plot(scalex, cfdsig,'green'); 
 % title('Signal and delayed signal'); pause
 %
 % Plot signal, cfd signal and cfd threshold
 clear figure(1);
 % if i==1;
 %     plot(scalex,fy(i,:)); hold; plot(scalex, cfdsig,'black');
 %     plot(scalex,thresh,'red');
 %    title ('Signal, cfd signal and threshold'); 
 % end
 %
 % Look for trigger 
 %
 tt=1;
 while fy(i,tt)<thresh
       tt=tt+1;
       if tt>=sfy
           disp('signal under threshold')
           max(fy(i,:))
           thresh
       end
 end
 % disp('time at threshold');
 tt;
 % Look for time
 while cfdsig(tt)<0
       tt=tt+1;
       if tt>=size(cfdsig)
           disp('signal under threshold')
           max(cfdsig)
           thresh
       end
 end
 if i==1
     tt0=tt(i);
 end
 tt(i)=tt;
 % disp('timeabs(i)=')
 timeabs(i)=tt(i)-t0(i)-tt0;
 % Checks
 timeabs;
end
% 
% Remove NaN's
%
p=1;
for i=1:stat
  if isnan(timeabs(i));
    time(i:stat-1)=timeabs(i+1:stat); timeabs=timeabs(1:stat-p);
    p=p+1;
  end
end
cumsn=0; sni=0;
for i=1:stat
       if snn~=inf
           cumsn=cumsn+snn(i);
           sni=sni+1;
       end;
end
% Checks 
disp('t0 residuals')
timeabs;
figure(10); 
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
disp('stdt t0 from cfd_stat')
stdt0=std(t0(2:stat))*tstep
disp('stda i0 from cfd_stat')
stda0=std(i0(2:stat))
disp('stdt from cfd_stat')
stdt=std(tstep*timeabs)
mean_snn=cumsn/(sni);
disp('mean snn from cfd_stat')
mean_snn
