function [peaknb] = sptFromPosList(fdir,fname,PathC_name,trackingInCells,...
    ImagingInterval, maxDisplacement, ToSave, Res_max, Int_min, std_max, raw_th, filter_th)

mkdir(strcat(fdir,'Tracks/',PathC_name,fname,'/negmask'))

if trackingInCells == 1
    n = ['.' filesep 'Tracks' filesep fname filesep];
    path_n = strcat('.', filesep, 'Tracks', filesep, PathC_name, fname, filesep);
else
    n = ['.' filesep 'Tracks' filesep fname filesep 'negmask' filesep];
    path_n = strcat('.', filesep, 'Tracks', filesep, PathC_name, fname, filesep, 'negmask', filesep);
end

fileList = dir([n,'*-*.mat']);
%dont load pathsC
name = {fileList.name};
a1 = ~contains(name,'pathsC');
name = name(a1);
%%
AllDcent = []; Allsteps = []; Alllens = []; AlldistC = [];Allangles = [];

StorePaths = [];
NB = 0;
KO = []; pathsNoTrafo = [];
fast = 0; slow = 0;
AnzPfade = 0; simSpotsCell = []; DiffAll = []; SpotsPerCell = []; TracksPerCell = [];
hl = jet(15);

% tracking parameters. (see trackCommentsshort for description)
param.mem = 0;
param.dim = 2;
param.good = 0;
param.quiet = 1;
param.stepLength = maxDisplacement;

TStep = ImagingInterval;

initial_stds = []; initial_fitted_intensity = []; initial_residuals = []; 
initial_raw_intensity = []; initial_raw_norm_intensity = []; initial_filter_intensity = []; 
final_stds = []; final_fitted_intensity = []; final_residuals = []; 
final_raw_intensity = []; final_raw_norm_intensity = []; final_filter_intensity = [];
median_stds = []; median_fitted_intensity = []; median_residuals = []; 
median_raw_intensity = []; median_raw_norm_intensity = []; median_filter_intensity = [];
cell_areas = []; dt = 50;

