function [Dall,len, other_len, nb315, time, Track_density, peak_density, sd_med, std_90, dir_R, mean_motion, fitp1, fitp2, spot_perc, RSSs, motions] = 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 = 0; % 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 = 3;

L = 4; 
maxL = Inf;
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; RSSs = []
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;
dir_R = []; mean_motion = []; fitp2 = []; fitp1 = []; spot_percs = []; motions = []
%%

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');
    
    track_pos = []
    
    %movie_name = fileList(i).name; movie_name = movie_name(max(end - 25,1):end - 11);
    %img_path = regexp(folder,'/','split');
    %img_path = strcat('/', img_path(2), '/', img_path(3), '/',img_path(4), '/',img_path(5), '/',img_path(6), '/',img_path(7), '/',img_path(8), '/',img_path(11), '/', movie_name,'/Pos0/');
    %img_dir = dir(strcat(img_path, '/img*'))
    
    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,L),:);
        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)

        len = [len leng];
        
%             idx = 0
%             for kk = path(:,3)'
%                 idx = idx + 1
%                 img = double(imread(strcat(img_path, img_dir(kk).name)));
%                 figure(15), hold on, imagesc(img)
%                 plot(path(:,2), path(:,1), '-or')
%                 x = path(idx,1); y = path(idx,2)
%                 scatter(y, x, 'og')
%                 xlim([path(1,2)-6,path(1,2)+6])
%                 ylim([path(1,1)-6,path(1,1)+6])
%                 path(idx,:)
%                 path(idx,9)/(path(idx,14)-path(idx,12))
%                 
%                 mean2(img(floor(x):ceil(x),floor(y):ceil(y)))
%                 
%                 caxis([1000 6000])
%                 colormap jet
%                 pause
%             end            
            
%%%

            NB = NB + 1;
            %len = [len leng];
            time = [time path(1,3)];
            eeDists = [eeDists sqrt(sum((path(1, 1:2) - path(end, 1:2)).^2))*80];
            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, msd] = calcMSD2(path, L, 1, 4, dt);
%             curve.m = mean(sum((path(5:end, 1:2)*0.08-path(1:end-4, 1:2)*0.08).^2, 2));
%             curve.m = mean(sum((path(4:end, 1:2)*0.08-path(1:end-3, 1:2)*0.08).^2, 2));
% 
%             % msd = [curvefit([0.2, 0.3], dts)]
%             
%             tvec = dt:dt:5*dt;
%             dts = [1:4];
%             curvefit = @(p, dts) p(1)*dts.^2 + p(2);
%             
%             p0 = [0.02, 0.02];
%             opts = optimset('Display','off');
%             [a ,RSS,residuals,exitflag,output] = lsqcurvefit(curvefit, p0, dts, msd, [],[],opts);
%             RSS = RSS^0.5;
%             
%             
%             
%             corrcoef(curvefit(a, dts), msd)
% 
%             r2 = corrcoef(curvefit(a, dts), msd)
%             r2 = r2(2,1)
%             
%             
%             
%             motions = [motions, a(1)*80 / 3.5];
%             RSSs = [RSSs, r2];
%             
% 
%             curve.b = 0;
%             %figure(3), plot(path(:,1), path(:,2))
%             
%             R = corrcoef(path(:,1:2));
%             dir_R = [dir_R, R(1,2)]; dir_R(end);
%             
            
            %mean_motion = [mean_motion, mean(sum((path(3:end, 1:2)-path(1:end-2, 1:2)).^2, 2))*0.080 - mean(sum((path(2:end, 1:2)-path(1:end-1, 1:2)).^2, 2))*0.080];
            mean_motion = [mean_motion; mean(sum((path(2:end, 1:2)-path(1:end-1, 1:2)).^2, 2)), mean(sum((path(3:end, 1:2)-path(1:end-2, 1:2)).^2, 2)), mean(sum((path(4:end, 1:2)-path(1:end-3, 1:2)).^2, 2)), mean(sum((path(5:end, 1:2)-path(1:end-4, 1:2)).^2, 2))];
            
            [p1, p2] = calcMSD3(path, L, min_dt, max_dt, dt);
            
            fitp1 = [fitp1; p1];
            fitp2 = [fitp2; p2];
            
            mean_motion(end);
            
            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
        
        if leng > 2 & sqrt(sum((path(1, 1:2) - path(end, 1:2)).^2)) < 2
            track_pos = [track_pos; [mean(path(:,1:2),1) , path(1,3), leng]];
        end
        
    end
    
    if length(track_pos) == 0, track_pos = [0, 0, 1, 1], end
    
    dist_mat = zeros(max(1,length(track_pos(:,1))));
    for dd = 1:length(track_pos(:,1))
        for ff = 1:length(track_pos(:,1))
            dist_mat(dd,ff) = sqrt(sum((track_pos(dd, 1:2) - track_pos(ff, 1:2)).^2));
            if dd ==ff
                dist_mat(dd,ff) = Inf;
            end
        end
    end
    
    mean(min(dist_mat)' < 2)
    density = length(track_pos(:,1))/(512*512)
    th_proba = 9*pi/(512*512)
    
    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
    
    nb_frames = max(paths(:,3))
    
    max_dist = 2
    loc_seq = []
    for j = 1:length(track_pos(:,1))
        x = track_pos(j,1); y = track_pos(j,2); lseq = zeros(1, nb_frames);
        for s = 1:length(track_pos(:,1))
            x_s = track_pos(s,1); y_s = track_pos(s,2);
            if sqrt((x - x_s)^2 + (y - y_s)^2) < max_dist
                lseq(track_pos(s,3):track_pos(s,3) + track_pos(s,4)) = ...
                lseq(track_pos(s,3):track_pos(s,3) + track_pos(s,4)) +1+ mod(j, 5);
            end
        end
        loc_seq = [loc_seq; lseq];
    end
    
    %pdist2(track_pos, track_pos);
    
    100*4*pi/500^2
    mask = loc_seq>0;
    
    [m ,max_idx] = max(mask.* (1:nb_frames),[], 2)
    [m ,min_idx] = max(-mask.* (-nb_frames:-1),[], 2)

    pi*2^2 / (512*512) * 200
    
    sparsness = (max_idx - min_idx + 1 - sum(mask, 2))% ./sum(mask, 2)
    
    %figure(7), histogram(sparsness, 3.5:10:60)
    spot_percs = [spot_percs, sparsness'>=3]

    %[x, y] = find(loc_seq>0)
    %[m ,idx] = min(,[], 2)
    
    figure(55), imagesc(loc_seq)
    colormap jet
    grid
    
    mean(sparsness'>=3)
    %pause
%     idx = 0
%     ii = 34
%     x = track_pos(ii,1); y = track_pos(ii,2)
%     for kk = min_idx(ii):max_idx(ii)
%         idx = idx + 1
%         img = double(imread(strcat(img_path, img_dir(kk).name)));
%         figure(15), hold on, imagesc(img)
%         ylim([x-6,x+6])
%         xlim([y-6,y+6])
%         
%         mean2(img(floor(x):ceil(x),floor(y):ceil(y)))
%         
%         caxis([1000 6000])
%         colormap jet
%         pause
%     end
end

spot_perc = mean(spot_percs)

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