<jupyter_start><jupyter_text>I think ideally, I should've removed the red geyser sample from datas right here, before looking for control sample.<jupyter_code>from astropy.table import Table
control_inds = np.array([]) #we'll save the indices of your control galaxes here
dm=0.1
for mass in np.log10(data['nsa_elpetro_mass']):
    #find the indices of galaxies which satisfy your criteria
    sel_control = np.where((np.log10(datas['nsa_elpetro_mass_2']) > (mass-0.1)) & (np.log10(datas['nsa_elpetro_mass_2']) < (mass+0.1)))
    #(you could add NUV-r > 5 reqiurement as well here )
    sel_control = sel_control[0]
    control_inds=np.concatenate([control_inds,sel_control])


res=[]
[res.append(x) for x in control_inds if x not in res] #getting rid of redundancies in indices

hmm=np.array(res) #turn list into array
print('The indices of the control galaxies:', len(res),'. I don"t know if the fact that this output is 50 less than 749 is a coincidence or not.')
t=Table(datas) #attempting to apply index slicing
t.add_index('nsa_elpetro_mass_2')
res2=list(map(int, res))
control=t[res2]

print(control)
control.write('controlll-take-22.fits', overwrite=True) #export this control table that still has some red geysers in it
cntrl = 'controlll-take-22.fits'
hdul6 = fits.open(cntrl)
data6 = hdul6[1].data

hdu6 = hdul6[1]
hdr6 = hdul6[0].header 
hdr6 = hdul6[1].header
#I'm pretty sure this cell is obsolete

last='control-geysers3.fits'

#hdul7 = fits.open(last)
#hdul7.info()
#data7 is the final control group without any of the red geyser sample, even though TOPCAT matching only removed 44 galaxies from controlll-take-3
#hdu7 = hdul7[1]
#hdr7 = hdul7[0].header 
#hdr7 = hdul7[1].header
#data7 = hdu7.data
#print(len(data7['nsa_elpetro_mass_1']))
#detections
mask7=(data6['LOGMHI_2'] > -999) & (data6['Z_2'] < 0.05)
mask10=(data6['LOGMHI_2'] < -999) & (data6['LOGHILIM200KMS_2'] > -999)
#print(len(data6['nsa_elpetro_mass_2'][mask10]))
xdetx=data6['nsa_elpetro_mass_2'][mask7]
ydety=(10**(data6['LOGMHI_2'][mask7]))/(xdetx)
print('Number of detected control galaxies =',len(xdetx))
#one galaxy (8082-6101) is neither non-detection nor detection
#nondetections
mask8=(data6['LOGHILIM200KMS_2'] > -999) & (data6['Z_2'] < 0.05)
xnondetx=data6['nsa_elpetro_mass_2'][mask8]
ynondety=(10**(data6['LOGHILIM200KMS_2'][mask8]))/(xnondetx)
print('Number of non-detected control galaxies =',len(xnondetx))
#print('Number of detected control galaxies =',len(ynondety))
import matplotlib.pyplot as plt
plt.figure(figsize=(8,8))
#ax = fig.add_axes([0.1, 0.1, 0.6, 0.75])
plt.scatter(xdetx, ydety, label='control detections', color='darkslategrey', marker='x',s=30)
plt.scatter(xnondetx, ynondety,label='control non-detections', color='silver', marker='o', s=30)
plt.scatter(detx,dety, label='geyser detections', color= 'fuchsia', marker='x', s=50)
plt.scatter(nondetx, nondety, label= 'geyser non-detections', color='red', marker='o', s=40)
plt.xscale('log')
plt.yscale('log')
plt.xlabel("log $M_* [M_o]$", fontsize=15)
plt.ylabel("log $M_{HI}/M_*$",fontsize=15)
plt.legend(loc='lower left', fontsize=15)
plt.savefig('all_mass_plot.png', bbox_inches='tight')
plt.title("Control Group and Geyser Sample $M_{HI}/M_*$ vs $M_*$", fontsize=15)
#combine red geyser table (data) with control sample table (controlll-take-22.fits)
file2='all-finall.fits'

hdul8 = fits.open(file2)
hdul8.info()

hdu8 = hdul8[1]
hdr08 = hdul8[0].header 
hdr8 = hdul8[1].header
alldata = hdu8.data
#make column for detection
hip = np.array(alldata['LOGMHI'])
sel = (alldata['LOGMHI'] > -999) 
sel = np.multiply(sel, 1)
print(sel)
#import sys
#np.set_printoptions(threshold=sys.maxsize)
#print(sel)
#print(alldata['LOGMHI'])
#make column for sample
bet=alldata['NUV_r']
bet[np.isnan(bet)] = 0
sel2 = bet > 0
sel2 = np.multiply(sel2, 1)

#where_are_NaNs = isnan(alldata['NUV_r'])
#alldata['NUV_r'][where_are_NaNs] = 0
bet[np.isnan(bet)] = 0
print(sel2)
#make new column for MHI/M* take two
ratio = np.zeros(len(alldata))
nondetections = alldata['LOGHILIM200KMS'] > -999
detections = (nondetections == False)
ratio[nondetections] = (10**alldata['LOGHILIM200KMS'][nondetections])/(alldata['nsa_elpetro_mass'][nondetections])
ratio[detections] = (10**alldata['LOGMHI'][detections])/(alldata['nsa_elpetro_mass'][detections])
print(len(ratio))
#table for Dave M*, MHI/M*, detection (1=True, 0=False), and sample (1=red geyser, 0=control)

tebl=Table()
tebl['MANGAID'] = alldata['MANGAID_1']
tebl['M*'] = alldata['nsa_elpetro_mass']
tebl['MHI/M*'] = ratio
tebl['Detection (0=non)'] = sel
tebl['Sample (0=red geyser)'] = sel2
tebl.remove_row(195)
tebl.show_in_notebook()
tebl.write('stattable2.fits', format='fits')<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup