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

Generalized a bit 1023, added background to SCO (latest)

parent c20554ab

scsearch_1023_multiobs.m

+51 −25
Original line number Diff line number Diff line
@@ -29,8 +29,8 @@ pathfi = [folname,'/']; 
pathfigu = [figfolname,'/'];





reset(gpuDevice(1));

gpuDevice(1);

reset(gpuDevice(2));

gpuDevice(2);



diary([pathfi,'loggerino_metrica_',char(datetime('now','Format','dd_MM')),'.log'])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -48,13 +48,17 @@ toc 


pathfi = [folname,'/']



Nyq = 2000 %cambiare in caso

Tseg = 128

dt_psd=1/(2*Nyq);

Nyqmin = 2000; %cambiare in caso

Tsegmin = 128;

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

f_supamax = 1550

f_supamin = 49 % Hz

f_supamax = 1550 % Hz

fr_int_width = 100 % Hz (-1)

n_fr_intervals = (f_supamax-f_supamin-1)/fr_int_width;

porb_min = 17115.3216 % s

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

a_fact = porb_max/(2.0*pi*299792458);
@@ -73,9 +77,9 @@ bruttone = zeros(lesupf,1,'single'); 
% 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))

for fifo = 1:n_fr_intervals

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

    f_max = f_supamax - fr_int_width*(n_fr_intervals - fifo)

    MJDREF = zeros(numfiles,1);

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

    lefgr = length(f_gr);
@@ -96,15 +100,18 @@ for fifo = 1:15 
        %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");

        t_raw = t_raw{1};

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

        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
@@ -124,6 +131,7 @@ for fifo = 1:15 
        save([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","tstart","nbin","mjdref","-v7.3","-nocompression");

    end

    

    Nphotot = sum(x);

    Mtot = M;

    MJDREF(1)=str2double(mjdref);

    % Norb=round((MJDREF(1)-Tasc_old)/(Porb/86400.0));
@@ -182,15 +190,18 @@ for fifo = 1:15 
            file=fits.openFile(filoc);

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

            fits.movAbsHDU(file,2);

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

            t_raw = fitsread(filoc,"binarytable");

            t_raw = t_raw{1};

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

            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
@@ -210,6 +221,7 @@ for fifo = 1:15 
            save([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","tstart","nbin","mjdref","-v7.3","-nocompression");

        end

        MJDREF(o)=str2double(mjdref);

        Nphotot = Nphotot + sum(x,"all");

        Mtot = Mtot + M;

        toff = tstart - tasc_mid;

        %%stepparino = delporb(a_max,porb_min,f_max,toff,Tseg,M,tol);
@@ -369,8 +381,9 @@ for fifo = 1:15 
    

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

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

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

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

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

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

                [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);
@@ -391,7 +404,7 @@ for fifo = 1:15 
            end

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

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

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

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

    

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

            save([pathfi,lamfiname],"Lambda","-v7.3","-nocompression");
@@ -544,7 +557,7 @@ for fifo = 1:15 
    %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;    

    bruttone((1+(fifo-1)*(fr_int_width/f_step)):((fifo*fr_int_width+1)/f_step)+1) = 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.']);
@@ -563,6 +576,19 @@ for fifo = 1:15 


end



[maxoni,maxindi] = maxk(bruttone,20);

if (maxoni(1) > sigmastar2) 

    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')

    disp(['A signal was found (total summed power ',num2str(maxoni(1)),') at a frequency = ',num2str(totfgr(maxindi(1)),12),' Hz']);

    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')

else

    pup = power_upper_limit(maxoni,Mtot,0.0027);

    amplicost = 1.61*(Nphotot)^(-0.5); %No background for now

    Auls = amplicost.*sqrt(pup)./(pi.*sinc(totfgr(maxindi)./2.*Nyq));

    [maxau,maxauind] = max(Auls);

    disp(['No signal was found, our UL on amplitude is ',num2str(maxau),' (p_ul at 3 sig) at a frequency = ',num2str(totfgr(maxindi(maxauind))),' Hz (no background)']);

end





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

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

scsearch_SCO_multiobs_multifreq.m

+25 −9
Original line number Diff line number Diff line
@@ -29,8 +29,8 @@ pathfi = [folname,'/']; 
pathfigu = [figfolname,'/'];





reset(gpuDevice(2));

gpuDevice(2);

reset(gpuDevice(1));

gpuDevice(1);



diary([pathfi,'loggerino_metrica_',char(datetime('now','Format','dd_MM')),'.log'])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -75,6 +75,13 @@ f_step = 1/Tseg 
lesupf = (f_supamax-f_supamin)/f_step + 1;

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

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

% BKG rates 2022/2023

backrates = [6354.0,7060.2,15510.6,6203.6];

% BKG rates 2018/2019

% backrates = [28929.8,15235.6,11246.5,12505.8,13026.6,16324.0,6496.1,12180.5,12072.7];



back = zeros(size(backrates));

Nphotot = zeros(size(backrates));



% tasc_mid = (tasc_max + tasc_min)/2

% tasc_step = tasc_max-tasc_min;
@@ -137,6 +144,7 @@ for fifo = 1:n_fr_intervals 
    end

    

    Nphotot(1) = sum(x,"all");

    back(1) = M*Tseg*backrates(1);

    Mtot = M;

    MJDREF(1)=str2double(mjdref);

    Norb=round((abs(MJDREF(1)-Tasc_old))/(Porb/86400.0));
@@ -227,6 +235,7 @@ for fifo = 1:n_fr_intervals 
        end

        MJDREF(o)=str2double(mjdref);

        Nphotot(o) = sum(x,"all");

        back(o) = M*Tseg*backrates(o);

        Mtot = Mtot + M;

        if (round(abs((MJDREF(o)-Tasc_old))/(Porb/86400.0)) > Norb)

            Norb = round(abs((MJDREF(o)-Tasc_old))/(Porb/86400.0));
@@ -307,9 +316,11 @@ for fifo = 1:n_fr_intervals 
    

    % Try s* and check ν^s range

    tic

    g_jj=(((pi*Tseg)^2)/3).*[1; (Tseg^2)/60; (Tseg^4)/1344; (Tseg^6)/172800]; %eq. 22 M2015 + calcoli da eq. 21 M2015

    mu_s=0.05;

    s_s = 4;

    while(1)

        if((nismax(s_s)-nismin(s_s))<delta_ni_new(s_s,Tseg)/f_max)

        if((nismax(s_s)-nismin(s_s))<delta_ni(mu_s,s_s,g_jj(s_s))/f_max)

            s_s=s_s-1;

        else

            break;
@@ -327,7 +338,7 @@ for fifo = 1:n_fr_intervals 
    % element is accessed as nis{s}(i)

    Nni = zeros(s_s,1);

    for s = 1:s_s

        Nni(s) = ceil((nismax(s)-nismin(s))*f_max/delta_ni_new(s,Tseg));

        Nni(s) = ceil((nismax(s)-nismin(s))*f_max/delta_ni(mu_s,s_s,g_jj(s)));

        if (mod(Nni(s),2)==0)

            Nni(s) = Nni(s)+1;

        end
@@ -335,6 +346,7 @@ for fifo = 1:n_fr_intervals 
        nis{s}=franco;

    end

    infax = gpuArray(1./factorial(1:s_s)); % MATLAB's factorial is slow

    lenin1 = length(nis{1});



    toc
@@ -421,7 +433,7 @@ for fifo = 1:n_fr_intervals 
            end

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

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

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

            promala = 1-(chi2cdf(mala,2))^(nino*(max(lesupf/lenin1,n_fr_intervals))*Mtot);

    

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

            save([pathfi,lamfiname],"Lambda","-v7.3","-nocompression");
@@ -601,13 +613,17 @@ if (maxoni(1) > sigmastar2) 
else

    pup = power_upper_limit(maxoni,Mtot,0.0027);

    if (exist('back','var')==1)

        amplicost = 1.61*(sum(Nphotot(o)/(Tseg*((Nphot/Tseg) + 15000)))*(Nphot*51)^(-0.5);

        amp = 0;

        for o = 1:length(mario)

            amp = amp + Nphotot(o) - back(o);

        end

        amplicost = 1.61*(sqrt(sum(Nphotot(o)))/amp);

        Auls = amplicost.*sqrt(pup)./(pi.*sinc(totfgr(maxindi)./2.*Nyq));

        [maxau,maxauind] = max(Auls);

        disp(['No signal was found, our UL on amplitude is ',num2str(maxau),' (p_ul at 3 sig) at a frequency = ',num2str(totfgr(maxindi(maxauind))),' Hz (no background)']);

        disp(['No signal was found, our UL on amplitude is ',num2str(maxau),' (p_ul at 3 sig) at a frequency = ',num2str(totfgr(maxindi(maxauind))),' Hz (back-subtracted)']);



    else

        amplicost = 1.61*(Nphotot)^(-0.5); %No background for now

        amplicost = 1.61*(sum(Nphotot(o)))^(-0.5); %No background for now

        Auls = amplicost.*sqrt(pup)./(pi.*sinc(totfgr(maxindi)./2.*Nyq));

        [maxau,maxauind] = max(Auls);

        disp(['No signal was found, our UL on amplitude is ',num2str(maxau),' (p_ul at 3 sig) at a frequency = ',num2str(totfgr(maxindi(maxauind))),' Hz (no background)']);