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

Reworked powdist to also plot titular power ditribution vs chi2 with 2 degrees... 

Reworked powdist to also plot titular power ditribution vs chi2 with 2 degrees of freedom, changed SCO_multiobs
parent 5993277e

functions/conf_int.m

+4 −4
Original line number Diff line number Diff line
@@ -5,7 +5,7 @@ function CP_ci = conf_int(x,n,alpha) 
%   trials. If x is a vector, CP_ci is a matrix with size (length(x),2). 

%

%   By default, the optional parameter alpha is set to 0.05 (corresponding 

%   to 95%-confidence intervals.

%   to 95%-confidence intervals).



% %--------------

if nargin < 3
@@ -17,7 +17,7 @@ if n<0 || n~=round(n) || isinf(n) || any(x>n) 
    disp('THE RESULTS WILL MAKE NO SENSE.')

end

if any(x<0)

    disp('All of x must not be negative.')

    disp('All of x must be non-negative.')

    disp('THE RESULTS WILL MAKE NO SENSE.')

end
@@ -35,8 +35,8 @@ maskx0 = logical(x==0); 
maskxn = logical(x==n);

xpu = x(maskxpu);



CP_ci(maskxpu,2) = icdf('Beta',1-alpha/2,xpu+1,n-xpu);

CP_ci(maskxpu,1) = icdf('Beta',alpha/2,xpu,n-xpu+1);

CP_ci(maskxpu,2) = betainv(1-alpha/2,xpu+1,n-xpu);

CP_ci(maskxpu,1) = betainv(alpha/2,xpu,n-xpu+1);

CP_ci(maskx0,2) = 1-((alpha./2).^(1./n));

CP_ci(maskxn,1) = ((alpha./2).^(1./n));

functions/powdist.m

+55 −2
Original line number Diff line number Diff line 
function [extrVoverE,CrossTalkProb,HalfMeanPow] = powdist(NyqFreq,LC,Binned,PSname,Noisename)

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

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

%POWDIST PS, crosstalk estimate, power distribution; graphs to check by eye

%

%   Input arguments:

%   NyqFreq is the Nyquist frequency in Hz
@@ -8,6 +8,10 @@ function [extrVoverE,CrossTalkProb,HalfMeanPow] = powdist(NyqFreq,LC,Binned,PSna 
%   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

%   PowDisname is the filename for the graph containing the comparison of 

%   a chi^2 distribution with two degrees of freedom vs our power

%   distribution, both with "input" powers and the PS renormalized by

%   2*HalfMeanPow

%   

%   Output arguments:

%   extrVoverE are the three values of the expected Var/E assuming 8-pixel,
@@ -211,6 +215,55 @@ exportgraphics(t,Noisename) 
extrVoverE = vsuetot;

CrossTalkProb = epscr;

HalfMeanPow = meanpow;





%% Part III: Power distribution at frequencies >= 10 Hz, comparison with chi^2_2 



PSgeq10 = Y0(ceil(Tseg*10)+1:end);

Fgeq10 = F0(ceil(Tseg*10)+1:end);

maxp = max(PSgeq10);

pedges = (0:1:ceil(maxp));

pmid = (0.5:1:ceil(maxp)-0.5);

Pcounts= fastcounts(PSgeq10,pedges);

sumpc = sum(Pcounts);

Pcountsresc = fastcounts(PSgeq10./(2*HalfMeanPow),pedges);

sumpcresc = sum(Pcounts);



t = tiledlayout(2,1);

nexttile



bintor = conf_int(Pcounts,sumpc);

histogram('BinEdges',pedges,'BinCounts',bintor(:,2),'EdgeColor','none','FaceAlpha',0.5)

hold on

yscale log

histogram('BinEdges',pedges,'BinCounts',bintor(:,1),'EdgeColor','none','FaceAlpha',1,'FaceColor','white')

ylim([0.01 max(Pcounts)*3])

histogram('BinEdges',pedges,'BinCounts',Pcounts,'EdgeColor','black','FaceColor','none','DisplayStyle','stairs')

semilogy(pmid,Pcounts,'o','Color',[0.8500 0.3250 0.0980],'MarkerFaceColor',[0.8500 0.3250 0.0980])

semilogy(pmid,(chi2pdf(pmid,2)).*(sumpc),'Color',[0.4660 0.6740 0.1880],'LineWidth',2)

xlabel('Power')

ylabel('Number of bins')

title('Distribution of powers for frequencies above 10 Hz')

hold off



nexttile 

bintor = conf_int(Pcountsresc,sumpcresc);

histogram('BinEdges',pedges,'BinCounts',bintor(:,2),'EdgeColor','none','FaceAlpha',0.5)

hold on

yscale log

histogram('BinEdges',pedges,'BinCounts',bintor(:,1),'EdgeColor','none','FaceAlpha',1,'FaceColor','white')

ylim([0.01 max(Pcountsresc)*3])

histogram('BinEdges',pedges,'BinCounts',Pcountsresc,'EdgeColor','black','FaceColor','none','DisplayStyle','stairs')

semilogy(pmid,Pcountsresc,'o','Color',[0.8500 0.3250 0.0980],'MarkerFaceColor',[0.8500 0.3250 0.0980])

semilogy(pmid,(chi2pdf(pmid,2)).*(sumpcresc),'Color',[0.4660 0.6740 0.1880],'LineWidth',2)

xlabel('Power')

ylabel('Number of bins')

title('Distribution of rescaled powers for frequencies above 10 Hz')

hold off

exportgraphics(t,PowDisname)







end



function res = ptotk(pncr_1,extrprim)

scsearch.m

+6 −4
Original line number Diff line number Diff line
@@ -97,11 +97,13 @@ toc 
pathfi = [folname,'/'];



%Rebin -------------------------------------------------------------------

Nyq = 2000 %cambiare in caso

Tseg = 256

Nyqmin = 2000; %cambiare in caso

Tsegmin = 256;

Nyq = Nyqmin



dt_psd=1/(2*Nyq); %risoluzione spettro di potenza, diverso dal dt dei segmenti

% nbin=round((t_raw(end)-t_raw(1))/dt_psd);

dt_psd=1/(2*Nyq); %

N = pow2(ceil(log2(Tsegmin/dt_psd)))

Tseg = N*dt_psd



countsname = sprintf('counts_check_Tseg_%d.mat',Tseg);

%%

scsearch_SCO_multiobs_multifreq.m

+23 −16
Original line number Diff line number Diff line
@@ -48,11 +48,13 @@ toc 


pathfi = [folname,'/']



Nyq = 2000 %cambiare in caso

Tseg = 512

dt_psd=1/(2*Nyq);

Nyqmin = 2000; %cambiare in caso

Tsegmin = 512;

Nyq = Nyqmin

dt_psd=1/(2*Nyq); %

N = ceil(Tsegmin/dt_psd); % for now

Tseg = N*dt_psd

dt = dt_psd;

N = fix(Tseg/dt);

f_supamin = 49 % Hz

f_supamax = 1550 % Hz

fr_int_width = 100 % Hz (-1)
@@ -63,7 +65,7 @@ porb_min = Porb - 3*porb_unc; % s 
porb_max = Porb + 3*porb_unc; % s

porb_gr = (Porb).';

Tasc_old = 56722.6303715306713

Tasc_old_unc = 3*33 + 4; % s (3sigma Killestein + 4 s GPS def. uncertainty)

Tasc_old_unc = 3*(33 + 4); % s (3sigma Killestein + 4 s GPS def. uncertainty)

a_fact = porb_max/(2.0*pi*299792458);

a_min = 1.44

a_max = 3.25
@@ -106,13 +108,15 @@ for fifo = 1:n_fr_intervals 
        t_raw = fitsread(filoc,"binarytable");

        t_raw = t_raw{1};

        % tstart = str2double(fits.readKey(file,'TSTART'));

        tstart = t_raw(1);

        fits.closeFile(file);

        nbin=ceil((t_raw(end)-t_raw(1))/dt_psd);

        [C,t]=(histcounts(t_raw,nbin));

        C=C.'; %conteggi

        t=(t(1:end-1)+(t(2)-t(1))/2).'; % vettore tempi rebinnato, prendo il centro del bin

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

        halfspli = 0.5*rem(nbin*dt,(t_raw(end)-t_raw(1)));

        tstart = t_raw(1)-halfspli;

        tedgeswhole = (tstart:dt:tstart+nbin*dt);

        C=(fastcounts(t_raw,tedgeswhole));

        %C=C.'; %conteggi

        t=(tedgeswhole(1:end-1)+dt/2).'; % vettore tempi rebinnato, prendo il centro del bin

        M=fix((nbin*dt)/Tseg);

        tm = zeros(M,N);

        tmid = zeros(M,1);

        for m = 1:M
@@ -193,13 +197,15 @@ for fifo = 1:n_fr_intervals 
            t_raw = fitsread(filoc,"binarytable");

            t_raw = t_raw{1};

            % tstart = str2double(fits.readKey(file,'TSTART'));

            tstart = t_raw(1);

            fits.closeFile(file);

            nbin=ceil((t_raw(end)-t_raw(1))/dt_psd);

            [C,t]=(histcounts(t_raw,nbin));

            C=C.'; %conteggi

            t=(t(1:end-1)+(t(2)-t(1))/2).'; % vettore tempi rebinnato, prendo il centro del bin

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

            halfspli = 0.5*rem(nbin*dt,(t_raw(end)-t_raw(1)));

            tstart = t_raw(1)-halfspli;

            tedgeswhole = (tstart:dt:tstart+nbin*dt);

            C=(fastcounts(t_raw,tedgeswhole));

            %C=C.'; %conteggi

            t=(tedgeswhole(1:end-1)+dt/2).'; % vettore tempi rebinnato, prendo il centro del bin

            M=fix((nbin*dt)/Tseg);

            tm = zeros(M,N);

            tmid = zeros(M,1);

            for m = 1:M
@@ -391,8 +397,9 @@ for fifo = 1:n_fr_intervals 
    

            graph1 = [pathfigu,'PS_',finame,'_Tseg',num2str(Tseg),'_m',num2str(m),'.pdf'];

            graph2 = [pathfigu,'CT_',finame,'_Tseg',num2str(Tseg),'_m',num2str(m),'.pdf'];

            graph3 = [pathfigu,'PD_',finame,'_Tseg',num2str(Tseg),'_m',num2str(m),'.pdf'];

            if (~(isfile(graph1)) || ~(isfile(graph2)))

                [VsuEtot,epsct,hmp] = powdist(Nyq,ttemp,false,graph1,graph2);

                [VsuEtot,epsct,hmp] = powdist(Nyq,ttemp,false,graph1,graph2,graph3);

            end

            norman = gather(2.*abs(Y0(1))./mean((abs(Y0(F1>1000 & F1<Nyq)).^2)));

            tau = ones(size(ttemp),"like",ttemp);