<jupyter_start><jupyter_text>MILLIQUAS Plots<jupyter_code>print(ra1_first.min(),ra1_first.max()) #smaller region
print(dec1_first.min(),dec1_first.max())
print(ra2_first.min(),ra2_first.max()) #larger region
print(dec2_first.min(),dec2_first.max())
ra_min1= Angle(f'{ra1_first.min()}d', unit='hour')
ra_max1= Angle(f'{ra1_first.max()}d', unit='hour')
ra_min2= Angle(f'{ra2_first.min()}d', unit='hour')
ra_max2= Angle(f'{ra2_first.max()}d', unit='hour')
dec_min1 = Angle(f'{dec1_first.min()}d', unit='degree')
dec_max1 = Angle(f'{dec1_first.max()}d', unit='degree')
dec_min2 = Angle(f'{dec2_first.min()}d', unit='degree')
dec_max2 = Angle(f'{dec2_first.max()}d', unit='degree')

print(ra_min1,'|', ra_max1)
print(dec_min1,'|', dec_max1)
print('\n')
print(ra_min2,'|', ra_max2)
print(dec_min2,'|', dec_max2)
MILLI_sub1 = pd.read_csv('/epyc/users/tw94/lsst-hits-bump/mrawls/MILLI_QSO_region1.1584579832', 
                         sep = '|',header=2, usecols = [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18], na_values = '     ')
MILLI_sub1.head()
for col in MILLI_sub1.columns: 
    print(col) 
MILLI_sub1 = MILLI_sub1.rename(columns={"ra         ":"ra", "dec        ":"dec", "bmag ":"bmag", "rmag ":"rmag",
                                        "radio_name            ":"radio_name", "xray_name             ":"xray_name",
                                        "lii      ":"lii", "bii     ":"bii", "class           ":"class"})
MILLI_sub2 = pd.read_csv('/epyc/users/tw94/lsst-hits-bump/mrawls/MILLI_QSO_region2.1584579992', 
                         sep = '|',header=2, usecols = [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18], na_values = '     ')
MILLI_sub2.head()
MILLI_sub2 = MILLI_sub2.rename(columns={"ra         ":"ra", "dec        ":"dec", "bmag ":"bmag", "rmag ":"rmag",
                                        "radio_name            ":"radio_name", "xray_name             ":"xray_name",
                                        "lii      ":"lii", "bii     ":"bii", "class           ":"class"})
MILLI_sub1['ra_deg']= (Angle(MILLI_sub1['ra'], unit ='hour')).deg
MILLI_sub1['dec_deg'] = (Angle(MILLI_sub1['dec'], unit ='degree')).deg
MILLI_sub1.head()
MILLI_sub2['ra_deg']= (Angle(MILLI_sub2['ra'], unit ='hour')).deg
MILLI_sub2['dec_deg'] = (Angle(MILLI_sub2['dec'], unit ='degree')).deg
MILLI_sub2.head()
cols = MILLI_sub2.columns.tolist()
cols_reorder = ['ra', 'dec', 'ra_deg', 'dec_deg', 'bmag', 'rmag', 'redshift', 'radio_name', 'xray_name', 
                'lii', 'bii', 'broad_type', 'optical_flag', 'red_psf_flag', 'blue_psf_flag', 'ref_name', 
                'ref_redshift', 'qso_prob', 'class'] 
MILLI_sub1 = MILLI_sub1[cols_reorder]
MILLI_sub2 = MILLI_sub2[cols_reorder]

MILLI_sub1.head()<jupyter_output><empty_output><jupyter_text>Preliminary Plots<jupyter_code>plt.figure(figsize = (10,10))
plt.scatter(MILLI_sub1['ra_deg'].values,MILLI_sub1['dec_deg'].values)
plt.figure(figsize = (10,10))
plt.scatter(MILLI_sub2['ra_deg'].values,MILLI_sub2['dec_deg'].values)<jupyter_output><empty_output><jupyter_text>Filtering to Visable and Plotting Objects on the Sky<jupyter_code>MILLI_sub1_visible = MILLI_sub1.query('rmag < 23 and bmag < 23')
MILLI_sub2_visible = MILLI_sub2.query('rmag < 23 and bmag < 23')
print(len(MILLI_sub1), len(MILLI_sub1_visible))
print(len(MILLI_sub2), len(MILLI_sub2_visible))
def plot_objects_on_sky2(ra1_first, dec1_first,
                        ra2_first, dec2_first,
                        ra1_second, dec1_second,
                        ra2_second, dec2_second,
                        ra1_third, dec1_third,
                        ra2_third, dec2_third):
    """This function takes three sets of RA and Dec and plots them
    both on the sky in different colors.
    
    It is customized to plot a specific region (three HiTS fields in two panels).
    
    "1" and "2" refer to the two panels in the plot.
    "first" and "second" refer to the two different datasets.
    """
    # Set up the figure object, title, and two axes
    fig = plt.figure(figsize=(12, 12))
    fig.suptitle('View of Objects from LSST, HiTS (QSOs), and MILLIQUAS on the Sky', x = 0.5, y= 0.9, fontsize=15)
    ax1 = plt.subplot2grid((100, 100), (0, 55), rowspan=50, colspan=45)
    ax2 = plt.subplot2grid((100, 100), (0, 0), rowspan=90, colspan=50)

    # Plot the first set of RAs and Decs in blue
    # This will be from the LSST database
    ax1.scatter(ra1_first, dec1_first, marker='.', s=0.5, alpha=0.3, c='white',label='LSST Objects')
    ax2.scatter(ra2_first, dec2_first, marker='.', s=0.5, alpha=0.3, c='white',label='LSST Objects')
    
    # Plot the second set of RAs and Decs in red
    # This will be from the HiTS DR1
    ax1.scatter(ra1_second, dec1_second, marker='.', s=17, alpha=1, c='cyan',label='HiTS QSO')
    ax2.scatter(ra2_second, dec2_second, marker='.', s=17, alpha=1, c='cyan',label='HiTS QSO')

    ax1.invert_xaxis()
    ax2.invert_xaxis()
    
    # Plot the third set of RAs and Decs in purple
    # This will be from the MILLIQUAS database
    ax1.scatter(ra1_third, dec1_third, marker = '.', s = 15, alpha = 0.6, c = 'blue', 
                label = 'MILLIQUAS QSO')
    ax2.scatter(ra2_third, dec2_third, marker = '.', s = 15, alpha = 0.6, c = 'blue', 
                label = 'MILLIQUAS QSO')
    
    
    # Label axes
    plt.xlabel('RA (deg)', fontsize=14)
    plt.ylabel('Dec (deg)', fontsize=14)
    
    # Increase x and y tick text size
    ax1.tick_params(axis='both', which='major', labelsize=11.5)
    ax1.tick_params(axis='both', which='minor', labelsize=11.5)

    ax2.tick_params(axis='both', which='major', labelsize=11.5)
    ax2.tick_params(axis='both', which='minor', labelsize=11.5)
    
    # Generate legends
    lgnd = ax1.legend(scatterpoints=1, fontsize=10, prop={'size': 11.5},bbox_to_anchor=(.69,-.38), loc = 'best')
    for lh in lgnd.legendHandles: 
        lh.set_alpha(1)
    lgnd.legendHandles[0]._sizes = [100]
    lgnd.legendHandles[0].set_color('white')
    lgnd.legendHandles[1]._sizes = [100]
    lgnd.legendHandles[2]._sizes = [100]
    
    #lgnd2 = ax2.legend(loc="lower right", scatterpoints=1, fontsize=10, prop={'size': 11.5})
    #lgnd2.legendHandles[0]._sizes = [60]
    #lgnd2.legendHandles[1]._sizes = [60]
ra1_third = MILLI_sub1_visible['ra_deg']
dec1_third = MILLI_sub1_visible['dec_deg']
ra2_third = MILLI_sub2_visible['ra_deg']
dec2_third = MILLI_sub2_visible['dec_deg']
plot_objects_on_sky2(ra1_first, dec1_first,
                    ra2_first, dec2_first,
                    ra1_second, dec1_second,
                    ra2_second, dec2_second,
                    ra1_third, dec1_third,
                    ra2_third, dec2_third)<jupyter_output><empty_output><jupyter_text>Combining Dataframes<jupyter_code>frames1 = [MILLI_sub1_visible, MILLI_sub2_visible]
MILLI_all_filtered = pd.concat(frames1)
MILLI_all_filtered.head()
ra1 = MILLI_all_filtered['ra_deg']
dec1 = MILLI_all_filtered['dec_deg'] 
ra2 = goodObj['ra']
dec2 = goodObj['decl']
c = SkyCoord(ra=ra1*u.degree, dec=dec1*u.degree)
catalog = SkyCoord(ra=ra2*u.degree, dec=dec2*u.degree)
idx, d2d, d3d = c.match_to_catalog_sky(catalog)
plt.figure(figsize = (10,10))
plt.scatter(idx,d2d)
x = [0,50000]
y = [0.02,0.02]
plt.plot(x,y, color = 'red')<jupyter_output><empty_output><jupyter_text>Choose cutoff for distance and dia obj with 5 or more sources<jupyter_code>goodSrc.head()
ccdVisitIds = goodSrc['ccdVisitId']
dataIds = [{'visit': int(str(dataId)[0:6]), 'ccd': int(str(dataId)[6:])} for dataId in ccdVisitIds]
idx
MILLI_all_filtered['idx'] = idx
MILLI_all_filtered['d2d_deg'] = d2d
MILLI_all_filtered.head()
a = MILLI_all_filtered.query('d2d_deg < 0.05')['ra'].count()
b = MILLI_all_filtered.query('d2d_deg < 0.02')['ra'].count()
c = MILLI_all_filtered.query('d2d_deg < 0.01')['ra'].count()
d = MILLI_all_filtered.query('d2d_deg < 0.001')['ra'].count()
e = MILLI_all_filtered.query('d2d_deg < 0.0001')['ra'].count()
f = MILLI_all_filtered.query('d2d_deg < 0.00001')['ra'].count()

print(f'The total number of visible MILLIQUAS QSOs is {len(MILLI_all_filtered)}.')
print(f'QSOs with distance < 0.05 deg = {a}')
print(f'QSOs with distance < 0.02 deg = {b}')
print(f'QSOs with distance < 0.01 deg = {c}')
print(f'QSOs with distance < 0.001 deg = {d}')
print(f'QSOs with distance < 0.0001 deg = {e}')
print(f'QSOs with distance < 0.00001 deg = {f}')
ra2 = MILLI_all_close['ra_deg']
dec2 = MILLI_all_close['dec_deg'] 
ra3 = goodObj['ra']
dec3 = goodObj['decl']
c2 = SkyCoord(ra=ra2*u.degree, dec=dec2*u.degree)
catalog2 = SkyCoord(ra=ra3*u.degree, dec=dec3*u.degree)
idx2, d2d2, d3d2 = c2.match_to_catalog_sky(catalog2)
goodObj.head()
idx2
diaObjectIdQSO = goodObj['diaObjectId'].iloc[idx] #Hopefully these are the right ones?
MILLI_all_filtered['diaObjectId'] = diaObjectIdQSO.values

MILLI_all_filtered
MILLI_all_filtered.query('d2d_deg < 0.00001')<jupyter_output><empty_output><jupyter_text>Light Curves<jupyter_code>def plotLSSTLightcurve(diaObjectId, dbPath='/epyc/users/mrawls/hits2015/rerun/cw_2020_02/association.db'):
    '''Plot a light curve for a DIA (Difference Image Analysis) Object
    from an LSST APDB (Alert Production database).
    
    Parameters
    ----------
    diaObjectId : Identifier for an LSST DIA object
    
    dbPath : Path on disk to an APDB we can load DIA Objects or DIA Sources from
        often the database is named `association.db`

    '''
        
    # Load all sources for this object only
    connection = sqlite3.connect(dbPath)
    forcedSrcTable = pd.read_sql_query(f'select diaObjectId, midPointTai, psFlux, psFluxErr, totFlux, totFluxErr \
                                       from DiaForcedSource where diaObjectId = {diaObjectId};', connection)
    srcTable = pd.read_sql_query(f'select diaObjectId, midPointTai, ccdVisitId, psFlux, psFluxErr, totFlux, totFluxErr \
                                 from DiaSource where diaObjectId = {diaObjectId};', connection)

    # Convert flux to magnitude
    #mags = (np.array(forcedSrcTable['totFlux'])*u.nJy).to(u.ABmag)
    #magErrs = 1.0857 / forcedSrcTable['totFluxErr'] / forcedSrcTable['totFlux']
    #forcedSrcTable['mags'] = mags
    #forcedSrcTable['magErrs'] = magErrs
    
    #mags = (np.array(srcTable['totFlux'])*u.nJy).to(u.ABmag)
    #magErrs = 1.0857 / srcTable['totFluxErr'] / srcTable['totFlux']
    #srcTable['mags'] = mags
    #srcTable['magErrs'] = magErrs
    
    plt.figure(figsize=(10,4))
    
  
  
    # Plot the LSST light curve (forced PSF photometry on calexp at DIASource position)
    plt.errorbar(srcTable['midPointTai'], srcTable['totFlux'], 
                 yerr=srcTable['totFluxErr'], ls=':', marker='o', label='LSST Source', color='C0')
    
    # Plot the LSST light curve (forced PSF photometry on calexp at DIAObject position)
    #plt.errorbar(forcedSrcTable['midPointTai'], forcedSrcTable['totFlux'], 
     #            yerr=forcedSrcTable['totFluxErr'], ls=':', marker='o', label='LSST', color='C0')
    
    # Plot the LSST light curve directly from source catalogs
    #if plotSourceCatLightCurve:
       # objTable = pd.read_sql_query('select diaObjectId, ra, decl, nDiaSources \
        #                      from DiaObject where validityEnd is NULL;', connection)
       # sourceCatTimes, sourceCatMags, sourceCatMagErrs = getSourceCatLightCurve(obj, objTable, srcTable, butler)
       # plt.errorbar(sourceCatTimes, sourceCatMags, yerr=sourceCatMagErrs, ls=':', marker='o', label='LSST SourceCat', color='C2')
    
    

    # Add some plot labels, and reverse the y-axis so bright stuff is at the top
    plt.xlabel('Time (MJD)')
    plt.ylabel('Instantaneous Flux (nJy)')
    #plt.legend(frameon=False, prop={'size': 16})
for obj in MILLI_all_filtered['diaObjectId'].loc[MILLI_all_filtered['d2d_deg'] < 0.00001]:
     plotLSSTLightcurve(obj)
MILLI_all_filtered.to_csv('/epyc/users/tw94/lsst-hits-bump/mrawls/MILLI_all_filtered.csv', index = False)
v = 1
d = 1
x = np.linspace(0,100)
y = v*(x/d)**(4/3)
y2 = v*(x/d)
plt.plot(x,y)
plt.plot(x,y2,linestyle = '--')
plt.xlim(0,100)
plt.ylim(0,500)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15612549800796813, 'fraction_numerical': 0.04282868525896414, 'mean_word_length': 3.825390468562275, 'pattern_counts': {'":': 18, '<': 43, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '382853', 'license': 'MIT', 'n_tokens_mistral': 4985, 'n_tokens_neox': 4678, 'n_words': 866, 'path': 'ASTR_499_NB4_QSOobjectplots.ipynb', 'repo_name': 'mrawls/qso-hitsvar', 'size': '12048'}