/mccountscripts.py Secret
Created
June 13, 2016 22:55
mccountscripts example
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/usr/bin/env python | |
| from __future__ import division,with_statement | |
| from numpy import * | |
| from numpy.random import * | |
| from matplotlib.pyplot import * | |
| import matplotlib.pyplot as plt | |
| from astropysics import phot,plotting | |
| from astropysics.utils import io | |
| import numpy as np | |
| import pyfits | |
| import mcpsscripts | |
| #All RC3 mags corrected (e.g. _0) | |
| LMCappmags = {'U':0.52,'B':0.57,'V':0.13} | |
| SMCappmags = {'U':2.05,'B':2.28,'V':1.87} | |
| M32appmags = {'U':9.12 ,'B':8.72 ,'V':7.84,'R':6.97} | |
| NGC205appmags = {'B':8.79 ,'V':7.97 } | |
| NGC185appmags = {'U':9.51,'B':9.28,'V':8.55} | |
| NGC147appmags = {'B':9.71,'V':8.93} | |
| M33appmags = {'U':5.59,'B':5.75,'V':5.28} | |
| LMCdm = 18.50 #OGLE/Araucaria | |
| SMCdm = 18.93 #?? | |
| M32dm = 24.41 | |
| NGC205dm = 24.57 | |
| NGC185dm = 24.07 | |
| NGC147dm = 24.15 | |
| M33dm = 24.70 | |
| LMCbothun = {'B-V':.56,'B':0.25} #foreground reddening corrected as in bothun, mags within 24.7 mag asec^-2 | |
| #LMCbothun['B-V'] = .61 #raw/no reddening | |
| LMCbothun['B-V'] = .61 -0.075 #schlegel extinction, LMC E(B=V)=0.075 | |
| LMCbothun['V'] = LMCbothun['B'] - LMCbothun['B-V'] | |
| LMCbothun['g-r'] = 1.02*LMCbothun['B-V']-.22 | |
| SMCbothun = {'B-V':.47,'B':2.58} #reddening corrected, mags within 24.7 mag asec^-2 | |
| #SMCbothun = {'B-V':.52} | |
| SMCbothun['V'] = LMCbothun['B'] - LMCbothun['B-V'] | |
| SMCbothun['g-r'] = 1.02*LMCbothun['B-V']-.22 | |
| LMCvtan,LMCvr = 367.0,89.0 #\pm 18,4 -> kav? | |
| M31appmags={'U':3.67,'B':3.36 ,'V':2.68 ,'I':0.73} #src? | |
| M31dm = 24.45 | |
| #where did I get this from? seems like the color is a bit crazy | |
| MWabsmags={'U':-20.74 , 'B':-20.74 , 'V':-21.26 , 'R':-21.78 , 'I':-22.3} | |
| for d in (LMCappmags,LMCbothun,SMCappmags,SMCbothun,M32appmags,NGC205appmags, | |
| NGC185appmags,NGC147appmags,M33appmags,M31appmags,MWabsmags): | |
| phot.transform_dict_ugriz_UBVRI(d) | |
| def load_VAGC(vagc='dr6',loadcat=True, | |
| loadimg=True, | |
| loadspec=False, | |
| loadkcnone=True, | |
| loadkcnear=True, | |
| loadcollnone=True, | |
| loadcollnear=True, | |
| memmap=1): | |
| from os.path import join | |
| if vagc == 'dr6': | |
| vagcdir = 'NYUVAGCdr6' | |
| kcorrtype = 'petro' | |
| elif vagc == 'dr7': | |
| vagcdir = 'NYUVAGCdr7' | |
| kcorrtype = 'petro' | |
| elif vagc == 'model' or vagc=='dr7model': | |
| vagcdir = 'NYUVAGCdr7' | |
| kcorrtype = 'model' | |
| else: | |
| raise ValueError('invalid VAGC') | |
| if loadcat: | |
| catf = pyfits.open(join(vagcdir,'object_catalog.fits'),memmap=memmap) | |
| catd = catf[1].data | |
| else: | |
| catd = None | |
| if loadimg: | |
| imsf = pyfits.open(join(vagcdir,'object_sdss_imaging.fits'),memmap=memmap) | |
| imd = imsf[1].data | |
| else: | |
| imd = None | |
| if loadspec: | |
| specsf = pyfits.open(join(vagcdir,'object_sdss_spectro.fits'),memmap=memmap) | |
| specd = specsf[1].data | |
| else: | |
| specd = None | |
| if loadkcnone: | |
| kcorrf = pyfits.open(join(vagcdir,'kcorrect.none.'+kcorrtype+'.z0.00.fits'),memmap=memmap) | |
| kcorrd = kcorrf[1].data | |
| else: | |
| kcorrd = None | |
| if loadkcnear: | |
| kcorrnearf = pyfits.open(join(vagcdir,'kcorrect.nearest.'+kcorrtype+'.z0.00.fits'),memmap=memmap) | |
| kcorrneard = kcorrnearf[1].data | |
| else: | |
| kcorrneard = None | |
| if loadcollnone: | |
| collnonef = pyfits.open(join(vagcdir,'collisions.none.fits'),memmap=memmap) | |
| collnonef = collnonef[1].data | |
| else: | |
| collnonef = None | |
| if loadcollnear: | |
| collnearf = pyfits.open(join(vagcdir,'collisions.nearest.fits'),memmap=memmap) | |
| collneard = collnearf[1].data | |
| else: | |
| collneard = None | |
| return {'catalog':catd, | |
| 'imaging':imd, | |
| 'spectro':specd, | |
| 'kcorr':kcorrd, | |
| 'collnearest':collneard, | |
| 'collnone':collnonef, | |
| 'dr':vagc} | |
| def load_cuts_dr6(): | |
| """ | |
| "hosts" are isolated galaxies w/ r < -20 | |
| "sats" are - 17.5 > r > -20 | |
| the volume is dictated by -17.5 completeness | |
| the output dictionary has entries for: | |
| * 'z500': hosts with satellites w/i 500 km/s | |
| * 'nozcut': all hosts with satellites | |
| * 'noz': hosts with no redshift (for the satellite?) | |
| * 'uniquenoz': hosts with no redshift and no other neighbor in noz or nozcut | |
| each of these is a recarray with `i` indicating the VAGC index of the host, | |
| `nn` indicating the number of neighbors, `nd` indicating the nearest one's | |
| distance, and `ni` giving the index of the nearest sat | |
| also, arrays of indecies of the following groups: | |
| *'pothosts': potential hosts | |
| *'potsats': potential sats with redshift | |
| *'potsatsnoz': potential sats without redshift | |
| """ | |
| names,types='i,nn,nd,ni','i,i,f,i' | |
| drdir = 'dr6dat' | |
| z500 = rec.fromarrays(loadtxt(drdir+'/hosts_M20_Dn250_M17.5_z500.dat').T,names=names,formats=types) | |
| nozcut = rec.fromarrays(loadtxt(drdir+'/hosts_M20_Dn250_M17.5_nozcut.dat').T,names=names,formats=types) | |
| noz = rec.fromarrays(loadtxt(drdir+'/hosts_M20_Dn250_M17.5_noz.dat').T,names=names,formats=types) | |
| uniquenoz = noz[~np.in1d(noz.i,nozcut.i)] | |
| #find those where the noz is closer than the z gal | |
| closernoz = noz[np.in1d(noz.i,nozcut.i)] | |
| nozcutclosei = searchsorted(np.sort(nozcut.i),closernoz.i) | |
| closernoz = closernoz[closernoz.nd < nozcut[nozcutclosei].nd] | |
| del names,types,nozcutclosei | |
| return locals() | |
| def load_pots_dr6(): | |
| drdir = 'dr6dat' | |
| pothosts=loadtxt(drdir+'/ph_M20_z17.5_Dn300.dat',dtype=int) | |
| potsats=loadtxt(drdir+'/pn_z_M17.5.dat',dtype=int) | |
| potsatsnoz=loadtxt(drdir+'/pn_noz_M17.5.dat',dtype=int) | |
| return locals() | |
| def load_cuts_dr7(): | |
| """ | |
| "hosts" are isolated galaxies w/ r < -20 | |
| "sats" are - 17.5 > r > -20 | |
| the volume is dictated by -17.5 completeness | |
| the output dictionary has entries for: | |
| * 'z500': hosts with satellites w/i 500 km/s | |
| * 'nozcut': hosts with satellites w/redshifts | |
| * 'zbad': hosts with satellites where dv>500 km/s | |
| * 'noz': hosts w/ satellites w/o redshifts | |
| * 'uniquenoz': noz objects where there is no sat w/z | |
| * 'closernoz':noz objects where the noz satellite is closer than a z sat | |
| each of these is a recarray with `i` indicating the VAGC index of the host, | |
| `nn` indicating the number of neighbors, `nd` indicating the nearest one's | |
| distance, and `ni` giving the index of the nearest sat | |
| also, arrays of indecies of the following groups: | |
| *'pothosts': potential hosts | |
| *'potsats': potential sats with redshift | |
| *'potsatsnoz': potential sats without redshift | |
| """ | |
| names,types='i,nn,nd,ni','i,i,f,i' | |
| drdir = 'dr7dat' | |
| z500=rec.fromarrays(loadtxt(drdir+'/hosts_M20_Dn250_M17.5_z500.dat').T,names=names,formats=types) | |
| nozcut=rec.fromarrays(loadtxt(drdir+'/hosts_M20_Dn250_M17.5_nozcut.dat').T,names=names,formats=types) | |
| noz=rec.fromarrays(loadtxt(drdir+'/hosts_M20_Dn250_M17.5_noz.dat').T,names=names,formats=types) | |
| zbad = nozcut[~np.in1d(nozcut.ni,z500.ni)] | |
| uniquenoz = noz[~np.in1d(noz.i,nozcut.i)] | |
| #find those where the noz is closer than the z gal | |
| closernoz = noz[np.in1d(noz.i,nozcut.i)] | |
| nozcutclosei = searchsorted(np.sort(nozcut.i),closernoz.i) | |
| closernoz = closernoz[closernoz.nd < nozcut[nozcutclosei].nd] | |
| del names,types,nozcutclosei | |
| return locals() | |
| def load_pots_dr7(): | |
| drdir = 'dr7dat' | |
| pothosts=loadtxt(drdir+'/ph_M20_z17.5.dat',dtype=int) | |
| potsats=loadtxt(drdir+'/pn_z_M17.5.dat',dtype=int) | |
| potsatsnoz=loadtxt(drdir+'/pn_noz_M17.5.dat',dtype=int) | |
| return locals() | |
| zentnerformat = [ | |
| 'id','int',# 1 original ID number from Zentner's output | |
| 'n_neighbors','int',# 2 within DLIM=0.7 Mpc/h, MASSIVE neighbors only (***above VcMIN***) | |
| 'halo_vin','f8',# 3 Vin correcting central galaxy for satellites above cut only | |
| 'z_obs','f8', | |
| 'quadrant','int',# 5 quadrant of the galaxy | |
| 'ra','f8', | |
| 'dec','f8', | |
| 'm_host','f8', | |
| 'minsep','f8', # 9 COMOVING separation of closest neighbor w/in 500 km/s ***above VcMIN***, in Mpc (NOT Mpc/h) | |
| 'vel_of_minsep','f8',#10 velocity difference to that neighbor***above VcMIN*** | |
| 'sorted_distances4','f8',# 11 fifth nearest neighbor distance, Mpc (h=0.7)***above VcMIN*** | |
| 'n_sat','int', # 12 number of satellites in host ***above VcMIN*** | |
| 'sorted_distances1','f8', # 13 second nearest neighbor distance | |
| 'sorted_distances2','f8', # 14 third nearest neighbor distance | |
| 'a_accreted','f8', # 15 scale factor at time accreted (0 for host) | |
| 't_encounter_last','f8', # 16 time since last close pass (if any) | |
| 'b_encounter_last','f8', # 17 b of last close pass, already in kpc/h | |
| 'v_encounter_last','int',# 18 v of last close pass | |
| 'ratio_encounter_last','f8',# 19 Vmax ratio of last close pass (Vmax ratios,not Vin) | |
| 'host_mass_encounter_last','f8',# 20 Vmax mass of host at last close pass | |
| 'n_encounter','int',# 21 number of previous close passes | |
| 'same_as_neighbor','int', # 22 last close pass same galaxy as nearest neighbor? 0=no/1=yes/9=no passes; | |
| 'ltidal','f8', # 23 log of tidal force in arbitrary scaling | |
| 'truesep','f8', # 24 true 3-D distance to apparent nearest neighbor *** above VcMin ***, MPC (NOT MPC/h) | |
| 'halo_orig','f8', # 25 original Vin without subtraction of satellites | |
| 'my_mass','f8', # 26 my final mass | |
| 'vin_host','f8', # 27 my host v_in | |
| 'avg_hostmass_cylinder','f8', # 28 average host mass in cylinder | |
| 'max_hostmass_neighbors','f8', # 29 maximum host mass in cylinder | |
| 'n_neighbors_1','int', # 30 within 1 Mpc/h ***above VcMIN*** | |
| 'n_neighbors_3','int', # 31 within 3 Mpc/h ***above VcMIN*** | |
| 'n_neighbors_6','int', # 32 within 6 Mpc/h ***above VcMIN*** | |
| 'id_host','int', # 33 id of host galaxy | |
| 'icount_host','int', # 34 icount number of host | |
| 'halo_vnow','f8', # 35 Vnow of halo | |
| 'icount','int', # 36 icount value of this galaxy | |
| 'icount_of_minsep','int', # 37 icount value of nearest apparent | |
| 'n_sat_inc_minor','int' # 38 number of satellites in host including minors | |
| ] | |
| def load_betsy_zentner(fn=1201,nisolimit=1,disolimit=250,minorvin=105,majorvin=196, | |
| minorrad=250,deltav=500): | |
| """ | |
| :param fn: filename to load | |
| :param nisolimit: number isolation limit (number of neighbors <= this is isolated) | |
| :param disolimit: distance isolation limit in kpc/h (nn >= this is isolated) | |
| :param minorvin: lower-limit on Vin for minors (km/s) | |
| :param majorvin: lower limit on Vin for major/minor boundary (km/s) | |
| :param minorrad: maximum distance of minor (kpc/h) | |
| :param deltav: upper limit of v_major-v_minor for lmclike (km/s) | |
| Note that h is taken to be .7 for zentner sim | |
| """ | |
| from astropysics.coords import cosmo_z_to_dist | |
| from astropysics.constants import asecperrad | |
| dtd = {'names':zentnerformat[::2],'formats':zentnerformat[1::2]} | |
| h = .7 | |
| if isinstance(fn,int): | |
| fn = 'out%i_700Kpc_cosmo_105_minor.dat'%fn | |
| #load raw zentner | |
| arr1 = np.loadtxt(fn,dtype=dtype(dtd)) | |
| #add ancillary data types | |
| dtd['names'].append('minsepkpch') #minsep in kpc/h | |
| dtd['formats'].append('float') | |
| dtd['names'].append('truesepkpch') #truesep in kpc/h | |
| dtd['formats'].append('float') | |
| dtd['names'].append('deltav') #km/s between minsep and self with sign information | |
| dtd['formats'].append('float') | |
| dtd['names'].append('deltav_host') #km/s between host and self with sign information | |
| dtd['formats'].append('float') | |
| dtd['names'].append('t_pass') #Gyr since last encounter | |
| dtd['formats'].append('float') | |
| dtd['names'].append('actual_sat') | |
| dtd['formats'].append('bool') | |
| dtd['names'].append('isolated') #meaningless unless also a major | |
| dtd['formats'].append('bool') | |
| dtd['names'].append('minor') #<mass cut | |
| dtd['formats'].append('bool') | |
| dtd['names'].append('major') #>mass cut | |
| dtd['formats'].append('bool') | |
| dtd['names'].append('minsep_major_iso') #minsep is a major isolated object | |
| dtd['formats'].append('bool') | |
| dtd['names'].append('minsep_not_major_iso') #minsep is not a major isolated object | |
| dtd['formats'].append('bool') | |
| dtd['names'].append('lmc_analog') #analog: sat of isolated host, mass and deltav cut | |
| dtd['formats'].append('bool') | |
| dtd['names'].append('lmc_analog_noz') #analog except no deltav cut: sat of isolated host, mass cut | |
| dtd['formats'].append('bool') | |
| dtd['names'].append('lmc_analog_noiso') #analog except no isolation for minsep | |
| dtd['formats'].append('bool') | |
| dtd['names'].append('nearest_lmc_analog') #nearest sat and analog | |
| dtd['formats'].append('bool') | |
| dtd['names'].append('nearest_lmc_analog_noz') #nearest sat and analog w/ noz | |
| dtd['formats'].append('bool') | |
| dtd['names'].append('nearest_lmc_analog_noiso') #nearest sat and analog w/no isolation | |
| dtd['formats'].append('bool') | |
| dtd['names'].append('analog_host') #hosts of an LMC analog | |
| dtd['formats'].append('bool') | |
| dtd['names'].append('angminsep') #angular seperation in arcsec | |
| dtd['formats'].append('float') | |
| dtd['names'].append('fibercoll') #if True, angminsep<55" | |
| dtd['formats'].append('bool') | |
| arr2 = np.empty(arr1.size,dtype=dtype(dtd)) | |
| for ns in zentnerformat[::2]: | |
| arr2[ns] = arr1[ns] | |
| minsepkpch = arr2['minsepkpch'] = arr1['minsep']*1000*h # -> kpc/h | |
| arr2['truesepkpch'] = arr1['truesep']*1000*h # -> kpc/h | |
| arr2['actual_sat'] = arr1['icount_host'] == arr1['icount_of_minsep'] | |
| arr2['deltav'] = 3e5*(arr1['z_obs']-arr1[arr1['icount_of_minsep']]['z_obs']) | |
| arr2['deltav_host'] = 3e5*(arr1['z_obs']-arr1[arr1['icount_host']]['z_obs']) | |
| arr2['t_pass'] = 13.48 - arr1['t_encounter_last'] | |
| arr2['isolated'] = ((arr1['n_neighbors'] <= nisolimit) &\ | |
| (minsepkpch >= (disolimit)) | |
| ) | (arr1['n_neighbors']==0) | |
| arr2['minor'] = (arr1['halo_vin']>minorvin)&(arr1['halo_vin']<majorvin) | |
| arr2['major'] = arr1['halo_vin'] >= majorvin | |
| arr2['minsep_major_iso'] = arr2['major'][arr2['icount_of_minsep']] &\ | |
| arr2['isolated'][arr1['icount_of_minsep']] #&\ | |
| #(arr1['id_host']!=arr1['id']) | |
| #TODO:replace this with check that you're not an apparent host of some kind? | |
| arr2['minsep_not_major_iso'] = arr2['major'][arr2['icount_of_minsep']] &\ | |
| (~arr2['isolated'][arr1['icount_of_minsep']]) | |
| arr2['lmc_analog'] = arr2['minor'] & arr2['minsep_major_iso'] &\ | |
| (minsepkpch<(minorrad)) &\ | |
| (arr1['vel_of_minsep']<deltav) | |
| arr2['lmc_analog_noz'] = arr2['minor'] & arr2['minsep_major_iso'] &\ | |
| (minsepkpch<(minorrad)) | |
| arr2['lmc_analog_noiso'] = arr2['minor'] & arr2['minsep_not_major_iso'] &\ | |
| (minsepkpch<(minorrad)) &\ | |
| (arr1['vel_of_minsep']<deltav) | |
| #figure out which of the lmc_analogs are nearest | |
| arrlana = arr2[arr2['lmc_analog']] | |
| lanahosts = arr2[np.unique(arrlana['icount_host'])] | |
| nsatsis = [] | |
| for host in lanahosts: | |
| #TODO: un-for if the data size gets larger | |
| sats = arrlana[host['icount']==arrlana['icount_host']] | |
| nsatsis.append(sats[sats['minsep']==sats['minsep'].min()][0]['icount']) | |
| hostnn = np.zeros(arr2['lmc_analog'].size,dtype=bool) | |
| hostnn[nsatsis] = True | |
| arr2['nearest_lmc_analog'] = arr2['lmc_analog']&hostnn | |
| #now repeat for the noz and noiso sets | |
| #noz | |
| arrlana = arr2[arr2['lmc_analog_noz']] | |
| lanahosts = arr2[np.unique(arrlana['icount_host'])] | |
| nsatsis = [] | |
| for host in lanahosts: | |
| #TODO: un-for if the data size gets larger | |
| sats = arrlana[host['icount']==arrlana['icount_host']] | |
| nsatsis.append(sats[sats['minsep']==sats['minsep'].min()][0]['icount']) | |
| nnnoz = np.zeros(arr2['lmc_analog_noz'].size,dtype=bool) | |
| nnnoz[nsatsis] = True | |
| arr2['nearest_lmc_analog_noz'] = arr2['lmc_analog_noz']&nnnoz | |
| #niso | |
| arrlana = arr2[arr2['lmc_analog_noiso']] | |
| lanahosts = arr2[np.unique(arrlana['icount_host'])] | |
| nsatsis = [] | |
| for host in lanahosts: | |
| #TODO: un-for if the data size gets larger | |
| sats = arrlana[host['icount']==arrlana['icount_host']] | |
| nsatsis.append(sats[sats['minsep']==sats['minsep'].min()][0]['icount']) | |
| nnnoiso = np.zeros(arr2['lmc_analog_noiso'].size,dtype=bool) | |
| nnnoiso[nsatsis] = True | |
| arr2['nearest_lmc_analog_noiso'] = arr2['lmc_analog_noiso']&nnnoiso | |
| ahostsi = np.unique(arr1['icount_host'][arr2['lmc_analog']]) | |
| arr2['analog_host'] = np.zeros_like(arr2['lmc_analog']) | |
| arr2['analog_host'][ahostsi] = True | |
| arr2['angminsep'] = asecperrad*arr1['minsep']/cosmo_z_to_dist(arr1['z_obs']) #comoving Mpc | |
| arr2['fibercoll'] = arr2['angminsep']<55 | |
| return arr2.view(np.recarray) | |
| class MSIICone(object): | |
| def __init__(self,fnorarr,deltavcut=500,niso=1): | |
| if isinstance(fnorarr,basestring): | |
| arrs = np.load(fnorarr) | |
| else: | |
| arrs = fnorarr | |
| self.arrs = arrs | |
| self.deltavcut = deltavcut | |
| self.niso = niso | |
| self._make_derived() | |
| self._make_masks() | |
| @property | |
| def n(self): | |
| return len(self.arrs) | |
| @property | |
| def fields(self): | |
| fds = list(self.arrs.dtype.names) | |
| fds.extend(self.derived.keys()) | |
| return tuple(fds) | |
| def _make_derived(self): | |
| from astropysics.models import NFWModel | |
| self.derived = d = {} | |
| a = self.arrs | |
| d['rvir_h'] = NFWModel.Vmax_to_Rvir(a['h_vmax'],z=0)*.7/1e3 #kpc>Mpc/h | |
| deld = a['dist_h']-a['dist_s'] #Mpc/h | |
| d['r2dto3d'] = (a['rproj']**2+deld**2)**0.5#pri/sec 3D seperation | |
| def _make_masks(self): | |
| a = self.arrs | |
| self.masks = masks = {} | |
| masks['conez'] = a['dist_h']<=100 | |
| masks['isolated'] = (a['n_close'] == 0) & (a['n_far']<=self.niso) | |
| masks['deltav'] = np.abs(a['delta_v']) < self.deltavcut | |
| masks['analog'] = masks['conez']&masks['isolated']&masks['deltav']&(a['rproj']<.25) | |
| masks['central_host'] = a['host_type']==0 | |
| masks['truesat_fof'] = masks['central_host']&(a['same_fof']==1) | |
| masks['truesat_3d_fof'] = masks['central_host']&(a['r3d_same_fof']==1) | |
| hrvir = self.derived['rvir_h'] | |
| psec3d = self.derived['r2dto3d'] | |
| masks['wivvir'] = psec3d < hrvir | |
| masks['twodisthreed'] = ((psec3d - a['r3d']) < .005) | |
| masks['truesat_spherical'] = masks['central_host']&masks['wivvir'] | |
| # locs=((dat['n_close'][i] == 0) & (dat['n_far'][i] < 2) & (numpy.abs(dat['delta_v'][i]) < 500) & (dat['dist_h'][i] <= 100)).nonzero()[0] | |
| # locs1=((dat['n_close'][i] == 0) & (dat['n_far'][i] < 2) & (dat['dist_h'][i] <= 100)).nonzero()[0] | |
| def get_masked_arr(self,arrname,masknames=None,i=slice(None)): | |
| """ | |
| gets the specified array name with the masks applied. `masknames` should | |
| be a 'name1,name2,...' string or None | |
| """ | |
| arr = getattr(self,arrname)[i] | |
| if masknames is None: | |
| masks = slice(None) | |
| else: | |
| masks = self.masks | |
| names = [n[1:] if n.startswith('~') else n for n in masknames.split(',')] | |
| invs = [n.startswith('~') for n in masknames.split(',')] | |
| masks = [(~masks[n][i]) if inv else masks[n][i] for n,inv in zip(names,invs)] | |
| masks = reduce(np.logical_and,masks) | |
| return arr[masks] | |
| def __getattr__(self,key): | |
| if key in self.masks: | |
| return self.masks[key] | |
| elif key in self.derived: | |
| return self.derived[key] | |
| elif key in self.arrs.dtype.names: | |
| return self.arrs[key] | |
| else: | |
| return object.__getattribute__(self,key) | |
| #<--------------------Plot support functions-----------------------------------> | |
| _plotreg = {} | |
| _mainfigs = [] | |
| def plotfunc(figtype='mpl',figsize=None,tweakbbox=None,mainfig=False): | |
| """ | |
| :param figtype: | |
| Specifies the type/format of the figure. Can be any of: | |
| * 'mpl' | |
| Matplotlib figure saving to both eps and pdf. | |
| * 'mpl<ext1>,<ext2>' | |
| Matplotlib figure saving each figure to filenames that have the | |
| requested extension. | |
| * 'mpltoeps' | |
| saves as png and then converts png->eps | |
| * 'mayavi' | |
| Mayavi figure, defaults to png | |
| * 'mayavi<ext1>,<ext2>' | |
| Mayavi figure, saving to the specified extension types | |
| :param figsize: | |
| Specifies the size of the figure as a (width,height) tuple or None for default. | |
| :param tweakbbox: | |
| Specifies the bounding box to be used for eps output as a (x1,y1,x2,y2) | |
| tuple or None for the default bounding box. | |
| :param bool mainfig: | |
| If True, this is a "main" figure - e.g. one that will actually appear in | |
| a paper (used for the command-line call) | |
| """ | |
| def deco(f): | |
| fname = f.func_name | |
| if 'mpl' in figtype: | |
| if figsize is None: | |
| fsz = (8,8) | |
| else: | |
| fsz = figsize | |
| elif 'mayavi' in figtype: | |
| if figsize is None: | |
| fsz = (800,600) | |
| else: | |
| fsz = figsize | |
| else: | |
| raise ValueError('unrecognized plot type %s'%figtype) | |
| _plotreg[fname]=(f,fsz,figtype,tweakbbox) | |
| if mainfig: | |
| _mainfigs.append(fname) | |
| return f | |
| if callable(figtype): | |
| f=figtype | |
| figtype='mpl' | |
| return deco(f) | |
| else: | |
| return deco | |
| def plot_names(): | |
| """ | |
| Returns the names of all registered plots for this file, as well as the | |
| names of just the 'main' plots. | |
| :returns: 2-tuple (plotnames,mainnames) | |
| """ | |
| return _plotreg.keys(),_mainfigs | |
| from os.path import abspath | |
| papersdir = abspath('papers') | |
| def make_plots(plots,argd=None,figs=None,save=papersdir,overwrite=True, | |
| showfigs=True,savetypes=['eps','pdf']): | |
| """ | |
| Generate plots and maybe save them. | |
| :param plots: A comma-seperated string of plots to show | |
| :param argd: | |
| The dictionary to use to find the arguments of the plot function. If | |
| None, global variables will be used. | |
| :param figs: | |
| A list of figure numbers to use for each of the figures in `plots`. If | |
| None, the current figure will be used as the starting point. | |
| :param str save: | |
| A directory in which to save the plots. If False, they will not be | |
| saved, if True, the current directory will be used. | |
| :param bool overwrite: | |
| If True and plots are being saved, any existing plot will be | |
| overwritten. Otherwise, the figure will not be saved. | |
| :parm bool showfigs: | |
| If True, the figures will be shown (in whatever GUI is the correct one | |
| for the type). | |
| """ | |
| import os,shutil,inspect,subprocess | |
| from os.path import exists | |
| if save is True: | |
| save='.' | |
| if type(save) is str and not save.endswith(os.sep): | |
| save+=os.sep | |
| if plots is None: | |
| plots = _plotreg.keys() | |
| elif isinstance(plots,basestring): | |
| if '*' in plots: | |
| import re | |
| plots = plots.replace('*','.*') | |
| regex = re.compile(plots) | |
| plots = [p for p in plot_names() if regex.match(p)] | |
| else: | |
| plots = plots.split(',') | |
| if figs is None: | |
| figs = arange(len(plots))+1 | |
| if len(figs) != len(plots): | |
| raise ValueError("plots don't match figs") | |
| if argd is None: | |
| argd = globals() | |
| plotd={} | |
| for p,fnum in zip(plots,figs): | |
| if p not in _plotreg: | |
| raise KeyError('plotname %s not found'%p) | |
| plotd[p] = (_plotreg[p][0],fnum,_plotreg[p][1],_plotreg[p][2],_plotreg[p][3]) | |
| isinter = plt.isinteractive() | |
| plt.ioff() | |
| anympl = False | |
| try: | |
| for fname,(figf,fignum,figsize,ftype,tweakbbox) in plotd.iteritems(): | |
| print 'figure',fname | |
| if 'mpl' in ftype: | |
| anympl = True | |
| if isinter: | |
| plt.ion() | |
| plt.close(fignum) | |
| plt.figure(fignum,figsize=figsize) | |
| plt.ioff() | |
| elif 'mayavi' in ftype: | |
| from enthought.mayavi import mlab as M | |
| if len(M.get_engine().scenes) > 0: | |
| M.gcf().parent.close_scene(M.gcf()) | |
| f=M.figure(fignum,size=figsize) | |
| #this is buggy | |
| #f.scene.set_size(figsize) | |
| #M.clf() | |
| else: | |
| raise ValueError('unrecognized figtype') | |
| argsp = inspect.getargspec(figf) | |
| if argsp[2] is not None: | |
| figf(**argd) | |
| else: | |
| args=[] | |
| for arg in argsp[0]: | |
| args.append(argd[arg]) | |
| figf(*args) | |
| if save: | |
| epsconv = False | |
| if ftype=='mpltoeps': | |
| epsconv = 'png' | |
| savefunc = plt.savefig | |
| saveexts = ['png'] | |
| elif ftype.startswith('mpl'): | |
| savefunc = plt.savefig | |
| saveexts = ftype[3:].strip() | |
| if saveexts=='': | |
| if isinstance(savetypes,basestring): | |
| saveexts = [savetypes] | |
| else: | |
| saveexts = savetypes | |
| else: | |
| saveexts = saveexts.split(',') | |
| elif ftype.startswith('mayavi'): | |
| savefunc = M.savefig | |
| saveexts = ftype[6:].strip() | |
| if saveexts=='': | |
| saveexts = ['png'] | |
| else: | |
| saveexts = saveexts.split(',') | |
| if 'eps' in saveexts: | |
| saveexts.remove('eps') | |
| if 'png' not in saveexts: | |
| saveexts.append('png') | |
| epsconv = 'png' | |
| else: | |
| raise ValueError('unrecognized figtype') | |
| for saveext in saveexts: | |
| fn = '%s%s.%s'%(save,fname,saveext) | |
| epsfns = None | |
| if not overwrite and exists(fn): | |
| print fn,'exists, skipping (use -o at command line to overwrite)' | |
| else: | |
| print 'saving',fn | |
| savefunc(fn) | |
| if saveext == epsconv: | |
| epsfns = (fn,'%s%s.%s'%(save,fname,'eps')) | |
| if epsfns is not None: | |
| from os import system | |
| fn,epsfn = epsfns | |
| print 'converting',fn,'to',epsfn | |
| retcode = subprocess.call('convert %s %s'%(fn,epsfn),shell=True) | |
| if retcode is not 0: | |
| raise ValueError('convert failed to convert %s to %s'%(fn,epsfn)) | |
| if tweakbbox is not None: | |
| print 'updating eps bounding box on',epsfn,'to',tweakbbox | |
| with open(epsfn+'.tmpwrite','w') as fw: | |
| with open(epsfn,'r') as fr: | |
| for l in fr: | |
| if 'BoundingBox' in l and 'Page' not in l: | |
| newline = l.split(':')[0]+(': %s %s %s %s\n'%tweakbbox) | |
| fw.write(newline) | |
| else: | |
| fw.write(l) | |
| shutil.move(epsfn+'.tmpwrite',epsfn) | |
| finally: | |
| if isinter: | |
| plt.ion() | |
| if showfigs and anympl: | |
| plt.draw() | |
| #<------------------functions here - decorating plots with @plotfunc-----------> | |
| def get_phot_mask(vagcd,rlim=(-12,None),gmrlim=(-2,5,10),finitecut=True): | |
| g,r = vagcd['kcorr'].ABSMAG[:,(1,2)].T | |
| gmr = g-r | |
| upperr = max(*rlim) | |
| lowerr = min(*rlim) | |
| if upperr == lowerr: | |
| if upperr<-20: | |
| rcut = lowerr<r | |
| else: | |
| rcut = r<upperr | |
| else: | |
| rcut = (lowerr<r)&(r<upperr) | |
| uppergmr = max(*gmrlim) | |
| lowergmr = min(*gmrlim) | |
| if uppergmr==lowergmr: | |
| if uppergmr < 1: | |
| gmrcut = lowergmr<r | |
| else: | |
| gmrcut = r < uppergmr | |
| else: | |
| gmrcut = (lowergmr<gmr)&(gmr<uppergmr) | |
| if finitecut: | |
| return rcut&gmrcut&isfinite(gmr) | |
| else: | |
| return rcut&gmrcut | |
| def get_skyserver_url(ind,vagc,cuts=None,mode='navigate',dr=7,openpage='tab'): | |
| """ | |
| `ind` can is an index into `cuts` or `vagc` if `cuts` is None | |
| `mode` can be 'navigate','explore', 'chart' or 'cutout'. | |
| `openpage` can be False, True, 'new',or 'tab' | |
| """ | |
| from operator import isMappingType | |
| if not (1<=dr<=7): | |
| raise ValueError('invalid data release') | |
| if isMappingType(vagc): | |
| vagc = vagc['kcorr'] | |
| if isMappingType(cuts): | |
| cuts = cuts['z500'] | |
| if cuts is None: | |
| ra = vagc.RA[ind] | |
| dec = vagc.DEC[ind] | |
| ra2 = dec2 = None | |
| z = vagc.Z[ind] | |
| else: | |
| ra = vagc.RA[cuts.i[ind]] | |
| dec = vagc.DEC[cuts.i[ind]] | |
| ra2 = vagc.RA[cuts.ni[ind]] | |
| dec2 = vagc.DEC[cuts.ni[ind]] | |
| z = vagc.Z[cuts.i[ind]] | |
| if 'navigate' in mode: | |
| opts = mode.replace('navigate','').strip() | |
| mode = 'navigate' | |
| else: | |
| opts = None | |
| url = _construct_url(ra,dec,ra2,dec2,z,dr,mode) | |
| if opts: | |
| url+='&opt='+opts | |
| if openpage: | |
| import webbrowser | |
| if openpage == 'tab': | |
| webbrowser.open_new_tab(url) | |
| if openpage == 'new': | |
| webbrowser.open_new(url) | |
| else: | |
| webbrowser.open(url) | |
| return url | |
| def _construct_url(ra,dec,ra2,dec2,z,dr,mode): | |
| from urllib import urlencode | |
| from astropysics import coords,constants | |
| h = constants.H0/100 | |
| if mode=='navigate': | |
| baseurl = 'http://cas.sdss.org/dr%i/en/tools/chart/navi.asp?'%dr | |
| data = {'ra':ra,'dec':dec} | |
| elif mode=='explore': | |
| baseurl = 'http://cas.sdss.org/dr%i/en/tools/explore/obj.asp?'%dr | |
| data = {'ra':ra,'dec':dec} | |
| elif mode=='cutout' or mode=='chart': | |
| if mode=='cutout': | |
| baseurl = 'http://casjobs.sdss.org/ImgCutoutDR%i/getjpeg.aspx?'%dr | |
| else: | |
| baseurl = 'http://cas.sdss.org/dr%i/en/tools/chart/chart.asp?'%dr | |
| imsz = 512 | |
| #this scale corresponds to dproj of 250 kpc/h | |
| dang = coords.cosmo_z_to_dist(z,disttype='angular')*1000*h | |
| theta = 2*250/dang | |
| scale = 206265*theta/imsz | |
| #this scale is for 1.5*sat-host distance | |
| #asecsep = 3600*((ra-ra2)**2+(dec-dec2)**2)**0.5 | |
| #scale = asecsep*3/imsz | |
| linestr='\r\n' | |
| radecs = ['RA+DEC','{0}+{1}'.format(ra,dec),'{0}+{1}'.format(ra2,dec2)] | |
| data = {'ra':ra,'dec':dec, | |
| 'width':imsz,'height':imsz, | |
| 'scale':scale, | |
| 'opt':'','query':linestr.join(radecs)} | |
| else: | |
| raise ValueError('invalid mode') | |
| url = baseurl+urlencode(data).replace('%2B','+') | |
| return url | |
| def generate_line_queries(vagcd,cutsrecarr,fnbase): | |
| """ | |
| Generates a pair of comma-seperated strings suitable for import into CasJobs, | |
| and write two files beginning with `fnbase` - the .radec file looks like: | |
| #VAGCidhost,VAGCidsat,RAh,Dech,RAs,Decs | |
| 1345656,213123,123.464814864,45.2456,123.462814864,45.1456 | |
| ... | |
| This file should by uploaded to MyDB on SDSS CasJobs (with the same name as | |
| the file's base), and then the returned query should be run to populate the | |
| objids for these objects. | |
| The second file - .query - has the sql queries to apply after the upload | |
| call as e.g. generate_radec_string(vagcd,cutsd['z500'],'lmcana_dr7z500') | |
| """ | |
| ras = vagcd['catalog'].RA[cutsrecarr.ni] | |
| decs = vagcd['catalog'].DEC[cutsrecarr.ni] | |
| rah = vagcd['catalog'].RA[cutsrecarr.i] | |
| dech = vagcd['catalog'].DEC[cutsrecarr.i] | |
| strlines = [','.join([str(ti) for ti in t]) for t in zip(cutsrecarr.i,cutsrecarr.ni,rah,dech,ras,decs)] | |
| strlines.insert(0,'#VAGCidhost,VAGCidsat,RAh,Dech,RAs,Decs') | |
| sradec = '\n'.join(strlines) | |
| with open(fnbase+'.radec','w') as f: | |
| f.write(sradec) | |
| queryids = """ | |
| SELECT VAGCidhost, RAh, Dech, dbo.fGetNearestSpecObjIdEq(RAh,Dech,3.0) as SpecObjIDhost, dbo.fGetNearestObjIdEq(RAh,Dech,3.0) as PriPhotObjIDhost, | |
| VAGCidsat, RAs, Decs, dbo.fGetNearestSpecObjIdEq(RAs,Decs,3.0) as SpecObjIDsat, dbo.fGetNearestObjIdEq(RAs,Decs,3.0) as PriPhotObjIDsat | |
| INTO mydb.{0}id | |
| FROM mydb.{0} | |
| """.format(fnbase) | |
| queryliness = """ | |
| SELECT | |
| o.VAGCidsat, o.RAh as RA, o.Dech as Dec, o.SpecObjIDsat as SpecObjID, dbo.fGetUrlSpecImg(o.SpecObjIDsat) as specimgurl,s.bestObjID as PhotObjID, | |
| L.ew as Ha_ew, L.ewErr as Ha_ewErr, L.continuum as Ha_cont, L.wave as Ha_wave, | |
| L2.ew as Hb_ew, L2.ewErr as Hb_ewErr, L2.continuum as Hb_cont, L2.wave as Hb_wave, | |
| L_Hd.ew as Hd_ew, L_Hd.ewErr as Hd_ewErr,L_Hd.continuum as Hd_cont, L_Hd.wave as Hd_wave, | |
| L_OII.ew as OII_ew,L_OII.ewErr as OII_ewErr,L_OII.continuum as OII_cont, L_OII.wave as OII_wave, | |
| L_MgIb.ew as MgIb_ew,L_MgIb.ewErr as MgIb_ewErr,L_MgIb.continuum as MgIb_cont,L_MgIb.wave as MgIb_wave, | |
| L_NaD.ew as NaD_ew,L_NaD.ewErr as NaD_ewErr,L_NaD.continuum as NaD_cont,L_NaD.wave as NaD_wave, | |
| L_N2.ew as N2_ew,L_N2.ewErr as N2_ewErr,L_N2.continuum as N2_cont,L_N2.wave as N2_wave, | |
| L_OIII.ew as OIII_ew,L_OIII.ewErr as OIII_ewErr,L_OIII.continuum as OIII_cont,L_OIII.wave as OIII_wave | |
| FROM mydb.{0}id o | |
| INTO mydb.{0}satlines | |
| JOIN SpecObj AS S ON o.SpecObjIDsat = S.SpecObjID | |
| JOIN SpecLine AS L ON o.SpecObjIDsat = L.SpecObjID | |
| JOIN SpecLine AS L2 ON o.SpecObjIDsat = L2.SpecObjID | |
| JOIN SpecLine AS L_OII ON o.SpecObjIDsat = L_OII.SpecObjID | |
| JOIN SpecLine AS L_Hd ON o.SpecObjIDsat = L_Hd.SpecObjID | |
| JOIN SpecLine AS L_MgIb ON o.SpecObjIDsat = L_MgIb.SpecObjID | |
| JOIN SpecLine AS L_NaD ON o.SpecObjIDsat = L_NaD.SpecObjID | |
| JOIN SpecLine AS L_N2 ON o.SpecObjIDsat = L_N2.SpecObjID | |
| JOIN SpecLine AS L_OIII ON o.SpecObjIDsat = L_OIII.SpecObjID | |
| WHERE | |
| L.LineId = 6565 | |
| and L2.LineId = 4863 | |
| and L_OII.LineId = 3727 | |
| and L_Hd.LineId = 4103 | |
| and L_MgIb.LineID = 5177 | |
| and L_NaD.LineID = 5896 | |
| and L_N2.LineID = 6585 | |
| and L_OIII.LineID = 5008 | |
| ORDER BY o.VAGCidhost | |
| """.format(fnbase) | |
| querylinesh = """ | |
| SELECT | |
| o.VAGCidhost, o.RAh as RA, o.Dech as Dec, o.SpecObjIDhost as SpecObjID, dbo.fGetUrlSpecImg(o.SpecObjIDhost) as specimgurl,s.bestObjID as PhotObjID, | |
| L.ew as Ha_ew, L.ewErr as Ha_ewErr, L.continuum as Ha_cont, L.wave as Ha_wave, | |
| L2.ew as Hb_ew, L2.ewErr as Hb_ewErr, L2.continuum as Hb_cont, L2.wave as Hb_wave, | |
| L_Hd.ew as Hd_ew, L_Hd.ewErr as Hd_ewErr,L_Hd.continuum as Hd_cont, L_Hd.wave as Hd_wave, | |
| L_OII.ew as OII_ew,L_OII.ewErr as OII_ewErr,L_OII.continuum as OII_cont, L_OII.wave as OII_wave, | |
| L_MgIb.ew as MgIb_ew,L_MgIb.ewErr as MgIb_ewErr,L_MgIb.continuum as MgIb_cont,L_MgIb.wave as MgIb_wave, | |
| L_NaD.ew as NaD_ew,L_NaD.ewErr as NaD_ewErr,L_NaD.continuum as NaD_cont,L_NaD.wave as NaD_wave, | |
| L_N2.ew as N2_ew,L_N2.ewErr as N2_ewErr,L_N2.continuum as N2_cont,L_N2.wave as N2_wave, | |
| L_OIII.ew as OIII_ew,L_OIII.ewErr as OIII_ewErr,L_OIII.continuum as OIII_cont,L_OIII.wave as OIII_wave | |
| FROM mydb.{0}id o | |
| INTO mydb.{0}hostlines | |
| JOIN SpecObj AS S ON o.SpecObjIDhost = S.SpecObjID | |
| JOIN SpecLine AS L ON o.SpecObjIDhost = L.SpecObjID | |
| JOIN SpecLine AS L2 ON o.SpecObjIDhost = L2.SpecObjID | |
| JOIN SpecLine AS L_OII ON o.SpecObjIDhost = L_OII.SpecObjID | |
| JOIN SpecLine AS L_Hd ON o.SpecObjIDhost = L_Hd.SpecObjID | |
| JOIN SpecLine AS L_MgIb ON o.SpecObjIDhost = L_MgIb.SpecObjID | |
| JOIN SpecLine AS L_NaD ON o.SpecObjIDhost = L_NaD.SpecObjID | |
| JOIN SpecLine AS L_N2 ON o.SpecObjIDhost = L_N2.SpecObjID | |
| JOIN SpecLine AS L_OIII ON o.SpecObjIDhost = L_OIII.SpecObjID | |
| WHERE | |
| L.LineId = 6565 | |
| and L2.LineId = 4863 | |
| and L_OII.LineId = 3727 | |
| and L_Hd.LineId = 4103 | |
| and L_MgIb.LineID = 5177 | |
| and L_NaD.LineID = 5896 | |
| and L_N2.LineID = 6585 | |
| and L_OIII.LineID = 5008 | |
| ORDER BY o.VAGCidhost | |
| """.format(fnbase) | |
| with open(fnbase+'.query','w') as f: | |
| f.write('#this query is to be applied to the original uploaded table and determines the SpecObjId of the targets\n') | |
| f.write(queryids) | |
| f.write('\n\n#this query finds the lines for the hosts\n') | |
| f.write(querylinesh) | |
| f.write('\n\n#this query finds the lines for the sats\n') | |
| f.write(queryliness) | |
| def match_line_query_to_cuts(hqueryfn,squeryfn,cutsd,cuttype='z500'): | |
| """ | |
| matches SDSS query results for spectral lines to the cuts by index (the SDSS | |
| query misses a few for unclear reasons). | |
| returns liness,scutsi,linesh,hcutsi | |
| (e.g. do cutsd[cuttype].i[hcutsi] to get the host indecies for hlines) | |
| """ | |
| cuts = cutsd[cuttype] | |
| linesh = pyfits.open(hqueryfn)[1].data | |
| liness = pyfits.open(squeryfn)[1].data | |
| hid = linesh.VAGCidhost | |
| sid = liness.VAGCidsat | |
| sorti = argsort(cuts.i) | |
| sortni = argsort(cuts.ni) | |
| sorti = sorti[np.searchsorted(cuts.i[sorti],hid)] | |
| sortni = sortni[np.searchsorted(cuts.ni[sortni],sid)] | |
| return liness,sortni,linesh,sorti | |
| @plotfunc | |
| def fullcmd(g,r): | |
| xl = (0.2,1.0) | |
| yl = (-17.75,-24.25) | |
| hexbin(g-r,r,bins='log',extent=(xl[0],xl[1],yl[0],yl[1]),marginals=False,reduce_C_function=sum) #marginals not compatible w/cb | |
| xlim(*xl) | |
| ylim(*yl) | |
| cb = colorbar() | |
| cb.set_label(r'$\log(1+N_{gal})$') | |
| xlabel('$g-r$') | |
| ylabel('$r$') | |
| @plotfunc | |
| def fullcmdBV(B,V): | |
| hexbin(B-V,B,bins='log') | |
| xlim(0.4,1.2) | |
| ylim(-17,-23) | |
| cb = colorbar() | |
| cb.set_label(r'$\log(1+N_{gal})$') | |
| xlabel('$B-V$') | |
| ylabel('$B$') | |
| @plotfunc | |
| def potscmd(g,r,potsd,galmask,h): | |
| subplots_adjust(top=.98,bottom=.08,left=.12,right=.97) | |
| pothmask = np.zeros_like(galmask) | |
| pothmask[potsd['pothosts']] = True | |
| pothmask = pothmask[galmask] | |
| potsmask = np.zeros_like(galmask) | |
| potsmask[potsd['potsats']] = True | |
| potsmask = potsmask[galmask] | |
| c = g - r | |
| xl = (0.2,1.0) | |
| yl = (-17.5 + 5*log10(h),-22.5 + 5*log10(h)) | |
| ycut = -20 + 5*log10(h) | |
| xsz = 50 | |
| ysz1 = int(10*(ycut - yl[1])) | |
| ysz2 = int(10*(yl[0] - ycut)) | |
| #marginals not compatible w/cb | |
| sarr = hexbin(c[potsmask],r[potsmask],gridsize=(xsz,ysz2),bins='log', | |
| extent=(xl[0],xl[1],yl[0],ycut)).get_array() | |
| vmin,vmax = sarr.min(),sarr.max() | |
| harr = hexbin(c[pothmask],r[pothmask],gridsize=(xsz,ysz1),bins='log', | |
| extent=(xl[0],xl[1],yl[1],ycut),vmin=vmin,vmax=vmax).get_array() | |
| hlines(ycut,xl[0],xl[1],colors='k') | |
| xlim(*xl) | |
| ylim(*yl) | |
| cb = colorbar() | |
| #cb.set_label(r'$\log(1+N_{\rm gal})$') | |
| cb.set_label(r'$N_{\rm gal}$') | |
| tcks = array([0,1,2,3,5,10,20,30,50,100,200]) | |
| cb.set_ticks(log10(tcks+1)) | |
| cb.set_ticklabels(['%i'%tck for tck in tcks]) | |
| xlabel('$g-r$') | |
| ylabel('$r$') | |
| @plotfunc | |
| def potscmdsec(g,r,potsd,galmask,h): | |
| subplots_adjust(top=.98,bottom=.08,left=.15,right=.97) | |
| pothmask = np.zeros_like(galmask) | |
| pothmask[potsd['pothosts']] = True | |
| pothmask = pothmask[galmask] | |
| potsmask = np.zeros_like(galmask) | |
| potsmask[potsd['potsats']] = True | |
| potsmask = potsmask[galmask] | |
| c = g - r | |
| xl = (0.2,1.0) | |
| yl = (-17.5 + 5*log10(h),-20 + 5*log10(h)) | |
| xsz = 50 | |
| ysz = int(10*(yl[0] - yl[1])) | |
| hexbin(c[potsmask],r[potsmask],gridsize=(xsz,ysz),bins='log',extent=(xl[0],xl[1],yl[0],yl[1])) #marginals not compatible w/cb | |
| xlim(*xl) | |
| ylim(*yl) | |
| cb = colorbar() | |
| #cb.set_label(r'$\log(1+N_{\rm gal})$') | |
| cb.set_label(r'$N_{\rm gal}$') | |
| tcks = array([0,1,2,3,5,10,20,30,50,100,200]) | |
| cb.set_ticks(log10(tcks+1)) | |
| cb.set_ticklabels(['%i'%tck for tck in tcks]) | |
| xlabel('$g-r$') | |
| ylabel('$r$') | |
| @plotfunc | |
| def potscmdBV(B,V,potsd,mask): | |
| pothmask = np.zeros_like(mask) | |
| pothmask[potsd['pothosts']] = True | |
| pothmask = pothmask[mask] | |
| potsmask = np.zeros_like(mask) | |
| potsmask[potsd['potsats']] = True | |
| potsmask = potsmask[mask] | |
| c = B-V | |
| xl = (0.4,1.2) | |
| yl = (-17.2,-22) | |
| ycut = -19 | |
| xsz = 50 | |
| ysz1 = int(10*(ycut - yl[1])) | |
| ysz2 = int(10*(yl[0] - ycut)) | |
| hexbin(c[pothmask],B[pothmask],gridsize=(xsz,ysz1),bins='log',extent=(xl[0],xl[1],yl[1],ycut)) #marginals not compatible w/cb | |
| hexbin(c[potsmask],B[potsmask],gridsize=(xsz,ysz2),bins='log',extent=(xl[0],xl[1],yl[0],ycut)) #marginals not compatible w/cb | |
| hlines(ycut,xl[0],xl[1],colors='k') | |
| xlim(*xl) | |
| ylim(*yl) | |
| cb = colorbar() | |
| cb.set_label(r'$\log(1+N_{gal})$') | |
| xlabel('$B-V$') | |
| ylabel('$B$') | |
| @plotfunc(mainfig=True) | |
| def matchcmd(g,r,potsd,cutsd,hiicutsd,galmask,LMCappmags,LMCbothun,h): | |
| cutsd = hiicutsd #comment to get non HII-removed | |
| potscmd(g,r,potsd,galmask,h) | |
| dcut = 250 | |
| xl = xlim() | |
| yl = ylim() | |
| rcut = -20 + 5*log10(h) | |
| dcutmask = cutsd['z500'].nd < dcut | |
| z500hmask = np.zeros_like(galmask) | |
| z500hmask[cutsd['z500'].i[dcutmask]] = True | |
| z500hmask = z500hmask[galmask] | |
| #z500hmask[r>rcut] = False | |
| z500smask = np.zeros_like(galmask) | |
| z500smask[cutsd['z500'].ni[dcutmask]] = True | |
| z500smask = z500smask[galmask] | |
| #z500smask[r<rcut] = False | |
| scatter((g-r)[z500hmask],r[z500hmask],c='r',edgecolor='k',lw=1,label='Primaries') | |
| scatter((g-r)[z500smask],r[z500smask],c='w',edgecolor='k',lw=1,label='Secondaries') | |
| print 'Median pri M_r =',median(r[z500hmask]) | |
| print 'Median sec M_r =',median(r[z500smask]) | |
| print 'Median blue M_r =',median(r[z500smask][(g-r)[z500smask] <.6]) | |
| # lmcr = LMCappmags['r']-LMCdm | |
| # scatter([LMCappmags['g']-LMCappmags['r']],[lmcr],c='r',marker=(7,1,0),s=250) | |
| # scatter([LMCbothun['g-r']],[lmcr],c='b',marker=(7,1,0),s=250) | |
| # #now get LMC track from Eskew 10 | |
| # lmchist = mcpsscripts.load_peg_data('MCPS/lmc_colors_spectra1011.dat')[:-5] #stoping @ 14.5 Gyr | |
| # lmchr = lmchist["V"] - lmchist["V-r'"] +2.5*log10(lmchist['M*'][-1]) - 2.5*log10(2.2e9) #H&Z09 SFH | |
| # print 'Eskew g-r=',lmchist["g'-r'"][-1],'r=',lmchr[-1] | |
| # lmchr -= lmchr[-1] - lmcr #forces it to match the final magnitude of the others | |
| # lmchgmr = lmchist["g'-r'"] | |
| # plot(lmchgmr,lmchr,'m') | |
| # scatter([lmchgmr[-1]],[lmchr[-1]],c='m',marker=(7,1,0),s=250) | |
| ylabel('$M_r$') | |
| legend(loc=2,prop={'size':18}) | |
| xlim(*xl) | |
| ylim(*yl) | |
| @plotfunc(mainfig=True) | |
| def matchcmdsec(g,r,potsd,cutsd,hiicutsd,galmask,LMCappmags,LMCbothun,h): | |
| cutsd = hiicutsd #comment to get non HII-removed | |
| potscmdsec(g,r,potsd,galmask,h) | |
| dcut = 250*h | |
| xl = xlim() | |
| yl = ylim() | |
| rcut = -20 + 5*log10(h) | |
| dcutmask = cutsd['z500'].nd < dcut | |
| z500smask = np.zeros_like(galmask) | |
| z500smask[cutsd['z500'].ni[dcutmask]] = True | |
| z500smask = z500smask[galmask] | |
| scatter((g-r)[z500smask],r[z500smask],c='w') | |
| lmcr = LMCappmags['r']-LMCdm | |
| lmcrbothun = -18.93 | |
| scatter([LMCappmags['g']-LMCappmags['r']],[lmcr],c='r',marker=(7,1,0),s=250,edgecolor='w') | |
| scatter([LMCbothun['g-r']],[lmcrbothun],c='b',marker=(7,1,0),s=250,edgecolor='w') | |
| #now get LMC track from Eskew 10 | |
| lmchist = mcpsscripts.load_peg_data('MCPS/lmc_colors_spectra1011.dat')[:-5] #stoping @ 14.5 Gyr | |
| lmchr = lmchist["V"] - lmchist["V-r'"] +2.5*log10(lmchist['M*'][-1]) - 2.5*log10(2.2e9) #H&Z09 SFH | |
| print 'Eskew g-r=',lmchist["g'-r'"][-1],'r=',lmchr[-1] | |
| lmchr -= lmchr[-1] - lmcr #forces it to match the final magnitude of the others | |
| lmchgmr = lmchist["g'-r'"] | |
| #plot(lmchgmr,lmchr,'m') | |
| scatter([lmchgmr[-1]],[lmchr[-1]],c='m',marker=(7,1,0),s=250,edgecolor='w') | |
| m33mag = -19.06 | |
| m33color = .44 | |
| scatter([m33color],[m33mag],c='k',marker=(7,1,0),s=250,edgecolor='w') | |
| ylabel('$M_r$') | |
| xlim(*xl) | |
| ylim(*yl) | |
| @plotfunc(mainfig=True) | |
| def colorhist(g,r,potsd,cutsd,hiicutsd,galmask,LMCappmags,LMCbothun,h): | |
| from scipy.stats import percentileofscore | |
| from matplotlib.transforms import blended_transform_factory | |
| cutsd = hiicutsd #comment to get non HII-removed | |
| crange = [.2,1] | |
| #SDSS Secondaries | |
| dcut = 250*h | |
| #dcut = 100*h | |
| dcutmask = cutsd['z500'].nd < dcut | |
| z500hmask = np.zeros_like(galmask) | |
| z500hmask[cutsd['z500'].i[dcutmask]] = True | |
| z500hmask = z500hmask[galmask] | |
| z500smask = np.zeros_like(galmask) | |
| z500smask[cutsd['z500'].ni[dcutmask]] = True | |
| z500smask = z500smask[galmask] | |
| pricolor = (g-r)[z500hmask] | |
| primag = r[z500hmask] | |
| seccolor = (g-r)[z500smask] | |
| secmag = r[z500smask] | |
| #LMCs | |
| lmcmag = LMCappmags['r']-LMCdm | |
| lmccRC3 = LMCappmags['g']-LMCappmags['r'] | |
| lmccbothun = LMCbothun['g-r'] | |
| #now get LMC track from Eskew 10 | |
| lmchist = mcpsscripts.load_peg_data('MCPS/lmc_colors_spectra1011.dat')[:-5] #stoping @ 14.5 Gyr | |
| lmchr = lmchist["V"] - lmchist["V-r'"] +2.5*log10(lmchist['M*'][-1]) - 2.5*log10(2.2e9) #H&Z09 SFH | |
| lmchr -= lmchr[-1] - lmcmag #forces it to match the final magnitude of the others | |
| lmcceskew = lmchist["g'-r'"][-1] | |
| #Potentials | |
| pothmask = np.zeros_like(galmask) | |
| pothmask[potsd['pothosts']] = True | |
| pothmask = pothmask[galmask] | |
| potsmask = np.zeros_like(galmask) | |
| potsmask[potsd['potsats']] = True | |
| potsmask = potsmask[galmask] | |
| c = g - r | |
| potpricolor = c[pothmask] | |
| potseccolor = c[potsmask] | |
| n,b,pa = hist(seccolor,fc=(.7,.7,1),ec='k',lw=1,bins=25,histtype='stepfilled',normed=True,range=crange,label='Clean Sample Secondaries',zorder=1) | |
| pn,pb,ppa = hist(potseccolor,color='g',lw=3,bins=100,histtype='step',normed=True,range=crange,label='Control Sample',zorder=2) | |
| # n,b,pa = hist(seccolor,color='g',lw=3,bins=25,histtype='step',normed=True,range=crange,label='Clean Sample Secondaries',zorder=2) | |
| # pn,pb,ppa = hist(potseccolor,fc=(.7,.7,1),ec='k',lw=1,bins=100,histtype='stepfilled',normed=True,range=crange,label='Potential Secondaries',zorder=1) | |
| print 'Clean f(red)=',sum(seccolor>.6)/seccolor.size | |
| print 'Control f(red)=',sum(potseccolor>.6)/potseccolor.size | |
| xdyaxtrans = blended_transform_factory(gca().transData,gca().transAxes) | |
| showp = True | |
| prc3pot = percentileofscore(potseccolor,lmccRC3) | |
| prc3sel = percentileofscore(seccolor,lmccRC3) | |
| rc3l = axvline(lmccRC3,0.7,0.77,linestyle='dashed',color='r',lw=2,label=r'LMC (RC3) $p=%.1f\%%$'%prc3sel,zorder=3) | |
| rc3l.set_label('_'+rc3l.get_label()) | |
| if showp: | |
| text(lmccRC3,.77,r' $%.1f\%%$'%prc3sel,ha='center',va='bottom',transform=xdyaxtrans) | |
| print 'RC3 percentile (pot,sel)',prc3pot,prc3sel | |
| peskewpot = percentileofscore(potseccolor,lmcceskew) | |
| peskewsel = percentileofscore(seccolor,lmcceskew) | |
| ekl = axvline(lmcceskew,0.7,0.77,linestyle='dashed',color='m',lw=2,label=r'LMC (Eskew 2010) $p=%.1f\%%$'%peskewsel,zorder=3) | |
| ekl.set_label('_'+ekl.get_label()) | |
| if showp: | |
| text(lmcceskew,.77,r' $%.1f\%%$'%peskewsel,ha='center',va='bottom',transform=xdyaxtrans) | |
| print 'Eskew percentile (pot,sel)',peskewpot,peskewsel | |
| pbothunpot = percentileofscore(potseccolor,lmccbothun) | |
| pbothunsel = percentileofscore(seccolor,lmccbothun) | |
| bol = axvline(lmccbothun,0.7,0.77,linestyle='dashed',color='b',lw=2,label=r'LMC (Bothun 1988) $p=%.1f\%%$'%pbothunsel,zorder=3) | |
| bol.set_label('_'+bol.get_label()) | |
| if showp: | |
| text(lmccbothun,.77,r' $%.1f\%%$'%pbothunsel,ha='center',va='bottom',transform=xdyaxtrans) | |
| print 'Bothun percentile (pot,sel)',pbothunpot,pbothunsel | |
| lmcmeancolor = average([lmccRC3,lmcceskew,lmccbothun]) | |
| meanp = percentileofscore(potseccolor,lmcmeancolor) | |
| means = percentileofscore(seccolor,lmcmeancolor) | |
| print 'LMC Mean,ppot,psel',lmcmeancolor,meanp,means | |
| arrow(lmcmeancolor,2.5,0,-.9,head_length=.5,head_width=.05,overhang=1, | |
| length_includes_head=True,lw=3,fc='k') | |
| text(lmcmeancolor,2.6,'LMC',ha='center',va='bottom',fontsize=24) | |
| m33color = .44 | |
| pm33sel = percentileofscore(seccolor,m33color) | |
| pm33pot = percentileofscore(potseccolor,m33color) | |
| m33l = axvline(m33color,0.88,0.95,linestyle='dashed',color='k',lw=2,label=r'M33 $p=%.1f\%%$'%pm33sel,zorder=3) | |
| m33l.set_label('_'+m33l.get_label()) | |
| if showp: | |
| text(m33color,.95,r' $%.1f\%%$'%pm33sel,ha='center',va='bottom',transform=xdyaxtrans) | |
| arrow(m33color,3.3,0,-.9,head_length=.5,head_width=.05,overhang=1, | |
| length_includes_head=True,lw=3,fc='k') | |
| text(m33color,3.4,'M33',ha='center',va='bottom',fontsize=24) | |
| xlim(*crange) | |
| ylim(0,4.2) | |
| yticks([]) | |
| xlabel('g-r') | |
| ylabel(r'$N_{\rm gal}$ (normalized)') | |
| legend(loc=0) | |
| def colorhist_table(g,r,potsd,cutsd,hiicutsd,galmask,LMCappmags,LMCbothun,h,fnclean=None,fncontrol=None): | |
| from scipy.stats import percentileofscore | |
| from matplotlib.transforms import blended_transform_factory | |
| cutsd = hiicutsd #comment to get non HII-removed | |
| crange = [.2,1] | |
| #SDSS Secondaries | |
| dcut = 250*h | |
| #dcut = 100*h | |
| dcutmask = cutsd['z500'].nd < dcut | |
| z500hmask = np.zeros_like(galmask) | |
| z500hmask[cutsd['z500'].i[dcutmask]] = True | |
| z500hmask = z500hmask[galmask] | |
| z500smask = np.zeros_like(galmask) | |
| z500smask[cutsd['z500'].ni[dcutmask]] = True | |
| z500smask = z500smask[galmask] | |
| pricolor = (g-r)[z500hmask] | |
| primag = r[z500hmask] | |
| seccolor = (g-r)[z500smask] | |
| secmag = r[z500smask] | |
| #LMCs | |
| lmcmag = LMCappmags['r']-LMCdm | |
| lmccRC3 = LMCappmags['g']-LMCappmags['r'] | |
| lmccbothun = LMCbothun['g-r'] | |
| #now get LMC track from Eskew 10 | |
| lmchist = mcpsscripts.load_peg_data('MCPS/lmc_colors_spectra1011.dat')[:-5] #stoping @ 14.5 Gyr | |
| lmchr = lmchist["V"] - lmchist["V-r'"] +2.5*log10(lmchist['M*'][-1]) - 2.5*log10(2.2e9) #H&Z09 SFH | |
| lmchr -= lmchr[-1] - lmcmag #forces it to match the final magnitude of the others | |
| lmcceskew = lmchist["g'-r'"][-1] | |
| #Potentials | |
| pothmask = np.zeros_like(galmask) | |
| pothmask[potsd['pothosts']] = True | |
| pothmask = pothmask[galmask] | |
| potsmask = np.zeros_like(galmask) | |
| potsmask[potsd['potsats']] = True | |
| potsmask = potsmask[galmask] | |
| c = g - r | |
| potpricolor = c[pothmask] | |
| potseccolor = c[potsmask] | |
| cleanstr = '\n'.join([str(e) for e in seccolor]) | |
| cleanstr = '#Clean sample secondaries\n#g-r\n' + cleanstr | |
| controlstr = '\n'.join([str(e) for e in potseccolor]) | |
| controlstr = '#Control sample\n#g-r\n' + controlstr | |
| if fnclean: | |
| with open(fnclean, 'w') as f: | |
| f.write(cleanstr) | |
| if fncontrol: | |
| with open(fncontrol, 'w') as f: | |
| f.write(controlstr) | |
| return cleanstr, controlstr | |
| @plotfunc | |
| def matchcmdBV(B,V,potsd,cutsd,mask): | |
| potscmdBV(B,V,potsd,mask) | |
| dcut = 100 | |
| xl = xlim() | |
| yl = ylim() | |
| dcutmask = cutsd['z500'].nd < dcut | |
| z500hmask = np.zeros_like(mask) | |
| z500hmask[cutsd['z500'].i[dcutmask]] = True | |
| z500hmask = z500hmask[mask] | |
| z500smask = np.zeros_like(mask) | |
| z500smask[cutsd['z500'].ni[dcutmask]] = True | |
| z500smask = z500smask[mask] | |
| c = B-V | |
| scatter(c[z500hmask],B[z500hmask],c='r') | |
| scatter(c[z500smask],B[z500smask],c='g') | |
| scatter([LMCappmags['B']-LMCappmags['V']],[LMCappmags['B']-LMCdm],c='r',marker=(7,1,0),s=250) | |
| ams = [M32appmags,NGC205appmags,NGC185appmags,NGC147appmags,M33appmags] | |
| dms = [M32dm,NGC205dm,NGC185dm,NGC147dm,M33dm] | |
| for am,dm in zip(ams,dms): | |
| scatter([am['B']-am['V']],[am['B']-dm],c='r',marker=(7,1,0),s=250) | |
| xlim(*xl) | |
| ylim(*yl) | |
| @plotfunc(mainfig=True) | |
| def fracrad(cutsd,hiicutsd,nhiicutsd,potsd,zsim,bolshoid,msiicone,LMCdm,zvol,h): | |
| from astropysics.plotting import cumulative_plot | |
| from astropysics.phot import distance_from_modulus | |
| from astropysics.coords import cosmo_z_to_dist | |
| from astropysics.constants import asecperrad | |
| from astropysics.models import NFWModel | |
| from scipy.interpolate import interp1d | |
| from matplotlib.transforms import blended_transform_factory | |
| cutsd = hiicutsd #comment this to go back to un-visually filtered for HII regions | |
| subplots_adjust(left=.15,right=.95,top=.95) | |
| LMCd = distance_from_modulus(LMCdm)*h/1000 | |
| zsimana = zsim[zsim.nearest_lmc_analog] #comoving kpc/h | |
| zsimsep = zsimana.minsepkpch | |
| fibcoll = zsimana.fibercoll | |
| redge = 250 #in kpc/h | |
| yl = [0,.41] | |
| nds = cutsd['z500'].nd | |
| ndsspec = cutsd['nozcut'].nd | |
| nds1 = cutsd['z500'].nd[cutsd['z500'].nn == 1] | |
| npots = len(potsd['pothosts']) | |
| ndsuniquenoz = cutsd['uniquenoz'].nd | |
| ndsphot = np.concatenate((ndsuniquenoz,ndsspec)) | |
| selfuncx = sort(ndsphot)#kpc/h | |
| selfuncphotf = (arange(ndsphot.size)+1)/npots | |
| selfuncz500f = interp(selfuncx,sort(nds),(arange(nds.size)+1)/npots) | |
| selfuncspecf = interp(selfuncx,sort(ndsspec),(arange(ndsspec.size)+1)/npots) | |
| selfunccorrf = selfuncphotf*(selfuncz500f/selfuncspecf) | |
| err = selfunccorrf*((selfuncphotf*npots)**-1+(selfuncz500f*npots)**-1+(selfuncspecf*npots)**-1)**0.5 | |
| plot(selfuncx/h,selfunccorrf,'b-',c=(0,0,0.5),lw=4, | |
| label='Clean Sample, Selec. Corr. (Observed) ',zorder=10) | |
| fill_between(selfuncx/h,selfunccorrf-err,selfunccorrf+err,color=(0,0,0.5,0.35),zorder=-1) | |
| print 'selfuncx@250:',interp(250,selfuncx/h,selfunccorrf) | |
| msiiinterps = [] | |
| nmsiilines = 50 | |
| for i in range(msiicone.n): | |
| d = msiicone.get_masked_arr('rproj','isolated,deltav,conez',i) #Mpc/h | |
| iso = msiicone.get_masked_arr('isolated','conez',i) | |
| n = np.sum(iso) | |
| d = d[d<(redge/1000)]*1e3/h #now in kpc, not Mpc/h, cut on 250 kpc/h | |
| res = cumulative_plot(d,edges=(0,redge/h),frac=d.size/n,doplot=i<nmsiilines, | |
| label=r'Simulated Clean Sample ($\Lambda$CDM)' if i==0 else '', | |
| zorder=9,ls='-',c='g',lw=1,alpha=.4) | |
| msiiinterps.append(interp1d(*res)) | |
| msiid = linspace(0,redge/h,50) | |
| msiifs = array([[intr(x) for intr in msiiinterps] for x in msiid]) | |
| msiimeans = mean(msiifs,axis=-1) | |
| msiistds = std(msiifs,axis=-1) | |
| #grey band of 1 sigma | |
| # plot(msiid,msiimeans,'-k',lw=3,zorder=1000) | |
| # fill_between(msiid,msiimeans+msiistds,msiimeans-msiistds, | |
| # facecolor=(0.75,0.75,0.75,1),lw=0) | |
| fiberd = cosmo_z_to_dist(zvol,disttype='angular')*(55/asecperrad)*1e3 #Mpc->kpc | |
| vlines(fiberd,yl[0],yl[1],zorder=1,linestyle='dashed',color='k')#,label=r'Fiber radius at $z_{\rm lim}$') | |
| xdyaxtrans = blended_transform_factory(gca().transData,gca().transAxes) | |
| text(fiberd+1,.95,r'Fiber radius at $z_{\rm lim}$',horizontalalignment='left',fontsize=16,transform=xdyaxtrans) | |
| matchrad = LMCd/h | |
| matchedx = linspace(matchrad*h,redge,100)/h | |
| matchedf = interp(matchedx,selfuncx/h,selfunccorrf) | |
| matchedf += interp(matchrad,msiid,msiimeans)-matchedf[0] | |
| # plot(matchedx,matchedf,'r--',lw=2, | |
| # label=r'SDSS Matched at $d_{\rm LMC}$',zorder=10) | |
| legend(loc=4) | |
| ylabel(r'$f_{\rm sat}(<d_{\rm proj})$') | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| ylim(*yl) | |
| xlim(0,250) | |
| @plotfunc(mainfig=True) | |
| def fracraddn(cutsd,hiicutsd,nhiicutsd,potsd,zsim,bolshoid,msiicone,LMCdm,zvol,h): | |
| from astropysics.plotting import cumulative_plot | |
| from astropysics.phot import distance_from_modulus | |
| from astropysics.coords import cosmo_z_to_dist | |
| from astropysics.constants import asecperrad | |
| from astropysics.models import NFWModel | |
| from scipy.interpolate import interp1d | |
| from matplotlib.transforms import blended_transform_factory | |
| from scipy.stats import chisqprob | |
| jackknifeobs = 2 | |
| jackknifesim = None | |
| cutsd = hiicutsd #comment this to go back to un-visually filtered for HII regions | |
| subplots_adjust(left=.17,right=.965,top=.95) | |
| LMCd = distance_from_modulus(LMCdm)*h/1000 | |
| zsimana = zsim[zsim.nearest_lmc_analog] #comoving kpc/h | |
| zsimsep = zsimana.minsepkpch | |
| fibcoll = zsimana.fibercoll | |
| redge = 250 #in kpc/h | |
| nbins = 40 | |
| binedges = linspace(0,redge,nbins+1) | |
| bincens = correlate(binedges,[.5,.5],mode='valid') | |
| dbin = bincens[1]-bincens[0] | |
| nds = cutsd['z500'].nd | |
| ndsspec = cutsd['nozcut'].nd | |
| nds1 = cutsd['z500'].nd[cutsd['z500'].nn == 1] | |
| ndsuniquenoz = cutsd['uniquenoz'].nd | |
| ndsphot = np.concatenate((ndsuniquenoz,ndsspec)) | |
| if jackknifeobs is not None: | |
| nds = nds[np.random.permutation(nds.size)[:nds.size/jackknifeobs]] | |
| ndsspec = ndsspec[np.random.permutation(ndsspec.size)[:ndsspec.size/jackknifeobs]] | |
| nds1 = nds1[np.random.permutation(nds1.size)[:nds1.size/jackknifeobs]] | |
| ndsuniquenoz = ndsuniquenoz[np.random.permutation(ndsuniquenoz.size)[:ndsuniquenoz.size/jackknifeobs]] | |
| ndsphot = ndsphot[np.random.permutation(ndsphot.size)[:ndsphot.size/jackknifeobs]] | |
| npots = len(potsd['pothosts']) | |
| selfuncphot = histogram(ndsphot/h,binedges)[0] | |
| selfuncz500 = histogram(nds/h,binedges)[0] | |
| selfuncspec = histogram(ndsspec/h,binedges)[0] | |
| selfunccorr = selfuncphot*(selfuncz500/selfuncspec) | |
| #plot(selfuncx/h,selfunccorr,'b-',c=(0,0,0.5),lw=4, | |
| # label='SDSS Projected',zorder=10) | |
| f = selfunccorr/npots | |
| if jackknifeobs is not None: | |
| f*=jackknifeobs | |
| obsx = bincens/h | |
| obsy = f | |
| obsyerr = f*(1/selfunccorr+1/npots)**0.5 | |
| l1 = errorbar(obsx,obsy,obsyerr,fmt='bo',ecolor='k', | |
| label='Clean Sample, Selec. Corr. (Observed)',zorder=9)[0] | |
| msiifs = [] | |
| msiinisos = [] | |
| for i in range(msiicone.n): | |
| d = msiicone.get_masked_arr('rproj','isolated,deltav,conez',i) #Mpc/h | |
| iso = msiicone.get_masked_arr('isolated','conez',i) | |
| niso = np.sum(iso) | |
| d = d[d<(redge/1000)]*1e3/h #now in kpc, cut on 250 kpc/h | |
| if jackknifesim is None: | |
| fs = histogram(d,binedges)[0]/niso | |
| else: | |
| d = d[np.random.permutation(d.size)[:d.size/jackknifesim]] | |
| fs = jackknifesim*histogram(d,binedges)[0]/niso | |
| # y = concatenate(([0],ns)) | |
| # plot(binedges/h,y,'g-',drawstyle='steps-pre',lw=1,alpha=.4, | |
| # label='MS-II Projected' if i==0 else '') | |
| msiifs.append(fs) | |
| msiinisos.append(niso) | |
| msiifs = array(msiifs) | |
| msiimeans = msiifs.mean(axis=0) | |
| msiistd = msiifs.std(axis=0) | |
| chi2mask = obsx<=250 | |
| chi2 = ((obsy - msiimeans)**2/(msiistd**2+obsyerr**2))[chi2mask].sum() | |
| dof = len(msiimeans) #NOT n-1, b/c the model is not the mean | |
| print 'rchi2',chi2/dof,'p',chisqprob(chi2,dof) | |
| text(180,.022,r'$\chi^2_{\rm red}=%.2f$'%(chi2/dof),fontsize=20) | |
| bar(binedges[:-1]/h,msiistd*2,dbin/h,msiimeans-msiistd, | |
| color=(0.25,.75,0.25),lw=0, | |
| label=r'Simulated Clean Sample $\pm 1\sigma$ ($\Lambda$CDM)') | |
| fiberd = cosmo_z_to_dist(zvol,disttype='angular')*(55/asecperrad)*1e3 #Mpc->kpc | |
| axvline(fiberd,zorder=1,linestyle='dashed',color='k')#,label=r'Fiber radius at $z_{\rm lim}$') | |
| xdyaxtrans = blended_transform_factory(gca().transData,gca().transAxes) | |
| text(fiberd+1,.95,r'Fiber radius at $z_{\rm lim}$',horizontalalignment='left',fontsize=16,transform=xdyaxtrans) | |
| l = legend(loc=4) | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| xlim(0,250) | |
| #ylim(0,.031) | |
| ylabel(r'$f_{\rm sat,bin}(d_{\rm proj})$') | |
| @plotfunc(mainfig=True) | |
| def fracraddn2(cutsd,hiicutsd,nhiicutsd,potsd,zsim,bolshoid,msiicone,LMCdm,zvol,h,b03tomsiiv,vagcd): | |
| from astropysics.plotting import cumulative_plot | |
| from astropysics.phot import distance_from_modulus | |
| from astropysics.coords import cosmo_z_to_dist | |
| from astropysics.constants import asecperrad | |
| from astropysics.models import NFWModel | |
| from scipy.interpolate import interp1d | |
| from matplotlib.transforms import blended_transform_factory | |
| from scipy.stats import chisqprob | |
| cutsd = hiicutsd #comment this to go back to un-visually filtered for HII regions | |
| magcut = -22 | |
| magcuth = magcut -5*log10(h) | |
| vmcut = b03tomsiiv(magcuth) | |
| subplots_adjust(left=.17,right=.965,top=.95) | |
| LMCd = distance_from_modulus(LMCdm)*h/1000 | |
| zsimana = zsim[zsim.nearest_lmc_analog] #comoving kpc/h | |
| zsimsep = zsimana.minsepkpch | |
| fibcoll = zsimana.fibercoll | |
| redge = 250 #in kpc/h | |
| nbins = 40 | |
| binedges = linspace(0,redge,nbins+1) | |
| bincens = correlate(binedges,[.5,.5],mode='valid') | |
| dbin = bincens[1]-bincens[0] | |
| z5mask = vagcd['kcorr'].ABSMAG[cutsd['z500'].i,2] > magcuth | |
| nzmask = vagcd['kcorr'].ABSMAG[cutsd['nozcut'].i,2] > magcuth | |
| unzmask = vagcd['kcorr'].ABSMAG[cutsd['uniquenoz'].i,2] > magcuth | |
| phmask = vagcd['kcorr'].ABSMAG[potsd['pothosts'],2] > magcuth | |
| nds = cutsd['z500'][z5mask].nd | |
| ndsspec = cutsd['nozcut'][nzmask].nd | |
| nds1 = cutsd['z500'][z5mask].nd[cutsd['z500'][z5mask].nn == 1] | |
| npots = len(potsd['pothosts'][phmask]) | |
| ndsuniquenoz = cutsd['uniquenoz'].nd[unzmask] | |
| ndsphot = np.concatenate((ndsuniquenoz,ndsspec)) | |
| selfuncphot = histogram(ndsphot/h,binedges)[0] | |
| selfuncz500 = histogram(nds/h,binedges)[0] | |
| selfuncspec = histogram(ndsspec/h,binedges)[0] | |
| selfunccorr = selfuncphot*(selfuncz500/selfuncspec) | |
| #plot(selfuncx/h,selfunccorr,'b-',c=(0,0,0.5),lw=4, | |
| # label='SDSS Projected',zorder=10) | |
| f = selfunccorr/npots | |
| obsx = bincens/h | |
| obsy = f | |
| obsyerr = f*(1/selfunccorr+1/npots)**0.5 | |
| l1 = errorbar(obsx,obsy,obsyerr,fmt='bo',ecolor='k', | |
| label='Clean Sample, $M_r>%i$, Selec. Corr. (Observed)'%magcut,zorder=9)[0] | |
| msiifs = [] | |
| msiinisos = [] | |
| for i in range(msiicone.n): | |
| d = msiicone.get_masked_arr('rproj','isolated,deltav,conez',i) #Mpc/h | |
| iso = msiicone.get_masked_arr('isolated','conez',i) | |
| #now get the ones in the vmaxcut | |
| vmh = msiicone.get_masked_arr('h_vmax','isolated,deltav,conez',i) | |
| vmha = msiicone.get_masked_arr('h_vmax','conez',i) | |
| d = d[vmh<vmcut] | |
| iso = iso[vmha<vmcut] | |
| niso = np.sum(iso) | |
| d = d[d<(redge/1000)]*1e3/h #now in kpc, cut on 250 kpc/h | |
| fs = histogram(d,binedges)[0]/niso | |
| # y = concatenate(([0],ns)) | |
| # plot(binedges/h,y,'g-',drawstyle='steps-pre',lw=1,alpha=.4, | |
| # label='MS-II Projected' if i==0 else '') | |
| msiifs.append(fs) | |
| msiinisos.append(niso) | |
| msiifs = array(msiifs) | |
| msiimeans = msiifs.mean(axis=0) | |
| msiistd = msiifs.std(axis=0) | |
| chi2mask = obsx<=250 | |
| chi2 = ((obsy - msiimeans)**2/(msiistd**2+obsyerr**2))[chi2mask].sum() | |
| dof = len(msiimeans) #NOT n-1, b/c the model is not the mean | |
| print 'rchi2',chi2/dof,'p',chisqprob(chi2,dof) | |
| text(180,.022,r'$\chi^2_{\rm red}=%.2f$'%(chi2/dof),fontsize=20) | |
| bar(binedges[:-1]/h,msiistd*2,dbin/h,msiimeans-msiistd, | |
| color=(0.25,.75,0.25),lw=0, | |
| label=r'Simulated Clean Sample, $v_{\rm max} <%i$ $\pm 1\sigma$ ($\Lambda$CDM)'%vmcut) | |
| fiberd = cosmo_z_to_dist(zvol,disttype='angular')*(55/asecperrad)*1e3 #Mpc->kpc | |
| axvline(fiberd,zorder=1,linestyle='dashed',color='k')#,label=r'Fiber radius at $z_{\rm lim}$') | |
| xdyaxtrans = blended_transform_factory(gca().transData,gca().transAxes) | |
| text(fiberd+1,.95,r'Fiber radius at $z_{\rm lim}$',horizontalalignment='left',fontsize=16,transform=xdyaxtrans) | |
| l = legend(loc=4) | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| xlim(0,250) | |
| ylim(0,.026) | |
| ylabel(r'$f_{\rm bin}(d_{\rm proj})$') | |
| @plotfunc(mainfig=True) | |
| def fracradsamps(cutsd,hiicutsd,potsd,zsim,zvol,msiicone,vagcd,h): | |
| from astropysics.plotting import cumulative_plot | |
| from astropysics.coords import cosmo_z_to_dist | |
| from astropysics.constants import asecperrad | |
| from scipy.interpolate import interp1d | |
| redge = 250 | |
| cutsd = hiicutsd | |
| subratio = 4 | |
| ax1 = subplot2grid((subratio,1),(1,0),rowspan=subratio-1) | |
| ax2 = subplot2grid((subratio,1),(0,0),rowspan=1) | |
| subplots_adjust(left=.12,right=.96,top=.98,bottom=.08,hspace=0) | |
| axes(ax1) | |
| msiiinterps = [] | |
| #for i,(d,n) in enumerate(zip(msiicone['d_proj'],msiicone['n_iso'])): #old cone | |
| for i in range(msiicone.n): | |
| d = msiicone.get_masked_arr('rproj','isolated,deltav,conez',i) #Mpc/h | |
| iso = msiicone.get_masked_arr('isolated','conez',i) | |
| n = np.sum(iso) | |
| d = d[d<(redge/1000)]*1e3/h #now in kpc, not Mpc/h, cut on 250 kpc/h | |
| res = cumulative_plot(d, | |
| edges=(0,redge/h), | |
| frac=d.size/n, | |
| label='MS-II Projected' if i==0 else '', | |
| zorder=9,ls='-',c='g',lw=1,alpha=.4,doplot=False) | |
| msiiinterps.append(interp1d(*res)) | |
| msiid = linspace(0,redge/h,50) | |
| msiifs = array([[intr(x) for intr in msiiinterps] for x in msiid]) | |
| msiimeans = mean(msiifs,axis=-1) | |
| msiistds = std(msiifs,axis=-1) | |
| zsimana = zsim[zsim.nearest_lmc_analog] #comoving kpc/h | |
| zsimsep = zsimana.minsepkpch | |
| fibcoll = zsimana.fibercoll | |
| nds = cutsd['z500'].nd | |
| nds1 = cutsd['z500'].nd[cutsd['z500'].nn == 1] | |
| uniquenoz = cutsd['uniquenoz'] | |
| ndsnozcut = cutsd['nozcut'].nd | |
| ndsupper = np.concatenate((uniquenoz.nd,cutsd['z500'].nd)) | |
| ndsphot = np.concatenate((uniquenoz.nd,cutsd['nozcut'].nd)) | |
| perr = True | |
| #phot | |
| resp = cumulative_plot(ndsphot/h, | |
| edges=(0,redge/h), | |
| frac=ndsphot.size/len(potsd['pothosts']), | |
| label='Full Sample', | |
| zorder=10,c='r',ls=':',perr=perr) | |
| # #non-z objects only | |
| # cumulative_plot(uniquenoz.nd, | |
| # edges=(0,redge), | |
| # frac=uniquenoz.size/len(potsd['pothosts']), | |
| # label='Cloesest Sat w/ no $z$', | |
| # zorder=10,c='r',ls=':',perr=perr) | |
| # #spec | |
| ress = cumulative_plot(ndsnozcut/h, | |
| edges=(0,redge/h), | |
| frac=ndsnozcut.size/len(potsd['pothosts']), | |
| label='Any $z$ Sample', | |
| zorder=10,c='g',ls='--',perr=perr) | |
| # #phot - spec>500 | |
| # cumulative_plot(ndsupper, | |
| # edges=(0,redge), | |
| # frac=ndsupper.size/len(potsd['pothosts']), | |
| # label='Phot - Spec, $\Delta v > 500$ (upper limit)', | |
| # zorder=10,c='g',ls='--',perr=perr) | |
| # #spec<500 | |
| resdv = cumulative_plot(nds/h, | |
| edges=(0,redge/h), | |
| frac=nds.size/len(potsd['pothosts']), | |
| label='Clean Sample', | |
| zorder=10,c='b',perr=perr) | |
| #phot - spec>500 | |
| # plot(selfuncx,selffunccorr,'c-', | |
| # label='$\Delta v < 500$+Phot Sel Func (Best)',zorder=10) | |
| # plot(selfuncx,selffunccorr-selfuncbkgfunc(selfuncx,1500),':', | |
| # label='Bkg $r_0=1500$',zorder=10,lw=3) | |
| # plot(selfuncx,selffunccorr-selfuncbkgfunc(selfuncx,1000),':', | |
| # label='Bkg $r_0=1000$',zorder=10,lw=3) | |
| # plot(selfuncx,selffunccorr-selfuncbkgfunc(selfuncx,750),':', | |
| # label='Bkg $r_0=750$',zorder=10,lw=3) | |
| # plot(selfuncx,selffunccorr-selfuncbkgfunc(selfuncx,650),':', | |
| # label='Bkg $r_0=650$',zorder=10,lw=3) | |
| # fp = fill_between(msiid,msiimeans+msiistds,msiimeans-msiistds, | |
| # facecolor=(0.5,0.5,0.5,1),lw=0) | |
| # bar(-1,.1,facecolor=fp.get_facecolor()[0],lw=0,label=r'MS-II $1\sigma$') | |
| #55" is minimum seperation of fibers | |
| fiberd = cosmo_z_to_dist(zvol,disttype='angular')*(55/asecperrad)*1e3 #Mpc->kpc | |
| #axvline(fiberd,zorder=1,label='Fiber radius at z=%.4f'%zvol) | |
| legend(loc=4) | |
| #ylabel(r'$f(<d_{\rm proj}) \equiv N_{\rm sec}(<d_{\rm proj})/N_{\rm pri}(<d_{\rm proj})$') | |
| ylabel(r'$f_{\rm sat}(<d_{\rm proj})$') | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| ylim(0,.5) | |
| yticks([.1,.2,.3,.4]) | |
| xlim(0,250) | |
| xl = xlim() | |
| xtl = xticks()[0] | |
| #selection function | |
| axes(ax2) | |
| nds = cutsd['z500'].nd | |
| ndsspec = cutsd['nozcut'].nd | |
| nds1 = cutsd['z500'].nd[cutsd['z500'].nn == 1] | |
| npots = len(potsd['pothosts']) | |
| ndsuniquenoz = cutsd['uniquenoz'].nd | |
| ndsphot = np.concatenate((ndsuniquenoz,ndsspec)) | |
| selfuncx = sort(ndsphot)#kpc/h | |
| selfuncphotf = (arange(ndsphot.size)+1)/npots | |
| selfuncz500f = interp(selfuncx,sort(nds),(arange(nds.size)+1)/npots) | |
| selfuncspecf = interp(selfuncx,sort(ndsspec),(arange(ndsspec.size)+1)/npots) | |
| selfunc = selfuncz500f/selfuncspecf | |
| selfuncerr = selfunc*((selfuncz500f*npots)**-1+(selfuncspecf*npots)**-1)**0.5 | |
| m = selfuncx>min(nds/h) | |
| fill_between(selfuncx[m],selfunc[m]+selfuncerr[m],selfunc[m]-selfuncerr[m],color=(.75,.75,.75)) | |
| plot(selfuncx[m],selfunc[m],'k-') | |
| ylabel(r'$F_{\rm sel}(<d_{\rm proj})$') | |
| xlim(*xl) | |
| xticks(xtl,['' for xt in xtl]) | |
| ylim(0.25,.95) | |
| yticks([.3,.5,.7,.9]) | |
| @plotfunc | |
| def fracradselfunconly(cutsd,hiicutsd,potsd,zsim,zvol,msiicone,vagcd,h): | |
| from astropysics.plotting import cumulative_plot | |
| from astropysics.coords import cosmo_z_to_dist | |
| from astropysics.constants import asecperrad | |
| from scipy.interpolate import interp1d | |
| cutsd = hiicutsd | |
| redge = 250 | |
| msiiinterps = [] | |
| #for i,(d,n) in enumerate(zip(msiicone['d_proj'],msiicone['n_iso'])): #old cone | |
| for i in range(msiicone.n): | |
| d = msiicone.get_masked_arr('rproj','isolated,deltav,conez',i) #Mpc/h | |
| iso = msiicone.get_masked_arr('isolated','conez',i) | |
| n = np.sum(iso) | |
| d = d[d<(redge/1000)]*1e3/h #now in kpc, not Mpc/h, cut on 250 kpc/h | |
| res = cumulative_plot(d, | |
| edges=(0,redge/h), | |
| frac=d.size/n,doplot=False) | |
| msiiinterps.append(interp1d(*res)) | |
| msiid = linspace(0,redge/h,50) | |
| msiifs = array([[intr(x) for intr in msiiinterps] for x in msiid]) | |
| msiimeans = mean(msiifs,axis=-1) | |
| msiistds = std(msiifs,axis=-1) | |
| subplots_adjust(left=.15,right=.95,top=.95) | |
| zsimana = zsim[zsim.nearest_lmc_analog] #comoving kpc/h | |
| zsimsep = zsimana.minsepkpch | |
| fibcoll = zsimana.fibercoll | |
| nds = cutsd['z500'].nd | |
| nds1 = cutsd['z500'].nd[cutsd['z500'].nn == 1] | |
| uniquenoz = cutsd['uniquenoz'] | |
| ndsnozcut = cutsd['nozcut'].nd | |
| ndsupper = np.concatenate((uniquenoz.nd,cutsd['z500'].nd)) | |
| ndsphot = np.concatenate((uniquenoz.nd,cutsd['nozcut'].nd)) | |
| #compute the lower line + fraction of the photometric sample | |
| selfuncx = concatenate((nds,ndsnozcut)) | |
| selfuncx = sort(selfuncx[selfuncx<250]) | |
| selfuncz500 = interp(selfuncx,np.sort(nds),(arange(nds.size)+1)/len(potsd['pothosts'])) | |
| selfuncnozcut = interp(selfuncx,np.sort(ndsnozcut),(arange(ndsnozcut.size)+1)/len(potsd['pothosts'])) | |
| selfuncfrac = selfuncz500/selfuncnozcut | |
| selfuncunoz = interp(selfuncx,np.sort(uniquenoz.nd),(arange(uniquenoz.size)+1)/len(potsd['pothosts'])) | |
| selfuncphot = interp(selfuncx,np.sort(ndsphot),(arange(ndsphot.size)+1)/len(potsd['pothosts'])) | |
| selffunccorr = selfuncz500+selfuncunoz*selfuncfrac | |
| selffunccorr = selfuncphot*selfuncfrac | |
| #bkgx0 = 1500 | |
| selfuncbkgfunc = lambda x,x0:1-np.exp(-(x/x0)**2) | |
| #selffunccorr -= selfuncbkgfunc(selfuncx,bkgx0) | |
| perr = True | |
| #phot | |
| # resp = cumulative_plot(ndsphot, | |
| # edges=(0,redge), | |
| # frac=ndsphot.size/len(potsd['pothosts']), | |
| # label='Phot Sample', | |
| # zorder=10,c='r',ls=':',perr=perr) | |
| # #non-z objects only | |
| # cumulative_plot(uniquenoz.nd, | |
| # edges=(0,redge), | |
| # frac=uniquenoz.size/len(potsd['pothosts']), | |
| # label='Cloesest Sat w/ no $z$', | |
| # zorder=10,c='r',ls=':',perr=perr) | |
| # #spec | |
| # ress = cumulative_plot(ndsnozcut, | |
| # edges=(0,redge), | |
| # frac=ndsnozcut.size/len(potsd['pothosts']), | |
| # label='Spec Sample', | |
| # zorder=10,c='m',ls='-.',perr=perr) | |
| #phot - spec>500 | |
| cumulative_plot(ndsupper/h, | |
| edges=(0,redge/h), | |
| frac=ndsupper.size/len(potsd['pothosts']), | |
| #label='Phot - Spec, $\Delta v > 500$ (upper limit)', | |
| label='SDSS Upper', | |
| zorder=10,c='m',ls='-.',perr=perr) | |
| #phot - spec>500 | |
| plot(selfuncx/h,selffunccorr,'b-', | |
| label='SDSS Best',zorder=10) | |
| #label='$\Delta v < 500$+Phot Sel Func (Best)',zorder=10) | |
| #spec<500 | |
| resdv = cumulative_plot(nds/h, | |
| edges=(0,redge/h), | |
| frac=nds.size/len(potsd['pothosts']), | |
| #label='Spec, $\Delta v < 500$ (lower limit)', | |
| label='SDSS Lower', | |
| zorder=10,c='g',ls='--',perr=perr) | |
| plot(selfuncx/h,selfuncbkgfunc(selfuncx,1500),':', | |
| label='Bkg $r_0=1500$',zorder=10,lw=3) | |
| plot(selfuncx/h,selfuncbkgfunc(selfuncx,1000),':', | |
| label='Bkg $r_0=1000$',zorder=10,lw=3) | |
| plot(selfuncx/h,selfuncbkgfunc(selfuncx,750),':', | |
| label='Bkg $r_0=750$',zorder=10,lw=3) | |
| plot(selfuncx/h,selfuncbkgfunc(selfuncx,650),':', | |
| label='Bkg $r_0=650$',zorder=10,lw=3) | |
| legend(loc=0) | |
| # plot(selfuncx/h,selffunccorr-selfuncbkgfunc(selfuncx,1500),':', | |
| # label='Bkg $r_0=1500$',zorder=10,lw=3) | |
| # plot(selfuncx/h,selffunccorr-selfuncbkgfunc(selfuncx,1000),':', | |
| # label='Bkg $r_0=1000$',zorder=10,lw=3) | |
| # plot(selfuncx/h,selffunccorr-selfuncbkgfunc(selfuncx,750),':', | |
| # label='Bkg $r_0=750$',zorder=10,lw=3) | |
| # plot(selfuncx/h,selffunccorr-selfuncbkgfunc(selfuncx,650),':', | |
| # label='Bkg $r_0=650$',zorder=10,lw=3) | |
| fill_between(msiid,msiimeans+msiistds,msiimeans-msiistds, | |
| facecolor=(0.75,0.75,0.75,1),lw=0) | |
| yl = [0,.62] | |
| # #55" is minimum seperation of fibers | |
| # fiberd = cosmo_z_to_dist(zvol,disttype='angular')*(55/asecperrad)*1e3 #Mpc->kpc | |
| # vlines(fiberd,yl[0],yl[1],zorder=1,label='Fiber radius at z=%.4f'%zvol) | |
| #legend(loc=4) | |
| ylabel(r'$f(<d)$') | |
| xlabel(r'$d/{\rm kpc}$') | |
| #yticks([.1,.2,.3,.3,.45]) | |
| ylim(*yl) | |
| xlim(0,redge/h) | |
| #@plotfunc(figsize=(12,6)) | |
| @plotfunc(mainfig=True) | |
| def fracisocosmo(msiicone,h,LMCdm): | |
| from astropysics.plotting import cumulative_plot | |
| from astropysics.models import NFWModel | |
| from astropysics.phot import distance_from_modulus | |
| from scipy import stats | |
| m2r,v2m = NFWModel.Mvir_to_Rvir,NFWModel.Vvir_to_Mvir | |
| subplots_adjust(left=.15,top=.98,bottom=.08,right=.96) | |
| redge = 250 | |
| nbins = 5 | |
| rpobslike = msiicone.get_masked_arr('rproj','isolated,deltav,conez')*1e3/h | |
| vmobslike = msiicone.get_masked_arr('h_vmax','isolated,deltav,conez') | |
| r3obslike = msiicone.get_masked_arr('r3dto2d','isolated,deltav,conez')*1e3/h | |
| iso = msiicone.get_masked_arr('isolated','conez') | |
| rptrue = msiicone.get_masked_arr('rproj','central_host,conez')*1e3/h | |
| vmtrue = msiicone.get_masked_arr('h_vmax','central_host,conez') | |
| r3true = msiicone.get_masked_arr('r3dto2d','central_host,conez')*1e3/h | |
| # val,N = cumulative_plot(r3obslike,edges=(0,redge),frac=np.sum(r3obslike<250)/np.sum(iso), | |
| # label='Simulated Clean 3D',lw=4,ls='-',c='b',zorder=12) | |
| # LMCd = distance_from_modulus(LMCdm)/1000 | |
| # print 'Frac within 3D d',LMCd,':',interp(LMCd,val,N) | |
| cumulative_plot(rpobslike,edges=(0,redge),frac=np.sum(rpobslike<250)/np.sum(iso), | |
| label='Simulated Clean',lw=4,ls='--',c='g',zorder=11) | |
| cumulative_plot(r3true,edges=(0,redge),frac=np.sum(r3true<250)/r3true.size, | |
| label=r'3D Sat (Any $v_{\rm max,host}$)',lw=4,ls='-.',c='k',zorder=10) | |
| edges = [stats.scoreatpercentile(vmobslike,pc) for pc in linspace(0,100,nbins+1)[::-1]] | |
| cs = ('r','y','m','c','b') | |
| meds = [] | |
| for e1,e2,c in zip(edges[1:],edges[:-1],cs): | |
| msk = (e1<vmtrue)&(vmtrue<e2) | |
| meds.append(np.median(vmtrue[msk])) | |
| cumulative_plot(r3true[msk],edges=(0,redge),frac=np.sum(r3true[msk]<250)/np.sum(msk), | |
| label=r'%i km/s'%meds[-1],lw=1,ls='-',c=c) | |
| # cumulative_plot(r3true[msk],edges=(0,redge),frac=np.sum(r3true[msk]<250)/np.sum(msk), | |
| # label=r'$%i<v_{\rm max,host}<%i$ km/s'%(e1,e2),lw=1,ls='-',c=c) | |
| text(150,.96,'Simulation',fontsize=24,va='top',ha='left') | |
| #hand tweak these | |
| # text(195,.19,r'$v_{\rm max,host}\sim%i$'%meds[-1],fontsize=16,ha='left',va='top') | |
| # text(195,.235,r'$v_{\rm max,host}\sim%i$'%meds[-2],fontsize=16,ha='left',va='bottom') | |
| # text(195,.33,r'$v_{\rm max,host}\sim%i$'%meds[-3],fontsize=16,ha='left',va='top') | |
| # text(195,.53,r'$v_{\rm max,host}\sim%i$'%meds[-4],fontsize=16,ha='left',va='bottom') | |
| # text(195,.8,r'$v_{\rm max,host}\sim%i$'%meds[-5],fontsize=16,ha='left',va='top') | |
| legend(loc=2) | |
| ylabel(r'$f_{\rm sat}(<d_{\rm proj})$') | |
| ylim(0,1) | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| xlim(0,250) | |
| @plotfunc | |
| def fracradtrueiso(msiicone,h): | |
| from astropysics.plotting import cumulative_plot | |
| from astropysics.models import NFWModel | |
| m2r,v2m = NFWModel.Mvir_to_Rvir,NFWModel.Vvir_to_Mvir | |
| subplots_adjust(left=.15,top=.98,bottom=.08,right=.96) | |
| redge = 250 | |
| rp = msiicone.get_masked_arr('rproj','isolated,conez')*1e3/h | |
| rptruesat = msiicone.get_masked_arr('rproj','isolated,conez,truesat_fof')*1e3/h | |
| rpni = msiicone.get_masked_arr('rproj','conez,~isolated')*1e3/h | |
| rpnitruesat = msiicone.get_masked_arr('rproj','conez,truesat_fof,~isolated')*1e3/h | |
| rpa = msiicone.get_masked_arr('rproj','conez')*1e3/h | |
| rpatruesat = msiicone.get_masked_arr('rproj','conez,truesat_fof')*1e3/h | |
| massbins = linspace | |
| cumulative_plot(rp,frac=rptruesat.size/rp.size,edges=(0,redge),c='b',label='Isolated',ls='-') | |
| cumulative_plot(rpa,frac=rpatruesat.size/rpa.size,edges=(0,redge),c='g',label='Any',ls=':') | |
| cumulative_plot(rpni,frac=rpnitruesat.size/rpni.size,edges=(0,redge),c='r',label='Not Isolated',ls='--') | |
| text(15,.40,'Simulation',fontsize=24,va='top',ha='left') | |
| legend(loc=4) | |
| ylabel(r'$f_{\rm sat}(<d_{\rm proj})$') | |
| ylim(0,0.42) | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| xlim(0,redge) | |
| @plotfunc | |
| def fracradtruebin(msiicone,h): | |
| from astropysics.plotting import cumulative_plot | |
| from astropysics.models import NFWModel | |
| m2r,v2m = NFWModel.Mvir_to_Rvir,NFWModel.Vvir_to_Mvir | |
| subplots_adjust(left=.15,top=.98,bottom=.08,right=.96) | |
| redge = 250 | |
| rp = msiicone.get_masked_arr('rproj','isolated,conez')*1e3/h | |
| rptruesat = msiicone.get_masked_arr('rproj','isolated,conez,truesat_fof')*1e3/h | |
| vmh = msiicone.get_masked_arr('h_vmax','isolated,conez') | |
| vmhtruesat = msiicone.get_masked_arr('h_vmax','isolated,conez,truesat_fof') | |
| binedges = (650,400,250,225,200,185,166) | |
| cs = ['r','y','m','g','c','b'] | |
| for e1,e2,c in zip(binedges[1:],binedges[:-1],cs): | |
| msk = (e1<vmh)&(vmh<e2) | |
| mskts = (e1<vmhtruesat)&(vmhtruesat<e2) | |
| cumulative_plot(rptruesat[mskts],frac=rptruesat[mskts].size/rp[msk].size,edges=(0,redge), | |
| c=c,ls='-',label=r'$%i < v_{\rm max,host} < %i$'%(e1,e2)) | |
| cumulative_plot(rptruesat,frac=rptruesat.size/rp.size,edges=(0,redge), | |
| c='k',ls=':',label=r'All') | |
| text(100,.97,'Simulation',fontsize=24,va='top',ha='left') | |
| ylabel(r'$f_{\rm truesat}(<d_{\rm proj})$') | |
| ylim(0,1) | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| xlim(0,redge) | |
| legend(loc=0,prop={'size':14}) | |
| @plotfunc(mainfig=True) | |
| def fractrue(msiicone,h): | |
| #purity= #true from sample/#sel from sample | |
| subplots_adjust(left=.12,top=.98,bottom=.08,right=.98) | |
| lw = 3 | |
| dana = msiicone.get_masked_arr('rproj','isolated,conez,deltav')*1e3 | |
| sorti = argsort(dana) | |
| dana = dana[sorti] | |
| truesat = msiicone.get_masked_arr('truesat_fof','isolated,conez,deltav')[sorti] | |
| ftrue = cumsum(truesat)/(arange(dana.size)+1) | |
| plot(dana/h,ftrue,'b-',label='Isolated Primary, $\Delta v < 500$ (Clean)',zorder=20,lw=lw) | |
| dananoz = msiicone.get_masked_arr('rproj','isolated,conez')*1e3 | |
| sortinoz = argsort(dananoz) | |
| dananoz = dananoz[sortinoz] | |
| truesatnoz = msiicone.get_masked_arr('truesat_fof','isolated,conez')[sortinoz] | |
| ftruenoz = cumsum(truesatnoz)/(arange(dananoz.size)+1) | |
| plot(dananoz/h,ftruenoz,'g:',label='Isolated Primary, Any $\Delta v$ (Any $z$)',lw=lw) | |
| dananoisodv = msiicone.get_masked_arr('rproj','conez,deltav')*1e3 | |
| sortinoisodv = argsort(dananoisodv) | |
| dananoisodv = dananoisodv[sortinoisodv] | |
| truesatnoisodv = msiicone.get_masked_arr('truesat_fof','conez,deltav')[sortinoisodv] | |
| ftruenoisodv = cumsum(truesatnoisodv)/(arange(dananoisodv.size)+1) | |
| plot(dananoisodv/h,ftruenoisodv,'r--',label='Any Primary, $\Delta v < 500$',lw=lw) | |
| dananoiso = msiicone.get_masked_arr('rproj','conez')*1e3 | |
| sortinoiso = argsort(dananoiso) | |
| dananoiso = dananoiso[sortinoiso] | |
| truesatnoiso = msiicone.get_masked_arr('truesat_fof','conez')[sortinoiso] | |
| ftruenoiso = cumsum(truesatnoiso)/(arange(dananoiso.size)+1) | |
| plot(dananoiso/h,ftruenoiso,'m-.',label='Any Primary, Any $\Delta v$',lw=lw) | |
| text(340,.96,'Simulation',fontsize=24,va='top',ha='right') | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| #ylabel(r'$F_{\rm pure}(<d_{\rm proj})$') | |
| ylabel(r'Purity') | |
| xlim(0,250/h) | |
| legend(loc=0) | |
| ylim(0.2,1.0) | |
| @plotfunc | |
| def fractrue_cv(msiicone,h): | |
| #purity= #true from sample/#sel from sample | |
| #estimate for cosmic variance effect | |
| subplots_adjust(left=.12,top=.98,bottom=.08,right=.98) | |
| samples = arange(1,250,10) | |
| lw = 1 | |
| al = .5 | |
| labed = False | |
| for i in samples: | |
| dana = msiicone.get_masked_arr('rproj','isolated,conez,deltav',i=i)*1e3 | |
| sorti = argsort(dana) | |
| dana = dana[sorti] | |
| truesat = msiicone.get_masked_arr('truesat_fof','isolated,conez,deltav',i=i)[sorti] | |
| ftrue = cumsum(truesat)/(arange(dana.size)+1) | |
| if labed: | |
| lab = '' | |
| else: | |
| lab = 'Isolated Primary, $\Delta v < 500$ (Clean), x'+str(len(samples)) | |
| plot(dana/h,ftrue,'b-',label=lab,zorder=20,lw=lw,alpha=al) | |
| dananoiso = msiicone.get_masked_arr('rproj','conez',i=i)*1e3 | |
| sortinoiso = argsort(dananoiso) | |
| dananoiso = dananoiso[sortinoiso] | |
| truesatnoiso = msiicone.get_masked_arr('truesat_fof','conez',i=i)[sortinoiso] | |
| ftruenoiso = cumsum(truesatnoiso)/(arange(dananoiso.size)+1) | |
| if labed: | |
| lab = '' | |
| else: | |
| lab = 'Any Primary, Any $\Delta v$, x'+str(len(samples)) | |
| labed = True | |
| plot(dananoiso/h,ftruenoiso,'r--',label=lab,lw=lw,alpha=al) | |
| #text(340,.96,'Simulation',fontsize=24,va='top',ha='right') | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| ylabel(r'Purity') | |
| xlim(0,250/h) | |
| legend(loc='lower right') | |
| ylim(0.2,1.0) | |
| @plotfunc | |
| def fractrue2(msiicone,h): | |
| from scipy import stats | |
| #purity= #true from sample/#sel from sample | |
| subplots_adjust(left=.12,top=.98,bottom=.08,right=.98) | |
| lw = 3 | |
| dana = msiicone.get_masked_arr('rproj','isolated,conez,deltav')*1e3 | |
| sorti = argsort(dana) | |
| dana = dana[sorti] | |
| truesat = msiicone.get_masked_arr('truesat_fof','isolated,conez,deltav')[sorti] | |
| ftrue = cumsum(truesat)/(arange(dana.size)+1) | |
| plot(dana/h,ftrue,'b-',label='Isolated Primary, $\Delta v < 500$ (Clean)',zorder=20,lw=lw) | |
| dananoz = msiicone.get_masked_arr('rproj','isolated,conez')*1e3 | |
| sortinoz = argsort(dananoz) | |
| dananoz = dananoz[sortinoz] | |
| truesatnoz = msiicone.get_masked_arr('truesat_fof','isolated,conez')[sortinoz] | |
| ftruenoz = cumsum(truesatnoz)/(arange(dananoz.size)+1) | |
| plot(dananoz/h,ftruenoz,'g:',label='Isolated Primary, Any $\Delta v$ (Any $z$)',lw=lw) | |
| zs = np.zeros((msiicone.n,msiicone.arrs[0][0].size)).astype(bool) | |
| zs[123] = True | |
| msiicone.masks['zs'] = zs | |
| d = msiicone.get_masked_arr('dist_h','isolated,conez,zs') | |
| dananoz = msiicone.get_masked_arr('rproj','isolated,conez,zs')*1e3 | |
| sortinoz = argsort(dananoz) | |
| dananoz = dananoz[sortinoz] | |
| truesatnoz = msiicone.get_masked_arr('truesat_fof','isolated,conez,zs')[sortinoz] | |
| ftruenoz = cumsum(truesatnoz)/(arange(dananoz.size)+1) | |
| del msiicone.masks['zs'] | |
| dperc = percentile(d,50) | |
| dmsk = d<dperc | |
| print 'dperc',dperc,sum(dmsk),dmsk.size | |
| plot(dananoz[dmsk]/h,ftruenoz[dmsk],'r-.',label='Isolated Primary, Any $\Delta v$ (Any $z$) no cone',lw=lw) | |
| # dananoisodv = msiicone.get_masked_arr('rproj','conez,deltav')*1e3 | |
| # sortinoisodv = argsort(dananoisodv) | |
| # dananoisodv = dananoisodv[sortinoisodv] | |
| # truesatnoisodv = msiicone.get_masked_arr('truesat_fof','conez,deltav')[sortinoisodv] | |
| # ftruenoisodv = cumsum(truesatnoisodv)/(arange(dananoisodv.size)+1) | |
| # plot(dananoisodv/h,ftruenoisodv,'r--',label='Any Primary, $\Delta v < 500$',lw=lw) | |
| # dananoiso = msiicone.get_masked_arr('rproj','conez')*1e3 | |
| # sortinoiso = argsort(dananoiso) | |
| # dananoiso = dananoiso[sortinoiso] | |
| # truesatnoiso = msiicone.get_masked_arr('truesat_fof','conez')[sortinoiso] | |
| # ftruenoiso = cumsum(truesatnoiso)/(arange(dananoiso.size)+1) | |
| # plot(dananoiso/h,ftruenoiso,'m-.',label='Any Primary, Any $\Delta v$',lw=lw) | |
| text(340,.96,'Simulation',fontsize=24,va='top',ha='right') | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| #ylabel(r'$F_{\rm pure}(<d_{\rm proj})$') | |
| ylabel(r'Purity') | |
| xlim(0,250/h) | |
| legend(loc=0) | |
| ylim(0.2,1.0) | |
| @plotfunc | |
| def fractruesimple(msiicone,h): | |
| #purity= #true from sample/#sel from sample | |
| subplots_adjust(left=.12,top=.98,bottom=.08,right=.98) | |
| lw = 3 | |
| dana = msiicone.get_masked_arr('rproj','isolated,conez,deltav')*1e3 | |
| sorti = argsort(dana) | |
| dana = dana[sorti] | |
| truesat = msiicone.get_masked_arr('truesat_fof','isolated,conez,deltav')[sorti] | |
| ftrue = cumsum(truesat)/(arange(dana.size)+1) | |
| plot(dana/h,ftrue,'b-',label='Isolated Primary, $\Delta v < 500$ (Clean)',zorder=20,lw=lw) | |
| dananoz = msiicone.get_masked_arr('rproj','isolated,conez')*1e3 | |
| sortinoz = argsort(dananoz) | |
| dananoz = dananoz[sortinoz] | |
| truesatnoz = msiicone.get_masked_arr('truesat_fof','isolated,conez')[sortinoz] | |
| ftruenoz = cumsum(truesatnoz)/(arange(dananoz.size)+1) | |
| # plot(dananoz/h,ftruenoz,'g:',label='Isolated Primary, Any $\Delta v$ (Any $z$)',lw=lw) | |
| dananoisodv = msiicone.get_masked_arr('rproj','conez,deltav')*1e3 | |
| sortinoisodv = argsort(dananoisodv) | |
| dananoisodv = dananoisodv[sortinoisodv] | |
| truesatnoisodv = msiicone.get_masked_arr('truesat_fof','conez,deltav')[sortinoisodv] | |
| ftruenoisodv = cumsum(truesatnoisodv)/(arange(dananoisodv.size)+1) | |
| # plot(dananoisodv/h,ftruenoisodv,'r--',label='Any Primary, $\Delta v < 500$',lw=lw) | |
| dananoiso = msiicone.get_masked_arr('rproj','conez')*1e3 | |
| sortinoiso = argsort(dananoiso) | |
| dananoiso = dananoiso[sortinoiso] | |
| truesatnoiso = msiicone.get_masked_arr('truesat_fof','conez')[sortinoiso] | |
| ftruenoiso = cumsum(truesatnoiso)/(arange(dananoiso.size)+1) | |
| # plot(dananoiso/h,ftruenoiso,'m-.',label='Any Primary, Any $\Delta v$',lw=lw) | |
| plot(dananoiso/h,ftruenoiso,'r--',label='Any Primary, Any $\Delta v$',lw=lw) | |
| text(340,.96,'Simulation',fontsize=24,va='top',ha='right') | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| #ylabel(r'$F_{\rm pure}(<d_{\rm proj})$') | |
| ylabel(r'Purity $(N_{\rm true\,sats}\,/N_{\rm secondary-like})$') | |
| xlim(0,250/h) | |
| ylim(0.3,1.0) | |
| # legend(loc=0) | |
| @plotfunc | |
| def completeness(msiicone,h): | |
| subplots_adjust(left=.12,top=.98,bottom=.08,right=.98) | |
| d = np.sort(msiicone.get_masked_arr('rproj','truesat_3d_fof,isolated,conez')*1e3/h) | |
| plot(d,(arange(d.size)+1)/d.size,'b-',label='Isolated Primary') | |
| d = np.sort(msiicone.get_masked_arr('rproj','truesat_3d_fof,conez')*1e3/h) | |
| plot(d,(arange(d.size)+1)/d.size,'r--',label='Any Primary') | |
| text(340,.96,'Simulation',fontsize=24,va='top',ha='right') | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| ylabel(r'$F_{\rm complete}(<d_{\rm proj})$') | |
| xlim(0,250/h) | |
| legend(loc=0) | |
| ylim(0.2,1.0) | |
| @plotfunc | |
| def compcorr(msiicone,h): | |
| subplots_adjust(left=.12,top=.98,bottom=.08,right=.9) | |
| dpure = msiicone.get_masked_arr('rproj','isolated,conez,deltav')*1e3/h | |
| sorti = argsort(dpure) | |
| dpure = dpure[sorti] | |
| truesat = msiicone.get_masked_arr('truesat_fof','isolated,conez,deltav')[sorti] | |
| fpure = cumsum(truesat)/(arange(dpure.size)+1) | |
| dcomp = np.sort(msiicone.get_masked_arr('rproj','truesat_3d_fof,isolated,conez')*1e3/h) | |
| fcomp = (arange(dcomp.size)+1)/dcomp.size | |
| x = linspace(0,1000,1000) | |
| y = interp(x,dpure,fpure)/interp(x,dcomp,fcomp) | |
| ccerr = y*(1/(interp(x,dpure,fpure)*dpure.size)+1/(interp(x,dcomp,fcomp)*dcomp.size))**0.5 | |
| plot(x,interp(x,dpure,fpure),':b',label='Purity') | |
| plot(x,interp(x,dcomp,fcomp),'-.g',label='Completeness') | |
| plot(x,10*interp(x,dpure,fpure)/interp(x,dcomp,fcomp),'r-',label='P/C') | |
| legend(loc=5) | |
| ylim(0,1) | |
| text(340,.96,'Simulation',fontsize=24,va='top',ha='right') | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| ylabel(r'$F_P$ or $F_C$') | |
| twinx() | |
| plot(x,y,'r-',label='P/C',zorder=2) | |
| fill_between(x,y-ccerr,y+ccerr,color=(1,.5,.5)) | |
| axhline(1,linestyle='-',color='k',lw=1,zorder=1) | |
| ylim(0,9) | |
| yticks([1,2,3,4,5,6,7,8,9]) | |
| ylabel('$P/C$ Correction') | |
| xlim(0,250/h) | |
| @plotfunc | |
| def fracradcompcorr(msiicone,h,cutsd,hiicutsd,potsd): | |
| from astropysics.phot import distance_from_modulus | |
| subplots_adjust(left=.12,top=.98,bottom=.08,right=.9) | |
| #obs | |
| cutsd = hiicutsd #comment this to go back to un-visually filtered for HII regions | |
| subplots_adjust(left=.15,right=.95,top=.95) | |
| LMCd = distance_from_modulus(LMCdm)*h/1000 | |
| redge = 250 #in kpc/h | |
| yl = [0,.41] | |
| nds = cutsd['z500'].nd | |
| ndsspec = cutsd['nozcut'].nd | |
| nds1 = cutsd['z500'].nd[cutsd['z500'].nn == 1] | |
| npots = len(potsd['pothosts']) | |
| ndsuniquenoz = cutsd['uniquenoz'].nd | |
| ndsphot = np.concatenate((ndsuniquenoz,ndsspec)) | |
| selfuncx = sort(ndsphot)#kpc/h | |
| selfuncphotf = (arange(ndsphot.size)+1)/npots | |
| selfuncz500f = interp(selfuncx,sort(nds),(arange(nds.size)+1)/npots) | |
| selfuncspecf = interp(selfuncx,sort(ndsspec),(arange(ndsspec.size)+1)/npots) | |
| selfunccorrf = selfuncphotf*(selfuncz500f/selfuncspecf) | |
| obsx = selfuncx/h | |
| obsy = selfunccorrf | |
| obserr = selfunccorrf*((selfuncphotf*npots)**-1+(selfuncz500f*npots)**-1+(selfuncspecf*npots)**-1)**0.5 | |
| #compcorr | |
| dpure = msiicone.get_masked_arr('rproj','isolated,conez,deltav')*1e3/h | |
| sorti = argsort(dpure) | |
| dpure = dpure[sorti] | |
| truesat = msiicone.get_masked_arr('truesat_fof','isolated,conez,deltav')[sorti] | |
| fpure = cumsum(truesat)/(arange(dpure.size)+1) | |
| dcomp = np.sort(msiicone.get_masked_arr('rproj','truesat_3d_fof,isolated,conez')*1e3/h) | |
| fcomp = (arange(dcomp.size)+1)/dcomp.size | |
| #ccx = linspace(0,1000,1000) | |
| compcorr = interp(obsx,dpure,fpure)/interp(obsx,dcomp,fcomp) | |
| ccerr = compcorr*(1/(interp(obsx,dpure,fpure)*dpure.size)+1/(interp(obsx,dcomp,fcomp)*dcomp.size))**-0.5 | |
| plot(obsx,obsy,'r--',label='Selection Corrected $(<d_{\rm proj})$') | |
| plot(obsx,obsy*compcorr,'g-',label='FOF halo fraction') | |
| fill_between(obsx,(obsy-obserr)*compcorr,(obsy+obserr)*compcorr,color=(.5,.75,.5)) | |
| legend(loc=0) | |
| xlim(0,250/h) | |
| ylim(0,.8) | |
| ylabel(r'$f_{\rm sat}(<d_{\rm proj})$') | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| @plotfunc | |
| def F1(msiicone,h): | |
| #F1 = clean hosts w/ nearest 3D w/i dproj(dproj)/clean(dproj) | |
| subplots_adjust(left=.14,top=.98,bottom=.08,right=.98) | |
| rsamedv = sort(msiicone.get_masked_arr('r3dto2d','isolated,conez,deltav')) | |
| ralldv = sort(msiicone.get_masked_arr('rproj','isolated,conez,deltav')) | |
| rdv = sort(concatenate((rsamedv,ralldv))) | |
| fdv = interp(rdv,rsamedv,arange(rsamedv.size)+1)/interp(rdv,ralldv,arange(ralldv.size)+1) | |
| rsamendv = sort(msiicone.get_masked_arr('r3dto2d','isolated,conez')) | |
| rallndv = sort(msiicone.get_masked_arr('rproj','isolated,conez')) | |
| rndv = sort(concatenate((rsamendv,rallndv))) | |
| fndv = interp(rndv,rsamendv,arange(rsamendv.size)+1)/interp(rndv,rallndv,arange(rallndv.size)+1) | |
| rsameadv = sort(msiicone.get_masked_arr('r3dto2d','conez,deltav')) | |
| ralladv = sort(msiicone.get_masked_arr('rproj','conez,deltav')) | |
| radv = sort(concatenate((rsameadv,ralladv))) | |
| fadv = interp(radv,rsameadv,arange(rsameadv.size)+1)/interp(radv,ralladv,arange(ralladv.size)+1) | |
| rsameandv = sort(msiicone.get_masked_arr('r3dto2d','conez')) | |
| rallandv = sort(msiicone.get_masked_arr('rproj','conez')) | |
| randv = sort(concatenate((rsameandv,rallandv))) | |
| fandv = interp(randv,rsameandv,arange(rsameandv.size)+1)/interp(randv,rallandv,arange(rallandv.size)+1) | |
| plot(rdv*1000/h,fdv,'b-',label='Isolated Host, $\Delta v < 500$') | |
| plot(rndv*1000/h,fndv,'g:',label='Isolated Host, Any $\Delta v$') | |
| plot(radv*1000/h,fadv,'r--',label='Any Host, $\Delta v < 500$') | |
| plot(randv*1000/h,fandv,'m-.',label='Any Host, Any $\Delta v$') | |
| text(340,.96,'Simulation',fontsize=24,va='top',ha='right') | |
| legend(loc=4) | |
| xlabel(r'$d_{\rm proj}$') | |
| xlim(0,250/.7) | |
| ylabel(r'$F_1(<d_{\rm proj})$') | |
| ylim(0.5,1) | |
| @plotfunc | |
| def F2(msiicone,h): | |
| #F2 = 2dis3d(r)/clean(r) | |
| subplots_adjust(left=.14,top=.98,bottom=.08,right=.98) | |
| rsamedv = sort(msiicone.get_masked_arr('rproj','twodisthreed,isolated,conez,deltav')) | |
| ralldv = sort(msiicone.get_masked_arr('rproj','isolated,conez,deltav')) | |
| rdv = sort(concatenate((rsamedv,ralldv))) | |
| fdv = interp(rdv,rsamedv,arange(rsamedv.size)+1)/interp(rdv,ralldv,arange(ralldv.size)+1) | |
| rsamendv = sort(msiicone.get_masked_arr('rproj','twodisthreed,isolated,conez')) | |
| rallndv = sort(msiicone.get_masked_arr('rproj','isolated,conez')) | |
| rndv = sort(concatenate((rsamendv,rallndv))) | |
| fndv = interp(rndv,rsamendv,arange(rsamendv.size)+1)/interp(rndv,rallndv,arange(rallndv.size)+1) | |
| rsameadv = sort(msiicone.get_masked_arr('rproj','twodisthreed,conez,deltav')) | |
| ralladv = sort(msiicone.get_masked_arr('rproj','conez,deltav')) | |
| radv = sort(concatenate((rsameadv,ralladv))) | |
| fadv = interp(radv,rsameadv,arange(rsameadv.size)+1)/interp(radv,ralladv,arange(ralladv.size)+1) | |
| rsameandv = sort(msiicone.get_masked_arr('rproj','twodisthreed,conez')) | |
| rallandv = sort(msiicone.get_masked_arr('rproj','conez')) | |
| randv = sort(concatenate((rsameandv,rallandv))) | |
| fandv = interp(randv,rsameandv,arange(rsameandv.size)+1)/interp(randv,rallandv,arange(rallandv.size)+1) | |
| plot(rdv*1000/h,fdv,'b-',label='Isolated Host, $\Delta v < 500$') | |
| plot(rndv*1000/h,fndv,'g:',label='Isolated Host, Any $\Delta v$') | |
| plot(radv*1000/h,fadv,'r--',label='Any Host, $\Delta v < 500$') | |
| plot(randv*1000/h,fandv,'m-.',label='Any Host, Any $\Delta v$') | |
| text(340,.96,'Simulation',fontsize=24,va='top',ha='right') | |
| legend(loc=4) | |
| xlabel(r'$d_{\rm proj}$') | |
| xlim(0,250/.7) | |
| ylabel(r'$F_2(<d_{\rm proj})$') | |
| ylim(0.4,1) | |
| @plotfunc | |
| def compopure(msiicone,h): | |
| #C/P | |
| subplots_adjust(left=.14,top=.98,bottom=.08,right=.98) | |
| #completeness | |
| rsamedv = sort(msiicone.get_masked_arr('r3dto2d','isolated,conez,deltav')) | |
| ralldv = sort(msiicone.get_masked_arr('rproj','isolated,conez,deltav')) | |
| rdv = sort(concatenate((rsamedv,ralldv))) | |
| fdvc = interp(rdv,rsamedv,arange(rsamedv.size)+1)/interp(rdv,ralldv,arange(ralldv.size)+1) | |
| rsamendv = sort(msiicone.get_masked_arr('r3dto2d','isolated,conez')) | |
| rallndv = sort(msiicone.get_masked_arr('rproj','isolated,conez')) | |
| rndv = sort(concatenate((rsamendv,rallndv))) | |
| fndvc = interp(rndv,rsamendv,arange(rsamendv.size)+1)/interp(rndv,rallndv,arange(rallndv.size)+1) | |
| rsameadv = sort(msiicone.get_masked_arr('r3dto2d','conez,deltav')) | |
| ralladv = sort(msiicone.get_masked_arr('rproj','conez,deltav')) | |
| radv = sort(concatenate((rsameadv,ralladv))) | |
| fadvc = interp(radv,rsameadv,arange(rsameadv.size)+1)/interp(radv,ralladv,arange(ralladv.size)+1) | |
| rsameandv = sort(msiicone.get_masked_arr('r3dto2d','conez')) | |
| rallandv = sort(msiicone.get_masked_arr('rproj','conez')) | |
| randv = sort(concatenate((rsameandv,rallandv))) | |
| fandvc = interp(randv,rsameandv,arange(rsameandv.size)+1)/interp(randv,rallandv,arange(rallandv.size)+1) | |
| #purity | |
| dana = msiicone.get_masked_arr('rproj','isolated,conez,deltav')*1e3 | |
| sorti = argsort(dana) | |
| dana = dana[sorti] | |
| truesat = msiicone.get_masked_arr('truesat_fof','isolated,conez,deltav')[sorti] | |
| fpure = cumsum(truesat)/(arange(dana.size)+1) | |
| fdv = fdvc/np.interp(rdv*1000/h,dana,fpure) | |
| dananoz = msiicone.get_masked_arr('rproj','isolated,conez')*1e3 | |
| sortinoz = argsort(dananoz) | |
| dananoz = dananoz[sortinoz] | |
| truesatnoz = msiicone.get_masked_arr('truesat_fof','isolated,conez')[sortinoz] | |
| fpurenoz = cumsum(truesatnoz)/(arange(dananoz.size)+1) | |
| fndv = fndvc/np.interp(rndv*1000/h,dananoz,fpurenoz) | |
| dananoisodv = msiicone.get_masked_arr('rproj','conez,deltav')*1e3 | |
| sortinoisodv = argsort(dananoisodv) | |
| dananoisodv = dananoisodv[sortinoisodv] | |
| truesatnoisodv = msiicone.get_masked_arr('truesat_fof','conez,deltav')[sortinoisodv] | |
| fpurenoisodv = cumsum(truesatnoisodv)/(arange(dananoisodv.size)+1) | |
| fadv = fadvc/np.interp(radv*1000/h,dananoisodv,fpurenoisodv) | |
| dananoiso = msiicone.get_masked_arr('rproj','conez')*1e3 | |
| sortinoiso = argsort(dananoiso) | |
| dananoiso = dananoiso[sortinoiso] | |
| truesatnoiso = msiicone.get_masked_arr('truesat_fof','conez')[sortinoiso] | |
| fpurenoiso = cumsum(truesatnoiso)/(arange(dananoiso.size)+1) | |
| fandv = fandvc/np.interp(randv*1000/h,dananoiso,fpurenoiso) | |
| #plots | |
| plot(rdv*1000/h,fdv,'b-',label='Isolated Host, $\Delta v < 500$') | |
| #plot(rndv*1000/h,fndv,'g:',label='Isolated Host, Any $\Delta v$') | |
| #plot(radv*1000/h,fadv,'r--',label='Any Host, $\Delta v < 500$') | |
| #plot(randv*1000/h,fandv,'m-.',label='Any Host, Any $\Delta v$') | |
| text(15,1.23,'Simulation',fontsize=24,va='top',ha='left') | |
| #legend(loc=3) | |
| xlabel(r'$d_{\rm proj}$') | |
| xlim(0,250/.7) | |
| ylabel(r'$f_{\rm complete}/f_{\rm pure}(<d_{\rm proj})$') | |
| ylim(0.9,1.25) | |
| @plotfunc | |
| def fractruefofsph(msiicone): | |
| dana = msiicone.get_masked_arr('rproj','isolated,conez,deltav')*1e3 | |
| sorti = argsort(dana) | |
| dana = dana[sorti] | |
| truesat = msiicone.get_masked_arr('truesat_fof','isolated,conez,deltav')[sorti] | |
| ftrue = cumsum(truesat)/(arange(dana.size)+1) | |
| plot(dana,ftrue,'b-',label='FOF Isolated Host, $\Delta v < 500$',zorder=20) | |
| dananoz = msiicone.get_masked_arr('rproj','isolated,conez')*1e3 | |
| sortinoz = argsort(dananoz) | |
| dananoz = dananoz[sortinoz] | |
| truesatnoz = msiicone.get_masked_arr('truesat_fof','isolated,conez')[sortinoz] | |
| ftruenoz = cumsum(truesatnoz)/(arange(dananoz.size)+1) | |
| plot(dananoz,ftruenoz,'g-.',label='FOF Isolated Host, Any $\Delta v$') | |
| dana = msiicone.get_masked_arr('rproj','isolated,conez,deltav')*1e3 | |
| sorti = argsort(dana) | |
| dana = dana[sorti] | |
| truesat = msiicone.get_masked_arr('truesat_spherical','isolated,conez,deltav')[sorti] | |
| ftrue = cumsum(truesat)/(arange(dana.size)+1) | |
| plot(dana,ftrue,'r--',label='Spherically Isolated Host, $\Delta v < 500$',zorder=20) | |
| dananoz = msiicone.get_masked_arr('rproj','isolated,conez')*1e3 | |
| sortinoz = argsort(dananoz) | |
| dananoz = dananoz[sortinoz] | |
| truesatnoz = msiicone.get_masked_arr('truesat_spherical','isolated,conez')[sortinoz] | |
| ftruenoz = cumsum(truesatnoz)/(arange(dananoz.size)+1) | |
| plot(dananoz,ftruenoz,'m:',label='Spherically Isolated sph Host, Any $\Delta v$') | |
| legend(loc=0) | |
| xlabel(r'$d_{\rm proj}/h$') | |
| ylabel(r'$f_{\rm actualsat}(<d_{\rm proj})$') | |
| xlim(0,250) | |
| ylim(.55,.95) | |
| @plotfunc | |
| def fracnearesthost(msiicone,h): | |
| binwidth = .1 #Mpc | |
| ulim = 10 | |
| nh = msiicone.get_masked_arr('d_3d_closest_host','conez')/h | |
| nhi = msiicone.get_masked_arr('d_3d_closest_host','conez,isolated')/h | |
| nhni = msiicone.get_masked_arr('d_3d_closest_host','conez,~isolated')/h | |
| if ulim is None: | |
| ulim = nh.max() | |
| xbs = linspace(nh.min(),ulim,250) | |
| xs = [] | |
| fs = [] | |
| for x,x0,x1 in zip(xbs,xbs-binwidth/2,xbs+binwidth/2): | |
| n = sum((x0<nh)&(nh<x1)) | |
| ni = sum((x0<nhi)&(nhi<x1)) | |
| if n!=0: | |
| xs.append(x) | |
| fs.append(float(ni)/float(n)) | |
| xs,fs = array(xs),array(fs) | |
| plot(xs,fs) | |
| axvline(.778,c='k') | |
| print 'LG-like f=',interp(.778,xs,fs) | |
| print '50% d=',interp(.5,sort(fs),xs[argsort(fs)]) | |
| xlabel(r'$d_{\rm nn,hosts}$') | |
| ylabel(r'$f_{\rm iso}$') | |
| xlim(0,ulim) | |
| @plotfunc | |
| def fractruehybrid(zsim): | |
| zsimana = zsim[zsim.nearest_lmc_analog] #comoving kpc/h | |
| sorti = argsort(zsimana.minsepkpch) | |
| zsimsep = zsimana.minsepkpch[sorti] | |
| zsimftrue = cumsum(zsimana.actual_sat[sorti])/(arange(zsimana.size)+1) | |
| plot(zsimsep,zsimftrue,'b-',label='LMC Analog',zorder=20) | |
| zsimana = zsim[zsim.nearest_lmc_analog_noz] #comoving kpc/h | |
| sorti = argsort(zsimana.minsepkpch) | |
| zsimsep = zsimana.minsepkpch[sorti] | |
| zsimftrue = cumsum(zsimana.actual_sat[sorti])/(arange(zsimana.size)+1) | |
| plot(zsimsep,zsimftrue,'y:',label='No $\Delta v$') | |
| zsimana = zsim[zsim.nearest_lmc_analog_noiso] #comoving kpc/h | |
| sorti = argsort(zsimana.minsepkpch) | |
| zsimsep = zsimana.minsepkpch[sorti] | |
| zsimftrue = cumsum(zsimana.actual_sat[sorti])/(arange(zsimana.size)+1) | |
| plot(zsimsep,zsimftrue,'r--',label='Non-isolated host') | |
| legend(loc=0) | |
| xlabel(r'$f(<d_{\rm proj})$') | |
| ylabel(r'$f_{\rm actualsat}$') | |
| @plotfunc | |
| def truesatdvdd(msiicone): | |
| subplots_adjust(left=.15,top=.97,bottom=.08,right=.95) | |
| m = msiicone.conez.ravel() | |
| dv = msiicone.delta_v.ravel()[m] | |
| dd = (msiicone.dist_h - msiicone.dist_s).ravel()[m] | |
| tfof = msiicone.truesat_fof.ravel()[m] | |
| tsph = msiicone.truesat_spherical.ravel()[m] | |
| both = tfof&tsph | |
| ns = (~tfof)&(~tsph) | |
| plot(dv[both],dd[both],'g.',ms=1,label='Both Satellites',zorder=2) | |
| plot(dv[tfof&(~both)],dd[tfof&(~both)],'b.',ms=1,label='Only FOF Satellites',zorder=4) | |
| plot(dv[tsph&(~both)],dd[tsph&(~both)],'c.',ms=1,label='Only Spherical Satellites',zorder=3) | |
| plot(dv[ns],dd[ns],'r.',ms=1,label='Non-Satellites',zorder=1) | |
| axvline(500,c='k',zorder=5) | |
| axvline(-500,c='k',zorder=5) | |
| legend(loc=4,markerscale=8) | |
| xlabel(r'$\Delta v$') | |
| xlim(-10000,10000) | |
| ylabel(r'$(d_{\rm pri} - d_{\rm sec} )/ ({\rm Mpc} h^{-1} )$') | |
| ylim(-100,100) | |
| @plotfunc | |
| def isodvdd(msiicone): | |
| subplots_adjust(left=.15,top=.97,bottom=.08,right=.95) | |
| m = msiicone.conez.ravel() | |
| dv = msiicone.delta_v.ravel()[m] | |
| dd = (msiicone.dist_h - msiicone.dist_s).ravel()[m] | |
| tfof = msiicone.truesat_fof.ravel()[m] | |
| iso = msiicone.isolated.ravel()[m] | |
| ns = ~tfof | |
| plot(dv[iso&tfof],dd[iso&tfof],'g.',ms=1,label='Satellites of Isolated Hosts',zorder=4) | |
| plot(dv[(~iso)&tfof],dd[(~iso)&tfof],'b.',ms=1,label='Satellites of Non-Isolated Hosts',zorder=3) | |
| plot(dv[ns],dd[ns],'m.',ms=1,label='Non-Satellites of Isolated Hosts',zorder=2,alpha=.8) | |
| plot(dv[ns],dd[ns],'r.',ms=1,label='Non-Satellites of Non-Isolated Hosts',zorder=1,alpha=.8) | |
| axvline(500,c='k',zorder=5) | |
| axvline(-500,c='k',zorder=5) | |
| legend(loc=4,markerscale=8) | |
| xlabel(r'$\Delta v$') | |
| xlim(-10000,10000) | |
| ylabel(r'$(d_{\rm pri} - d_{\rm sec} )/ ({\rm Mpc} h^{-1} )$') | |
| ylim(-100,100) | |
| @plotfunc | |
| def isovelfunc(zsim): | |
| from astropysics.plotting import plot_histogram | |
| mj = zsim[zsim.major] | |
| mji = mj[mj.isolated] | |
| mjni = mj[~mj.isolated] | |
| fiso = np.sum(mji.id==mji.id_host)/mji.size | |
| fniso = np.sum(mjni.id==mjni.id_host)/mjni.size | |
| fall = np.sum(mj.id==mj.id_host)/mj.size | |
| #vmin = np.min(mj.vin_host) | |
| #for some reason there's a single inner outlier, so go to second... | |
| vmin = mj.vin_host[np.argsort(mj.vin_host)[1]] | |
| vmax = 600 | |
| nb = 50 | |
| b = linspace(vmin,vmax,nb) | |
| n = histogram(mj.vin_host,bins=b)[0] | |
| ni = histogram(mji.vin_host,bins=b)[0] | |
| nni = histogram(mjni.vin_host,bins=b)[0] | |
| plot_histogram(nni/np.max(nni),b,color='r',ls='--',label=r'Non-Isolated $f_{\rm host}=%.2f$'%fniso) | |
| plot_histogram(n/np.max(n),b,color='g',ls=':',label=r'All $f_{\rm host}=%.2f$'%fall) | |
| plot_histogram(ni/np.max(ni),b,color='b',label=r'Isolated $f_{\rm host}=%.2f$'%fiso) | |
| legend(loc=0) | |
| xlabel(r'$v_{\rm in}$') | |
| xlim(vmin,vmax) | |
| ylabel(r'$N/N_{\rm max}$') | |
| ylim(0,1.02) | |
| @plotfunc | |
| def simcorr(varvals_sim,varnames_sim,varobslast): | |
| subplots_adjust(left=.28,bottom=.25) | |
| jet() | |
| n = len(varnames_sim) | |
| cc = corrcoef(varvals_sim) | |
| imshow(cc,interpolation='nearest',norm=Normalize(-1,1)) | |
| hlines(varobslast+.5,-0.5,n-.5) | |
| vlines(varobslast+.5,-0.5,n-.5) | |
| cb = colorbar() | |
| cb.set_label('Correlation Coefficient') | |
| xticks(arange(n),varnames_sim,rotation=90) | |
| yticks(arange(n),varnames_sim) | |
| ylim(-0.5,n-.5) | |
| @plotfunc | |
| def simcorrbw(varvals_sim,varnames_sim,varobslast): | |
| subplots_adjust(left=.28,bottom=.25) | |
| gray() | |
| n = len(varnames_sim) | |
| acc = abs(corrcoef(varvals_sim)) | |
| imshow(acc,interpolation='nearest',norm=Normalize(0,1)) | |
| hlines(varobslast+.5,-0.5,n-.5) | |
| vlines(varobslast+.5,-0.5,n-.5) | |
| cb = colorbar() | |
| cb.set_label('Absolute Correlation Coefficient') | |
| xticks(arange(n),varnames_sim,rotation=90) | |
| yticks(arange(n),varnames_sim) | |
| ylim(-0.5,n-.5) | |
| def tilted_ks(d1,frac1,d2,frac2): | |
| from scipy.stats import ks_2samp,ksprob | |
| n1,n2 = len(d1),len(d2) | |
| if frac1>frac2: | |
| nadd = int(round(d2.size*(frac1/frac2-1))) | |
| d1s = d1 | |
| d2s = np.concatenate((d2,ones(nadd)*d2.max())) | |
| else: | |
| nadd = int(round(d1.size*(frac2/frac1-1))) | |
| d1s = np.concatenate((d1,ones(nadd)*d1.max())) | |
| d2s = d2 | |
| D,p = ks_2samp(d1s,d2s) | |
| en = (n1*n2/(n1+n2))**0.5 | |
| try: | |
| return D,ksprob((en+0.12+0.11/en)*D) #why not ksprob(en*D)? | |
| except: | |
| return D,1.0 | |
| @plotfunc | |
| def sepcolor(sep,us,gs,rs,i_s,z_s,Js,Hs,Ks): | |
| from astropysics.plotting import cumulative_plot | |
| cs = gs-rs | |
| cedge = [0,.5,.64,1] | |
| sepedge = 250 | |
| mb,mg,mr = [(cedge[i]<cs)&(cs<cedge[i+1]) for i in range(3)] | |
| redlabel = '%.2f<g-r<%.2f'%tuple(cedge[2:4])+', N=%i'%sum(sep[mr]<=sepedge) | |
| greenlabel = '%.2f<g-r<%.2f'%tuple(cedge[1:3])+', N=%i'%sum(sep[mg]<=sepedge) | |
| bluelabel = '%.2f<g-r<%.2f'%tuple(cedge[0:2])+', N=%i'%sum(sep[mb]<=sepedge) | |
| # bins=linspace(0,sepedge,11) | |
| # hist(sep[mr],histtype='step',label=redlabel,bins=bins,color='r',normed=True) | |
| # hist(sep[mg],histtype='step',label=greenlabel,bins=bins,color='g',normed=True) | |
| # hist(sep[mb],histtype='step',label=bluelabel,bins=bins,color='b',normed=True) | |
| cumulative_plot(sep[mr],label=redlabel,color='r',edges=(0,sepedge),frac=True) | |
| cumulative_plot(sep[mg],label=greenlabel,color='g',edges=(0,sepedge),frac=True) | |
| cumulative_plot(sep[mb],label=bluelabel,color='b',edges=(0,sepedge),frac=True) | |
| xlabel('$r_{\\rm sep}/({\\rm kpc}/h)$') | |
| ylabel('$f(<r)$') | |
| legend(loc=2) | |
| @plotfunc | |
| def sepcolorfrac(cutsd,potsd,vagcd): | |
| from astropysics.plotting import cumulative_plot | |
| u,g,r,i,z,J,H,K = vagcd['kcorr'].ABSMAG.T | |
| sep = cutsd['z500'].nd | |
| pinds = potsd['potsats'] | |
| sinds = cutsd['z500'].ni | |
| cs = g[sinds]-r[sinds] | |
| cp = g[pinds]-r[pinds] | |
| cedges = [0,.5,.64,1] | |
| cmid = (cedges[1]+cedges[2])/2 | |
| Nall = sum((cedges[0]<cp)&(cp<cedges[-1])) | |
| #red/all including green valley galaxies | |
| Nr = sum((cedges[2]<cp)&(cp<cedges[3])) | |
| Nb = sum((cedges[0]<cp)&(cp<cedges[1])) | |
| rbpots = Nr/Nall | |
| rbpe = (rbpots**2/Nr+rbpots**2/Nall)**0.5 | |
| brpots = Nb/Nall | |
| brpe = (brpots**2/Nb+brpots**2/Nall)**0.5 | |
| #red/all with green valley removed from sample | |
| Nrg = sum((cmid<cp)&(cp<cedges[3])) | |
| rbgpots = Nrg/Nall | |
| rbgpe = (rbgpots**2/Nrg+rbgpots**2/Nall)**0.5 | |
| brgpots = 1/rbgpots | |
| brgpe = rbgpe*brgpots**2 | |
| #green/all ratio | |
| Ng = sum((cedges[1]<cp)&(cp<cedges[2])) | |
| gpots = Ng/Nall | |
| gpe = (gpots**2/Ng+gpots**2/Nall)**0.5 | |
| sorti = argsort(sep)#small->large | |
| sepsort = sep[sorti] | |
| Nall = cumsum((cedges[0]<cs[sorti])&(cs[sorti]<cedges[-1])) | |
| Nr = cumsum((cedges[2]<cs[sorti])&(cs[sorti]<cedges[3])) | |
| Nb = cumsum((cedges[0]<cs[sorti])&(cs[sorti]<cedges[1])) | |
| rbsats = Nr/Nall | |
| rbse = (rbsats**2/Nr + rbsats**2/Nall)**0.5 | |
| rbf = rbsats/rbpots | |
| rbfe = ((rbse/rbpots)**2 + (rbf*rbpe/rbpots)**2)**0.5 | |
| brsats = Nb/Nall | |
| brse = (brsats**2/Nb + brsats**2/Nall)**0.5 | |
| brf = brsats/brpots | |
| brfe = ((brse/brpots)**2 + (brf*brpe/brpots)**2)**0.5 | |
| rbfe[~isfinite(rbfe)] = 0 | |
| brfe[~isfinite(brfe)] = 0 | |
| Nrg = cumsum((cmid<cs[sorti])&(cs[sorti]<cedges[3])) | |
| rbgsats = Nrg/Nall | |
| rbgf = rbgsats/rbgpots | |
| rbgse = (rbgsats**2/Nrg+rbgsats**2/Nall)**0.5 | |
| rbgfe = ((rbgse/rbgpots)**2 + (rbgf*rbgpe/rbgpots)**2)**0.5 | |
| brgsats = 1/rbgsats | |
| brgse = rbgse*brgsats**2 | |
| brgf = brgsats/brgpots | |
| brgfe = ((brgse/brgpots)**2 + (brgf*brgpe/brgpots)**2)**0.5 | |
| rbgfe[~isfinite(rbgfe)] = 0 | |
| brgfe[~isfinite(brgfe)] = 0 | |
| Ng = cumsum((cedges[1]<cs[sorti])&(cs[sorti]<cedges[2])) | |
| gsats = Ng/Nall | |
| gse = (gsats**2/Ng+gsats**2/Nall)**0.5 | |
| gf = gsats/gpots | |
| gfe = ((gse/gpots)**2 + (gf*gpe/gpots)**2)**0.5 | |
| gfe[~isfinite(gfe)] = 0 | |
| plot(insert(sepsort,0,0),insert(rbgf,0,0),'r',drawstyle='steps',label='$f=N_r/N_{\\rm all}$ including green') | |
| #fill_between(repeat(sepsort,2)[1:],repeat(rbgf-rbgfe,2)[:-1],repeat(rbgf+rbgfe,2)[:-1],color=(.8,.7,.7)) | |
| plot(insert(sepsort,0,0),insert(brgf,0,0),'b',drawstyle='steps',label='$f=N_b/N_{\\rm all}$ including green') | |
| #fill_between(repeat(sepsort,2)[1:],repeat(brgf-brgfe,2)[:-1],repeat(brgf+brgfe,2)[:-1],color=(.7,.7,.8)) | |
| plot(insert(sepsort,0,0),insert(rbf,0,0),'r:',drawstyle='steps',label='$f=N_r/N_{\\rm all}$') | |
| fill_between(repeat(sepsort,2)[1:],repeat(rbf-rbfe,2)[:-1],repeat(rbf+rbfe,2)[:-1],color=(.8,.7,.7)) | |
| plot(insert(sepsort,0,0),insert(brf,0,0),'b:',drawstyle='steps',label='$f=N_b/N_{\\rm all}$') | |
| fill_between(repeat(sepsort,2)[1:],repeat(brf-brfe,2)[:-1],repeat(brf+brfe,2)[:-1],color=(.7,.7,.8)) | |
| plot(insert(sepsort,0,0),insert(gf,0,0),'g:',drawstyle='steps',label='$f=N_g/N_{\\rm all}$') | |
| fill_between(repeat(sepsort,2)[1:],repeat(gf-gfe,2)[:-1],repeat(gf+gfe,2)[:-1],color=(.7,.8,.7)) | |
| legend(loc=0) | |
| hlines(1,0,250,'k','dashed',lw=1) | |
| xlim(0,250) | |
| ylim(0,1.75) | |
| xlabel('$r_{\\rm sep}/({\\rm kpc}/h)$') | |
| ylabel('$f_{\\rm sat}(<r)/f_{\\rm field}$') | |
| @plotfunc | |
| def sepfracmstar(cutsd,potsd,vagcd): | |
| from astropysics.plotting import cumulative_plot | |
| u,g,r,i,z,J,H,K = vagcd['kcorr'].ABSMAG.T | |
| lmstar = log10(vagcd['kcorr'].MASS) | |
| sep = cutsd['z500'].nd | |
| pinds = potsd['potsats'] | |
| sinds = cutsd['z500'].ni | |
| cs = g[sinds]-r[sinds] | |
| cp = g[pinds]-r[pinds] | |
| cedges = [0,.5,.64,1] | |
| cmid = (cedges[1]+cedges[2])/2 | |
| lmstarssort = lmstar[sinds] | |
| lmstarssort.sort() | |
| logmedge = [2*lmstarssort[0]-lmstarssort[1], | |
| lmstarssort[int(lmstarssort.size/3)], | |
| lmstarssort[int(lmstarssort.size*2/3)], | |
| 2*lmstarssort[-1]-lmstarssort[-2]] | |
| colors = ['b','g','r'] | |
| styles = ['-','--',':'] | |
| print logmedge | |
| for color,sty,e1,e2 in zip(colors,styles,logmedge[:-1],logmedge[1:]): | |
| label = '$10^{%.1f}< M_{*{\\rm sat}} < 10^{%.1f}$'%(e1,e2) | |
| ms = (e1<lmstar[sinds])&(lmstar[sinds]<e2) | |
| mp = (e1<lmstar[pinds])&(lmstar[pinds]<e2) | |
| sorti = argsort(sep[ms])#small->large | |
| sepsort = sep[ms][sorti] | |
| rbsats = cumsum((cmid<cs[ms][sorti])&(cs[ms][sorti]<cedges[3])) / \ | |
| cumsum((cedges[0]<cs[ms][sorti])&(cs[ms][sorti]<cedges[-1])) | |
| rbpots = sum((cmid<cp[mp])&(cp[mp]<cedges[3]))/sum((cedges[0]<cp[mp])&(cp[mp]<cedges[-1])) | |
| plot(insert(sepsort,0,0),insert(rbsats,0,0)/rbpots,color=color,linestyle=sty,drawstyle='steps',label=label) | |
| legend(loc=0) | |
| hlines(1,0,250,'k','dashed',lw=1) | |
| xlim(0,250) | |
| ylim(.9,2.5) | |
| xlabel('$r_re{\\rm sep}/({\\rm kpc}/h)$') | |
| ylabel('$f_{r_{\\rm sat}}(<r)/f_{r_{\\rm field}}$') | |
| @plotfunc | |
| def sepssfr300(sep,b3s): | |
| from astropysics.plotting import cumulative_plot | |
| cs = log10(b3s) | |
| cs[argsort(cs)][cs.size/3] | |
| cedge = [cs.min(),cs[argsort(cs)][cs.size/3],cs[argsort(cs)][2*cs.size/3],cs.max()] | |
| sepedge = 250 | |
| mr,mg,mb = [(cedge[i]<=cs)&(cs<=cedge[i+1]) for i in range(3)] | |
| bluelabel = r'$%.4f<\log(b300)<%.4f$'%tuple(cedge[2:4])+', N=%i'%sum(sep[mr]<=sepedge) | |
| greenlabel = r'$%.4f<\log(b300)<%.4f$'%tuple(cedge[1:3])+', N=%i'%sum(sep[mg]<=sepedge) | |
| redlabel = r'$%.4f<\log(b300)<%.4f$'%tuple(cedge[0:2])+', N=%i'%sum(sep[mb]<=sepedge) | |
| cumulative_plot(sep[mr],label=redlabel,color='r',edges=(0,sepedge),frac=True) | |
| cumulative_plot(sep[mg],label=greenlabel,color='g',edges=(0,sepedge),frac=True) | |
| cumulative_plot(sep[mb],label=bluelabel,color='b',edges=(0,sepedge),frac=True) | |
| xlabel('$r_{\\rm sep}/({\\rm kpc}/h)$') | |
| ylabel('$f(<r)$') | |
| legend(loc=2) | |
| @plotfunc | |
| def sepssfr1000(sep,b10s): | |
| from astropysics.plotting import cumulative_plot | |
| cs = log10(b10s) | |
| cs[argsort(cs)][cs.size/3] | |
| cedge = [cs.min(),cs[argsort(cs)][cs.size/3],cs[argsort(cs)][2*cs.size/3],cs.max()] | |
| sepedge = 250 | |
| mr,mg,mb = [(cedge[i]<=cs)&(cs<=cedge[i+1]) for i in range(3)] | |
| bluelabel = r'$%.4f<\log(b1000)<%.4f$'%tuple(cedge[2:4])+', N=%i'%sum(sep[mr]<=sepedge) | |
| greenlabel = r'$%.4f<\log(b1000)<%.4f$'%tuple(cedge[1:3])+', N=%i'%sum(sep[mg]<=sepedge) | |
| redlabel = r'$%.4f<\log(b1000)<%.4f$'%tuple(cedge[0:2])+', N=%i'%sum(sep[mb]<=sepedge) | |
| cumulative_plot(sep[mr],label=redlabel,color='r',edges=(0,sepedge),frac=True) | |
| cumulative_plot(sep[mg],label=greenlabel,color='g',edges=(0,sepedge),frac=True) | |
| cumulative_plot(sep[mb],label=bluelabel,color='b',edges=(0,sepedge),frac=True) | |
| xlabel('$r_{\\rm sep}/({\\rm kpc}/h)$') | |
| ylabel('$f(<r)$') | |
| legend(loc=2) | |
| @plotfunc | |
| def Dnnenc(zsim): | |
| anas = zsim[zsim.lmc_analog] | |
| m0 = anas.n_encounter==0 | |
| m1 = anas.n_encounter==1 | |
| m2 = anas.n_encounter==2 | |
| m3 = anas.n_encounter>2 | |
| # m2 = (2<=anas.n_encounter)&(anas.n_encounter<=4) | |
| # m3 = 4<anas.n_encounter | |
| bins=linspace(0,100,11) | |
| hist(zsim.minsepkpch[m0],histtype='step',label='First Infall',bins=bins,normed=True) | |
| hist(zsim.minsepkpch[m1],histtype='step',label='$N_{\\rm enc} = 1$',bins=bins,normed=True) | |
| hist(zsim.minsepkpch[m2],histtype='step',label='$N_{\\rm enc} = 2$',bins=bins,normed=True) | |
| hist(zsim.minsepkpch[m3],histtype='step',label='$N_{\\rm enc} > 2$',bins=bins,normed=True) | |
| xlabel('$D_n$') | |
| #ylabel('N') | |
| legend(loc=2) | |
| @plotfunc | |
| def Dntpass(zsim): | |
| anas = zsim[zsim.lmc_analog] | |
| m0 = anas.t_pass < 0.5 | |
| m1 = (0.5 < anas.t_pass)&(anas.t_pass < 1) | |
| m2 = (1 < anas.t_pass)&(anas.t_pass < 2) | |
| m3 = anas.t_pass > 2 | |
| # m2 = (2<=anas.n_encounter)&(anas.n_encounter<=4) | |
| # m3 = 4<anas.n_encounter | |
| bins=linspace(0,100,11) | |
| hist(zsim.minsepkpch[m0],histtype='step',label=r'$t_{\rm pass} < 0.5 \; {\rm Gyr}$',bins=bins,normed=True) | |
| hist(zsim.minsepkpch[m1],histtype='step',label=r'$0.5<t_{\rm pass}<1 \; {\rm Gyr}$',bins=bins,normed=True) | |
| hist(zsim.minsepkpch[m2],histtype='step',label=r'$1<t_{\rm pass}<2 \; {\rm Gyr}$',bins=bins,normed=True) | |
| hist(zsim.minsepkpch[m3],histtype='step',label=r'$2<t_{\rm pass} \; {\rm Gyr}$',bins=bins,normed=True) | |
| xlabel('$D_n$') | |
| #ylabel('N') | |
| legend(loc=2) | |
| @plotfunc(figsize=(8,4)) | |
| def blast(zsim): | |
| subplots_adjust(bottom=.18) | |
| zana = zsim[zsim.lmc_analog] | |
| x,y = zana.minsepkpch,zana.b_encounter_last | |
| edges = [0,20,50,100,150,250] | |
| ms = [(e1<x)&(x<e2) for e1,e2 in zip(edges[:-1],edges[1:])] | |
| mns=array([(mean(x[m]),mean(y[m])) for m in ms]).T | |
| stds=array([(std(x[m]),std(y[m])) for m in ms]).T | |
| errorbar(mns[0],mns[1],stds[1],stds[0],'o') | |
| ylim(0,160) | |
| #errorbar(mns[0],mns[1],None,stds[0],'o') | |
| #ylim(40,100) | |
| xlabel('$d_{\\rm proj}/(\\rm kpc/h)$') | |
| ylabel('$r_{\\rm last}/(\\rm kpc/h)$') | |
| @plotfunc | |
| def emmdist(liness,scutsi,cutsd): | |
| from astropysics.utils import median_absolute_deviation,biweight_midvariance,interquartile_range | |
| d = cutsd['z500'].nd | |
| si = cutsd['z500'].ni | |
| li = liness.VAGCidsat | |
| sortedvagcinds = argsort(si) | |
| si = si[sortedvagcinds] | |
| d = d[sortedvagcinds] | |
| dli = d[np.searchsorted(si,li)] | |
| linenames = ('Ha_ew','OII_ew','OIII_ew') | |
| colors = ('r','b','g') | |
| for linename,color in zip(linenames,colors): | |
| currliness = liness.field(linename) | |
| msk = currliness>-9999 | |
| scatter(dli[msk],currliness[msk],color=color,label=linename) | |
| xlim(0,250) | |
| ylim(-50,185) | |
| xlabel('$d_{\\rm proj}/(\\rm kpc/h)$') | |
| ylabel('EW') | |
| legend(loc=0) | |
| @plotfunc | |
| def emmdistbin(liness,scutsi,cutsd): | |
| from astropysics.utils import median_absolute_deviation,biweight_midvariance,interquartile_range | |
| d = cutsd['z500'].nd | |
| si = cutsd['z500'].ni | |
| li = liness.VAGCidsat | |
| sortedvagcinds = argsort(si) | |
| si = si[sortedvagcinds] | |
| d = d[sortedvagcinds] | |
| dli = d[np.searchsorted(si,li)] | |
| dbinedges = [0,20,40,60,100,150,200,250] | |
| linenames = ('Ha_ew','OII_ew','OIII_ew') | |
| colors = ('r','b','g') | |
| for linename,color in zip(linenames,colors): | |
| currliness = liness.field(linename) | |
| bincens = [] | |
| binhszs = [] | |
| binstds = [] | |
| binmeds = [] | |
| linesinbins = [] | |
| for l,u in zip(dbinedges[:-1],dbinedges[1:]): | |
| mask = (l <= dli)&(dli < u) &(currliness>-9999) | |
| bincens.append((u+l)/2) | |
| binhszs.append((u-l)/2) | |
| linesinbins.append(currliness[mask]) | |
| #binstds.append(biweight_midvariance(currliness[mask])[0]**0.5) | |
| binstds.append(std(currliness[mask])) | |
| #binstds.append(np.std(currliness[mask])) | |
| binmeds.append(np.median(currliness[mask])) | |
| bincens,binmeds,binstds,binhszs = np.array((bincens,binmeds,binstds,binhszs)) | |
| maxscale = max(abs(binmeds)) | |
| label = linename.replace('_ew','')+'/%.3f'%maxscale | |
| errorbar(bincens,binmeds/maxscale,binstds/maxscale,binhszs,fmt='o',color=color,label=label) | |
| xlim(0,250) | |
| ylim(-1.5,3) | |
| xlabel('$d_{\\rm proj}/(\\rm kpc/h)$') | |
| ylabel('scaled EW') | |
| legend(loc=0) | |
| @plotfunc | |
| def zdist(vagcd,cutsd,potsd): | |
| from astropysics.coords import cosmo_z_to_dist | |
| zs = vagcd['kcorr'].Z | |
| zh = zs[potsd['pothosts']] | |
| zdv = zs[cutsd['z500'].i] | |
| zedge = max(zh.max(),zdv.max()) | |
| dedge = cosmo_z_to_dist(zedge) | |
| hist(zdv,bins=25,histtype='step',normed=False,label=r'Clean Sample',lw=2) | |
| hist(zh,bins=25,histtype='step',normed=False,label=r'All Potential',lw=2,ls='dashed') | |
| xticks([0,.01,.02,.03]) | |
| xlim(0,zedge) | |
| xlabel(r'$z_{\rm pri}$') | |
| ylim(0,130) | |
| ylabel('$N$') | |
| legend(loc=2) | |
| twiny() | |
| xlim(0,dedge) | |
| xlabel(r'$d_{\rm pri}/({\rm comoving \, Mpc})$') | |
| @plotfunc(mainfig=True) | |
| def dvdists(vagcd,cutsd,msiicone): | |
| from scipy.stats import ks_2samp,kstest | |
| lw = 4 | |
| dvobs = 3e5*(vagcd['kcorr'].Z[cutsd['z500'].i]-vagcd['kcorr'].Z[cutsd['z500'].ni]) | |
| dvsim = msiicone.get_masked_arr('delta_v','conez,isolated,deltav') | |
| dvsimnoiso = msiicone.get_masked_arr('delta_v','conez,deltav') | |
| ks1 = ks_2samp(dvobs,dvsim) | |
| ks1o = kstest(dvobs,'norm') | |
| ks1s = kstest(dvsim,'norm') | |
| print '2-KS',ks1,'obs-norm',ks1o,'sim-norm',ks1s | |
| print 'obs sig',dvobs.std() | |
| print 'sim sig',dvsim.std() | |
| hist(dvobs,histtype='step',label='Clean Sample',ls='solid',lw=lw,normed=True,bins=25) | |
| hist(dvsim,histtype='step',label='Simulated Clean Sample',ls='dashed',lw=lw,normed=True,bins=25) | |
| legend(loc=0) | |
| xlabel(r'$\Delta v/({\rm km/s})$') | |
| xlim(-500,500) | |
| ylabel(r'$N_{\rm gal}$ (Linear normalized )') | |
| ylim(0,.004) | |
| ytcks = yticks() | |
| yticks(ytcks[0],['' for i in ytcks[1]]) | |
| @plotfunc | |
| def twodthreed(msiicone,h): | |
| subplots_adjust(left=.15,top=.97,right=.97) | |
| r2 = msiicone.get_masked_arr('r3dto2d','conez')*1000/h | |
| r3 = msiicone.get_masked_arr('r3d','conez')*1000/h | |
| r2iso = msiicone.get_masked_arr('r3dto2d','conez,isolated')*1000/h | |
| r3iso = msiicone.get_masked_arr('r3d','conez,isolated')*1000/h | |
| x = r2 | |
| y = (r3/r2) | |
| xs = np.linspace(0,250,101)[1:] | |
| ys = [percentile(y[x<xi],50) for xi in xs] | |
| yus = [percentile(y[x<xi],16) for xi in xs] | |
| yls = [percentile(y[x<xi],84) for xi in xs] | |
| plot(xs,ys,'o-b') | |
| fill_between(xs,yls,yus,color='b',alpha=.5) | |
| text(240,4.8,'Simulation',fontsize=24,va='top',ha='right') | |
| xlabel(r'$d_{\rm proj}$') | |
| xlim(0,250) | |
| ylabel(r'$d_{\rm 3D}/d_{\rm proj}$') | |
| ylim(1,5) | |
| @plotfunc | |
| def twodthreeddist(msiicone,h): | |
| subplots_adjust(left=.15,top=.97,right=.97) | |
| r2 = msiicone.get_masked_arr('r3dto2d','conez')*1000/h | |
| r3 = msiicone.get_masked_arr('r3d','conez')*1000/h | |
| r2iso = msiicone.get_masked_arr('r3dto2d','conez,isolated')*1000/h | |
| r3iso = msiicone.get_masked_arr('r3d','conez,isolated')*1000/h | |
| # plot(r2,r3,'.',ms=1,label='Satelluites of All Hosts') | |
| # plot(r2iso,r3iso,'.',ms=1,label='Satellites of Isolated Hosts') | |
| msk2 = r2<250 | |
| msk3 = r3<250 | |
| hist((r2/r3)[msk2],bins=50,histtype='step') | |
| axvline(median(r2[msk2]/r3[msk2]),ls='dashed',color='k') | |
| print 'Median r2/r3',median(r2/r3) | |
| print 'Median r2/r3 r2<250',median(r2[msk2]/r3[msk2]) | |
| print 'Median r2/r3 r3<250',median(r2[msk3]/r3[msk3]) | |
| xlabel(r'$d_{\rm 3D}/d_{\rm proj}$') | |
| xlim(0,1.03) | |
| ylabel(r'N') | |
| @plotfunc(mainfig=True) | |
| def lfuncs(cutsd,hiicutsd,msiicone,galmask,g,r,h,msiivtob03,b03tomsiiv,b05lf,potsd,vagcd): | |
| subplots_adjust(left=.10,bottom=.08,top=.97,right=.97) | |
| cutsd = hiicutsd #comment to get non HII-removed | |
| dcutmask = (cutsd['z500'].nd/h) < 250 | |
| z500hmask = np.zeros_like(galmask) | |
| z500hmask[cutsd['z500'].i[dcutmask]] = True | |
| z500hmask = z500hmask[galmask] | |
| z500smask = np.zeros_like(galmask) | |
| z500smask[cutsd['z500'].ni[dcutmask]] = True | |
| z500smask = z500smask[galmask] | |
| prim = r[z500hmask] #primary mags | |
| secm = r[z500smask] #sec mags | |
| #prim = vagcd['kcorr'].ABSMAG[potsd['pothosts'],2]+5*log10(h) | |
| vtom,mtov = msiivtob03,b03tomsiiv | |
| pbins = np.linspace(np.min(prim),-20+5*log10(h),16) | |
| pcens = np.convolve(pbins,[.5,.5],mode='valid') | |
| sbins = np.linspace(-20+5*log10(h),-17.5+5*log10(h),16) | |
| scens = np.convolve(sbins,[.5,.5],mode='valid') | |
| hhist = [] | |
| shist = [] | |
| for realization in range(len(msiicone.arrs)): | |
| hsrproj = msiicone.get_masked_arr('rproj','conez,isolated,deltav',realization) | |
| hostv = msiicone.get_masked_arr('h_vmax','conez,isolated,deltav',realization) | |
| satv = msiicone.get_masked_arr('s_vmax','conez,isolated,deltav',realization) | |
| hostm = vtom(hostv[(hsrproj*1e3/h)<250]) + 5*log10(h) | |
| hhist.append(histogram(hostm,bins=pbins,normed=True)[0]) | |
| satm = vtom(satv[(hsrproj*1e3/h)<250]) + 5*log10(h) | |
| shist.append(histogram(satm,bins=sbins,normed=True)[0]) | |
| hhist,shist = array(hhist),array(shist) | |
| #pri/host | |
| npri = histogram(prim,bins=pbins)[0] | |
| errorbar(pcens,npri,npri**0.5,fmt='o',mfc='m',ecolor='k',label='Clean (Observed) Pris') | |
| hhist *= np.mean(npri)/np.mean(hhist) | |
| hmed = percentile(hhist,50,axis=0) | |
| herr = array([percentile(hhist,15.85,axis=0),percentile(hhist,84.15,axis=0)]) | |
| hvar = var(hhist,axis=0) | |
| bar(pbins[:-1],herr[1]-herr[0],np.diff(pbins[:-1]).mean(),herr[0], | |
| color=(1,.5,.5),linewidth=0,label='Simulated Clean Pris') | |
| rc2h = ((hmed-npri)**2/(hvar+npri)).sum()/(pbins.size-1) | |
| text(-22.5,55,r'$\chi^2_{\rm red}=%.1f$'%rc2h,ha='right',va='center',fontsize=18,color='r') | |
| print 'Rchi2 w/o faintest pri bin:',(((hmed-npri)**2/(hvar+npri))[:-1]).sum()/(pbins.size-1) | |
| plfx = linspace(pbins[0],pbins[-1],50) | |
| plfy = interp(plfx,b05lf['m'][::-1] + 5*log10(h),b05lf['p1'][::-1]) | |
| plfy *= np.mean(npri)/np.mean(plfy) | |
| plot(plfx,plfy,'--k') | |
| #sec/sat | |
| nsec = histogram(secm,bins=sbins)[0] | |
| errorbar(scens,nsec,nsec**0.5,fmt='o',mfc='b',ecolor='k',label='Clean (Observed) Secs') | |
| shist *= np.mean(nsec)/np.mean(shist) | |
| smed = percentile(shist,50,axis=0) | |
| serr = array([percentile(shist,15.85,axis=0),percentile(shist,84.15,axis=0)]) | |
| svar = var(shist,axis=0) | |
| bar(sbins[:-1],serr[1]-serr[0],np.diff(sbins[:-1]).mean(),serr[0], | |
| color=(.5,.75,.5),linewidth=0,label='Simulated Clean Secs') | |
| rc2s = ((smed-nsec)**2/(svar+nsec)).sum()/(sbins.size-1) | |
| text(-19.75,45,r'$\chi^2_{\rm red}=%.1f$'%rc2s,ha='left',va='center',fontsize=18,color='g') | |
| slfx = linspace(sbins[0],sbins[-1],50) | |
| slfy = interp(slfx,b05lf['m'][::-1] + 5*log10(h),b05lf['p1'][::-1]) | |
| slfy *= np.mean(nsec)/np.mean(slfy) | |
| plot(slfx,slfy,'--k',label='Blanton 2005') | |
| axvline(-20+5*log10(h),c='k') | |
| 'Clean Sample, Selec. Corr. (Observed)' 'Simulated Clean Sample' | |
| xlim(sbins[-1],pbins[0]) | |
| xlabel('$M_r$') | |
| ylabel(r'$N_{\rm gal,clean}$') | |
| ylim(0,70) | |
| secmid = (-20-17.5)/2+5*log10(h) | |
| primid = (-20+pbins[0]+5*log10(h))/2 | |
| text(secmid,66,'Secondaries',fontsize=20,va='top',ha='center') | |
| text(primid,66,'Primaries',fontsize=20,va='top',ha='center') | |
| # h,l = gca().get_legend_handles_labels() | |
| # #blanton to end | |
| # h.insert(4,h.pop(2)) | |
| # l.insert(4,l.pop(2)) | |
| # #pris and secs together | |
| # h.insert(2,h.pop(1)) | |
| # l.insert(2,l.pop(1)) | |
| # legend(h,l,loc=2) | |
| @plotfunc | |
| def lfuncsks(cutsd,hiicutsd,msiicone,galmask,g,r,h,msiivtob03,b03tomsiiv,msiivtob05,b05tomsiiv): | |
| from scipy.stats import ks_2samp | |
| cutsd = hiicutsd #comment to get non HII-removed | |
| vtom,mtov = msiivtob03,b03tomsiiv | |
| realization = 5#slice(None) | |
| hsrproj = msiicone.get_masked_arr('rproj','conez,isolated,deltav',realization) | |
| hostv = msiicone.get_masked_arr('h_vmax','conez,isolated,deltav',realization) | |
| hostm = vtom(hostv[(hsrproj*1e3/h)<250]) + 5*log10(h) | |
| satv = msiicone.get_masked_arr('s_vmax','conez,isolated,deltav',realization) | |
| satm = vtom(satv[(hsrproj*1e3/h)<250]) + 5*log10(h) | |
| dcutmask = (cutsd['z500'].nd/h) < 250 | |
| z500hmask = np.zeros_like(galmask) | |
| z500hmask[cutsd['z500'].i[dcutmask]] = True | |
| z500hmask = z500hmask[galmask] | |
| z500smask = np.zeros_like(galmask) | |
| z500smask[cutsd['z500'].ni[dcutmask]] = True | |
| z500smask = z500smask[galmask] | |
| prim = r[z500hmask] | |
| secm = r[z500smask] | |
| print 'nhost/npri',hostm.size/prim.size,hostm.size,prim.size | |
| hist(prim,bins=15,histtype='step',ls='solid',color='r',normed=True) | |
| hist(hostm,bins=15,histtype='step',ls='dotted',color='m',normed=True) | |
| hsks = ks_2samp(prim,hostm) | |
| print 'host/pri KS',hsks | |
| brightest = min(prim.min(),hostm.min()) | |
| text(-22,.9,r'$p= %0.1f \%%$'%(hsks[1]*100),ha='left',va='center',fontsize=18,color='r') | |
| print 'nsat/nsec',satm.size/secm.size,satm.size,secm.size | |
| hist(secm,bins=15,histtype='step',ls='solid',color='g',normed=True) | |
| hist(satm,bins=15,histtype='step',ls='dotted',color='b',normed=True) | |
| ssks = ks_2samp(secm,satm) | |
| print 'sec/sat KS',ssks | |
| faintest = max(satm.max(),secm.max()) | |
| text(-18.5,.75,r'$p= %0.1f \%%$'%(ssks[1]*100),ha='left',va='center',fontsize=18,color='g') | |
| xlim(faintest,brightest) | |
| xlabel('$M_r$') | |
| yticks([]) | |
| ylabel('Abundance (Linear normalized)') | |
| def _do_abel(obsx,obsy,nbins): | |
| obssig = obsy/obsx #surfden w/o normalization | |
| dsig = np.diff(obssig) | |
| xdsig = np.convolve(obsx,[.5,.5],mode='valid') | |
| r = np.linspace(0,xdsig.max(),nbins) | |
| rs = np.tile(r,dsig.size).reshape(dsig.size,r.size).T | |
| xs = np.tile(xdsig,r.size).reshape(r.size,xdsig.size) | |
| ig = (-1/pi)*np.tile(dsig,r.size).reshape(r.size,dsig.size)*(xs**2-rs**2)**-0.5 | |
| ig[xs<rs]=0 | |
| rho = np.sum(ig,axis=1) | |
| return r,rho | |
| @plotfunc(mainfig=True) | |
| def threedprof(cutsd,hiicutsd,nhiicutsd,potsd,msiicone,h): | |
| from astropysics.models import LinearModel,BurkertModel,TwoPowerModel,PowerLawModel | |
| from astropysics.plotting import cumulative_plot | |
| cutsd = hiicutsd #comment this to go back to un-visually filtered for HII regions | |
| subplots_adjust(left=.17,right=.965,top=.95) | |
| redge = 250 #in kpc/h | |
| nbins1 = 42 | |
| nbins2 = 20 | |
| frac250 = 0.425594195247 | |
| binedges = linspace(0,redge,nbins1+1) | |
| bincens = correlate(binedges,[.5,.5],mode='valid') | |
| dbin = bincens[1]-bincens[0] | |
| nds = cutsd['z500'].nd | |
| ndsspec = cutsd['nozcut'].nd | |
| nds1 = cutsd['z500'].nd[cutsd['z500'].nn == 1] | |
| npots = len(potsd['pothosts']) | |
| ndsuniquenoz = cutsd['uniquenoz'].nd | |
| ndsphot = np.concatenate((ndsuniquenoz,ndsspec)) | |
| #sel corr clean | |
| selfuncphot = histogram(ndsphot/h,binedges)[0] | |
| selfuncz500 = histogram(nds/h,binedges)[0] | |
| selfuncspec = histogram(ndsspec/h,binedges)[0] | |
| selfunccorr = selfuncphot*(selfuncz500/selfuncspec) | |
| f = selfunccorr/npots | |
| obsx = bincens/h | |
| obsy = f | |
| obsyerr = f*(1/selfunccorr+1/npots)**0.5 | |
| #abel inversion | |
| r, rho = _do_abel(obsx,obsy,nbins2) | |
| lr,lrho = log10(r),log10(rho) | |
| #make linear fit to be able to do offset | |
| mfin = np.isfinite(lr)&np.isfinite(lrho) | |
| lr = lr[mfin] | |
| lrho = lrho[mfin] | |
| m = LinearModel() | |
| print 'Lin Fit',m.fitBasic(lr,lrho) | |
| m.fitData(lr,lrho) | |
| #compute the power law over the range shown (assume 0 inside, basically) | |
| minlr = 10**lr.min() | |
| pm = PowerLawModel(A=1,p=m.m) | |
| pm.A = frac250/pm.integrateSpherical(minlr,250) | |
| print pm.integrateSpherical(minlr,250),log10(pm.A) | |
| lrhooffset = m.b - log10(pm.A) | |
| m.b = log10(pm.A) | |
| lrho -= lrhooffset | |
| #m.fitData(lr,lrho) | |
| scatter(lr,lrho,s=30,c='b',label='Clean Sample (Abel Inverted)',zorder=10) | |
| m.plot(1.2,2.5,ls='--',c='k',data=None,clf=False,label='Power Law fit to Observed') | |
| print 'Power law',m.pardict,10**m.b | |
| #now do random errors in y | |
| scat = [] | |
| for rv in randn(100): | |
| newy = obsy+obsyerr*rv | |
| scat.append(_do_abel(obsx,newy,nbins2)) | |
| rhos = np.log10(scat)[:,1] | |
| rhos = np.ma.masked_where(~np.isfinite(rhos),rhos) | |
| stds = rhos.std(axis=0) | |
| #try bootstrap | |
| bootstrap = [] | |
| for inds in randint(obsy.size,size=(100,obsy.size)): | |
| newy = obsy[inds] | |
| bootstrap.append(_do_abel(obsx,newy,nbins2)) | |
| rhos = np.log10(bootstrap)[:,1] | |
| rhos = np.ma.masked_where(~np.isfinite(rhos),rhos) | |
| stds = rhos.std(axis=0) | |
| #both | |
| botherrs = [] | |
| rvs = randn(100) | |
| bsinds = randint(obsy.size,size=(100,obsy.size)) | |
| for rv,inds in zip(rvs,bsinds): | |
| newy = obsy[inds] + obsyerr[inds]*rv | |
| botherrs.append(_do_abel(obsx,newy,nbins2)) | |
| rhos = np.log10(botherrs)[:,1] | |
| rhos = np.ma.masked_where(~np.isfinite(rhos),rhos) | |
| stds = rhos.std(axis=0) | |
| errorbar(lr,lrho,stds.data[mfin],fmt=None,ecolor='k') | |
| #now add MS-II | |
| simclean3d = msiicone.get_masked_arr('r2dto3d','isolated,deltav,conez')*1e3/h | |
| simb = logspace(lr.min(),3,30) | |
| nbs = np.histogram(simclean3d,bins=simb)[0] | |
| simbcen = np.convolve(simb,(.5,.5),mode='valid') | |
| simrho = nbs/((4*pi/3)*(simb[1:]**3-simb[:-1]**3)) | |
| lsimx,lsimy = log10(simbcen),log10(simrho) | |
| lxmid = np.mean(lsimx) | |
| simoffset = m(lxmid) - interp(lxmid,lsimx,lsimy) | |
| plot(lsimx,lsimy+simoffset,'r-',label='Simulated Clean Sample') | |
| bm1 = BurkertModel(rho0=1,r0=67) | |
| bm1.setCall(xtrans='pow',ytrans='log') | |
| bm1.fitData(lr,lrho) | |
| #now integrate it to match the fraction @ 250 | |
| # bm1.setCall(xtrans=None,ytrans=None) | |
| # bm1.rho0 /= frac250/bm1.integrateSpherical(minlr,250) | |
| bm1.setCall(xtrans='pow',ytrans='log') | |
| bm1.plot(1.25,2.5,data=None,clf=False,label='Burkert Model of Observed',ls=':',color='g') | |
| print 'Burkert core',bm1.r0,'log',bm1.rho0 | |
| xlabel(r'$\log(r/{\rm kpc})$') | |
| ylabel(r'$\log(n/[{\rm kpc}^{-3}])$') | |
| xlim(1.25,2.5) | |
| ylim(-10,-6) | |
| h,l = gca().get_legend_handles_labels() | |
| #move clean to top | |
| h.insert(0,h.pop(-1)) | |
| l.insert(0,l.pop(-1)) | |
| #sim clean to bottom | |
| h.append(h.pop(2)) | |
| l.append(l.pop(2)) | |
| legend(h,l,loc=0) | |
| #return lr,lrho | |
| @plotfunc | |
| def threedprof2(cutsd,hiicutsd,nhiicutsd,potsd,msiicone,h): | |
| from astropysics.phot import distance_from_modulus | |
| from astropysics.coords import cosmo_z_to_dist | |
| cutsd = hiicutsd #comment this to go back to un-visually filtered for HII regions | |
| subplots_adjust(left=.15,right=.95,top=.95) | |
| nds = cutsd['z500'].nd | |
| ndsspec = cutsd['nozcut'].nd | |
| nds1 = cutsd['z500'].nd[cutsd['z500'].nn == 1] | |
| npots = len(potsd['pothosts']) | |
| ndsuniquenoz = cutsd['uniquenoz'].nd | |
| ndsphot = np.concatenate((ndsuniquenoz,ndsspec)) | |
| selfuncx = sort(ndsphot)#kpc/h | |
| selfuncphotf = (arange(ndsphot.size)+1)/npots | |
| selfuncz500f = interp(selfuncx,sort(nds),(arange(nds.size)+1)/npots) | |
| selfuncspecf = interp(selfuncx,sort(ndsspec),(arange(ndsspec.size)+1)/npots) | |
| selfunccorrf = selfuncphotf*(selfuncz500f/selfuncspecf) | |
| err = selfunccorrf*((selfuncphotf*npots)**-1+(selfuncz500f*npots)**-1+(selfuncspecf*npots)**-1)**0.5 | |
| x = selfuncx/h | |
| y = selfunccorrf | |
| xd = np.linspace(0,250,22) | |
| sig = np.diff(interp(xd,x,y))/(pi*(xd[1:]**2-xd[:-1]**2)) | |
| plot(xd[:-1],sig) | |
| dsig = np.diff(sig) | |
| xd[1:-1] | |
| plot(xd[1:-1],dsig) | |
| #legend(loc=4) | |
| ylabel(r'$f_{\rm sat}(<d_{\rm proj})$') | |
| xlabel(r'$d_{\rm proj}/{\rm kpc}$') | |
| #ylim(0,1) | |
| xlim(0,250) | |
| def cc_plots(data,varnames,vartoplot,doclf=True): | |
| cc = corrcoef(data) | |
| if isinstance(vartoplot,basestring): | |
| i = varnames.index(vartoplot) | |
| varstr = vartoplot | |
| else: | |
| i = vartoplot | |
| varstr = varnames[i] | |
| cce=list(cc[i]) | |
| cce.insert(0,cce[i]) | |
| cce=array(cce) | |
| if doclf: | |
| clf() | |
| subplots_adjust(bottom=.28,top=.93) | |
| plot(arange(cce.size),cce,ls='steps') | |
| hlines(0,0,len(varnames),colors='k',linestyles='dashed') | |
| xticks(arange(len(varnames))+.5,varnames,rotation=90) | |
| xlim(0,len(varnames)) | |
| ylim(-1,1) | |
| title(varstr) | |
| def visual_check(vagcd,cutsd,cuttype='z500',fnbase='vischeck',dr=7,randomorder=False,secctype=None): | |
| """ | |
| run this to visually check the | |
| """ | |
| import os,time,random | |
| fn = fnbase+'.'+'dr%i'%dr+'.'+cuttype | |
| hostsvisited = [] | |
| if os.path.exists(fn): | |
| with open(fn,'r') as f: | |
| for l in f: | |
| if l.strip() != '': | |
| hosti,sati,val = l.split()[:3] | |
| hostsvisited.append(int(hosti)) | |
| if len(hostsvisited)>0: | |
| print 'Already visited',len(hostsvisited),'Pairs' | |
| with open(fn,'a') as f: | |
| quit = False | |
| cutzip = list(enumerate(zip(cutsd[cuttype].i,cutsd[cuttype].ni,cutsd[cuttype].nd))) | |
| if randomorder: | |
| random.shuffle(cutzip) | |
| ivis = 0 | |
| for i,(hi,ni,nd) in cutzip: | |
| if hi in hostsvisited: | |
| continue | |
| zh = vagcd['kcorr'].Z[hi] | |
| zs = vagcd['kcorr'].Z[ni] | |
| rh = vagcd['kcorr'].ABSMAG[hi,2] | |
| rs = vagcd['kcorr'].ABSMAG[ni,2] | |
| gh = vagcd['kcorr'].ABSMAG[hi,1] | |
| gs = vagcd['kcorr'].ABSMAG[ni,1] | |
| # get_skyserver_url(hi,vagcd,mode='navigate',dr=dr,openpage=True) | |
| # time.sleep(.5) | |
| # get_skyserver_url(ni,vagcd,mode='navigate',dr=dr,openpage=True) | |
| # time.sleep(1) | |
| get_skyserver_url(i,vagcd,cuts=cutsd[cuttype],mode='cutout',dr=dr,openpage=True) | |
| time.sleep(1) | |
| get_skyserver_url(ni,vagcd,mode='navigateST',dr=dr,openpage=True) | |
| time.sleep(.5) | |
| print ivis+1,'of',cutsd[cuttype].size-len(hostsvisited),'Sep='+str(nd),'r=%.2f,%.2f'%(rh,rs),'g-r=%.2f,%.2f'%(gh-rh,gs-rs),'z=%.5f,%.5f'%(zh,zs) | |
| if secctype: | |
| m = where(cutsd[secctype].i==hi)[0] | |
| if len(m)==0: | |
| print 'no',secctype,'matches' | |
| else: | |
| for mi in m: | |
| match = cutsd[secctype][mi] | |
| secra = vagcd['kcorr'].RA[match.ni] | |
| secdec = vagcd['kcorr'].DEC[match.ni] | |
| print secctype,'matches host',match.i,'sep:',match.nd,'RA',secra,'Dec',secdec | |
| get_skyserver_url(match.ni,vagcd,mode='navigateST',dr=dr,openpage=True) | |
| time.sleep(.5) | |
| needinput = True | |
| while needinput: | |
| resri = raw_input('Good=1,Warning=2,Unk=0,Bad sat=-1,Bad host=-2,Possble Grp=-3,quit=q:') | |
| needinput = False | |
| resri = resri.split(',') | |
| comment = resri[1] if len(resri)>1 else '' | |
| resri = resri[0] | |
| if resri.strip() == '1': | |
| f.write('%i\t%i\t1\tGood %s\n'%(hi,ni,comment)) | |
| elif resri.strip() == '2': | |
| f.write('%i\t%i\t2\tWarning %s\n'%(hi,ni,comment)) | |
| elif resri.strip() == '0': | |
| f.write('%i\t%i\t0\tUnknown %s\n'%(hi,ni,comment)) | |
| elif resri.strip() == '-1': | |
| f.write('%i\t%i\t-1\tBadsat %s\n'%(hi,ni,comment)) | |
| elif resri.strip() == '-2': | |
| f.write('%i\t%i\t-2\tBadhost %s\n'%(hi,ni,comment)) | |
| elif resri.strip() == '-3': | |
| f.write('%i\t%i\t-3\tPossiblegrp %s\n'%(hi,ni,comment)) | |
| elif resri.strip().lower() == 'q': | |
| quit = True | |
| else: | |
| needinput = True | |
| hostsvisited.append(hi) | |
| ivis +=1 | |
| if quit: | |
| break | |
| def build_vischeck_table(vagcd,cutsd,cuttype='z500',fn=None,order=None,mask=None): | |
| """ | |
| vagcd and cutsd are from the loaded files, cuttype is a key of cutsd, fn is | |
| the filename, `order` is 'random','increasing',or 'decreasing', and mask must | |
| match the cutsd[cuttype]'s shape or be None. | |
| """ | |
| import random | |
| if fn is None: | |
| orderstr = '' if order is None else ('_'+str(order)) | |
| fn = 'vischeck_%s_%s%s'%(cuttype,vagcd['dr'],orderstr) | |
| vagc = vagcd['kcorr'] | |
| cuts = cutsd[cuttype] | |
| if mask is not None: | |
| cuts = cuts[mask] | |
| print 'Making table at',fn | |
| with open(fn,'w')as f: | |
| f.write('#cut:%s,vagcdr:%s,count:%i'%(cuttype,vagcd['dr'],len(cuts))) | |
| if mask is not None: | |
| f.write(',mask present') | |
| f.write('\n') | |
| f.write('#cuti pri_vagci sec_vagci pri_ra pri_dec sec_ra sec_dec pri_z sec_z sep pri_r sec_r pri_g-r sec_g-r') | |
| f.write('\n') | |
| cutzip = list(enumerate(zip(cuts.i,cuts.ni,cuts.nd))) | |
| if order=='random': | |
| random.shuffle(cutzip) | |
| elif order=='distance' or order=='increasing' or order=='decreasing': | |
| cutzip.sort(cmp=lambda x,y:cmp(x[1][2],y[1][2]),reverse=order=='decreasing') | |
| elif order is not None: | |
| raise ValueError('invalid order '+str(order)) | |
| for i,(pi,si,d) in cutzip: | |
| rap,decp = vagc.RA[pi],vagc.DEC[pi] | |
| ras,decs = vagc.RA[si],vagc.DEC[si] | |
| zp,zs = vagc.Z[pi],vagc.Z[si] | |
| gp,gs = vagc.ABSMAG[(pi,si),1] | |
| rp,rs = vagc.ABSMAG[(pi,si),2] | |
| gmrp = gp-rp | |
| gmrs = gs-rs | |
| line = '{i} {pi} {si} {rap} {decp} {ras} {decs} {zp} {zs} {d} {rp} {rs} {gmrp} {gmrs}' | |
| f.write(line.format(**locals())) | |
| f.write('\n') | |
| def do_vischeck_table(fn,opts='',offset=0): | |
| import webbrowser,time | |
| from shutil import move | |
| if opts!='': | |
| opts = '&opt='+opts | |
| wlines = [] | |
| with open(fn,'r+') as f: | |
| l1 = f.readline()[:-1] | |
| firstlineelems = l1[1:].split(',') | |
| if len(firstlineelems)>3: | |
| print 'Mask used to create this vischeck table' | |
| cut,vagc,cnt = firstlineelems[:3] | |
| cut = cut.replace('cut:','') | |
| dr = int(vagc.replace('vagcdr:','').replace('dr','')) | |
| cnt = cnt.replace('count:','') | |
| l2 = f.readline()[:-1] | |
| wlines.append(l1) | |
| if 'qcode comment' in l2: | |
| wlines.append(l2) | |
| else: | |
| wlines.append(l2+' qcode comment') | |
| quit = False | |
| locali = 0 | |
| ioff = 0 | |
| prompt = 'Good=1,Warning=2,Unknown=0,Bad sat=-1,Bad host=-2,Possble Grp=-3,repeat=r,skip=s,quit=q:' | |
| lastell = f.tell() | |
| for l in f: | |
| locali += 1 | |
| if not quit: | |
| try: | |
| ls = l[:-1].split() | |
| if len(ls)!=14 : | |
| valid = 'start' | |
| elif ioff<offset: | |
| valid = 'offset' | |
| ioff+=1 | |
| else: | |
| valid = False | |
| while not valid: | |
| i,pi,si,rap,decp,ras,decs,zp,zs,sep,rp,rs,gmrp,gmrs = ls | |
| print '{locali} of {cnt},z={zp},r_pri={rp},r_sec={rs},sep={sep}'.format(**locals()) | |
| webbrowser.open_new_tab(_construct_url(ras,decs,None,None,float(zs if float(zs)>=-1 else zp),dr,'navigate')+opts) | |
| time.sleep(.5) | |
| webbrowser.open_new_tab(_construct_url(rap,decp,ras,decs,float(zp),dr,'cutout')+opts) | |
| time.sleep(1) | |
| res = raw_input(prompt) | |
| sres = res.split(',') | |
| if len(sres)==1: | |
| q = sres[0] | |
| c = None | |
| elif len(sres)==2: | |
| q = int(sres[0]) | |
| c = sres[1] | |
| else: | |
| print 'Invalid number of inputs, repeating!' | |
| continue | |
| if q=='r': | |
| valid = False | |
| elif q=='s': | |
| valid = 'skip' | |
| elif q=='q': | |
| quit = True | |
| valid = 'quit' | |
| else: | |
| try: | |
| q = int(q) | |
| valid = True | |
| except ValueError: | |
| print 'Invalid input, repeating!' | |
| if valid is True: | |
| l2 = l[:-1]+' '+str(q) | |
| if c is not None: | |
| l2 += ' '+c | |
| wlines.append(l2) | |
| else: | |
| wlines.append(l[:-1]) | |
| except Exception,e: | |
| wlines.append(l[:-1]) | |
| quit = e | |
| else: | |
| wlines.append(l[:-1]) | |
| lastell = f.tell() | |
| with open(fn,'w') as f: | |
| for l in wlines: | |
| f.write(l) | |
| f.write('\n') | |
| if quit and quit is not True: | |
| raise quit | |
| def get_vischeck_filtered(fn,filter=None,includenonvis=True,matchcuts=None): | |
| """ | |
| Filter can be a callable f(code,comment) returning a bool (if True, the | |
| object will be included), a '<#' string, or a number (will filter everything | |
| less than that). | |
| nmatchd specifies a cuts array that can be used to get 'nn' | |
| returns filteredcuts,ncut | |
| """ | |
| from operator import isSequenceType | |
| dt = {'i':int,'ni':int,'nd':float}.items() | |
| if matchcuts is not None: | |
| dt.append(('nn',int)) | |
| dt = dtype(dt) | |
| pris = [] | |
| secis = [] | |
| seps = [] | |
| if filter is None: | |
| filter = lambda code,comment:True | |
| elif callable(filter): | |
| pass | |
| elif isinstance(filter,basestring): | |
| filter = filter.strip() | |
| if filter.startswith('<'): | |
| filternum = float(filter.replace('<','')) | |
| filter = lambda code,comment:float(code)<filternum | |
| elif filter.startswith('>'): | |
| filternum = float(filter.replace('>','')) | |
| filter = lambda code,comment:float(code)>filternum | |
| elif filter.startswith('<='): | |
| filternum = float(filter.replace('<=','')) | |
| filter = lambda code,comment:float(code)<=filternum | |
| elif filter.startswith('>='): | |
| filternum = float(filter.replace('>=','')) | |
| filter = lambda code,comment:float(code)>=filternum | |
| elif filter.startswith('=='): | |
| filternum = float(filter.replace('==','')) | |
| filter = lambda code,comment:float(code)==filternum | |
| else: | |
| filternum = float(filter) | |
| filter = lambda code,comment:float(code)>filternum | |
| elif isSequenceType(filter): | |
| filterwords = filter | |
| filter = lambda code,comment:all([s not in comment for s in filterwords]) | |
| else: | |
| filternum = filter | |
| filter = lambda code,comment:float(code)>filternum | |
| ncut = 0 | |
| with open(fn) as f: | |
| for l in f: | |
| ls = l.strip() | |
| if ls and not ls.startswith('#'): | |
| lss = ls.split() | |
| pri = int(lss[1]) | |
| seci = int(lss[2]) | |
| sep = float(lss[9]) | |
| code = int (lss[14]) if len(lss)>14 else None | |
| comment = ' '.join(lss[15:]) if len(lss)>15 else '' | |
| if code is None: | |
| add = includenonvis | |
| else: | |
| add = filter(code,comment) | |
| if add: | |
| pris.append(pri) | |
| secis.append(seci) | |
| seps.append(sep) | |
| else: | |
| ncut +=1 | |
| arr = np.empty((len(pris),),dtype=dt) | |
| arr['i'] = pris | |
| arr['ni'] = secis | |
| arr['nd'] = seps | |
| if matchcuts is not None: | |
| sorti = argsort(matchcuts['i']) | |
| matches = searchsorted(matchcuts['i'][sorti],arr['i']) | |
| arr['nn'] = matchcuts['nn'][matches] | |
| return arr.view(np.recarray),ncut | |
| def _load_vischeck_cutsds(cutsd,filter,d2fn): | |
| filtcuts = {} | |
| ncutd = {} | |
| for k,c in cutsd.iteritems(): | |
| if k in d2fn: | |
| cuts,ncut = get_vischeck_filtered(d2fn[k],filter,matchcuts=c) | |
| filtcuts[k] = cuts | |
| ncutd[k] = ncut | |
| else: | |
| filtcuts[k] = c | |
| ncutd[k] = 0 | |
| return filtcuts,ncutd | |
| def _copy_vischecks(fn1='test',fn2='vischeck_z500_dr7_increasing',wfn=None): | |
| """ | |
| use `fn1` as base, add in `fn2` codes/comments, write to `fn` | |
| """ | |
| strs = [] | |
| with open(fn1) as f1: | |
| for l in f1: | |
| if not l.startswith('#') and not l.strip()=='': | |
| i = int (l.split()[0]) | |
| with open(fn2) as f2: | |
| for l2 in f2: | |
| l2s = l2.split() | |
| if len(l2s)>0 and not l2s[0].startswith('#') and int(l2s[0])==i: | |
| strs.append(l2) | |
| break | |
| else: | |
| strs.append(l) | |
| else: | |
| strs.append(l) | |
| if wfn: | |
| with open(wfn,'w') as fw: | |
| fw.writelines(strs) | |
| return strs | |
| def abund_match_samples(n1,n2,gtr1=True,gtr2=True,doplot=False): | |
| """ | |
| Does abundance matching, and returns functions 1->2 and 2->1. | |
| `gtr1` and `gtr2` specify if the matching should be for n>val (True) | |
| or n<val (False) for each of the abundances. | |
| use as abund_match(mag,vmax,False,True) | |
| """ | |
| from scipy.interpolate import interp1d | |
| n1 = np.sort(n1) | |
| if gtr1: | |
| n1 = n1[::-1] | |
| n2 = np.sort(n2) | |
| if gtr2: | |
| n2 = n2[::-1] | |
| p1 = linspace(0,1,n1.size) | |
| p2 = linspace(0,1,n2.size) | |
| if doplot: | |
| if doplot is True: | |
| doplot = 'cdf' | |
| if not ('noclf' in doplot): | |
| clf() | |
| if 'cdf' in doplot: | |
| l1 = plot(n1,p1,'-b',label='PDF 1')[0] | |
| xlabel('PDF 1') | |
| ylabel('p') | |
| twiny() | |
| l2 = plot(n2,p2,'--g',label='PDF 2')[0] | |
| xlabel('PDF 2') | |
| legend([l1,l2],['PDF1','PDF2']) | |
| elif 'pdf' in doplot: | |
| doplot.replace('pdf','').replace('noclf','') | |
| l1 = hist(n1,color='b',linestyle='solid',lw=2,histtype='step')[-1] | |
| xlabel('PDF 1') | |
| ylabel('N') | |
| twiny() | |
| l2 = hist(n2,color='g',linestyle='dashed',lw=2,histtype='step')[-1] | |
| xlabel('PDF 2') | |
| legend([l1,l2],['PDF1','PDF2']) | |
| else: | |
| raise ValueError('unrecognized doplot') | |
| subplots_adjust(top=.9) | |
| n2grid1 = np.interp(p1,p2,n2) #if gtr2 else np.interp(p1[::-1],p2[::-1],n2[::-1]) | |
| i12 = interp1d(n1[::-1],n2grid1[::-1]) if gtr1 else interp1d(n1,n2grid1) | |
| n1grid2 = np.interp(p2,p1,n1) #if gtr1 else np.interp(p2[::-1],p1[::-1],n1[::-1]) | |
| i21 = interp1d(n2[::-1],n1grid2[::-1]) if gtr2 else interp1d(n2,n1grid2) | |
| return i12,i21 | |
| def abund_match_sim_to_sch(simvals,simvol,phistar,alpha,Mstar,nsamples=100,doplot=False,h4mag=1): | |
| """ | |
| Abundance matches a sampled simulation to a Schechter function (luminosity function) | |
| retruns (sim->lf,lf->sim) callables | |
| """ | |
| from astropysics.models import SchechterMagModel | |
| from scipy.interpolate import interp1d | |
| mod = SchechterMagModel(phistar=phistar,alpha=alpha,Mstar=Mstar) | |
| mod.setCall('integrate',lower=Mstar-20) #20 is fairly arbitrary but more than enough | |
| simden = (np.arange(simvals.size)+1)/simvol | |
| simvals = np.sort(simvals)[::-1] | |
| dens = logspace(log10(simden[0]),log10(simden[-1]),nsamples) | |
| simvalsamps = interp(dens,simden,simvals) | |
| lfsamps = array([mod.inv(d) for d in dens]) + 5*log10(h4mag) | |
| if doplot: | |
| clf() | |
| l1 = semilogy(lfsamps,dens,'-b',label='LF')[0] | |
| xlabel('mags') | |
| xlim(*xlim()[::-1]) | |
| ylabel('$n$') | |
| twiny() | |
| l2 = semilogy(simvalsamps,dens,'--g',label='Simulation')[0] | |
| xlabel('sim units') | |
| legend((l1,l2),('LF','Sim')) | |
| subplots_adjust(top=.9) | |
| sorti1 = argsort(simvalsamps) | |
| sorti2 = argsort(lfsamps) | |
| return (interp1d(simvalsamps[sorti1],lfsamps[sorti1]), | |
| interp1d(lfsamps[sorti2],simvalsamps[sorti2])) | |
| def abund_match_sim_to_data(simvals,simvol,lfmags,lfphi,doplot=False,h4mag=1): | |
| """ | |
| Abundance matches a sampled simulation to a sampled LF | |
| returns (sim->lf,lf->sim) callables | |
| """ | |
| from scipy.interpolate import interp1d | |
| from scipy.integrate import cumtrapz | |
| sortlf = argsort(lfmags) | |
| lfmags,lfphi = lfmags[sortlf],lfphi[sortlf] | |
| lfden = cumtrapz(lfphi,lfmags) | |
| lfden = np.insert(lfden.astype(float),0,0) | |
| simden = (np.arange(simvals.size)+1)/simvol | |
| simvals = np.sort(simvals)[::-1] | |
| if simvals.size < lfmags.size: | |
| dens = simden | |
| simvalsamps = simvals | |
| #need to resample lf | |
| sorti = argsort(lfden) | |
| lfsamps = interp(dens,lfden[sorti],lfmags[sorti]) + 5*log10(h4mag) | |
| else: | |
| dens = lfden | |
| lfsamps = lfmags + 5*log10(h4mag) | |
| #need to resample sim | |
| sorti = argsort(simden) | |
| simvalsamps = interp(dens,simden[sorti],simvals[sorti]) | |
| if doplot: | |
| clf() | |
| l1 = semilogy(lfsamps,dens,'-b',label='LF')[0] | |
| xlabel('mags') | |
| xlim(*xlim()[::-1]) | |
| ylabel('$n$') | |
| twiny() | |
| if 'logsim'==doplot: | |
| l2 = loglog(simvalsamps,dens,'--g',label='Simulation')[0] | |
| else: | |
| l2 = semilogy(simvalsamps,dens,'--g',label='Simulation')[0] | |
| xlabel('sim units') | |
| legend((l1,l2),('LF','Sim')) | |
| subplots_adjust(top=.9) | |
| sorti1 = argsort(simvalsamps) | |
| sorti2 = argsort(lfsamps) | |
| return (interp1d(simvalsamps[sorti1],lfsamps[sorti1]), | |
| interp1d(lfsamps[sorti2],simvalsamps[sorti2])) | |
| def abund_match_data_to_data(data1,data2,cf1=1,cf2=1,reverse1=False,reverse2=False): | |
| """ | |
| Abundance matches two data sets assuming the same volume for each. Always | |
| matches cumulatively from *largest* value if `reverse` is False. | |
| `cf`s are completeness factors - e.g. if .5, that data set will have its | |
| abundance doubled, 2 halved, etc. | |
| """ | |
| from scipy.interpolate import interp1d | |
| sd1 = np.sort(data1) | |
| sd2 = np.sort(data2) | |
| c1 = (arange(sd1.size)[::-1]+1)/cf1 | |
| c2 = (arange(sd2.size)[::-1]+1)/cf2 | |
| if reverse1: | |
| d2c1 = interp1d(sd1,c1[::-1]) | |
| c2d1 = interp1d(c1[::-1],sd1) | |
| else: | |
| d2c1 = interp1d(sd1,c1) | |
| c2d1 = interp1d(c1[::-1],sd1[::-1]) | |
| if reverse2: | |
| d2c2 = interp1d(sd2,c2[::-1]) | |
| c2d2 = interp1d(c2[::-1],sd2) | |
| else: | |
| d2c2 = interp1d(sd2,c2) | |
| c2d2 = interp1d(c2[::-1],sd2[::-1]) | |
| #return d2c1,c2d1,d2c2,c2d2 | |
| one2two = lambda d1:c2d2(d2c1(d1)) | |
| two2one = lambda d2:c2d1(d2c2(d2)) | |
| return one2two,two2one | |
| def abund_match_w_errs(mfunc,mfsyserr,simvol,lfmag,lfphi,lferr,magval,n=100): | |
| """ | |
| abundance matches the given mass function to the given lum func with a | |
| systematic (fractional) shift in the MF. | |
| """ | |
| vels = [] | |
| for i in range(n): | |
| print 'Doing abundance matching iteration #',i+1 | |
| #mfuncjk = np.random.permutation(mfunc)[:-mfuncerr*mfunc.size]#jackknife | |
| mfuncwerr = mfunc*(1+mfsyserr*randn(1)[0]) | |
| lfoff = lferr*randn(lfphi.size) #lferr should match lfphi | |
| sim2lf,lf2sim = abund_match_sim_to_data(mfuncwerr,simvol,lfmag,lfphi+lfoff) | |
| vels.append(lf2sim(magval)) | |
| return array(vels) | |
| def abund_match_shuffle(mfunc,simvol,lfmag,lfphi,lferr,lscatter,magval=None,n=100): | |
| """ | |
| does abundance matching with an added (logarithmic) scatter in the luminosity | |
| function as mapped | |
| """ | |
| import gc | |
| vels = [] | |
| sim2lfs = [] | |
| lf2sims = [] | |
| for i in range(n): | |
| print 'doing',i+1,'of',n | |
| lfoff = lferr*randn(lfphi.size) | |
| sim2lf,lf2sim = abund_match_sim_to_data(mfunc,simvol,lfmag,lfphi+lfoff) | |
| mfunccut = mfunc[(sim2lf.x.min()<mfunc)&(mfunc<sim2lf.x.max())] | |
| lfofsim = sim2lf(mfunccut)#note that LF is in abs mags | |
| #lfofsim = 10**(log10(lfofsim)+lscatter*randn(lfofsim.size)) #this is for lum | |
| lfofsim = lfofsim + 2.5*lscatter*randn(lfofsim.size) | |
| sim2lf,lf2sim = abund_match_data_to_data(mfunccut,lfofsim,reverse2=True) | |
| if magval is None: | |
| sim2lfs.append(sim2lf) | |
| lf2sims.append(lf2sim) | |
| else: | |
| vels.append(lf2sim(magval)) | |
| gc.collect() | |
| if magval is not None: | |
| return array(vels) | |
| else: | |
| return sim2lfs,lf2sims | |
| def get_hiiex_cutouts(vagcd,vcfn='vischeck_z500_dr7_increasing', | |
| inds=[26,58,262,0,3,1],imsz=(256,256), | |
| scales=[.7,.5,5,.25,.5,.35],showpobj=[0,0,0,1,1,1], | |
| showsat=[0,0,0,1,1,1],action='dl'): | |
| from operator import isSequenceType | |
| from math import ceil | |
| from urllib import urlencode,urlretrieve,urlcleanup | |
| import webbrowser | |
| sorti = argsort(inds) | |
| hiis = get_vischeck_filtered(vcfn,lambda code,comment:'HII' in comment or 'overlapping' in comment,includenonvis=False)[0] | |
| hiis = get_vischeck_filtered(vcfn,includenonvis=True)[0] | |
| rap = vagcd['kcorr'].RA[hiis.i] | |
| decp = vagcd['kcorr'].DEC[hiis.i] | |
| ras = vagcd['kcorr'].RA[hiis.ni] | |
| decs = vagcd['kcorr'].DEC[hiis.ni] | |
| zs = vagcd['kcorr'].Z[hiis.i] | |
| dr = 7 | |
| mode = 'chart' | |
| if isscalar(scales) or scales is None: | |
| scales = [scales for i in ras] | |
| scales = [0.396127 if s is None else s for s in scales] | |
| showpobj = showpobj if isSequenceType(showpobj) else [showpobj for s in scales] | |
| showsat = showsat if isSequenceType(showsat) else [showsat for s in scales] | |
| urls = [] | |
| j=0 | |
| for i,(ra,dec,ra2,dec2,z) in enumerate(zip(rap,decp,ras,decs,zs)): | |
| if i is None or i in inds: | |
| print i,sorti | |
| scale = np.array(scales)[sorti][j] | |
| showp = np.array(showpobj)[sorti][j] | |
| shows = np.array(showsat)[sorti][j] | |
| j+=1 | |
| baseurl = 'http://casjobs.sdss.org/ImgCutoutDR%i/getjpeg.aspx?'%dr | |
| data = {'ra':ra,'dec':dec, | |
| 'width':imsz[0],'height':imsz[1], | |
| 'scale':scale, | |
| 'opt':'P' if showp else ''} | |
| if shows: | |
| linestr='\r\n' | |
| radecs = ['RA+DEC','{0}+{1}'.format(ra2,dec2)] | |
| data['query'] = linestr.join(radecs) | |
| url = baseurl+urlencode(data).replace('%2B','+') | |
| urls.append(url) | |
| urls = np.array(urls)[np.argsort(sorti)] | |
| if action=='dl' or 'combine' in action: | |
| fns = [] | |
| urlcleanup() | |
| for i,url in enumerate(urls): | |
| fns.append('hii%i.jpg'%(i+1)) | |
| print 'dling',url,'to',fns[-1] | |
| urlretrieve(url,fns[-1]) | |
| if 'combine' in action: | |
| lstr = action.replace('combine','') | |
| try: | |
| edge = int(lstr) | |
| except ValueError: | |
| edge = 1 | |
| from PIL import Image as pilim | |
| ims = [pilim.open(fn) for fn in fns] | |
| ncol = len(ims) if len(ims)<=3 else 3 | |
| nrow = 1 if len(ims)<=3 else int(ceil(len(ims)/3)) | |
| widths = [np.max([ims[i*nrow+j].size[0] for j in range(nrow)]) for i in range(ncol)] | |
| heights = [np.max([ims[i*ncol+j].size[1] for j in range(ncol)]) for i in range(nrow)] | |
| width = sum(widths)+(ncol-1)*edge | |
| height = sum(heights)+(nrow-1)*edge | |
| mainim = pilim.new('RGB',(width,height),'white') | |
| upper = -ims[0].size[1]-edge | |
| for i,im in enumerate(ims): | |
| if (i%ncol)==0: | |
| left = 0 | |
| upper += im.size[1]+edge | |
| mainim.paste(im,(left,upper)) | |
| left += im.size[0] | |
| left += edge | |
| print 'saving',mainim.size,'image to hiiex.png,.pdf,.eps' | |
| mainim.save('hiiex.png') | |
| mainim.save('hiiex.pdf') | |
| mainim.save('hiiex.eps') | |
| elif action=='open': | |
| for url in urls: | |
| webbrowser.open_new_tab(url) | |
| return list(urls) | |
| def select_match(dataval,datamatch,compval,compmatch,nmatch=20): | |
| """ | |
| data should be single val, `comp` should be matching length n, and | |
| `compmatch` must be sorted | |
| `nmatch` can be a number, or a fraction of the compmatch size. | |
| :returns: data-mean(compval in nmatch range of compmatch) | |
| """ | |
| if compval.shape != compmatch.shape: | |
| raise ValueError('compval and compmatch do not match') | |
| dataval = array(dataval,ndmin=1,copy=False) | |
| datamatch = array(datamatch,ndmin=1,copy=False) | |
| if dataval.shape != datamatch.shape: | |
| raise ValueError('dataval and datamatch do not match') | |
| if nmatch<1 and nmatch is not 1.0: | |
| nmatch = int(nmatch*compmatch.size) | |
| if nmatch>compmatch.size: | |
| raise ValueError('nmatch too large') | |
| lnmatch = int(nmatch/2) | |
| rnmatch = lnmatch+1 if lnmatch != nmatch/2 else lnmatch | |
| cens = searchsorted(compmatch,datamatch) | |
| ledge,redge = cens-lnmatch,cens+rnmatch | |
| redge[ledge<0] -= ledge[ledge<0] #makes the right edge *biger* if the left edge is too long | |
| ledge[ledge<0] = 0 | |
| ledge[redge>=compmatch.size] -= (redge[redge>=compmatch.size]-compmatch.size) | |
| redge[redge>=compmatch.size] = compmatch.size | |
| return dataval - array([mean(compval[li:ri]) for li,ri in zip(ledge,redge)]) | |
| if __name__ == '__main__': | |
| from optparse import OptionParser | |
| from os import uname | |
| op = OptionParser() | |
| op.usage = '%prog [options] [plot1 plot2 ...]' | |
| op.add_option('-m','--main-figs',help='Shows all main figures even if they are not included in command line', | |
| dest='mainfigs',action='store_true',default=False) | |
| op.add_option('-s','--save',help='save in the specified location',default=None) | |
| op.add_option('-o','--overwrite',help='Overwrite existing figure',action='store_true',default=False) | |
| op.add_option('-d','--dr',help='data release (6 or 7)',type='int',default=7) | |
| ops,args = op.parse_args() | |
| dr = ops.dr | |
| #on-run actions here | |
| from astropysics.phot import ugriz_to_UBVRI | |
| from astropysics.constants import H0 | |
| from astropysics.utils import funpickle | |
| h = H0/100 | |
| if dr == 6: | |
| cutsd = load_cuts_dr6() | |
| potsd = load_pots_dr6() | |
| elif dr == 7: | |
| cutsd = load_cuts_dr7() | |
| potsd = load_pots_dr7() | |
| else: | |
| raise ValueError('invalid data release '+str(dr)) | |
| _cut2fn = {'z500':'vischeck_z500_dr7_increasing', | |
| 'uniquenoz':'vischeck_uniquenoz_dr7_increasing', | |
| 'closernoz':'vischeck_closernoz_dr7_random', | |
| 'noz':'vischeck_noz_dr7_increasing'} | |
| hiicutsd,nhiicutsd = _load_vischeck_cutsds(cutsd,['HII','fg star'],_cut2fn) | |
| #a silly hack to not have to copy everything over on my home computer | |
| if 'teal' in uname()[1]: | |
| vagcd = load_VAGC('dr6' if dr == 6 else 'dr7',loadcat=True, | |
| loadimg=False, | |
| loadspec=False, | |
| loadkcnone=True, | |
| loadkcnear=False, | |
| loadcollnone=False, | |
| loadcollnear=False) | |
| else: | |
| vagcd = load_VAGC('dr6' if dr == 6 else 'dr7') | |
| zsim = load_betsy_zentner(disolimit=400) | |
| sanas = zsim[zsim.nearest_lmc_analog] | |
| sanah = zsim[zsim.analog_host] | |
| measured_varnames_sim = ['minsepkpch','truesepkpch','vel_of_minsep','deltav'] | |
| inferred_varnames_sim = ['halo_vin','vin_host','a_accreted','t_encounter_last', | |
| 'b_encounter_last','v_encounter_last','n_encounter', | |
| 'n_sat'] | |
| varnames_sim = [] | |
| varvals_sim = [] | |
| varnames_sim.extend(measured_varnames_sim) | |
| varvals_sim.extend([sanas[p] for p in measured_varnames_sim]) | |
| varnames_sim.append('ms+v') | |
| varvals_sim.append(((sanas.minsep/250)**2+(sanas.vel_of_minsep/500)**2)**0.5) | |
| varnames_sim.append('ms*v') | |
| varvals_sim.append(sanas.minsep*sanas.deltav) | |
| varnames_sim.append('ms*abs(v)') | |
| varvals_sim.append(sanas.minsep*sanas.vel_of_minsep) | |
| varnames_sim.append('abs(v)/ms') | |
| varvals_sim.append(sanas.vel_of_minsep/sanas.minsep) | |
| varobslast = len(varvals_sim)-1 | |
| varnames_sim.extend(inferred_varnames_sim) | |
| varvals_sim.extend([sanas[p] for p in inferred_varnames_sim]) | |
| varvals_sim = array(varvals_sim) | |
| #use as cc_plots(varvals_sim,varnames_sim,0) | |
| galmask = get_phot_mask(vagcd,rlim=(-17,-23.5),gmrlim=(.2,1),finitecut=True) | |
| # hostmask = zeros(vagcd['kcorr'].size,dtype=bool) | |
| # hostmask[cutsd['z500'].i] = True | |
| # satmask = zeros(vagcd['kcorr'].size,dtype=bool) | |
| # satmask[cutsd['z500'].ni] = True | |
| sep = cutsd['z500'].nd | |
| zvol = 0.0342 #chris says this is the volume limit | |
| ra,dec = vagcd['kcorr'].RA[galmask],vagcd['kcorr'].DEC[galmask] | |
| rah,dech = vagcd['kcorr'].RA[cutsd['z500'].i],vagcd['kcorr'].DEC[cutsd['z500'].i] | |
| ras,decs = vagcd['kcorr'].RA[cutsd['z500'].ni],vagcd['kcorr'].DEC[cutsd['z500'].ni] | |
| #"raw" - unfiltered by galmask | |
| ur,gr,rr,ir,zr,Jr,Hr,Kr = vagcd['kcorr'].ABSMAG.T + 5*np.log10(h) | |
| u,g,r,i,z,J,H,K = vagcd['kcorr'].ABSMAG[galmask].T + 5*np.log10(h) | |
| uh,gh,rh,ih,zh,Jh,Hh,Kh = vagcd['kcorr'].ABSMAG[cutsd['z500'].i].T + 5*np.log10(h) | |
| us,gs,rs,i_s,z_s,Js,Hs,Ks = vagcd['kcorr'].ABSMAG[cutsd['z500'].ni].T + 5*np.log10(h) | |
| b3h = vagcd['kcorr'].B300[cutsd['z500'].i] | |
| b3s = vagcd['kcorr'].B300[cutsd['z500'].ni] | |
| b10h = vagcd['kcorr'].B1000[cutsd['z500'].i] | |
| b10s = vagcd['kcorr'].B1000[cutsd['z500'].ni] | |
| mstarh = vagcd['kcorr'].MASS[cutsd['z500'].i] | |
| mstars = vagcd['kcorr'].MASS[cutsd['z500'].ni] | |
| gmr = g - r | |
| gmrh = gh - rh | |
| gmrs = gs - rs | |
| U,B,V,R,I = ugriz_to_UBVRI(u,g,r,i,z) | |
| Uh,Bh,Vh,Rh,Ih = ugriz_to_UBVRI(uh,gh,rh,ih,zh) | |
| Us,Bs,Vs,Rs,Is = ugriz_to_UBVRI(us,gs,rs,i_s,z_s) | |
| zs = vagcd['kcorr'].Z[galmask] | |
| zsh = vagcd['kcorr'].Z[cutsd['z500'].i] | |
| zss = vagcd['kcorr'].Z[cutsd['z500'].ni] | |
| #load bolshoi dists | |
| bolshoid = io.loadtxt_text_fields ('erik_first_sats.416.dat') | |
| msiicone = MSIICone('cone_B03_rotated.2011.01.20.npy') | |
| #for vmax>50 - note MSII is 100 h^-1 Mpc | |
| #do msiimassfuncs['Mvir']/h to get physical vir mass | |
| msiimassfuncs = np.load('z0infall_vmax50.npy') | |
| msiivtob03,b03tomsiiv = abund_match_sim_to_sch(msiimassfuncs['vmax'],1e6,.0149,-1.05,-20.44) #Blanton 03 LF | |
| b05lf = np.load('blanton05lf.npy') #r-band | |
| b05lfugiz = np.load('blanton05lfugiz.npy') | |
| msiivtob05,b05tomsiiv = abund_match_sim_to_data(msiimassfuncs['vmax'],1e6,b05lf['m'],b05lf['p1']) #Blanton 05 dwarf LF | |
| msiivtodr6lf,dr6lftomsiiv = abund_match_sim_to_sch(msiimassfuncs['vmax'],1e6,.0090,-1.23,-20.73)#Montero-Dorta DR6 LF | |
| #abund_match_w_errs(msiimassfuncs['vmax'],.1,1e6,b05lf['m'],b05lf['p1'],b05lf['e1'],#mag#) | |
| #abund_match_shuffle(msiimassfuncs['vmax'],1e6,b05lf['m'],b05lf['p1'],b05lf['e1'],lscatter,magval=None,n=100): | |
| if dr==7: | |
| liness,scutsi,linesh,hcutsi = match_line_query_to_cuts('lmcana_dr7z500hostlines_eteq.fits','lmcana_dr7z500satlines_eteq.fits',cutsd) | |
| if ops.mainfigs: | |
| figstomake = _mainfigs[:] | |
| figstomake.extend(args) | |
| else: | |
| figstomake = args | |
| if len(figstomake)>0: | |
| make_plots(set(figstomake),save=ops.save,overwrite=ops.overwrite) | |
| #DR6 errors in redshift: hosts median=9.4205233835964464e-06 | |
| # sats median=1.0891547390201595e-05 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment