% localization example 3: localization with dot detection
%
% this script (i) generates simulated data with a single diffusing
% fluorophore in it, (ii) runs a set of localization analyses on each
% image, (iii) saves the results, and (iv) deletes the image data. 
%
% Requirements: 
% - matlab R2015a or later, with 
%   image_toolbox
%   optimization_toolbox
%   statistics_toolbox
% - uncertainSPT, with all folders and subfolders on the matlab path, e.g., 
%       addpath(genpath(fullfile(pwd,'..','..','uncertainSPT'))); 
% - the ../psfModels/ folder must be in the matlab path, e.g.,
%       addpath(genpath(fullfile(pwd,'..','psfModels')));
% - SMeagol must be installed and on the matlab path. See 
%   http://smeagol.sourceforge.net/
% To verify that the mesoRD simulation works in matlab, add SMeagol to your
% matlab path and run SM_mesoRD_tracking_test()
%
% The active parameters below constitute a minimal test set. Commented out
% is a set of parameters to reproduce the numerical experiments behind figs
% 2,3 in the paper. Other settings are hardcoded, either in the SMeagol
% runinput file runinput0_SMeagol.m, or in the lower level production
% function localization_analysis.m (not that the latter overrides some of
% the parameters in the runinput file). 

% production runs: The number of spots per set of conditions is set by the
% Nspots parameter. Each call to localization_analysis.m (below) creates
% its own temporary working directory, and writes results to its own .mat
% file in EX1results/, in file names that are partly random and do not
% clash. To stack up statistics fast, one can therefore run multiple
% independent matlab processes in parallell.

addpath('FRST')

%% parameters 
% for test runs
r_Nphotons=[100 320] 
r_bg      =[1]  
dt=30*1e-3;
r_tE=[6e-3] 
r_psfIndex=[2];
Nspots=500
EMgain=50;
readRMS=10;

% parameters for our production runs
%%%r_Nphotons=[100 140 200 320 600]
%%%r_bg      =[1 3] % 
%%%r_tE=[1e-3 6e-3 10e-3 22e-3 3e-3 16e-3 ]% 

% psfindex:
%  1 : RW_514nm_NA149
%  2 : RW_639nm_NA14
%  3 : RW_680nm_NA14
%% check likelihood object
cMax=EMgain*(max(r_Nphotons)/4+max(r_bg));
Emax=cMax/EMgain;
try
    LL=load('logL_EX1.mat');
    missingEX1=false;
catch
    missingEX1=true;
end
if( missingEX1 || LL.EMgain ~= EMgain || LL.readRMS ~= readRMS || LL.cMax ~= cMax || LL.Emax ~= Emax )
    disp('constructing new logL object and saving in logL_EX1.mat ...')
    logLobj=EMCCDfit.logL_EMCCD_lookup(EMgain,readRMS,cMax,Emax);
    save logL_EX1.mat logLobj cMax Emax EMgain readRMS
    clear logLobj
    disp('... done!')
else
    disp('Found pre-computed logL_EX1.mat with correct parameters.')
end
clear LL cMax Emax
%% production
rng('shuffle')
for kN=1:numel(r_Nphotons)
    Nphotons=r_Nphotons(kN);
    for kb=1:numel(r_bg)
        bg=r_bg(kb);
        for kE=1:numel(r_tE)
            tE=r_tE(kE);
            for kp=1:numel(r_psfIndex)
                psfIndex=r_psfIndex(kp);
                trj_movie_FRST_localization(Nphotons,bg,dt,tE,psfIndex,Nspots,EMgain,readRMS)
            end
        end
    end
end
