function[fy,snn,sn_mcp,sn_elec,t0,i0,maxsig,tstep]=database(stat);
return
% Create a database of MCP signals randomly distributed in  
% delays and amplitudes
% warning off
% clear fy snn t0 i0 maxsig
% [tstep,trise,tfall,gain,NPE,Qshot,G,nadc,bw,tau,taus,T,R,fs,ts,apjit]=params()
disp('math time step from database')
tstep
% disp('Bandwidth')
bw=.35/(tau*tstep)
%
% Central distribution for MCPs with NPE photoelectrons. 
% Charge signal q0=imax*(trise+tfall)*tstep/2
imax=2*NPE*gain*1.6e-19/((trise+tfall)*tstep);
disp('(tr+tf)/tstep')
(trise+tfall)*tstep
i0=random('normal',imax,imax/sqrt(NPE),1,stat);
%
% Fixed amplitude
% i0=imax*ones(1,stat);
% Non central
% i0=imax*random('nct',10,1000,1,stat);
% hist(random('nct',10,10000,1,20000),50);
%
t0=-round(random('uniform',-.5*tau,.5*tau,1,stat));
%
% First normalized run
t0(1)=0;
i0(1)=imax;
% Check i0
disp('i0, t0 (database)');
i0
t0
%
if 4*tau-t0-2<=0 
    disp('t0>=4*tau-2'); disp(t0); disp(tau); return
if n<4*tau+t0
    disp('n<4*tau+t0');disp(n); disp(t0); disp(tau); return
end
end
    % Calculate signal
% t0max=floor(3*tau);
t0max=floor(1*tau);
n=round(20*tau);
t=1:n;
% sh(1:n)=(2.718.*(t-(t0max+t0))/tau).*exp(-(t-(t0max+t0))/tau);
% sh(1:t0max+t0-1)=0;
% sh=conv(sh,sh);
% sh=sh(1:2:size(sh,2));
for i=1:stat;
    disp('stat')
    i
 % MCP signal 
    [fy(i,:),snn(i),sn_mcp(i),sn_elec(i),maxsig(i)]=mcp1(tau,taus,i0(i),t0(i));
end
%
% Signal references generation
%
i0_0=i0(1);
t0_0=t0(1);
ref=ref_gen(stat,i0(1),t0(1));
% fy=[fy zeros(1,size(fy,2)) zeros(1,size(fy,2))];
sfy=size(fy,2);
fy(i+1,:)=ref(1,1:sfy);
fy(i+2,:)=ref(2,1:sfy);
%
% plot signals and references
%
disp('mean snn_mcp');
mean(sn_mcp)
disp('mean snn_elec');
mean(sn_elec)
disp('mean snn');
mean(snn)
%
% Check reference
% figure(2); subplot(1,2,1); plot(fy(stat+1,:)); 
% subplot(1,2,2); plot(fy(stat+2,:));
% hold off
dataplot(fy);
