content: import sys, os, re, yaml, subprocess, shlex, FITS_tools
  import pandas as pd
  import numpy as np
  
  import pickle
  
  import matplotlib
  import matplotlib.pyplot as plt
  from matplotlib import image
  import matplotlib.cm as cm
  import matplotlib.image as mpimg
  
  from scipy.optimize import curve_fit
  import scipy.integrate as integrate
  from scipy import interpolate
  from scipy.interpolate import interp1d
  import scipy.stats
  
  from astropy.io import fits
  from astropy.table import Table, join
  from astropy import units as u
  from astropy.coordinates import SkyCoord
  from astropy.cosmology import FlatLambdaCDM
  cosmo = FlatLambdaCDM(H0=71, Om0=0.3, Tcmb0=2.725)
  
  import extra_program as ex
  from PIL import Image as Image_PIL
  import ebvpy  #Galactic Reddening
  
  """
  Example: 
  python pisco_pipeline/pisco_photometry_all_2019.py PKS1353 psf allslr 2mass
  python pisco_pipeline/pisco_photometry_all_2019.py PKS1353 psf allslr no2mass
  python pisco_pipeline/pisco_photometry_all_2019.py PKS1353 psf noslr no2mass
  python pisco_pipeline/pisco_photometry_all_2019.py PKS1353 auto noslr no2mass
  
  field: name of the fields
  mode: psf, auto, aper, hybrid, model
  allslr:
    - allslr: run everything including photometry_v4, cut_frame, SLR
    - slr: run just SLR and update the color
    - noslr: don't run slr, just update the color with different modes
  2mass
    - 2mass: run SLR with 2MASS to match
    - no2mass: run SLR without 2MASS
  """
  
  ###--------------------------------------------------------------------------###
  def find_seeing(field,band):
      df_see=pd.read_csv('/Users/taweewat/Documents/red_sequence/total_chips_field_seeing.csv',index_col=0)
      if field[0:5]=='CHIPS':
          seeing = df_see[df_see.chips==field]['seeing_q25_%s'%band].values[0] #_%s'%band
          return seeing
      elif (field[0:5]=='Field')|(field[0:3]=='PKS')|(field[0:4]=='SDSS'):
          seeing = df_see[df_see.name==field]['seeing_q25_%s'%band].values[0] #_%s'%band
          return seeing
  
  def find_seeing_new(dir,field):
      myReg3=re.compile(r'(CHIPS)[^\_]*\_[^\_]*')
      seeing = float(fits.open(list_file_name(dir,myReg3.search(field).group())[0])[0].header['FWHM1'])
      return seeing
  
  def find_seeing_fits(field):
      home='/Users/taweewat/Documents/pisco_code/'
      dirs=['ut170103/','ut170104/','ut170619/','ut170621/','ut170624/','ut171208/',\
      'ut171209/','ut171212/','ut190412/','ut190413/']
      myReg=re.compile(r'(%s_A).*'%field)
      for di in dirs:
          dir=home+di
          for text in os.listdir(dir):
              if myReg.search(text) != None:
                  seeing=float(fits.open(dir+myReg.search(text).group())[0].header['FWHM1'])
      return seeing
  def read_param():
      with open("pisco_pipeline/params.yaml", 'r') as stream:
          try:
              param=yaml.load(stream, Loader=yaml.FullLoader)
              return param
          except yaml.YAMLError as exc:
              print(exc)
  def read_param_izp(mode):
      if mode=='psf':
          mode_izp=''
      elif mode=='model':
          mode_izp='' #'_model'
      else:
          mode_izp=''
      # print "/Users/taweewat/Documents/pisco_code/pisco_pipeline/params_izeropoint%s.yaml" % mode_izp
      with open("/Users/taweewat/Documents/pisco_code/pisco_pipeline/params_izeropoint%s.yaml"%mode_izp, 'r') as stream:
          try:
              param=yaml.load(stream, Loader=yaml.FullLoader)
              return param
          except yaml.YAMLError as exc:
              print(exc)
  
  def star_galaxy_bleem(field):
      sg_dir = 'star_galaxy'
      if not os.path.exists(sg_dir):
          os.makedirs(sg_dir)
  
      param=read_param()
      # seeing=find_seeing(field,'i')
      # seeing=find_seeing_fits(field)
      seeing = 1.0
      # seeing=1.5
      # seeing=0.95
      minarea=1.7
      data, header = fits.getdata('final/coadd_c%s_i.fits'%field, header=True)
      data2=data**2
      # pxscale=0.11
      pxscale=0.22
      
      fits.writeto('final/coadd_c%s_sq_i.fits'%field, data2, header=header, overwrite=True)
      cmd='sex final/coadd_c%s_i.fits -c pisco_pipeline/config.sex -PARAMETERS_NAME pisco_pipeline/%s -CATALOG_NAME %s -CATALOG_TYPE FITS_1.0 -SEEING_FWHM %s -SATUR_LEVEL %s -PHOT_APERTURES 15 -PIXEL_SCALE %s -DETECT_MINAREA %s -CHECKIMAGE_NAME checki.fits,segmenti.fits'%\
      (field,'sex.param',sg_dir+'/%s_catalog.fits'%(field),str(seeing),str(param['satur_level_i_psf']),str(pxscale),str(1.1/minarea*np.pi*(seeing/pxscale)**2)); print cmd
      subprocess.check_call(shlex.split(cmd))
      cmd='sex final/coadd_c%s_i.fits,final/coadd_c%s_sq_i.fits -c pisco_pipeline/config.sex -PARAMETERS_NAME pisco_pipeline/%s -CATALOG_NAME %s -CATALOG_TYPE FITS_1.0 -SEEING_FWHM %s -SATUR_LEVEL %s -PHOT_APERTURES 15 -PIXEL_SCALE %s -DETECT_MINAREA %s'%\
      (field,field,'sex.param',sg_dir+'/%s_sq_catalog.fits'%(field),str(seeing),str(param['satur_level_i_sq_psf']),str(pxscale),str(1.1/minarea*np.pi*(seeing/pxscale)**2)); print cmd
      subprocess.check_call(shlex.split(cmd))
  
  def pisco_photometry_v4(field):
      def aperature_proj(field,band):
          param=read_param()
          seeing=find_seeing(field,band)
          # seeing=find_seeing_fits(field)
          # seeing = 1.1
          # seeing=1.5
  
          slrdir = 'slr_output'
          to_be_projected = 'final/coadd_c%s_%s.fits'%(field,band)
          reference_fits  = 'final/coadd_c%s_i.fits'%field
          im1,im2, header = FITS_tools.match_fits(to_be_projected,reference_fits,return_header=True)
          outname = 'final/proj_coadd_c%s_%s.fits'%(field,band)
          print 'projecting from %s band to i band the fits file '%band + outname
          fits.writeto(outname, im1, header, overwrite=True)
  
          minarea=1.7 #1.7
          pxscale=0.22
          # pxscale=0.11
          cmd='sex final/coadd_c%s_%s.fits -c pisco_pipeline/config.sex -PARAMETERS_NAME pisco_pipeline/%s -CATALOG_NAME %s -SEEING_FWHM %s -SATUR_LEVEL %s -PHOT_APERTURES 23 -PIXEL_SCALE %s -DETECT_MINAREA %s -CHECKIMAGE_NAME check_psf_%s.fits,segment_psf_%s.fits'%\
          (field,band,'sex_psf.param','psfex_output/psf_%s_%s.fits'%(field,band),str(seeing),str(param['satur_level_%s_psf'%band]),str(pxscale),str(1.1/minarea*np.pi*(seeing/pxscale)**2), band, band) 
          print cmd
          subprocess.check_call(shlex.split(cmd))
  
          Tf=Table(fits.open('psfex_output/psf_%s_%s.fits'%(field,band))[2].data)
          # Tfcut = Tf[(Tf['CLASS_STAR'] > 0.97) & (Tf['FLAGS'] == 0)].copy()  #0.97 Field292
          if len(Tf[(Tf['CLASS_STAR'] > 0.95) & (Tf['FLAGS'] < 5)]) > 0:
              Tfcut = Tf[(Tf['CLASS_STAR'] > 0.95) & (Tf['FLAGS'] < 5)].copy()
          else: 
              Tfcut = Tf[(Tf['CLASS_STAR'] > 0.9) & (Tf['FLAGS'] < 5)].copy()
          # Tfcut = Tf[(Tf['CLASS_STAR'] > 0.9) & (Tf['FLAGS'] < 5)].copy()  #0.97 Field292
          Tfcut_edge=Tfcut[(Tfcut['XWIN_IMAGE']<np.max(Tfcut['XWIN_IMAGE'])-60)&(Tfcut['XWIN_IMAGE']>np.min(Tfcut['XWIN_IMAGE'])+60)&\
                  (Tfcut['YWIN_IMAGE']<np.max(Tfcut['YWIN_IMAGE'])-60)&(Tfcut['YWIN_IMAGE']>np.min(Tfcut['YWIN_IMAGE'])+60)].copy()
          Tfcut_more=Tfcut_edge[(np.abs(Tfcut_edge['FLUX_RADIUS']-np.mean(Tfcut_edge['FLUX_RADIUS']))<2*np.std(Tfcut_edge['FLUX_RADIUS']))]
          Tfcut_more2=Tfcut_more[(np.abs(Tfcut_more['ELONGATION']-np.mean(Tfcut_more['ELONGATION']))<2*np.std(Tfcut_more['ELONGATION']))].copy()
          print "length of Tf: all: {}, CS>0.97: {}, edges: {}, flux_radius: {}, elong: {}".format(len(Tf), len(Tfcut), len(Tfcut_edge), len(Tfcut_more), len(Tfcut_more2))
          hdu = fits.open('psfex_output/psf_%s_%s.fits'%(field,band))
          hdu[2].data = hdu[2].data[Tfcut_more2['NUMBER']-1]
          # hdu[2].data = hdu[2].data[Tfcut['NUMBER']-1]
          hdu.writeto('psfex_output/psf_%s_%s.fits'%(field,band), overwrite=True)
  
          cmd='psfex %s -c pisco_pipeline/pisco.psfex' % ('psfex_output/psf_%s_%s.fits'%(field,band))
          print cmd
          subprocess.check_call(shlex.split(cmd))
  
          # minarea=3.0
          cmd='sex final/coadd_c%s_i.fits,final/proj_coadd_c%s_%s.fits -c pisco_pipeline/config.sex -PSF_NAME %s -PARAMETERS_NAME pisco_pipeline/%s -CATALOG_NAME %s -SEEING_FWHM %s -SATUR_LEVEL %s -PIXEL_SCALE %s -CATALOG_TYPE FITS_1.0 -PHOT_APERTURES 23 -DETECT_MINAREA %s -CHECKIMAGE_NAME check%s.fits,segment%s.fits'%\
              (field, field, band, 'psfex_output/psf_%s_%s.psf' % (field, band), 'sex_after_psf.param', '%s/a_psf_%s_%s.fits' % (slrdir, field, band),
               str(seeing), str(param['satur_level_%s_psf' % band]), str(pxscale), str(1.1 / minarea * np.pi * (seeing / pxscale)**2), band, band)
          print cmd
          subprocess.check_call(shlex.split(cmd))
  
          table=Table.read('%s/a_psf_%s_%s.fits'%(slrdir,field,band))
          for name in table.colnames[:]:
              table.rename_column(name, name + '_%s' % band)
          return table
      
      slrdir = 'slr_output'
      if not os.path.exists(slrdir):
          os.makedirs(slrdir)
  
      tableg=aperature_proj(field,'g')
      tablei=aperature_proj(field,'i')
      tabler=aperature_proj(field,'r')
      tablez=aperature_proj(field,'z')
  
      print 'len of all table', len(tableg), len(tablei), len(tabler), len(tablez)
  
      ci=SkyCoord(ra=np.array(tablei['ALPHA_J2000_i'])*u.degree, dec=np.array(tablei['DELTA_J2000_i'])*u.degree)# print len(ci)
      cg=SkyCoord(ra=np.array(tableg['ALPHA_J2000_g'])*u.degree, dec=np.array(tableg['DELTA_J2000_g'])*u.degree)# print len(cg)
      cr=SkyCoord(ra=np.array(tabler['ALPHA_J2000_r'])*u.degree, dec=np.array(tabler['DELTA_J2000_r'])*u.degree)# print len(cr)
      cz=SkyCoord(ra=np.array(tablez['ALPHA_J2000_z'])*u.degree, dec=np.array(tablez['DELTA_J2000_z'])*u.degree)# print len(cz)
  
      idxn, d2dn, d3dn=cg.match_to_catalog_sky(ci)
      # Table_I=tablei[idxn][['NUMBER_i','XWIN_IMAGE_i','YWIN_IMAGE_i','ALPHA_J2000_i','DELTA_J2000_i','MAG_APER_i','MAGERR_APER_i','MAG_AUTO_i','MAGERR_AUTO_i','MAG_HYBRID_i','MAGERR_HYBRID_i',\
      #               'CLASS_STAR_i','FLAGS_i','MAG_PSF_i','MAGERR_PSF_i','MAG_MODEL_i','MAGERR_MODEL_i','SPREAD_MODEL_i']]
      Table_I=tablei[idxn][['NUMBER_i','XWIN_IMAGE_i','YWIN_IMAGE_i','ALPHA_J2000_i','DELTA_J2000_i','MAG_APER_i','MAGERR_APER_i','MAG_AUTO_i','MAGERR_AUTO_i','MAG_SPHEROID_i','MAGERR_SPHEROID_i',\
                    'CLASS_STAR_i','FLAGS_i','MAG_PSF_i','MAGERR_PSF_i','MAG_MODEL_i','MAGERR_MODEL_i','SPREAD_MODEL_i','SPREADERR_MODEL_i','MAG_ISO_i','MAGERR_ISO_i']]
      Table_I.rename_column('ALPHA_J2000_i','ALPHA_J2000')
      Table_I.rename_column('DELTA_J2000_i','DELTA_J2000')
  
      idxn, d2dn, d3dn=cg.match_to_catalog_sky(cr)
      # Table_R=tabler[idxn][['NUMBER_r','ALPHA_J2000_r','DELTA_J2000_r','MAG_APER_r','MAGERR_APER_r','MAG_AUTO_r','MAGERR_AUTO_r','MAG_HYBRID_r','MAGERR_HYBRID_r',\
      #               'CLASS_STAR_r','FLAGS_r','MAG_PSF_r','MAGERR_PSF_r','MAG_MODEL_r','MAGERR_MODEL_r','SPREAD_MODEL_r']]
      Table_R=tabler[idxn][['NUMBER_r','ALPHA_J2000_r','DELTA_J2000_r','MAG_APER_r','MAGERR_APER_r','MAG_AUTO_r','MAGERR_AUTO_r','MAG_SPHEROID_r','MAGERR_SPHEROID_r',\
                    'CLASS_STAR_r','FLAGS_r','MAG_PSF_r','MAGERR_PSF_r','MAG_MODEL_r','MAGERR_MODEL_r','SPREAD_MODEL_r','SPREADERR_MODEL_r','MAG_ISO_r','MAGERR_ISO_r']]
      Table_R.rename_column('ALPHA_J2000_r','ALPHA_J2000')
      Table_R.rename_column('DELTA_J2000_r','DELTA_J2000')
  
      idxn, d2dn, d3dn=cg.match_to_catalog_sky(cz)
      # Table_Z=tablez[idxn][['NUMBER_z','ALPHA_J2000_z','DELTA_J2000_z','MAG_APER_z','MAGERR_APER_z','MAG_AUTO_z','MAGERR_AUTO_z','MAG_HYBRID_z','MAGERR_HYBRID_z',\
      #               'CLASS_STAR_z','FLAGS_z','MAG_PSF_z','MAGERR_PSF_z','MAG_MODEL_z','MAGERR_MODEL_z','SPREAD_MODEL_z']]
      Table_Z=tablez[idxn][['NUMBER_z','ALPHA_J2000_z','DELTA_J2000_z','MAG_APER_z','MAGERR_APER_z','MAG_AUTO_z','MAGERR_AUTO_z','MAG_SPHEROID_z','MAGERR_SPHEROID_z',\
                    'CLASS_STAR_z','FLAGS_z','MAG_PSF_z','MAGERR_PSF_z','MAG_MODEL_z','MAGERR_MODEL_z','SPREAD_MODEL_z','SPREADERR_MODEL_z','MAG_ISO_z','MAGERR_ISO_z']]
      Table_Z.rename_column('ALPHA_J2000_z','ALPHA_J2000')
      Table_Z.rename_column('DELTA_J2000_z','DELTA_J2000')
  
      # Table_G=tableg[['NUMBER_g','ALPHA_J2000_g','DELTA_J2000_g','MAG_APER_g','MAGERR_APER_g','MAG_AUTO_g','MAGERR_AUTO_g','MAG_HYBRID_g','MAGERR_HYBRID_g',\
      #               'CLASS_STAR_g','FLAGS_g','MAG_PSF_g','MAGERR_PSF_g','MAG_MODEL_g','MAGERR_MODEL_g','SPREAD_MODEL_g']]
      Table_G = tableg[['NUMBER_g', 'ALPHA_J2000_g', 'DELTA_J2000_g', 'MAG_APER_g', 'MAGERR_APER_g', 'MAG_AUTO_g', 'MAGERR_AUTO_g', 'MAG_SPHEROID_g', 'MAGERR_SPHEROID_g',
                    'CLASS_STAR_g','FLAGS_g','MAG_PSF_g','MAGERR_PSF_g','MAG_MODEL_g','MAGERR_MODEL_g','SPREAD_MODEL_g','SPREADERR_MODEL_g','MAG_ISO_g','MAGERR_ISO_g']]
      Table_G.rename_column('ALPHA_J2000_g','ALPHA_J2000')
      Table_G.rename_column('DELTA_J2000_g','DELTA_J2000')
  
      print 'len of all new table', len(Table_G), len(Table_I), len(Table_R), len(Table_Z)
  
      total=join(join(join(Table_I,Table_G,keys=['ALPHA_J2000','DELTA_J2000']),Table_R,keys=['ALPHA_J2000','DELTA_J2000']),\
           Table_Z,keys=['ALPHA_J2000','DELTA_J2000'])
  
      # total=join(join(join(mag_ii,mag_ig,keys='NUMBER'), mag_ir,keys='NUMBER'),\
      #            mag_iz,keys='NUMBER')
  
      # total2=total[['ALPHA_J2000','DELTA_J2000','NUMBER_i','NUMBER_r','NUMBER_g','XWIN_IMAGE_i','YWIN_IMAGE_i',\
      #               'MAG_APER_i','MAGERR_APER_i','MAG_APER_g','MAGERR_APER_g','MAG_APER_r',\
      #               'MAGERR_APER_r','MAG_APER_z','MAGERR_APER_z','MAG_AUTO_i','MAGERR_AUTO_i',\
      #               'MAG_AUTO_g','MAGERR_AUTO_g','MAG_AUTO_r','MAGERR_AUTO_r','MAG_AUTO_z',\
      #               'MAGERR_AUTO_z','MAG_HYBRID_i','MAGERR_HYBRID_i','MAG_HYBRID_g',\
      #               'MAGERR_HYBRID_g','MAG_HYBRID_r','MAGERR_HYBRID_r','MAG_HYBRID_z',\
      #               'MAGERR_HYBRID_z','CLASS_STAR_i','CLASS_STAR_g','CLASS_STAR_r',\
      #               'CLASS_STAR_z','FLAGS_g','FLAGS_r','FLAGS_i','FLAGS_z','MAG_PSF_g',\
      #               'MAG_PSF_r','MAG_PSF_i','MAG_PSF_z','MAGERR_PSF_g','MAGERR_PSF_r',\
      #               'MAGERR_PSF_i','MAGERR_PSF_z','MAG_MODEL_g','MAG_MODEL_r',\
      #               'MAG_MODEL_i','MAG_MODEL_z','MAGERR_MODEL_g','MAGERR_MODEL_r',\
      #               'MAGERR_MODEL_i','MAGERR_MODEL_z','SPREAD_MODEL_g','SPREAD_MODEL_r',\
      #               'SPREAD_MODEL_i','SPREAD_MODEL_z',]]
  
      total.write(os.path.join(slrdir, 'total0_psf_%s.csv' % field), overwrite=True)
  
      total2=total[['ALPHA_J2000','DELTA_J2000','NUMBER_i','NUMBER_r','NUMBER_g','XWIN_IMAGE_i','YWIN_IMAGE_i',\
                  'MAG_APER_i','MAGERR_APER_i','MAG_APER_g','MAGERR_APER_g','MAG_APER_r',\
                  'MAGERR_APER_r','MAG_APER_z','MAGERR_APER_z','MAG_AUTO_i','MAGERR_AUTO_i',\
                  'MAG_AUTO_g','MAGERR_AUTO_g','MAG_AUTO_r','MAGERR_AUTO_r','MAG_AUTO_z',\
                  'MAGERR_AUTO_z','MAG_ISO_g','MAGERR_ISO_g','MAG_ISO_r','MAGERR_ISO_r',\
                  'MAG_ISO_i','MAGERR_ISO_i','MAG_ISO_z','MAGERR_ISO_z',\
                  'MAG_SPHEROID_i','MAGERR_SPHEROID_i','MAG_SPHEROID_g',\
                  'MAGERR_SPHEROID_g','MAG_SPHEROID_r','MAGERR_SPHEROID_r','MAG_SPHEROID_z',\
                  'MAGERR_SPHEROID_z','CLASS_STAR_i','CLASS_STAR_g','CLASS_STAR_r',\
                  'CLASS_STAR_z','FLAGS_g','FLAGS_r','FLAGS_i','FLAGS_z','MAG_PSF_g',\
                  'MAG_PSF_r','MAG_PSF_i','MAG_PSF_z','MAGERR_PSF_g','MAGERR_PSF_r',\
                  'MAGERR_PSF_i','MAGERR_PSF_z','MAG_MODEL_g','MAG_MODEL_r',\
                  'MAG_MODEL_i','MAG_MODEL_z','MAGERR_MODEL_g','MAGERR_MODEL_r',\
                  'MAGERR_MODEL_i','MAGERR_MODEL_z','SPREAD_MODEL_g','SPREAD_MODEL_r',\
                  'SPREAD_MODEL_i','SPREAD_MODEL_z','SPREADERR_MODEL_g','SPREADERR_MODEL_r',\
                  'SPREADERR_MODEL_i','SPREADERR_MODEL_z']]
  
      total2.write(os.path.join(slrdir, 'total_psf_%s.csv' % field), overwrite=True)
      # total2.write(slrdir+'/all_psf_%s.fits' % field, overwrite=True)
  
  def pisco_cut_star(field,c_a,c_b,c_d,c_delta):
      seeing=find_seeing_fits(field)
      true_seeing=find_seeing(field,'i')
  
      df_i=Table(fits.open('/Users/taweewat/Documents/pisco_code/star_galaxy/%s_catalog.fits'%field)[1].data).to_pandas()
      df_isq=Table(fits.open('/Users/taweewat/Documents/pisco_code/star_galaxy/%s_sq_catalog.fits'%field)[1].data).to_pandas()
  
      #cut the object out so that it has the same number of object between the sq catalog list and the psf mag list.
      fname = "/Users/taweewat/Documents/pisco_code/slr_output/total_psf_%s.csv"%field
      df0 = pd.read_csv(fname)
      df0['NUMBER'] = np.arange(0, len(df0), 1).tolist()
      cf_i=SkyCoord(ra=np.array(df_i['ALPHA_J2000'])*u.degree, dec=np.array(df_i['DELTA_J2000'])*u.degree)
      cf_isq=SkyCoord(ra=np.array(df_isq['ALPHA_J2000'])*u.degree, dec=np.array(df_isq['DELTA_J2000'])*u.degree)
      cf0=SkyCoord(ra=np.array(df0['ALPHA_J2000'])*u.degree, dec=np.array(df0['DELTA_J2000'])*u.degree)
      df0.rename(columns={'ALPHA_J2000': 'ALPHA_J2000_i'}, inplace=True)
      df0.rename(columns={'DELTA_J2000': 'DELTA_J2000_i'}, inplace=True)
  
      idxn, d2dn, d3dn=cf0.match_to_catalog_sky(cf_i)
      df_i_cut0=df_i.loc[idxn].copy()
      df_i_cut0['NUMBER']=np.arange(0,len(df0),1).tolist()
      df_i_cut=pd.merge(df_i_cut0,df0,on='NUMBER')
  
      idxn, d2dn, d3dn=cf0.match_to_catalog_sky(cf_isq)
      df_isq_cut0=df_isq.loc[idxn].copy()
      df_isq_cut0['NUMBER']=np.arange(0,len(df0),1).tolist()
      df_isq_cut=pd.merge(df_isq_cut0,df0,on='NUMBER')
  
      fig,ax=plt.subplots(2,3,figsize=(15,10))
      df_i0=df_i_cut[(df_i_cut.MAG_APER<0)&(df_isq_cut.MAG_APER<0)]
      df_isq0=df_isq_cut[(df_i_cut.MAG_APER<0)&(df_isq_cut.MAG_APER<0)]# print len(df_i), len(df_isq)
  
      # c_d=-7.5
      df_i2=df_i0[(df_i0.CLASS_STAR>c_a) & (df_i0.MAG_APER<c_d)]# & (df_i0.MAG_APER>c_c)]
      df_isq2=df_isq0[(df_i0.CLASS_STAR>c_a) & (df_i0.MAG_APER<c_d)]# & (df_i0.MAG_APER>c_c)];# print len(df_i2), len(df_isq2)
  
      icut_per=np.percentile(df_i2.MAG_APER,35) #35
      df_i3=df_i2[df_i2.MAG_APER>icut_per]
      df_isq3=df_isq2[df_i2.MAG_APER>icut_per]
  
      fit=np.polyfit(df_i3.MAG_APER, df_i3.MAG_APER-df_isq3.MAG_APER, 1)
      f=np.poly1d(fit)
      ax[0,0].plot(df_i2.MAG_APER,f(df_i2.MAG_APER),'--')
  
      res=(df_i3.MAG_APER-df_isq3.MAG_APER)-f(df_i3.MAG_APER)
      aa=np.abs(res)<1.5*np.std(res)
      # outl=np.abs(res)>=1.5*np.std(res)
      fit=np.polyfit(df_i3.MAG_APER[aa], df_i3.MAG_APER[aa]-df_isq3.MAG_APER[aa], 1)
      f=np.poly1d(fit)
  
      ax[0,0].axvline(icut_per,color='blue',label='35th quantile')
      ax[0,0].errorbar(df_i2.MAG_APER,df_i2.MAG_APER-df_isq2.MAG_APER,yerr=np.sqrt(df_i2.MAGERR_APER**2+df_isq2.MAGERR_APER**2),fmt='o')
      ax[0,0].set_title('only for star')
      ax[0,0].plot(df_i2.MAG_APER,f(df_i2.MAG_APER),'--',label='no outlier')
      ax[0,0].set_ylabel('MAG_APER-MAG_APER_sq')
      ax[0,0].set_xlabel('MAG APER i')
  
      #--->  #0.1 default, 0.2
      c_c=df_i2[f(df_i2.MAG_APER)-(df_i2.MAG_APER-df_isq2.MAG_APER)<0.1]['MAG_APER'].values\
      [np.argmin(df_i2[f(df_i2.MAG_APER)-(df_i2.MAG_APER-df_isq2.MAG_APER)<0.1]['MAG_APER'].values)] #edit10/30 (previous 0.1)
      #--->
  
      ax[0,0].axvline(c_c,color='red',label='new upper cut')
      ax[0,0].legend(loc='best')
  
      # color_axis='CLASS_STAR'
      color_axis='SPREAD_MODEL_i'
  
  
      ax[0,1].scatter(df_i0.MAG_APER,df_i0.MAG_APER-df_isq0.MAG_APER,marker='.',c=df_i0[color_axis],vmin=0., vmax=0.005)
      ax[0,1].plot(df_i3.MAG_APER,df_i3.MAG_APER-df_isq3.MAG_APER,'x')
      ax[0,1].set_title('for all objects')
      ax[0,1].set_ylabel('MAG_APER-MAG_APER_sq')
      ax[0,1].set_xlabel('MAG APER i')
      ax[0,1].axvline(c_b,ls='--')
      ax[0,1].axvline(c_c,ls='--')
  
      delta=(df_i0.MAG_APER-df_isq0.MAG_APER) - f(df_i0.MAG_APER)
  
      ax[0,2].scatter(df_i0.MAG_APER,delta,marker='.',c=df_i0[color_axis],vmin=0., vmax=0.005)
      ax[0,2].axhline(0,ls='--')
      ax[0,2].axvline(c_c,ls='--')
      ax[0,2].axvline(c_b,ls='--')
      ax[0,2].set_ylabel('Delta')
      ax[0,2].set_xlabel('MAG APER i')
      ax[0,2].set_ylim(0.5,-1.2)
  
      df_i1=df_i0[(df_i0.MAG_APER>c_c)&(df_i0.MAG_APER<c_b)].copy()
      df_isq1=df_isq0[(df_i0.MAG_APER>c_c)&(df_i0.MAG_APER<c_b)].copy()
      delta1=(df_i1.MAG_APER-df_isq1.MAG_APER) - f(df_i1.MAG_APER)
      ax[1,0].scatter(df_i1.MAG_APER, delta1, marker='o', c=df_i1[color_axis],vmin=0., vmax=0.005)
      ax[1,0].axhline(0,ls='--')
      ax[1,0].axhline(c_delta, ls='--')
      ax[1,0].set_ylabel('Delta')
      ax[1,0].set_xlabel('MAG APER i')
      ax[1,0].set_ylim(0.5,-2)
  
      # deltag=delta1[delta1<c_delta] #galaxy  0.1, 0.2 (0.005), 0.5 ()
      deltas=delta1[(delta1>=c_delta)&(delta1<3.)] #star
  
      def gauss(x, *p):
          A, mu, sigma = p
          return A*np.exp(-(x-mu)**2/(2.*sigma**2))
      p0 = [1., 0., 0.1]
      # def gauss(x, *p):
      #     A, sigma = p
      #     return A*np.exp(-(x-0)**2/(2.*sigma**2))
      # p0 = [1., 0.1]
  
      #galaxy
      # hist, bin_edges = np.histogram(deltag,bins=np.arange(-1.2,0.5,0.02))
      hist, bin_edges = np.histogram(delta1,bins=np.arange(-1.2,0.5,0.02))
      bin_centres = (bin_edges[:-1] + bin_edges[1:])/2
      ax[1,1].plot(bin_centres, hist, label='galaxies',linestyle='steps')
  
      #stars
      hist, bin_edges = np.histogram(deltas,bins=np.arange(-1,0.5,0.02)) #(0 vs -1,0.5,0.02)
      # hist, bin_edges = np.histogram(delta1, bins=np.arange(c_delta, 0.5, 0.02))
      bin_centres = (bin_edges[:-1] + bin_edges[1:])/2
      coeff2, var_matrix = curve_fit(gauss, bin_centres, hist, p0=p0)
      ax[1,1].plot(bin_centres, hist, label='stars',linestyle='steps')
  
      # hist, bin_edges = np.histogram(delta1,bins=np.arange(-1.2,0.5,0.02)) #added for right gaussian fitting
      # bin_centres = (bin_edges[:-1] + bin_edges[1:])/2  # added for right gaussian fitting
      x=np.arange(-1.25,0.5,0.02)
      # hist_fit2 = gauss(x, *coeff2)
      hist_fit2 = gauss(x, *coeff2)
      hist_fit3 = gauss(x, *coeff2)/np.max(gauss(x, *coeff2)) #added for right gaussian fitting
      ax[1,1].plot(x, hist_fit2, label='stars_fit')
      ax[1,1].plot(x, hist_fit3, label='stars_fit_norm') #added for right gaussian fitting
      ax[1,1].axvline(x[hist_fit3>star_cut][0],c='tab:pink',label='cut:%.3f'%x[hist_fit3>star_cut][0])  #added for right gaussian fitting
      ax[1,1].legend(loc='best')
      ax[1,1].set_xlabel('Delta')
      ax[1,1].set_ylabel('Histogram')
  
      ax[0,2].axhline(x[hist_fit3>star_cut][0],c='tab:pink') #added for right gaussian fitting
      ax[1,0].axhline(x[hist_fit3>star_cut][0],c='tab:pink') #added for right gaussian fitting
      ax[1,2].axhline(star_cut, c='tab:red')  # added for right gaussian fitting
  
      maxi=np.max(gauss(delta,*coeff2))
      def prob_SG(delta,maxi,*coeff2):
          if delta>0.:
              return 0.
          elif delta<=0.:
              return 1. - (gauss(delta, *coeff2) / maxi)
  
      vprob_SG= np.vectorize(prob_SG)
      SG=1.-vprob_SG(delta1,maxi,*coeff2)
      df_i1.loc[:,'SG']=SG
      param_izp=read_param_izp('psf')
      mag0=param_izp['i_zp_day%i'%dir_dict[find_fits_dir(field)[-9:]]]
      axi = ax[1, 2].scatter(df_i1.MAG_APER + mag0, SG,
                             marker='.', c=df_i1[color_axis], vmin=0., vmax=0.005)
      ax[1,2].axvline(aper_cut, ls='--', c='tab:blue')
      ax[1,2].axhline(SG_upper, ls='--', c='tab:blue')
      ax[1,2].set_ylim(-0.02,1.02)
      ax[1,2].set_xlabel('MAG APER i')
      ax[1,2].set_ylabel('SG (probability to be a star)')
      plt.suptitle(field+' seeing vs true_seeing: '+str(seeing)+','+str(true_seeing))
      fig.colorbar(axi)
      plt.tight_layout(rect=[0, 0., 1, 0.98])
      plt.savefig('/Users/taweewat/Documents/red_sequence/pisco_color_plots/star_galaxy_sep_12_all%s.png' % field, dpi=120)
      plt.close(fig)
      return df_i_cut, df_i1
  
  def pisco_cut_frame(field):
      # df_i=Table(fits.open('/Users/taweewat/Documents/pisco_code/star_galaxy/'+
      #                      '%s_catalog.fits'%field)[1].data).to_pandas()
      """
      c_a: CLASS_STAR lower limit for stars used for the linear fit
      c_b, c_c: upper and lower limit for all objects selection
      c_c can be moved with the for loop to include more objects until the confusion limit
      c_d: Faintest magnitude for stars used for the linear fit
      c_delta: lower limit for Delta to consider stars before fitting the gaussian and find SG (Star/Galaxy) factor 
      """
      seeing=find_seeing_fits(field)
      true_seeing=find_seeing(field,'i')
  
      ##Using SPREAD_MODEL to seperate star/galaxies
      fname = "/Users/taweewat/Documents/pisco_code/slr_output/total_psf_%s.csv"%field
      df0 = pd.read_csv(fname)
      df0['NUMBER'] = np.arange(0, len(df0), 1).tolist()
      df0.rename(columns={'ALPHA_J2000': 'ALPHA_J2000_i'}, inplace=True)
      df0.rename(columns={'DELTA_J2000': 'DELTA_J2000_i'}, inplace=True)
  
      #EXTENDED_COADD: 0 star, 1 likely star, 2 mostly galaxies, 3 galaxies
      # df0['EXTENDED_COADD']=np.array(((df0['SPREAD_MODEL_i']+ 3*df0['SPREADERR_MODEL_i'])>0.005).values, dtype=int)+\
      # np.array(((df0['SPREAD_MODEL_i']+df0['SPREADERR_MODEL_i'])>0.003).values, dtype=int)+\
      # np.array(((df0['SPREAD_MODEL_i']-df0['SPREADERR_MODEL_i'])>0.003).values, dtype=int)
      # dff=df0[df0['EXTENDED_COADD']>1]
      # dff_star=df0[df0['EXTENDED_COADD']<2]
  
      dfi=df0[df0['MAG_AUTO_i']<-16]
      x=dfi['MAG_AUTO_i']
      y=dfi['SPREAD_MODEL_i']
      p_spread=np.poly1d(np.polyfit(x,y,1))
      xs=np.arange(np.min(df0['MAG_AUTO_i']),np.max(df0['MAG_AUTO_i']),0.01)
  
      df0['SPREAD_MODEL_i2']=df0['SPREAD_MODEL_i']-p_spread(df0['MAG_AUTO_i'])
  
  
      dff=df0[(df0['SPREAD_MODEL_i']>0.005)]
      # dff_star=df0[np.abs(df0['SPREAD_MODEL_i'])<0.004] #+5/3.*df0['SPREADERR_MODEL_i'] <0.002
      dff_star=df0[(df0['SPREAD_MODEL_i']<0.004)]#&(df0['FLAGS_i']<4)]
  
      fig=plt.figure(figsize=(4,4))
      plt.plot(df0['MAG_AUTO_i'],df0['SPREAD_MODEL_i'],'.',c='grey',alpha=0.1)
      plt.plot(dff['MAG_AUTO_i'],dff['SPREAD_MODEL_i'],'.',alpha=1,label='galaxies')
      plt.plot(dff_star['MAG_AUTO_i'],dff_star['SPREAD_MODEL_i'],'.',alpha=1,label='stars')
      plt.ylim(-0.08,0.08)
      plt.xlim(-19,-10.5)
      plt.axhline(0.005,color='tab:orange')
      plt.legend(loc='best')
      plt.savefig('/Users/taweewat/Documents/red_sequence/pisco_color_plots/spread_model_real_i_fit_%s_%s.png' %
          (mode, field), dpi=120)
      plt.close(fig)
  
  
      dff0=dff
      dff0.to_csv("/Users/taweewat/Documents/pisco_code/slr_output/"+\
                  "galaxy_psf_total_%s.csv"%field)
      # dff_star0=pd.merge(dff_star, df0, on='NUMBER') # for non-SPREAD_MODEL
      dff_star0=dff_star #for SPREAD_MODEL
      dff_star0.to_csv("/Users/taweewat/Documents/pisco_code/slr_output/"+\
                       "star_psf_total_%s.csv"%field)
  
  def pisco_photometry_psf_v4(field, mode='psf', mode2mass='', slr=True):  #mode2mass: '' vs '_no2mass'
      def slr_running_psf(field, infile="None", mode="psf", mode2mass='', bigmacs="pisco_pipeline/big-macs-calibrate-master"):
          """
          slr_running: running SLR script from github.com/patkel/big-macs-calibrate to get a calibrated magnitude
          INPUT:
          - field: object of interset e.g., 'Field026'
          - bigmacs: the location for "big-macs-calibrate" directoty
          OUTPUT:
          - a new table with added columns with name MAG_g,...,MAGERR_g,...
          """
          slrdir = 'slr_output'
          pyfile = os.path.join(bigmacs, 'fit_locus.py')
          # cmd = "python %s --file %s --columns %s --extension 1 --bootstrap 15 -l -r ALPHA_J2000_i -d DELTA_J2000_i -j --plot=PLOTS_%s_%s" \
          #     % (pyfile, infile, os.path.join(bigmacs, "coadd_mag_sex_%s%s.columns"%(mode,'')), mode, field) 
          if mode2mass=='':
              cmd = "python %s --file %s --columns %s --extension 1 --bootstrap 15 -l -r ALPHA_J2000_i -d DELTA_J2000_i -j --plot=PLOTS_%s_%s" \
                  % (pyfile, infile, os.path.join(bigmacs, "coadd_mag_sex_%s%s.columns"%(mode,mode2mass)), mode, field)  #'' vs '_no2mass'
          elif mode2mass=='_no2mass':
              cmd = "python %s --file %s --columns %s --extension 1 --bootstrap 15 -l -r ALPHA_J2000_i -d DELTA_J2000_i --plot=PLOTS_%s_%s" \
                      % (pyfile, infile, os.path.join(bigmacs, "coadd_mag_sex_%s%s.columns"%(mode,mode2mass)), mode, field)  #'' vs '_no2mass'
          print cmd
          subprocess.check_call(shlex.split(cmd))
  
      def update_color(fname, table, mode='psf'):
          """
          update_color: using the output from SLR, update to the correct magnitude
          INPUT:
          - fname: input file from SLR output (...offsets.list)
          - table: the table that we want to update the value (from column magg,etc to MAG_g,etc)
          OUTPUT:
          - a new table with added columns with name MAG_g,...,MAGERR_g,...
          """
          print fname
          with open(fname) as f:
              content = f.readlines()
          content = [x.strip() for x in content]
          # print content
          if len(content)==8:
              red_content=content[4:]
          elif len(content)==10:    
              red_content=content[5:-1]
          # if len(content)==7:
          #     red_content=content[4:]
          # elif len(content)==9:    
          #     red_content=content[5:-1]
          band = [x.split(' ')[0][-1] for x in red_content]
          corr = [float(x.split(' ')[1]) for x in red_content]
          ecorr = [float(x.split(' ')[3]) for x in red_content]    
          print 'bands = ', band
  
          if mode=='psf':
              MODE1='PSF'
          elif mode=='model':
              MODE1='MODEL'
          elif mode=='auto':
              MODE1='AUTO'
          elif mode=='aper':
              MODE1='APER'
          elif mode=='hybrid':
              MODE1='HYBRID'
          elif mode=='iso':
              MODE1='ISO'
  
  
          table['MAG_' + band[0]] = table['MAG_%s_'%MODE1 + band[0]] + corr[0]
          table['MAG_' + band[1]] = table['MAG_%s_'%MODE1 + band[1]] + corr[1]
          table['MAG_' + band[2]] = table['MAG_%s_'%MODE1 + band[2]] + corr[2]
          table['MAG_' + band[3]] = table['MAG_%s_'%MODE1 + band[3]] + corr[3]
          table['MAGERR_' + band[0]] = (table['MAGERR_%s_'%MODE1 + band[0]]**2)**0.5# + ecorr[0]**2)**0.5
          table['MAGERR_' + band[1]] = (table['MAGERR_%s_'%MODE1 + band[1]]**2)**0.5# + ecorr[1]**2)**0.5
          table['MAGERR_' + band[2]] = (table['MAGERR_%s_'%MODE1 + band[2]]**2)**0.5# + ecorr[2]**2)**0.5
          table['MAGERR_' + band[3]] = (table['MAGERR_%s_'%MODE1 + band[3]]**2)**0.5# + ecorr[3]**2)**0.5
  
          return table
  
      slrdir = 'slr_output'
      df0=pd.read_csv("/Users/taweewat/Documents/pisco_code/slr_output/star_psf_total_%s.csv" % field,index_col=0)
      # if field=='SDSS603':
      #     df0=df0.drop([399,258,357,157,381,310,86,81,31,66,422,232,208,19,10])
      # elif field=='SDSS501':
      #     df0=df0.drop([265,108,196,213,160])
      # elif field=='SDSS123':
          # df0=df0.drop([68,5,61])
      # else:
      #     df0=df0    
      total3 = Table.from_pandas(df0)
      total3=total3[['NUMBER','ALPHA_J2000_i','DELTA_J2000_i','XWIN_IMAGE_i','YWIN_IMAGE_i',\
              'MAG_APER_i','MAGERR_APER_i','MAG_APER_g','MAGERR_APER_g','MAG_APER_r',\
              'MAGERR_APER_r','MAG_APER_z','MAGERR_APER_z','MAG_AUTO_i','MAGERR_AUTO_i',\
              'MAG_AUTO_g','MAGERR_AUTO_g','MAG_AUTO_r','MAGERR_AUTO_r','MAG_AUTO_z',\
              'MAGERR_AUTO_z','MAG_ISO_i','MAGERR_ISO_i','MAG_ISO_g','MAGERR_ISO_g','MAG_ISO_r',\
              'MAGERR_ISO_r','MAG_ISO_z','MAGERR_ISO_z',\
              'MAG_SPHEROID_i','MAGERR_SPHEROID_i','MAG_SPHEROID_g',\
              'MAGERR_SPHEROID_g','MAG_SPHEROID_r','MAGERR_SPHEROID_r','MAG_SPHEROID_z',\
              'MAGERR_SPHEROID_z','CLASS_STAR_i','CLASS_STAR_g','CLASS_STAR_r',\
              'CLASS_STAR_z','FLAGS_g','FLAGS_r','FLAGS_i','FLAGS_z','MAG_PSF_g',\
              'MAG_PSF_r','MAG_PSF_i','MAG_PSF_z','MAGERR_PSF_g','MAGERR_PSF_r',\
              'MAGERR_PSF_i','MAGERR_PSF_z','MAG_MODEL_g','MAG_MODEL_r',\
              'MAG_MODEL_i','MAG_MODEL_z','MAGERR_MODEL_g','MAGERR_MODEL_r',\
              'MAGERR_MODEL_i','MAGERR_MODEL_z','SPREAD_MODEL_g','SPREAD_MODEL_r',\
              'SPREAD_MODEL_i','SPREAD_MODEL_z','SPREADERR_MODEL_g','SPREADERR_MODEL_r',\
              'SPREADERR_MODEL_i','SPREADERR_MODEL_z']]
  
      print 'number of stars =', len(total3)
  
      if (mode2mass==''):
          starpsfmode = '_psf'
      elif (mode2mass=='_no2mass'):
          starpsfmode ='_no2mass'
  
      # total3.write(slrdir+'/star_psf%s_%s_%i.fits' % ('_psf',field,0), overwrite=True) #with 2MASS stars: star_psf_psf_%s_%i.fits
      total3.write(slrdir+'/star_psf%s_%s_%i.fits'%(starpsfmode,field,0),overwrite=True)  
      # no 2MASS star mode vs , '_psf' vs '_no2mass'
  
      if slr:
          slr_running_psf(field, infile=slrdir + '/star_psf%s_%s_%i.fits' %
                          (starpsfmode, field, 0), mode='psf', mode2mass=mode2mass)  # '_psf' vs '_no2mass'
  
      print 'mode=', mode, '/star_psf%s_%s_%i.fits.offsets.list' % (starpsfmode, field, 0)
  
      total_gal=Table.from_pandas(pd.read_csv("/Users/taweewat/Documents/pisco_code/slr_output/galaxy_psf_total_%s.csv"%(field)))
      ntotal_gal = update_color(slrdir+'/star_psf%s_%s_%i.fits.offsets.list' %
                                (starpsfmode, field, 0), total_gal, mode=mode)
      ntotal_gal.write(os.path.join(
          slrdir, 'galaxy_%s%s_ntotal_%s.csv'%(mode,mode2mass,field)), overwrite=True)
  
  
      total_star=Table.from_pandas(pd.read_csv("/Users/taweewat/Documents/pisco_code/slr_output/star_psf_total_%s.csv"%(field)))
      ntotal_star = update_color(slrdir+'/star_psf%s_%s_%i.fits.offsets.list'%
                                (starpsfmode, field, 0), total_star, mode=mode)
      ntotal_star.write(os.path.join(
          slrdir, 'star_%s%s_ntotal_%s.csv'%(mode,mode2mass,field)), overwrite=True)
  
  
  def make_images(field,ax=None):
      dir='/Users/taweewat/Documents/pisco_code/Chips_images/'
      try:
          ax.imshow(image.imread(dir+"aplpy4_%s_img4.jpeg"%field))
      except:
          ax.imshow(image.imread(dir+"aplpy4_%s_img.jpeg"%field))
      # ax.imshow(image.imread(dir+"aplpy4_%s_img4.jpeg"%field))
      ax.axes.get_xaxis().set_visible(False)
      ax.axes.get_yaxis().set_visible(False)
      ax.axis('off')
      return None
  
  # def sur_pro(r): #Mpc
  #     def fn(x):
  #         if x>=1:
  #             return 1.-(2/np.sqrt(x**2-1)*np.arctan(np.sqrt((x-1.)/(x+1.))))
  #         elif x<1:
  #             return 1.-(2/np.sqrt(1-x**2)*np.arctanh(np.sqrt((1.-x)/(x+1.))))
  #     rs=0.15/0.71 #Mpc
  #     if r>=(0.1/0.71):
  #         return 1/((r/rs)**2-1)*fn(r/rs)
  #     elif r<(0.1/0.71):
  #         return 1./(((0.1/0.71)/rs)**2-1)*fn((0.1/0.71)/rs)
  # def k_NFW():
  #     def integrated(y):
  #         return 1./integrate.quad(lambda r: 2*np.pi*r*sur_pro(r),0,y)[0]
  #     xy=np.logspace(-3,3,num=30)
  #     X = np.log(xy)
  #     Y = np.log([integrated(np.e**(y)) for y in X])
  #     Z=np.polyfit(X,Y,6)
  #     k_NFW = np.poly1d(Z)
  #     return k_NFW
  # def sur_pro_prob(r,rc,k_NFW): #(Mpc,Mpc) # Weighted based on the distance from the center (Rykoff+12)
  #     return np.e**(k_NFW(np.log(rc)))*sur_pro(r)
  
  
  name=['z','dist','age','mass','Abs_g','App_g','kcorr_g','Abs_r',\
        'App_r','kcorr_r','Abs_i','App_i','kcorr_i','Abs_z','App_z','kcorr_z']
  df=pd.read_csv('/Users/taweewat/Documents/red_sequence/rsz/model/'+\
  # 'ezmodel2_bc03_zf2.5_chab_0.016_exp_0.1.txt',
  'ezmodel2_bc03_zf2.5_chab_0.02_exp_0.1.txt',
              #    'ezmodel2_c09_zf3.0_chab_0.02_exp_0.1.txt',
                 skiprows=27,delim_whitespace=True,names=name)
  df=df[(df.z>=0.1) & (df.z<1.)]
  z_new=np.arange(0.1, 0.95, 0.0025)
  Appi_new = interpolate.splev(z_new, interpolate.splrep(df.z, df.App_i, s=0), der=0)
  Appi_f = interpolate.interp1d(df.z, df.App_i, kind='cubic')
  
  #all extra options
  extra_name= 'gnorm_zf2.5_bc03_noebv_auto_bin1.0_root15_sur0.25' #'gremove_lum_silk_zf2.5_c09_11', 'gremove_silk_zf3_c09_noebv_model_complete_no2mass'
  core_radius=0.25
  gremove = False # remove non-detect g objects from the list
  duplicate = False # remove duplicate redshift (uncertain)
  colorerr = True  # add redshift with color_error taken into account
  transparent = True # make transparent plot for flip book
  img_filp = False  # make image flip from transparent
  img_redshift = True # make image with redshift for each object
  
  def linear_rmi(x0,redshift):
      x=df.z[:-11] #-12
      y=(df.App_r-df.App_i)[:-11] #-12
      yhat = np.polyfit(x, y, 5) #5 vs 9
      f_rmi = np.poly1d(yhat)
      slope=-0.0222174237562*1.007
      # Appi0=Appi_new[np.where(abs(z_new-redshift)<=1e-9)[0][0]]
      Appi0=Appi_f(redshift)
      return slope*(x0-Appi0)+f_rmi(redshift)
  
  def linear_gmr(x0,redshift):
      x=df.z[:-24] #-25
      y=(df.App_g-df.App_r)[:-24] #-25
      yhat = np.polyfit(x, y, 5)
      f_gmr = np.poly1d(yhat)
      slope=-0.0133824600874*1.646
      # Appi0=Appi_new[np.where(abs(z_new-redshift)<=1e-9)[0][0]]
      Appi0=Appi_f(redshift)
      return slope*(x0-Appi0)+f_gmr(redshift)
  
  def linear_gmi(x0,redshift):
      x=df.z[:-9] 
      y=(df.App_g-df.App_i)[:-9] 
      yhat = np.polyfit(x, y, 5)
      f_gmi = np.poly1d(yhat) 
      Appi0=Appi_f(redshift)
      slope = -0.04589707934164738 * 1.481
      return slope*(x0-Appi0)+f_gmi(redshift)
  
  def find_fits_dir(field):
      home = '/Users/taweewat/Documents/pisco_code/'
      dirs = ['ut170103/', 'ut170104/', 'ut170619/', 'ut170621/',\
              'ut170624/', 'ut171208/', 'ut171209/', 'ut171212/']
      myReg = re.compile(r'(%s_A).*' % field)
      for di in dirs:
          diri = home + di
          for text in os.listdir(diri):
              if myReg.search(text) != None:
                  # filename = myReg.search(text).group()
                  allfilename = diri
      return allfilename
  
  dir_dict = dict(zip(['ut170103/','ut170104/','ut170619/',\
  'ut170621/','ut170624/','ut171208/','ut171209/','ut171212/'], np.arange(1, 9)))
  
  def find_ra_dec(field):
      if field == 'PKS1353':
          RA = 209.0225
          DEC = -34.3530556
          redshift = 0.223
      elif field == 'CHIPS2249-2808': #CHIPS2227-4333
          # RA = 336.99975202151825
          # DEC = -43.57623068466675
          RA = 336.98001
          DEC = -43.56472
          redshift = -1
      elif field == 'CHIPS2246-2854': #'CHIPS2223-3455'
          # RA = 335.7855174238757
          # DEC = -34.934569299688185
          RA = 335.78
          DEC = -34.9275
          redshift = -1
      elif field[0:5] == 'Field':
          base = pd.read_csv(
              '/Users/taweewat/Dropbox/Documents/MIT/Observation/2017_1/all_objs.csv')
          RA = base[base.name == field].ra.values[0]
          DEC = base[base.name == field].dec.values[0]
          redshift = base[base.name == field].redshift.values[0]
      elif field[0:5] == 'CHIPS':
          base = pd.read_csv(
              '/Users/taweewat/Documents/red_sequence/chips_all_obj.csv', index_col=0)
          RA = base[base.chips == field].ra.values[0]
          DEC = base[base.chips == field].dec.values[0]
          redshift = base[base.chips == field].redshift.values[0]
      elif field[0:4] == 'SDSS':
          base = pd.read_csv(
              '/Users/taweewat/Documents/xray_project/ned-result/final_sdss_cut5.csv', index_col=0)
          RA = base[base.name == field].RA.values[0]
          DEC = base[base.name == field].DEC.values[0]
          redshift = base[base.name == field].redshift.values[0]
      return RA, DEC, redshift
  
  
  def pisco_tilt_resequence(field, mode='psf', mode2mass=''):
      RA, DEC, redshift = find_ra_dec(field)
  
      if redshift!=-1:
          qso_redshift=redshift
      else:
          qso_redshift=0.2
  
      print 'RA', RA
      print 'DEC', DEC
  
      ebv = ebvpy.calc_ebv(ra=[RA],dec=[DEC]); print 'ebv:', ebv[0]
      # ebv_g=ebvpy.calc_color_correction('g', ebv)[0]
      # ebv_r=ebvpy.calc_color_correction('r', ebv)[0]
      # ebv_i=ebvpy.calc_color_correction('i', ebv)[0]
      # ebv_z=0.0
      ebv_g,ebv_r,ebv_i,ebv_z=0.0,0.0,0.0,0.0  #no longer use reddening correction because it is already included in SLR
      print 'ebv_g:', ebv_g, 'ebv_r:', ebv_r, 'ebv_i:', ebv_i
  
      param_izp=read_param_izp(mode) #i zero point
  
      # fname = "/Users/taweewat/Documents/pisco_code/slr_output/galaxy_ntotal_%s.csv"%field
      dir_slrout='/Users/taweewat/Documents/pisco_code/slr_output/'
      fname = dir_slrout+"galaxy_%s%s_ntotal_%s.csv" % (
          mode, mode2mass, field)  # '' vs '_no2mass'
      df0 = pd.read_csv(fname,index_col=0)
      if gremove:
          nog=len(df0[df0['MAG_PSF_g'] >= 50.]); print "no g detected:", nog
          df0 = df0[df0['MAG_PSF_g'] < 50.].copy()  # cut out not detected objects in g band
      else: 
          nog=0
  
      c5 = SkyCoord(ra=df0['ALPHA_J2000_i'].values*u.degree, dec=df0['DELTA_J2000_i'].values*u.degree)
      c0 = SkyCoord(ra=RA*u.degree, dec=DEC*u.degree)
      sep = c5.separation(c0)
      df0['sep(deg)']=sep
      df0['sep(Mpc)']=sep*60.*cosmo.kpc_proper_per_arcmin(qso_redshift).value/1e3
      cut=df0
  
      dfi = cut#.drop_duplicates(subset=['XWIN_WORLD', 'YWIN_WORLD'], keep='first').copy()
      print 'duplicates:', len(df0), len(dfi)
      
      # Added Galactic Reddening (6/16/18)
      if mode2mass == '':
          dfi['MAG_i']=dfi['MAG_i']-ebv_i
          dfi['MAG_g']=dfi['MAG_g']-ebv_g
          dfi['MAG_r']=dfi['MAG_r']-ebv_r
      # Use i Zero Point from each day and g,r zero point fron the color (6/22/18)
      elif mode2mass == '_no2mass':
          mag0 = param_izp['i_zp_day%i'%dir_dict[find_fits_dir(field)[-9:]]]
          print 'i_zp_day', find_fits_dir(field), mag0
          dfi['MAG_i']=dfi['MAG_i']-ebv_i+mag0
          dfi['MAG_g']=dfi['MAG_g']-ebv_g+mag0
          dfi['MAG_r']=dfi['MAG_r']-ebv_r+mag0
          dfi['MAG_z']=dfi['MAG_z']-ebv_z+mag0
  
      dfi.to_csv(dir_slrout+"galaxy_%s_final_%s.csv"%(mode,field))
      
      fname=dir_slrout+"star_%s%s_ntotal_%s.csv" % (mode, mode2mass, field)
      df0=pd.read_csv(fname,index_col=0)
      dfi=df0
      # Added Galactic Reddening (6/16/18)
      if mode2mass == '':
          dfi['MAG_i']=dfi['MAG_i']-ebv_i
          dfi['MAG_g']=dfi['MAG_g']-ebv_g
          dfi['MAG_r']=dfi['MAG_r']-ebv_r
      # Use i Zero Point from each day and g,r zero point fron the color (6/22/18)
      elif mode2mass == '_no2mass':
          mag0 = param_izp['i_zp_day%i'%dir_dict[find_fits_dir(field)[-9:]]]
          print 'i_zp_day', find_fits_dir(field), mag0
          dfi['MAG_i']=dfi['MAG_i']-ebv_i+mag0
          dfi['MAG_g']=dfi['MAG_g']-ebv_g+mag0
          dfi['MAG_r']=dfi['MAG_r']-ebv_r+mag0
          dfi['MAG_z']=dfi['MAG_z']-ebv_z+mag0
      dfi.to_csv(dir_slrout+"star_%s_final_%s.csv"%(mode,field))
          
      
      
      return None
  
      # dfi=dfi[dfi['MAG_i']<21.5].copy()
      # dfi=dfi[dfi.MAGERR_g<0.5]
      # dfi=dfi[(dfi.MAG_g<100)&(dfi.MAG_i<100)&(dfi.MAG_r<100)]
      # dfi=dfi[(dfi.FLAGS_g<5)&(dfi.FLAGS_r<5)&(dfi.FLAGS_i<5)&(dfi.FLAGS_z<5)]
  
  def xxx(x):
      dfi=dfi[np.sqrt(dfi['MAGERR_r']**2+dfi['MAGERR_i']**2)<0.3].copy() #0.5
      x=dfi['MAG_i']
      y=np.sqrt(dfi['MAGERR_r']**2+dfi['MAGERR_i']**2)
      p=np.poly1d(np.polyfit(x,np.log(y),1, w=np.sqrt(y)))
      Mag_cut=(p-np.log(0.067*1.5)).roots;
      print "Mag_cut: %.2f"%(Mag_cut)
      xs=np.arange(np.min(x),np.max(x),0.01)
  
      fig,ax=plt.subplots(figsize=(5,5))
      plt.plot(x,y,'.',label='r-i')
      plt.plot(xs,np.exp(p(xs)),label='exp({:.1f}+{:.1f}x)'.format(p[0],p[1]))
      plt.xlabel('Mag_i'); plt.ylabel('$\Delta r-i$ err')
      plt.ylim(-0.05,0.35)
      plt.axvline(Mag_cut,label='Mag_cut')
      plt.legend(loc='best')
      plt.savefig('/Users/taweewat/Documents/red_sequence/pisco_color_plots/uncer_%s_%s.png' %
          (mode, field), dpi=120)
      # plt.tight_layout()
      plt.close(fig)
      #Magnitude cut
  
      print field, qso_redshift, df0.shape, cut.shape, dfi.shape, dfi['sep(deg)'].max(), dfi['sep(Mpc)'].max()
  
      norm = matplotlib.colors.Normalize(vmin=0.10,vmax=0.675)
      c_m = matplotlib.cm.cool
      s_m = matplotlib.cm.ScalarMappable(cmap=c_m, norm=norm)
      s_m.set_array([])
  
      I=np.arange(16,24,0.01)
      dfi.loc[:,"z_gmr"] = np.nan
      dfi.loc[:,"z_rmi"] = np.nan
      dfi.loc[:,"w_gmr"] = np.nan
      dfi.loc[:,"w_rmi"] = np.nan
      dfi.loc[:,"w_col_gmr"] = np.nan
      dfi.loc[:,"w_col_rmi"] = np.nan
      # dfi.loc[:,"z_gmi"] = np.nan
      # dfi.loc[:,"w_gmi"] = np.nan
      # dfi.loc[:,"w_col_gmi"] = np.nan
  
      # k_NFW0=k_NFW()
  
      bin_width=0.035 #0.025
  
      # bins_gmr_cen = np.arange(0.15815-0.0175, 0.33315-0.0175+0.01, bin_width)
      # bins_gmr_edge = np.arange(0.14065-0.0175, 0.35065-0.0175+0.01, bin_width)
  
      # bins_gmr_cen = np.arange(0.12315+0.0175, 0.33315+0.0175+0.01, bin_width)
      # bins_gmr_edge = np.arange(0.10565+0.0175, 0.35065+0.0175+0.01, bin_width)
      # bins_rmi_cen = np.arange(0.36815-0.0175, 0.64815-0.0175+0.01, bin_width) 
      # bins_rmi_edge = np.arange(0.35065-0.0175, 0.66565-0.0175+0.01, bin_width)
  
      #new one: combine last gmr with "new" rmi
      bins_gmr_cen = np.arange(0.12315, 0.33315+0.01, bin_width)
      bins_gmr_edge = np.arange(0.10565, 0.35065+0.01, bin_width)
      bins_rmi_cen = np.arange(0.36815-bin_width, 0.64815+0.01, bin_width) 
      bins_rmi_edge = np.arange(0.35065-bin_width, 0.66565+0.01, bin_width) 
  
      z_rmi,w_rmi,w_col_rmi=[],[],[]
      for i, row in dfi.iterrows():
          for z in bins_rmi_cen:
              # if row['MAG_i'] < -18+5.*np.log10(ex.d_L(z)*1e6)-5.:
               # if row['MAG_i'] < magi_cut_rmi:
              # if np.sqrt(row['MAGERR_r']**2+row['MAGERR_i']**2)<0.134: #np.mean(f_rmi(x+0.07)-f_rmi(x))
              # if np.sqrt(row['MAGERR_r']**2+row['MAGERR_i']**2)<0.067*1.5: #0.067*1.5
              if row['MAG_i'] < Mag_cut:
                  rmi=row['MAG_r']-row['MAG_i']
                  # rmierr=np.sqrt(row['MAGERR_r']**2+row['MAGERR_i']**2)
                  low_edge=linear_rmi(row['MAG_i'],round(z-0.0175,4)) #0.0125
                  high_edge=linear_rmi(row['MAG_i'],round(z+0.0175,4)) #0.0125
                  if (rmi > low_edge) & (rmi <= high_edge):
                      # if (np.sqrt(row['MAGERR_r']**2+row['MAGERR_i']**2) < 3.5*(high_edge-low_edge)):
                      z_rmi.append(round(z,3))
                      # wrmi0=sur_pro_prob(row['sep(Mpc)'],1.,k_NFW0)
                      wrmi0=ex.sur_pro_prob_ang(row['sep(deg)']*60, core_radius); w_rmi.append(wrmi0) #arcmin
                      # w_col_rmi0=scipy.stats.norm(rmi,rmierr).cdf(high_edge)-scipy.stats.norm(rmi,rmierr).cdf(low_edge); w_col_rmi.append(w_col_rmi0)
                      w_col_rmi0=1.; w_col_rmi.append(w_col_rmi0)
                      dfi.loc[i,"z_rmi"]=z
                      dfi.loc[i,"w_rmi"]=wrmi0
                      dfi.loc[i,"w_col_rmi"]=w_col_rmi0
  
      z_gmr,w_gmr,w_col_gmr=[],[],[]
      for i, row in dfi.iterrows():
          for z in bins_gmr_cen:
              # if row['MAG_i'] < -18+5.*np.log10(ex.d_L(z)*1e6)-5.:
              # if row['MAG_i'] < magi_cut_gmr:
              # if np.sqrt(row['MAGERR_g']**2+row['MAGERR_r']**2)<0.165: #np.mean(f_gmr(x+0.07)-f_gmr(x))
              # if np.sqrt(row['MAGERR_g']**2+row['MAGERR_r']**2)<0.0825*1.5: #0.0825*1.5
              if row['MAG_i'] < Mag_cut:
                  gmr=row['MAG_g']-row['MAG_r']
                  # gmrerr=np.sqrt((row['MAGERR_g'])**2+row['MAGERR_r']**2) #add factor 2.2 to reduce the g error to be similar to other bands
                  low_edge=linear_gmr(row['MAG_i'],round(z-0.0175,4)) #0.0125
                  high_edge=linear_gmr(row['MAG_i'],round(z+0.0175,4)) #0.0125
                  if (gmr > low_edge) & (gmr <= high_edge):
                      # if (np.sqrt(row['MAGERR_g']**2+row['MAGERR_r']**2) < 3.5*(high_edge-low_edge)):
                      z_gmr.append(round(z,3))
                      # w_col_gmr0=scipy.stats.norm(gmr,gmrerr).cdf(high_edge)-scipy.stats.norm(gmr,gmrerr).cdf(low_edge); w_col_gmr.append(w_col_gmr0)
                      w_col_gmr0=1.; w_col_gmr.append(w_col_gmr0)
                      # wgmr0=sur_pro_prob(row['sep(Mpc)'],1.,k_NFW0); w_gmr.append(wgmr0)
                      wgmr0 = ex.sur_pro_prob_ang(row['sep(deg)'] * 60, core_radius); w_gmr.append(wgmr0)  # arcmin
                      dfi.loc[i,"z_gmr"]=z
                      dfi.loc[i,"w_gmr"]=wgmr0
                      dfi.loc[i,"w_col_gmr"]=w_col_gmr0
  
      ns1,xs1=np.histogram(z_gmr,bins=bins_gmr_edge,weights=np.array(w_gmr)*np.array(w_col_gmr)) #0.15-0.325
      bin_cen1 = (xs1[:-1] + xs1[1:])/2
      ns2,xs2=np.histogram(z_rmi,bins=bins_rmi_edge,weights=np.array(w_rmi)*np.array(w_col_rmi)) #0.36-0.675
      bin_cen2 = (xs2[:-1] + xs2[1:])/2
      # z_total=np.append(bin_cen1, bin_cen2)
      # n_total=np.append(ns1,ns2)
      z_total=np.append(bin_cen1, bin_cen2[1:])
      n_total=np.append(np.append(ns1[:-1],np.array(ns1[-1]+ns2[0])),np.array(ns2[1:]))  
      z_max=z_total[np.where(n_total==np.max(n_total))[0][0]]
      n_median = np.median(n_total[n_total != 0])
      n_mean = np.mean(n_total)
  
      n_bkg = np.mean(sorted(n_total)[2:-2]);
      z_total_added = np.insert(
          np.append(z_total, z_total[-1] + bin_width), 0, z_total[0] - bin_width)
      n_total_added = np.insert(np.append(n_total, 0), 0, 0) - n_bkg
      # print 'n_total_added', n_total_added
      
      lumfn=pd.read_csv('/Users/taweewat/Documents/red_sequence/coma_cluster_luminosity_function/schecter_fn.csv',\
      names=['M_r','theta(M)Mpc^-3']) 
      h=0.7
      x=lumfn['M_r']+5*np.log10(h);
      y=lumfn['theta(M)Mpc^-3']*(h**3)
      f1d=interp1d(x, y,kind='cubic')
      def lum_function(M):
          alpha = -1.20
          Nb = np.log(10) / 2.5 * 0.002 * (70 / 50.)**3
          Mb_s = -21. + 5 * np.log10(70 / 50.)
          return Nb * (10.**(0.4 * (alpha + 1) * (Mb_s - M))) * np.exp(-10.**(0.4 * (Mb_s - M)))
  
      lum_fn = lambda z: integrate.quad( f1d, -23.455, ex.abs_mag(22.25, z))[0]
      lum_vfn = np.vectorize(lum_fn)
      dense_fn = lambda z: integrate.quad(ex.NFW_profile,0.001,cosmo.kpc_proper_per_arcmin(z).value/1e3)[0]
      dense_vfn = np.vectorize(dense_fn)
      n_total_adj=n_total_added #/(lum_vfn(z_total_added)*dense_vfn(z_total_added)) (adjusted the peak before picking it)
      print 'n_total_added:', n_total_added
      print 'n_total_adj:', n_total_adj
  
      indi = np.where(n_total_adj == np.max(n_total_adj))[0][0]
      # indi = np.where(n_total_added == np.max(n_total_added))[0][0]
      z_fit = z_total_added[[indi - 1, indi, indi + 1]]; print 'z_fit', z_fit
      n_fit = n_total_added[[indi - 1, indi, indi + 1]]; print 'n_fit', n_fit
      def gaussian_func(x, a, mu):
          sigma=0.035
          return a * np.exp(-(x-mu)**2/(2*(sigma**2)))
      
      if (n_fit[0]<0.) and (n_fit[2]<0.):
          popt, pcov = curve_fit(gaussian_func, z_fit, [0,n_fit[1],0], p0=[n_fit[1],z_fit[1]])
      else:
          popt, pcov = curve_fit(gaussian_func, z_fit,
                                 n_fit, p0=[n_fit[1], z_fit[1]])
      # signal=tuple(popt)[0]
  
      # def v_func(z):
      #     return (z**2+2*z)/(z**2+2*z+2)
      # signal=((np.max(n_total)-np.mean(n_total))*(v_func(z_max)*(4000))**2)/5.3e6 #normalization for r~1 at z~0.15
      # signal = (
      #     (tuple(popt)[0]) * (cosmo.luminosity_distance(tuple(popt)[1]).value)**1.5) / 5.3e5  # normalization for r~1 at z~0.15
  
      def lum_function(M):
          alpha = -1.20
          Nb = np.log(10) / 2.5 * 0.002 * (70 / 50.)**3
          Mb_s = -21. + 5 * np.log10(70 / 50.)
          return Nb * (10.**(0.4 * (alpha + 1) * (Mb_s - M))) * np.exp(-10.**(0.4 * (Mb_s - M)))
  
      lumfn=pd.read_csv('/Users/taweewat/Documents/red_sequence/coma_cluster_luminosity_function/schecter_fn.csv',\
      names=['M_r','theta(M)Mpc^-3']) 
      h=0.7
      x=lumfn['M_r']+5*np.log10(h);
      y=lumfn['theta(M)Mpc^-3']*(h**3)
      f1d=interp1d(x, y,kind='cubic')
  
  
      z_max_fit = tuple(popt)[1]
      # lum_factor = integrate.quad(lum_function, -24, abs_mag(21.60, tuple(popt)[1]))[0]
      # lum_factor = cosmo.luminosity_distance(tuple(popt)[1]).value**-1.5*100
      lum_factor = integrate.quad( f1d, -23.455, ex.abs_mag(22.25, z_max_fit))[0]  
      #-23.455: min abs Mag from schecter_fn.csv, 22.25: median of Mag r
      density_factor=integrate.quad(ex.NFW_profile, 0.001, core_radius*cosmo.kpc_proper_per_arcmin(z_max_fit).value/1e3)[0]
      signal = tuple(popt)[0] / (lum_factor * density_factor)
      
  
      print 'z_max_fit', z_max_fit
      print 'lum_factor:', lum_factor
      print 'density_factor', density_factor
  
      # duplicate=False ## set duplication
      n_total_dup=0
  
      ## Plot the figure
      cmap=matplotlib.cm.RdYlGn
      if duplicate or colorerr:
          fig,ax=plt.subplots(1,5,figsize=(25,5))
      else: 
          fig,ax=plt.subplots(1,4,figsize=(20,5))
      make_images(field,ax[0])
  
      norm = matplotlib.colors.Normalize(vmin=0.01,vmax=2)
      dfi_ri=dfi.loc[dfi['z_rmi'].dropna().index]
      ax[1].scatter(dfi['MAG_i'],dfi['MAG_r']-dfi['MAG_i'],c='black',alpha=0.1)#dfi['w_rmi'],cmap=cmap)
      ax[1].scatter(dfi_ri['MAG_i'],dfi_ri['MAG_r']-dfi_ri['MAG_i'],c=dfi_ri['w_rmi'],cmap=cmap)#,norm=norm)
      ax[1].errorbar(dfi_ri['MAG_i'],dfi_ri['MAG_r']-dfi_ri['MAG_i'],xerr=dfi_ri['MAGERR_i'],yerr=np.sqrt(dfi_ri['MAGERR_r']**2+dfi_ri['MAGERR_i']**2),fmt='none',c='k',alpha=0.05)
      # plt.plot(df.App_i,df.App_r-df.App_i,'.')
      # ax[1].axhline(xs[:-1][(xs[:-1]<1.33) & (xs[:-1]>0.6)][0],lw=0.7,color='green')
      for z in bins_rmi_cen:
          ax[1].plot(I,linear_rmi(I,round(z,4)),color=s_m.to_rgba(z))
      ax[1].set_ylim(0.25,1.5)
      ax[1].set_xlim(16,24)
      # cbar=plt.colorbar(s_m)
      ax[1].set_xlabel('I')
      ax[1].set_ylabel('R-I')
      ax[1].set_title('z=0.35-0.675')#, icut:'+str(magi_cut_rmi))
  
      # plt.plot([corr_f(z) for z in df.z.values[5:-12]],df.App_r[5:-12]-df.App_i[5:-12],'-')
      dfi_gr=dfi.loc[dfi['z_gmr'].dropna().index]
      ax[2].scatter(dfi['MAG_i'],dfi['MAG_g']-dfi['MAG_r'],c='black',alpha=0.1)#,c=dfi['w_gmr'],cmap=cmap)
      ax[2].scatter(dfi_gr['MAG_i'],dfi_gr['MAG_g']-dfi_gr['MAG_r'],c=dfi_gr['w_gmr'],cmap=cmap)#,norm=norm)
      ax[2].errorbar(dfi_gr['MAG_i'],dfi_gr['MAG_g']-dfi_gr['MAG_r'],xerr=dfi_gr['MAGERR_i'],yerr=np.sqrt(dfi_gr['MAGERR_g']**2+dfi_gr['MAGERR_r']**2),fmt='none',c='k',alpha=0.05)
      # plt.plot(df.App_i,df.App_g-df.App_r,'.')
      # ax[2].axhline(xs[:-1][(xs[:-1]<1.65) & (xs[:-1]>np.min(x2))][0],lw=0.7,color='green')
      for z in bins_gmr_cen:
          ax[2].plot(I,linear_gmr(I,round(z,4)),color=s_m.to_rgba(z))
      ax[2].set_ylim(0.75,2)
      ax[2].set_xlim(16,24)
      # cbar=plt.colorbar(s_m)
      ax[2].set_xlabel('I')
      ax[2].set_ylabel('G-R')
      ax[2].set_title('z=0.15-0.325')
      # plt.plot([corr_f(z) for z in df.z.values[:-25]],df.App_g[:-25]-df.App_r[:-25],'-')
  
      xs=np.arange(np.min(z_fit)-0.1,np.max(z_fit)+0.1,0.001)
      ax[3].bar(bin_cen2, ns2, width=bin_width, color='#1f77b4', alpha=1.0)
      ax[3].bar(bin_cen1, ns1, width=bin_width, color='#ff7f0e', alpha=1.0)
      ax[3].bar(z_total, n_total, width=bin_width, color='grey', alpha=0.5)
      ax[3].axvline(0.3525,ls='--')
      ax[3].axvline(z_max,ls='--',color='purple',label='z_max:%.2f'%z_max)
      ax[3].axvline(redshift,color='red',label='z:%.2f'%redshift)
      ax[3].plot(z_fit,n_fit+n_bkg,'o',c='tab:purple')
      ax[3].plot(xs, gaussian_func(xs, *popt)+n_bkg, c='tab:green', ls='--', label='fit: a=%.2f, mu=%.4f'% tuple(popt))
      ax[3].axhline(n_median,color='tab:green',label='median:%.2f'%n_median)
      ax[3].axhline(n_mean,color='tab:red',label='mean:%.2f'%n_mean)
      ax[3].legend(loc='best')
      ax[3].set_xlabel('z')
      ax[3].set_xlim(0.1,0.7)
      ax[3].set_title('ebv:%.3f,ebv_g-r:-%.3f,ebv_r-i:-%.3f'%(ebv[0],ebv_g-ebv_r,ebv_r-ebv_i))
      if np.max(n_total)<30:
          ax[3].set_ylim(0,30)
  
      if duplicate:
          xs = np.arange(np.min(z_fit_dup) - 0.1, np.max(z_fit_dup) + 0.1, 0.001)
          ax[4].bar(bin_cen2, ns2-ns_dup2, width=bin_width, color='#1f77b4') #widht = 0.025
          ax[4].bar(bin_cen1, ns1-ns_dup1, width=bin_width, color='#ff7f0e') #width = 0.025
          ax[4].axvline(z_max,ls='--',color='purple',label='z_max:%.2f'%z_max)
          ax[4].axvline(redshift,color='red',label='z:%.2f'%redshift)
          ax[4].plot(z_fit_dup,n_fit_dup+n_bkg_dup,'o',c='tab:purple')
          ax[4].plot(xs, gaussian_func(xs, *popt_dup)+n_bkg_dup, c='tab:green', ls='--', label='fit: a=%.2f, mu=%.4f'% tuple(popt))
          ax[4].legend(loc='best')
          ax[4].set_xlabel('z')
          ax[4].set_xlim(0.1,0.7)
          if np.max(n_total)<30:
              ax[4].set_ylim(0,30)
          
      if colorerr:
          dfi_rmi = dfi[~np.isnan(dfi['z_rmi'])]
          dfi_gmr = dfi[~np.isnan(dfi['z_gmr'])]
          zs_gmr = np.arange(0.11, 0.3425, 0.002)
          zs_rmi = np.arange(0.3425, 0.65, 0.002)
  
          ntot_rmi = np.repeat(0, len(zs_rmi))
          ntot_gmr = np.repeat(0, len(zs_gmr))
          for i, row in dfi_rmi.iterrows():
              # for i, row in dfi.iterrows():
              i0 = row['MAG_i']
              rmi = row['MAG_r'] - row['MAG_i']
              rmierr = np.sqrt((row['MAGERR_r'])**2 + row['MAGERR_i']**2)
              ntot_rmi0 = scipy.stats.norm(rmi, rmierr).pdf(
                  linear_rmi(i0, zs_rmi))
              ntot_rmi = ntot_rmi + ntot_rmi0 * row['w_rmi']
              ax[4].plot(zs_rmi,ntot_rmi0*row['w_rmi'],'-',color='tab:red',alpha=0.2)
  
          for i, row in dfi_gmr.iterrows():
              # for i, row in dfi.iterrows():
              i0 = row['MAG_i']
              gmr = row['MAG_g'] - row['MAG_r']
              gmrerr = np.sqrt((row['MAGERR_g'])**2 + row['MAGERR_r']**2)
              ntot_gmr0 = scipy.stats.norm(gmr, gmrerr).pdf(
                  linear_gmr(i0, zs_gmr))
              ntot_gmr = ntot_gmr + ntot_gmr0 * row['w_gmr']
              ax[4].plot(zs_gmr,ntot_gmr0*row['w_gmr'],'-',color='tab:cyan',alpha=0.2)
          
          ax[4].plot(zs_gmr, ntot_gmr, '.')
          ax[4].plot(zs_rmi, ntot_rmi, '.')
          ax[4].axvline(z_max,ls='--',color='purple',label='z_max:%.2f'%z_max)
          ax[4].axvline(redshift,color='red',label='z:%.2f'%redshift)
          ax[4].legend(loc='best')
          ax[4].set_xlabel('z')
          ax[4].set_xlim(0.1, 0.7)
          if np.max(np.append(ntot_gmr,ntot_rmi)) < 200:
              ax[4].set_ylim(0, 200)
          n_total_cerr=np.append(ntot_gmr,ntot_rmi)
      else:
          n_total_cerr=0
  
      signal_final = signal_dup if duplicate else signal
  
      plt.tight_layout(rect=[0, 0., 1, 0.98])
      purge('/Users/taweewat/Documents/red_sequence/pisco_color_plots/',
            'redsq_richg%s_%s_all_.*_%s_tilted.png' % ('', mode, field))
      plt.savefig('/Users/taweewat/Documents/red_sequence/pisco_color_plots/redsq_richg%s_%s_all_%.3f_%s_tilted.png' % ('',mode,signal_final,field), dpi=120)
      plt.close(fig)
  
      # fig,ax=plt.subplots(1,4,figsize=(20,5))
      # make_images(field,ax[0])
      # dfi_gmi=dfi[~np.isnan(dfi['z_gmi'])]
      # zs_gmi=np.arange(0.115,0.69,0.002)
      # ntot_gmi=np.repeat(0,len(zs_gmi))
      # for i, row in dfi_gmi.iterrows():
      #     i0 = row['MAG_i']
      #     gmi = row['MAG_g'] - row['MAG_i']
      #     gmierr = np.sqrt((row['MAGERR_g'])**2 + row['MAGERR_i']**2)
      #     ntot_gmi0 = scipy.stats.norm(gmi, gmierr).pdf(
      #         linear_gmi(i0, zs_gmi))
      #     ntot_gmi = ntot_gmi + ntot_gmi0 * row['w_gmi']
      #     ax[3].plot(zs_gmi,ntot_gmi0*row['w_gmi'],'-',color='tab:cyan',alpha=0.2)
      # ax[1].scatter(dfi['MAG_i'],dfi['MAG_g']-dfi['MAG_i'],c='black',alpha=0.1)#dfi['w_rmi'],cmap=cmap)
      # ax[1].scatter(dfi_gmi['MAG_i'],dfi_gmi['MAG_g']-dfi_gmi['MAG_i'],c=dfi_gmi['w_gmi'],cmap=cmap)
      # ax[1].errorbar(dfi_gmi['MAG_i'], dfi_gmi['MAG_g'] - dfi_gmi['MAG_i'], xerr=dfi_gmi['MAGERR_i'],
      #                yerr=np.sqrt(dfi_gmi['MAGERR_g']**2 + dfi_gmi['MAGERR_i']**2), fmt='none', c='k', alpha=0.05)
      # for z in np.arange(0.15, 0.71, bin_width):
      #     ax[1].plot(I,linear_gmi(I,z),color=s_m.to_rgba(z))
      # ax[1].set_ylim(1.0,3.5)   
      # ax[1].set_xlim(16,24)   
      # ax[1].set_xlabel('I')
      # ax[1].set_ylabel('G-I')
      # ax[1].set_title('z=0.15-0.675')
      # ns3,xs3=np.histogram(z_gmi,bins=np.arange(0.1325,0.7,0.035),weights=np.array(w_gmi)*np.array(w_col_gmi))
      # bin_cen3 = (xs3[:-1] + xs3[1:])/2
      # z_max_gmi = bin_cen3[np.where(ns3 == np.max(ns3))[0][0]]
      # n_bkg = np.mean(sorted(ns3)[2:-2]);
      # z_total_added = np.insert(
      #     np.append(bin_cen3, bin_cen3[-1] + bin_width), 0, bin_cen3[0] - bin_width)
      # n_total_added = np.insert(np.append(ns3, 0), 0, 0) - n_bkg
      # indi = np.where(n_total_added == np.max(n_total_added))[0][0]
      # z_fit = z_total_added[[indi - 1, indi, indi + 1]]; print 'z_fit', z_fit
      # n_fit = n_total_added[[indi - 1, indi, indi + 1]]; print 'n_fit', n_fit    
      # if (n_fit[0]<0.) and (n_fit[2]<0.):
      #     popt_gmi, pcov_gmi = curve_fit(gaussian_func, z_fit, [0,n_fit[1],0], p0=[n_fit[1],z_fit[1]])
      # else:
      #     popt_gmi, pcov_gmi = curve_fit(gaussian_func, z_fit,
      #                            n_fit, p0=[n_fit[1], z_fit[1]])
      # lum_factor2 = integrate.quad( f1d, -23.455, abs_mag(22.25, tuple(popt_gmi)[1]))[0]
      # density_factor2=integrate.quad(NFW_profile,0.001,cosmo.kpc_proper_per_arcmin(tuple(popt_gmi)[1]).value/1e3)[0]
      # signal_gmi = tuple(popt_gmi)[0] / (lum_factor2 * density_factor2)
      # z_max_fit_gmi = tuple(popt_gmi)[1]    
      # ax[2].bar(bin_cen3, ns3, width = 0.035, color='#1f77b4')#, alpha=0.5)
      # ax[2].axvline(z_max_gmi, ls='--', color='purple',
      #               label='z_max=%.3f'%z_max_gmi)
      # ax[2].axvline(z_max_fit_gmi, ls='--', color='tab:green',
      #               label='z_max_fit=%.3f'%z_max_fit_gmi)
      # ax[2].axvline(redshift,color='red',label='z:%.3f'%redshift)    
      # ax[2].plot(z_fit,n_fit+n_bkg,'o',c='tab:purple')
      # xs=np.arange(np.min(z_fit)-0.1,np.max(z_fit)+0.1,0.001)
      # ax[2].plot(xs, gaussian_func(xs, *popt_gmi) + n_bkg, c='tab:green',
      #            ls='--', label='fit: a=%.2f, mu=%.4f' % tuple(popt_gmi))
      # ax[2].legend(loc='best')
      # ax[2].set_xlabel('z')
      # ax[2].set_xlim(0.1,0.7)
      # if np.max(n_total)<30:
      #     ax[2].set_ylim(0,30)
      # ax[3].plot(zs_gmi,ntot_gmi,'.')
      # ax[3].set_xlabel('z')
      # ax[3].set_xlim(0.1,0.7)
      # ax[3].axvline(z_max_fit_gmi,ls='--',color='purple',label='z_max_fit:%.2f'%z_max_fit_gmi)
      # ax[3].axvline(redshift,color='red',label='z:%.2f'%redshift)
      # if np.max(ntot_gmi)<70:
      #     ax[3].set_ylim(0,70)
      # ntot_gmi_max=np.max(ntot_gmi)
      # zs_gmi_max=zs_gmi[np.argmax(ntot_gmi)]
      # ax[3].axvline(zs_gmi_max,ls='--',color='pink',label='zs_gmi_max:%.2f'%zs_gmi_max)
      # plt.tight_layout(rect=[0, 0., 1, 0.98])
      # plt.savefig('/Users/taweewat/Documents/red_sequence/pisco_color_plots/redsq_gmi_%s_all_%.3f_%s_tilted.png' %
      #             (mode, signal_gmi, field), dpi=120)
      # plt.close(fig)
  
      # transparent=False
      if transparent:
          fig,ax=plt.subplots(figsize=(7,4))
          ax.bar(bin_cen2, ns2, width=0.035, color='#1f77b4') #widht = 0.025
          ax.bar(bin_cen1, ns1, width = 0.035, color='#ff7f0e') #width = 0.025
          ax.axvline(z_max,ls='--',color='purple',label='z_max:%.2f'%z_max)
          ax.set_xlabel('z')
          ax.set_xlim(0.1,0.7)
          if np.max(n_total)<30:
              ax.set_ylim(0,30)
          for axp in ax.spines:
              ax.spines[axp].set_color('white') 
          ax.xaxis.label.set_color('white')
          ax.yaxis.label.set_color('white')
          ax.tick_params(axis='x', colors='white')
          ax.tick_params(axis='y', colors='white')
          purge('/Users/taweewat/Documents/red_sequence/pisco_color_plots/',
              'redsq_transparent_%.3f_%s_tilted.png' % (signal_final,field)) 
          plt.savefig('/Users/taweewat/Documents/red_sequence/pisco_color_plots/redsq_transparent_%.3f_%s_tilted.png' % (signal_final,field), dpi=120, transparent=True)
          plt.close(fig)
  
      # red_dir='/Users/taweewat/Documents/red_sequence/'
      # rich_filename = 'all_richness_%s.csv'%extra_name
      # if not os.path.isfile(red_dir + rich_filename):
      #     os.system("cp %s %s"%(red_dir+'all_richness_gremove_lum_silk_zf2.5.csv',red_dir+rich_filename))
      #     df_richness=pd.read_csv(red_dir+rich_filename)
      #     df_richness[['Nmax','Nbkg_mean','Nbkg_median','zmax','amp','zmax_fit','gremove','lum_factor','density_factor']]=np.nan
      #     df_richness.to_csv(red_dir+rich_filename)
  
      # df_richness=pd.read_csv(red_dir+rich_filename)
      # df_richness=df_richness.copy()
      # df_richness.loc[df_richness['name'] == field, 'Nmax'] = np.max(n_total)
      # df_richness.loc[df_richness['name'] == field, 'Nbkg_mean'] = np.mean(n_total)
      # df_richness.loc[df_richness['name'] == field, 'Nbkg_median'] = np.median(n_total)
      # df_richness.loc[df_richness['name'] == field, 'zmax'] = z_max
      # df_richness.loc[df_richness['name'] == field, 'amp'] = signal_final
      # df_richness.loc[df_richness['name'] == field, 'zmax_fit'] = z_max_fit
      # df_richness.loc[df_richness['name'] == field, 'gremove'] = nog
      # df_richness.loc[df_richness['name'] == field, 'lum_factor'] = lum_factor
      # df_richness.loc[df_richness['name'] == field, 'density_factor'] = density_factor
      # df_richness.to_csv(red_dir+rich_filename,index=0)
  
      red_dir='/Users/taweewat/Documents/red_sequence/'
      rich_filename = 'all_richness_%s.csv'%extra_name
      if not os.path.isfile(red_dir + rich_filename):
          df_richness=pd.DataFrame(columns=['name','Nmax','Nbkg_mean','Nbkg_median','zmax','amp','zmax_fit','gremove','lum_factor','density_factor'])
          df_richness.to_csv(red_dir+rich_filename)
      
      df_richness=pd.read_csv(red_dir+rich_filename,index_col=0)
      dic={'name':field, 'Nmax':np.max(n_total), 'Nbkg_mean':np.mean(n_total), 'Nbkg_median':np.median(n_total), 'zmax':z_max,\
      'amp':signal_final, 'zmax_fit':z_max_fit, 'gremove':nog, 'lum_factor':lum_factor, 'density_factor':density_factor}
      if field in df_richness['name'].values: 
          df_richness=df_richness[df_richness['name']!=field]
      df_richness=df_richness.append(pd.Series(dic),ignore_index=True).copy()
      df_richness.to_csv(red_dir+rich_filename)
  
      # get member redshfit in the figure
      if img_redshift:
          image_redshift(field,signal,tuple(popt)[1],mode)
      # get total images with red-sequence
      if img_filp:
          image_flip(field,signal,tuple(popt)[1],mode)
  
      if colorerr:
          return z_total, n_total, n_total_cerr
      else:
          return z_total, n_total, n_total_dup
  
  
  def pisco_combine_imgs(fields, mode='psf', mode2mass=''):
      dir1='/Users/taweewat/Documents/red_sequence/pisco_color_plots/psf_est/'
      dir2='/Users/taweewat/Documents/red_sequence/pisco_color_plots/'
      dir3='/Users/taweewat/Documents/red_sequence/pisco_color_plots/'
      dirout='/Users/taweewat/Documents/red_sequence/pisco_all/'
  
      myReg = re.compile(r'(redsq_richg_%s_all_.*%s.*png)' % (mode, field))
      myReg2=re.compile(r'(\d{1,3}\.\d{1,3})')
      names=[]
      for text in os.listdir(dir3):
          if myReg.search(text) != None:
              names.append(myReg.search(text).group())
  
      if names==[]:
          print 'no files', field
  
      signal=myReg2.search(names[0]).group()
  
      img1=dir1+'psf_est3_'+field+'_i.png'
      img15='/Users/taweewat/Documents/red_sequence/pisco_color_plots/uncer_%s_%s.png'%(mode,field)
      # img2=dir2+'star_galaxy_sep_12_all'+field+'.png'
      img2='/Users/taweewat/Documents/red_sequence/pisco_image_redshift/img_redshift_%s_%.3f_%s.png' %(mode,float(signal),field)
      print img2
      img3=dir3+names[0]
  
      images_list=[img1, img2, img3, img15]
      imgs=[]
      try:
          imgs = [ Image_PIL.open(i) for i in images_list ]
      except:
          print 'no image file', field
      mw = imgs[2].width/2
      h = imgs[0].height+imgs[1].height/1+imgs[2].height/2
      result = Image_PIL.new("RGBA", (mw, h))
      y,index=0,0
      for i in imgs:
          if index<3:
              if (index==2):# or (index==1):
                  i=i.resize((i.width/2,i.height/2))
              result.paste(i, (0, y))
              y += i.size[1]
              index+=1
          elif index==3: 
              i=i.resize((i.width/2,i.height/2))
              result.paste(i, (imgs[0].width,0))
      result.save(dirout + 'all_combine_%s_%s_%s_%s_%s.png' %
                  (extra_name, mode2mass, myReg2.search(names[0]).group(), mode, field)) 
  
  def purge(dir, pattern):
      for f in os.listdir(dir):
          if re.search(pattern, f):
              print 'remove', f
              os.remove(os.path.join(dir, f))
  
  def image_redshift(field,signal,redshift,mode):
      df_total=pd.read_csv('/Users/taweewat/Documents/pisco_code/slr_output/galaxy_%s_final_%s.csv'%(mode,field),index_col=0)
      df_star=pd.read_csv('/Users/taweewat/Documents/pisco_code/slr_output/star_psf_total_%s.csv'%field,index_col=0)
      # df_star=df_star[df_star['SG']>0.95]
      hdu=fits.open('/Users/taweewat/Documents/pisco_code/final/coadd_c%s_i.fits'%field)
      img=hdu[0].data.astype(float)
      img -= np.median(img)    
      df_total['redshift_m']=df_total.apply(lambda row: ex.redshift_f(row), axis=1)
      def size_f(row):
          if not np.isnan(row['w_gmr']):
              size=row['w_gmr']
          if not np.isnan(row['w_rmi']):
              size=row['w_rmi']
          if np.isnan(row['w_rmi']) and np.isnan(row['w_gmr']):
              size=0
          return size
      df_total['size_m']=df_total.apply(lambda row: size_f(row), axis=1)
  
      df_total=df_total[df_total['redshift_m'] > 0].copy()
  
      norm = matplotlib.colors.Normalize(vmin=0,vmax=500)
      c_m = matplotlib.cm.Greys_r
      s_m = matplotlib.cm.ScalarMappable(cmap=c_m, norm=norm)
      s_m.set_array([])
  
      normalize = matplotlib.colors.Normalize(vmin=0.1, vmax=0.7)
  
  
      fig, (a0, a1) = plt.subplots(1,2, figsize=(30,18), gridspec_kw = {'width_ratios':[0.8, 1]})
      # a0.imshow(img, cmap=c_m, norm=norm, origin='lower')
      # a0.scatter(df_star['XWIN_IMAGE_i'].values,df_star['YWIN_IMAGE_i'].values,s=100, marker='*', facecolors='none', edgecolors='yellow', label='star')
      # df1i=df_total[df_total['w_rmi']>0.1]
      # df2i=df_total[df_total['w_rmi']<=0.1]
      # # a0.scatter(df1i['XWIN_IMAGE_i'].values,df1i['YWIN_IMAGE_i'].values,s=100, facecolors='none', edgecolors='blue')
      # a0.scatter(df1i['XWIN_IMAGE_i'].values, df1i['YWIN_IMAGE_i'].values, s=100, c=df1i['size_m'].values, cmap='RdYlGn')
      # a0.scatter(df2i['XWIN_IMAGE_i'].values,df2i['YWIN_IMAGE_i'].values, s=100, facecolors='none', edgecolors='white')
      # a0.set_xlim(0,1600)
      # a0.set_ylim(0, 2250)
      try:
          img2 = mpimg.imread('/Users/taweewat/Documents/pisco_code/Chips_images/aplpy4_%s_img4.jpeg' % field)
      except:
          img2 = mpimg.imread('/Users/taweewat/Documents/pisco_code/Chips_images/aplpy4_%s_img.jpeg' % field)
      imgplot = a0.imshow(img2)
      a0.axis('off')
      a0.annotate('Redshift: %.3f\nRichness: %.2f' %
                  (redshift, signal), xy=(150, 100), color='white')
  
      a1.imshow(img, cmap=c_m, norm=norm, origin='lower')
      a1.scatter(df_star['XWIN_IMAGE_i'].values,df_star['YWIN_IMAGE_i'].values, s=300,edgecolor='orange', facecolor='none',lw=3)
      #,s=100, marker='*', facecolors='none', edgecolors='yellow', label='star')
      axi = a1.scatter(df_total['XWIN_IMAGE_i'].values, df_total['YWIN_IMAGE_i'].values,
                       s=(df_total['size_m'].values * 200)+30, c=df_total['redshift_m'].values, cmap='tab20b', norm=normalize)
      plt.colorbar(axi)  # df_total['size_m'].values*300
      a1.set_xlim(0, 1600)
      a1.set_ylim(0, 2250)
      plt.tight_layout()
  
      left, bottom, width, height = [0.05, 0.24, 0.3, 0.2]
      ax2 = fig.add_axes([left, bottom, width, height])
      ax2.imshow(mpimg.imread(
          '/Users/taweewat/Documents/red_sequence/pisco_color_plots/redsq_transparent_%.3f_%s_tilted.png' % (signal, field)))
      ax2.axes.get_xaxis().set_visible(False)
      ax2.axes.get_yaxis().set_visible(False)
      ax2.axis('off')
  
      plt.savefig('/Users/taweewat/Documents/red_sequence/pisco_image_redshift/img_redshift_%s_%.3f_%s.png' %
                  (mode,signal,field), dpi=50)
      plt.close(fig)
  
  
  def image_flip(field, signal, redshift, mode):
      img = mpimg.imread(
          '/Users/taweewat/Documents/pisco_code/Chips_images/aplpy4_%s_img.jpeg' % field)
      fig, ax = plt.subplots(figsize=(7, 7))
      imgplot = ax.imshow(img)
      ax.axis('off')
      ax.annotate('Redshift: %.3f\nRichness: %.2f' %
                  (redshift, signal), xy=(150, 100), color='white')
  
      left, bottom, width, height = [0.2, 0.18, 0.3, 0.2]
      ax2 = fig.add_axes([left, bottom, width, height])
      ax2.imshow(mpimg.imread(
          '/Users/taweewat/Documents/red_sequence/pisco_color_plots/redsq_transparent_%.3f_%s_tilted.png' % (signal, field)))
      ax2.axes.get_xaxis().set_visible(False)
      ax2.axes.get_yaxis().set_visible(False)
      ax2.axis('off')
  
      # plt.tight_layout()
      plt.savefig('/Users/taweewat/Documents/red_sequence/pisco_image_redshift/image_flip_%s_%.3f_%s.png' %
                  (mode, signal, field), dpi=200)
      plt.close(fig)
  
  
  if __name__ == "__main__":
      """
      execute:
      python pisco_pipeline/pisco_photometry_all.py CHIPS111 psf slr
      #updated version with no2mass option for no more comparison with known 2mass stars
      python pisco_pipeline/pisco_photometry_all.py CHIPS111 psf allslr no2mass
      """
      print 'Number of arguments:', len(sys.argv), 'arguments.'
      print 'Argument List:', str(sys.argv)
  
      field = str(sys.argv[1])
      mode = str(sys.argv[2])   #aper, psf, auto, hybrid
      all_argv=sys.argv[3:]  #allslr, slr, noslr
      if (all_argv[0]=='allslr') | (all_argv[0]=='slr'):
          slr=str(all_argv[0])
          slr_param=True
      elif all_argv[0]=='noslr':
          slr='no_slr'
          slr_param=False
      if all_argv[1]=='2mass':
          mode2mass=''
      elif all_argv[1]=='no2mass':
          mode2mass='_no2mass'
  
      home='/Users/taweewat/Documents/pisco_code/' #09, 171208
      # dirs=['ut170103/','ut170104/','ut170619/','ut170621/','ut170624/','ut171208/','ut171209/','ut171212/']
      dirs=['ut190412','ut190413']
      # 'ut171208/', 'ut171209/','ut171212/', 'ut170621/', 'ut170624/'
      # dirs = ['ut170621/','ut170624/']
      # dirs = ['ut170619/']
      # dirs = ['ut170103/']
      names=[]
      myReg=re.compile(r'(CHIPS\d{4}[+-]\d{4})|(Field\d{3})')
      for di in dirs:
          dir=home+di
          for text in os.listdir(dir):
              if myReg.search(text) != None:
                  names.append(myReg.search(text).group())
      all_fields=list(set(names))
      # print all_fields
  
      infile = open('/Users/taweewat/Documents/xray_project/code_github/allremove_chips.txt', 'r')
      exception = [i.strip() for i in infile.readlines()]
  
      all_fields_cut = all_fields[:]
      all_fields_cut = ['SDSS603','SDSS501','SDSS123']
      print all_fields_cut
      # all_fields_cut = ['CHIPS1422-2728']
  
  
      notgoflag=True
      z_total_all,n_total_all,n_total_dup_all=[],[],[]
      for index, field in enumerate(all_fields_cut):
          print field, '%i/%i' % (index, len(all_fields_cut))
  
          # if field == 'CHIPS0122-2646':
          #     notgoflag = False; continue
          # if notgoflag:
          #     continue 
          if field in exception:
              continue
          if field in ['CHIPS1933-1511']:
              continue
  
          if slr=='allslr':
              print 'allslr'
              pisco_photometry_v4(field)
          elif slr=='slr':
              # star_galaxy_bleem(field)
              pisco_cut_frame(field)
              pisco_photometry_psf_v4(field, mode=mode, mode2mass=mode2mass, slr=slr_param)
  
              purge('/Users/taweewat/Documents/red_sequence/pisco_color_plots/'\
                  ,r'(redsq_%s_all_.*%s.*png)'%(mode,field))
              
              # pisco_tilt_resequence(field, mode='psf', mode2mass='')
              z_total=pisco_tilt_resequence(field, mode=mode, mode2mass=mode2mass)
              # z_total_all.append(z_total)
              # n_total_all.append(n_total)
              # n_total_dup_all.append(n_total_dup)
              # pisco_combine_imgs(field, mode=mode, mode2mass=mode2mass)
              # pickle.dump( [z_total_all,n_total_all,n_total_dup_all], open( "pickle_all_richness_%s.pickle"%extra_name, "wb" ) )
              # print 'save pickle fie at', "pickle_all_richness_%s.pickle" % extra_name
          elif slr == 'no_slr':
              pisco_cut_frame(field)
              pisco_photometry_psf_v4(field, mode=mode, mode2mass=mode2mass, slr=slr_param)
              purge('/Users/taweewat/Documents/red_sequence/pisco_color_plots/'\
                  ,r'(redsq_%s_all_.*%s.*png)'%(mode,field))
              z_total=pisco_tilt_resequence(field, mode=mode, mode2mass=mode2mass)
              # z_total_all.append(z_total)
              # n_total_all.append(n_total)
              # n_total_dup_all.append(n_total_dup)
              # pisco_combine_imgs(field, mode=mode, mode2mass=mode2mass)
              # pickle.dump( [z_total_all,n_total_all,n_total_dup_all], open( "pickle_all_richness_%s.pickle"%extra_name, "wb" ) )
              # print 'save pickle fie at', "pickle_all_richness_%s.pickle" % extra_name
      
          purge('final', "proj_coadd_c%s_.*\.fits" % field)
          purge('.', "proto_psf_%s_.*\.fits" % field)
          purge('.', "samp_psf_%s_.*\.fits" % field)
          purge('.', "resi_psf_%s_.*\.fits" % field)
          purge('.', "snap_psf_%s_.*\.fits" % field)
          purge('.', "chi_psf_%s_.*\.fits" % field)
          # purge('psfex_output', "psf_%s_.*\.fits" % field)
          # purge('slr_output', "a_psf_%s_.*\.fits" % field)
          purge('final', "coadd_c%s_sq_.*\.fits" % field)
  
          
  