Commit 645a9abf authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

Finished powdist_f, now to be tested and added.

parent 71dff050

powdist_f.m

+124 −116
Original line number Diff line number Diff line 
function [outputArg1,outputArg2] = powdist_f(NyqFreq,LC,Binned,PSname,Noisename)

function [extrVoverE,CrossTalkProb,HalfMeanPow] = powdist_f(NyqFreq,LC,Binned,PSname,Noisename)

%POWDIST_F PS, rebinning and crosstalk estimate; graphs for visual inspection 

%   Detailed explanation goes here

%   Input arguments:

%   NyqFreq is the Nyquist frequency in Hz 

%   LC is the lightcurve 

%   If Binned is true, we assume uniform binning and dt = 1/2*NyqFreq,

%   otherwise we assume LC are actually the times of arrival

%   PSname and Noisename are the filenames (or complete paths) of graphs to

%   be saved, resp. PSs+lightcurve and counts dist. with crosstalk estimate

%   Output arguments:

%   extrVoverE are the three values of the expected Var/E assuming 8-pixel,

%   4-pixel, and 1-pixel crosstalk probability distributions

%   CrossTalkProb is the agnostic crosstalk porb. (only assumes underlying

%   poissonian noise)

%   HalfMeanPow is 0.5 the average high-freq. power in the power spectrum

%   (this is what extrVoverE has to be compared to)



% import matlab.io.* prob./ unneeded



%% Part I: Power spectrum, quasi-log. rebinning and light-curve printing

dt = 1/(2*NyqFreq);

Hfreq = 1000;



if (~Binned)

    N = fix((LC(end)-LC(1))/dt);
@@ -58,12 +66,12 @@ xlabel('Frequency [Hz]') 
ylabel('Leahy-normalized power (no corr.)')



nexttile

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

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

semilogx(Fqlmid,Yqlbin)

xlim([10 3000])

xlim([10 NyqFreq+1000])

ylim([yme*0.8 yme*1.2])

yline(yme,'--','Avg. pow. for f > 100 Hz')

xline(100,':')

yline(yme,'--',['Avg. pow. for f > ',num2str(fix(Hfreq)),' Hz'])

xline(Hfreq,':')

xlabel('Frequency [Hz]')

ylabel('Leahy-normalized power (no corr.)')

title('Quasi-logarithmically-rebinned PS')
@@ -87,7 +95,7 @@ if (isempty(xtemp(xtemp ==0))) 
        ' and no bin with 0 counts.']);

    if(~Binned)

        disp('Will rebin the times of arrival on a finer grid')

        while (length(xtemp(xtemp ==0))/N<exp(-1))

        while (length(xtemp(xtemp ==0))<N*exp(-1))

            dt = dt/2;

            ted = (0:dt:Tseg);

            [xtemp,~]=(histcounts(LC,ted));
@@ -98,6 +106,10 @@ if (isempty(xtemp(xtemp ==0))) 
    else

        disp('Data were provided already binned, nothing we can do')

        disp('Try with times of arrival instead, and set Binned to false')

        extrVoverE = NaN;

        CrossTalkProb = NaN;

        HalfMeanPow = mean(Y0(Hfreq*Tseg+1:end))/2;

        return

    end



end
@@ -111,8 +123,11 @@ end 
mi = -log(ptot0);

% mivalues(i,m) = mi;

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

    ncr = 8;

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

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

vsuedat = var(xtemp)/mean(xtemp);



ncr = 8;

% epsvalues(i,m) = epscr;

% epsmedio = epsmedio + epscr/M;

pcr = 1- (1-epscr)^(1/ncr);
@@ -128,28 +143,22 @@ end 
for k = 6:maxcount

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

end

    expec1 = p8_1(5)*(1+4*r8_1)/r8_1/r8_1 + sum((1:4).*p8_1(1:4));

expec1 = sum((1:maxcount).*p8_1(1:maxcount));

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

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

expectot = expec1*mi;

vartot = var1*mi + mi*expec1^2;

    vard = var(xtemp);

    vsuetot = vartot/expectot;

    vsuedat = vard/mean(xtemp);

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

vsuetot(1) = vartot/expectot;

% varvalues(i,m,1) = 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;

% 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);
@@ -161,50 +170,54 @@ end 
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;

