<jupyter_start><jupyter_text>Generate a mock redshift efficiency fileThis one is meant to approximate the DEEP2 map, shown in Issue 1<jupyter_code>Rab = np.arange(18, 25, 0.05)
#prob_specz = []
#for i in range(1,5):
# frac_i = np.e

def prob_specz(xsamples, psamples):
    """Do a cubic polynomial fit to the given input data
    
    xsamples : array-like
        samples along the 'x' dimension (maybe Rab or SFR or Mass)
    psamples : array-like, same length as xparam
        probability of getting a spec-z for a galaxy, at each xsample value
        
    returns a function to deliver probability for any given value.
    """
    f = interpolate.interp1d(xsamples, psamples, kind='cubic', 
                             axis=- 1, copy=True, bounds_error=False, 
                             fill_value=([0.99], [0.]), 
                             assume_sorted=True)      
    return(f)
rab1    = [19,  20,   21,  22,   23,    24, 24.5]
pspecz1 = [0.9, 0.85, 0.8, 0.75, 0.65, 0.45, 0.15]

p1 = prob_specz(rab1, pspecz1)

plt.plot(Rab, p1(Rab), ls='-', marker='None', color='r')
plt.xlabel(r"R$_{\rm AB}$")
plt.ylabel('prob of getting a reliable specz')<jupyter_output><empty_output><jupyter_text>Read in an existing SNANA HOSTLIB fileThe one we are using here is based on the CANDELS+3DHST galaxy catalogs.<jupyter_code># Big file, 100k galaxies.  Use, e.g. data_end=300 for quick testing
hostlib = table.Table.read("data/candels_hostlib_v4.txt", format='ascii.basic', data_start=1, header_start=0) 

hostlib = table.Table.read("data/Akari_Hosts_subset_SNR_v7.HOSTLIB", 
                           format='ascii.basic', data_start=1, header_start=0) 
#histout = plt.hist2d(hostlib['logssfr'], hostlib['logmass'], bins=20, range=[[-15,-7],[6,13]])

histout = plt.hexbin(hostlib['logssfr'], hostlib['logmass'], gridsize=20, extent=[-15,-7,6,13])
ax = plt.gca()
plt.xlabel('log(sSFR)')
plt.ylabel('log(Mass)')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Connect the (mock) map of redshift efficiency to this map of host galaxy propertiesand fold in the probability of a galaxy hosting a SN to get the probability of getting a specz for SN host galaxies vs z.<jupyter_code>def ssnr_ah17_smooth(logssfr):
    """ Returns the Type Ia specific SN rate per Tyr 
    (number of SN Ia exploding per 10^12 yr per solar mass)
    for a galaxy, using the model of Andersen & Hjorth 2017, which is based 
    on the specific star formation rate, given as log10(SSFR).
    """
    a = 1.2e-13
    b = 0.73
    k = 0.49
    ssfr0 = 1.665e-10
    # logssfr0 = -9.778585762157661    # log10(ssfr0)
    ssfr = 10**logssfr
    ssnr = (a + (a/k) * np.log10(ssfr/ssfr0 + b)) * 1e12
    return(ssnr)
hostlib[:3]<jupyter_output><empty_output><jupyter_text>This is all very kludgy... but as a quick proof-of-concept it seems to be doing the job.  Need to capture this into a module, write some tests, make it more general, then extend it to other parameter dimensions.<jupyter_code>### 1. predict the AB mag of each galaxy in the catalog
rfit = hostlib['sdssr_fit']
logssfr = hostlib['logssfr']
logmass = hostlib['logmass']


### 2. predict the probability of getting a spec-z
pspecz = p1(rfit)
pspecz /= pspecz.sum()

### 3. apply the SN rate estimation to predict SN yields
ssnr = ssnr_ah17_smooth(logssfr)
snr = ssnr * 10**logmass
netsnr = snr.sum()
psnhost = snr / netsnr

# Probability that galaxy is a SN host AND gets a spec-z
wgt = pspecz * psnhost


### 4. count up the fraction of spec-z's for SN hosts vs z, vs Rab, vs H, 
normalization_specz = pspecz.sum()
normalization_snhost = psnhost.sum()

all_binvals, all_binedges = np.histogram(hostlib['ZTRUE'], density=True,
                                         weights=None, bins = np.arange(0,4,0.2))

host_binvals, host_binedges = np.histogram(hostlib['ZTRUE'], 
                                           weights=psnhost, bins = np.arange(0,4,0.2))

specz_binvals, specz_binedges = np.histogram(hostlib['ZTRUE'], 
                                             weights=pspecz, bins = np.arange(0,4,0.2))

hostspecz_binvals, hostspecz_binedges = np.histogram(hostlib['ZTRUE'], 
                                                     weights=pspecz*psnhost, 
                                                     bins = np.arange(0,4,0.2))
