<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Compare photometry in the new  Stripe82 catalog
to Gaia DR2 photometry and derive corrections for
gray systematics using Gmag photometry
input
Step1:  <a id='dataReading'></a>
Define paths and catalogs
Step2:  Simple positional match using ra/dec
Step3:  apply standard cuts as in old catalog
Step4:  now match to Gaia DR2...
Step5:  Select good matches and compare both catalogs to Gaia DR2
Step6:  Final Figures for the Paper
with Karun's 2D histogram implementation
Step7:  Final Gmag-based recalibration
Recalibrate R.A. residuals
Step8:  Recalibrate Dec residuals
Step9:  Now save correction arrays, then apply to original file, and then test
Step10:  paper plot showing the jump in Gaia Gmag
Step11:  TEMP code for color corrections to go from 3.1 to 3.2 and 3.3
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
from astropy.table import Table
from astropy.coordinates import SkyCoord
from astropy import units as u
from astropy.table import hstack
import matplotlib.pyplot as plt 
import numpy as np
from astroML.plotting import hist
# for astroML installation see https://www.astroml.org/user_guide/installation.html
## automatically reload any modules read below that might have changed (e.g. plots)
%load_ext autoreload
%autoreload 2
# importing ZI and KT tools: 
import ZItools as zit
import KTtools as ktt
Explanation: Compare photometry in the new  Stripe82 catalog
to Gaia DR2 photometry and derive corrections for
gray systematics using Gmag photometry
input: N2020_stripe82calibStars.dat
output: stripe82calibStars_v3.1.dat
files with RA/Dec corrections: ZPcorrectionsRA_v3.1_final.dat and ZPcorrectionsDec_v3.1_final.dat
makes paper plots:
GmagCorrection_RA_Hess.png  GmagCorrection_Dec_Hess.png
GmagCorrectionTest_Gmag_Hess.png
End of explanation
ZIdataDir = "/Users/ivezic/Work/Science/CalibrationV2/SDSS_SSC/Data"
# the original SDSS catalog from 2007
sdssOldCat = ZIdataDir + "/" + "stripe82calibStars_v2.6.dat"
# INPUT: Karun's new catalog from 2020
sdssNewCatIn = ZIdataDir + "/" + "N2020_stripe82calibStars.dat"
readFormat = 'csv'
# OUTPUT: with Gmag-based gray corrections 
sdssNewCatOut = ZIdataDir + "/" + "stripe82calibStars_v3.1.dat" 
# Gaia DR2 
GaiaDR2Cat = ZIdataDir + "/" + "Stripe82_GaiaDR2.dat"
# Gaia DR2 with BP and RP data
GaiaDR2CatBR = ZIdataDir + "/" + "Stripe82_GaiaDR2_BPRP.dat" 
# both new and old files use identical data structure
colnamesSDSS = ['calib_fla', 'ra', 'dec', 'raRMS', 'decRMS', 'nEpochs', 'AR_val', 
                'u_Nobs', 'u_mMed', 'u_mMean', 'u_mErr', 'u_rms_scatt', 'u_chi2',
                'g_Nobs', 'g_mMed', 'g_mMean', 'g_mErr', 'g_rms_scatt', 'g_chi2',
                'r_Nobs', 'r_mMed', 'r_mMean', 'r_mErr', 'r_rms_scatt', 'r_chi2',
                'i_Nobs', 'i_mMed', 'i_mMean', 'i_mErr', 'i_rms_scatt', 'i_chi2',
                'z_Nobs', 'z_mMed', 'z_mMean', 'z_mErr', 'z_rms_scatt', 'z_chi2']
%%time
# old
sdssOld = Table.read(sdssOldCat, format='ascii', names=colnamesSDSS) 
np.size(sdssOld)
%%time
# new 
sdssNew = Table.read(sdssNewCatIn, format=readFormat, names=colnamesSDSS)
np.size(sdssNew)
Explanation: <a id='dataReading'></a>
Define paths and catalogs
End of explanation
sdssOld_coords = SkyCoord(ra = sdssOld['ra']*u.degree, dec= sdssOld['dec']*u.degree) 
sdssNew_coords = SkyCoord(ra = sdssNew['ra']*u.degree, dec= sdssNew['dec']*u.degree) 
# this is matching sdssNew to sdssOld, so that indices are into sdssNew catalog
# makes sense in this case since the sdssOld catalog is (a little bit) bigger 
# than sdssNew (1006849 vs 1005470)
idx, d2d, d3d = sdssNew_coords.match_to_catalog_sky(sdssOld_coords)  
# object separation is an object with units, 
# I add that as a column so that one can 
# select based on separation to the nearest matching object
new_old = hstack([sdssNew, sdssOld[idx]], table_names = ['new', 'old'])
new_old['sep_2d_arcsec'] = d2d.arcsec
# good matches between the old and new catalogs
MAX_DISTANCE_ARCSEC = 0.5
sdss = new_old[(new_old['sep_2d_arcsec'] < MAX_DISTANCE_ARCSEC)]
print(np.size(sdss))
Explanation: Simple positional match using ra/dec
End of explanation
mOK3 = sdss[sdss['ra_new']<1]
mOK3 = zit.selectCatalog(sdss, mOK3)
print(996147/1006849)
print(993774/1006849)
print(991472/1006849)
Explanation: apply standard cuts as in old catalog:
End of explanation
colnamesGaia = ['ra', 'dec', 'nObs', 'Gmag', 'flux', 'fluxErr', 'pmra', 'pmdec']
colnamesGaia = colnamesGaia + ['BPmag', 'BPeI', 'RPmag', 'RPeI', 'BRef']
gaia = Table.read(GaiaDR2CatBR, format='ascii', names=colnamesGaia)
gaia['raG'] = gaia['ra']
gaia['decG'] = gaia['dec']    
gaia['GmagErr'] = gaia['fluxErr'] / gaia['flux'] 
gaia['BR'] = gaia['BPmag'] - gaia['RPmag'] 
gaia['GBP'] = gaia['Gmag'] - gaia['BPmag']
gaia['GRP'] = gaia['Gmag'] - gaia['RPmag']
sdss_coords = SkyCoord(ra = sdss['ra_old']*u.degree, dec= sdss['dec_old']*u.degree) 
gaia_coords = SkyCoord(ra = gaia['raG']*u.degree, dec= gaia['decG']*u.degree) 
# this is matching gaia to sdss, so that indices are into sdss catalog
# makes sense in this case since the sdss catalog is bigger than gaia
idxG, d2dG, d3dG = gaia_coords.match_to_catalog_sky(sdss_coords)  
# object separation is an object with units, 
# I add that as a column so that one can 
# select based on separation to the nearest matching object
gaia_sdss = hstack([gaia, sdss[idxG]], table_names = ['gaia', 'sdss'])
gaia_sdss['sepSG_2d_arcsec'] = d2dG.arcsec
### code for generating new quantities, such as dra, ddec, colors, differences in mags, etc
def derivedColumns(matches):
    matches['dra'] = (matches['ra_new']-matches['ra_old'])*3600
    matches['ddec'] = (matches['dec_new']-matches['dec_old'])*3600
    matches['ra'] = matches['ra_old']
    ra = matches['ra'] 
    matches['raW'] = np.where(ra > 180, ra-360, ra) 
    matches['dec'] = matches['dec_old']
    matches['u'] = matches['u_mMed_old']
    matches['g'] = matches['g_mMed_old']
    matches['r'] = matches['r_mMed_old']
    matches['i'] = matches['i_mMed_old']
    matches['z'] = matches['z_mMed_old']
    matches['ug'] = matches['u_mMed_old'] - matches['g_mMed_old']
    matches['gr'] = matches['g_mMed_old'] - matches['r_mMed_old']
    matches['ri'] = matches['r_mMed_old'] - matches['i_mMed_old']
    matches['gi'] = matches['g_mMed_old'] - matches['i_mMed_old']
    matches['du'] = matches['u_mMed_old'] - matches['u_mMed_new']
    matches['dg'] = matches['g_mMed_old'] - matches['g_mMed_new']
    matches['dr'] = matches['r_mMed_old'] - matches['r_mMed_new']
    matches['di'] = matches['i_mMed_old'] - matches['i_mMed_new']
    matches['dz'] = matches['z_mMed_old'] - matches['z_mMed_new']
    # Gaia 
    matches['draGold'] = -3600*(matches['ra_old'] - matches['raG']) 
    matches['draGnew'] = -3600*(matches['ra_new'] - matches['raG']) 
    matches['ddecGold'] = -3600*(matches['dec_old'] - matches['decG']) 
    matches['ddecGnew'] = -3600*(matches['dec_new'] - matches['decG']) 
    # photometric
    matches['gGr_old'] = matches['Gmag'] - matches['r_mMed_old']
    matches['gGr_new'] = matches['Gmag'] - matches['r_mMed_new']
    matches['gRPr_new'] = matches['RPmag'] - matches['r_mMed_new']
    return
derivedColumns(gaia_sdss) 
Explanation: now match to Gaia DR2...
End of explanation
# doGaiaAll(mOK)
def doGaiaGmagCorrection(d, Cstr, Gmax=20.0, yMax=0.03):
    # Cstr = 'gGr_old' or 'gGr_new'  
    gi = d['gi']
    Gr = d[Cstr]
    Gmag = d['Gmag']
    zit.qpBM(d, 'gi', -1, 4.5, Cstr, -2, 1.0, 56) 
    xBin, nPts, medianBin, sigGbin = zit.fitMedians(gi, Gr, -0.7, 4.0, 47, 0)  
    data = np.array([xBin, medianBin, sigGbin])
    Ndata = xBin.size
    ### HERE WE ARE FITTING 7-th ORDER POLYNOMIAL TO Gmag-rSDSS vs. g-i ###
    # get best-fit parameters  
    thetaCloc = zit.best_theta(data,7)
    # generate best fit lines on a fine grid 
    xfit = np.linspace(-1.1, 4.3, 1000)
    yfit = zit.polynomial_fit(thetaCloc, xfit) 
    ## added "Poly" because switched to piecewise linear interpolation below
    d['gGrFitPoly'] = zit.polynomial_fit(thetaCloc, gi)
    d['dgGrPoly'] = d[Cstr] - d['gGrFitPoly']  
    ### PIECEWISE LINEAR INTERPOLATION (AS FOR ALL OTHER COLORS AND SURVEYS)
    d['gGrFit'] = np.interp(gi, xBin, medianBin)
    d['dgGr'] = d[Cstr] - d['gGrFit']   
    
    # SELECT FOR RECALIBRATION wrt RA and Dec
    giMin = 0.4
    giMax = 3.0 
    Dc = d[(d['gi']>giMin)&(d['gi']<giMax)]
    print('N before and after color cut:', np.size(d), np.size(Dc))
    DcB = Dc[(Dc['Gmag']>14.5)&(Dc['Gmag']<Gmax)]
    DcB['GrResid'] = DcB['dgGr'] - np.median(DcB['dgGr'])
    zit.printStats(DcB['dgGr'])
    DcBok = DcB[np.abs(DcB['dgGr'])<0.1]
    print(np.size(DcB), np.size(DcBok))
    zit.qpBM(DcBok, 'Gmag', 14.5, Gmax, 'GrResid', -1*yMax, yMax, 56) 
    zit.qpBM(DcBok, 'dec', -1.3, 1.3, 'GrResid', -1*yMax, yMax, 126) 
    zit.qpBM(DcBok, 'raW', -51.5, 60, 'GrResid', -1*yMax, yMax, 112) 
    
    return thetaCloc, DcBok 
