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

Working tests

parent 16182096

power_dist_test.m

+148 −4
Original line number Diff line number Diff line
@@ -32,8 +32,9 @@ for i = 1:length(mario) 
    % fits.closeFile(file);

    t_raw = fitsread(filoc,"binarytable");

    t_raw = t_raw{1};

    M=fix((t_raw(end)-t_raw(1))/Tseg);

    M=fix((t_raw(end)-t_raw(1))/Tseg)

    for m = 1:M

        disp(['Segment ',num2str(m),' of ',num2str(M),' in file',filoc,' '])

        tsplit = t_raw((t_raw>(t_raw(1)+Tseg*(m-1)))&(t_raw<(t_raw(1)+Tseg*m)));

        % t_raw = t_raw(t_raw<(t_raw(1)+Tseg));

        tsplit = tsplit - tsplit(1);
@@ -49,17 +50,28 @@ for i = 1:length(mario) 
        end

        Yqlbin(:) = Yqlbin(:)./(Fqlbin(:).*Tseg);



        tiledlayout(1,1)

        nexttile



        yme = mean(Y0(100*Tseg+1:end));

        semilogx(Fqlmid,Yqlbin)

        xlim([100 3000])

        xlim([10 3000])

        ylim([yme*0.8 yme*1.2])

        yline(yme)

        xline(100)



        pause





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

            maxcount = max(xtemp);

            datadist = zeros(1,maxcount);

            for bu = 1:maxcount

                datadist(bu) = length(xtemp(xtemp == bu))/N;

            end

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

            mi = -log(ptot0);

            extrpoiss = poisspdf([1:maxcount],mi);

            ncr = 8;

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

            pcr = 1- (1-epscr)^(1/ncr);
@@ -69,18 +81,25 @@ for i = 1:length(mario) 
            expectot = expec1*mi

            vartot = var1*mi + mi*expec1^2

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

            p8_1 = zeros(1,maxcount);

            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)));

            for k = 6:maxcount

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

            end

            vard = var(xtemp)

            vsuetot = vartot/expectot

            vsuedat = vard/mean(xtemp)

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

            meanpow/vsuetot

            meanpow/vsuedat

            extrtot8 = ptotk(p8_1,extrpoiss);





            ncr = 4;

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

            pcr = 1- (1-epscr)^(1/ncr);
@@ -90,6 +109,7 @@ for i = 1:length(mario) 
            expectot = expec1*mi

            vartot = var1*mi + mi*expec1^2

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

            p4_1 = zeros(1,maxcount);

            p4_1(1) = (1-epscr);

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

            p4_1(3) = (18*(pcr^2)*(qcr^8));
@@ -97,7 +117,6 @@ for i = 1:length(mario) 
            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)
@@ -106,6 +125,16 @@ for i = 1:length(mario) 
            meanpow = mean(Y0(100*Tseg+1:end))/2

            meanpow/vsuetot

            meanpow/vsuedat

            extrtot4 = ptotk(p4_1,extrpoiss);

            tiledlayout(2,1)

            nexttile

            semilogy([1:maxcount],extrtot4,[1:maxcount],extrtot8,[1:maxcount],datadist,'+')

            legend('extrtot4','extrtot8','datadist')

            nexttile

            semilogy([1:maxcount],abs(extrtot4-datadist)./datadist,[1:maxcount],abs(extrtot8-datadist)./datadist)

            ylim([0.0001 100])

            legend('abs(extrtot4-datadist)/datadist','abs(extrtot8-datadist)/datadist')

            yline([0.01 0.1 0.5 1.0])

        else

            disp('Too many counts!')

            disp(['Avg. counts per bin :',num2str(mean(xtemp)),'.']);
@@ -114,7 +143,122 @@ for i = 1:length(mario) 
    end

end



%% 

filoc = ['/data/Sorgenti/power_dist_test/dark.mat']

load(filoc)

% file = fits.openFile(filoc)

% fits.closeFile(file);

% t_raw = fitsread(filoc,"binarytable");

t_raw = ToAs2;

M=fix((t_raw(end)-t_raw(1))/Tseg)

for m = 1:M

    disp(['Segment ',num2str(m),' of ',num2str(M),' in file',filoc,' '])

    tsplit = t_raw((t_raw>(t_raw(1)+Tseg*(m-1)))&(t_raw<(t_raw(1)+Tseg*m)));

    % t_raw = t_raw(t_raw<(t_raw(1)+Tseg));

    tsplit = tsplit - tsplit(1);

    Nphot = length(tsplit);

    [xtemp,~]=(histcounts(tsplit,ted));

    xtemp = gpuArray(xtemp);

    Y0 = abs(fft(xtemp));

    Y0 = 2.0.*(Y0(1:N/2+1).^2)./Y0(1);

    % Y02 = Y0(end/2:end);

    Yqlbin = zeros(size(Fqlbin));

    for nino = 1:nbqlog

        Yqlbin(nino) = sum(Y0((Fqlog(nino)<=F0)&(F0<Fqlog(nino+1))));

    end

    Yqlbin(:) = Yqlbin(:)./(Fqlbin(:).*Tseg);



    tiledlayout(1,1)

    nexttile



    yme = mean(Y0(100*Tseg+1:end));

    semilogx(Fqlmid,Yqlbin)

    xlim([10 3000])

    ylim([yme*0.8 yme*1.2])

    yline(yme)

    xline(100)



    pause





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

        maxcount = max(xtemp);

        datadist = zeros(1,maxcount);

        for bu = 1:maxcount

            datadist(bu) = length(xtemp(xtemp == bu))/N;

        end

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

        mi = -log(ptot0);

        extrpoiss = poisspdf([1:maxcount],mi);

        ncr = 8;

        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);

        p8_1 = zeros(1,maxcount);

        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)));

        for k = 6:maxcount

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

        end

        vard = var(xtemp)

        vsuetot = vartot/expectot

        vsuedat = vard/mean(xtemp)

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

        meanpow/vsuetot

        meanpow/vsuedat

        extrtot8 = ptotk(p8_1,extrpoiss);





        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 = zeros(1,maxcount);

        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

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

        end

        vard = var(xtemp)

        vsuetot = vartot/expectot

        vsuedat = vard/mean(xtemp)

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

        meanpow/vsuetot

        meanpow/vsuedat

        extrtot4 = ptotk(p4_1,extrpoiss);

        tiledlayout(2,1)

        nexttile

        semilogy([1:maxcount],extrtot4,[1:maxcount],extrtot8,[1:maxcount],datadist,'+')

        legend('extrtot4','extrtot8','datadist')

        nexttile

        semilogy([1:maxcount],abs(extrtot4-datadist)./datadist,[1:maxcount],abs(extrtot8-datadist)./datadist)

        ylim([0.0001 100])

        legend('abs(extrtot4-datadist)/datadist','abs(extrtot8-datadist)/datadist')

        yline([0.01 0.1 0.5 1.0])

    else

        disp('Too many counts!')

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

    end

    pause

end

% % xtemp = x(:,1);

% % length(xtemp(xtemp == 0))

% % length(xtemp(xtemp == 0))/N
@@ -222,7 +366,7 @@ end 
% datadist(i) = length(xtemp(xtemp == i))/N;

% end

%%

extrtot = ptotk(p4_1,extrpoiss)

% extrtot = ptotk(p4_1,extrpoiss)



function res = ptotk(pncr_1,extrprim)

    %Calculate crosstalk prob. as in Gallego+13