Nall = all_binvals.sum()
Nhost = host_binvals.sum()
Nspecz = specz_binvals.sum()
Nhostspecz = hostspecz_binvals.sum()
Nall, Nhost, Nspecz, Nhostspecz
plt.plot(all_binedges[:-1], all_binvals, drawstyle='steps-mid', color='g', label='all gal.')
plt.plot(host_binedges[:-1], host_binvals*2, drawstyle='steps-mid', color='b', label='sn hosts')
#plt.plot(specz_binedges[:-1], specz_binvals, drawstyle='steps-mid', color='r')
plt.plot(hostspecz_binedges[:-1], hostspecz_binvals*4e4, drawstyle='steps-mid', color='k', label='SN host, got specz')
plt.xlabel('redshift')
plt.ylabel('relative count')

plt.legend()<jupyter_output><empty_output><jupyter_text>OLD CODE FOLLOWS BELOW HERE<jupyter_code>#sim1 = table.Table.read("WFIRST_sim_3.8sqdeg_v0.0.fits")
simcatA = table.Table.read("WFIRST_CANDELSbased_sim_0.2sqdeg_v0.0.fits")
simcatB = table.Table.read("fiducial_lensing_galaxia_truth_gal.fits")
simcatBid = np.array(simcatA['id'], dtype='int')
np.max(simcatB['pind']), np.max(simcatBid), len(np.unique(simcatBid))

fits1 = fits.open("WFIRST_sim_3.8sqdeg_v0.0.fits")
fits2 = fits.open("WFIRST_CANDELSbased_sim_0.2sqdeg_v0.0.fits")
fits2dat = table.Table(fits2[1].data)
fits2dat
# Define a dict that translates a fieldname into an integer that matches the 
# 3DHST 'ifield' entry.  Note that AEGIS and EGS are the same field.
ifield_dict = {'aegis':1, 'egs':1, 'cosmos':2, 'goodsn':3, 'goodss':4, 'uds':5}
# Set limits for the allowable Star formation rate and mass ranges
# to be included in the catalogs and weight maps.
_LOGSSFR_MIN=-50
_LOGSSFR_MAX=-5
_LOGMASS_MIN=3
_LOGMASS_MAX=13<jupyter_output><empty_output><jupyter_text>Preamble: computing HST+WFC3 to WFIRST+WFI filter color terms<jupyter_code># read in the filter data, downloaded from http://svo2.cab.inta-csic.es/svo/theory//fps3/
# also available at: https:// wfirst.gsfc.nasa.gov/science/sdt_public/wps/references/ instrument/WFIRST-WFI-Transmission_160720.xlsm
J_hst = ascii.read("HST_WFC3_IR.F125W.dat")
H_hst = ascii.read("HST_WFC3_IR.F160W.dat")
J_wfirst = ascii.read("WFIRST_WFI.J129.dat")
H_wfirst = ascii.read("WFIRST_WFI.H158.dat")

# Define the AB mag zero-points 
# from HST WFC3 data tables www.stsci.edu/hst/wfc3/phot_zp_lbn
# and from Hounsell+ 2018  (which references the wfirst spreadsheet linked above)
zpt_J_hst = 26.25
zpt_H_hst = 25.96
zpt_J_wfirst = 26.35
zpt_H_wfirst = 26.41
# make a plot showing the transmission curves
for filt, filtname in zip([J_hst, H_hst, J_wfirst, H_wfirst],
                          ['HST WFC3 J (F125W)', 'HST WFC3 H (F125W)', 
                          'WFIRST WFI J (J129)', 'WFIRST WFI H (H158)']):
    plt.plot(filt['col1'], filt['col2']/filt['col2'].max(), label=filtname)
plt.legend(loc='lower right')
ax = plt.gca()
ax.set_xlabel(r'Wavelength ($\AA$)')
ax.set_ylabel('Filter Transmission')
# compute the HST to WFIRST filter transformation offsets
def integrate_bandpass(band):
    wave = band['col1']
    transmission = band['col2']/band['col2'].max()
    return(integrate.trapz(transmission, wave))

J_hst_tot = integrate_bandpass(J_hst)
J_wfi_tot = integrate_bandpass(J_wfirst)
H_hst_tot = integrate_bandpass(H_hst)
H_wfi_tot = integrate_bandpass(H_wfirst)

# The magnitude value to add to an HST observed magnitude to change it
# into the equivalent WFIRST magnitude
delta_mAB_J_hst_to_wfirst = (-2.5*np.log10(J_wfi_tot/J_hst_tot) 
                             + zpt_J_wfirst - zpt_J_hst)

delta_mAB_H_hst_to_wfirst = (-2.5*np.log10(H_wfi_tot/H_hst_tot) 
                             + zpt_H_wfirst - zpt_H_hst)<jupyter_output><empty_output><jupyter_text>Read in the master 3DHST catalogThis catalog has two columns that define a unique object identifier:  - ifield : integer number for each field  - 1 : aegis (egs)  - 2 : cosmos  - 3 : goodsn  - 4 : goodss  - 5 : uds- phot_id  = runs from 1 to N for each field, counting phot objects (galaxies)   - Note: max N is about 50,000  We make a unique object ID ('obj_id') by multiplying the ifield by 10^5 and adding the phot_id.<jupyter_code>cat3dhst = table.Table.read("3dhst.v4.1.5.master.fits")
objid_column = table.Column(data=cat3dhst['ifield']*100000 + cat3dhst['phot_id'], name='obj_id')
cat3dhst.add_column(objid_column, index=0)<jupyter_output><empty_output><jupyter_text>Define functions for merging EazyPy and GalFit parameters with the 3DHST data For each field we add EazyPy results (SED fitting parameters) and then GalFit results (Sersic profile parameters) by defining a matching obj_id column and doing an inner join with the 3DHST catalog.<jupyter_code>def get_eazypy_data(fieldname):
    """Read in the EazyPy data (best-fit SED coefficients) for the given field"""
    filepath_data = 'eazypy/{0}_3dhst.v4.1.eazypy.data.fits'.format(fieldname)
    assert os.path.exists(filepath_data)
    datahdu = fits.open(filepath_data)

    # make an object ID that carries the field in the highest-value digit    
    objid_data = datahdu['ID'].data + ifield_dict[fieldname] * 100000
    objid_column = table.Column(data=objid_data, name='obj_id')

    # make a column with the sum of coefficients (sum=0 means no usable fit)
    coeff_sum_data = np.sum(datahdu['COEFFS'].data, axis=1)
    coeff_sum_column = table.Column(data=coeff_sum_data, name='coeff_sum')

    coeff_table = table.Table(datahdu['COEFFS'].data)
    for i in range(len(coeff_table.columns)):
        coeff_table.columns[i].name = 'eazy_coeff{0:02d}'.format(i)
    coeff_table.add_column(coeff_sum_column, index=0)
    coeff_table.add_column(objid_column, index=0)
    return(coeff_table)
def mk_galfit_filepath(fieldname, filtername):
    """Construct the 3DHST catalog file path for the given field"""
    filepath = '{0}/{0}_3dhst.v4.1_{1}.galfit'.format(fieldname, filtername)
    assert os.path.exists(filepath)
    return(filepath)

def mk_merged_catalog(fieldname, trim_unusable=True, checkradec=True,
                     verbose=True):
    """Do an inner merge of a GalFit catalog with the 3DHST catalog, then do an 
    inner merge of an EazyPy coefficients table with the merged GalFit+3DHST catalog.
    
    trim_unusable : exclude galaxies with bad sersic 
                    fits, no redshift, or no SSFR measure
                    
    checkradec : do a consistency check on the merge by verifying 
                 consistent coordinate positions.
    """
    ifield = ifield_dict[fieldname]

    # Get the EazyPy coefficients data
    cateazypy = get_eazypy_data(fieldname)
    
    # Get the GalFit data
    galfit_f125w_filepath = mk_galfit_filepath(fieldname, 'f125w')
    galfit_f160w_filepath = mk_galfit_filepath(fieldname, 'f160w')
    
    catgalfit = table.Table.read(galfit_f160w_filepath, format='ascii.commented_header')
    objid_column = table.Column(data=ifield*100000 + catgalfit['NUMBER'], name='obj_id')
    catgalfit.add_column(objid_column, index=0)
    catgalfit.rename_column('RA', 'ra_galfit')
    catgalfit.rename_column('DEC', 'dec_galfit')
    catgalfit.rename_column('mag', 'mag_f160w')
    catgalfit.rename_column('dmag', 'dmag_f160w')

    # Add a column the observed F125W mag from the Sersic fits
    # (we already have the F160W mags)
    cat125 = table.Table.read(galfit_f125w_filepath, format='ascii.commented_header')
    cat125.rename_column('mag', 'mag_f125w')
    cat125.rename_column('dmag', 'dmag_f125w')
       
    i160 = catgalfit.index_column('mag_f160w')
    catgalfit.add_columns([cat125['mag_f125w'],cat125['dmag_f125w']], indexes=[i160, i160] )
    

    # Add columns with the WFIRST J129 and H158 mags, 
    # using the HST-to-WFIRST offsets computed above
    col129 = table.Column(data=cat125['mag_f125w']+delta_mAB_J_hst_to_wfirst,
                         name='mag_j129')
    col158 = table.Column(data=catgalfit['mag_f160w']+delta_mAB_H_hst_to_wfirst,
                         name='mag_h158')
    catgalfit.add_columns([col129, col158], indexes=[i160, i160] )

    
    # Merge the 3dHST and GalFit and EazyPy catalogs with an inner join, matching 
    # on the obj_id column, which is the Skelton 2014 galaxy id. 
    catmerged0 = table.join(cat3dhst, catgalfit, join_type='inner', keys=['obj_id'])
    catmerged = table.join(catmerged0, cateazypy, join_type='inner', keys=['obj_id'])
    
    if trim_unusable:
        # limit to only galaxies that have good galfit data (flag=0) 
        # good galaxy properties (sfr, mass) and good
        # redshifts (zbest>0), so they'll deliver a good SN rate.
        ivalid = np.where((catmerged['f']<1) & 
                          (catmerged['z_best']>0) & 
                          (catmerged['lssfr']>_LOGSSFR_MIN) &
                          (catmerged['lssfr']<_LOGSSFR_MAX) &
                          (catmerged['lmass']>_LOGMASS_MIN) &
                          (catmerged['lmass']<_LOGMASS_MAX) & 
                          (catmerged['coeff_sum']>0)
                         )[0]
        if verbose:
            print("Selected {0} of {1} galaxies".format(
                len(ivalid), len(catmerged)) + 
              "  with ok redshifts, Sersic fits, sSFR, and mass data.")
        catmerged = catmerged[ivalid]
        
    if checkradec:
        # If all of our objects are positional matches, then none should 
        # be offset by more than 0.01 arcsec in either direction
        delta_ra_max_arcsec = np.max( 
            np.abs(catmerged['ra'] - catmerged['ra_galfit']) 
            * np.cos(catmerged['dec'])) * 3600.
        delta_dec_max_arcsec = np.max(
            np.abs(catmerged['dec'] - catmerged['dec_galfit'])) * 3600.
        assert delta_ra_max_arcsec < 0.06
        assert delta_dec_max_arcsec < 0.06
        if verbose:
            print("{} galfit catalog merged with 3DHST catalog successfully.".format(fieldname))
            print("Max RA discrepancy = {} arcsec".format(delta_ra_max_arcsec))
            print("Max Dec discrepancy = {} arcsec".format(delta_dec_max_arcsec))
            
    return(catmerged)<jupyter_output><empty_output><jupyter_text>Define functions for Computing SN rates for each galaxyAnd show us a plot of the function<jupyter_code># Show the SN Rate model from Andersen & Hjorth 2017
