Commit 16182096 authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

More strides

parent b06386cd

power_dist_test.m

+89 −11
Original line number Diff line number Diff line
@@ -54,10 +54,10 @@ for i = 1:length(mario) 
        xlim([100 3000])

        ylim([yme*0.8 yme*1.2])

        yline(yme)

        pause





        if (~isempty(xtemp(xtemp ==0)))

            maxcount = max(xtemp);

            ptot0 = length(xtemp(xtemp == 0))/N;

            mi = -log(ptot0);

            ncr = 8;
@@ -69,14 +69,36 @@ for i = 1:length(mario) 
            expectot = expec1*mi

            vartot = var1*mi + mi*expec1^2

            qcr = (1-epscr)^(1/ncr);

            p8l(1) = (1-epscr);

            p8l(2) = (8*pcr*qcr^14);

            p8l(3) = (84*(pcr^2)*(qcr^20));

            p8l(4) = 24*(pcr^3)*(qcr^24)*(1+3*qcr+38*qcr^2);

            p8l(5) = 4*(pcr^4)*(qcr^30)*(180+540*qcr+2521*qcr^2);

            r8l = p8l(5)/(1-sum(p8l(1:4)));

            for k = 6:10

                p8l(k) = p8l(5)*(1-r8l)^(k-5);

            p8_1(1) = (1-epscr);

            p8_1(2) = (8*pcr*qcr^14);

            p8_1(3) = (12*(pcr^2)*(qcr^18)*(1+2*qcr+4*qcr^2));

            p8_1(4) = 4*(pcr^3)*(qcr^20)*(1+3*qcr+14*qcr^2+30*qcr^3+61*qcr^4+59*qcr^5+72*qcr^6);

            p8_1(5) = 5*(pcr^4)*(qcr^24)*(9+36*qcr+98*qcr^2+188*qcr^3+310*qcr^4+372*qcr^5+520*qcr^6+396*qcr^7+341*qcr^8);

            r8_1 = p8_1(5)/(1-sum(p8_1(1:4)));

            vard = var(xtemp)

            vsuetot = vartot/expectot

            vsuedat = vard/mean(xtemp)

            meanpow = mean(Y0(100*Tseg+1:end))/2

            meanpow/vsuetot

            meanpow/vsuedat

            ncr = 4;

            epscr = 1 - (length(xtemp(xtemp == 1))/N)/(mi*ptot0);

            pcr = 1- (1-epscr)^(1/ncr);

            expec1 = 1+ epscr + (3*(ncr-1)/(2*ncr))*epscr^2;

            var1 = epscr + (epscr^2)*(7-9/ncr)/2;

            enf = 1 + var1/(expec1^2);

            expectot = expec1*mi

            vartot = var1*mi + mi*expec1^2

            qcr = (1-epscr)^(1/ncr);

            p4_1(1) = (1-epscr);

            p4_1(2) = (4*pcr*qcr^6);

            p4_1(3) = (18*(pcr^2)*(qcr^8));

            p4_1(4) = 4*(pcr^3)*(qcr^8)*(1+3*qcr+18*qcr^2);

            p4_1(5) = 5*(pcr^4)*(qcr^10)*(8+24*qcr+55*qcr^2);

            r4_1 = p4_1(5)/(1-sum(p4_1(1:4)));

            for k = 6:maxcount

                p8_1(k) = p8_1(5)*(1-r8_1)^(k-5);

                p4_1(k) = p4_1(5)*(1-r4_1)^(k-5);

            end

            vard = var(xtemp)

            vsuetot = vartot/expectot
@@ -86,7 +108,7 @@ for i = 1:length(mario) 
            meanpow/vsuedat

        else

            disp('Too many counts!')

            disp(['Avg. counts per bin :',num2str(mean(xtemp)/N),'.']);

            disp(['Avg. counts per bin :',num2str(mean(xtemp)),'.']);

        end

        pause

    end
@@ -171,3 +193,59 @@ end 
% % [histnosub,~] = histcounts(Y0sub2,powedges,'normalization','pdf');

% % powmid(1:52) = (powedges(2:53)+powedges(1:52))./2;

% % semilogy(powmid,histno,'o',powmid,histnodiv,'+',powmid,chi2pdf(powmid,(2.0)),powmid,histnosub,'*')



%% A function to calculate crosstalk probabilities given any distribution 

%% of primary events and the pncr_1 (see Gallego+2013)



% if (length(pncr_1)< maxcount) 

%     disp('p8_1 was not calculated for the whole distribution of observed counts')

%     disp(['Setting the size of p8_m to that (',num2str(length(pncr_1)),') instead of ',num2str(maxcount),'.'])

%     maxcount = length(pncr_1);

% end

% pncr_m = zeros(maxcount);

% pncr_m(1,:) = pncr_1(:); 

% for m = 1:9 

%     for k = 1:maxcount

%         for i = 1:k-m

%             pncr_m(m+1,k) = pncr_m(m+1,k) + pncr_m(m,k-i)*pncr_1(i);

%         end

%         % pncr_m(m+1,k) = sum(pncr_m(m,k-(1:k-m)).*pncr_1(1:k-m));

%     end

% end

% extrtot = zeros(size(pncr_1));

% for k = 1:maxcount

%     for m = 1:k

%         extrtot(k) = extrtot(k) + extrprim(m).*pncr_m(m,k);

%     end

% end

% for i = 1:maxcount

% datadist(i) = length(xtemp(xtemp == i))/N;

% end

%%

extrtot = ptotk(p4_1,extrpoiss)



function res = ptotk(pncr_1,extrprim)

    %Calculate crosstalk prob. as in Gallego+13

    % if (length(pncr_1)<maxcount)

    %     disp('p8_1 was not calculated for the whole distribution of observed counts')

    %     disp(['Setting the size of p8_m to that (',num2str(length(pncr_1)),') instead of ',num2str(maxcount),'.'])

    %     maxcount = length(pncr_1);

    % end

    maxcount = length(pncr_1);

    pncr_m = zeros(maxcount);

    pncr_m(1,:) = pncr_1(:);

    for m = 1:maxcount-1

        for k = 1:maxcount

            for i = 1:k-m

                pncr_m(m+1,k) = pncr_m(m+1,k) + pncr_m(m,k-i)*pncr_1(i);

            end

            % pncr_m(m+1,k) = sum(pncr_m(m,k-(1:k-m)).*pncr_1(1:k-m));

        end

    end

    res = zeros(size(pncr_1));

    for k = 1:maxcount

        for m = 1:k

            res(k) = res(k) + extrprim(m).*pncr_m(m,k);

        end

    end

end
 
 No newline at end of file