function[fy,sn,sn_mcp,sn_elec,signal]=mcp(tau,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)
% sn=signal/noise ratio; i0=pulse amplitude, t0=pulse delay 
%
[tstep,trise,tfall,gain,NPE,Qshot,G,nadc,bw,tau,T,R,fs,ts,apjit]=params();
%
% new
% t0max=floor(3*tau);
t0max=floor(1*tau);
% n=round(50*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
% figure(1); plot(y,'r'); title('Detector input'); return
%
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);
% Scale to speedup
% sh=sh(1:2:size(sh,2));
% plot(sh); return
%
% figure(3);  sz=size(sh,2); scalex=(1:sz)*tstep;
% plot(scalex,sh); title('Shaper impulse response'); return
% Shot noise current from MCP. 
%
% Current noise in 3.75 kOhms = voltage noise
% 
fs=1/(ts*tstep);
bw=.35/(tau*tstep); 
Ishot=Qshot*gain*1.6e-19/(3*tau*tstep)
shotnoise=5e4*sqrt(2*1.6e-19*Ishot*(.35*fs/2))*randn(1,n);
%                                                        % Shot noise
% tau/tstep=physical time
% shotnoise=1000*sqrt(2*1.6e-19*gain*Ishot*(2/tstep))*randn(1,n);
% figure(4);  sz=size(shotnoise,2); scalex=(1:sz)*tstep;
% plot(scalex,shotnoise); title('Shotnoise'); return
% disp('Input Shot noise'); std(shotnoise)    
% disp('signal, shot noise  signal/shot noise before convolutions');
%
insig=y+shotnoise;
%
% Check input noisy signal; 
% figure(5); sz=size(y,2); scalex=(1:sz)*tstep;
% subplot(4,1,1); plot(scalex,y); subplot(4,1,2); plot(scalex,shotnoise);
% subplot(4,1,3); plot(scalex,insig); return
%                   
%                                       Conv MCP signal fy
% MCP Signal
fy=conv(insig,sh);
fy=fy(t0max+1:n);
% figure(777)
% plot(fy)
%
% Anti-alias filter
%
nfy=size(fy,2); 
h=hamming(nfy);
fy=fy.*h';
qf=1/(nfy*tstep); nfs=floor(fs/qf);
%
% Shannon window
%
w=zeros(1,nfy); % w(1:floor(nfs/2)-10)=1; w=w+fliplr(w); 
wslope=1; 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));
%
% Conv noiseless MCP signal
%                                       (for signal to noise)
% MCP signal without noise 
iy=conv(y,sh);                         
iy=iy(t0max+1:n);
s=max(abs(iy));
%
% Anti-aliased signal
%
niy=size(iy,2); qf=1/(nfy*tstep); nfs=floor(fs/qf); 
fe=fft(iy); p=fe.*w; 
iy=real(ifft(p));
%
% Normalize signal and noiseless signal to i0 whatever the convolution
% factor
iy=iy*i0/s;
fy=fy*i0/s;
% figure(31); plot(fy); plot(fy);
%
% npy=size(fy,2);
%
% Check shot noisy signal; 
% figure(7);
% sz=size(fy,2); scalex=(1:sz)*tstep;
% plot(scalex,fy); title('Shot noisy signal'); 
%
% Add electronics noise (Thermal up to fs/2)          %  Elec_noise
%
elec_noise=sqrt(4*1.38e-23*300*R*(fs/2));
e_noise=elec_noise*randn(1,n-t0max);
e_noise=e_noise(1:nfy);
%
% Anti-aliased elec noise
%
% nfy=size(e_noise,2); 
qf=1/(nfy*tstep); nfs=floor(fs/qf);
% w=zeros(1,nfy); w(1:floor(nfs/2)-10)=1; w=w+fliplr(w);
fe=fft(e_noise); p=fe.*w; 
e_noise=real(ifft(p));
%
% figure(1); plot(e_noise); figure(2); plot(abs(fft(e_noise)));
%
% Normalize to elec_noise whatever the convolution
e_noise=elec_noise*e_noise/std(e_noise);
%
disp('elec noise');
elec_noise=std(e_noise)
e_noise=e_noise(1:nfy);
%
% e_noise=zeros(1,npy);                               % No elec_noise
%
% figure(8); subplot(2,1,1); plot(fy(1,:)); 
% subplot(2,1,2); plot(e_noise);
% Signal + electronics noise
fy=fy+e_noise;
%
% Check electronics noisy signal; 
% figure(9); sz=size(fy,2); scalex=(1:sz)*tstep;
% plot(scalex,fy); title('Electronics noisy signal'); return
%  
% Processed shot noise                  Conv MCP noise        
%
nsy=conv(shotnoise,sh);
nsy=nsy(t0max+1:n);
%
% Anti-aliased shot-noise
%
nnsy=size(nsy,2); qf=1/(nfy*tstep); nfs=floor(fs/qf);
% w=zeros(1,nnsy); w(1:floor(nfs/2)-10)=1; w=w+fliplr(w);
fe=fft(nsy); p=fe.*w; 
nsy=real(ifft(p));
%
% Same scale factor as fy and iy (78 79)
nsy=nsy*i0/s;
% Convolution tstep factor;
nsy=nsy/(sqrt(fstep)/2);
%
% figure(77)
% plot(nsy);
disp('mcp processed noise');
mcp_noise=std(nsy)
% nsy=zeros(1,npy);                                   % No MCP noise
%
% Total noise
ny=nsy+e_noise;
%
disp('Total output noise')
std(ny)
%
% Spectra
%
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(1:nsh),fiy(1:nsh),'green');
hold;
% disp('MCP shot noise')
nsy=nsy(1:nfy);
hnsy=h.*nsy';
fnsy=abs(fft(hnsy)); % fnsy=fnsy(2:nf);
% loglog(scalef(2:nf),fnsy(2:nf),'red');
outfnsy=fnsy(1:nsh);
loglog(scalef(1:nsh),fnsy(1:nsh),'red');
% 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');
hold off;
%
signal=i0;
stdnoise=std(ny);
% Output Noise spectrum
% fny=abs(fft(ny.*hamming(1,size(ny,2))));
% df=1/(size(fny,2)*tstep);
% scalef=(2*df:df:floor(1/(2*tstep)));
% figure(10); semilogx(scalef,fny(1:size(scalef,2)));
%
% Check processed noise
% figure(11);
% sz=size(ny,2); scalex=(1:sz)*tstep;
% plot(scalex,ny); title('Total noise'); return
% figure(12); plot(iy(1,:))
%
% 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'); 