## first limit astrometric distance and 
## require at least 4 epochs as in the old catalog
MAX_DISTANCE_ARCSEC = 0.5
m1 = gaia_sdss[(gaia_sdss['sepSG_2d_arcsec'] < MAX_DISTANCE_ARCSEC)]
a1 = m1['g_Nobs_new']
a2 = m1['r_Nobs_new']
a3 = m1['i_Nobs_new']
mOK = m1[(a1>3)&(a2>3)&(a3>3)]
print(len(new_old))
print(len(m1))
print(len(mOK))
def plotAstro2Ddiagrams(d):
    ### plots 
    plotNameRoot = 'astroVSpm_RA_pm'
    plotName = plotNameRoot + '.png'    
    kw = {"Xstr":'pmra', "Xmin":-40, "Xmax":40, "Xlabel":'R.A. proper motion (mas/yr)', \
          "Ystr":'draGnew', "Ymin":-0.5, "Ymax":0.5, "Ylabel":'raw SDSS R.A. - Gaia R.A. (arcsec)', \
          "XminBin":-35, "XmaxBin":35, "nBin":70, \
          "plotName":plotName, "Nsigma":0, "offset":-0.1, "symbSize":0.05}
    kw["nBinX"] = 90
    kw["nBinY"] = 40
    kw["cmap"] = 'plasma'
    ktt.plotdelMagBW_KT(d, kw)
    print('made plot', plotName)
    # need to fit draGnew vs. pmra and correct for the mean trend, then plot vs. r mag
    pmra = d['pmra']
    draGnew = d['draGnew']
    xBin, nPts, medianBin, sigGbin = zit.fitMedians(pmra, draGnew, -60, 60, 120, 0)  
    ### PIECEWISE LINEAR INTERPOLATION 
    d['draGnewFit'] = np.interp(d['pmra'], xBin, medianBin)
    draCorr = d['draGnew'] - d['draGnewFit'] 
    draCorrOK = np.where(np.abs(draCorr) < 0.25, draCorr, 0)
    d['draGnewCorr'] = draCorrOK   
    plotNameRoot = 'astroVSpm_RA_r'
    plotName = plotNameRoot + '.png'  
    kw = {"Xstr":'r_mMed_new', "Xmin":14, "Xmax":21, "Xlabel":'SDSS r magnitude', \
          "Ystr":'draGnewCorr', "Ymin":-0.12, "Ymax":0.12, "Ylabel":'corr. SDSS R.A. - Gaia R.A. (arcsec)', \
          "XminBin":14, "XmaxBin":21, "nBin":30, \
          "plotName":plotName, "Nsigma":0, "offset":0.050, "symbSize":0.05}
    kw["nBinX"] = 30
    kw["nBinY"] = 24
    kw["cmap"] = 'plasma'
    ktt.plotdelMagBW_KT(d, kw)
    print('made plot', plotName)
    plotNameRoot = 'astroVSpm_Dec_pm'
    plotName = plotNameRoot + '.png'    
    kw = {"Xstr":'pmdec', "Xmin":-40, "Xmax":40, "Xlabel":'Dec. proper motion (mas/yr)', \
          "Ystr":'ddecGnew', "Ymin":-0.5, "Ymax":0.5, "Ylabel":'raw SDSS Dec. - Gaia Dec. (arcsec)', \
          "XminBin":-35, "XmaxBin":35, "nBin":70, \
          "plotName":plotName, "Nsigma":0, "offset":-0.1, "symbSize":0.05}
    kw["nBinX"] = 90
    kw["nBinY"] = 40
    kw["cmap"] = 'plasma'
    ktt.plotdelMagBW_KT(d, kw)
    print('made plot', plotName)
### produce astrometric plots showing correlation with proper motions
plotAstro2Ddiagrams(mOK)
# print(np.std(mOK['draGnew']), np.std(mOK['ddecGnew'])) 
#mOK
x = mOK['draGnewCorr']
xOK = x[np.abs(x)<0.25]
print(np.std(xOK), zit.sigG(xOK))
zit.qpBM(mOK, 'pmra', -50, 50, 'draGnew', -0.6, 0.6, 50) 
zit.qpBM(mOK, 'pmdec', -50, 50, 'ddecGnew', -0.6, 0.6, 50) 
theta, mOKc = doGaiaGmagCorrection(mOK, 'gGr_new')
thetaLoc = theta
## for zero point calibration, in addition to color cut in doGaiaAll, take 16 < G < 19.5 
mOKcB = mOKc[(mOKc['Gmag']>16)&(mOKc['Gmag']<19.5)]
mOKcB['GrResid'] = mOKcB['dgGr'] - np.median(mOKcB['dgGr'])
mOKcBok = mOKcB[np.abs(mOKcB['dgGr'])<0.1]
print(np.size(mOKc), np.size(mOKcB), np.size(mOKcBok))
print(np.std(mOKcBok['GrResid']), zit.sigG(mOKcBok['GrResid'])) 
zit.qpBM(mOKcBok, 'dec', -1.3, 1.3, 'GrResid', -0.03, 0.03, 260) 
zit.qpBM(mOKcBok, 'raW', -51.5, 60, 'GrResid', -0.03, 0.03, 112) 
Explanation: Select good matches and compare both catalogs to Gaia DR2
End of explanation
def plotGmag2Ddiagrams(d):
    ### plots 
    plotNameRoot = 'GrVSgi'
    plotName = plotNameRoot + '.png'    
    kw = {"Xstr":'gi', "Xmin":0.0, "Xmax":3.5, "Xlabel":'SDSS g-i', \
          "Ystr":'gGr_new', "Ymin":-1.25, "Ymax":0.25, "Ylabel":'Gaia Gmag - SDSS r', \
          "XminBin":-0.5, "XmaxBin":4.0, "nBin":90, \
          "plotName":plotName, "Nsigma":3, "offset":0.0, "symbSize":0.05}
    kw["nBinX"] = 90
    kw["nBinY"] = 40
    kw["cmap"] = 'plasma'
    ktt.plotdelMagBW_KT(d, kw)
    print('made plot', plotName)
def plotGmag2DdiagramsX(d, kw):
    # Gaia G
    print('-----------')
    print(' stats for SDSS r binning medians:')
    plotName = plotNameRoot + '_Gmag.png'    
    kwOC = {"Xstr":'Gmag', "Xmin":14.3, "Xmax":21.01, "Xlabel":'Gaia G (mag)', \
          "Ystr":kw['Ystr'], "Ymin":-0.06, "Ymax":0.06, "Ylabel":Ylabel, \
          "XminBin":14.5, "XmaxBin":21.0, "nBin":130, \
          "plotName":plotName, "Nsigma":3, "offset":0.01, "symbSize":kw['symbSize']}
    zit.plotdelMag(goodC, kwOC)
    plotName = plotNameRoot + '_Gmag_Hess.png'    
    kwOC["plotName"] = plotName
    kwOC["nBinX"] = 130
    kwOC["nBinY"] = 50
    kwOC["cmap"] = 'plasma'
    ktt.plotdelMagBW_KT(goodC, kwOC)
    print('made plot', plotName)
    print('------------------------------------------------------------------')
def plotGmagCorrections(d, kw):
    ### REDEFINE residuals to correspond to "SDSS-others", as other cases
    d['redef'] = -1*d[kw['Ystr']] 
    kw['Ystr'] = 'redef'
    goodC = d[np.abs(d['redef'])<0.1]
    
    ### plots 
    plotNameRoot = kw['plotNameRoot']
    # RA
    print(' stats for RA binning medians:')
    plotName = plotNameRoot + '_RA.png'
    Ylabel = 'residuals for (Gmag$_{SDSS}$ - Gmag$_{GaiaDR2}$)  '
    kwOC = {"Xstr":'raW', "Xmin":-52, "Xmax":60.5, "Xlabel":'R.A. (deg)', \
          "Ystr":kw['Ystr'], "Ymin":-0.07, "Ymax":0.07, "Ylabel":Ylabel, \
          "XminBin":-51.5, "XmaxBin":60, "nBin":112, \
          "plotName":plotName, "Nsigma":3, "offset":0.01, "symbSize":kw['symbSize']}
    zit.plotdelMag(goodC, kwOC)
    plotName = plotNameRoot + '_RA_Hess.png'    
    kwOC["plotName"] = plotName
    kwOC["nBinX"] = 112
    kwOC["nBinY"] = 50
    kwOC["cmap"] = 'plasma'
    ktt.plotdelMagBW_KT(goodC, kwOC)  
    print('made plot', plotName)
    # Dec
    print('-----------')
    print(' stats for Dec binning medians:')
    plotName = plotNameRoot + '_Dec.png'
    kwOC = {"Xstr":'dec', "Xmin":-1.3, "Xmax":1.3, "Xlabel":'Declination (deg)', \
          "Ystr":kw['Ystr'], "Ymin":-0.07, "Ymax":0.07, "Ylabel":Ylabel, \
          "XminBin":-1.266, "XmaxBin":1.264, "nBin":252, \
          "plotName":plotName, "Nsigma":3, "offset":0.01, "symbSize":kw['symbSize']}
    zit.plotdelMag(goodC, kwOC)
    plotName = plotNameRoot + '_Dec_Hess.png'    
    kwOC["plotName"] = plotName
    kwOC["nBinX"] = 252  
    kwOC["nBinY"] = 50
    kwOC["cmap"] = 'plasma'
    ktt.plotdelMagBW_KT(goodC, kwOC)
    print('made plot', plotName)   
    
    # Gaia G
    print('-----------')
    print(' stats for SDSS r binning medians:')
    plotName = plotNameRoot + '_Gmag.png'    
    kwOC = {"Xstr":'Gmag', "Xmin":14.3, "Xmax":21.01, "Xlabel":'Gaia G (mag)', \
          "Ystr":kw['Ystr'], "Ymin":-0.06, "Ymax":0.06, "Ylabel":Ylabel, \
          "XminBin":14.5, "XmaxBin":21.0, "nBin":130, \
          "plotName":plotName, "Nsigma":3, "offset":0.01, "symbSize":kw['symbSize']}
    zit.plotdelMag(goodC, kwOC)
    plotName = plotNameRoot + '_Gmag_Hess.png'    
    kwOC["plotName"] = plotName
    kwOC["nBinX"] = 130
    kwOC["nBinY"] = 50
    kwOC["cmap"] = 'plasma'
    ktt.plotdelMagBW_KT(goodC, kwOC)
    print('made plot', plotName)
    print('------------------------------------------------------------------')
 
mOK['GrResid'] = mOK['dgGr'] - np.median(mOK['dgGr']) + 0.006
mOKok = mOK[np.abs(mOK['dgGr'])<0.1]
print(np.size(mOK), np.size(mOKok))
keywords = {"Ystr":'GrResid', "plotNameRoot":'GmagCorrection', "symbSize":0.05}
plotGmagCorrections(mOKok, keywords) 
!cp GmagCorrection_Gmag_Hess.png GmagCorrectionTest_Gmag_Hess.png
mOKokX = mOKok[(mOKok['Gmag']>15)&(mOKok['Gmag']<15.5)]
print(np.median(mOKokX['GrResid']))
mOKokX = mOKok[(mOKok['Gmag']>16)&(mOKok['Gmag']<16.2)]
print(np.median(mOKokX['GrResid']))
keywords = {"Ystr":'GrResid', "plotNameRoot":'GmagCorrection', "symbSize":0.05}
plotGmagCorrections(mOKcBok, keywords) 
# for calibration: giMin = 0.4 & giMax = 3.0 
mOKB = mOK[(mOK['Gmag']>16)&(mOK['Gmag']<19.5)]
plotGmag2Ddiagrams(mOKB)
mOKB
Explanation: Final Figures for the Paper
with Karun's 2D histogram implementation
End of explanation
RAbin, RAnPts, RAmedianBin, RAsigGbin = zit.fitMedians(mOKcBok['raW'], mOKcBok['GrResid'], -51.5, 60.0, 112, 1)
Explanation: Final Gmag-based recalibration
Recalibrate R.A. residuals
End of explanation
decOK = mOKcBok['dec_new']
GrResid =  mOKcBok['GrResid']
fig,ax = plt.subplots(1,1,figsize=(8,6))
ax.scatter(decOK, GrResid, s=0.01, c='blue')
ax.set_xlim(-1.3,1.3)
ax.set_ylim(-0.06,0.06)
ax.set_ylim(-0.04,0.04)
ax.set_xlabel('Declination (deg)')
ax.set_ylabel('Gaia G - SDSS G')
xBin, nPts, medianBin, sigGbin = zit.fitMedians(decOK, GrResid, -1.266, 1.264, 252, 0)
ax.scatter(xBin, medianBin, s=30.0, c='black', alpha=0.9)
ax.scatter(xBin, medianBin, s=15.0, c='yellow', alpha=0.5)
TwoSigP = medianBin + 2*sigGbin
TwoSigM = medianBin - 2*sigGbin 
ax.plot(xBin, TwoSigP, c='yellow')
ax.plot(xBin, TwoSigM, c='yellow')
xL = np.linspace(-100,100)
ax.plot(xL, 0*xL+0.00, c='yellow')
ax.plot(xL, 0*xL+0.01, c='red')
ax.plot(xL, 0*xL-0.01, c='red')
dCleft = -1.3
ax.plot(0*xL+dCleft, xL, c='red')
alltheta = []
for i in range(0,12):
    decCol = -1.2655 + (i+1)*0.2109
    ax.plot(0*xL+decCol, xL, c='red')
    xR = xBin[(xBin>dCleft)&(xBin<decCol)]
    yR = medianBin[(xBin>dCleft)&(xBin<decCol)]
    dyR = sigGbin[(xBin>dCleft)&(xBin<decCol)]
    data = np.array([xR, yR, dyR])
    theta2 = zit.best_theta(data,5)
    alltheta.append(theta2)
    yfit = zit.polynomial_fit(theta2, xR)
    ax.plot(xR, yfit, c='cyan', lw=2)
    dCleft = decCol
    rrr = yR - yfit
    # print(i, np.median(rrr), np.std(rrr))  # 2 milimag scatter 
    # print(i, theta2)
plt.savefig('GmagDecCorrections.png')
# let's now correct all mags with this correction
thetaRecalib = alltheta
decLeft = -1.3
for i in range(0,12):
    decRight = -1.2655 + (i+1)*0.2109
    decArr = np.linspace(decLeft, decRight, 100)
    thetaBin = thetaRecalib[i] 
    ZPfit = zit.polynomial_fit(thetaBin, decArr)
    if (i==0):
        decCorrGrid = decArr
        ZPcorr = ZPfit
    else: 
        decCorrGrid = np.concatenate([decCorrGrid, decArr]) 
        ZPcorr = np.concatenate([ZPcorr, ZPfit])
    decLeft = decRight
mOKtest = mOK[mOK['r_Nobs_new']>3]
# Dec correction
decGrid2correct = mOKtest['dec_new']
ZPcorrectionsDec = np.interp(decGrid2correct, decCorrGrid, ZPcorr)
# RA correction 
raWGrid2correct = mOKtest['raW'] 
ZPcorrectionsRA = np.interp(raWGrid2correct, RAbin, RAmedianBin)
print(np.std(ZPcorrectionsDec), np.std(ZPcorrectionsRA))
fig,ax = plt.subplots(1,1,figsize=(8,6))
ax.scatter(decGrid2correct, ZPcorrectionsDec, s=0.01, c='blue')
ax.plot(decCorrGrid, ZPcorr, c='red')
ax.set_xlim(-1.3,1.3)
ax.set_ylim(-0.02,0.02)
ax.set_xlabel('Declination (deg)')
ax.set_ylabel('Correction')
fig,ax = plt.subplots(1,1,figsize=(8,6))
ax.scatter(raWGrid2correct, ZPcorrectionsRA, s=0.01, c='blue')
ax.plot(RAbin, RAmedianBin, c='red')
ax.set_xlim(-52,61)
ax.set_ylim(-0.05,0.05)
ax.set_xlabel('RA (deg)')
ax.set_ylabel('Correction') 
np.min(ZPcorrectionsRA)
mOKtest['ZPcorrectionsRA'] = ZPcorrectionsRA
mOKtest['ZPcorrectionsDec'] = ZPcorrectionsDec
mOKtest['r_mMed_new'] = mOKtest['r_mMed_new'] + mOKtest['ZPcorrectionsRA'] + mOKtest['ZPcorrectionsDec']
mOKtest['gGr_new'] = mOKtest['Gmag'] - mOKtest['r_mMed_new']
mOKtest['gGrFit'] = zit.polynomial_fit(thetaCloc, mOKtest['gi'])
mOKtest['dgGr'] = mOKtest['gGr_new'] - mOKtest['gGrFit']
d = mOKtest
gi = d['gi']
Gr = d['gGr_new']
Gmag = d['Gmag']
zit.qpBM(d, 'gi', -1, 4.5, 'gGr_new', -2, 1.0, 56) 
thetaCtest, DcBokTest_new = doGaiaGmagCorrection(mOKtest, 'gGr_new')
keywords = {"Ystr":'gGr_new', "plotNameRoot":'GmagCorrectionTest', "symbSize":0.05}
mOKtest2 = mOKtest[(mOKtest['gi']>0.4)&(mOKtest['gi']<3.0)]
x = mOKtest2[(mOKtest2['Gmag']>14.5)&(mOKtest2['Gmag']<15.5)]
mOKtest2['gGr_new'] = mOKtest2['gGr_new'] - np.median(x['gGr_new'])  
plotGmagCorrections(mOKtest2, keywords) 
Explanation: Recalibrate Dec residuals
End of explanation
# final refers to the July 2020 analysis, before the paper submission
#np.savetxt('ZPcorrectionsRA_v3.1_final.dat', (RAbin, RAmedianBin)) 
#np.savetxt('ZPcorrectionsDec_v3.1_final.dat', (decCorrGrid, ZPcorr))
sdssOut = sdss[sdss['ra_new']<1]
sdssOut = zit.selectCatalog(sdss, sdssOut)
sdssOut.sort('calib_fla_new') 
# read back gray zero point recalibration files 
zpRAgrid, zpRA = np.loadtxt('ZPcorrectionsRA_v3.1_final.dat')  
zpDecgrid, zpDec = np.loadtxt('ZPcorrectionsDec_v3.1_final.dat')  
sdssOut
# Dec correction
decGrid2correct = sdssOut['dec_new']
ZPcorrectionsDec = np.interp(decGrid2correct, zpDecgrid, zpDec)
# RA correction 
ra = sdssOut['ra_new'] 
raWGrid2correct = np.where(ra > 180, ra-360, ra) 
ZPcorrectionsRA = np.interp(raWGrid2correct, zpRAgrid, zpRA)
print('gray std RA/Dec:', np.std(ZPcorrectionsRA), np.std(ZPcorrectionsDec)) 
for b in ('u', 'g', 'r', 'i', 'z'):
    for mtype in ('_mMed_new', '_mMean_new'):
        mstr = b + mtype
        # applying here gray corrections  
        sdssOut[mstr] = sdssOut[mstr] + ZPcorrectionsRA + ZPcorrectionsDec 
SSCindexRoot = 'CALIBSTARS_'
outFile = ZIdataDir + "/" + "stripe82calibStars_v3.1_noheader_final.dat"
newSSC = open(outFile,'w')
df = sdssOut
Ngood = 0
for i in range(0, np.size(df)):
    Ngood += 1
    NoldCat = df['calib_fla_new'][i]
    strNo = f'{Ngood:07}'
    SSCindex = SSCindexRoot + strNo 
    SSCrow = zit.getSSCentry(df, i)
    zit.SSCentryToOutFileRow(SSCrow, SSCindex, newSSC) 
newSSC.close()
print(Ngood, 'rows in file', outFile)
Explanation: Now save correction arrays, then apply to original file, and then test
End of explanation
np.size(zpDec)
Explanation: paper plot showing the jump in Gaia Gmag
End of explanation
### need to figure out where were ZPcorrections2_rz_Dec.dat etc produced ... 
## color corrections 
for mtype in ('_mMed', '_mMean'):
    ## u band from u-r color
    color = 'ur'
    zpcFilename = 'ZPcorrections_' + color + '_RA.dat'
    zpcRAgrid, zpcRA = np.loadtxt(zpcFilename)   
    zpcFilename = 'ZPcorrections_' + color + '_Dec.dat'
    zpcDecgrid, zpcDec = np.loadtxt(zpcFilename)   
    ZPcorrectionsRA = np.interp(raWGrid2correct, zpcRAgrid, zpcRA) 
    ZPcorrectionsDec = np.interp(decGrid2correct, zpcDecgrid, zpcDec) 
    print('u-r std RA/Dec:', np.std(ZPcorrectionsRA), np.std(ZPcorrectionsDec)) 
    mstr = 'u' + mtype
    sdssOut[mstr] = sdssOut[mstr] - ZPcorrectionsRA - ZPcorrectionsDec 
    ## g band from g-r color
    color = 'gr'
    zpcFilename = 'ZPcorrections_' + color + '_RA.dat'
    zpcRAgrid, zpcRA = np.loadtxt(zpcFilename)   
    zpcFilename = 'ZPcorrections_' + color + '_Dec.dat'
    zpcDecgrid, zpcDec = np.loadtxt(zpcFilename)   
    ZPcorrectionsRA = np.interp(raWGrid2correct, zpcRAgrid, zpcRA) 
    ZPcorrectionsDec = np.interp(decGrid2correct, zpcDecgrid, zpcDec) 
    print('g-r std RA/Dec:', np.std(ZPcorrectionsRA), np.std(ZPcorrectionsDec)) 
    mstr = 'g' + mtype
    sdssOut[mstr] = sdssOut[mstr] - ZPcorrectionsRA - ZPcorrectionsDec 
    ## i band from r-i color
    color = 'ri'
    zpcFilename = 'ZPcorrections_' + color + '_RA.dat'
    zpcRAgrid, zpcRA = np.loadtxt(zpcFilename)   
    zpcFilename = 'ZPcorrections_' + color + '_Dec.dat'
    zpcDecgrid, zpcDec = np.loadtxt(zpcFilename)   
    ZPcorrectionsRA = np.interp(raWGrid2correct, zpcRAgrid, zpcRA) 
    ZPcorrectionsDec = np.interp(decGrid2correct, zpcDecgrid, zpcDec) 
    mstr = 'i' + mtype
    print('r-i std RA/Dec:', np.std(ZPcorrectionsRA), np.std(ZPcorrectionsDec)) 
    sdssOut[mstr] = sdssOut[mstr] + ZPcorrectionsRA + ZPcorrectionsDec 
    ## i band from r-z color
    color = 'rz'
    zpcFilename = 'ZPcorrections_' + color + '_RA.dat'
    zpcRAgrid, zpcRA = np.loadtxt(zpcFilename)   
    zpcFilename = 'ZPcorrections_' + color + '_Dec.dat'
    zpcDecgrid, zpcDec = np.loadtxt(zpcFilename)   
    ZPcorrectionsRA = np.interp(raWGrid2correct, zpcRAgrid, zpcRA) 
    ZPcorrectionsDec = np.interp(decGrid2correct, zpcDecgrid, zpcDec) 
    mstr = 'z' + mtype
    print('r-z std RA/Dec:', np.std(ZPcorrectionsRA), np.std(ZPcorrectionsDec)) 
    sdssOut[mstr] = sdssOut[mstr] + ZPcorrectionsRA + ZPcorrectionsDec 
Explanation: TEMP code for color corrections to go from 3.1 to 3.2 and 3.3
End of explanation
<END_TASK>