function [Rout,BWmax,MSRout]=optimalBW05(Rstin,Tin)

%Tbin is the window to run the optimal bin on, such as Tbin=start:step:end i.e. Tbin=-1:0.005:2
%Tin the PETH duration
%T0 is window w/ the smallest binsize
%
global Tbin

BW=[0.01:0.005:1];
NUMTRL=length(Rstin);PTH0=[];

Bsize=Tbin(2)-Tbin(1);

Rst_all=sort(cell2mat(Rstin'));
T0=[Tbin(1)-Bsize:Bsize:Tbin(end)+Bsize];L0=length(Tbin)-2;
R0=hist(Rst_all,T0)/(Bsize*NUMTRL);R0=R0(2:end-1);SST=sum((R0-mean(R0)).^2);

if isempty(Rst_all)
    R0=R0';
end

AIC=[];
for i=1:length(BW)
    T=[Tbin(1)-BW(i):BW(i):Tbin(end)+BW(i)];P=length(T)-2;PTH=[];
    PTH_all=hist(Rst_all,T)/(BW(i)*NUMTRL);
    lambda_est=interp1(T(2:end-1),PTH_all(2:end-1),T0(2:end-1),'nearest','extrap');
    
    SSE=sum((lambda_est-R0).^2);
    %AIC 
    AIC(i)=L0*log(SSE/L0)+2*P;
      % AIC(i)=2*mean(lambda_est)/(NUMTRL*BW(i))-var(lambda_est)
end

AICF=filtfilt(0.1*ones(1,10),1,AIC);
[qwe,Imin]=min(AICF);BWmax=BW(Imin);
% figure(1);plot(BW,AIC);hold on;plot(BW,AICF,'r')
% fprintf(1,'final BW %1.4f\n',BWmax);

MSRout=AIC;
Tn=[-BWmax/2:-BWmax:Tin(1)-BWmax]; Tp=[BWmax/2:BWmax:Tin(end)+BWmax];T=[Tn(end:-1:1),Tp];%making PSTH centered on zero
%T=[T0(1)-BWmax:BWmax:T0(end)+BWmax];
PTH_opt=hist(Rst_all,T)/(BWmax*NUMTRL);

%Assigning the output
Rout=interp1(T(2:end-1),PTH_opt(2:end-1),Tin,'nearest','extrap');