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

Two files deleted 

Added new way of properly counting number of trials (still a slight overestimate, but lower)
parent 11c37c47

J1023_B_2017_Bary.fits

deleted100644 → 0
−111 MiB

File deleted.

View file

scsearch.m

+57 −69
Original line number Diff line number Diff line
@@ -3,10 +3,12 @@ close all; 


import matlab.io.*



pathfi = '/data/Sorgenti/3FGLJ1544.6-1125/';

% pathfi = '/data/Sorgenti/SCORPIUS-X-1/';

maxNumCompThreads(48);



% pathfi = '/data/Sorgenti/3FGLJ1544.6-1125/';

pathfi = '/data/Sorgenti/SCORPIUS-X-1/';

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

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

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

if (exist(folname,"dir")==7)

    disp(['Folder ',folname,' already exists, previous Lambda files (if present) may be lost']);

else
@@ -25,10 +27,11 @@ disp(daterella); clear daterella 


tic



finame = 'SiFAP2_20220429_3FGLJ1544_N1_Bary.fits';

% finame = 'ScoX1_20230423_Bary.fits';

% finame = 'J1023_B_2017_Bary.fits';

% finame = 'EPN_0744840201_bary.fits';

% finame = 'SiFAP2_20220429_3FGLJ1544_N1_Bary.fits'

% finame = 'ScoX1_20230423_Bary.fits'

% finame = 'J1023_B_2017_Bary.fits'

% finame = 'EPN_0744840201_bary.fits'

finame = 'ScoX1_20220628_N2_Targ_Bary.fits'



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

file=fits.openFile(filoc);
@@ -39,10 +42,7 @@ tstart = str2double(fits.readKey(file,'TSTART')); 
fits.closeFile(file);



t_raw = fitsread(filoc,"binarytable");

%t_raw = fitsread('C:\Users\Filippo\Desktop\J1544_2022\SiFAP2_20220429_3FGLJ1544_N1_Bary.fits','binarytable');

t_raw = t_raw{1};

% MJDREF 1544

% MJDREF=59698.943055555602768;





% info = fitsinfo(finame).BinaryTable.Keywords;
@@ -63,7 +63,7 @@ pathfi = [folname,'/']; 
t0 = (t_raw(end)+t_raw(1))/2;

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

Nyq = 2000 %cambiare in caso

Tseg = 490

Tseg = 512



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

nbin=round((t_raw(end)-t_raw(1))/dt_psd);
@@ -168,26 +168,26 @@ end 


% Parameters for Sco X-1

 % s 

f_max = 800

f_min = 750

f_max = 750

f_min = 700

f_step = 1/Tseg

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

% Porb = 68023.919 % s

% porb_gr = (Porb).'; % s

% porb_unc = 0.017; % s

Porb = 68023.919 % s

porb_gr = (Porb).'; porb_max = Porb; porb_min = Porb; % s

porb_unc = 2*0.017; % s

% porb_max = 28800

% porb_min = 14400

% porb_step = 60

% porb_gr = (porb_min:porb_step:porb_max);

% Tasc_old = 56722.63037153067129629630; % MJD

Tasc_old = 57508.3323; % MJD

Porb = 0.2415361*86400.0

Tasc_old = 56722.6303715306713; % MJD

% Tasc_old = 57508.3323; % MJD

% Porb = 0.2415361*86400.0

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

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

Tasc_propagato = Tasc_old+Norb*(Porb/86400.0)

porb_unc =  0.0000036*86400.0;

% Tasc_old_unc = 37; % in sec

Tasc_old_unc = 0.0053*86400.0; % in sec

Tasc_old_unc = 2*33 + 4; % in sec

% 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)*86400

tasc_min = tasc_mid - 2.0*Tasc_new_unc
@@ -197,18 +197,18 @@ tasc_max = tasc_mid + 2.0*Tasc_new_unc 
% tasc_gr = (((58106.9494:8e-4:58107.1595) - MJDREF).*86400).';

% tasc_min = min(tasc_gr)

% tasc_max = max(tasc_gr)

tasc_mid = 0.5*(tasc_max+tasc_min)

a_min = 0.008

a_max = 0.15

% tasc_mid = 0.5*(tasc_max+tasc_min)

a_min = 1.44

a_max = 3.25

nrot = ceil(f_max*Tseg);

toff = tstart - tasc_mid;

porb_max = 21600

porb_min = 18000

% porb_max = 21600

% porb_min = 18000

% porb_step = 60

tol = 0.9

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)

porb_gr = (porb_min:porb_step:porb_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)

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

stepparino = delasenc(porb_min,f_max,toff,Tseg,M,tol)

a_step = (a_max-a_min)/ceil((a_max-a_min)/stepparino)

% a_step = 1.0/(2.0*f_max*0.9)
@@ -368,7 +368,7 @@ nino = length(nibank); 
s_s

M

avetime = 0.0;

guess = nino*M*5614.2461/(5376.0*48);

guess = nino*M*3071.5583/(3100.0*28);

disp(['Rough estimate of the time necessary for the first loop = ',num2str(guess),'']);  disp(' ');

for m=1:M
@@ -389,7 +389,7 @@ for m=1:M 
    % ttemp = gpuArray(tm(m,:));

    % xtemp = gpuArray(x(:,m).');

    % ttemp = tm(m,:);

    ttemp = gpuArray(t_raw((t_raw>=tm(m,1)-0.5*dt) & t_raw<=tm(m,N)+0.5*dt));

    ttemp = gpuArray(t_raw((t_raw>=gather(tm(m,1))-0.5*dt) & t_raw<=gather(tm(m,N))+0.5*dt));

    Nphot = length(ttemp);

    ttempdif = (ttemp(:)-tmid(m)).';

    ttempdifl(1:N) = (tm(m,1:N)-tmid(m)-0.5*dt).';
@@ -451,7 +451,7 @@ for m=1:M 
    % disp(length(f_gr)*i);

    % disp(m);



    lamfiname = sprintf('Lambda_seg_%d.mat', m);

    lamfiname = sprintf('Lambda_fmax_%d_seg_%d.mat', f_max, m);

    save([pathfi,lamfiname],"Lambda","-v7.3","-nocompression");

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

    % Lambda = zeros(length(f_gr),length(nibank),'single','gpuArray');
@@ -480,47 +480,29 @@ parbank = combinations(porb_gr,a_gr,tasc_gr).Variables; 
%bestpar = zeros(1,5);

toc



parno = length(parbank);

if (nino<intmax('uint16'))

    niwalk = ones(M,parno,'uint16');

elseif (nino<intmax('uint32'))

    niwalk = ones(M,parno,'uint32');

elseif (nino<intmax('uint64'))

    niwalk = ones(M,parno,'uint64');

else

    disp(['The length of nibank in greater than ',num2str(intmax('uint64')),', something definitely makes no sense here'])

end



disp(['length(parbank) = ',num2str(length(parbank)),'']);

nimask = false(length(nibank),1);

% totlam = zeros(length(f_gr),length(parbank),'single');

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

guess = length(parbank)*M*29761.2957/(48.0*2180530.0);

disp(['Probably terrible estimate of the time necessary for the second loop = ',num2str(guess),'']);  disp(' ');

% tic

% for n=1:length(f_gr)

%     for i=1:length(parbank)

%         curpar = [f_gr(n),parbank(i,:)];

%         % Move from ν,P,a,T_asc to ν,Ω,a,γ

%         curpar(2) = 2*pi/curpar(2);

%         for m = 1:M

%             curpar(4) = curpar(2)*(tmid(m) - parbank(i,3));

%             curni=zeros(1,s_s);

%             for s=1:s_s

%                 curni(s) = (curpar(2)^s)*sin(curpar(4)+0.5*s*pi);

%             end

%             curni = -curni.*curpar(3);

%             curni(1) = curni(1) + 1;

% 

%             lamfiname = sprintf('Lambda_seg_%d.mat', m);

%             % load(lamfiname);

%             load([pathfi,lamfiname]);

% 

%             [Idx,D] = knnsearch(nisearcher,curni);

%             totlam(i,n) = totlam(i,n)+Lambda(Idx,n);

% 

%             clear Lambda

% 

%         end

%         % % if (totlam > bestpar(1))

%         % %     bestpar = [totlam,f_gr(n),parbank(i,:)];

%         % % end

%     end

% end

disp(['(Probably terrible) Estimate of the time necessary for the second loop = ',num2str(guess),'']);  disp(' ');



bles = cospi(1/2);

avetime = 0.0;

for m = 1:M

    segstart = tic;

    lamfiname = sprintf('Lambda_seg_%d.mat', m);

    lamfiname = sprintf('Lambda_fmax_%d_seg_%d.mat', f_max, m);

    load([pathfi,lamfiname]);

    % toc

    % disp(m);
@@ -537,10 +519,13 @@ for m = 1:M 
        end

        curni = -curni.*curpar(3);

        curni(1) = curni(1) + 1;



        [Idx,D] = knnsearch(nisearcher,curni);

        % if (~nimask(Idx))

        nimask(Idx) = 1;

        % end

        niwalk(m,i) = Idx;



        totlam(:,i) = totlam(:,i)+Lambda(:,Idx);



    end
@@ -600,8 +585,11 @@ disp(['t_asc in days is = ',num2str((MJDREF+bestpar(4)/86400)),'.']); 
disp(['Its detection statistics is Σ_max = ',num2str(maxtotlam),',']);

disp(['which considering all parameters, 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),'.']);

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

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

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

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

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

disp(['When considering only the nis walked instead, Σ_max corresponds to a false alarm probability of ',num2str(pfamax2),','])

disp(['while the multi-trial threshold for a false alarm probability = ',num2str(pfa),' is eq. to ',num2str(sigmastar2),'.']);
@@ -620,8 +608,8 @@ diary off 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%% Fai il grafico delle totlam a posteriori %%%%%%%%%%%%%%%%



% plot(f_gr,max(totlam))

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

plot(f_gr,bru)

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



% for n = 1:length(f_gr)

%     hold on

scsearch_J1544_test.m

deleted100644 → 0
+0 −486

File deleted.

Preview size limit exceeded, changes collapsed.