Commit 7648f6fc authored by Philip Roy Charlton's avatar Philip Roy Charlton

Merge branch 'O3Development' of git.ligo.org:stochastic/stochastic into O3Development

parents 062eb3b7 339c1e5c
......@@ -677,7 +677,7 @@ for I = 1:params.numIntervalsTotal
% This is useful for the stationarity veto which excludes
% segments for which the naive sigma differs too much from
% the one calculated with the sliding PSD average.
[naiQ, result.naiVar, naiSensInt] ...
[naiQ, result.naiVar, result.naiSensInt] ...
= calOptimalFilter(params.segmentDuration, params.gamma, ...
params.fRef, params.alphaExp, ...
calPSD1, calPSD2, ...
......
......@@ -85,6 +85,7 @@ function writeOutputToMatFile(params, I, K, segmentStartTime, result)
if (params.doAllSkyComparison & params.writeSensIntsToFiles)
% write sensitivity integrands to output files
params.matFile.sensIntAllSky(I, K) = result.sensIntAllSky;
params.matFile.naiSensInt(I, K) = result.naiSensInt;
end
% Q
......
function makeDeltaSigmaCutsInPostProc(jobfile, stochOutputPrefix, dsc, doOverlap, segmentDuration, deltaF, outname, alpha)
%
% This script goes through the *_naivesigmas* output files for all of the
% jobs processed by stochastic.m, calculates a user-specified delta sigma
% cut (taking into account bias factor), and writes to file
% sGc 02/15/16
%
% Input
% - jobfile = jobfile used for particular run of stochastic.m
% - stochOutputPrefix = prefix for output of stochastic.m
% - dsc = desired delta sigma cut
% - doOverlap = 0 for no stochastic.m overlap of time segments (non-standard),
% 1 for overlap
% - segmentDuration = segment duration used for stochastic.m run
% - deltaF = frequency bin size used in stochastic.m
% - outname = name of output bad GPS times based on delta sigma cut
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%USER-SPECIFIED INPUT%%%
%Input
%jobfile='/home/sgwynne.crowder/sgwb/O1_zerolag/C02/input/jobfiles/JOB-FILE-1126623617-1136649617-master_C02-events-removed.dat';
%stochOutputPrefix='/home/sgwynne.crowder/sgwb/O1_zerolag/C02/output/a3/a3';
%Desired delta sigma cut
%dsc=0.2;
%Info for calculating bias factor for naive and sliding sigmas (info about bias factor on ilog & Stefan Ballmer's thesis)
%doOverlap=1;
%segmentDuration=192;
%deltaF=1/32;
if ~exist('alpha','var')
alpha=0;
end
fref=1;
if doOverlap
segs=segmentDuration*deltaF*2-1;
else
segs=segmentDuration*deltaF;
end
nn=2*9/11*segs; %9/11 for Welch factor
bf_ss=nn/(nn-1); %sliding sigma bias factor
nn=9/11*segs;
bf_ns=nn/(nn-1); %naive sigma bias factor
%bf_ss=1; %trying getting rid of bias factor
%bf_ns=1; %trying getting rid of bias factor
%Output
%outname='/home/sgwynne.crowder/sgwb/O1_zerolag/C02/output/a0_postproc/badGPSTimes_a3.txt';
%%%%%%%%%%%%%%%%%%%%%%%%%%
%Parameters
jobs=load(jobfile);
num_jobs=size(jobs,1);
%Process
fid=fopen(outname,'w+');
for ii=1:num_jobs %step through jobs
%sigmasfilename=[stochOutputPrefix '_naivesigmas.job' num2str(ii) '.trial1.dat'];
disp(['Processing job ' num2str(ii) ' of ' num2str(num_jobs)]);
sigmasfilename=[stochOutputPrefix '.job' num2str(ii) '.mat']
if exist(sigmasfilename)
tmp=load(sigmasfilename); %column 1 = start sec, column 2 = naive sigma, column 3 = sliding sigma
if ~isfield(tmp,'segmentStartTime')
disp('No data for this job')
continue
end
naiSensInt=tmp.naiSensInt;
sensInt=tmp.sensInt;
segmentStartTime=tmp.segmentStartTime;
% define frequency array
flow=naiSensInt(1).flow;
deltaF=naiSensInt(1).deltaF;
nbins=length(naiSensInt(1).data);
fmax=flow+(nbins-1)*deltaF;
freqs=transpose(flow:deltaF:fmax);
% Define H(f)=((f/fref)^alpha
Hf=(freqs/fref).^alpha;
% Rescale the naive snsiviity integrand and integrate to get naive sigma
nsegs=length(naiSensInt)
naiSigmaAlpha=zeros(nsegs,1);
ccSigmaAlpha=zeros(nsegs,1);
for jj=1:nsegs
naiSensIntAlpha=naiSensInt(jj).data.*Hf.^2;
naiSigmaAlpha(jj)=1/sqrt(sum(naiSensIntAlpha*deltaF));
sensIntAlpha=sensInt(jj).data.*Hf.^2;
ccSigmaAlpha(jj)=1/sqrt(sum(sensIntAlpha*deltaF));
end
sigmas=[segmentStartTime,naiSigmaAlpha,ccSigmaAlpha];
[m,n]=size(sigmas);
if m~=0
for jj=1:m %write times exceeding delta sigma cut to file
if abs(sigmas(jj,3)*bf_ss-sigmas(jj,2)*bf_ns)/(sigmas(jj,3)*bf_ss)>dsc %apply bias factor
%fprintf(fid,'%i %e %e\n',sigmas(jj,1),sigmas(jj,2),sigmas(jj,3));
fprintf(fid,'%i\n',sigmas(jj,1));
end
end
end
end
clear naiSigmaAlpha;
clear ccSigmaAlpha;
clear freqs;
clear Hf;
clear tmp;
end %jobs
fclose(fid)
return;
% Parameters which are the same for all runs of makeDeltaSigmaCutsInPostProc
jobfile='../input/jobs/JOB-FILE-1238166018-1239458471.dat';
dsc=0.2;
doOverlap=1;
segmentDuration=192;
deltaF=0.03125;
% For a=-5
stochOutputPrefix_aminus5='../output/a0/a0/H1L1';
outname_aminus5='./tmp_a-5.txt';
makeDeltaSigmaCutsInPostProc_naiSensInt(jobfile,stochOutputPrefix_aminus5,dsc,doOverlap,segmentDuration,deltaF,outname_aminus5,-5);
% For a=0
stochOutputPrefix_a0='../output/a0/a0/H1L1';
outname_a0='./tmp_a0.txt';
makeDeltaSigmaCutsInPostProc_naiSensInt(jobfile,stochOutputPrefix_a0,dsc,doOverlap,segmentDuration,deltaF,outname_a0,0);
% For a=3
stochOutputPrefix_a3='../output/a0/a0/H1L1';
outname_a3='./tmp_a3.txt';
makeDeltaSigmaCutsInPostProc_naiSensInt(jobfile,stochOutputPrefix_a3,dsc,doOverlap,segmentDuration,deltaF,outname_a3,3);
% Super cut
outname_supercut='./badGPSTimes_naiSensInt.dat'
produceSuperCut(outname_aminus5,outname_a0,outname_a3,outname_supercut);
% cleanup
delete(outname_aminus5,outname_a0,outname_a3);
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment