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

Added multifrequency

parent 96379c7b

functions/powdist.m

+1 −0
Original line number Diff line number Diff line
@@ -191,6 +191,7 @@ t = tiledlayout(2,1); 
nexttile

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 distribution'},'Location','southwest')

text(0.8,0.870,['$\epsilon_{ct} = $',num2str(epscr)])

xlabel('Counts')

ylabel('Number of bins with x counts')

title('Extrapolated generalized Poisson distributions vs data')

scsearch_1023_multiobs.m

+468 −393
Original line number Diff line number Diff line
@@ -4,9 +4,10 @@ close all; 
import matlab.io.*

addpath('functions/')



maxNumCompThreads(48);

maxNumCompThreads(96);



pathfi = '/data/Sorgenti/J1023/03_2016/';

% pathfi = '/data/Sorgenti/J1023/03_2016/';

pathfi = '/data/Sorgenti/J1023/30-31_01_2020/';

folname = [pathfi,'SP_search_',char(datetime('now','Format','dd_MM'))];

figfolname = [folname,'/figures'];

if (exist(folname,"dir")==7)
@@ -48,12 +49,12 @@ toc 
pathfi = [folname,'/']



Nyq = 2000 %cambiare in caso

Tseg = 256

Tseg = 128

dt_psd=1/(2*Nyq);

dt = dt_psd;

N = fix(Tseg/dt);

f_min = 550

f_max = 600

f_supamin = 49

f_supamax = 1550

porb_min = 17115.3216 % s

porb_max = 17115.3216; porb_gr = (porb_max).'; % s

a_fact = porb_max/(2.0*pi*299792458);
@@ -62,24 +63,38 @@ a_max = a_fact*280000*0.72 
% tasc_min = (58878.9911009 - MJDREF).*86400 - porb_max/2

% tasc_max = (58878.9911009 - MJDREF).*86400 + porb_max/2

f_step = 1/Tseg

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

lesupf = (f_supamax-f_supamin)/f_step + 1;

totfgr = (f_supamin:f_step:f_supamax).';

bruttone = zeros(lesupf,1,'single');



% tasc_mid = (tasc_max + tasc_min)/2

% tasc_step = tasc_max-tasc_min;

% porb_step = porb_max;

% a_step = a_max;

tol = 0.9