expec1 = sum((1:maxcount).*p4_1(1:maxcount));

expectot = expec1*mi;

vartot = var1*mi + mi*expec1^2;

vsuetot(2) = vartot/expectot;

% varvalues(i,m,2) = meanpow/vsuetot;

% meanpow/vsuedat

extrtot4 = ptotk(p4_1,extrpoiss);



% The case for n = 1 is especially simple: it's a Polya–Aeppli dist.;

% [see Univ. Disc. Dist. - Johnson, Kemp, and Kotz, sec 9.7, same

% notation, while in their chap. 5 their p is our q and vice versa]

ncr = 1;

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

pcr = epscr;

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

qcr = (1-epscr);

p1_1 = zeros(1,maxcount);

p1_1(1:maxcount) = (pcr.^(0:maxcount-1)).*(1-pcr);

    vard = var(xtemp);

    vsuetot = vartot/expectot;

    vsuedat = vard/mean(xtemp);

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

    varvalues(i,m,3) = meanpow/vsuetot;

expec1 = 1/qcr;

var1 = pcr/(qcr^2);

expectot = expec1*mi;

vartot = var1*mi + mi*expec1^2;

vsuetot(3) = vartot/expectot;

% varvalues(i,m,3) = meanpow/vsuetot;

% meanpow/vsuedat

extrtot1 = ptotk(p1_1,extrpoiss);







    if (i == 7)

        tiledlayout(2,1)

        nexttile

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

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

blindpoiss = poisspdf((1:maxcount),mean(xtemp));

t = tiledlayout(2,1);

nexttile

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

        ylim([0.0001 100])

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

legend('1-pixel crosstalk','4-pixel crosstalk','8-pixel crosstalk','Pure Poissonian','Data distr.')

xlabel('Counts')

ylabel('Number of bins with x counts')

title('Extrapolated gener. Poiss. dist. vs data')

babbo = nexttile;

semilogy([1:maxcount],abs(extrtot1-datadist)./datadist,'x',[1:maxcount],abs(extrtot4-datadist)./datadist,'o',[1:maxcount],abs(extrtot8-datadist)./datadist,'*',[1:maxcount],abs(blindpoiss-datadist)./datadist,'square')

ylim([0.0001 2])

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

        yline([0.01 0.1 0.5 1.0])

        pause

    end

yline([0.01 0.1 0.5 1.0],'--')

babbo.YGrid = 'on';

xlabel('Counts')

ylabel('|Extrapolated - data|/data')

title('Relative errors')

exportgraphics(t,Noisename)



outputArg1 = inputArg1;

outputArg2 = inputArg2;

extrVoverE = vsuetot;

CrossTalkProb = epscr;

HalfMeanPow = meanpow;

end



function res = ptotk(pncr_1,extrprim)
@@ -225,10 +238,5 @@ function res = ptotk(pncr_1,extrprim) 
            % 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

    res = extrprim*pncr_m;

end
 
 No newline at end of file

power_dist_test.m

+9 −8
Original line number Diff line number Diff line
@@ -164,7 +164,7 @@ for i = 1:length(mario) 
            





            if (i == 7)

            % if (i == 7)

            tiledlayout(2,1)

            nexttile

            semilogy([1:maxcount],extrtot1,[1:maxcount],extrtot4,[1:maxcount],extrtot8,[1:maxcount],datadist,'+')
@@ -175,7 +175,7 @@ for i = 1:length(mario) 
            % legend('abs(extrtot1-datadist)/datadist','abs(extrtot4-datadist)/datadist','abs(extrtot8-datadist)/datadist')

            yline([0.01 0.1 0.5 1.0])

            pause

            end

            % end

        else

            disp('Too many counts!')

            disp(['Avg. counts per bin :',num2str(mean(xtemp)),'.']);
@@ -464,10 +464,11 @@ function res = ptotk(pncr_1,extrprim) 
            % 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

%     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

    res = extrprim*pncr_m;

end
 
 No newline at end of file

scsearch.m

+2 −2
Original line number Diff line number Diff line
@@ -195,8 +195,8 @@ end 


% Parameters for Sco X-1

 % s 

f_max = 1050

f_min = 948

f_max = 1500

f_min = 1448

f_step = 1/Tseg

f_gr = (f_min:f_step:f_max).';

Porb = 68023.919 % s