% if the loop breaks start from i+1
for i = 1:length(name)
    try
        clear paths angles vels steps res Dcent PosList fitP peak_details
        
        load([n name{i}])

        cell_areas = [cell_areas, mean2(BW)];
        
        initial_stds = [initial_stds; fitP(:,3)];
        initial_fitted_intensity = [initial_fitted_intensity; fitP(:,2)];
        initial_residuals = [initial_residuals; fitP(:,4)./fitP(:,2)];
        initial_raw_intensity = [initial_raw_intensity; peak_details(:, 5)];
        initial_raw_norm_intensity = [initial_raw_norm_intensity; peak_details(:, 5)./peak_details(:, 2)];
        initial_filter_intensity = [initial_filter_intensity; peak_details(:, 6)];
        
        residuals_mask = fitP(:,4)./fitP(:,2) < Res_max; %fitP(:,4)./fitP(:,2) = residuals normalised by peak intensity
        std_mask = fitP(:,3)<std_max; %std = peak width
        raw_th_mask = peak_details(:, 5) > raw_th*peak_details(:, 2);
        filter_th_mask = peak_details(:, 6) > filter_th*peak_details(:, 4);
        
        global_mask = residuals_mask & std_mask & raw_th_mask & filter_th_mask;
        
        PosList = PosList(global_mask,:);
        fitP = fitP(global_mask,:);
        peak_details = peak_details(global_mask,:);
        
        if Int_min == "relative"
            Modes = []; % calculate the mode for different binnings and take the median one
            for ii =5:50
                for j = 0:0.2:1
                    Modes = [Modes, (ii*mode(floor(fitP(:,2)/ii + j)))*0.5];
                end
            end
            adjusted_Int_min = median(Modes); strcat('minimal intensity = ', string(' '), string(adjusted_Int_min))
        
            PosList = PosList(fitP(:,2)>adjusted_Int_min,:); % intensity higher than Int_min
            peak_details = peak_details(fitP(:,2)>adjusted_Int_min,:);
            fitP = fitP(fitP(:,2)>adjusted_Int_min,:);
            
            PosList = PosList(fitP(:,2)<adjusted_Int_min*10,:); % intensity higher than Int_min
            peak_details = peak_details(fitP(:,2)<adjusted_Int_min*10,:);
            fitP = fitP(fitP(:,2)<adjusted_Int_min*10,:);
        elseif Int_min == "percentage"
            perc = 0.01; % remove the outlayers to keep the closer to the mean values 
            lims = [quantile(fitP(:,2), perc), quantile(fitP(:,2), 1-perc)];
            Int_mask = (fitP(:,2) >  lims(1)) & (fitP(:,2) <  lims(2));
            
            strcat(string('Int bounds = ['), string(round(lims(1))), string(', '), string(round(lims(2))), ']')
            
            PosList = PosList(Int_mask,:); % keep the peaks with intensities belonging to the 
            fitP = fitP(Int_mask,:);
            peak_details = peak_details(Int_mask,:);
        else
            PosList = PosList(fitP(:,2)>Int_min,:); % intensity higher than Int_min
            peak_details = peak_details(fitP(:,2)>Int_min,:);
            fitP = fitP(fitP(:,2)>Int_min,:);
        end
        
        final_stds = [final_stds; fitP(:,3)];
        final_fitted_intensity = [final_fitted_intensity; fitP(:,2)];
        final_residuals = [final_residuals; fitP(:,4)./fitP(:,2)];
        final_raw_intensity = [final_raw_intensity; peak_details(:, 5)];
        final_raw_norm_intensity = [final_raw_norm_intensity; peak_details(:, 5)./peak_details(:, 2)];
        final_filter_intensity = [final_filter_intensity; peak_details(:, 6)];

        peaknb = length(fitP(:,3))
        
        time = PosList(:,3);
        
        median_stds = [median_stds; time_median(fitP(:,3), time, dt)];
        median_fitted_intensity = [median_fitted_intensity; time_median(fitP(:,2), time, dt)];
        median_residuals = [median_residuals; time_median(fitP(:,4)./fitP(:,2), time, dt)];
        median_raw_intensity = [median_raw_intensity; time_median(peak_details(:, 5), time, dt)];
        median_raw_norm_intensity = [median_raw_norm_intensity; time_median(peak_details(:, 5)./peak_details(:, 2), time, dt)];
        median_filter_intensity = [median_filter_intensity; time_median(peak_details(:, 6), time, dt)];
        
        TStep = ImagingInterval;
        dummyCell = [];
        for ceI = 1:length(cells)
            if not(isempty(cells(ceI).Area))
                dummyCell = [dummyCell; cells(ceI)];
            end
        end
        clear cells
        cells = dummyCell;
        
        str = [fdir refimgMreB];
        str = strrep(str,'/',filesep); str = strrep(str,'\',filesep);
        if (exist(str, 'file') == 2)       
            IMG = double(imread(str));
        else
            IMG = zeros([512, 512]); str
        end
        
        if size(BW, 1) > 512 || size(BW, 2) > 512
            BW = BW(1:512, 1:512);
        end
        IMG = IMG ./ max(max(IMG(BW > 0)));

        % Add a fourth column to the position list with the cell number as
        % found in the mask
        POSLIST = [];
        FITP = [];
        for j = 1:max(PosList(:,3))
            ind = find(PosList(:,3)==j);
            if isempty(ind)
                POSLIST = [POSLIST; 0 0 j 0];
                FITP = [FITP; 0 0 0 0];
            else
                for k = ind'
                    cellNB = BW(round(PosList(k,1)), round(PosList(k,2)));
                    if cellNB == 0
                        cellNB = BW(ceil(PosList(k,1)), ceil(PosList(k,2)));
                        if cellNB == 0
                            cellNB = BW(floor(PosList(k,1)), floor(PosList(k,2)));
                        end
                    end
                    if cellNB == 0
                        POSLIST = [POSLIST; 0 0 j 0];
                        FITP = [FITP; 0 0 0 0];
                    else
                        POSLIST = [POSLIST; [PosList(k,1:3) cellNB]];
                        FITP = [FITP; fitP(k,:)];
                    end
                end
            end
        end
        PosList = POSLIST;
        fitP = FITP;

        PosListGem = [];
        for zeit = 1:max(POSLIST(:,3)) 
            subPosList = POSLIST(POSLIST(:,3)==zeit, :);
            zellen = unique(subPosList(:,4))'; 
            zellen = zellen(zellen > 0) ;
            if ~isempty(zellen)
                for k = zellen
                    if (sum(subPosList(:,4) == k) > 1)
                        PosListGem = [PosListGem; subPosList(subPosList(:,4) == k,:)];
                    end
                end
            end
        end
        if ~isempty(PosListGem)
            simSpotsCell = [];
            for j = unique(PosListGem(:,3))'
                subPosListGem = PosListGem(PosListGem(:,3) == j, :);
                for k = unique(subPosListGem(:,4))'
                    simSpotsCell = [simSpotsCell; sum(subPosListGem(:,4) == k)];
                end
            end
        end

        zellen = unique(PosList(:,4))';
        zellen = zellen(zellen > 0) ;
        if ~isempty(zellen)
            for k = zellen
                SpotsPerCell = [SpotsPerCell; sum(PosList(:,4) == k)];
            end
        end

        idx = find(PosList(:,4) ~= 0);
        PosList = PosList(idx,:); 
        fitP = fitP(idx,:);
        for iT = 1:floor(max(PosList(:,3))/10)*10/10+1 
            TimeGrob = floor(PosList(:,3)/10)*10/10+1;
            idxTime = find(TimeGrob == iT);

        end
        
        % Parameters of the cells   %%from regionprops in mask.mat
        alpha= cat(1, cells.Orientation);
        AxMin = cat(1, cells.MinorAxisLength);
        AxMaj = cat(1, cells.MajorAxisLength);
        Centr = cat(1, cells.Centroid);
        
        % tracking parameters. (see trackCommentsshort for description)
        param.mem = 0;
        param.dim = 2;
        param.good = 0;
        param.quiet = 1;
        param.stepLength = maxDisplacement;
        % if two peaks present in the maxDisplacement px area, track is stopped.
        
        % The tracking algorithm expects a time vector that has constant
        % increments. Therefore, if there are frames with no peaks detected, we
        % have to add am empty row in the position list. The 'fake' peaks have
        % coordinates of (0, 0) to not confuse them with real peaks and they
        % are sorted out again later.
        PList = [];
        tempfitP = [];
        for tim = 1:max(PosList(:,3))
            idxT = find(PosList(:,3) == tim);
            if ~isempty(idxT)
                PList = [PList; PosList(idxT,:)];
                tempfitP = [tempfitP; fitP(idxT,:)];
            else
                PList = [PList; 0 0 tim 0];
                tempfitP = [tempfitP; 0 0 0 0];
            end
        end
        
        % Call tracking function, based on Crocker Grier
        [res, fitP] = trackCommentsshort(PList(:,1:3),param.stepLength, param, tempfitP);
        ind = find(res(:,1) ~= 0);
        res = res(ind,:);
        fitP = fitP(ind,:);
        % res: x-Coord | y-Coord | time [frame] | path number

        meanResiduals = [];
        meanStd = [];
        for trackID = 1:max(res(:,4))
            idx = find(res(:,4) == trackID);
            meanResiduals(idx,1) = mean(fitP(idx,4));
            meanStd(idx,1) = mean(fitP(idx,3));
        end
        
        Dind = 1;
        resalt = res;

        for iconv = unique(res(:,4))'
            clear  a b PathX PathY indices cellNB
            indices = find(res(:,4) == iconv);
            PathX = res(find(res(:,4) == iconv) ,1);
            PathY = res(find(res(:,4) == iconv) ,2);
            cellNB = PosList(find((PathX(1) == PosList(:,1)) & (PathY(1) == PosList(:,2))),4); % find Cell number
            res(indices, 5) = cellNB(1);
        end
        ce = cat(1, cells.Centroid);
        MjAx = cat(1, cells.MajorAxisLength);

        %%
        paths = []; vels = []; angles = []; Dangles = []; steps = []; nbPath = 0; Dcents = []; lens = []; Int = [];

        for j = 1:max(res(:,4))
            clear path angle vel step Dangle Dcent len
            path_idx = find(res(:,4) == j);
            path = res(path_idx,:);
            pathalt = resalt(path_idx,:);
            leng = length(path(:,1));
            if leng >= 1
                NB = NB + 1;
                nbPath = nbPath + 1;

                if leng > 1
                    for k = 1:length(path(:,1))-1
                        PosList(find((PathX(1) == PosList(:,1)) & (PathY(1) == PosList(:,2))),4);
                        step(k) = sqrt((path(k,1)-path(k+1,1))^2 + (path(k,2)-path(k+1,2))^2);
                        cellNB = res(k,5);
                        Dcent(k) = sqrt((path(k,1)-ce(cellNB,1))^2 + (path(k,2)-ce(cellNB,2))^2)/(MjAx(cellNB)/2);
                        vel(k) = step(k)/(path(k+1,3) - path(k,3));
                        angle(k) = acos((path(k,1)-path(k+1,1))/step(k))/pi*180;
                        len(k) = leng;
                        
                        if k > 1
                            Dangle(k-1) = angle(k) - angle(k-1);
                        else
                            Dangle = [];
                        end
                    end
                    angles = [angles angle];
                    Dangles = [Dangles Dangle];
                    Dcents = [Dcents Dcent];
                    vels = [vels vel];
                    steps = [steps step];
                    lens = [lens len];
                end
                
                paths = [paths; path];
                pathsNoTrafo = [pathsNoTrafo; pathalt];
                
                if leng >= 5 % only show paths with more than 4 steps.
                    clear Path PathNEW Steps MreBInt MreBIntS Pathneu PathneuNEW msd
                    Path = resalt(path_idx,:);
                    Pathneu = res(path_idx,1:4);
                    AnzPfade = AnzPfade + 1;
                    SD = [];
                    for p = 1:size(Path,1)
                        PathNEW(p,:) = [Path(p,1) Path(p,2) 1] ;
                        PathneuNEW(p,:) = [Pathneu(p,1) Pathneu(p,2) 1] ;
                    end
                    
                    msd = calcMSD(Path);
                    limit = length(msd);
                    x = (1:limit)*(TStep*0.001);
                    fo = fitoptions('Method','NonlinearLeastSquares','Robust','Off', 'Algorithm',...
                        'Trust-Region', 'Lower', [0 0], 'Upper',[0.0014*3 4],'StartPoint',...
                        [0 0], 'MaxFunEvals', 600, 'MaxIter', 400);
                    ft = fittype('m*x + b','options',fo);
                    [curve,gof] = fit(x',msd'*0.075^2,ft);
                    DiffAll = [DiffAll curve.m/4];
                    
                    Cx = mean(PathNEW(:,1));
                    Cy = mean(PathNEW(:,2));
                    for p = 1:size(Path,1)
                        distC(p) = (PathNEW(p,1)-Cx)^2 + (PathNEW(p,2)-Cy)^2;
                    end
                    AlldistC = [AlldistC distC];
                    if DiffAll(end) > 0.025%0.025
                        fast = fast + 1;
                    else
                        slow  = slow + 1;
                    end
                    StorePaths = [StorePaths; PathNEW(:, 1:2) ones(size(PathNEW, 1), 1)*...
                        AnzPfade ones(size(PathNEW, 1), 1)*DiffAll(end)];        
                end
                if leng > 4
                end
            end
        end

        AllDcent  = [AllDcent Dcents];
        Allsteps  = [Allsteps steps];
        Alllens   = [Alllens lens];
        Allangles = [Allangles angles];

        disp([max(size(PList)) length(paths)])
    if(ToSave)
        na = strrep(name{i},'.mat','');
        na = strcat(path_n, na, '-pathsC.mat');
        save('/home/mcm/temp', 'paths','param', 'TStep', 'nbPath', 'cells', 'refimg', 'pathsNoTrafo', 'fitP', 'BW', 'peaknb')
        movefile(['/home/mcm/temp.mat'], na)
    end
    catch
        disp('no peaks found on movie')
        paths = []; nbPath = 0; 
    end
end

save('/home/mcm/temp', 'initial_stds', 'initial_fitted_intensity', 'initial_raw_intensity', 'initial_raw_norm_intensity',...
'initial_filter_intensity', 'final_stds', 'final_fitted_intensity', 'final_raw_intensity', ...
'final_raw_norm_intensity', 'final_filter_intensity', 'median_stds', 'median_fitted_intensity',...
'median_raw_intensity', 'median_raw_norm_intensity', 'median_filter_intensity', 'cell_areas', 'dt')

movefile(['/home/mcm/temp.mat'], strcat(path_n, 'peak_details.mat'))
