Commit 91e2f76d authored by Michal Was's avatar Michal Was
Browse files

Merge branch 'master' of ligo-vcs.phys.uwm.edu:/usr/local/git/lalsuite

Original: 4a851e1f6d309d1790ac414cf431dd21ccc2d2aa
parents 900d929a e36a6b99
......@@ -96,7 +96,7 @@ try:
# perform sky localization
log.info("starting sky localization")
sky_map, epoch, elapsed_time, instruments = gracedb_sky_map(coinc_file, psd_file, "TaylorF2threePointFivePN", 10, prior_distance_power=-1, reference_frequency=120)
sky_map, epoch, elapsed_time, instruments = gracedb_sky_map(coinc_file, psd_file, "TaylorF2threePointFivePN", 10, reference_frequency=120)
log.info("sky localization complete")
# upload FITS file
......
......@@ -48,8 +48,10 @@ typedef struct{
BOOLEAN inclAutoCorr; /**< include auto-correlation terms (an SFT with itself) */
REAL8 fStart; /**< start frequency in Hz */
REAL8 fBand; /**< frequency band to search over in Hz */
REAL8 fdotStart; /**< starting value for first spindown in Hz/s*/
REAL8 fdotBand; /**< range of first spindown to search over in Hz/s */
/* REAL8 fdotStart; /\**< starting value for first spindown in Hz/s*\/ */
/* REAL8 fdotBand; /\**< range of first spindown to search over in Hz/s *\/ */
REAL8 alphaRad; /**< right ascension in radians */
REAL8 deltaRad; /**< declination in radians */
REAL8 refTime; /**< reference time for pulsar phase definition */
CHAR *sftLocation; /**< location of SFT data */
CHAR *ephemYear; /**< range of years for ephemeris file */
......@@ -85,11 +87,15 @@ int main(int argc, char *argv[]){
/* sft related variables */
MultiSFTVector *inputSFTs = NULL;
MultiPSDVector *psd = NULL;
SFTIndexList *sftIndices;
SFTPairIndexList *sftPairs;
MultiLIGOTimeGPSVector *gpsTimes = NULL;
SFTIndexList *sftIndices = NULL;
SFTPairIndexList *sftPairs = NULL;
REAL8 fMin, fMax; /* min and max frequencies read from SFTs */
REAL8 deltaF; /* Tsft; frequency resolution and time baseline of SFTs */
REAL8 deltaF, Tsft; /* frequency resolution and time baseline of SFTs */
REAL8 roughFreq; /* Approximate frequency to use for metric calculations */
REAL8Vector *sigmaUnshifted = NULL;
/* initialize and register user variables */
if ( XLALInitUserVars( &uvar ) != XLAL_SUCCESS ) {
......@@ -113,7 +119,7 @@ int main(int argc, char *argv[]){
}
deltaF = config.catalog->data[0].header.deltaF;
/*Tsft = 1. / deltaF;*/
Tsft = 1. / deltaF;
/* now read the data */
/* FIXME: need to correct fMin and fMax for Doppler shift, rngmedian bins and spindown range */
......@@ -128,6 +134,12 @@ int main(int argc, char *argv[]){
return 1;
}
/* read the timestamps from the SFTs*/
if ((gpsTimes = XLALExtractMultiTimestampsFromSFTs ( inputSFTs )) == NULL){
LogPrintf ( LOG_CRITICAL, "%s: XLALExtractMultiTimestampsFromSFTs() failed with errno=%d\n", __func__, xlalErrno );
return 1;
}
/* calculate the psd and normalize the SFTs */
if (( psd = XLALNormalizeMultiSFTVect ( inputSFTs, uvar.rngMedBlock )) == NULL){
LogPrintf ( LOG_CRITICAL, "%s: XLALNormalizeMultiSFTVect() failed with errno=%d\n", __func__, xlalErrno );
......@@ -167,6 +179,24 @@ int main(int argc, char *argv[]){
XLAL_ERROR( XLAL_EFUNC );
}
roughFreq = uvar.fStart + 0.5 * uvar.fBand;
/* Get weighting factors for calculation of metric */
if ( ( XLALCalculateCrossCorrSigmaUnshifted( &sigmaUnshifted, sftPairs, sftIndices, psd, roughFreq, Tsft) != XLAL_SUCCESS ) ) {
/* XLALDestroySFTPairIndexList(sftPairs) */
/* XLALDestroySFTIndexList(sftIndices) */
XLALDestroyMultiSFTVector ( inputSFTs );
XLALDestroyMultiPSDVector ( psd );
XLALDestroySFTCatalog (config.catalog );
XLALFree( config.edat->ephemE );
XLALFree( config.edat->ephemS );
XLALFree( config.edat );
/* de-allocate memory for user input variables */
XLALDestroyUserVars();
XLAL_ERROR( XLAL_EFUNC );
}
/* /\* get SFT parameters so that we can initialise search frequency resolutions *\/ */
/* /\* calculate deltaF_SFT *\/ */
/* deltaF_SFT = catalog->data[0].header.deltaF; /\* frequency resolution *\/ */
......@@ -260,8 +290,10 @@ int XLALInitUserVars (UserInput_t *uvar)
uvar->inclAutoCorr = FALSE;
uvar->fStart = 100.0;
uvar->fBand = 0.1;
uvar->fdotStart = 0.0;
uvar->fdotBand = 0.0;
/* uvar->fdotStart = 0.0; */
/* uvar->fdotBand = 0.0; */
uvar->alphaRad = 0.0;
uvar->deltaRad = 0.0;
uvar->rngMedBlock = 50;
/* default for reftime is in the middle */
......@@ -281,8 +313,10 @@ int XLALInitUserVars (UserInput_t *uvar)
XLALregBOOLUserStruct ( inclAutoCorr, 0, UVAR_OPTIONAL, "Include auto-correlation terms (an SFT with itself)");
XLALregREALUserStruct ( fStart, 0, UVAR_OPTIONAL, "Start frequency in Hz");
XLALregREALUserStruct ( fBand, 0, UVAR_OPTIONAL, "Frequency band to search over in Hz ");
XLALregREALUserStruct ( fdotStart, 0, UVAR_OPTIONAL, "Start value of spindown in Hz/s");
XLALregREALUserStruct ( fdotBand, 0, UVAR_OPTIONAL, "Band for spindown values in Hz/s");
/* XLALregREALUserStruct ( fdotStart, 0, UVAR_OPTIONAL, "Start value of spindown in Hz/s"); */
/* XLALregREALUserStruct ( fdotBand, 0, UVAR_OPTIONAL, "Band for spindown values in Hz/s"); */
XLALregREALUserStruct ( alphaRad, 0, UVAR_OPTIONAL, "Right ascension for directed search (radians)");
XLALregREALUserStruct ( deltaRad, 0, UVAR_OPTIONAL, "Declination for directed search (radians)");
XLALregSTRINGUserStruct( ephemYear, 0, UVAR_OPTIONAL, "String Ephemeris year range");
XLALregSTRINGUserStruct( sftLocation, 0, UVAR_REQUIRED, "Filename pattern for locating SFT data");
XLALregINTUserStruct ( rngMedBlock, 0, UVAR_OPTIONAL, "Running median block size for PSD estimation");
......
......@@ -126,12 +126,24 @@ def from_files(filelist, etg, columns=None, start=None, end=None,
end is not None and end or segments.PosInfinity)
filelist = filelist.sieve(segment=span).pfnlist()
if len(filelist) == 0:
return utils.new_ligolw_table(etg, columns=columns)
if (kwargs.has_key('channel') and
(not isinstance(kwargs['channel'], basestring) and
kwargs['channel'] is not None)):
channels = kwargs['channel']
return dict((c, utils.new_ligolw_table(etg, columns=columns)) for
c in channels)
else:
return utils.new_ligolw_table(etg, columns=columns)
if verbose:
print_verbose(0)
out = from_file(filelist[0], etg, columns=columns, start=start, end=end,
**kwargs)
extend = out.extend
if isinstance(out, dict):
def extend(in_):
for key,tab in in_.iteritems():
out[key].extend(tab)
else:
extend = out.extend
if verbose:
print_verbose(1)
for i,fp in enumerate(filelist[1:]):
......
......@@ -82,6 +82,10 @@ def ascii_trigger(line, columns=KLEINEWELLE_COLUMNS):
channel = channel.rsplit("_", 2)[0]
if 'channel' in columns:
t.channel = channel
if 'ifo' in columns and re.match("[A-Z]\d:", channel):
t.ifo = channel[:2]
elif 'ifo' in columns:
t.ifo = None
if len(dat) == 8:
(start, stop, peak, freq, energy,
amplitude, n_pix, sig) = list(map(float, dat))
......@@ -195,14 +199,19 @@ def from_ascii(filename, columns=None, start=None, end=None,
@returns a `LIGO_LW` table containing the triggers
"""
if channel and re.match("[A-Z]\d:", channel):
ifo = channel[:2]
if isinstance(channel, basestring):
channels = [channel]
elif channel is None:
channels = [None]
else:
ifo = None
channels = channel
if not columns:
columns = KLEINEWELLE_COLUMNS
if columns:
if channel:
columns.append('channel')
columns.append('ifo')
columns = set(columns)
if 'peak_time' not in columns and (start or end):
columns.add("peak_time")
......@@ -219,27 +228,31 @@ def from_ascii(filename, columns=None, start=None, end=None,
else:
check_time = False
# generate table
out = lsctables.New(lsctables.SnglBurstTable, columns=columns)
columns = out.columnnames
out = dict()
for ch in channels:
# generate table
out[ch] = lsctables.New(lsctables.SnglBurstTable, columns=columns)
columns = out[ch].columnnames
# read file and generate triggers
append = out.append
next_id = out.get_next_id
append = dict((ch, out[ch].append) for ch in channels)
next_id = out[ch].get_next_id
with open(filename, 'r') as f:
for line in f:
if _comment.match(line):
continue
t = ascii_trigger(line, columns=columns)
if channel and t.channel not in channels:
continue
t.event_id = next_id()
if channel:
t.ifo = ifo
t.channel = channel
if not check_time or (float(t.get_peak()) in span):
append(t)
append[t.channel](t)
return out
if isinstance(channel, basestring) or channel is None:
return out[channel]
else:
return out
# TODO: remove following lines when a good trigfind method is in place
......
......@@ -240,7 +240,7 @@ class TaylorF2RedSpinTemplate(Template):
0, df, self.m1 * LAL_MSUN_SI, self.m2 * LAL_MSUN_SI, self.chi,
self.bank.flow, 0, 1000000 * LAL_PC_SI, 7, 3)
# have to resize wf to next pow 2 for FFT plan caching
lal.ResizeCOMPLEX16FrequencySeries( wf, 0, ceil_pow_2(wf.data.length) )
wf = lal.ResizeCOMPLEX16FrequencySeries( wf, 0, ceil_pow_2(wf.data.length) )
return wf
def metric_match(self, other, df, **kwargs):
......@@ -399,7 +399,7 @@ class SEOBNRv1Template(Template):
self.bank.flow, 1e6*LAL_PC_SI, 0., self.spin1z, self.spin2z)
# zero-pad up to 1/df
N = int(sample_rate / df)
lal.ResizeREAL8TimeSeries(hplus, 0, N)
hplus = lal.ResizeREAL8TimeSeries(hplus, 0, N)
# taper
lalsim.SimInspiralREAL8WaveTaper(hplus.data, lalsim.LAL_SIM_INSPIRAL_TAPER_START)
......@@ -485,7 +485,7 @@ class EOBNRv2Template(Template):
self.bank.flow, 1e6*LAL_PC_SI, 0.)
# zero-pad up to 1/df
N = int(sample_rate / df)
lal.ResizeREAL8TimeSeries(hplus, 0, N)
hplus = lal.ResizeREAL8TimeSeries(hplus, 0, N)
# taper
lalsim.SimInspiralREAL8WaveTaper(hplus.data, lalsim.LAL_SIM_INSPIRAL_TAPER_START)
......
......@@ -194,3 +194,50 @@ int XLALCreateSFTPairIndexList
return XLAL_SUCCESS;
}
/** Construct vector of sigma_alpha values for each SFT pair */
/* This version uses a single frequency rather than the doppler-shifted ones */
/* Allocates memory as well */
int XLALCalculateCrossCorrSigmaUnshifted
(
REAL8Vector **sigma_alpha, /* Output: vector of sigma_alpha values */
SFTPairIndexList *pairIndexList, /* Input: list of SFT pairs */
SFTIndexList *indexList, /* Input: list of SFTs */
MultiPSDVector *psds, /* Input: PSD estimate (Sn*Tsft/2) for each SFT */
REAL8 freq, /* Frequency to extract from PSD */
REAL8 Tsft /**< SFT duration */
)
{
UINT8 j, numPairs, numSFTs, freqInd;
REAL8Vector *psdData;
REAL8FrequencySeries *psd;
SFTIndex sftIndex;
REAL8Vector *ret = NULL;
REAL8 Tsft4;
Tsft4 = SQUARE(SQUARE(Tsft));
numPairs = pairIndexList->length;
numSFTs = indexList->length;
XLAL_CHECK ( ( psdData = XLALCreateREAL8Vector ( numSFTs ) ) != NULL, XLAL_EFUNC, "XLALCreateREAL8Vector ( %d ) failed.", numSFTs );
for (j=0; j < numSFTs; j++) {
sftIndex = indexList->data[j];
psd = &(psds->data[sftIndex.detInd]->data[sftIndex.sftInd]);
freqInd = (UINT8) floor( (freq - psd->f0)/psd->deltaF + 0.5 );
psdData->data[j] = psd->data->data[freqInd];
}
XLAL_CHECK ( ( ret = XLALCreateREAL8Vector ( numPairs ) ) != NULL, XLAL_EFUNC, "XLALCreateREAL8Vector ( %d ) failed.", numPairs );
for (j=0; j < numPairs; j++) {
ret->data[j] = ( psdData->data[pairIndexList->data[j].sftNum[0]]
* psdData->data[pairIndexList->data[j].sftNum[1]]
) / Tsft4 ;
}
(*sigma_alpha) = ret;
XLALDestroyREAL8Vector ( psdData );
return XLAL_SUCCESS;
}
......@@ -120,6 +120,7 @@ int XLALCreateSFTIndexListFromMultiSFTVect
MultiSFTVector *sfts
)
;
int XLALCreateSFTPairIndexList
(
SFTPairIndexList **pairIndexList,
......@@ -130,6 +131,17 @@ int XLALCreateSFTPairIndexList
)
;
int XLALCalculateCrossCorrSigmaUnshifted
(
REAL8Vector **sigma_alpha,
SFTPairIndexList *pairIndexList,
SFTIndexList *indexList,
MultiPSDVector *psds,
REAL8 freq,
REAL8 Tsft
)
;
/*@}*/
#ifdef __cplusplus
......
......@@ -2186,7 +2186,7 @@ int XLALSimInspiralChooseFDWaveform(
ABORT_NONZERO_TIDES(waveFlags);
/* Call the waveform driver routine */
ret = XLALSimInspiralTaylorF2ReducedSpin(hptilde, phiRef, deltaF,
m1, m2, XLALSimIMRPhenomBComputeChi(m1, m2, S1z, S2z),
m1, m2, XLALSimInspiralTaylorF2ReducedSpinComputeChi(m1, m2, S1z, S2z),
f_min, f_max, r, phaseO, amplitudeO);
if (ret == XLAL_FAILURE) XLAL_ERROR(XLAL_EFUNC);
/* The above returns h(f) for optimal orientation (i=0, Fp=1, Fc=0)
......
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