for fifo = 1:15

    f_min = f_supamin + (100*(fifo-1))

    f_max = f_supamax - 1400 + (100*(fifo-1))

    MJDREF = zeros(numfiles,1);

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

    lefgr = length(f_gr);

    

    tic

    %% First step in par., def. also max min tasc

    Mold = 0;

    finame = char(mario(1))

    filoc = [pathfi,'../',finame];

    countsname = ['counts_',finame,'_Tseg_',num2str(Tseg),'.mat'];

    if(exist([pathfi,countsname],"file")==2)

        load([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","tstart","nbin","mjdref");

        % clear t_raw

    else

        file=fits.openFile(filoc);

        mjdref=fits.readKey(file,'MJDREF');

MJDREF(1)=str2double(mjdref);

tasc_min = (57449.72579 - MJDREF(1)).*86400 - porb_max/2

tasc_max = (57449.72579 - MJDREF(1)).*86400 + porb_max/2

tasc_mid = (tasc_max + tasc_min)/2

        %tasc_min = (57449.72579 - MJDREF(1)).*86400 - porb_max/2

        %tasc_max = (57449.72579 - MJDREF(1)).*86400 + porb_max/2

        %tasc_mid = (tasc_max + tasc_min)/2

        fits.movAbsHDU(file,2);

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

        t_raw = fitsread(filoc,"binarytable");
@@ -90,8 +105,6 @@ nbin=ceil((t_raw(end)-t_raw(1))/dt_psd); 
        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);

Mtot = M;

Mold = 0;

        tm = zeros(M,N);

        tmid = zeros(M,1);

        for m = 1:M
@@ -108,11 +121,34 @@ aux1 = M*N; 
        aux2 = C(1:aux1); clear aux1

        x = reshape(aux2,[N,M]); clear aux2

        clear C

countsname = ['counts_',finame,'_Tseg_',num2str(Tseg),'.mat'];

save([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","tstart","nbin","-v7.3","-nocompression");

        save([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","tstart","nbin","mjdref","-v7.3","-nocompression");

    end

    

    Mtot = M;

    MJDREF(1)=str2double(mjdref);

    % Norb=round((MJDREF(1)-Tasc_old)/(Porb/86400.0));

    % Tasc_propagato = Tasc_old+Norb*(Porb/86400.0)

    % porb_unc =  0.0000036*86400.0;

    % Tasc_old_unc = 0.0053*86400.0; % in sec

    % Tasc_new_unc = sqrt(Tasc_old_unc^2 + (Norb*porb_unc)^2); % s 

    % tasc_mid = (Tasc_propagato - MJDREF(1))*86400

    % tasc_min = tasc_mid - 2.0*Tasc_new_unc

    % tasc_max = tasc_mid + 2.0*Tasc_new_unc

    tasc_min = (57449.72579 - MJDREF(1)).*86400 - porb_max/2

    tasc_max = (57449.72579 - MJDREF(1)).*86400 + porb_max/2

    tasc_mid = (tasc_max + tasc_min)/2

    

    toff = tstart - tasc_mid;

% delporb(a_max,porb_min,f_max,toff,Tseg,M,tol);

    % a_min = 0.008

    % a_max = 0.7

    % porb_max = Porb + 3*porb_unc

    % porb_min = Porb - 3*porb_unc

    %%stepparino = delporb(a_max,porb_min,f_max,toff,Tseg,M,tol);

    %%if (((porb_max-porb_min)/stepparino) > 0.5)

    %%    porb_step = (porb_max-porb_min)/ceil((porb_max-porb_min)/stepparino)

    %%else

    %%    porb_step = stepparino

    %%end

    % (porb_min^2)*sqrt(0.1*(1-tol))/(a_max*2*pi*Tseg*f_max);

    % % stepparino = max(delporb(a_max,porb_min,f_max,toff,Tseg,M,tol),(porb_min^2)*sqrt(0.1*(1-tol))/(a_max*2*pi*Tseg*f_max))

    % % porb_step = (porb_max-porb_min)/ceil((porb_max-porb_min)/stepparino)
@@ -138,9 +174,13 @@ for o = 2:numfiles 
        %% WILL have to choose whether to do only the parameter space check here

        

        filoc = [pathfi,'../',finame];

        countsname = ['counts_',finame,'_Tseg_',num2str(Tseg),'.mat'];

        if(exist([pathfi,countsname],"file")==2)

            load([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","tstart","nbin","mjdref");

            % clear t_raw

        else

            file=fits.openFile(filoc);

            mjdref=fits.readKey(file,'MJDREF');

    MJDREF(o)=str2double(mjdref);

            fits.movAbsHDU(file,2);

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

            t_raw = fitsread(filoc,"binarytable");
@@ -151,7 +191,6 @@ for o = 2:numfiles 
            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);

    Mtot = Mtot + M;

            tm = zeros(M,N);

            tmid = zeros(M,1);

            for m = 1:M
@@ -168,11 +207,16 @@ for o = 2:numfiles 
            aux2 = C(1:aux1); clear aux1

            x = reshape(aux2,[N,M]); clear aux2

            clear C

    countsname = ['counts_',finame,'_Tseg_',num2str(Tseg),'.mat'];

    save([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","tstart","nbin","-v7.3","-nocompression");



            save([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","tstart","nbin","mjdref","-v7.3","-nocompression");

        end

        MJDREF(o)=str2double(mjdref);

        Mtot = Mtot + M;

        toff = tstart - tasc_mid;

    delporb(a_max,porb_min,f_max,toff,Tseg,M,tol);

        %%stepparino = delporb(a_max,porb_min,f_max,toff,Tseg,M,tol);

        %%stepparino = (porb_max-porb_min)/ceil((porb_max-porb_min)/stepparino)

        %%if ((stepparino < porb_step) && (stepparino>0.5*(porb_max-porb_min)))

        %%    porb_step = stepparino

        %%end

        % (porb_min^2)*sqrt(0.1*(1-tol))/(a_max*2*pi*Tseg*f_max)

        % % stepparino = max(delporb(a_max,porb_min,f_max,toff,Tseg,M,tol),(porb_min^2)*sqrt(0.1*(1-tol))/(a_max*2*pi*Tseg*f_max))

        % % porb_step = (porb_max-porb_min)/ceil((porb_max-porb_min)/stepparino)
@@ -200,12 +244,17 @@ end 
    %% Finally make grid of orb. parameters with smallest steps (and define trupars)

    a_gr = (a_min:a_step:a_max).';

    tasc_gr = (tasc_min:tasc_step:tasc_max).';

% porb_gr = (porb_min:porb_step:porb_max).';

    %%if (porb_step < 6*porb_unc)

    %%    porb_gr = Porb;

    %%else

    %%    porb_gr = (porb_min:porb_step:porb_max).';

    %%end

    

    trupar = [5.924214682724856e+02,1.71155216592e+04,0.343356,58878.9911009];

    toc

    

    %% Just getting an idea

nodesnum = length(porb_gr)*length(a_gr)*length(tasc_gr)*length(f_gr)

    nodesnum = length(porb_gr)*length(a_gr)*length(tasc_gr)*lefgr

    

    %% Aux. param. and coherent grid generation

    % Aux par
@@ -270,10 +319,10 @@ s_s 
    Mtot

    % Lambda = zeros(length(f_gr),length(nibank),'single','gpuArray'); % MATLAB

    % stores array elements in a column-major layout

Lambda = ones(length(f_gr),length(nibank),'single');

    Lambda = ones(lefgr,length(nibank),'single');

    F1=gpuArray((0:N-1)./(N*dt_psd)).';

f_ind = zeros(1,length(f_gr),'uint32');

for n = 1:length(f_gr)

    f_ind = zeros(1,lefgr,'uint32');

    for n = 1:lefgr

        [~,f_ind(n)] = min(abs(F1-f_gr(n)));

    end
@@ -282,7 +331,7 @@ for o = 1:numfiles 
            finame = char(mario(o))

            countsname = ['counts_',finame,'_Tseg_',num2str(Tseg),'.mat'];

            load([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","tstart","nbin");

        Lambda = ones(length(f_gr),length(nibank),'single');

            Lambda = ones(lefgr,length(nibank),'single');

        end

        tm=gpuArray(tm);

        tedge = ones(1,N+1,'gpuArray');
@@ -341,6 +390,8 @@ for o = 1:numfiles 
    

            end

            Lambda = (Lambda.^2).*norman;

            mala = double(max(Lambda,[],'all'));

            promala = 1-(chi2cdf(mala,2))^(nino*lefgr);

    

            lamfiname = ['Lambda_',finame,'_fmax_',num2str(f_max),'_seg_',num2str(m),'.mat'];

            save([pathfi,lamfiname],"Lambda","-v7.3","-nocompression");
@@ -349,6 +400,10 @@ for o = 1:numfiles 
            segend = toc(segstart);

            avetime = avetime + segend;

            disp(['Just finished segment num.: ',num2str(m),'.']);

            disp(['The maximum Lambda for this segment (',num2str(mala),') has probability ',num2str(promala),' to arise from noise.']);

            if (promala<0.0028)

                disp(['---!!!---Go check it out!---!!!---']);

            end

            disp(['Time elapsed in this loop = ',num2str(avetime),' s, estimated time left = ',num2str(((M-m)*avetime/m)),' s.']);  disp(' ');

    

        end
@@ -356,10 +411,10 @@ for o = 1:numfiles 
    

        %% The first time in this loop we need to set some general variables such as niwalk and totlam 

        if (o == 1)

        nodesnum = length(porb_gr)*length(a_gr)*length(tasc_gr)*length(f_gr)

            nodesnum = length(porb_gr)*length(a_gr)*length(tasc_gr)*lefgr

            nodesmem = (nodesnum*8.0)/1024^3; %->to bytes -> to GBs

    

        memthres = 600.0 - (length(f_gr)*nino*8.0)/1024^3;

            memthres = 600.0 - (lefgr*nino*8.0)/1024^3;

            if (nodesmem > memthres)

                % THIS DISCLAIMER IS STILL HERE AFTER MONTHS

                % FIX IIIIT
@@ -380,7 +435,7 @@ for o = 1:numfiles 
            end

            % Writing on a matrix of ones is faster in MATLAB than on one made 

            % of zeroes, we shall subtract one after all calculations are done

        totlam = ones(length(f_gr),parno,'single');

            totlam = ones(lefgr,parno,'single');

        end

    

        if (o == 1)
@@ -389,7 +444,7 @@ for o = 1:numfiles 
            mjdiff = (MJDREF(o) - MJDREF(1))*86400.0;

        end

    

    guess = parno*M*16192.4539/(4.0*7546032.0);

        guess = parno*((nino/66560)^1.3)*6.0*(M/7.0);

        disp(['(Probably terrible) Estimate of the time needed for the second loop in the current file (',num2str(o),' of ',num2str(numfiles),') =']);

        if (guess > 60)

            guess = guess/60;
@@ -450,14 +505,14 @@ end 
    totlam = totlam - 1;

    %% Calculate detection threshold for a multi-trial false alarm prob. = pfa

    pfa = 0.01;

sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(parno*length(f_gr)))),Mtot,'upper');

%sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(sum(nimask)*length(f_gr)))),M,'upper');

    sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(parno*lesupf))),Mtot,'upper');

    %sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(sum(nimask)*lefgr))),M,'upper');

    [bru,tto] = max(totlam,[],2);

    [maxtotlam,maxinde] = max(bru);

    bestpar = [f_gr(maxinde), parbank(tto(maxinde),:), double(maxtotlam)]; 

    bparch = string(bestpar);

    tascmjd = MJDREF(1)+tasc_gr./86400;

pfamax = 1.0-( (gammainc(0.5*double(maxtotlam),Mtot,'lower'))^(parno*length(f_gr)));

    pfamax = 1.0-( (gammainc(0.5*double(maxtotlam),Mtot,'lower'))^(parno*lesupf));

    disp('Our most significant peak has parameters:');

    disp(['f_spin = '+bparch(1)+' Hz, p_orb = '+bparch(2)+' s, a*sin(i)/c = '+bparch(3)+' s, t_asc = '+bparch(4)+' s.']);

    disp(['t_asc in days is = ',num2str((MJDREF(1)+bestpar(4)/86400)),'.']);
@@ -465,8 +520,8 @@ disp(['Its detection statistics is Σ_max = ',num2str(maxtotlam),',']); 
    disp(['which considering all parameters (',num2str(parno),') corresponds to a false alarm probability of ',num2str(pfamax),','])

    disp(['while the multi-trial threshold for a false alarm probability = ',num2str(pfa),' is eq. to ',num2str(sigmastar),'.']);

    truwalk = size(unique(niwalk.', 'rows'),1);

sigmastar2 = 2*gammaincinv(((1-pfa)^(1/(truwalk*length(f_gr)))),Mtot,'lower');

pfamax2 = 1.0-( (gammainc(0.5*double(maxtotlam),Mtot,'lower'))^(truwalk*length(f_gr)));

    sigmastar2 = 2*gammaincinv(((1-pfa)^(1/(truwalk*lesupf))),Mtot,'lower');

    pfamax2 = 1.0-( (gammainc(0.5*double(maxtotlam),Mtot,'lower'))^(truwalk*lesupf));

    disp(['The number of combinations of ni combinations actually walked is ',num2str(truwalk),',']);

    % sigmastar2 = 2*gammaincinv((1-(1-pfa)^(1/(sum(nimask)*length(f_gr)))),M,'upper');

    % pfamax2 = 1.0-( (1.0-gammainc(0.5*double(maxtotlam),M,'upper'))^(sum(nimask)*length(f_gr)));
@@ -489,7 +544,7 @@ inGB = summina/1024^3; 
    %fprintf("The workspace will need %.5f KB or %.5f MB or %.5f GB of memory", inKB, inMB, inGB)

    

    timerapp = 24.0*summina/518260001575;



    bruttone(1+(fifo-1)*(lefgr):(fifo)*(lefgr)) = bru;    

    % disp(bestpar);

    % save('C:\Users\Filippo\Desktop\XMM_Jxxx\risultelli.mat');

    disp(['I am about to save around ',num2str(round(inGB)),' GB of data (+ headers), it should take roughly ',num2str(ceil(timerapp)),' minutes.']);
@@ -498,6 +553,26 @@ disp(['Results saved in ',pathfi,'risultelli_tot_',char(datetime('now','Format', 
    %% 

    



    figu = figure('visible','on');

    % figu = figure('visible','on');

    rina = [pathfigu,'f_sig_',num2str(f_max),'.m'];

    plot(f_gr,bru)

    yline(sigmastar,'--','Σ*')

    yline(sigmastar2,'--','Σ*2')

    saveas(figu,rina); clear figu



end





figu = figure('visible','on');

% figu = figure('visible','on');

rina = [pathfigu,'f_sig_tot.m'];

plot(totfgr,bruttone)

yline(sigmastar,'--','Σ*')

yline(sigmastar2,'--','Σ*2')

saveas(figu,rina); clear figu





daterella = ['Date and time now ',char(datetime('now')),' (UTC), maybe'];

disp(daterella); clear daterella
@@ -558,7 +633,7 @@ function res = delporb(asenc,porb,fspin,toff,tseg,M,eqtol) 
    res = 10.0^10;

    for m = 1:M

        fun = @(step) ((sum(sin(omespin.*pert((0:nrot-1),step,m))).^2 + sum(cos(omespin.*pert((0:nrot-1),step,m))).^2)./(nrot).^2.0 - eqtol);

        root = fzero(fun,[1 porb]);

        root = fzero(fun,[0.01 porb]);

        if (root<res)

            res = root;

        end