#cat = aegis_cat
#cat = cosmos_cat
#lssfrcat = cat['lssfr']
#ssnrcat = snr_ah17_smooth(lssfrcat)

lssfr_testvals = np.linspace(-13., -5., 1000)
ssnr = ssnr_ah17_smooth(lssfr_testvals)

#plt.semilogy(lssfrcat, ssnrcat, marker='.', ms=2, alpha=0.5, ls=' ', color='r')
plt.semilogy(lssfr_testvals, ssnr, marker=' ', ls='-', lw=1, color='k')
ax = plt.gca()
ax.set_ylabel("specific SN Rate [Tyr-1 Msun-1]")
ax.set_xlabel("log10( sSFR [yr-1] )")<jupyter_output><empty_output><jupyter_text>Define a function to add a SN rate column to a catalogThis predicted SN Rate per galaxy will define the weight in the WGTMAP for a SNANA HOSTLIB file.We also trim out bad entries from the catalog (missing redshift, Sersic fit, or info required for the SN rate).<jupyter_code>def mk_weight_column(catalog, snr_model='AH17', verbose=True):
    """Add a column to the given catalog that gives the relative probability
    for each galaxy hosting a SN in any given observer-frame year.  This is
    computed based on the predicted SN rate (number of SN explosions per 
    observer-frame year) of each galaxy, adopting the specified
    SN rate model.  
        
    snr_model='A+B' : SNR = A*M + B*SFR   (Scannapieco & Bildsten 2005)
             ='AH17' : the smooth logarithmic sSFR model (Andersen & Hjorth 2017)
             
    """
    if snr_model.lower()=='a+b':
        # Note: adopting the A and B values from Andersen & Hjorth 2017
        # but dividing by 1e-4 (so the SNR below actually counts the number of 
        # SN explodiing per 10000 yrs)
        A = 4.66 * 1e-10
        B = 4.88  
        snr = A * 10**catalog['lmass'] + B * 10**catalog['lsfr']
        # divide by the total snr to get relative probabilities
        snr /= np.nanmax(snr)
        snrcol = table.Column(data=snr, name='snr_A+B')
        if 'snr_A+B' in catalog.colnames:
            catalog['snr_A+B'] = snr
        else:
            catalog.add_column(snrcol)
    elif snr_model.lower().startswith('ah17'):
        a = 1.2e-13
        b = 0.73
        k = 0.49
        ssfr0 = 1.665e-10
        ssfr = 10**(catalog['lsfr']-catalog['lmass'])
        ssnr = a + (a/k) * np.log10(ssfr/ssfr0 + b)
        snr = ssnr * 10**catalog['lmass']
        snr /= np.nanmax(snr)
        snrcolname = 'snr_AH17_smooth'
        snrcol = table.Column(data=snr, name=snrcolname)
        if snrcolname in catalog.colnames:
            catalog[snrcolname] = snr
        else:
            catalog.add_column(snrcol)
   
    if verbose:
        print("Added/updated SN rate weight column using {} model".format(snr_model))      

    return(catalog)<jupyter_output><empty_output><jupyter_text>Converting to a HOSTLIBExtract just the useful information (useful to SN sims) for each galaxy and write it out to a SNANA HOSTLIB file. From the 3DHST catalog we get redshift, position and magnitude information: - z : best redshift from the 3DHST catalog- RA, Dec- observed AB mags in SDSS u,g,r,i,z and 2MASS J, HThe GalFit catalogs from Arjen Van Der Weil contain the following columns:- COL 1,  NUMBER: ID from the Skelton et al. (2014) photometric catalogs from the 3D-HST team- COL 2,  RA: from SExtractor, measured in F160W- COL 3,  DEC: idem.- COL 4,  f: FLAG value (0: good fit; 1: suspicious fit; 2: bad fit; 3: no fit -- see van der Wel et al. 2012)- COL 5,  mag: total AB magnitude from best-fitting Sersic model (GALFIT)- COL 6,  dmag: 1-sigma uncertainty in mag- COL 7,  re: semi-major axis in arcsec of the ellipse that contains half of the total light in the best fitting Sersic model- COL 8,  1-sigma uncertainty on re- COL 9,  n: Sersic index of the best-fitting Sersic model- COL 10, dn: 1-sigma uncertainty on n- COL 11, q: axis ratio of the best fitting Sersic model- COL 12, dq: 1-sigma uncertainty on q- COL 13, pa: position angle in degrees (0: North; 90: East)- COL 14, dpa: 1-sigma uncertainty on pa- COL 15, sn: S/N as measured in the F???W filter, integrated over the F160W segmentation regionFor the HOSTLIB file we get Sersic fit info, and HST J and H band magnitudes from this catalog. There are two sets of functions:  - one set for defining the HOSTLIB entries (a header and a line for each galaxy)- another set for making a weight map in the HOSTLIB (defining the weight based on the predicted SN rate per galaxy).<jupyter_code>_HOSTLIB_HEADER = """
VARNAMES: GALID RA_GAL  DEC_GAL      ZTRUE      ZERR     ZPHOT  ZPHOTERR    logsfr   logmass  logssfr hstj_obs  hsth_obs n0_Sersic a0_Sersic b0_Sersic a_rot eazy_coeff00 eazy_coeff01 eazy_coeff02 eazy_coeff03 eazy_coeff04 eazy_coeff05 eazy_coeff06 eazy_coeff07 eazy_coeff08 eazy_coeff09 eazy_coeff10 eazy_coeff11 eazy_coeff12

"""

    
def mk_hostlib_line(linedat):
    """ Extract galaxy info  that is useful for SN simulations from a 
    merged 3DHST+GalFit catalog.  Return a string in the SNANA HOSTLIB format
    """
    # question for <NAME>: can a SNANA HOSTLIB handle asymmetric z errors? 
    zbesterr = (linedat['z_best_u68']-linedat['z_best_l68'])/2.   # assuming 'No.'
    zphoterr = (linedat['z_phot_u68']-linedat['z_phot_l68'])/2.   # assuming 'No.'

    # Galfit: pa: position angle in degrees (0: North; 90: East)
    # SNANA:  “a_rot” is the rotation angle (in degrees) of the major axis w.r.t. the +RA coordinate. 
    # If North is up and East is to the right, a_rot is measured clockwise, from the East toward the South.
    a_rot = linedat['pa'] - 90 
    if a_rot < 0 : 
        a_rot += 180

    outline = ("GAL: {:1} " + 2*"{:10.5f}"+ 7*"{:10.4f}" + 
               2*"{:8.2f}"+ 4*"{:10.4f}" + 
               13*"{:12.5e}" + "\n").format(
        linedat['obj_id'], linedat['ra'], linedat['dec'], 
        linedat['z_best'], zbesterr, linedat['z_peak_phot'], zphoterr, 
        linedat['lsfr'], linedat['lmass'], linedat['lssfr'],
        #-2.5*np.log10(linedat['L156'])+25, 
        #-2.5*np.log10(linedat['L157'])+25, 
        #-2.5*np.log10(linedat['L158'])+25, 
        #-2.5*np.log10(linedat['L159'])+25, 
        #-2.5*np.log10(linedat['L160'])+25, 
        #-2.5*np.log10(linedat['L161'])+25, 
        #-2.5*np.log10(linedat['L162'])+25, 
        linedat['mag_j129'], linedat['mag_h158'],
        linedat['n'], linedat['re'], linedat['re']*linedat['q'], a_rot, 
        linedat['eazy_coeff00'], linedat['eazy_coeff01'], linedat['eazy_coeff02'], 
        linedat['eazy_coeff03'], linedat['eazy_coeff04'], linedat['eazy_coeff05'], 
        linedat['eazy_coeff06'], linedat['eazy_coeff07'], linedat['eazy_coeff08'], 
        linedat['eazy_coeff09'], linedat['eazy_coeff10'], linedat['eazy_coeff11'], 
        linedat['eazy_coeff12']
    )
    return(outline)

# First option for a weight map: 
#   use the sSFR to define the weight and host mass to define the SN mag shift. 
#  The weight map is defined with a few dozen points and SNANA is left to
# interpolate between them.

_WGTMAP_HEADER_SSFR = """
VARNAMES_WGTMAP: logssfr logmass  WGT SNMAG_SHIFT

"""

def mk_wgtmap_block(snr='AH17_smooth', 
                    logssfr_stepsize=0.5, 
                    logmass_stepsize=0.5):
    """ 
    Construct the HOSTLIB weight map: a block of text, with 
    each line including N-2 observable host galaxy parameters, and the 
    last two giving WGT and SNMAGSHIFT.  
    
    WGT is the weight (relative probability of hosting a SN Ia)
    assigned for the preceding set of galaxy parameters.  
    SNMAGSHIFT is the magnitude shift applied toa SN with the matching
    set of host galaxy parameters. 
    
    SNANA will do interpolation between the host galaxy parameter 
    values given in the WGTMAP to assign a WGT and SNMAGSHIFT to each 
    simulated SN host galaxy in the HOSTLIB section, which follows 
    below the WGTMAP. 
    """
    wgtmap_str = '\n\n'
    logssfr_gridpoints = np.arange(_LOGSSFR_MIN,
                                   _LOGSSFR_MAX + logssfr_stepsize, 
                                   logssfr_stepsize)
    logmass_gridpoints = np.arange(_LOGMASS_MIN,
                                   _LOGMASS_MAX + logmass_stepsize, 
                                   logmass_stepsize)
    massstep_values = np.where(logmass_gridpoints<10, 0.0, 0.04)
    for lssfr in logssfr_gridpoints:
        ssnr = ssnr_ah17_smooth(lssfr)
        for imass in range(len(logmass_gridpoints)):
            wgt = ssnr * np.power(10,logmass_gridpoints[imass])
            wgtmap_str += "WGT: {:8.3f} {:8.3f} {:8.3f} {:6.2f}\n".format(
                lssfr, logmass_gridpoints[imass], wgt, massstep_values[imass]
            )
            
    wgtmap_str += '\n\n'
    return(wgtmap_str)


   
# The second option : specify a weight and mag shift for each galaxy individually
# so the weight map is defined point by point for each unique GALID 
_WGTMAP_HEADER_GALID = """
VARNAMES_WGTMAP: GALID WGT SNMAG_SHIFT

"""  
    
def mk_wgtmap_galid(linedat, snrcolname='snr_AH17_smooth'):
    """ Construct a weight line for the given single line of galaxy data.  
    Return a string in the SNANA HOSTLIB WGTMAP format
    """
    outline = ("WGT: {:1} {:14.5e}  {:10.5f} \n").format(
        linedat['obj_id'], linedat[snrcolname], 0.0)
    return(outline)<jupyter_output><empty_output><jupyter_text>These two functions bring it all together1. construct a weight map file that can apply to any CANDELS/3DHST-based HOSTLIB2. constructing a HOSTLIB file from an input merged catalog<jupyter_code>def mk_wgtmap_file(wgtmapfilename):
    """Make a SNANA HOSTLIB WGTMAP file that defines a parameterized
    weightmap using the Andersen & Hjorth 2017 SN rate model, which 
    predicts the SN rate for each galaxy based on that galaxy's sSFR.
    """    
    fout = open(wgtmapfilename, mode='w')
    fout.write(_WGTMAP_HEADER_SSFR)
    wgtmap_block = mk_wgtmap_block()
    fout.write(wgtmap_block)
    fout.close()
    return

def mk_hostlib_file(cat, hostlib_filename, clobber=True, verbose=True):
    """Construct a SNANA HOSTLIB file from the given catalog, which must be the 
    result of merging a 3DHST and a GalFit catalog. 
    """
    if not clobber and os.path.exists(hostlib_filename):
        print("{} exists. Not clobbering.".format(hostlib_filename))
        return(-1)    
        
    fout = open(hostlib_filename, mode='w')
    fout.write(_HOSTLIB_HEADER)
    for i in range(len(cat)):
        fout.write(mk_hostlib_line(cat[i]))
    fout.close()
    return<jupyter_output><empty_output><jupyter_text>Interlude:  checking on some of the peculiarly large high-z galaxies<jupyter_code>fieldname = 'cosmos'
cat_trim = mk_merged_catalog(fieldname)
cat_trim = mk_weight_column(cat_trim, snr_model='AH17')<jupyter_output>Selected 18499 of 33879 galaxies  with ok redshifts, Sersic fits, sSFR, and mass data.
cosmos galfit catalog merged with 3DHST catalog successfully.
Max RA discrepancy = -0.0 arcsec
Max Dec discrepancy = 0.0018000000002516003 arcsec
Added/updated SN rate weight column using AH17 model<jupyter_text>identify galaxies that are unusually large<jupyter_code>FLCDM = FlatLambdaCDM(H0=70, Om0=0.3)
z = np.array(cat_trim['z_best'])
mu = FLCDM.distmod(z).value
a0 = np.array(cat_trim['re'])  # Sersic fit major axis length in arcsec
M = cat_trim['mag_h158']-mu
a0_thresh = 3
ibig = np.where(a0>a0_thresh)[0]
def plot_hexbin_scatter(ax, xval, yval, xlabel='x values', ylabel='y values'):
    """Plot yval vs xval as a combined hexbin and scatter plot, highlighting in red 
    the galaxies that are unusually large"""
    ax.plot(xval, yval, color='k', ms=2, marker='.', ls=' ', zorder=-10)
    hexbinout = ax.hexbin(xval, yval, gridsize=30, cmap=cm.Greys_r, mincnt=20, zorder=10)
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)
    ax.plot(xval[ibig], yval[ibig], color='r', ms=3, 
            marker='.', ls=' ', zorder=10)
    return(hexbinout)
fig = plt.figure(figsize=[12,4])
ax1 = fig.add_subplot(121)

xval = M
yval = np.log10(a0)
plot_hexbin_scatter(ax1, xval, yval, xlabel='abs. mag.', ylabel='log10(a0)')
ax1.axhline(np.log10(a0_thresh), color='r', lw=1)

ax2 = fig.add_subplot(122)
xval = z
yval = M
hexbinout = plot_hexbin_scatter(ax2, xval, yval, xlabel='redshift', ylabel='abs. mag.')
ax.plot(np.log10(z), M-np.log10(z)*6, marker='.', ls=' ',ms=2)

plt.colorbar(hexbinout)<jupyter_output><empty_output><jupyter_text>Show thumbnails for some unusually large galaxies<jupyter_code>from astropy import coordinates
from astropy import units as u
# Read in the COSMOS mosaic
mosaic = fits.open("mosaics/hlsp_candels_hst_wfc3_cos-tot_f160w_v1_0_drz.fits")
wcsobj = WCS(mosaic[0].header)
ra_big = np.array(cat_trim['ra'][ibig])*u.degree
dec_big = np.array(cat_trim['dec'][ibig])*u.degree
x_big, y_big = wcsobj.wcs_world2pix(ra_big, dec_big, 1)


ra = np.array(cat_trim['ra'])*u.degree
dec = np.array(cat_trim['dec'])*u.degree
x_all, y_all = wcsobj.wcs_world2pix(ra, dec, 1)
i=ibig[10]
x,y = int(x_all[i]), int(y_all[i])
dxy=50 # half of the stamp size
mosaic_stamp = mosaic[0].data[y-dxy:y+dxy,x-dxy:x+dxy]
plt.imshow(mosaic_stamp, vmin=-0.01, vmax=0.03)<jupyter_output><empty_output><jupyter_text>Define functions for Summmary plots from a trimmed catalog<jupyter_code>def plot_weighted_redshift_hist(z, weights, axnum=1, binsize=0.2, label='', 
                                ncol=4, nrow=2, fig=None, **kwargs):    
    if fig is None:
        fig = plt.gcf()
    ax = fig.add_subplot(nrow, ncol, axnum)
    histout = ax.hist(z, bins=np.arange(0,5,binsize), range=(0,5), weights=weights, 
                      **kwargs)
    ax.set_xlabel('redshift')
    ax.set_ylabel(label)
    return(ax)

def plot_catalog_summary(cat, catalog_name):
    """
    Make some plots showing summaries of the SN host galaxy catalog data
    """
    fig = plt.figure(figsize=[12,7])
    fig.subplots_adjust(wspace=0.4, hspace=0.4)

    z = cat['z_best']
    snr = cat['snr_AH17_smooth']
    sfr = 10**cat['lsfr']
    mass = 10**cat['lmass']    

    # total SNR, binned by redshift
    ax2 = plot_weighted_redshift_hist(z, snr, 
                                      axnum=2, color='b', histtype='step', lw=2, 
                                      label = 'SN Rate [yr-1]')

    # total SFR, binned by redshift
    ax3 = plot_weighted_redshift_hist(z, sfr, axnum=3, color='darkorange', histtype='step', lw=2, 
                                      label = 'SFR [Msun yr-1]')


    # total mass, binned by redshift
    ax4 = plot_weighted_redshift_hist(z, mass, axnum=4, color='g', histtype='step', lw=2, 
                                      label = 'Galaxy Mass [Msun]')

    # Histogram of redshifts
    ax5 = fig.add_subplot(2, 4, 5)
    ax5.set_ylabel('Galaxy Count')
    histout = ax5.hist(z, bins=np.arange(0,5,0.2), color='r')
    ax5.set_xlabel('redshift')


    # histogram of galaxy masses ( log10(mass) )
    ax5 = fig.add_subplot(2, 4, 8)
    histout = ax5.hist(cat['lmass'], bins=np.arange(5,13,0.2), color='g')
    ax5.set_xlabel('log10(Mass/Msun)')

    # histogram of Sersic index
    ax6 = fig.add_subplot(2, 4, 7)
    histout = ax6.hist(cat['n'], color='m')#, bins=np.arange(5,13,0.2), color='0.5')
    ax6.set_xlabel('Sersic Index')

    ax7 = fig.add_subplot(2, 4, 6)
    histout = ax7.hist(cat['snr_AH17_smooth'], bins=np.linspace(0,0.05,50), color='b')
    ax7.set_yscale('log')
    ax7.set_xlabel('SN Rate')


    catalog_info = """
{0}\n
N = {1:d}
""".format(catalog_name, len(cat))
    fig.text(0.1, 0.9, catalog_info, transform=fig.transFigure, ha='left', va='top',
            fontsize=18)
def plot_mag_hist_at_z(cat, zval, band='V_rest', dz=0.1, ax=None):
    if ax is None:
        ax = plt.gca()
    iz = np.where((cat['z_best']<zval+dz) & (cat['z_best']>zval-dz))[0]

    if band=='U_rest':
        m = 25. - 2.5*np.log10(cat['L135'][iz])
    elif band=='V_rest':
        m = 25. - 2.5*np.log10(cat['L137'][iz])
    elif band=='I_rest':
        I = 25. - 2.5*np.log10(cat['L139'][iz])
    elif band=='u':
        m = 25. - 2.5*np.log10(cat['L156'])[iz]
    elif band=='r':
        r = 25. - 2.5*np.log10(cat['L158'])[iz]
    elif band=='z':
        z = 25. - 2.5*np.log10(cat['L160'])[iz]

        #plt.scatter(u-r, r-z, marker='.', color='k', alpha=0.2)
    histout = ax.hist(m, bins=np.arange(18,30,0.5))
    ax.text(0.05, 0.95, "z={:.2f}+-{:.2f}".format(zval,dz), transform=ax.transAxes,
           ha='left', va='top', fontsize=14)
    return
    
def plot_restV_histograms(cat, field_name):
    fig = plt.figure(figsize=[12,8])
    #fig.subplots_adjust(wspace=0.4, hspace=0.2)
    ax1 = fig.add_subplot(2, 4, 1)

    for i in range(1,9):
        if i>1:
            ax = fig.add_subplot(2, 4, i, sharey=ax1)
        else:
            ax = ax1
        plot_mag_hist_at_z(cat, zval=0.4*i, dz=0.2, band='V_rest', ax=ax)
    #ax1.set_ylim(0,800)
    fig.text(0.5, -0.05, "Apparent Magnitude in rest-frame V band", 
             transform=fig.transFigure, ha='center', va='bottom', fontsize=18)
    fig.text(-0.05, 0.5, "Number of Galaxies in {} Catalog".format(field_name),
             transform=fig.transFigure, ha='left', va='center', 
             fontsize=18, rotation=90)
    plt.tight_layout()
    return<jupyter_output><empty_output><jupyter_text>Make the files for SNANA and summary plots for All 5 CANDELS Fields * First make a weight map file (common for all fields)* then for each  merged catalog:   * compute SN rates and trim bad entries  * construct the HOSTLIB file<jupyter_code>mk_wgtmap_file('wgtmap_for_candels3dhst_hostlibs.txt')
catlist = []
for fieldname in ['aegis', 'goodss', 'goodsn', 'cosmos', 'uds']:
    cat_trim = mk_merged_catalog(fieldname)
    cat_trim = mk_weight_column(cat_trim, snr_model='AH17')
    catlist.append(cat_trim)
    mk_hostlib_file(cat_trim, '{}_hostlib.txt'.format(fieldname))
    print('\n')
    plot_catalog_summary(cat_trim, fieldname.upper())
    plt.savefig("{}_catalog_summary.png".format(fieldname))

    plot_restV_histograms(cat_trim, fieldname)
    plt.savefig("{}_catalog_mag_hist.png".format(fieldname))
cat_trim = table.vstack(catlist, join_type='exact')
fieldname = 'candels'

mk_hostlib_file(cat_trim, '{}_hostlib.txt'.format(fieldname))
print('\n')

plot_catalog_summary(cat_trim, fieldname.upper())
plt.savefig("{}_catalog_summary.png".format(fieldname))

plot_restV_histograms(cat_trim, fieldname)
plt.savefig("{}_catalog_mag_hist.png".format(fieldname))<jupyter_output><jupyter_text>Write out the EazyPy spectral templatesExtract the EazyPy template data from any of the EazyPy catalog files, and write out as separate ascii data files.<jupyter_code>def mk_eazypy_template_datafile(outfilename='eazy_13_spectral_templates.dat'):
    fieldname = 'uds'  # any field will do. 
    filepath_data = 'eazypy/{0}_3dhst.v4.1.eazypy.data.fits'.format(fieldname)
    assert os.path.exists(filepath_data)
    hdu = fits.open(filepath_data)

    templ_wave = hdu['TEMPL'].data # rest-frame wavelength, Angstroms
    templ_flux = hdu['TEMPF'].data # template flux density for AB zeropoint = 25

    wave_column = table.Column(data=templ_wave, name='wavelength_angstroms')
    flux_columns = [table.Column(data=templ_flux[i], 
                                 name='eazyspectemp{:02d}'.format(i))
                    for i in range(len(templ_flux))]

    column_list = [wave_column] + flux_columns
    templ_data = table.Table(column_list)
    templ_data.write(outfilename, format='ascii.commented_header', overwrite=True)
    print("Wrote out Eazy spectral templates to {}".format(outfilename))
    return(templ_data)
templ_data = mk_eazypy_template_datafile()<jupyter_output>Wrote out Eazy spectral templates to eazy_13_spectral_templates.dat
{'fraction_non_alphanumeric': 0.1225169466548777, 'fraction_numerical': 0.0350132625994695, 'mean_word_length': 3.4112064482579303, 'pattern_counts': {'":': 0, '<': 79, '<?xml version=': 0, '>': 79, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '111864', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 13183, 'n_tokens_neox': 12162, 'n_words': 3000, 'path': 'romanz/roman_sn_hostlib.ipynb', 'repo_name': 'AnomandarisPurake/roman-sn-redshifts', 'size': '33936'}