diff --git a/zdm/MC_sample/loading.py b/zdm/MC_sample/loading.py index 99dfd0f5..84221e9d 100644 --- a/zdm/MC_sample/loading.py +++ b/zdm/MC_sample/loading.py @@ -91,11 +91,15 @@ def set_state(alpha_method=1, cosmo=Planck18): def survey_and_grid(survey_name:str='CRAFT/CRACO_1_5000', + opdir='', + bPosNum=0, init_state=None, state_dict=None, iFRB:int=0, alpha_method=1, NFRB:int=100, lum_func:int=2,sdir=None,nz=500,ndm=1400, - nbins=5,edir=''): + nbins=5, cluster=False, + clusterRedshift=np.nan, + lensing=False): """ Load up a survey and grid for a CRACO mock dataset Args: @@ -108,10 +112,7 @@ def survey_and_grid(survey_name:str='CRAFT/CRACO_1_5000', 0=power-law, 1=gamma, 2=gamma+spline. Defaults to 0. state_dict (dict, optional): Used to init state instead of alpha_method, lum_func parameters - sdir (string, optional): Directory containing surveys - edir (string, optional): - Directory containing efficiency files if using FRB-specific responses - + Raises: IOError: [description] @@ -139,18 +140,31 @@ def survey_and_grid(survey_name:str='CRAFT/CRACO_1_5000', datdir=resource_filename('zdm', 'GridData'), zlog=False,nz=nz,ndm=ndm) + ############## Initialise surveys ############## if sdir is not None: print("Searching for survey in directory ",sdir) else: sdir = os.path.join(resource_filename('zdm', 'craco'), 'MC_Surveys') - isurvey = survey.load_survey(survey_name, state, dmvals, - NFRB=NFRB, sdir=sdir, nbins=nbins, - iFRB=iFRB, edir=edir) + isurvey = survey.load_survey( + survey_name = survey_name, + state = state, + opdir = opdir, + bPosNum = bPosNum, + dmvals = dmvals, + cluster = cluster, + clusterRedshift = clusterRedshift, + zvals = zvals, + lensing = lensing, + NFRB=NFRB, + sdir=sdir, + nbins=nbins, + iFRB=iFRB + ) # generates zdm grid grids = misc_functions.initialise_grids( - [isurvey], zDMgrid, zvals, dmvals, state, wdist=True) + [isurvey], opdir, bPosNum, zDMgrid, zvals, dmvals, state, wdist=True, cluster=cluster, clusterRedshift = clusterRedshift) print("Initialised grid") # Return Survey and Grid diff --git a/zdm/data/Surveys/CHIME_SynthBeam.ecsv b/zdm/data/Surveys/CHIME_SynthBeam.ecsv new file mode 100644 index 00000000..1b7a6631 --- /dev/null +++ b/zdm/data/Surveys/CHIME_SynthBeam.ecsv @@ -0,0 +1,28 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: '{"observing": {"NORM_FRB": 0, "TOBS": 1, "MAX_DM": 7000}, +# "telescope": {"BMETHOD": 0, "DIAM": 80.0, "NBEAMS": 1, +# "NBINS": 100, "THRESH":5.0, "TRES":1.0, +# "FRES":0.025,"FBAR":600, "BTHRESH":0.001, "BW":400}}'} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z +DUMMY 400 300 30.0 600 0.025 "" "" 8 8. 5. 1.0 1 1.0 "" "" 0.3 diff --git a/zdm/data/Surveys/CHORD.ecsv b/zdm/data/Surveys/CHORD.ecsv new file mode 100644 index 00000000..76f09dcb --- /dev/null +++ b/zdm/data/Surveys/CHORD.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 0.625,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD1_5.ecsv b/zdm/data/Surveys/CHORD1_5.ecsv new file mode 100644 index 00000000..b0d82927 --- /dev/null +++ b/zdm/data/Surveys/CHORD1_5.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 420.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 240,\n \"BTHRESH\": 0.001,\n \"THRESH\": 11.18,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD2_5.ecsv b/zdm/data/Surveys/CHORD2_5.ecsv new file mode 100644 index 00000000..3094c02e --- /dev/null +++ b/zdm/data/Surveys/CHORD2_5.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 660.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 240,\n \"BTHRESH\": 0.001,\n \"THRESH\": 11.18,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD3_5.ecsv b/zdm/data/Surveys/CHORD3_5.ecsv new file mode 100644 index 00000000..dfa5f5ad --- /dev/null +++ b/zdm/data/Surveys/CHORD3_5.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 240,\n \"BTHRESH\": 0.001,\n \"THRESH\": 11.18,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD4_5.ecsv b/zdm/data/Surveys/CHORD4_5.ecsv new file mode 100644 index 00000000..aed3ac0c --- /dev/null +++ b/zdm/data/Surveys/CHORD4_5.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 1140.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 240,\n \"BTHRESH\": 0.001,\n \"THRESH\": 11.18,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD5_5.ecsv b/zdm/data/Surveys/CHORD5_5.ecsv new file mode 100644 index 00000000..17576b73 --- /dev/null +++ b/zdm/data/Surveys/CHORD5_5.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 1380.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 240,\n \"BTHRESH\": 0.001,\n \"THRESH\": 11.18,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORDSB1_4.ecsv b/zdm/data/Surveys/CHORDSB1_4.ecsv new file mode 100644 index 00000000..783a813f --- /dev/null +++ b/zdm/data/Surveys/CHORDSB1_4.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 6.0,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0019,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5.0,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_0.ecsv b/zdm/data/Surveys/CHORD_BeamPos_0.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_0.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_1.ecsv b/zdm/data/Surveys/CHORD_BeamPos_1.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_1.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_10.ecsv b/zdm/data/Surveys/CHORD_BeamPos_10.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_10.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_11.ecsv b/zdm/data/Surveys/CHORD_BeamPos_11.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_11.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_12.ecsv b/zdm/data/Surveys/CHORD_BeamPos_12.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_12.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_13.ecsv b/zdm/data/Surveys/CHORD_BeamPos_13.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_13.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_14.ecsv b/zdm/data/Surveys/CHORD_BeamPos_14.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_14.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_15.ecsv b/zdm/data/Surveys/CHORD_BeamPos_15.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_15.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_16.ecsv b/zdm/data/Surveys/CHORD_BeamPos_16.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_16.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_17.ecsv b/zdm/data/Surveys/CHORD_BeamPos_17.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_17.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_18.ecsv b/zdm/data/Surveys/CHORD_BeamPos_18.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_18.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_19.ecsv b/zdm/data/Surveys/CHORD_BeamPos_19.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_19.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_2.ecsv b/zdm/data/Surveys/CHORD_BeamPos_2.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_2.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_3.ecsv b/zdm/data/Surveys/CHORD_BeamPos_3.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_3.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_4.ecsv b/zdm/data/Surveys/CHORD_BeamPos_4.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_4.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_5.ecsv b/zdm/data/Surveys/CHORD_BeamPos_5.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_5.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_6.ecsv b/zdm/data/Surveys/CHORD_BeamPos_6.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_6.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_7.ecsv b/zdm/data/Surveys/CHORD_BeamPos_7.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_7.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_8.ecsv b/zdm/data/Surveys/CHORD_BeamPos_8.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_8.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORD_BeamPos_9.ecsv b/zdm/data/Surveys/CHORD_BeamPos_9.ecsv new file mode 100644 index 00000000..8624f79c --- /dev/null +++ b/zdm/data/Surveys/CHORD_BeamPos_9.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 5,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CHORDdiv.ecsv b/zdm/data/Surveys/CHORDdiv.ecsv new file mode 100644 index 00000000..321e1568 --- /dev/null +++ b/zdm/data/Surveys/CHORDdiv.ecsv @@ -0,0 +1,27 @@ +# %ECSV 1.0 +# --- +# datatype: +# - {name: TNS, datatype: string} +# - {name: BW, datatype: float64} +# - {name: DM, datatype: float64} +# - {name: DMG, datatype: float64} +# - {name: FBAR, datatype: float64} +# - {name: FRES, datatype: float64} +# - {name: Gb, datatype: string, subtype: 'float64[null]'} +# - {name: Gl, datatype: string, subtype: 'float64[null]'} +# - {name: SNR, datatype: float64} +# - {name: SNRTHRESH, datatype: float64} +# - {name: THRESH, datatype: float64} +# - {name: TRES, datatype: float64} +# - {name: WIDTH, datatype: float64} +# - {name: XC, datatype: string} +# - {name: XDec, datatype: string} +# - {name: XRA, datatype: string} +# - {name: Z, datatype: float64} +# meta: !!omap +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 0,\n \"TOBS\": 1\n },\n \ +# \"telescope\": {\n \"BEAM\": \"IGNORE\",\n \"DIAM\": 48,\n \"BMETHOD\": 0,\n \ +# \"FBAR\": 900.0,\n \"NBEAMS\": 1,\n \"TRES\": 1.0,\n \"FRES\": 0.0375,\n \ +# \"BW\": 1200,\n \"BTHRESH\": 0.001,\n \"THRESH\": 11.18,\n \"NBINS\": 10\n }\n}"} +# schema: astropy-2.0 +TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XC XDec XRA Z diff --git a/zdm/data/Surveys/CRAFT_class_I_and_II.ecsv b/zdm/data/Surveys/CRAFT_class_I_and_II.ecsv index 7a43ee21..a55e7da0 100644 --- a/zdm/data/Surveys/CRAFT_class_I_and_II.ecsv +++ b/zdm/data/Surveys/CRAFT_class_I_and_II.ecsv @@ -25,7 +25,7 @@ # - {name: XRA, datatype: string} # - {name: Z, datatype: float64} # meta: !!omap -# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 20,\n \"TOBS\": 1274.6,\n \"MAX_IDT\": 4096\n},\n \"telescope\": {\n \ +# - {survey_data: "{\n \"observing\": {\n \"NORM_FRB\": 20,\n \"TOBS\": 1274.6\n },\n \"telescope\": {\n \ # \ \"BEAM\": \"lat50_log\",\n \"DIAM\": 12.0,\n \"NBEAMS\": 36,\n \"NBINS\": 5\n }\n}"} # schema: astropy-2.0 TNS BW DM DMG FBAR FRES Gb Gl SNR SNRTHRESH THRESH TRES WIDTH XA XD XDec XE XF XG XH XI XRA Z diff --git a/zdm/energetics.py b/zdm/energetics.py index f47b0b8a..100d7033 100644 --- a/zdm/energetics.py +++ b/zdm/energetics.py @@ -39,6 +39,8 @@ import numpy as np from scipy import interpolate import mpmath +from astropy.cosmology import Planck18 as cosmo +from pathlib import Path from IPython import embed @@ -199,7 +201,7 @@ def vector_cum_power_law(Eth, *params): ndarray Fraction of bursts with E > Eth. Returns 1 for Eth < Emin, 0 for Eth > Emax. """ - params=np.array(params) + #params=np.array(params) Emin=params[0] Emax=params[1] gamma=params[2] @@ -255,7 +257,7 @@ def vector_diff_power_law(Eth,*params): low=np.where(Eth < Emin)[0] if len(low) > 0: - result[low]=0. + result[low]=0. high=np.where(Eth > Emax)[0] if len(high) > 0: result[high]=0. @@ -328,8 +330,6 @@ def vector_cum_gamma_spline(Eth: np.ndarray, *params): ----- Automatically initializes splines for new gamma values if needed. """ - global SplineLog - params=np.array(params) Emin=params[0] Emax=params[1] @@ -340,10 +340,7 @@ def vector_cum_gamma_spline(Eth: np.ndarray, *params): Eth_Emax = Eth/Emax if gamma not in igamma_splines.keys(): init_igamma_splines([gamma]) - if SplineLog: - numer = 10**interpolate.splev(np.log10(Eth_Emax), igamma_splines[gamma]) - else: - numer = interpolate.splev(Eth_Emax, igamma_splines[gamma]) + numer = interpolate.splev(Eth_Emax, igamma_splines[gamma]) result=numer/norm # Low end @@ -463,6 +460,82 @@ def vector_diff_gamma(Eth, *params): result= (Eth/Emax)**(gamma-1) * np.exp(-Eth/Emax) / norm low= Eth < Emin - result[low]=0. - + result[low]=0. + return result + +########### lensing modified ####### + +def distanceFraction(zD, zS): + Dds = cosmo.angular_diameter_distance_z1z2(zD,zS) + Ds = cosmo.angular_diameter_distance(zS) + return Dds/Ds + +def lensingPDF(mu, zD, zS, bPosNum, beami, opdir): + formatted_redshift = "{:03.2f}".format(zS) + parent = Path(opdir).parent + x = np.load(str(parent)+'/mux.npy') + yFull = np.load(str(parent)+'/pmux_BP_'+str(bPosNum)+str(formatted_redshift)+'.npy') + if np.sum(np.isnan(yFull[:,beami]))==len(yFull[:,beami]): + return np.ones(len(mu))*np.nan + y = yFull[:,beami] + interpFunc = interpolate.interp1d(x,y, bounds_error=False, fill_value=0) + return interpFunc(np.log10(mu)) + +def vector_cum_lensed_power_law(Eth,*params): + """ Calculates the fraction of bursts above a certain power law + for a given Eth. + """ + #params=np.array(params) + Emin=params[0] + Emax=params[1] + gamma=params[2] + zvals=params[4] + beami = params[5] + surveyName = params[6] + zD = params[7] + opdir = params[8] + bPosNum = params[9] + #print(Eth, Emin, Emax, gamma) + logEn = np.log(Emin) + logEx = np.log(Emax) + logSpacing = 0.01 + logERange = np.arange(logEn, logEx+logSpacing, logSpacing) + muNum = int((10+2)/logSpacing)+1 + #print(muNum) + logMu = np.arange(-2, muNum*logSpacing, logSpacing) + #print(len(logERange), len(logMu)) + result = np.zeros(Eth.shape) + for i in range(len(zvals)): + if zD < np.round(zvals[i],2): + probGrid = lensingPDF(np.e**logMu, zD, zvals[i], bPosNum, beami, opdir) + if np.sum(np.isnan(probGrid))==len(probGrid): + result[i,:]=vector_cum_power_law(Eth[i,:],*params) + else: + phiGrid = vector_diff_power_law(np.e**logERange, *params) + phiL = np.convolve(probGrid*(np.e**logMu), phiGrid)*logSpacing + logE_muRange = np.arange(logEn-2, np.amax(logERange)+np.amax(logMu),logSpacing) + #print(np.amin(np.e**logE_muRange), np.amax(np.e**logE_muRange)) + phiLCumConv = np.cumsum(np.flip((np.e**logE_muRange)*phiL*logSpacing)) + interpFunc = interpolate.interp1d(np.flip(np.e**logE_muRange), phiLCumConv, bounds_error=False, fill_value=(1.0,0.0)) + #iprint(interpFunc(1e31),interpFunc(1e32)) + result[i,:] = interpFunc(Eth[i,:]) + else: + result[i,:]=vector_cum_power_law(Eth[i,:],*params) + + return result + + + +def array_cum_lensed_power_law(Eth,*params): + """ Calculates the fraction of bursts above a certain power law + for a given Eth, where Eth is an N-dimensional array + """ + dims=Eth.shape + #if gamma >= 0: #handles crazy dodgy cases. Or just return 0? + # result=np.zeros([Eth.size]) + # result[np.where(Eth < Emax)]=1. + # result=result.reshape(dims) + # Eth=Eth.reshape(dims) + # return result + result=vector_cum_lensed_power_law(Eth,*params) return result diff --git a/zdm/grid.py b/zdm/grid.py index 5e940546..f4dc82c7 100644 --- a/zdm/grid.py +++ b/zdm/grid.py @@ -14,8 +14,9 @@ ------------ - Builds normalized 2D probability grids for expected FRB rates - Handles beam response and detection efficiency -- Supports multiple luminosity functions (power-law, gamma) +- Supports multiple luminosity functions (power-law, gamma, lensed power-law) - Efficient updating for MCMC parameter exploration +- Cluster/lensing support: per-beam smearing and lensed luminosity functions Example ------- @@ -65,9 +66,19 @@ class Grid: Parameter state used for grid calculation. survey : survey.Survey Associated survey object. + cluster : bool + If True, use per-beam smearing and lensed luminosity function (LF4). + clusterRedshift : float + Redshift of the cluster/lens (used when cluster=True and LF=4). + bPosNum : int + Beam position number passed to the lensed luminosity function. + opdir : str + Output directory passed to the lensed luminosity function. """ - def __init__(self, survey, state, zDMgrid, zvals, dmvals, smear_mask, wdist=None, prev_grid=None): + def __init__(self, survey, state, zDMgrid, zvals, dmvals, smear_mask, + wdist=None, prev_grid=None, + cluster=False, clusterRedshift=np.nan, bPosNum=0, opdir=''): """Initialize the Grid for a survey and parameter state. Parameters @@ -84,12 +95,22 @@ def __init__(self, survey, state, zDMgrid, zvals, dmvals, smear_mask, wdist=None dmvals : ndarray DM bin centers in pc/cm^3. Bins span [DM - dDM/2, DM + dDM/2]. smear_mask : ndarray - 1D convolution kernel for host DM smearing. + Convolution kernel for host DM smearing. For cluster mode this + must be 3D (nz, ndm, n_beam). wdist : bool, optional If True, include width distribution effects. prev_grid : Grid, optional Another Grid with same z/DM values but different survey. Allows reusing pre-computed cosmological quantities. + cluster : bool, optional + If True, activate per-beam smearing and lensed luminosity function + support (LF=4). Default False. + clusterRedshift : float, optional + Redshift of the cluster/lens. Used when cluster=True and LF=4. + bPosNum : int, optional + Beam position number passed to the lensed LF. Default 0. + opdir : str, optional + Output directory passed to the lensed LF. Default ''. """ self.grid = None self.survey = survey @@ -102,6 +123,12 @@ def __init__(self, survey, state, zDMgrid, zvals, dmvals, smear_mask, wdist=None # State self.state = state self.MCinit = False + # Cluster / lensing parameters + self.cluster = cluster + self.clusterRedshift = clusterRedshift + self.bPosNum = bPosNum + self.opdir = opdir + self.source_function = cos.choose_source_evolution_function( state.FRBdemo.source_evolution ) @@ -133,7 +160,7 @@ def __init__(self, survey, state, zDMgrid, zvals, dmvals, smear_mask, wdist=None weights = survey.wplist # Warning -- THRESH could be different for each FRB, but we don't treat it that way self.calc_thresholds(survey.meta["THRESH"], - efficiencies,weights=weights) + efficiencies, weights=weights) else: # this is called when the grid is not iterating over widths internally efficiencies = survey.mean_efficiencies # one dimension @@ -142,11 +169,11 @@ def __init__(self, survey, state, zDMgrid, zvals, dmvals, smear_mask, wdist=None # Calculate self.calc_pdv() - self.set_evolution() # sets star-formation rate scaling with z - here, no evoltion... + self.set_evolution() # sets star-formation rate scaling with z - here, no evolution... self.calc_rates() # includes sfr smearing factors and pdv mult def init_luminosity_functions(self): - """ Set the luminsoity function for FRB energetics """ + """ Set the luminosity function for FRB energetics """ if self.luminosity_function == 0: # Power-law self.array_cum_lf = energetics.array_cum_power_law self.vector_cum_lf = energetics.vector_cum_power_law @@ -168,6 +195,9 @@ def init_luminosity_functions(self): self.vector_cum_lf = energetics.vector_cum_gamma_linear self.array_diff_lf = energetics.array_diff_gamma self.vector_diff_lf = energetics.vector_diff_gamma + elif self.luminosity_function == 4: # Lensed power-law (cluster mode) + self.array_cum_lf = energetics.array_cum_lensed_power_law + self.vector_cum_lf = energetics.vector_cum_lensed_power_law else: raise ValueError( "Luminosity function must be 0, not ", self.luminosity_function @@ -241,7 +271,7 @@ def get_z_coeffs(self,zlist): return izs1, izs2, dkzs1, dkzs2 - def check_grid(self,TOLERANCE = 1e-6): + def check_grid(self, TOLERANCE=1e-6): """ Check that the grid values are behaving as expected @@ -292,7 +322,7 @@ def check_grid(self,TOLERANCE = 1e-6): maxoff ** 0.5, "detected, aborting", ) - + # Ensures that log-spaced bins are truly bin centres if not self.zlog and np.abs(self.zvals[0] - self.dz/2.) > TOLERANCE*self.dz: raise ValueError( @@ -300,8 +330,7 @@ def check_grid(self,TOLERANCE = 1e-6): " first value ",self.zvals[0]," expected to be half of spacing ", self.dz,", aborting..." ) - - + expectation = self.ddm * np.arange(0, self.ndm) + self.dmvals[0] diff = self.dmvals - expectation maxoff = np.max(diff ** 2) @@ -357,7 +386,8 @@ def set_evolution(self): # ,n,alpha=None): self.nuObs / self.nuRef ) ** -self.state.energy.alpha # alpha positive, nuObs 0.: - self.smear_zgrid = self.smear_z(self.rates,self.zsigma) - self.rates=self.smear_zgrid - + if self.cluster: + # --- Cluster mode: compute rates per-beam then sum --- + # smear_grid is 3D [nz, ndm, n_beam]; b_fractions is 3D [nz, ndm, n_beam] + tempRates = np.zeros([self.grid.shape[0], self.grid.shape[1], len(self.beam_b)]) + for i in range(len(self.beam_b)): + sfr_smear_i = np.multiply(self.smear_grid[:, :, i].T, self.sfr).T + tempRates[:, :, i] = (self.b_fractions[:, :, i].T * self.dV).T * sfr_smear_i + self.sfr_smear = np.multiply(self.smear_grid[:, :, 0].T, self.sfr).T # kept for compatibility + self.rates = np.sum(tempRates, axis=2) + else: + # --- Standard mode --- + # zfraction describes the fraction of host galaxies estimated to be + # visible at a given redshift. + if self.survey.survey_data.observing.Z_FRACTION is not None: + fdir = str(resources.files('zdm').joinpath('data/optical')) + ffile = fdir + "/fz_"+str(self.survey.survey_data.observing.Z_FRACTION)+".npy" + zfile = fdir + "/z_"+str(self.survey.survey_data.observing.Z_FRACTION)+".npy" + self.construct_fz(ffile, zfile) + self.sfr *= self.fz + + self.sfr_smear = np.multiply(self.smear_grid.T, self.sfr).T + + # below could pass more parameters internally, but this may not + # be the final implementation + self.rates = self.pdv * self.sfr_smear + self.zsigma = self.survey.survey_data.observing.Z_PHOTO + if self.zsigma > 0.: + self.smear_zgrid = self.smear_z(self.rates, self.zsigma) + self.rates = self.smear_zgrid + def get_rates(self): """ Returns rates, multiplied by the relevant constant, @@ -587,13 +684,14 @@ def calc_thresholds(self, F0:float, Args: F0 (float): base survey threshold - eff_table ([type]): table of efficiencies corresponding to DM-values. 1, 2, or 3 dimensions! - bandwidth ([type], optional): [description]. Defaults to 1e9. - nuObs ([float], optional): survey frequency (affects sensitivity via alpha - only for alpha method) - Defaults to 1.3e9. - nuRef ([float], optional): reference frequency we are calculating thresholds at - Defaults to 1.3e9. - weights ([type], optional): [description]. Defaults to None. + eff_table: table of efficiencies corresponding to DM-values. + 1D → (NDM,) — single width, standard + 2D → (nW, NDM) — multiple widths, standard + 3D → (nW, NDM, NZ) — multiple widths, z-dependent + 4D → (nW, NDM, NZ, N_beam) — cluster/per-beam mode + bandwidth (float, optional): Defaults to 1e9. + nuRef (float, optional): reference frequency. Defaults to 1.3e9. + weights: width weights. Required for multi-width tables. Raises: ValueError: [description] @@ -608,7 +706,7 @@ def calc_thresholds(self, F0:float, self.nthresh = 1 self.eff_weights = np.array([1]) self.eff_table = np.array([eff_table]) # make it an extra dimension - else: # multiple FRB widths: dimensions nW x NDM + else: # multiple FRB widths: 2D, 3D, or 4D # check that the weights dimensions check out self.nthresh = eff_table.shape[0] # number of width bins. if weights is not None: @@ -628,59 +726,80 @@ def calc_thresholds(self, F0:float, self.nw = self.eff_weights.shape[0] - # now two or three dimensions + # now two, three, or four dimensions Eff_thresh = F0 / self.eff_table self.EF(self.state.energy.alpha, bandwidth) # sets FtoE values - could have been done *WAY* earlier - self.thresholds = np.zeros([self.nthresh, self.zvals.size, self.dmvals.size]) - - # Performs an outer multiplication of conversion from fluence to energy. - # The FtoE array has one value for each redshift. - # The effective threshold array has one value for each combination of - # FRB width (nthresh) and DM. - # We loop over nthesh and generate a NDM x Nz array for each - for i in np.arange(self.nthresh): - if self.eff_table.ndim == 2: - self.thresholds[i,:,:] = np.outer(self.FtoE, Eff_thresh[i,:]) - else: - self.thresholds[i,:,:] = ((Eff_thresh[i,:,:]).T * self.FtoE).T + if eff_table.ndim == 4: + # --- Cluster / per-beam mode --- + # thresholds shape: [nthresh, n_beam, nz, ndm] + n_beam = self.beam_b.size + self.thresholds = np.zeros([self.nthresh, n_beam, self.zvals.size, self.dmvals.size]) + for i in np.arange(self.nthresh): + for j in range(n_beam): + # Eff_thresh[i, :, :, j] has shape (NDM, NZ) or similar + self.thresholds[i, j, :, :] = (self.FtoE * (Eff_thresh[i, :, :, j])).T + else: + # --- Standard mode: thresholds shape [nthresh, nz, ndm] --- + self.thresholds = np.zeros([self.nthresh, self.zvals.size, self.dmvals.size]) + for i in np.arange(self.nthresh): + if self.eff_table.ndim == 2: + self.thresholds[i,:,:] = np.outer(self.FtoE, Eff_thresh[i,:]) + else: + self.thresholds[i,:,:] = ((Eff_thresh[i,:,:]).T * self.FtoE).T def smear_dm(self, smear:np.ndarray): # ,mean:float,sigma:float): """ Smears DM using the supplied array. Example use: DMX contribution - smear_grid is created in place + smear_grid is created in place. + + For cluster mode (self.cluster=True), smear must be 3D + (nz, ndm, n_beam) and smear_grid will be 3D (nz, ndm, n_beam). Args: smear (np.ndarray): Smearing array """ # just easier to have short variables for this - - ls = smear.size lz, ldm = self.grid.shape - if not hasattr(self, "smear_grid"): - self.smear_grid = np.zeros([lz, ldm]) - self.smear = smear - - # this method is O~7 times faster than the 'brute force' above for large arrays - for i in np.arange(lz): - # we need to get the length of mode='same', BUT - # we do not want it 'centred', hence must make cut on full - if smear.ndim == 1: - self.smear_grid[i, :] = np.convolve( - self.grid[i, :], smear, mode="full" - )[0:ldm] - elif smear.ndim == 2: - self.smear_grid[i, :] = np.convolve( - self.grid[i, :], smear[i, :], mode="full" - )[0:ldm] - else: + if self.cluster: + # --- Cluster mode: per-beam smearing --- + if smear.ndim != 3: raise ValueError( - "Wrong number of dimensions for DM smearing ", smear.shape + "Wrong number of dimensions for cluster DM smearing ", smear.shape ) + self.smear = smear + self.smear_grid = np.zeros([lz, ldm, len(self.beam_b)]) + for j in range(len(self.beam_b)): + for i in np.arange(lz): + self.smear_grid[i, :, j] = np.convolve( + self.grid[i, :], smear[i, :, j], mode="full" + )[0:ldm] + else: + # --- Standard mode --- + if not hasattr(self, "smear_grid"): + self.smear_grid = np.zeros([lz, ldm]) + self.smear = smear + + # this method is O~7 times faster than the 'brute force' above for large arrays + for i in np.arange(lz): + # we need to get the length of mode='same', BUT + # we do not want it 'centred', hence must make cut on full + if smear.ndim == 1: + self.smear_grid[i, :] = np.convolve( + self.grid[i, :], smear, mode="full" + )[0:ldm] + elif smear.ndim == 2: + self.smear_grid[i, :] = np.convolve( + self.grid[i, :], smear[i, :], mode="full" + )[0:ldm] + else: + raise ValueError( + "Wrong number of dimensions for DM smearing ", smear.shape + ) def get_p_zgdm(self, DMs: np.ndarray): """ Calcuates the probability of redshift given a DM @@ -750,7 +869,7 @@ def GenMCSample(self, N, Poisson=False): def initMC(self): """ - Initialises the MC sample, if it has not been doen already + Initialises the MC sample, if it has not been done already This uses a great deal of RAM - hence, do not do this lightly! """ @@ -838,25 +957,22 @@ def GenMCFRB(self, Emax_boost): NOTE: currently, the actual FRB widths are not part of 'grid' only the relative probabilities of any given width. Hence, this routine only returns the integer of the width bin - not the width itelf. + not the width itself. Args: - pwb (optional): probability(width,beam) - Emax_boost (float, optional): + Emax_boost (float): Allow for larger energies than Emax The factor is logarithmic, i.e. Emax_boost = 2. allows for 10**2 higher energies than Emax Returns: - tuple: FRBparams=[MCz,MCDM,MCb,j,MCs], pwb values + list: FRBparams=[MCz, MCDM, MCb, MCs, MCw] These are: MCz: redshift MCDM: dispersion measure (extragalactic) MCb: beam value - j: MCs: SNR/SNRth value of FRB MCw: width value of FRB - [MCz, MCDM, MCb, j, MCs, MCw] """ # shorthand @@ -1224,7 +1340,7 @@ def update(self, vparams: dict, ALL=False, prev_grid=None): if calc_pdv or ALL: self.calc_pdv() if set_evol or ALL: - self.set_evolution() # sets star-formation rate scaling with z - here, no evoltion... + self.set_evolution() # sets star-formation rate scaling with z - here, no evolution... if new_sfr_smear or ALL: self.calc_rates() # includes sfr smearing factors and pdv mult elif new_pdv_smear: @@ -1262,7 +1378,7 @@ def chk_upd_param(self, param: str, vparams: dict, update=False): # return updated - def smear_z(self,array,zsigma): + def smear_z(self, array, zsigma): """ Smear a 2-D z-DM grid along the redshift axis to account for photometric redshift uncertainty. @@ -1304,25 +1420,25 @@ def smear_z(self,array,zsigma): ``self.survey.survey_data.observing.Z_PHOTO`` and applied to ``self.rates`` after the FRB rate grid has been computed. """ - r,c=array.shape + r, c = array.shape # get sigma in grid units - sigma=zsigma/(self.dz) - smear_size=int(self.state.photo.sigma_width*sigma) - smear_size=smear_size-smear_size%2+1 - smear_arr=np.linspace(-(smear_size-1)/2,(smear_size-1)//2,smear_size) + sigma = zsigma / (self.dz) + smear_size = int(self.state.photo.sigma_width * sigma) + smear_size = smear_size - smear_size % 2 + 1 + smear_arr = np.linspace(-(smear_size-1)/2, (smear_size-1)//2, smear_size) # makes the approximation of taking the central value in the bin. - smear_arr=np.exp(-(smear_arr**2)/(2*(sigma**2))) - #normalise - smear_arr/=np.sum(smear_arr) + smear_arr = np.exp(-(smear_arr**2) / (2*(sigma**2))) + # normalise + smear_arr /= np.sum(smear_arr) - smear_zgrid=np.zeros([r,c]) + smear_zgrid = np.zeros([r, c]) for i in range(c): - smear_zgrid[:,i]=np.convolve(array[:,i],smear_arr,mode="same") + smear_zgrid[:, i] = np.convolve(array[:, i], smear_arr, mode="same") return smear_zgrid - def construct_fz(self,ffile,zfile): + def construct_fz(self, ffile, zfile): """ linearly interpolates passed fz values onto own zvals array @@ -1334,8 +1450,8 @@ def construct_fz(self,ffile,zfile): z = np.load(zfile) from scipy.interpolate import interp1d - f=interp1d(z,fz,kind="linear",bounds_error=False) - newfz=f(self.zvals) + f = interp1d(z, fz, kind="linear", bounds_error=False) + newfz = f(self.zvals) # check for unphysical values toolow = np.where(newfz < 0.) @@ -1344,8 +1460,3 @@ def construct_fz(self,ffile,zfile): newfz[toohigh] = 1. self.fz = newfz - - - #np.save(path+"/"+name+"_fz",newfz) - #np.save(path+"/"+name+"_z",newz) - diff --git a/zdm/loading.py b/zdm/loading.py index a51c1145..54d34900 100644 --- a/zdm/loading.py +++ b/zdm/loading.py @@ -219,7 +219,12 @@ def surveys_and_grids(init_state=None, alpha_method=1, sdir=None, edir=None, rand_DMG=False, discard_empty=False, state_dict=None, - survey_dict=None,verbose=False): + survey_dict=None, verbose=False, + opdir='', + bPosNum=0, + cluster=False, + clusterRedshift=np.nan, + lensing=False): """ Load up a survey and grid for a real dataset Args: @@ -247,6 +252,12 @@ def surveys_and_grids(init_state=None, alpha_method=1, discard_empty (bool, optional): If true, does not calculate empty surveys (mostly for after latitude cuts) survey_dict (dict,None): list of survey metadata and values to apply + cluster (bool, optional): + If True, apply cluster (host-galaxy cluster) modelling. Defaults to False. + clusterRedshift (float, optional): + Redshift of the cluster. Defaults to np.nan (not used unless cluster=True). + lensing (bool, optional): + If True, apply lensing corrections. Defaults to False. Raises: IOError: [description] @@ -285,8 +296,11 @@ def surveys_and_grids(init_state=None, alpha_method=1, for survey_name in survey_names: # print(f"Initializing {survey_name}") s = survey.load_survey(survey_name, state, dmvals, zvals, - NFRB=NFRB, sdir=sdir, edir=edir, - rand_DMG=rand_DMG,survey_dict=survey_dict) + NFRB=NFRB, sdir=sdir, edir=edir, + rand_DMG=rand_DMG, survey_dict=survey_dict, + opdir=opdir,bPosNum=bPosNum, + cluster=cluster, clusterRedshift=clusterRedshift, + lensing=lensing) if discard_empty == False or s.NFRB != 0: # Check necessary parameters exist if considering repeaters @@ -302,7 +316,8 @@ def surveys_and_grids(init_state=None, alpha_method=1, # generates zdm grid grids = misc_functions.initialise_grids( - surveys, zDMgrid, zvals, dmvals, state, wdist=True, repeaters=repeaters) + surveys, zDMgrid, zvals, dmvals, state, wdist=True, repeaters=repeaters, + opdir=opdir, bPosNum=bPosNum, cluster=cluster, clusterRedshift=clusterRedshift) if verbose: print("Initialised grids") diff --git a/zdm/misc_functions.py b/zdm/misc_functions.py index 4fdfec75..b1b48672 100644 --- a/zdm/misc_functions.py +++ b/zdm/misc_functions.py @@ -41,6 +41,7 @@ from zdm import repeat_grid as zdm_repeat_grid from zdm import pcosmic from zdm import parameters +import magnificationMapper def get_w_tau_dist(grid,norm=True): @@ -1515,6 +1516,10 @@ def initialise_grids( state: parameters.State, wdist=True, repeaters=False, + opdir: str=None, + bPosNum: str=None, + cluster=False, + clusterRedshift = np.nan, ): """ For a list of surveys, construct a zDMgrid object wdist indicates a distribution of widths in the survey, @@ -1538,21 +1543,27 @@ def initialise_grids( # generates a DM mask # creates a mask of values in DM space to convolve with the DM grid - mask = pcosmic.get_dm_mask( - dmvals, (state.host.lmean, state.host.lsigma), zvals, plot=False - ) + if cluster: + hostMask = pcosmic.get_dm_mask( + dmvals, (state.host.lmean, state.host.lsigma), zvals, plot=True + ) + else: + mask = pcosmic.get_dm_mask( + dmvals, (state.host.lmean, state.host.lsigma), zvals, plot=True + ) grids = [] for survey in surveys: prev_grid = None - # print(f"Working on {survey.name}") + if cluster: + mask = pcosmic.get_cluster_dm_mask(survey, opdir, bPosNum, dmvals, zvals, hostMask, clusterRedshift) if repeaters: grid = zdm_repeat_grid.repeat_Grid( - survey, copy.deepcopy(state), zDMgrid, zvals, dmvals, mask, wdist, prev_grid=prev_grid + survey, copy.deepcopy(state), zDMgrid, zvals, dmvals, mask, wdist, prev_grid, cluster, clusterRedshift, bPosNum, opdir ) else: grid = zdm_grid.Grid( - survey, copy.deepcopy(state), zDMgrid, zvals, dmvals, mask, wdist, prev_grid=prev_grid + survey, copy.deepcopy(state), zDMgrid, zvals, dmvals, mask, wdist, prev_grid, cluster, clusterRedshift, bPosNum, opdir ) grids.append(grid) diff --git a/zdm/pcosmic.py b/zdm/pcosmic.py index 7f52dd56..26b62ba7 100644 --- a/zdm/pcosmic.py +++ b/zdm/pcosmic.py @@ -45,7 +45,11 @@ from frb.dm import igm from zdm import cosmology as cos from zdm import parameters - +from astropy.io import fits +from astropy import wcs +from astropy import units as u +import magnificationMapper +import scipy # from zdm import c_code import scipy as sp @@ -529,6 +533,25 @@ def plot_mean(zvals, saveas, title="Mean DM"): plt.show() plt.close() +def get_cluster_dm_mask(survey, opdir, bPosNum, dmvals, zvals, mask, clusterRedshift): + + DMThresh = np.load(opdir+'DMThresh.npy') + + new_mask = np.zeros([mask.shape[0], mask.shape[1], survey.meta["NBINS"]]) + for j in range(mask.shape[0]): + formatted_redshift = "{:03.2f}".format(zvals[j]) + pdms = np.load(opdir+'pdms_BP_'+str(bPosNum)+str(formatted_redshift)+'.npy') + for i in range(survey.meta["NBINS"]): + if np.round(zvals[j],2) > clusterRedshift: + if np.sum(np.isnan(pdms[:,i]))==len(pdms[:,i]): + new_mask[j,:,i] = mask[j,:] + else: + interpFunc = scipy.interpolate.interp1d(DMThresh[:-1], pdms[:,i], bounds_error=False, fill_value=0) + clusterConv = interpFunc(dmvals) + new_mask[j,:,i] = np.convolve(mask[j,:],clusterConv/np.sum(clusterConv), mode='full')[:mask.shape[1]] + else: + new_mask[j,:,i] = mask[j,:] + return new_mask def get_dm_mask(dmvals, params, zvals=None, plot=False): """Generate a convolution kernel for host galaxy DM contribution. diff --git a/zdm/scripts/beamMagDist.py b/zdm/scripts/beamMagDist.py new file mode 100644 index 00000000..91e13d6d --- /dev/null +++ b/zdm/scripts/beamMagDist.py @@ -0,0 +1,38 @@ +from magnificationMapper import normalisedLensFuncsAcrossBeam +from astropy.io import fits +from astropy import wcs +import astropy +from astropy import units as u +from astropy import constants as const + +import numpy as np +from zdm import survey +from matplotlib import pyplot as plt + + + + +magni = fits.getdata('hlsp_frontier_model_macs0717_bradac_v1_z01-magnif.fits') +info = fits.getheader('hlsp_frontier_model_macs0717_bradac_v1_z01-magnif.fits') +proj = wcs.WCS(info) + +xMagni = np.meshgrid(np.arange(0,len(magni[:,0]),1), np.arange(0,len(magni[0,:]),1)) +tempCoords = proj.array_index_to_world_values(xMagni[0], xMagni[1]) +relBeamPositions = np.load('relBeamPos.npy') +surveyName = 'CHORD' +ratesArr=np.zeros(len(relBeamPositions[:,0])) + +for i in range(len(relBeamPositions[:,0])): + print('---Beam Pos:', i) + name = 'CHORD_BeamPos_'+str(i) + mux, pmux = normalisedLensFuncsAcrossBeam(48*u.m, 900*u.MHz, 1e-3, 10, np.array([np.mean(tempCoords[0])+relBeamPositions[i,0], np.mean(tempCoords[1])+relBeamPositions[i,1]]), proj, magni, name) + np.save('mux_BP_'+str(i), np.log10(mux)) + np.save('pmux_BP_'+str(i), pmux) + np.save('mus', np.log10(mux)) + np.save('pmus', pmux) + fig, ax = plt.subplots() + ax.plot(np.log10(mux), np.log10(pmux)) + ax.set_xlabel('$\log_{10}\\mu$') + ax.set_ylabel('$\log_{10}p(\\mu)$') + fig.savefig('beamMagDist/probs'+format(i, '02d')) + plt.close() diff --git a/zdm/scripts/for_mawson.py b/zdm/scripts/for_mawson.py new file mode 100644 index 00000000..f74090b5 --- /dev/null +++ b/zdm/scripts/for_mawson.py @@ -0,0 +1,87 @@ +""" +This script creates zdm grids and plots localised FRBs + +It can also generate a summed histogram from all CRAFT data + +""" +import os + +from zdm import cosmology as cos +from zdm import misc_functions +from zdm import parameters +from zdm import survey +from zdm import pcosmic +from zdm import iteration as it +from zdm.craco import loading +from zdm import io + +import numpy as np +from zdm import survey +from matplotlib import pyplot as plt + + +def main(): + + # in case you wish to switch to another output directory + #opdir = "Localised_FRBs/" + opdir = "CHORD/" + + if not os.path.exists(opdir): + os.mkdir(opdir) + + # Initialise surveys and grids + + # The below is for private, unpublished FRBs. You will NOT see this in the repository! + sdir = "../data/Surveys/" + name = 'CHORD' + + # specifies state, updates variables according to H0 + # best fit, but with Emax extended as per Ryder et al +# state = parameters.State() +# state.energy.lEmax = 41.63 +# state.energy.gamma = -0.948 +# state.energy.alpha = -1.03 +# state.FRBdemo.sfr_n = 1.15 +# state.host.lsigma = 0.57 +# state.host.lmean = 2.22 +# state.FRBdemo.lC = 1.443 + state = parameters.State() + state.energy.lEmax = 41.38 + state.energy.gamma = -1.3 + state.energy.alpha = -1.39 + state.FRBdemo.sfr_n = 1.0 + state.host.lsigma = 0.57 + state.host.lmean = 2.22 + state.FRBdemo.lC = 4.86 + state.energy.luminosity_function=0 + + + s,g = loading.survey_and_grid(survey_name=name, + NFRB=None,sdir=sdir,init_state=state) + + np.save('zvals', g.zvals) + np.save('dmvals', g.dmvals) + + np.save('ratesUnlensed', g.rates) + + FRB_rate_per_day = np.sum(g.rates) * 10**g.state.FRBdemo.lC + print("Rate of FRBs per day is ",FRB_rate_per_day) + FRB_rate_per_day = np.sum(g.rates[g.zvals>1.0,:]) * 10**g.state.FRBdemo.lC + print("Rate of FRBs per day with z > 1.0 is ",FRB_rate_per_day) + + misc_functions.plot_grid_2( + g.rates, + g.zvals, + g.dmvals, + name=opdir + name + ".pdf", + norm=3, + log=True, + label="$\\log_{10} p({\\rm DM}_{\\rm EG},z)$ [a.u.]", + project=False, + zmax=4, + Aconts=[0.01, 0.1, 0.5], + DMmax=4000 + ) # + + +main() diff --git a/zdm/scripts/forecastZDM.py b/zdm/scripts/forecastZDM.py new file mode 100644 index 00000000..0bdb318b --- /dev/null +++ b/zdm/scripts/forecastZDM.py @@ -0,0 +1,122 @@ +""" +This script creates zdm grids and plots localised FRBs + +It can also generate a summed histogram from all CRAFT data + +""" +import os + +from zdm import cosmology as cos +from zdm import misc_functions +from zdm import parameters +from zdm import survey +from zdm import pcosmic +from zdm import iteration as it +from zdm import loading +from zdm import io +from magnificationMapper import normalisedLensFuncsAcrossBeam +from astropy.io import fits +from astropy import wcs +import astropy +from astropy import units as u +from astropy import constants as const + +import numpy as np +from zdm import survey, figures +from matplotlib import pyplot as plt + + +def pullTrigger(clusterRedshift, name, gamma, n, opdir=''): + + formatted_cluster_redshift = "{:03.2f}".format(clusterRedshift) + formatted_energy_index = "{:03.2f}".format(gamma) + opdir = opdir+"/"+name+'/z'+formatted_cluster_redshift+'/'+formatted_energy_index+'/' + print(opdir) + #opdir = "/arc/projects/chime_frb/msammons/CHIME/ClusterLensed/"+name+'/z'+formatted_cluster_redshift+'/' + + # Initialise surveys and grids + + # The below is for private, unpublished FRBs. You will NOT see this in the repository! + sdir = "../data/Surveys/" + + # specifies state, updates variables according to H0 + # best fit, but with Emax extended as per Ryder et al +# state = parameters.State() +# state.energy.lEmax = 41.63 +# state.energy.gamma = -0.948 +# state.energy.alpha = -1.03 +# state.FRBdemo.sfr_n = 1.15 +# state.host.lsigma = 0.57 +# state.host.lmean = 2.22 +# state.FRBdemo.lC = 1.443 + state = parameters.State() + state.energy.lEmax = 41.63 + state.energy.lEmin = 30.0 + state.energy.gamma = gamma + state.energy.alpha = 1.03 + state.FRBdemo.sfr_n = n + state.host.lsigma = 0.57 + state.host.lmean = 2.22 + state.FRBdemo.lC = 2.3-9 + state.energy.luminosity_function=4 + state.FRBdemo.alpha_method=1 + #state.FRBdemo.source_evolution=0 +# state = parameters.State() +# state.energy.lEmax = 41.42 +# state.energy.gamma = -1.16 +# state.energy.alpha = 0.92 +# state.FRBdemo.sfr_n = 0.91 +# state.host.lsigma = 0.46 +# state.host.lmean = 2.02 +# state.FRBdemo.lC = 2.0 +# state.energy.luminosity_function=2 +# state.FRBdemo.alpha_method=1 + + + cluster=True + lensing =True + + #relBeamPositions = np.load('relBeamPos.npy') #relative to magni + relBeamPositions = np.array([[0,0]]) + ratesArr=np.zeros([len(relBeamPositions[:,0]), 2]) + + + for i in range(len(relBeamPositions[:,0])): + print('---Beam Pos:', i) + formatted_number = "{:02d}".format(i) + surveyName = 'CHIME_SynthBeam' + bPosNum = "{:02d}".format(i) + s,gSet = loading.surveys_and_grids(survey_names=[surveyName], opdir=opdir, bPosNum=bPosNum, + NFRB=None,sdir=sdir,init_state=state, cluster=cluster, + clusterRedshift=clusterRedshift, lensing=lensing) + + print('gSet lengths:', len(gSet)) + g = gSet[0] + np.save(opdir+'rates_BP_'+str(formatted_number)+'SourceFunc1thresh5NEWZDMSFRn'+str("{:.2f}".format(n)), g.rates*10**g.state.FRBdemo.lC) + + + FRB_rate_per_day = np.sum(g.rates) * 10**g.state.FRBdemo.lC + print("Rate of FRBs per day is ",FRB_rate_per_day) + ratesArr[i,0] = FRB_rate_per_day + FRB_rate_per_day = np.sum(g.rates[g.zvals>1.0,:]) * 10**g.state.FRBdemo.lC + print("Rate of FRBs per day with z > 1.0 is ",FRB_rate_per_day) + + ratesArr[i,1] = FRB_rate_per_day + np.save(opdir+surveyName+'RatesArrSourceFunc1thresh5NEWZDMSFRn'+str("{:.2f}".format(n)), ratesArr) + + print('SAVING TO', opdir + surveyName + 'SourceFunc1thresh5NEWZDMSFRn'+str("{:.2f}".format(n))) + figures.plot_grid( + g.rates, + g.zvals, + g.dmvals, + name=opdir + surveyName + 'SourceFunc1thresh5NEWZDMSFRn'+str("{:.2f}".format(n))+'.pdf', + norm=3, + log=True, + label="$\\log_{10} p({\\rm DM}_{\\rm EG},z)$ [a.u.]", + project=False, + zmax=4, + Aconts=[0.01, 0.1, 0.5], + DMmax=4000 + ) # + + diff --git a/zdm/scripts/initialiseClusterContributions.py b/zdm/scripts/initialiseClusterContributions.py new file mode 100644 index 00000000..220e0e08 --- /dev/null +++ b/zdm/scripts/initialiseClusterContributions.py @@ -0,0 +1,186 @@ +""" +This script creates zdm grids and plots localised FRBs + +It can also generate a summed histogram from all CRAFT data + +""" +import os +from magnificationMapper import normalisedLensFuncsAcrossBeam, clusterDMFuncAcrossBeam, mapRescaler +from astropy.io import fits +from astropy import wcs +import astropy +from astropy import units as u +from astropy import constants as const +from astropy.cosmology import Planck18 as cosmo + +import numpy as np +from matplotlib import pyplot as plt + + +def initialise_Magni(opdir,clusterRedshift, name, clusterNeFile, trueClusterRedshift, sourceMagniArr=False): + + # in case you wish to switch to another output directory + #opdir = "Localised_FRBs/" + formatted_cluster_redshift = "{:03.2f}".format(clusterRedshift) + opdir = opdir+"/"+name+'/z'+formatted_cluster_redshift+'/' + print(opdir) + + if not os.path.exists(opdir): + os.makedirs(opdir) + + dishDiam = 80*u.m + fbar = 600*u.MHz + bThresh = 1e-3 + bbins = 100 + + fileList = [name+'_kappa.fits', name+'_gamma.fits', clusterNeFile] + mapRescaler(opdir, fileList, trueClusterRedshift, clusterRedshift) + + kappa = fits.getdata(opdir+fileList[0]) + gamma = fits.getdata(opdir+fileList[1]) + info = fits.getheader(opdir+fileList[0]) + proj = wcs.WCS(info) + xMagni = np.meshgrid(np.arange(0,len(kappa[:,0]),1), np.arange(0,len(kappa[0,:]),1)) + tempCoords = proj.array_index_to_world_values(xMagni[0], xMagni[1]) + + infoNe = fits.getheader(opdir+fileList[2]) + projNe = wcs.WCS(infoNe) + ne = fits.getdata(opdir+fileList[2]) + + cluster=True + lensing =True + + #relBeamPositions = np.load('relBeamPos.npy') #relative to magni + relBeamPositions = np.array([[0,0]]) + ratesArr=np.zeros(len(relBeamPositions[:,0])) + zvals = np.load('zvals.npy') + mux = 10**(np.arange(-3,2,0.02)+0.05) + np.save(opdir+'mux', np.log10(mux)) + scatThresh = 10**np.arange(-4,3,0.02) + xProbScat = scatThresh[:-1]*10**(np.diff(np.log10(scatThresh))[0]/2) + np.save(opdir+'xProbScat', (xProbScat)) + DMThresh = np.arange(0,15000,200) + np.save(opdir+'DMThresh', DMThresh) + + + for i in range(len(relBeamPositions[:,0])): + bPos = np.array([np.mean(tempCoords[0])+relBeamPositions[i,0], np.mean(tempCoords[1])+relBeamPositions[i,1]]) + for j in range(len(zvals)): + formatted_number = "{:02d}".format(i) + formatted_redshift = "{:03.2f}".format(zvals[j]) + surveyName = 'CHIME_BeamPos_'+str(formatted_number) + + if zvals[j] > clusterRedshift: + if sourceMagniArr: + pixCoordsArr= np.load(opdir+("{:03.2f}".format(zvals[j]))+'_sourceCoordsArr.npy') + magni = np.load(opdir+("{:03.2f}".format(zvals[j]))+'_sourceCoordsArr.npy') + + else: + rescaleFactor = ((cosmo.angular_diameter_distance(clusterRedshift)**2/(cosmo.angular_diameter_distance(trueClusterRedshift)**2))*cosmo.angular_diameter_distance_z1z2(clusterRedshift, zvals[j])*cosmo.angular_diameter_distance(trueClusterRedshift)/cosmo.angular_diameter_distance(zvals[j])/cosmo.angular_diameter_distance(clusterRedshift)).value + pixCoordsArr = xMagni + magni = 1/np.abs((1-kappa*rescaleFactor)**2-(gamma*rescaleFactor)**2) + magni[magni>=100] = 100 + + + pmux, magni, xMagni = normalisedLensFuncsAcrossBeam(dishDiam, fbar, bThresh, bbins, bPos, proj, xMagni, magni, pixCoordsArr, opdir+surveyName, sourceMagniArr=sourceMagniArr, muThresh = mux) + + np.save(opdir+'pmux_BP_'+str(formatted_number)+str(formatted_redshift), pmux) + + +def initialise_DM_Scattering(opdir, clusterRedshift, name, clusterNeFile, trueClusterRedshift, energyIndex, sourceMagniArr=False): + + # in case you wish to switch to another output directory + #opdir = "Localised_FRBs/" + formatted_cluster_redshift = "{:03.2f}".format(clusterRedshift) + formatted_energy_index = "{:03.2f}".format(energyIndex) + opdir = opdir+"/"+name+'/z'+formatted_cluster_redshift+'/'+formatted_energy_index+'/' + print(opdir) + + if not os.path.exists(opdir): + os.makedirs(opdir) + + dishDiam = 80*u.m + fbar = 600*u.MHz + bThresh = 1e-3 + bbins = 100 + + fileList = [name+'_kappa.fits', name+'_gamma.fits', clusterNeFile] + mapRescaler(opdir, fileList, trueClusterRedshift, clusterRedshift) + + kappa = fits.getdata(opdir+fileList[0]) + gamma = fits.getdata(opdir+fileList[1]) + info = fits.getheader(opdir+fileList[0]) + proj = wcs.WCS(info) + xMagni = np.meshgrid(np.arange(0,len(kappa[:,0]),1), np.arange(0,len(kappa[0,:]),1)) + tempCoords = proj.array_index_to_world_values(xMagni[0], xMagni[1]) + + infoNe = fits.getheader(opdir+fileList[2]) + projNe = wcs.WCS(infoNe) + ne = fits.getdata(opdir+fileList[2]) + + cluster=True + lensing =True + + #relBeamPositions = np.load('relBeamPos.npy') #relative to magni + relBeamPositions = np.array([[0,0]]) + ratesArr=np.zeros(len(relBeamPositions[:,0])) + zvals = np.load('zvals.npy') + mux = 10**(np.arange(-3,2,0.02)+0.05) + np.save(opdir+'mux', np.log10(mux)) + scatThresh = 10**np.arange(-4,3,0.02) + xProbScat = scatThresh[:-1]*10**(np.diff(np.log10(scatThresh))[0]/2) + np.save(opdir+'xProbScat', (xProbScat)) + DMThresh = np.arange(0,15000,200) + np.save(opdir+'DMThresh', DMThresh) + + + for i in range(len(relBeamPositions[:,0])): + bPos = np.array([np.mean(tempCoords[0])+relBeamPositions[i,0], np.mean(tempCoords[1])+relBeamPositions[i,1]]) + for j in range(len(zvals)): + formatted_number = "{:02d}".format(i) + formatted_redshift = "{:03.2f}".format(zvals[j]) + surveyName = 'CHIME_BeamPos_'+str(formatted_number) + + if zvals[j] > clusterRedshift: + if sourceMagniArr: + pixCoordsArr= np.load(opdir+("{:03.2f}".format(zvals[j]))+'_sourceCoordsArr.npy') + magni = np.load(opdir+("{:03.2f}".format(zvals[j]))+'_sourceCoordsArr.npy') + + else: + rescaleFactor = (cosmo.angular_diameter_distance_z1z2(clusterRedshift, zvals[j])*cosmo.angular_diameter_distance(trueClusterRedshift)/cosmo.angular_diameter_distance(zvals[j])/cosmo.angular_diameter_distance(clusterRedshift)).value + pixCoordsArr = xMagni + magni = 1/np.abs((1-kappa*rescaleFactor)**2-(gamma*rescaleFactor)**2) + magni[magni>=100] = 100 + + if sourceMagniArr: + rawWeights = (1/magni)**(energyIndex) + else: + rawWeights = 1/magni*(1/magni)**(energyIndex) + + tempCoords = proj.array_index_to_world_values(xMagni[0], xMagni[1]) + + probScat, fractionUnscattered, pdms = clusterDMFuncAcrossBeam( + D = dishDiam, + freq = fbar, + thresh = bThresh, + nbins = bbins, + bPos = bPos, + proj = projNe, + clusterRedshift = clusterRedshift, + z = zvals[j], + ne = ne, + name = opdir+'DM_BP_'+str(formatted_number), + weights = rawWeights, + imageProj = proj, + imageCoords = tempCoords, + DMThresh = DMThresh, + scatThresh = scatThresh + ) + np.save(opdir+'probScat_BP_'+str(formatted_number)+str(formatted_redshift), probScat) + np.save(opdir+'fractionUnscattered_BP_'+str(formatted_number)+str(formatted_redshift), fractionUnscattered) + np.save(opdir+'pdms_BP_'+str(formatted_number)+str(formatted_redshift),pdms) + else: + # fill in based on other else conditions + np.save(opdir+'probScat_BP_'+str(formatted_number)+str(formatted_redshift), np.zeros([len(xProbScat),bbins])) + np.save(opdir+'fractionUnscattered_BP_'+str(formatted_number)+str(formatted_redshift), np.ones(bbins)) + np.save(opdir+'pdms_BP_'+str(formatted_number)+str(formatted_redshift),np.zeros([len(DMThresh[:-1]),bbins])) diff --git a/zdm/scripts/magnificationMapper.py b/zdm/scripts/magnificationMapper.py new file mode 100644 index 00000000..f10816b7 --- /dev/null +++ b/zdm/scripts/magnificationMapper.py @@ -0,0 +1,245 @@ +import scipy.signal +import scipy.interpolate +import scipy.integrate +import matplotlib.pyplot as plt +import numpy as np +from astropy.cosmology import Planck18 as cosmo +from astropy.cosmology import LambdaCDM +from astropy import units as u +from astropy import constants as const +import matplotlib.pyplot as plt +import numpy as np +#import dynspec +#import pygedm +import astropy.coordinates as c +from astropy.io import fits +from astropy import wcs +import astropy +from astropy.convolution import Gaussian2DKernel + +def offSetBeamGains(bPos, imageCoords, beamSigma): + xOffset = imageCoords[0] - bPos[0] + yOffset = imageCoords[1] - bPos[1] + radiusOffset = (xOffset**2+yOffset**2)**0.5 + gains = BFG(radiusOffset, beamSigma) + return gains + +def BFG(x, sigma): + return np.exp(-1/2*(x**2/(sigma**2))) + + +def logSpaceIntegrand(logmu, func, funcArgs, base): + """Useful for evaluating integrals of func in log space""" + return func(base**logmu, funcArgs)*base**logmu*np.log(base) + +def unnormalisedLensFuncAtSubBeam(log10b, dlog10b, OmegaB, imagePlaneBGains, bGains, pixRes, magniArr, sourceMagniArr, muThresh): + #OmegaB in arcminutes^2, same as pixRes + inBeam = np.abs(np.log10(bGains)-log10b)0: + print(log10b, dlog10b, 'ALERT: somethings dead wrong !!!!!!!!!!!!!!!!!!') + return np.nan + if np.sum(inBeam)>0: + gtrMu = np.zeros(len(muThresh)) + for i in range(len(muThresh)): + gtrMu[i] = np.sum(magniArr[inBeam]>=muThresh[i]) + modelledArea = np.sum(planeInBeam)*(pixRes[0]*pixRes[1]) + numUnmodelledCells = (OmegaB - modelledArea)/(pixRes[0]*pixRes[1]) + extra1SWhere = muThresh<1 + gtrMu[extra1SWhere] = gtrMu[extra1SWhere]+numUnmodelledCells + if sourceMagniArr: + probUN = (-1*np.diff((gtrMu))/np.diff(muThresh)) + else: + probUN = (-1*np.diff((gtrMu))/np.diff(muThresh)/muThresh[:-1]) + + # smoothingKernel = scipy.signal.windows.gaussian(len(muThresh[:-1]),0.05/(np.mean(np.diff(np.log10(muThresh))))) + # probUNSmooth= np.convolve(probUN,smoothingKernel, mode='same')/np.sum(smoothingKernel) + interpFunc = scipy.interpolate.interp1d(np.log10(muThresh[:-1])+np.diff(np.log10(muThresh))[0], probUN, bounds_error=False, fill_value=0) + else: + interpFunc = None + return interpFunc + +def mapRescaler(opdir, fileList, zTrue, zNew): + scale_factor = (cosmo.angular_diameter_distance(zNew)/cosmo.angular_diameter_distance(zTrue)).value + for i in range(len(fileList)): + hdulist = fits.open(fileList[i]) + header = hdulist[0].header + header['CDELT1'] *= scale_factor + header['CDELT2'] *= scale_factor + if 'CD1_1' in header: + header['CD1_1'] *= scale_factor + header['CD1_2'] *= scale_factor + header['CD2_1'] *= scale_factor + header['CD2_2'] *= scale_factor + fits.writeto(opdir+fileList[i],hdulist[0].data, header, overwrite=True) + hdulist.close() + + +def mapWidener(magni, wideningFrac): + smoothKernel = Gaussian2DKernel(int(magni.shape[0]/100)) + smoothMagni = scipy.signal.fftconvolve(np.log10(magni), smoothKernel, mode='valid') + croppedEachSide = np.floor((np.asarray(magni.shape) - np.asarray(smoothMagni.shape))/2) + interpFunc = scipy.interpolate.RegularGridInterpolator((np.arange(int(croppedEachSide[0]), (magni.shape[0] - int(croppedEachSide[0])),1), np.arange(int(croppedEachSide[1]), (magni.shape[1] - int(croppedEachSide[1])),1)), (smoothMagni), bounds_error=False, fill_value=None) + wA = int(magni.shape[0]*wideningFrac/2) + x = np.meshgrid(np.arange(-wA,magni.shape[0]+wA,1), np.arange(-wA,magni.shape[1]+wA,1)) + expandedMagni = interpFunc((x[1].flatten(), x[0].flatten())) + finalMagni = (expandedMagni.reshape(np.asarray(magni.shape)+wA*2)) + finalMagni[finalMagni<0] = 0 + finalMagni[:wA,:wA] = np.mean(finalMagni[wA:-wA,:wA]) + finalMagni[-wA:, :wA] = np.mean(finalMagni[-wA:, wA:-wA]) + finalMagni[-wA:,-wA:] = np.mean(finalMagni[wA:-wA, -wA:]) + finalMagni[:wA, -wA:] = np.mean(finalMagni[:wA, wA:-wA]) + completeMagni = finalMagni + completeMagni[wA:-wA,wA:-wA] = np.log10(magni) + return 10**completeMagni, x + +def normalisedLensFuncsAcrossBeam(D, freq, thresh, nbins, bPos, proj, x, magniArr, pixCoordsArr, name, sourceMagniArr=False, muThresh = 10**(np.arange(-3,2,0.02)+0.05)): + FWHM = 1.22*(const.c/(freq))/D + beamSigma=(FWHM/2.)*(2*np.log(2))**-0.5 + dlnb=-np.log(thresh)/nbins + log10min=np.log10(thresh) + dlog10b=log10min/nbins + log10b=(np.arange(nbins)+0.5)*dlog10b + OmegaB= (2*np.pi*dlnb*(beamSigma*180/np.pi*60)**2).decompose().value + pixRes = np.abs(np.diag(proj.pixel_scale_matrix*60)) + + if sourceMagniArr: + pixCoords = [] + pixCoords.append(pixCoordsArr[:,0]) + pixCoords.append(pixCoordsArr[:,1]) + else: + dataEdge = np.mean(np.concatenate((magniArr[:,0], magniArr[:,-1], magniArr[0,:], magniArr[-1,:]))) + count = 0 + tempMagni = magniArr.copy() + #while (dataEdge -1) > 0.1: + # tempMagni,x = mapWidener(magniArr, 0.5+0.1*count) + # dataEdge = np.mean(np.concatenate((tempMagni[:,0], tempMagni[:,-1], tempMagni[0,:], tempMagni[-1,:]))) + # count = count+1 + # print('trapped forever', count) + magniArr = tempMagni + pixCoords = x + + imageCoords = proj.array_index_to_world_values(pixCoords[0], pixCoords[1]) + bGains = offSetBeamGains(bPos, imageCoords, beamSigma.decompose().value*180/np.pi) + imagePlaneCoords = proj.array_index_to_world_values(x[0], x[1]) + imagePlaneBGains = offSetBeamGains(bPos, imagePlaneCoords, beamSigma.decompose().value*180/np.pi) + + + +# fig = plt.figure() +# ax = plt.subplot(111, projection=proj) +# tower = np.zeros(bGains.shape) +# for i in range(len(log10b)): +# gainLevel = np.abs(np.log10(bGains)-log10b[i])0: + gtrDM = np.zeros(len(DMThresh)) + gtrScat = np.zeros([len(scatThresh)]) + probScat = np.zeros([len(scatThresh)-1]) + for i in range(len(DMThresh)): + gtrDM[i] = np.sum((neWeightedHist[0]*inBeam)*((1e6/(1+clusterRedshift)*ne)>=DMThresh[i])) + for i in range(len(scatThresh)): + if z>clusterRedshift: + scat = (4.1e-5/(1+clusterRedshift)*(lam/1)**4*((cosmo.angular_diameter_distance(clusterRedshift)*cosmo.angular_diameter_distance_z1z2(clusterRedshift,z)/cosmo.angular_diameter_distance(z)).value/1e3)*(8.4e-13*(ne/1e-4)**2*3.08567758e+22/((1+clusterRedshift)**2)/1e12)*(2.06264806e+9)**(1/3)*1e3) + if(np.amin(scat)0): + print('WARNING: Scattering outside threshold') + print('z = ', z, np.amin(scat), np.amin(scatThresh)) + break + gtrScat[i] = np.sum((scat>=scatThresh[i])*neWeightedHist[0]*inBeam) + if i==0: + gtrScat[0]=np.sum((scat>=0)*neWeightedHist[0]*inBeam) + + probScat[:] = (-1*np.diff(gtrScat[:])/(gtrScat[0])) + else: + probScat[:] = 0 + + modelledArea = np.sum(inBeam_2)*(pixRes[0]*pixRes[1]) + + + numUnmodelledCells = (OmegaB - modelledArea)/(pixRes[0]*pixRes[1]) + if numUnmodelledCells < 0: + numUnmodelledCells = 0 + + fractionUnscattered = (numUnmodelledCells+DMLessWeights)/(np.sum(weights*inBeam_2)+numUnmodelledCells) + gtrDM[0] = gtrDM[0]+numUnmodelledCells+DMLessWeights + probUN = (-1*np.diff((gtrDM))/np.diff(DMThresh)) + #interpFunc = scipy.interpolate.interp1d((DMThresh[:-1]), probUN, bounds_error=False, fill_value=0) + else: + gtrScat = np.ones(len(scatThresh[:-1]))*np.nan + numUnmodelledCells = np.nan + #interpFunc = None + probUN = np.ones(len(DMThresh[:-1]))*np.nan + probScat = np.zeros([len(scatThresh[:-1])]) + fractionUnscattered = 1 + return probUN, probScat, fractionUnscattered + diff --git a/zdm/scripts/runClusterInitialise.py b/zdm/scripts/runClusterInitialise.py new file mode 100644 index 00000000..262d7024 --- /dev/null +++ b/zdm/scripts/runClusterInitialise.py @@ -0,0 +1,13 @@ +import initialiseClusterContributions as iCC +import os + + + +print('initialising magnifications') +iCC.initialise_Magni('/arc/projects/chime_frb/msammons/clusterLensing/testZDM', 0.18, 'Abell2218', 'RadialABELL_2218.fits', 0.18) + +print('initialising DM and Scattering') +iCC.initialise_DM_Scattering('/arc/projects/chime_frb/msammons/clusterLensing/testZDM', 0.18, 'Abell2218', 'RadialABELL_2218.fits', 0.18 ,-1.0) + + + diff --git a/zdm/scripts/run_forecastZDM.py b/zdm/scripts/run_forecastZDM.py new file mode 100644 index 00000000..0dec8ad3 --- /dev/null +++ b/zdm/scripts/run_forecastZDM.py @@ -0,0 +1,13 @@ +from forecastZDM import pullTrigger +import numpy as np + +args = [] +gammas = np.array([-1]) +ns = np.array([1]) +print('pulling the trigger in serial') +for i in range(len(gammas)): + for j in range(len(ns)): + print('----->>> Executing: gamma:'+str(gammas[i])+', n:'+str(ns[j])) + pullTrigger(0.18, 'Abell2218', gammas[i], ns[j], opdir='/arc/projects/chime_frb/msammons/clusterLensing/testZDM') + +print('fin') diff --git a/zdm/scripts/sampleZDE.py b/zdm/scripts/sampleZDE.py new file mode 100644 index 00000000..99609823 --- /dev/null +++ b/zdm/scripts/sampleZDE.py @@ -0,0 +1,155 @@ +""" +This script creates zdm grids and plots localised FRBs + +It can also generate a summed histogram from all CRAFT data + +""" +import os + +from zdm import cosmology as cos +from zdm import misc_functions +from zdm import parameters +from zdm import survey +from zdm import pcosmic +from zdm import iteration as it +from zdm.craco import loading +from zdm import io +import astropy +from astropy import units as u +from astropy import constants as const +import pickle +import numpy as np +from zdm import survey +from matplotlib import pyplot as plt +import scipy.stats +from astropy.cosmology import Planck18 as cosmo +from scipy.interpolate import interp1d + +def renormalise(enorm,emin,emax,gamma): + oldNorm = (emin**gamma-emax**gamma) + newNorm = (enorm**gamma-emax**gamma) + renormFactor = oldNorm/newNorm #multiply existing rates by this + return renormFactor + +def F_to_E(F,z,alpha=0, bandwidth=1e9, Fobs=1.3e9, Fref=1.3e9): + """ Converts a fluence in Jy ms to an energy in erg + Formula from Macquart & Ekers 2018 + Works with an array of z. + + Arguments are: + Fluence: of an FRB [Jy ms] + + Redshift: assumed redshift of an FRB producing the fluence F. + Standard cosmological definition [unitless] + + alpha: F(\nu)~\nu^-\alpha. Note that this is an internal definition. + The paper uses ^alpha, not ^-alpha. [unitless] + + Bandwidth: over which to integrate fluence [Hz] + + Fobs: the observation frequency [Hz] + + Fref: reference frequency at which FRB energies E are normalised. + It defaults to 1.3 GHz (ASKAP lat50, Parkes). + + Return value: energy [erg] + + """ + E=F*4*np.pi*(cosmo.luminosity_distance(z).value)**2/(1.+z)**(2.-alpha) + # now convert from dl in MPc and F in Jy ms + # 10^-26 from Jy to W per m2 per Hz + # 1e-3 from Jy ms to J per m2 per Hz + # (3.086e16 m in 1 pc x 10^6 Mpc)^2 for dl in m + # 1e7 from J to erg + # total factor is 9.523396e22 + E *= 9.523396e22*bandwidth + + # now corrects for reference frequency + # according to value of alpha + # effectively: if fluence was X at F0, it was X*(F0/Fref)**alpha at Fref + # i.e. if alpha is positive (stronger at low frequencies), we reduce E + # This acts to reduce the telescope threshold at higher frequencies + E *= (Fobs/Fref)**alpha + + return E + + +def rollFRBSample(Ugrid, N,renormEnergy=1e39): + renormF = renormalise(renormEnergy, 10**Ugrid.state.energy.lEmin,10**Ugrid.state.energy.lEmax, np.asarray([Ugrid.state.energy.gamma])) + zdmMesh = np.meshgrid(Ugrid.dmvals,Ugrid.zvals) + tempBase = Ugrid.rates*10**Ugrid.state.FRBdemo.lC*renormF[0]*1024 + tempInts = np.random.choice(np.arange(len(tempBase.flatten())), p=tempBase.flatten()/np.sum(tempBase),size=N) + es = np.linspace(Ugrid.state.energy.lEmin, Ugrid.state.energy.lEmax+2, 10000) + efracs = Ugrid.array_cum_lf(10**es, 10**Ugrid.state.energy.lEmin, 10**Ugrid.state.energy.lEmax, Ugrid.state.energy.gamma, Ugrid.use_log10) + efunc = interp1d(np.log10(efracs),es) + subfractions = np.random.uniform(size=N) + randZs = (zdmMesh[1].flatten())[tempInts] + randDMs = (zdmMesh[0].flatten())[tempInts] + Fth=5 + eth = F_to_E(Fth,randZs) + ethFrac = Ugrid.array_cum_lf(eth, 10**Ugrid.state.energy.lEmin, 10**Ugrid.state.energy.lEmax, Ugrid.state.energy.gamma) + randEs = efunc(np.log10((ethFrac)*subfractions)) + return randZs, randDMs, randEs + +def main(gamma, n, emax, N_frbs, save_grid=False): + # in case you wish to switch to another output directory + #opdir = "Localised_FRBs/" + #clusterRedshift = 0.38 + opdir = "./" + + # Initialise surveys and grids + + sdir = "../data/Surveys/" +# state = parameters.State() +# state.energy.lEmax = 42.63 +# state.energy.lEmin = 30 +# state.energy.gamma = gamma +# state.energy.alpha = 1.03 +# state.FRBdemo.sfr_n = n +# state.host.lsigma = 0.57 +# state.host.lmean = 2.22 +# #state.FRBdemo.lC = -0.49 +# state.FRBdemo.lC = 2.3-9 +# state.energy.luminosity_function=0 +# state.FRBdemo.alpha_method=1 + + state = parameters.State() + state.energy.lEmax = emax + state.energy.lEmin = 39 + state.energy.gamma = gamma + state.energy.alpha = 0 + state.FRBdemo.sfr_n = n + state.host.lsigma = 0.41 + state.host.lmean = 1.93 + state.FRBdemo.lC = 2.3-9 + #state.FRBdemo.lC = 2.3-9+1.31 + state.energy.luminosity_function=2 + state.FRBdemo.alpha_method=0 + + + surveyName = 'CHIME_SynthBeam' + s,g = loading.survey_and_grid(survey_name=surveyName, opdir=opdir, + NFRB=None,sdir=sdir,init_state=state) + + z,d,e = rollFRBSample(g, N_frbs) + + np.save(opdir+'KaitUnlensedthresh5Gamma'+str("{:.2f}".format(gamma))+'SFRn'+str("{:.2f}".format(n))+'Emax'+str("{:.2f}".format(emax))+'N'+str(N_frbs)+'sampleRedshift',z) + np.save(opdir+'KaitUnlensedthresh5Gamma'+str("{:.2f}".format(gamma))+'SFRn'+str("{:.2f}".format(n))+'Emax'+str("{:.2f}".format(emax))+'N'+str(N_frbs)+'sampleDM',d) + np.save(opdir+'KaitUnlensedthresh5Gamma'+str("{:.2f}".format(gamma))+'SFRn'+str("{:.2f}".format(n))+'Emax'+str("{:.2f}".format(emax))+'N'+str(N_frbs)+'sampleEnergy',e) + + + + # Save + if save_grid: + print('saving grid') + with open(opdir+'KaitGridUnlensedthresh5Gamma'+str("{:.2f}".format(gamma))+'SFRn'+str("{:.2f}".format(n))+'Emax'+str("{:.2f}".format(emax)), "wb") as f: + pickle.dump(g, f) + print('fin') + +if __name__ == "__main__": + gamma=-0.01 + nsfr = 0.96 + N_frbs=1000 + emax = 41.38 + main(gamma, nsfr, emax, N_frbs) + diff --git a/zdm/survey.py b/zdm/survey.py index aa27905a..162f5d3e 100644 --- a/zdm/survey.py +++ b/zdm/survey.py @@ -11,6 +11,8 @@ - Detection efficiency as a function of DM and width - DM budget calculations (Galactic, halo, host contributions) - Threshold calculations for FRB detection +- Optional cluster scattering modelling (via pre-computed per-beam scattering + probability distributions) Survey Definition Files ----------------------- @@ -82,6 +84,8 @@ class Survey: Survey metadata from file header. efficiencies : ndarray Detection efficiency grid as function of DM and width. + cluster : bool + Whether cluster scattering is applied to this survey. """ def __init__(self, state, survey_name: str, @@ -92,7 +96,12 @@ def __init__(self, state, survey_name: str, iFRB: int = 0, edir=None, rand_DMG=False, - survey_dict=None): + survey_dict=None, + opdir=None, + bPosNum=None, + cluster=False, + clusterRedshift=None, + lensing=False): """Initialize an FRB Survey. Parameters @@ -117,6 +126,22 @@ def __init__(self, state, survey_name: str, If True, randomize Galactic DM within uncertainty. Default False. survey_dict : dict, optional Override survey metadata parameters. + opdir : str, optional + Directory containing pre-computed cluster scattering files + (xProbScat.npy, probScat_BP_*.npy, fractionUnscattered_BP_*.npy). + Required if cluster=True. + bPosNum : int or str, optional + Beam position number used to select cluster scattering files. + Required if cluster=True. + cluster : bool, optional + If True, load and apply cluster scattering distributions. + Default False. + clusterRedshift : float, optional + Redshift of the cluster. Informational; not used directly in + calculations. + lensing : bool, optional + If True, apply lensing magnification (reserved for future use). + Default False. """ # Proceed self.state = state @@ -127,14 +152,13 @@ def __init__(self, state, survey_name: str, if zvals is not None: self.NZ = zvals.size self.edir = edir + self.lensing = lensing + self.clusterRedshift = clusterRedshift # Load up self.process_survey_file(filename, NFRB, iFRB, min_lat=state.analysis.min_lat, - dmg_cut=state.analysis.DMG_cut,survey_dict = survey_dict) + dmg_cut=state.analysis.DMG_cut, survey_dict=survey_dict) # Check if repeaters or not and set relevant parameters - # Now done in loading - # self.repeaters=False - # self.init_repeaters() # DM EG self.init_halo_coeffs() if rand_DMG: @@ -155,18 +179,57 @@ def __init__(self, state, survey_name: str, plot=False, thresh=beam_thresh) # tells the survey to use the beam file + # --- Cluster scattering setup --- + # Store cluster flag and load per-(z, beam) scattering distributions + # if cluster scattering is requested. + self.cluster = cluster + if cluster: + if opdir is None or bPosNum is None: + raise ValueError("opdir and bPosNum must be provided when cluster=True") + if zvals is None: + raise ValueError("zvals must be provided when cluster=True") + xProbScat = np.load(os.path.join(opdir, 'xProbScat.npy')) + NBINS = self.meta['NBINS'] + NZ = len(zvals) + probScat = np.zeros([len(xProbScat), NZ, NBINS]) + fractionUnscattered = np.zeros([NZ, NBINS]) + for j in range(NZ): + formatted_redshift = "{:03.2f}".format(zvals[j]) + probScat[:, j, :] = np.load( + os.path.join(opdir, + 'probScat_BP_' + str(bPosNum) + + str(formatted_redshift) + '.npy')) + fractionUnscattered[j, :] = np.load( + os.path.join(opdir, + 'fractionUnscattered_BP_' + str(bPosNum) + + str(formatted_redshift) + '.npy')) + self.xProbScat = xProbScat + self.probScat = probScat # shape [Nsc, NZ, NBINS] + self.fractionUnscattered = fractionUnscattered # shape [NZ, NBINS] + else: + self.xProbScat = None + self.probScat = None + self.fractionUnscattered = None + # -------------------------------- + # initialise scattering/width and resulting efficiences self.init_widths() self.calc_max_dm() - self.init_frb_bvals() #initial;ise weights for FRBs with known beam values + self.init_frb_bvals() #initialise weights for FRBs with known beam values self.init_frb_wvals() self.calc_max_dm() - def init_widths(self,state=None): + def init_widths(self, state=None): """ - Performs initialisation of width and scattering distributions + Performs initialisation of width and scattering distributions. + + When ``self.cluster`` is True and ``width_method`` is 2 or 3, a beam + loop is added so that each beam position receives a separate width + distribution computed with its own cluster-scattering parameters. + The beam-solid-angle-weighted average of those distributions is stored + in ``self.wplist`` for use in efficiency calculations. Args: state (parameters.state object., optional): if set, assume @@ -183,7 +246,7 @@ def init_widths(self,state=None): self.thresh = self.state.width.Wthresh self.wlogmean = self.state.width.Wlogmean self.wlogsigma = self.state.width.Wlogsigma - if self.meta['WBIAS'] == 'Quadrature' or self.meta['WBIAS'] == 'Sammons' or self.meta['WBIAS']=="StdDev": + if self.meta['WBIAS'] == 'Quadrature' or self.meta['WBIAS'] == 'Sammons' or self.meta['WBIAS'] == "StdDev": # these allow width-dependent sensitivities self.width_method = self.state.width.Wmethod else: @@ -195,104 +258,116 @@ def init_widths(self,state=None): if self.width_method == 3 and self.zvals is None: raise ValueError("Width method 3 requires z-values to be set") - self.NInternalBins=self.state.width.WNInternalBins + self.NInternalBins = self.state.width.WNInternalBins # records scattering information, scaling # according to frequency - self.slogmean=self.state.scat.Slogmean \ - + self.state.scat.Sfpower*np.log10( - self.meta['FBAR']/self.state.scat.Sfnorm + self.slogmean = self.state.scat.Slogmean \ + + self.state.scat.Sfpower * np.log10( + self.meta['FBAR'] / self.state.scat.Sfnorm ) - self.slogsigma=self.state.scat.Slogsigma - self.maxsigma=self.state.scat.Smaxsigma - self.scatdist=self.state.scat.ScatDist - self.backproject=self.state.scat.Sbackproject + self.slogsigma = self.state.scat.Slogsigma + self.maxsigma = self.state.scat.Smaxsigma + self.scatdist = self.state.scat.ScatDist + self.backproject = self.state.scat.Sbackproject # sets internal functions WF = self.state.width.WidthFunction - if WF ==0: + if WF == 0: self.WidthFunction = constant elif WF == 1: self.WidthFunction = lognormal elif WF == 2: self.WidthFunction = halflognormal else: - raise ValueError("state parameter scat.WidthFunction ",WF," not implemented, use 0-2 only") + raise ValueError("state parameter scat.WidthFunction ", WF, " not implemented, use 0-2 only") SF = self.state.scat.ScatFunction - if SF ==0: + if SF == 0: self.ScatFunction = constant elif SF == 1: self.ScatFunction = lognormal elif SF == 2: self.ScatFunction = halflognormal else: - raise ValueError("state parameter scat.ScatFunction ",SF," not implemented, use 0-2 only") + raise ValueError("state parameter scat.ScatFunction ", SF, " not implemented, use 0-2 only") # sets n width bins equal to zero for this survey if self.width_method == 0 or self.width_method == 4: self.NWbins = 1 ###### calculate width bins. We fix these here ###### - # Unless w and tau are explicitly being fit, it is not actually necessary - # to have constant bin values over z and DM. But best to do so! - # Here, wbins are the bin edges, w list the midpoint values used for calculations - # Nbins describes the number of bins, so Nedges is Nbins+1 - wbins = np.zeros([self.NWbins+1]) + wbins = np.zeros([self.NWbins + 1]) if self.NWbins > 1: - wbins = np.logspace(np.log10(self.WMin),np.log10(self.WMax),self.NWbins+1) - dlogw = np.log10(wbins[2]/wbins[1]) - #wbins[0] = wbins[1]-dlogw # no longer tint: 1.e-10 # set to a tiny value, to ensure we capture all small widths - # offsets the mean by half the log-spacing for each - wlist = np.logspace(np.log10(self.WMin)+dlogw/2.,np.log10(self.WMax)-dlogw/2.,self.NWbins) + wbins = np.logspace(np.log10(self.WMin), np.log10(self.WMax), self.NWbins + 1) + dlogw = np.log10(wbins[2] / wbins[1]) + wlist = np.logspace(np.log10(self.WMin) + dlogw / 2., np.log10(self.WMax) - dlogw / 2., self.NWbins) wbins[0] -= 3 # ensures we capture low values! else: wbins[0] = self.WMin wbins[1] = self.WMax - dlogw = np.log10(wbins[1]/wbins[0]) - wlist = np.array([(self.WMax*self.WMin)**0.5]) + dlogw = np.log10(wbins[1] / wbins[0]) + wlist = np.array([(self.WMax * self.WMin) ** 0.5]) self.wbins = wbins self.wlist = wlist self.dlogw = dlogw - ####### generates internal width values of numerical calculation purposes ##### - #minval = np.min([self.wlogmean - self.maxsigma*self.wlogsigma, - # self.slogmean - self.maxsigma*self.slogsigma, - # np.log10(self.WMin)]) - minval = np.log10(self.WMin)-3 + ####### generates internal width values for numerical calculation purposes ##### + minval = np.log10(self.WMin) - 3 maxval = np.log10(self.WMax) - # I haven't decided yet where to put the value of 1000 internal bins - # in terms of parameters. - self.internal_logwvals = np.linspace(minval,maxval,self.NInternalBins) + self.internal_logwvals = np.linspace(minval, maxval, self.NInternalBins) + NBINS = self.meta['NBINS'] + # initialise probability bins if self.width_method == 3: # evaluate efficiencies at each redshift - self.efficiencies = np.zeros([self.NWbins,self.NZ,self.NDM]) - self.wplist = np.zeros([self.NWbins,self.NZ]) - self.DMlist = np.zeros([self.NZ,self.NDM]) - self.mean_efficiencies = np.zeros([self.NZ,self.NDM]) - - if self.backproject: - self.pws = np.zeros([self.NZ,self.internal_logwvals.size,self.NWbins]) #[iz,:,:] = pw - self.ptaus = np.zeros([self.NZ,self.internal_logwvals.size,self.NWbins]) #[iz,:,:] = ptau - - # we have a z-dependent scattering and width model - for iz,z in enumerate(self.zvals): - self.make_widths(iz) - #_ = self.get_efficiency_from_wlist(self.wlist,self.wplist[:,iz], - # model=self.meta['WBIAS'], edir=self.edir, iz=iz) + self.efficiencies = np.zeros([self.NWbins, self.NZ, self.NDM]) + self.DMlist = np.zeros([self.NZ, self.NDM]) + self.mean_efficiencies = np.zeros([self.NZ, self.NDM]) + + if self.cluster: + # Per-(z, beam) width distributions; shape [NWbins, NZ, NBINS] + self._wplist_cluster = np.zeros([self.NWbins, self.NZ, NBINS]) + if self.backproject: + raise NotImplementedError( + "backproject=True is not supported with cluster=True") + for iz, z in enumerate(self.zvals): + for ib in range(NBINS): + self.make_widths(iz, ib) + # Beam-solid-angle-weighted average -> [NWbins, NZ] + self.wplist = np.average( + self._wplist_cluster, weights=self.beam_o, axis=2) + else: + self.wplist = np.zeros([self.NWbins, self.NZ]) + if self.backproject: + self.pws = np.zeros([self.NZ, self.internal_logwvals.size, self.NWbins]) + self.ptaus = np.zeros([self.NZ, self.internal_logwvals.size, self.NWbins]) + for iz, z in enumerate(self.zvals): + self.make_widths(iz) else: - self.wplist = np.zeros([self.NWbins]) - self.make_widths() - if self.backproject: - self.pws = np.zeros([self.internal_logwvals.size,self.NWbins]) #[iz,:,:] = pw - self.ptaus = np.zeros([self.internal_logwvals.size,self.NWbins]) #[iz,:,:] = ptau - #_ = self.get_efficiency_from_wlist(self.wlist,self.wplist, - # model=self.meta['WBIAS'], edir=self.edir, iz=None) + if self.cluster and self.width_method in [2]: + # Per-beam width distributions for non-z-dependent case; + # shape [NWbins, NBINS] + self._wplist_cluster = np.zeros([self.NWbins, NBINS]) + if self.backproject: + raise NotImplementedError( + "backproject=True is not supported with cluster=True") + # For method 2 without z-dependence we use iz=0 as a proxy + for ib in range(NBINS): + self.make_widths(iz=None, ib=ib) + # Beam-weighted average -> [NWbins] + self.wplist = np.average( + self._wplist_cluster, weights=self.beam_o, axis=1) + else: + self.wplist = np.zeros([self.NWbins]) + self.make_widths() + if self.backproject: + self.pws = np.zeros([self.internal_logwvals.size, self.NWbins]) + self.ptaus = np.zeros([self.internal_logwvals.size, self.NWbins]) self.do_efficiencies() - def process_dm_mask(self,edir=None): + def process_dm_mask(self, edir=None): """ Processes a DM mask into units of local DM """ @@ -303,7 +378,7 @@ def process_dm_mask(self,edir=None): filename = os.path.expanduser(os.path.join(edir, self.meta['DMMASK'])) if not os.path.exists(filename): - raise ValueError("Could not find DM mask " + filename + " in directory "+edir) + raise ValueError("Could not find DM mask " + filename + " in directory " + edir) # Should contain DM in the first row and efficiencies in the second row sensitivity_array = np.load(filename) @@ -311,7 +386,7 @@ def process_dm_mask(self,edir=None): effective_vals = self.dmvals + self.state.MW.DMhalo + np.median(self.DMGs) else: effective_vals = self.dmvals + self.state.MW.DMhalo + 30 - dm_mask = np.interp(effective_vals, sensitivity_array[0,:], sensitivity_array[1,:], left=1.,right=0) + dm_mask = np.interp(effective_vals, sensitivity_array[0, :], sensitivity_array[1, :], left=1., right=0) self.dm_mask = dm_mask def do_efficiencies(self): @@ -327,13 +402,14 @@ def do_efficiencies(self): self.dm_mask = None if self.width_method == 3: - for iz,z in enumerate(self.zvals): - _ = self.get_efficiency_from_wlist(self.wlist,self.wplist[:,iz], + for iz, z in enumerate(self.zvals): + _ = self.get_efficiency_from_wlist(self.wlist, self.wplist[:, iz], model=self.meta['WBIAS'], edir=self.edir, iz=iz) else: - _ = self.get_efficiency_from_wlist(self.wlist,self.wplist, - model=self.meta['WBIAS'], edir=self.edir, iz=None) - def make_widths(self,iz=None): + _ = self.get_efficiency_from_wlist(self.wlist, self.wplist, + model=self.meta['WBIAS'], edir=self.edir, iz=None) + + def make_widths(self, iz=None, ib=None): """ This routine calculates width distributions of FRBs, which is required to estimate detection efficiency, and (potentially) @@ -352,29 +428,33 @@ def make_widths(self,iz=None): this with a lognormal scattering distribution. NOTE: if wmethod=2, iz must be called with iz=None. 3: As above, but allows for redshift dependence of both intrinsic - width and scattering. Distributions both scale width redshift: + width and scattering. Distributions both scale with redshift: width with 1+z, scattering with (1+z)^-3 4: Use a specific width of a particular FRB. Used for detailed calcs for individual FRB-by-FRB analyses. + + Args: + iz (int, optional): redshift index. Required for width_method=3. + ib (int, optional): beam index. When not None (cluster mode), the + cluster scattering distribution for beam bin ib is passed to + ``quadrature_convolution`` as a third term, and the result is + stored in ``self._wplist_cluster``. """ if self.width_method == 0: # do not take a distribution, just use 1ms for everything - # this is done for tests, for complex surveys such as CHIME, - # or for estimating the properties of a single FRB self.wlist[0] = 1. self.wplist[0] = 1. elif self.width_method == 1: # take intrinsic width function only for i in np.arange(self.NWbins): - norm=(2.*np.pi)**-0.5/self.wlogsigma - args=(self.wlogmean,self.wlogsigma,norm) - weight,err=quad(self.WidthFunction, - np.log10(self.wbins[i]),np.log10(self.wbins[i+1]),args=args) - self.wplist[i]=weight + norm = (2. * np.pi) ** -0.5 / self.wlogsigma + args = (self.wlogmean, self.wlogsigma, norm) + weight, err = quad(self.WidthFunction, + np.log10(self.wbins[i]), np.log10(self.wbins[i + 1]), args=args) + self.wplist[i] = weight elif self.width_method == 2 or self.width_method == 3: # include scattering distribution. 3 means include z-dependence - #gets cumulative hist and bin edges if iz is not None: if self.width_method == 2: @@ -383,51 +463,69 @@ def make_widths(self,iz=None): exit() z = self.zvals[iz] else: - z=0. + z = 0. # performs z-scaling (does nothing when z=0.) - wlogmean = self.wlogmean + np.log10(1+z) + wlogmean = self.wlogmean + np.log10(1 + z) wlogsigma = self.wlogsigma - slogmean = self.slogmean - 3.*np.log10(1+z) + slogmean = self.slogmean - 3. * np.log10(1 + z) slogsigma = self.slogsigma - # we need these for normalisation purposes. Otherwise, the - # scattering distribution gets diluted with a constant max value - # these are only used for normalisation purposes of otherwise - # unnormalisable functions. - smax = np.log10(self.WMax) - 3.*np.log10(1+z) - wmax = np.log10(self.WMax) + np.log10(1+z) + # we need these for normalisation purposes. + smax = np.log10(self.WMax) - 3. * np.log10(1 + z) + wmax = np.log10(self.WMax) + np.log10(1 + z) - WidthArgs = (wlogmean,wlogsigma,wmax) - ScatArgs = (slogmean,slogsigma,smax) + WidthArgs = (wlogmean, wlogsigma, wmax) + ScatArgs = (slogmean, slogsigma, smax) - if self.backproject: - dist,pw,ptau = quadrature_convolution(self.WidthFunction, WidthArgs, self.ScatFunction, ScatArgs, - self.internal_logwvals, self.wbins,backproject=self.backproject) + # --- Cluster scattering: select per-(iz,ib) parameters --- + if self.cluster and ib is not None: + xps = self.xProbScat + ps = self.probScat[:, iz, ib] if iz is not None else self.probScat[:, 0, ib] + fu = self.fractionUnscattered[iz, ib] if iz is not None else self.fractionUnscattered[0, ib] else: - dist = quadrature_convolution(self.WidthFunction, WidthArgs, self.ScatFunction, ScatArgs, - self.internal_logwvals, self.wbins,backproject=self.backproject) + xps = None + ps = None + fu = 1.0 + # ---------------------------------------------------------- + if self.backproject: + dist, pw, ptau = quadrature_convolution( + self.WidthFunction, WidthArgs, self.ScatFunction, ScatArgs, + self.internal_logwvals, self.wbins, backproject=self.backproject, + xProbScat=xps, probScat=ps, fractionUnscattered=fu) + else: + dist = quadrature_convolution( + self.WidthFunction, WidthArgs, self.ScatFunction, ScatArgs, + self.internal_logwvals, self.wbins, backproject=self.backproject, + xProbScat=xps, probScat=ps, fractionUnscattered=fu) - if iz is not None: - self.wplist[:,iz] = dist + # Store result in the right place + if ib is not None and self.cluster: + # cluster mode: store in per-beam array + if iz is not None: + self._wplist_cluster[:, iz, ib] = dist + else: + self._wplist_cluster[:, ib] = dist + elif iz is not None: + self.wplist[:, iz] = dist if self.backproject: - self.pws[iz,:,:] = pw - self.ptaus[iz,:,:] = ptau + self.pws[iz, :, :] = pw + self.ptaus[iz, :, :] = ptau else: self.wplist[:] = dist if self.backproject: - self.pws=pw - self.ptaus=ptau + self.pws = pw + self.ptaus = ptau elif self.width_method == 4: # use specific width of FRB. This requires there to be only a single FRB in the survey - if s.meta['NFRB'] != 1: - raise ValueError("If width method in make_widths is 4 only one FRB should be specified in the survey but ", str(s.meta['NFRB']), " FRBs were specified") + if self.meta['NFRB'] != 1: + raise ValueError("If width method in make_widths is 4 only one FRB should be specified in the survey but ", str(self.meta['NFRB']), " FRBs were specified") else: self.wplist[0] = 1. - self.wlist[0] = s.frbs['WIDTH'][0] + self.wlist[0] = self.frbs['WIDTH'][0] else: - raise ValueError("Width method in make_widths must be 0-4, not ",width_method) + raise ValueError("Width method in make_widths must be 0-4, not ", self.width_method) def init_repeaters(self): @@ -442,7 +540,7 @@ def init_repeaters(self): if not "NREP" in self.frbs: print("Warning, no information of repetition provided") print("Assuming all FRBs are once-off bursts") - self.frbs["NREP"] = np.full([self.NFRB],1,dtype='int') + self.frbs["NREP"] = np.full([self.NFRB], 1, dtype='int') # Set repeater/singles list self.replist = np.where(self.frbs["NREP"] > 1)[0] @@ -496,7 +594,6 @@ def init_repeaters(self): # Case of NORM_FRB not set else: self.NORM_FRB = self.NORM_REPS + self.NORM_SINGLES - # print("NORM_FRB set to NORM_REPS + NORM_SINGLES = " + str(self.NORM_FRB)) elif self.NORM_REPS is None: # Case invalid if self.NORM_SINGLES is None or self.NORM_SINGLES > self.NORM_FRB: @@ -504,7 +601,6 @@ def init_repeaters(self): # Case NORM_REPS not set else: self.NORM_REPS = self.NORM_FRB - self.NORM_SINGLES - # print("NORM_REPS set to NORM_FRB - NORM_SINGLES = " + str(self.NORM_REPS)) elif self.NORM_SINGLES is None: # Do not need to consider NORM_REPS == None as that is done in the previous elif # Case invalid @@ -513,7 +609,6 @@ def init_repeaters(self): # Case NORM_SINGLES not set else: self.NORM_SINGLES = self.NORM_FRB - self.NORM_REPS - # print("NORM_SINGLES set to NORM_FRB - NORM_REPS = " + str(self.NORM_SINGLES)) # Case all 3 are set else: # Common sense check @@ -524,7 +619,7 @@ def init_repeaters(self): self.init_zs_reps() - def get_internal_coeffs(self,wlist): + def get_internal_coeffs(self, wlist): """ Returns indices and coefficients for linear interpolation between intrinsic width values @@ -540,20 +635,20 @@ def get_internal_coeffs(self,wlist): # convert to log-widths - the bins are in log10 space logwlist = np.log10(wlist) - dinternal = self.internal_logwvals[1]-self.internal_logwvals[0] - kws=(logwlist-self.internal_logwvals[0])/dinternal + dinternal = self.internal_logwvals[1] - self.internal_logwvals[0] + kws = (logwlist - self.internal_logwvals[0]) / dinternal Bin0 = np.where(kws < 0.)[0] kws[Bin0] = 0. - iws1=kws.astype('int') - iws2=iws1+1 - dkws2=kws-iws1 # applies to izs2 + iws1 = kws.astype('int') + iws2 = iws1 + 1 + dkws2 = kws - iws1 # applies to izs2 dkws1 = 1. - dkws2 # checks for values which are too large toobigw = np.where(logwlist > self.internal_logwvals[-1])[0] if len(toobigw) > 0: - raise ValueError("Width value ",wlist[toobigw], - " too large for max internal log value of ",10**self.internal_logwvals[-1]) + raise ValueError("Width value ", wlist[toobigw], + " too large for max internal log value of ", 10 ** self.internal_logwvals[-1]) return iws1, iws2, dkws1, dkws2 @@ -567,33 +662,32 @@ def init_frb_bvals(self): """ # contains beam-dependent weights for each FRB - frb_bweights = np.zeros([self.NFRB,self.meta["NBINS"]]) + frb_bweights = np.zeros([self.NFRB, self.meta["NBINS"]]) lbs = np.log(self.beam_b) - for i,B in enumerate(self.frbs["B"]): + for i, B in enumerate(self.frbs["B"]): if B == -1: # code for "ignore it" - # still have to decide what to do here. Likely give equal weights? - frb_bweights[i,:] = 1./self.meta["NBINS"] + frb_bweights[i, :] = 1. / self.meta["NBINS"] elif B > self.beam_b[-1]: # greater value than max, just use max - frb_bweights[i,-1] = 1. + frb_bweights[i, -1] = 1. elif B < self.beam_b[0]: - frb_bweights[i,0] = 1. + frb_bweights[i, 0] = 1. else: # at least one value of beam_b will be greater and one lesser than B iB2 = np.where(self.beam_b > B)[0][0] iB1 = iB2 - 1 # do log-scaling lB = np.log(B) - kB2 = (lB- lbs[iB1])/(lbs[iB2]-lbs[iB1]) - kB1 = 1.-kB2 - frb_bweights[i,iB1] = kB1 - frb_bweights[i,iB2] = kB2 + kB2 = (lB - lbs[iB1]) / (lbs[iB2] - lbs[iB1]) + kB1 = 1. - kB2 + frb_bweights[i, iB1] = kB1 + frb_bweights[i, iB2] = kB2 self.frb_bweights = frb_bweights # speedups when iterating through 1D and 2D likelihoods - self.frb_zbweights = frb_bweights[self.zlist,:] - self.frb_nozbweights = frb_bweights[self.nozlist,:] + self.frb_zbweights = frb_bweights[self.zlist, :] + self.frb_nozbweights = frb_bweights[self.nozlist, :] def init_frb_wvals(self): @@ -603,28 +697,26 @@ def init_frb_wvals(self): This is for a slightly different purpose than the init widths routine """ - - #wlist = survey.WIDTHs # measured total widths nw = self.wlist.size - frb_wweights = np.zeros([self.NFRB,nw]) + frb_wweights = np.zeros([self.NFRB, nw]) OKw = np.where(self.WIDTHs > 0.) notOKw = np.where(self.WIDTHs <= 0.) # equal weights for all FRBs with no measured width - frb_wweights[notOKw,:] = 1./nw + frb_wweights[notOKw, :] = 1. / nw - # iterature through the list - iws1,iws2,dkws1,dkws2 = self.get_w_coeffs(self.WIDTHs[OKw]) - frb_wweights[OKw,iws1] = dkws1 - frb_wweights[OKw,iws2] = dkws2 + # iterate through the list + iws1, iws2, dkws1, dkws2 = self.get_w_coeffs(self.WIDTHs[OKw]) + frb_wweights[OKw, iws1] = dkws1 + frb_wweights[OKw, iws2] = dkws2 self.frb_wweights = frb_wweights # speedups when iterating through 1D and 2D likelihoods - self.frb_zwweights = self.frb_wweights[self.zlist,:] - self.frb_nozwweights = self.frb_wweights[self.nozlist,:] + self.frb_zwweights = self.frb_wweights[self.zlist, :] + self.frb_nozwweights = self.frb_wweights[self.nozlist, :] - def get_w_coeffs(self,wlist): + def get_w_coeffs(self, wlist): """ Returns indices and coefficients for linear interpolation between width values Bin edges run from [small~1e-10, self.WMin, self.WMin + self.dlowg, ..., self.Wmax] @@ -644,54 +736,50 @@ def get_w_coeffs(self,wlist): if self.NWbins == 1: # only when there is a single width bin nfrb = logwlist.size - iws1 = np.full([nfrb],0,dtype='int') + iws1 = np.full([nfrb], 0, dtype='int') iws2 = iws1 - dkws1 = np.full([nfrb],1.,dtype='float') + dkws1 = np.full([nfrb], 1., dtype='float') dkws2 = dkws1 - # dkws2 is identical to 1. This over-writes 1, but ensures - # that order of implementation of 1 and 2 does not matter return iws1, iws2, dkws1, dkws2 - # the below assumes that - kws=(logwlist-np.log10(self.WMin))/self.dlogw # now will assume it begins at Wmin+dlogw + 0.5 + kws = (logwlist - np.log10(self.WMin)) / self.dlogw # forces any low values to zero Bin0 = np.where(kws < 0.)[0] kws[Bin0] = 0. - iws1=kws.astype('int') - iws2=iws1+1 - dkws2=kws-iws1 # applies to izs2 + iws1 = kws.astype('int') + iws2 = iws1 + 1 + dkws2 = kws - iws1 # applies to izs2 dkws1 = 1. - dkws2 - # in case iws1 is in the largets bin, then - # iws2 will be too large + # in case iws1 is in the largest bin, then iws2 will be too large toobig1 = np.where(iws2 >= self.wlist.size)[0] - iws2[toobig1]=self.wlist.size-1 + iws2[toobig1] = self.wlist.size - 1 # checks for values which are too large toobigw = np.where(wlist > self.WMax)[0] if len(toobigw) > 0: - raise ValueError("Width value ",wlist[toobigw], - " too large for Wmax of ",self.WMax) + raise ValueError("Width value ", wlist[toobigw], + " too large for Wmax of ", self.WMax) return iws1, iws2, dkws1, dkws2 def randomise_DMG(self, uDMG=0.5): """ Change the DMG_ISM values to a random value within uDMG Gaussian uncertainty """ - new_DMGs = np.random.normal(self.DMGs, uDMG*self.DMGs) + new_DMGs = np.random.normal(self.DMGs, uDMG * self.DMGs) neg = np.where(new_DMGs < 0)[0] while len(neg) != 0: - new_DMGs[neg] = np.random.normal(self.DMGs[neg], uDMG*self.DMGs[neg]) + new_DMGs[neg] = np.random.normal(self.DMGs[neg], uDMG * self.DMGs[neg]) neg = np.where(new_DMGs < 0)[0] self.DMGs = new_DMGs - def init_DMEG(self,DMhalo,halo_method=0): + def init_DMEG(self, DMhalo, halo_method=0): """ Calculates extragalactic DMs assuming halo DM """ - self.DMhalo=DMhalo + self.DMhalo = DMhalo self.process_dmhalo(halo_method) - self.DMEGs=self.DMs-self.DMGs - self.DMhalos + self.DMEGs = self.DMs - self.DMGs - self.DMhalos def process_dmhalo(self, halo_method): """ @@ -704,7 +792,6 @@ def process_dmhalo(self, halo_method): # Constant halo if halo_method == 0: self.DMhalos = np.ones(self.DMs.shape) * self.DMhalo - # self.DMGals = self.DMhalos + self.DMGs # Yamasaki and Totani 2020 elif halo_method == 1: @@ -766,8 +853,6 @@ def process_dmhalo(self, halo_method): self.DMhalo = np.median(self.DMhalos) print(self.DMhalos) - # self.DMGal = np.median(self.DMGals) - def init_halo_coeffs(self): """ Initialise coefficients for Yamasaki and Totani 2020 implementation of @@ -811,15 +896,15 @@ def init_zs(self): self.ignored_Zlist = [] # Pandas resolves None to Nan - if len(self.frbs["Z"])>0: + if len(self.frbs["Z"]) > 0: - self.Zs=np.array(self.frbs["Z"].values).astype('float') + self.Zs = np.array(self.frbs["Z"].values).astype('float') # checks for any redhsifts identically equal to zero #exactly zero can be bad... only happens in MC generation # 0.001 is chosen as smallest redshift in original fit zeroz = np.where(self.Zs == 0.)[0] - if len(zeroz) >0: - self.Zs[zeroz]=0.001 + if len(zeroz) > 0: + self.Zs[zeroz] = 0.001 # checks to see if there are any FRBs which are localised self.zlist = np.where(self.Zs > 0.)[0] @@ -827,19 +912,19 @@ def init_zs(self): if len(self.zlist) < self.NFRB: self.nozlist = np.where(self.Zs < 0.)[0] if len(self.nozlist) == len(self.Zs): - self.nD=1 # they all had -1 as their redshift! - self.zlist=None + self.nD = 1 # they all had -1 as their redshift! + self.zlist = None else: - self.nD=3 # code for both + self.nD = 3 # code for both else: self.nozlist = None - self.nD=2 + self.nD = 2 else: - self.nD=1 - self.Zs=None - self.nozlist=np.arange(self.NFRB) - self.zlist=None + self.nD = 1 + self.Zs = None + self.nozlist = np.arange(self.NFRB) + self.zlist = None def init_zs_reps(self): """ @@ -909,22 +994,22 @@ def init_zs_reps(self): self.nDs = 3 # initialise rep-dependent beam and width weights - self.frb_zbweights_singles = self.frb_bweights[self.zsingles,:] - self.frb_zbweights_reps = self.frb_bweights[self.zreps,:] - self.frb_nozbweights_singles = self.frb_bweights[self.nozsingles,:] - self.frb_nozbweights_reps = self.frb_bweights[self.nozreps,:] - - self.frb_zwweights_singles = self.frb_wweights[self.zsingles,:] - self.frb_zwweights_reps = self.frb_wweights[self.zreps,:] - self.frb_nozwweights_singles = self.frb_wweights[self.nozsingles,:] - self.frb_nozwweights_reps = self.frb_wweights[self.nozreps,:] - - def process_survey_file(self,filename:str, - NFRB:int=None, - iFRB:int=0, + self.frb_zbweights_singles = self.frb_bweights[self.zsingles, :] + self.frb_zbweights_reps = self.frb_bweights[self.zreps, :] + self.frb_nozbweights_singles = self.frb_bweights[self.nozsingles, :] + self.frb_nozbweights_reps = self.frb_bweights[self.nozreps, :] + + self.frb_zwweights_singles = self.frb_wweights[self.zsingles, :] + self.frb_zwweights_reps = self.frb_wweights[self.zreps, :] + self.frb_nozwweights_singles = self.frb_wweights[self.nozsingles, :] + self.frb_nozwweights_reps = self.frb_wweights[self.nozreps, :] + + def process_survey_file(self, filename: str, + NFRB: int = None, + iFRB: int = 0, min_lat=None, dmg_cut=None, - survey_dict = None): + survey_dict=None): """ Loads a survey file, then creates dictionaries of the loaded variables @@ -948,20 +1033,20 @@ def process_survey_file(self,filename:str, frb_tbl.meta['survey_data']) - # Meta -- for convenience for now;  best to migrate away from this + # Meta -- for convenience for now; best to migrate away from this default_telescope = survey_data.Telescope() for key in self.survey_data.params: DC = self.survey_data.params[key] if DC == "telescope": - value = getattr(self.survey_data[DC],key) - if value == getattr(default_telescope,key): + value = getattr(self.survey_data[DC], key) + if value == getattr(default_telescope, key): # using default value - check if the FRBs have this if key in frb_tbl.columns: value = np.mean(frb_tbl[key]) self.meta[key] = value else: - self.meta[key] = getattr(self.survey_data[DC],key) + self.meta[key] = getattr(self.survey_data[DC], key) # Get default values from default frb data default_frb = survey_data.FRB() @@ -994,12 +1079,34 @@ def process_survey_file(self,filename:str, # Cut down? # NFRB if self.NFRB is not None: - self.NFRB=min(len(self.frbs), NFRB) - if self.NFRB < NFRB+iFRB: + self.NFRB = min(len(self.frbs), NFRB) + if self.NFRB < NFRB + iFRB: raise ValueError("Cannot return sufficient FRBs, did you mean NFRB=None?") - # Not sure the following linematters given the Error above - themax = max(NFRB+iFRB,self.NFRB) - self.frbs=self.frbs[iFRB:themax] + # Not sure the following line matters given the Error above + themax = max(NFRB + iFRB, self.NFRB) + self.frbs = self.frbs[iFRB:themax] + + # fills in missing coordinates if possible + # also converts RA and Dec strings to floats + self.fix_coordinates(verbose=False) + + # Min latitude + if min_lat is not None and min_lat > 0.0: + tot = len(self.frbs) + mask = [(Gb is None) or (np.abs(Gb) > min_lat) for Gb in self.frbs['Gb'].values] + self.frbs = self.frbs[mask] + included = len(self.frbs) + + print("Using minimum galactic latitude of " + str(min_lat) + ". Excluding " + str(tot - included) + " FRBs") + # Max DM + if dmg_cut is not None: + tot = len(self.frbs) + self.frbs = self.frbs[np.abs(self.frbs['DMG'].values) < dmg_cut] + included = len(self.frbs) + print("Using maximum DMG of " + str(dmg_cut) + ". Excluding " + str(tot - included) + " FRBs") + + # Get new number of FRBs + self.NFRB = len(self.frbs) # fills in missing coordinates if possible # also converts RA and Dec strings to floats @@ -1050,26 +1157,26 @@ def process_survey_file(self,filename:str, if survey_dict is not None: for key in survey_dict: self.meta[key] = survey_dict[key] - print(self.meta[key], "overidden by ",survey_dict[key]) + print(self.meta[key], "overidden by ", survey_dict[key]) ### processes galactic contributions self.process_dmg() ### get pointers to correct results ,for better access - self.DMs=self.frbs['DM'].values - self.DMGs=self.frbs['DMG'].values - self.SNRs=self.frbs['SNR'].values - self.WIDTHs=self.frbs['WIDTH'].values - self.TRESs=self.frbs['TRES'].values - self.FRESs=self.frbs['FRES'].values - self.FBARs=self.frbs['FBAR'].values - self.BWs=self.frbs['BW'].values - self.THRESHs=self.frbs['THRESH'].values - self.SNRTHRESHs=self.frbs['SNRTHRESH'].values - self.Ss=self.SNRs/self.SNRTHRESHs - self.TOBS=self.meta['TOBS'] - self.NORM_FRB=self.meta['NORM_FRB'] - self.Gbs=self.frbs['Gb'].values - self.Gls=self.frbs['Gl'].values + self.DMs = self.frbs['DM'].values + self.DMGs = self.frbs['DMG'].values + self.SNRs = self.frbs['SNR'].values + self.WIDTHs = self.frbs['WIDTH'].values + self.TRESs = self.frbs['TRES'].values + self.FRESs = self.frbs['FRES'].values + self.FBARs = self.frbs['FBAR'].values + self.BWs = self.frbs['BW'].values + self.THRESHs = self.frbs['THRESH'].values + self.SNRTHRESHs = self.frbs['SNRTHRESH'].values + self.Ss = self.SNRs / self.SNRTHRESHs + self.TOBS = self.meta['TOBS'] + self.NORM_FRB = self.meta['NORM_FRB'] + self.Gbs = self.frbs['Gb'].values + self.Gls = self.frbs['Gl'].values # calculates intrinsic widths # Uses the model of James et al 2025 @@ -1077,74 +1184,74 @@ def process_survey_file(self,filename:str, # if scattering dominates total width, expect tau = 0.816 w tscale = 1.225 # scale scattering time to total width at +- 1 sigma - TEMP = self.frbs['WIDTH'].values**2 - (tscale*self.frbs['TAU'].values)**2 + TEMP = self.frbs['WIDTH'].values ** 2 - (tscale * self.frbs['TAU'].values) ** 2 self.OKTAU = np.where(self.frbs['TAU'].values != -1.)[0] # code for non-existent toolow = np.where(TEMP <= 0.) - TEMP[toolow] = 0.01*self.frbs['TAU'].values[toolow]**2 # 10% of scattering width - iwidths = TEMP**0.5 # scale to SNR max width assuming Gaussian shape + TEMP[toolow] = 0.01 * self.frbs['TAU'].values[toolow] ** 2 # 10% of scattering width + iwidths = TEMP ** 0.5 # scale to SNR max width assuming Gaussian shape self.IWIDTHs = iwidths self.TAUs = self.frbs['TAU'].values # sets the 'beam' values to unity by default - self.beam_b=np.array([1]) - self.beam_o=np.array([1]) - self.NBEAMS=1 + self.beam_b = np.array([1]) + self.beam_o = np.array([1]) + self.NBEAMS = 1 - # checks for incorrectSNR values + # checks for incorrect SNR values toolow = np.where(self.Ss < 1.)[0] if len(toolow) > 0: - raise ValueError("FRBs ",toolow," have SNR < SNRTHRESH!!! Please correct this. Exiting...") + raise ValueError("FRBs ", toolow, " have SNR < SNRTHRESH!!! Please correct this. Exiting...") - print("FRB survey sucessfully initialised with ",self.NFRB," FRBs starting from", self.iFRB) + print("FRB survey sucessfully initialised with ", self.NFRB, " FRBs starting from", self.iFRB) - def fix_coordinates(self,verbose=False): + def fix_coordinates(self, verbose=False): """ Takes and FRB, and fills out missing coordinate values Note that now, RA, DEC, Gl, and Gb will be present But their default values are None """ # converts to float if in string. Will do nothing if None - if isinstance(self.frbs['RA'][0],str): + if isinstance(self.frbs['RA'][0], str): RAs = np.zeros([len(self.frbs['RA'])]) - for i,RA in enumerate(self.frbs['RA']): - RAs[i] = misc_functions.coord_string_to_deg(self.frbs['RA'][i],hr=True) + for i, RA in enumerate(self.frbs['RA']): + RAs[i] = misc_functions.coord_string_to_deg(self.frbs['RA'][i], hr=True) self.frbs['RA'] = RAs - if isinstance(self.frbs['DEC'][0],str): + if isinstance(self.frbs['DEC'][0], str): DECs = np.zeros([len(self.frbs['DEC'])]) - for i,DEC in enumerate(self.frbs['DEC']): - DECs[i] = misc_functions.coord_string_to_deg(self.frbs['DEC'][i],hr=False) + for i, DEC in enumerate(self.frbs['DEC']): + DECs[i] = misc_functions.coord_string_to_deg(self.frbs['DEC'][i], hr=False) self.frbs['DEC'] = DECs - if isinstance(self.frbs['Gb'][0],str): + if isinstance(self.frbs['Gb'][0], str): Gbs = np.zeros([len(self.frbs['Gb'])]) - for i,Gb in enumerate(self.frbs['Gb']): - Gbs[i] = misc_functions.coord_string_to_deg(self.frbs['Gb'][i],hr=True) + for i, Gb in enumerate(self.frbs['Gb']): + Gbs[i] = misc_functions.coord_string_to_deg(self.frbs['Gb'][i], hr=True) self.frbs['Gb'] = Gbs - if isinstance(self.frbs['Gl'][0],str): + if isinstance(self.frbs['Gl'][0], str): Gls = np.zeros([len(self.frbs['Gl'])]) - for i,Gl in enumerate(self.frbs['Gl']): - Gls[i] = misc_functions.coord_string_to_deg(self.frbs['Gl'][i],hr=False) + for i, Gl in enumerate(self.frbs['Gl']): + Gls[i] = misc_functions.coord_string_to_deg(self.frbs['Gl'][i], hr=False) self.frbs['Gl'] = Gls - for i,gl in enumerate(self.frbs['Gl']): + for i, gl in enumerate(self.frbs['Gl']): if gl is None or self.frbs['Gb'][i] is None: # test RA if self.frbs['RA'][i] is None or self.frbs['DEC'][i] is None: if verbose: - print("WARNING: no coordinates calculable for FRB ",i) + print("WARNING: no coordinates calculable for FRB ", i) else: - Gb,Gl = misc_functions.j2000_to_galactic(self.frbs['RA'][i], self.frbs['DEC'][i]) - self.frbs[i,'Gb'] = Gb - self.frbs[i,'Gl'] = Gl + Gb, Gl = misc_functions.j2000_to_galactic(self.frbs['RA'][i], self.frbs['DEC'][i]) + self.frbs[i, 'Gb'] = Gb + self.frbs[i, 'Gl'] = Gl elif self.frbs['RA'][i] is None or self.frbs['DEC'][i] is None: - RA,Dec = misc_functions.galactic_to_j2000(self.frbs['Gl'][i], self.frbs['Gb'][i]) - self.frbs[i,'RA'] = RA - self.frbs[i,'DEC'] = Dec + RA, Dec = misc_functions.galactic_to_j2000(self.frbs['Gl'][i], self.frbs['Gb'][i]) + self.frbs[i, 'RA'] = RA + self.frbs[i, 'DEC'] = Dec def process_dmg(self): """ Estimates galactic DM according to @@ -1158,47 +1265,47 @@ def process_dmg(self): it as DMG') print("Calculating DMG from NE2001. Please record this, it takes a while!") ne = density.ElectronDensity() - DMGs=np.zeros([self.NFRB]) - for i,l in enumerate(self.frbs["Gl"]): - b=self.frbs["Gb"][i] + DMGs = np.zeros([self.NFRB]) + for i, l in enumerate(self.frbs["Gl"]): + b = self.frbs["Gb"][i] ismDM = ne.DM(l, b, 100.) - print(i,l,b,ismDM) - DMGs=np.array(DMGs) - self.frbs["DMG"]=DMGs - self.DMGs=DMGs + print(i, l, b, ismDM) + DMGs = np.array(DMGs) + self.frbs["DMG"] = DMGs + self.DMGs = DMGs - def init_beam(self,plot=False, - method=1,thresh=1e-3): + def init_beam(self, plot=False, + method=1, thresh=1e-3): """ Initialises the beam """ # Gaussian beam if method == 0 - if method==0: - b,omegab=beams.gauss_beam(thresh=thresh, + if method == 0: + b, omegab = beams.gauss_beam(thresh=thresh, nbins=self.meta["NBINS"], - freq=self.meta["FBAR"],D=self.meta["DIAM"]) - self.beam_b=b - self.beam_o=omegab*self.meta["NBEAMS"] - self.orig_beam_b=self.beam_b - self.orig_beam_o=self.beam_o + freq=self.meta["FBAR"], D=self.meta["DIAM"]) + self.beam_b = b + self.beam_o = omegab * self.meta["NBEAMS"] + self.orig_beam_b = self.beam_b + self.orig_beam_o = self.beam_o elif self.meta["BEAM"] is not None: - logb,omegab=beams.load_beam(self.meta["BEAM"]) - self.orig_beam_b=10**logb - self.orig_beam_o=omegab + logb, omegab = beams.load_beam(self.meta["BEAM"]) + self.orig_beam_b = 10 ** logb + self.orig_beam_o = omegab if plot: - savename='Plots/Beams/'+self.name+'_'+self.meta["BEAM"]+'_'+str(method)+'_'+str(thresh)+'_beam.pdf' + savename = 'Plots/Beams/' + self.name + '_' + self.meta["BEAM"] + '_' + str(method) + '_' + str(thresh) + '_beam.pdf' else: - savename=None - b2,o2=beams.simplify_beam(logb,omegab,self.meta["NBINS"], - savename=savename,method=method,thresh=thresh) - # there is a chance that this method alters the expected number of bins. Reset it!~ + savename = None + b2, o2 = beams.simplify_beam(logb, omegab, self.meta["NBINS"], + savename=savename, method=method, thresh=thresh) + # there is a chance that this method alters the expected number of bins. Reset it! self.meta["NBINS"] = len(o2) - self.beam_b=b2 - self.beam_o=o2 - self.do_beam=True + self.beam_b = b2 + self.beam_o = o2 + self.do_beam = True # sets the 'beam' values to unity by default - self.NBEAMS=b2.size + self.NBEAMS = b2.size else: print("No beam found to initialise...") @@ -1206,23 +1313,19 @@ def init_beam(self,plot=False, def calc_max_dm(self): ''' Calculates the maximum searched DM. - - Calculates bandwidth using ''' - fbar=self.meta['FBAR'] - t_res=self.meta['TRES'] - nu_res=self.meta['FRES'] - max_idt=self.meta['MAX_IDT'] - max_dm=self.meta['MAX_DM'] + fbar = self.meta['FBAR'] + t_res = self.meta['TRES'] + nu_res = self.meta['FRES'] + max_idt = self.meta['MAX_IDT'] + max_dm = self.meta['MAX_DM'] if max_dm is None and max_idt is not None: - k_DM=4.149 #ms GHz^2 pc^-1 cm^3 - #f_low = fbar - (Nchan/2. - 1)*nu_res - #f_high = fbar + (Nchan/2. - 1)*nu_res - f_low = fbar - self.meta['BW']/2. # bottom of lowest band - f_high = fbar + self.meta['BW']/2. # top of highest band + k_DM = 4.149 #ms GHz^2 pc^-1 cm^3 + f_low = fbar - self.meta['BW'] / 2. # bottom of lowest band + f_high = fbar + self.meta['BW'] / 2. # top of highest band max_dt = t_res * max_idt - max_dm = max_dt / (k_DM * ((f_low/1e3)**(-2) - (f_high/1e3)**(-2))) + max_dm = max_dt / (k_DM * ((f_low / 1e3) ** (-2) - (f_high / 1e3) ** (-2))) self.max_dm = max_dm @@ -1236,7 +1339,7 @@ def calc_max_dm(self): self.max_idm = None self.max_dmeg = None - def get_efficiency_from_wlist(self,wlist,plist, + def get_efficiency_from_wlist(self, wlist, plist, model="Quadrature", addGalacticDM=True, edir=None, iz=None): @@ -1267,18 +1370,16 @@ def get_efficiency_from_wlist(self,wlist,plist, """ DMlist = self.dmvals - efficiencies=np.zeros([wlist.size,DMlist.size]) + efficiencies = np.zeros([wlist.size, DMlist.size]) if addGalacticDM: - # toAdd = self.DMhalo + self.meta['DMG'] toAdd = np.median(self.DMhalos + self.DMGs) - # toAdd = self.DMGal else: toAdd = 0. - for i,w in enumerate(wlist): - efficiencies[i,:]=calc_relative_sensitivity( - None,DMlist+toAdd,w, + for i, w in enumerate(wlist): + efficiencies[i, :] = calc_relative_sensitivity( + None, DMlist + toAdd, w, self.meta['FBAR'], self.meta['TRES'], self.meta['FRES'], @@ -1288,20 +1389,20 @@ def get_efficiency_from_wlist(self,wlist,plist, dsmear=False, edir=edir, max_iw=self.meta['MAX_IW'], - max_meth = self.meta['MAXWMETH']) + max_meth=self.meta['MAXWMETH']) # keep an internal record of this if iz is None: - self.efficiencies=efficiencies - self.DMlist=DMlist - mean_efficiencies=np.mean(efficiencies,axis=0) - self.mean_efficiencies=mean_efficiencies #be careful here!!! This may not be what we want! + self.efficiencies = efficiencies + self.DMlist = DMlist + mean_efficiencies = np.mean(efficiencies, axis=0) + self.mean_efficiencies = mean_efficiencies else: - self.efficiencies[:,iz,:]=efficiencies - self.DMlist[iz,:]=DMlist - mean_efficiencies=np.mean(efficiencies,axis=0) - self.mean_efficiencies[iz,:]=mean_efficiencies #be careful here!!! This may not be what we want! + self.efficiencies[:, iz, :] = efficiencies + self.DMlist[iz, :] = DMlist + mean_efficiencies = np.mean(efficiencies, axis=0) + self.mean_efficiencies[iz, :] = mean_efficiencies return efficiencies @@ -1317,9 +1418,9 @@ def __repr__(self): # implements something like Mawson's formula for sensitivity # t_res in ms -def calc_relative_sensitivity(DM_frb,DM,w,fbar,t_res,nu_res,Nchan=336,max_idt=None, - max_dm=None,model='Quadrature',dsmear=True,edir=None,max_iw=None, - max_meth = 0): +def calc_relative_sensitivity(DM_frb, DM, w, fbar, t_res, nu_res, Nchan=336, max_idt=None, + max_dm=None, model='Quadrature', dsmear=True, edir=None, max_iw=None, + max_meth=0): """ Calculates DM-dependent sensitivity This function adjusts sensitivity to a given burst as a function of DM. @@ -1353,7 +1454,7 @@ def calc_relative_sensitivity(DM_frb,DM,w,fbar,t_res,nu_res,Nchan=336,max_idt=No # this model returns the parameterised CHIME DM-dependent sensitivity # it is independent of width - if model=='CHIME': + if model == 'CHIME': # polynomial coefficients for fit to CHIME DM bias data (4th order poly) coeffs = np.array([ 7.79309074e-03, -2.09210057e-01, 1.93122752e+00, -7.05813760e+00, 8.93355593e+00]) @@ -1361,73 +1462,67 @@ def calc_relative_sensitivity(DM_frb,DM,w,fbar,t_res,nu_res,Nchan=336,max_idt=No coeffs /= 1.118694423940629 # fit is to natural log of DM values ldm = np.log(DM) - rate = np.polyval(coeffs,ldm) + rate = np.polyval(coeffs, ldm) # scale rate by assumed Cartesian logN-logS - sensitivity = rate**(2./3.) + sensitivity = rate ** (2. / 3.) # calculates relative sensitivity to bursts as a function of DM # Check for Quadrature and Sammons - elif model == 'Quadrature' or model == 'Sammons' or model=="StdDev": + elif model == 'Quadrature' or model == 'Sammons' or model == "StdDev": # constant of DM - k_DM=4.149 #ms GHz^2 pc^-1 cm^3 + k_DM = 4.149 #ms GHz^2 pc^-1 cm^3 # total smearing factor within a channel - dm_smearing=2*(nu_res/1.e3)*k_DM*DM/(fbar/1e3)**3 #smearing factor of FRB in the band + dm_smearing = 2 * (nu_res / 1.e3) * k_DM * DM / (fbar / 1e3) ** 3 #smearing factor of FRB in the band # this assumes that what we see are measured widths including all the smearing factors # hence we must first adjust for this prior to estimating the DM-dependence # for this we use the *true* DM at which the FRB was observed - if dsmear==True: + if dsmear == True: # width is the total width - measured_dm_smearing=2*(nu_res/1.e3)*k_DM*DM_frb/(fbar/1e3)**3 #smearing factor of FRB in the band + measured_dm_smearing = 2 * (nu_res / 1.e3) * k_DM * DM_frb / (fbar / 1e3) ** 3 if model == "StdDev": - uw = w**2 - dm_smearing**2/3. - t_res**2/3. + uw = w ** 2 - dm_smearing ** 2 / 3. - t_res ** 2 / 3. else: - uw = w**2-measured_dm_smearing**2-t_res**2 # uses the quadrature model to calculate intrinsic width uw + uw = w ** 2 - measured_dm_smearing ** 2 - t_res ** 2 if uw < 0: - uw=1e-2 # replace this with some fraction of minimum width? + uw = 1e-2 else: - uw=uw**0.5 + uw = uw ** 0.5 else: # w represents the intrinsic width - uw=w + uw = w if model == "StdDev": - # 2* to be +- one standard deviation. But we're now fixing Tres to be unity - totalw = (uw**2 + dm_smearing**2/3. + t_res**2)**0.5 - nosmearw = (uw**2 + t_res**2)**0.5 + # 2* to be +- one standard deviation. + totalw = (uw ** 2 + dm_smearing ** 2 / 3. + t_res ** 2) ** 0.5 + nosmearw = (uw ** 2 + t_res ** 2) ** 0.5 else: - totalw = (uw**2 + dm_smearing**2 + t_res**2)**0.5 - nosmearw = (uw**2 + t_res**2)**0.5 + totalw = (uw ** 2 + dm_smearing ** 2 + t_res ** 2) ** 0.5 + nosmearw = (uw ** 2 + t_res ** 2) ** 0.5 # calculates relative sensitivity to bursts as a function of DM - if model=='Quadrature' or model=="StdDev": - sensitivity=totalw**-0.5 - elif model=='Sammons': - sensitivity=0.75*(0.93*dm_smearing + uw + 0.35*t_res)**-0.5 + if model == 'Quadrature' or model == "StdDev": + sensitivity = totalw ** -0.5 + elif model == 'Sammons': + sensitivity = 0.75 * (0.93 * dm_smearing + uw + 0.35 * t_res) ** -0.5 # implements max integer width cut. if max_meth != 0 and max_iw is not None: - max_w = t_res*(max_iw+0.5) + max_w = t_res * (max_iw + 0.5) if max_meth == 1 or max_meth == 2: - # NOTE: for CRAFT, dm smearing already accounted for prior to width search - # this means that the smearing cut is DM independent if nosmearw > max_w: toolong = np.arange(dm_smearing.size) else: toolong = [] - #toolong = np.where(nosmearw > max_w)[0] elif max_meth == 3 or max_meth == 4: - # NOTE: for CRAFT, dm smearing already accounted for prior to width search toolong = np.where(totalw > max_w)[0] if max_meth == 1 or max_meth == 3: - sensitivity[toolong] = MIN_THRESH # something close to zero + sensitivity[toolong] = MIN_THRESH elif max_meth == 2 or max_meth == 4: - # we have already reduced it by \sqrt{t} - # we thus add a further sqrt{t} factor - sensitivity[toolong] *= (max_w / totalw[toolong])**0.5 + sensitivity[toolong] *= (max_w / totalw[toolong]) ** 0.5 # If model not CHIME, Quadrature or Sammons assume it is a filename else: @@ -1440,22 +1535,27 @@ def calc_relative_sensitivity(DM_frb,DM,w,fbar,t_res,nu_res,Nchan=336,max_idt=No # Should contain DM in the first row and efficiencies in the second row sensitivity_array = np.load(filename) - sensitivity = np.interp(DM, sensitivity_array[0,:], sensitivity_array[1,:], right=1e-2) + sensitivity = np.interp(DM, sensitivity_array[0, :], sensitivity_array[1, :], right=1e-2) return sensitivity -def load_survey(survey_name:str, state:parameters.State, - dmvals:np.ndarray, - zvals:np.ndarray=None, - sdir:str=None, NFRB:int=None, - nbins=None, iFRB:int=0, +def load_survey(survey_name: str, state: parameters.State, + dmvals: np.ndarray, + zvals: np.ndarray = None, + sdir: str = None, NFRB: int = None, + nbins=None, iFRB: int = 0, dummy=False, edir=None, rand_DMG=False, - survey_dict = None, - verbose=False): + survey_dict=None, + verbose=False, + opdir=None, + bPosNum=None, + cluster=False, + clusterRedshift=None, + lensing=False): """Load a survey Args: @@ -1478,6 +1578,11 @@ def load_survey(survey_name:str, state:parameters.State, survey_dict (dict, optional): dictionary of survey metadata to over-ride values in file verbose (bool): print output + opdir (str, optional): directory with pre-computed cluster scattering files + bPosNum (int or str, optional): beam position number for cluster files + cluster (bool, optional): if True, apply cluster scattering. Default False. + clusterRedshift (float, optional): redshift of the cluster. + lensing (bool, optional): if True, apply lensing (reserved). Default False. Raises: IOError: [description] @@ -1522,12 +1627,17 @@ def load_survey(survey_name:str, state:parameters.State, zvals=zvals, NFRB=NFRB, iFRB=iFRB, edir=edir, rand_DMG=rand_DMG, - survey_dict = survey_dict) + survey_dict=survey_dict, + opdir=opdir, + bPosNum=bPosNum, + cluster=cluster, + clusterRedshift=clusterRedshift, + lensing=lensing) return srvy -def vet_frb_table(frb_tbl:pandas.DataFrame, - mandatory:bool=False, - fill:bool=False): +def vet_frb_table(frb_tbl: pandas.DataFrame, + mandatory: bool = False, + fill: bool = False): """ This should not be necessary anymore, since all required FRB data should be populated with @@ -1555,10 +1665,19 @@ def vet_frb_table(frb_tbl:pandas.DataFrame, # These all return p(w) dlogw, and must take as arguments np.log10(widths) def quadrature_convolution(width_function, width_args, scat_function, scat_args, - internal_logvals, bins, backproject = False): + internal_logvals, bins, backproject=False, + xProbScat=None, probScat=None, fractionUnscattered=1.0): ''' - Numerically evaluates the resulting distribution of y=\sqrt{x1^2+x2^2}, - where x1 is the width distribution, and x2 is the scattering distribution. + Numerically evaluates the resulting distribution of y=\\sqrt{x1^2+x2^2[+x3^2]}, + where x1 is the intrinsic width distribution, x2 is the intrinsic scattering + distribution, and the optional x3 is a discrete cluster scattering distribution. + + When ``xProbScat`` is provided (and ``fractionUnscattered < 1``), a fraction + ``(1 - fractionUnscattered)`` of FRBs receive an additional cluster-scattering + term x3 drawn from the discrete distribution defined by ``(xProbScat, probScat)``. + The full histogram is a weighted sum of the unscattered contribution + ``y = sqrt(x1^2 + x2^2)`` and, for each discrete cluster value x3k, the + contribution ``y = sqrt(x1^2 + x2^2 + x3k^2)``. Args: width_function (float function(float,args)): function to call giving p(logw) dlogw @@ -1569,7 +1688,14 @@ def quadrature_convolution(width_function, width_args, scat_function, scat_args, values of log dw to use for internal calculation purposes. bins (np.ndarray([NBINS+1],dtype='float')): bin edges for final width distribution backproject (bool, optional): if True, calculates p(tau|totalw) and p(w|totalw) - for this redshift, and returns additional values. + for this redshift, and returns additional values. Not supported with + cluster scattering (xProbScat is not None). + xProbScat (np.ndarray, optional): 1-D array of cluster scattering values [ms]. + When None (default) no cluster term is added. + probScat (np.ndarray, optional): 1-D array of probabilities corresponding to + each entry of xProbScat. Must sum to 1. Required if xProbScat is given. + fractionUnscattered (float, optional): fraction of FRBs that are *not* scattered + by the cluster (i.e. receive no x3 contribution). Default 1.0 (no cluster). Returns: hist (np.ndarray): histogram of probability within bins wfracs (np.ndarray,only if backproject): p(tau|tw) @@ -1583,99 +1709,82 @@ def quadrature_convolution(width_function, width_args, scat_function, scat_args, # these functions should *not* be normalsid, since some true distribution # may extend beyond the range of interest. But it means that the below functions # absolutely should be correctly normalised - - pw = width_function(internal_logvals, *width_args)*logbinwidth - ptau = scat_function(internal_logvals, *scat_args)*logbinwidth + pw = width_function(internal_logvals, *width_args) * logbinwidth + ptau = scat_function(internal_logvals, *scat_args) * logbinwidth # adds extra bits onto the lowest bin. Does this by integrating in # log space. Assumes exp(-10) is small enough! - lowest = internal_logvals[0] - logbinwidth/2. - extrapw,err = quad(width_function,lowest-10,lowest,args=width_args) - extraptau,err = quad(scat_function,lowest-10,lowest,args=scat_args) + lowest = internal_logvals[0] - logbinwidth / 2. + extrapw, err = quad(width_function, lowest - 10, lowest, args=width_args) + extraptau, err = quad(scat_function, lowest - 10, lowest, args=scat_args) pw[0] += extrapw ptau[0] += extraptau - linvals = 10**internal_logvals + linvals = 10 ** internal_logvals # calculate total widths - all done in linear domain - Nbins = bins.size-1 + Nbins = bins.size - 1 hist = np.zeros([Nbins]) - for i,x1 in enumerate(linvals): - totalwidths = (x1**2 + linvals**2)**0.5 - probs = pw[i]*ptau - h,b = np.histogram(totalwidths,bins=bins,weights=probs) + + # Determine whether cluster scattering is active + use_cluster = (xProbScat is not None) and (fractionUnscattered < 1.0) + + for i, x1 in enumerate(linvals): + # --- Two-term contribution: sqrt(x1^2 + x2^2) --- + # weighted by fractionUnscattered (= 1 if no cluster) + totalwidths = (x1 ** 2 + linvals ** 2) ** 0.5 + probs = pw[i] * ptau * fractionUnscattered + h, b = np.histogram(totalwidths, bins=bins, weights=probs) hist += h + + # --- Three-term contribution: sqrt(x1^2 + x2^2 + x3^2) --- + # Added for the cluster-scattered fraction only + if use_cluster: + scattered_weight = (1.0 - fractionUnscattered) + for k, x3 in enumerate(xProbScat): + totalwidths_sc = (x1 ** 2 + linvals ** 2 + x3 ** 2) ** 0.5 + probs_sc = pw[i] * ptau * scattered_weight * probScat[k] + h_sc, _ = np.histogram(totalwidths_sc, bins=bins, weights=probs_sc) + hist += h_sc # calculate p(w) and p(tau) for each w for this z if backproject: + if use_cluster: + raise NotImplementedError( + "backproject=True is not supported simultaneously with cluster " + "scattering (xProbScat is not None).") # generate arrays to hold probabilities, so values of tau and # w can be fit - wfracs = np.zeros([internal_logvals.size,Nbins]) - taufracs = np.zeros([internal_logvals.size,Nbins]) - - # maps the probabilities as a function of intrinsic - # width to generate a p(w|total_width) and p(tau|total_width) - # note that these are p(observed) values, i.e. after z-correction - # arrays have dimensions(Nw,Ntotal_width) so for each total - # width, we get the probability - for i,x1 in enumerate(linvals): - # total widths corresponding to linvals for tau/iw x1 - totalwidths = (x1**2 + linvals**2)**0.5 + wfracs = np.zeros([internal_logvals.size, Nbins]) + taufracs = np.zeros([internal_logvals.size, Nbins]) + + for i, x1 in enumerate(linvals): + totalwidths = (x1 ** 2 + linvals ** 2) ** 0.5 - # ptau is the probability of achieving linvals, so combined - # probability is pw[i]*ptau - probs = pw[i]*ptau - h,b = np.histogram(totalwidths,bins=bins,weights=probs) - wfracs[i,:] = h + probs = pw[i] * ptau + h, b = np.histogram(totalwidths, bins=bins, weights=probs) + wfracs[i, :] = h - # pw is the probability of achieving linvals, so combined - # probability is ptau[i] * pw - probs = ptau[i]*pw - h,b = np.histogram(totalwidths,bins=bins,weights=probs) - taufracs[i,:] = h + probs = ptau[i] * pw + h, b = np.histogram(totalwidths, bins=bins, weights=probs) + taufracs[i, :] = h # we need p(tau|w). This means sum for a given w must be 1! - wnorm = np.sum(wfracs,axis=0) - # where is p(tau=0 for a given w?) + wnorm = np.sum(wfracs, axis=0) bad = np.where(wnorm == 0.)[0] - wfracs[:,bad] = 1./internal_logvals.size # if a particular width value has no iw, equalise probability over all iw + wfracs[:, bad] = 1. / internal_logvals.size wnorm[bad] = 1. - tnorm = np.sum(taufracs,axis=0) + tnorm = np.sum(taufracs, axis=0) bad = np.where(tnorm == 0.)[0] - taufracs[:,bad] = 1./internal_logvals.size # if a particular width value has no possible tau, equalise probability over all tau + taufracs[:, bad] = 1. / internal_logvals.size tnorm[bad] = 1. - # normalise probabilities for each intrinsic w - #wfracs = (wfracs.T/wnorm).T - #taufracs = (taufracs.T/tnorm).T - wfracs = wfracs/wnorm - taufracs = taufracs/tnorm - - # plot some examples. This code is kept here for internal analysis purposes - if False: - plt.figure() - plt.plot(internal_logvals,ptau,label="Intrinsic p(tau)") - plt.plot(internal_logvals,pw,label="Intrinsic p(w)") - for ib in np.arange(Nbins): - plt.plot(internal_logvals,taufracs[:,ib],label="p(tau) width "+str(ib)) - plt.plot(internal_logvals,wfracs[:,ib],label="p(w) width "+str(ib)) - plt.plot(np.log10([bins[ib],bins[ib]]),[0,1],linestyle=":",color=plt.gca().lines[-1].get_color()) - plt.plot(np.log10([bins[ib+1],bins[ib+1]]),[0,1],linestyle=":",color=plt.gca().lines[-1].get_color()) - break - plt.yscale("log") - #plt.xscale("log") - plt.xlabel("log10 Tau [ms]") - plt.ylabel("p(tau |w)") - plt.legend() - plt.tight_layout() - plt.show() - plt.close() - # exit now, to prevent very many such plots being generated - exit() - - return hist,wfracs,taufracs + wfracs = wfracs / wnorm + taufracs = taufracs / tnorm + + return hist, wfracs, taufracs else: return hist @@ -1696,12 +1805,12 @@ def lognormal(log10w, *args): """ logmean = args[0] logsigma = args[1] - norm = (2.*np.pi)**-0.5/logsigma + norm = (2. * np.pi) ** -0.5 / logsigma result = norm * np.exp(-0.5 * ((log10w - logmean) / logsigma) ** 2) return result -def halflognormal(log10w, *args):#logmean,logsigma,minw,maxw,nbins): +def halflognormal(log10w, *args): """ Generates a parameterised half-lognormal distribution. This acts as a lognormal in the lower half, but @@ -1723,36 +1832,29 @@ def halflognormal(log10w, *args):#logmean,logsigma,minw,maxw,nbins): logmean = args[0] logsigma = args[1] logmax = args[2] - norm = (2.*np.pi)**-0.5/logsigma #Currently no normalisation - if hasattr(log10w,"__len__"): + norm = (2. * np.pi) ** -0.5 / logsigma + if hasattr(log10w, "__len__"): large = np.where(log10w > logmean)[0] - modlogw = np.copy(log10w) # ensures we don't change the original values - modlogw[large] = logmean # subs mean value in for values larger than the mean + modlogw = np.copy(log10w) + modlogw[large] = logmean else: if log10w > logmean: modlogw = logmean else: modlogw = log10w - result = lognormal(modlogw,logmean,logsigma) + result = lognormal(modlogw, logmean, logsigma) - # normalises the distribution. We note that the lower half - # is correctly normalised to 0.5 via the lognormal function - # the upper half spans the range [logmax-logmean] at amplitude - # norm. Hence, the total integral is - # 0.5 + [logmax-logmean]*norm - result /= (0.5 + (logmax-logmean)*norm) + # normalises the distribution. + result /= (0.5 + (logmax - logmean) * norm) return result -def constant(log10w,*args): +def constant(log10w, *args): """ Dummy function that returns a constant of unity, down to a certain minimum, below which it is zero. - NOTE: to include 1+z scaling here, one will need to - reduce the minimum width argument with z. Feature - to be added. Maybe have args also contain min and max values? Args: log10w: log base 10 of widths @@ -1763,10 +1865,10 @@ def constant(log10w,*args): """ width = args[2] - args[0] - if hasattr(log10w,"__len__"): + if hasattr(log10w, "__len__"): good = np.where(log10w > args[0])[0] result = np.zeros([log10w.size]) - result[good] = 1./width + result[good] = 1. / width else: if log10w < args[0]: result = np.array([0])