function [Dall,len, other_len, nb315, time, Track_density, peak_density, sd_med, std_90] = ComputeDiffConstants3(expName,minT,dminT,dt, Res_max)
% First section of script will create the file with all diffusion
% coefficients of single trajectories from one data set.
% Diffusion constants are calculated based on the first 4 step interval
% (defined by the variable 'limit')
% 
%expName = strcat(PathC_name,expName)

steps_5 = 1; % takes only the 5 first steps of each track if steps_5 = 1
double_dt = 0 % take only the odd steps to multiply the Delta t by 2
min_dt = 1; max_dt = 4;

L = 5; 
maxL = 12;
limit = 4; % time points for linear fit to data

folder = strcat(string(cd), filesep, 'Tracks', filesep, expName, filesep);
fileList = dir(strcat(folder,'*pathsC.mat'));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fileList = dir('run/user/1000/gvfs/smb-share:server=gaia.pasteur.fr,share=%40mcm/Francois')01-All/Tracks/Res0.15_Intpercentage_std3_MaxDisp4_rawth1.65');

NB = 0; NBimmo = 0;
Xdispls=[];   Ydispls = []; len = [];   eeDists = []; 
Dall = [];    ball = [];    immo = [];     NBdispls = []; 
R2all =[];    immoTraj=[];  DDisplVec= []; AllPaths = []; other_len = [];
nb315 = 0; time =[]; Track_density = 0; movie_NB = 0; sample_peak_density = 0;

%%
for i = 1:length(fileList) % loop over all path files
    clear paths fitP BW peaknb
    
    load(strcat(folder,fileList(i).name),'paths', 'fitP', 'BW', 'peaknb');
    
    sd_med = median(fitP(:,3)); std_90 = quantile(fitP(:,3), 0.9)
    
    if exist('peaknb') ~= 1
        peaknb = 0
    end
        
    paths = [paths, fitP];
    
    if(isfinite(minT)&& sum(isfinite(paths(:))))
        paths = paths(paths(:,3)>= minT & paths(:,3)< dminT+minT, :);
    end
    %figure, histogram(paths(:,9)./paths(:,7))
    if(isempty(paths))
        pathnb=[];
    else
        pathnb = unique(paths(:, 4));
    end
    
    for j = pathnb' % loop over all paths in a file
        
        path = paths(paths(:, 4) == j, :);
        
        
        leng = size(path, 1);
        if double_dt == 1
           path = path(1:2:end,:);
           leng = size(path, 1);
        end
        if steps_5 == 1 
            path = path(1:min(end,5),:);
        end
        avg_Res = mean(path(:,9)./path(:,7));
        if (leng >= 3) && (leng <= 15), nb315 = nb315 + 1; end
        if (leng >= L) && (leng <= maxL) && (avg_Res < Res_max)
            NB = NB + 1;
            len = [len leng];
            time = [time path(1,3)];
            eeDists = [eeDists sqrt(sum((path(1, 1:2) - path(end, 1:2)).^2))*75];
            clear xdispl ydispl
%             for k = 1:leng-1
%                 xdispl(k) = (path(k, 1) - path(k+1, 1))*75;
%                 ydispl(k) = (path(k, 2) - path(k+1, 2))*75;
%             end
%             Xdispls = [Xdispls xdispl];
%             Ydispls = [Ydispls ydispl];
%            NBdispls = [NBdispls ones(size(ydispl))*NB];
            
            % Calculate MSD and fit linear curve to it
%             msd = calcMSD(path)*75^2;
%             x = (1:length(msd))*(TStep*0.001);
%             fo = fitoptions('Method','NonlinearLeastSquares','Robust','Off', 'Algorithm', ...
%                 'Trust-Region', 'Lower', [-Inf -Inf], 'Upper',[Inf Inf],'StartPoint',[0.0016*4 0], ...
%                 'MaxFunEvals', 600, 'MaxIter', 400);
%             
%             ft = fittype('m*x + b','options',fo);
%             [curve,gof] = fit(x(1:limit)', msd(1:limit)'/(1000^2),ft);

            curve.m = calcMSD2(path, L, min_dt, max_dt, dt);
            curve.b = 0;
            
            Dall = [Dall curve.m/4];
            ball = [ball curve.b];
            %R2all = [R2all gof.rsquare];
            %DDisplVec = [DDisplVec ones(size(ydispl))*curve.m/4];
            
%             AllPaths = [AllPaths; [path ones(size(path, 1), 1)*NB...
%                 ones(size(path, 1), 1)*curve.m/4 ones(size(path, 1), 1)*curve.b ones(size(path, 1), 1)*gof.rsquare]];
%             
%             if curve.m/4 < 0.01%Inf%
%                 immo = [immo ones(size(ydispl))];
%                 immoTraj = [immoTraj 1];
%                 NBimmo = NBimmo + 1;
%             else
%                 immo = [immo zeros(size(ydispl))];
%                 immoTraj = [immoTraj 0];
%             end
        else
            other_len = [other_len, leng];
        end
    
    end
    movie_NB = NB - movie_NB;
    movie_Track_density = movie_NB/mean2(BW)
    sample_peak_density = sample_peak_density + peaknb/mean2(BW)

    Track_density = Track_density + movie_Track_density; % number of tracks if the field of view was full
end

Track_density = round(Track_density/length(fileList))
peak_density = round(sample_peak_density/length(fileList))
