function[fy,sn,sn_mcp,sn_elec,signal,outfiy,outfnsy,outfe_noise,outffy]=mcp1(tau,taus,i0,t0);
% Calculates MCP detector + electronics output signal
% with amplitude i0 and time t0. 
% tstep=math time step;   n=number of time steps
% tau=oscilloscope time constant (number of samples)
%
[tstep,trise,tfall,gain,NPE,Qshot,G,nadc,bw,tau,taus,T,R,fs,ts,apjit]=params();
%
% tau/tstep=physical time
% t0max=floor(3*tau);
t0max=floor(1*tau);
% n=round(20*tau);
n=round(10*tau);
% tstep scaling factor for convolutions
fstep=tstep/500e-15;
clear x y; 
clear figure;
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
%
[y]=wave(tstep,trise,tfall,t0max,i0);
% plot(y); return
% disp('i0 from mcp'); i0
%
% Check detector signal
% Total input charge (e-)
% disp('q0 in electrons (mcp2)'); q0=tstep*sum(y)/1.6e-19
%
t=1:n;
sh(1:n)=(2.718.*(t-(t0max+t0))/tau).*exp(-(t-(t0max+t0))/tau);
sh(1:t0max+t0-1)=0;
% plot(sh); return
sh=conv(sh,sh); sh=sh(1:2:size(sh,2)); sh=sh/max(sh);
%
% figure(3);  sz=size(sh,2); scalex=(1:sz)*tstep;
% plot(scalex,sh); title('Shaper impulse response'); return
% Shot noise current from MCP. 
% 
fs=1/(ts*tstep); bw=.35/(tau*tstep); 
Ishot=Qshot*gain*1.6e-19/(3*tau*tstep)
%
%  Shot noise up to 1/2*tstep
%% Maximum shot noise  
% Current shot noise in 10 kOhms (arbitrary !) = voltage noise
%
shotnoise=10e3*sqrt(2*1.6e-19*Ishot*(1/(2*tstep)))*randn(1,n);
%
% Low pass shotnoise at taus=.35/bws
t=1:n;
shs(1:n)=(2.718.*(t-(t0max+t0))/taus).*exp(-(t-(t0max+t0))/taus);
shs(1:t0max+t0-1)=0;
shs=conv(shs,shs); shs=shs(1:2:size(shs,2)); shs=shs/max(shs);
% Process shotnoise
shotnoise=conv(shs,shotnoise); shotnoise=shotnoise(t0max+1:n);
%                   % Shot noise
% figure(4);  sz=size(shotnoise,2); scalex=(1:sz)*tstep;
% plot(scalex,shotnoise); title('Processed Shot noise');
% disp('Input Shot noise'); std(shotnoise)
%
% shotnoise=1e-10*zeros(1,size(shotnoise));                  % No shot noise    
%
% Check input; 
% figure(5); sz=size(y,2); scalex=(1:sz)*tstep; plot(scalex,y); return           
%
% MCP Signal
fy=conv(y,sh); fy=fy(t0max+1:n); nfy=size(fy,2);                          
fy=fy*i0/max(abs(fy));                                  % Normalize fy to i0
%                                                   
% Check signal; 
sz=size(fy,2); scalex=(1:sz)*tstep;
% figure(6);
% plot(scalex,fy); title('MCP Signal');
%
% New Shotnoise 2q Id: sqrt((Q2/T)*(1/T)) = Q/T = I
shotnoise=shotnoise.*sqrt(fy);
%
% figure(7); plot(scalex,shotnoise); title('Shotnoise');
%
fy=fy+shotnoise;
%
% figure(8); plot(scalex,fy); Title('Noisy signal');
%
% Add electronics noise (Thermal up to 1/(2*tstep))        Elec_noise
%
elec_noise=sqrt(4*1.38e-23*300*R*(1/(2*tstep)));
%
% Elec_noise does not depend upon tau before processing, R only
e_noise=randn(1,sz);
% 
% Normalize to elec_noise 
e_noise=elec_noise*e_noise/std(e_noise);
%
% No elec_noise
% e_noise=1e-10*zeros(1,sz);    
% No shot noise
% shotnoise=1e-20*zeros(1,sz);
%
% figure(8); subplot(2,1,1); plot(fy(1,:)); 
% subplot(2,1,2); plot(e_noise);
%
% Signal + shotnoise + electronics noise 
%
fy=fy+e_noise;
%
% Check noisy signal; 
% figure(9); scalex=(1:sz)*tstep;
% plot(scalex,fy); title('Noisy signal');                 return
%
% Noiseless MCP signal     for signal/noise (useless)
%
iy=conv(y,sh);                         
iy=iy(t0max+1:n);
%
s=max(abs(iy))
% Normalize noiseless signal to i0 
iy=iy*i0/s;
%
%  Anti-alias everyone
%
% Smoothed Shannon window
nfy=sz;
qf=1/(nfy*tstep); nfs=floor(fs/qf); w=zeros(1,nfy);  
wslope=5; w(1:floor(nfs/2)-wslope)=1; 
w(floor(nfs/2)-(wslope-1):floor(nfs/2))=1-[0:1/wslope:1-(1/wslope)];
w=w+fliplr(w);
%
% Anti-aliased signal+noise
fe=fft(fy); p=fe.*w; 
fy=real(ifft(p));
%
% Anti-aliased signal
fe=fft(iy); p=fe.*w; 
iy=real(ifft(p));
%
% Anti-aliased elec noise
fe=fft(e_noise); p=fe.*w; 
e_noise=real(ifft(p));
%
% Anti-aliased shot-noise
fe=fft(shotnoise); p=fe.*w; 
nsy=real(ifft(p));
%
% Spectra
% More frequency points
%
% z=zeros(1,size(fy,2)); iy=[z iy z];
% fy=[z fy(1,:) z]; nsy=[z nsy z];
% e_noise=[z e_noise z]; nfy=size(fy,2);
%
figure(10);
h=hamming(nfy);
nf=round(nfy/2)
nsh=(5e8+fs/2)*tstep*nfy;
scalef=(1:nf)*(1/(tstep*nfy));
% disp('MCP without noise')
hiy=h.*iy(1:nfy)';
fiy=abs(fft(hiy));
% loglog(scalef(2:nf),fiy(2:nf),'green');
outfiy=fiy(1:nsh);
% loglog(scalef(2:nsh),fiy(2:nsh),'green');
% disp('MCP shot noise')
nsy=nsy(1:nfy);
hnsy=h.*nsy';
fnsy=abs(fft(hnsy));
qqf=1/(nfy*tstep)
loglog(abs(fnsy(1:nfy/2)))
% loglog(scalef(2:nf),fnsy(2:nf),'red');
outfnsy=fnsy(1:nsh);
loglog(scalef(1:nsh),fnsy(1:nsh),'red');
hold;
% disp('Elec noise')
he_noise=h.*e_noise';
fe_noise=abs(fft(he_noise));fe_noise=fe_noise(1:nf);
% loglog(scalef(2:nf),fe_noise(2:nf),'blue');
outfe_noise=fe_noise(1:nsh);
loglog(scalef(1:nsh),fe_noise(1:nsh),'blue');
% disp('MCP signal spectrum ');
hfy=h.*fy';
ffy=abs(fft(hfy)); % ffy=ffy(2:nf);
% loglog(scalef(2:nf),ffy(2:nf),'black');
outffy=ffy(1:nsh);
loglog(scalef(1:nsh),ffy(1:nsh),'black');
loglog(scalef(1:nsh),fiy(1:nsh),'green');
hold off;
%
signal=i0;
%
% Noises stds
%
disp('mcp processed noise');
% figure(200);
% plot(1:3000,fy(1:3000),round(2.5*tau+t0),0,'r+',round(3.5*tau+t0),0,'r+'); 
mcp_noise=std(shotnoise(2.5*tau+t0:3.5*tau+t0))
% mcp_noise=std(shotnoise)
%
disp('elec noise');
elec_noise=std(e_noise)
%
% Total noise
ny=shotnoise+e_noise;
%
disp('Total output noise')
stdnoise=std(ny)
%
%
% Actual signal/noise
disp('signal to noises');
sn_mcp=signal/mcp_noise
sn_elec=signal/elec_noise
sn=signal/stdnoise 
% sz=size(fy,2); scalex=(1:sz)*tstep;
% figure(20); plot(scalex,fy); title('Last signal from mcp'); 
