function trj_movie_localization(Nphotons,bg,dt,tE,psfIndex,Nspots,EMgain,readRMS)
% trj_movie_localization(Nphotons,bg,dt,tE,psfIndex,Nspots,EMgain,readRMS)
% Run of localization/precision analysis on simulated data, using a single
% set of parameters:
% Nphotons  : number of photons/spot
% bg        : uniform background intensity, photons/pixel
% dt        : time step
% tE        : exposure time
% psfIndex  : choice of simulated point-spread function, from the list 
%  1 : RW_514nm_NA149 (514 nm, NA=1.49, Richards&Wolf model)
%  2 : RW_639nm_NA14  (639 nm, NA=1.4,  Richards&Wolf model)
%  3 : RW_680nm_NA14  (680 nm, NA=1.4,  Richards&Wolf model)
% these are parameterized in the ../psfModel/ folder, and 
% Nspots    : number of spots to simulate (number of frames in the simulated movie) 
% EMgain    : EMCCD gain parameter for the simulated camera
% readRMS   : std. of camera readout noise (Gaussian)
%
% other simulation parameters are either hardcoded below, or in the SMeagol
% runinput file runinput0_SMeagol.m

%% hardcoded parameters
if(~exist('EX1results','dir'))
   !mkdir EX1results 
end

% localization parameters
ROIwidth=9;
Nquad=5;
saveErrors=false; 

% temporary write folder
d0=pwd;
tmpNum=ceil(9e7*rand);
tmpDir=fullfile(d0,['tmp_' int2str(tmpNum) ]);
try
    system(['rm -rf ' tmpDir])
catch
end
system(['mkdir  ' tmpDir])

% output file name & PSF model
psfFile_rel={...
'../psfModels/psf_RW_514nm_NA149_ni1515.mat',...
'../psfModels/psf_RW_639nm_NA14_ni1515.mat',...
'../psfModels/psf_RW_680nm_NA14_ni1515.mat'};
psfString={...
'RW_514nm_NA149',...
'RW_639nm_NA140',...
'RW_680nm_NA140'};

lambda=[514  639  680 ];
NA    =[1.49 1.4  1.4 ];
PSFfile=fullfile(pwd,psfFile_rel{psfIndex});
lambda=lambda(psfIndex);
NA=NA(psfIndex);

resFile=fullfile(pwd,'EX1results',sprintf('resultEX1_bg%03.1f_N%04d_%s_dt%05dms_tE%05dms_%08d.mat',...
    bg,Nphotons,psfString{psfIndex},dt/1e-3,tE/1e-3,tmpNum));
%% run mesord 
% mesord parameters
modelPath =fullfile(pwd,'..','mesoRDmodel_D100_longCells','Ecoli_1state_D100.xml');
mesordOptions=['  ' modelPath ' -K -g -x 0 -c 1 -t ' num2str(Nspots*dt+0.2) ' -i 1 -I 50 -w 0.0037 -q 0.020 um -c 5 '];
SM_runMesoRD(mesordOptions,tmpDir)
cd(d0)
tmpDir_relative=SM_relative_path(d0,tmpDir);
%% run SMeagol simulations
SMopt=SM_getOptions('runinput0_SMeagol.m');

SMopt.background.photons_per_pixel=bg;

SMopt.sample.dt=dt;
SMopt.sample.tE=tE;
SMopt.baseIntensity.intensity=Nphotons/tE;
SMopt.camera.alpha=1/EMgain;
SMopt.camera.sigmaReadout=readRMS;
SMopt.psf.type = 'SM_psf_numericalPSF_Ruzp';
SMopt.psf.PSFfile=PSFfile;

SMopt.trj.trajectoryFile=fullfile(tmpDir_relative,'trajectories.txt');
SMopt.trj.reactionFile  =fullfile(tmpDir_relative,'reactions.txt');

SMopt.output.resultFile=fullfile(tmpDir_relative,'SMeagol_result.mat');
SMopt.output.maxFrames=Nspots;

SMopt=SM_runsimulation(SMopt);
R=SM_getResults(SMopt);
%% extract true coordinates
emTrj_nm=R.emissionAverage{1}; % x y z frame cell  
nm2px=1/R.opt.camera.pixLength;
emTrj_px=[emTrj_nm(:,1:3)*nm2px emTrj_nm(:,4:end)]; % the initial 'guess'
T=size(emTrj_px,1);
clear emTrj_nm
%% load MLEobject
load logL_EX1.mat
%% set up localizations
fluoOffset=SMopt.camera.offset*ones(SMopt.camera.yrange_px,SMopt.camera.xrange_px);
MV=EMCCDfit.ML_loadStack2(R.movieName{1});
sigMin=0.21*lambda/NA/SMopt.camera.pixLength;

% set up symmetric Gaussian PSF objects and coordinate lists
PSFsymGauss=cell(0);
% straight MLE
PSFsymGauss{1}=PSF.SymGauss_MLE(...
    'initialGuess',[0 0 log([bg Nphotons 1.1*sigMin])]);   

% hard lower limit on S, ds=S/S0-1~exp(1), lnB~N(0.ln30)
PSFsymGauss{end+1}=PSF.SymGaussS0_logNormBN_expS0('lambda',lambda*nm2px,'NA',NA,...
    'initialGuess',[0 0 log([bg Nphotons 1.1*sigMin])],...
    'priorParameters',[log([bg 30 Nphotons inf]) 1]);

for np=1:numel(PSFsymGauss)
    SG_coord{np}=zeros(T,2);
    SG_bNSdr{np}=zeros(T,4);
    SG_cov{np}  =zeros(T,3);
    SG_xyErr{np}=zeros(T,2);
    SG_c{np}    =zeros(T,1);
    SG_CRLB{np} =zeros(T,1);
end
% asymmetric Gauss MLE
PSFasymGauss=cell(0);
% asymmetric Gauss MLE
PSFasymGauss=cell(0);
PSFasymGauss{1}=PSF.AsymGauss_angle_MLE(...
    'initialGuess',[0 0 log([bg Nphotons 1.1*sigMin*[1 1]]) 0]);
% hard limit on S1,S2, lnB~N(0,ln30), dS1,dS2~exp(0.5)
PSFasymGauss{end+1}=PSF.AsymGaussS0_BNlnN_expS0('lambda',lambda*nm2px,'NA',NA,...
          'initialGuess',[0 0 log([bg Nphotons 1.1*sigMin*[1 1]]) 0],...
          'priorParameters',[log([bg 30 Nphotons inf]) 0.5 ]);

for np=1:numel(PSFasymGauss)
    AG_coord{np}=zeros(T,2);
    AG_bNS12vdr{np}=zeros(T,6);
    AG_cov{np}  =zeros(T,3);
    AG_xyErr{np}=zeros(T,2);
    AG_c{np}    =zeros(T,1);
    %AG_LAB{np}  =zeros(T,3);
end
%% run localizations
fitTime=zeros(2,max(numel(PSFsymGauss),numel(PSFasymGauss)));
fprintf('\n\n')
fprintf('%05d of %05d ',0,size(MV,3));
for frame=1:size(MV,3)
    ind=find(emTrj_px(:,4)==frame);
    if(~isempty(ind));
        currentFluoFrame=double(MV(:,:,frame));
        frameCoord0=emTrj_px(ind,:);
        
        % symmetric Gaussian with fixed initial guesses
        for np=1:numel(PSFsymGauss)
            t0=tic;
            
            [frameCoord,frameCov,framePar,~,dotPar]=EMCCDfit.refineSingleFrame(...
                frameCoord0,currentFluoFrame,fluoOffset,ROIwidth,Nquad,logLobj,PSFsymGauss{np},[],saveErrors);
            
            fitTime(1,np)=fitTime(1,np)+toc(t0);
            SG_coord{np}(ind,:)=frameCoord(:,1:2);
            SG_bNSdr{np}(ind,:)= [dotPar [framePar(:).dr]'];
            SG_cov{np}(ind,:)   =frameCov;
            SG_xyErr{np}(ind,:) =frameCoord(:,1:2)-frameCoord0(:,1:2);
        end
        
        % asymmetric Gauss with fixed initial Guess
        for np=1:numel(PSFasymGauss)
            t0=tic;
            
            [frameCoord,frameCov,framePar,~,dotPar]=EMCCDfit.refineSingleFrame(...
                frameCoord0,currentFluoFrame,fluoOffset,ROIwidth,Nquad,logLobj,PSFasymGauss{np},[],saveErrors);
            
            fitTime(2,np)=fitTime(2,np)+toc(t0);
            AG_coord{np}(ind,:)=frameCoord(:,1:2);
            AG_bNS12vdr{np}(ind,:)= [dotPar [framePar(:).dr]'];
            AG_cov{np}(ind,:)   =frameCov;
            AG_xyErr{np}(ind,:) =frameCoord(:,1:2)-frameCoord0(:,1:2);
        end
        
        
        fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b')
        fprintf('%05d of %05d ',frame,size(MV,3));
    end
end
fprintf('\n')
%% compute CRLBs
a=1; % pixel width, in pixel units
for np=1:numel(PSFsymGauss)
    s2a=SG_bNSdr{np}(:,3).^2+a^2/12;
    N  =SG_bNSdr{np}(:,2);
    b  =SG_bNSdr{np}(:,1);
    tau =2*pi.*b.*s2a./N/a^2;
    F0=1+4*tau+sqrt(2*tau./(1+4*tau)); % approximation from Rieger&Stallinga
        
    SG_CRLB{np} = 2*s2a./N.*F0;
    SG_c{np}=isfinite(SG_coord{np}(:,1)) & SG_bNSdr{np}(:,4)<ROIwidth/2; % basic convergence indicator
end
%% save results
clear MV currentFluoFrame fluoOffset R s2a N b tau F0  X 
disp(['Saving results to ' resFile ])
save(resFile)
%% delete movies etc
system(['rm ' tmpDir filesep '*'])
system(['rm -rf ' tmpDir])
