Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix quickquasars targetid truth #457

Merged
merged 9 commits into from Dec 4, 2018
45 changes: 24 additions & 21 deletions py/desisim/scripts/quickquasars.py
Expand Up @@ -359,24 +359,22 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
for ii in range(len(metadata)):

# quasars with z < min_z will not have any DLA in spectrum
if min_lya_z > metadata[ii]['Z']:
continue
if min_lya_z>metadata['Z'][ii]: continue

# quasar ID
idd=metadata['MOCKID'][ii]
dlas=[]

if args.dla=='file':

for dla in dla_info[dla_info['MOCKID']==idd]:

# Adding only DLAs with z < zqso
if (dla['Z_DLA']< metadata[ii]['Z']):
dlas.append(dict(z=dla['Z_DLA']+dla['DZ_DLA'],N=dla['N_HI_DLA']))
if dla['Z_DLA']>=metadata['Z'][ii]: continue
dlas.append(dict(z=dla['Z_DLA']+dla['DZ_DLA'],N=dla['N_HI_DLA']))
transmission_dla = dla_spec(trans_wave,dlas)

elif args.dla=='random':

dlas, transmission_dla = insert_dlas(trans_wave, metadata[ii]['Z'], rstate=random_state_just_for_dlas)
dlas, transmission_dla = insert_dlas(trans_wave, metadata['Z'][ii], rstate=random_state_just_for_dlas)

# multiply transmissions and store information for the DLA file
if len(dlas)>0:
Expand All @@ -388,9 +386,9 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
# write file with DLA information
if len(dla_id)>0:
dla_meta=Table()
dla_meta['NHI']=dla_NHI
dla_meta['z']=dla_z
dla_meta['ID']=dla_id
dla_meta['NHI'] = dla_NHI
dla_meta['Z'] = dla_z
dla_meta['TARGETID'] = dla_id

hdu_dla = pyfits.convenience.table_to_hdu(dla_meta)
hdu_dla.name="DLA_META"
Expand Down Expand Up @@ -462,6 +460,8 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
lyaforest=False, nocolorcuts=True,
noresample=True, seed=seed, south=issouth)

_meta['TARGETID'][:] = metadata['MOCKID']
_qsometa['TARGETID'][:] = metadata['MOCKID']
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Those two lines are what fixes the original bug?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, it is. Plus some other for DLA and BAL.

meta[these] = _meta
qsometa[these] = _qsometa
tmp_qso_flux[these, :] = _tmp_qso_flux
Expand Down Expand Up @@ -492,6 +492,9 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
# restore random state to get the same random numbers later
# as when we don't insert BALs
np.random.set_state(rnd_state)
meta_bal['TARGETID'] = metadata['MOCKID']
w = meta_bal['TEMPLATEID']!=-1
meta_bal = meta_bal[:][w]
londumas marked this conversation as resolved.
Show resolved Hide resolved
hdu_bal=pyfits.convenience.table_to_hdu(meta_bal); hdu_bal.name="BAL_META"
del meta_bal
else:
Expand Down Expand Up @@ -606,7 +609,7 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
meta=specmeta,seed=seed,fibermap_columns=fibermap_columns,use_poisson=False) # use Poisson = False to get reproducible results.


##Adedd to write the truth file, includen metadata for DLA's and BALs
##Adedd to write the truth file, includen metadata for DLA's and BALs
log.info('Writing a truth file {}'.format(truth_filename))
meta.rename_column('REDSHIFT','TRUEZ')
meta.add_column(Column(Z_noFOG,name='TRUEZ_noFOG'))
Expand All @@ -624,9 +627,9 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
hduqso.header['EXTNAME'] = 'QSO_META'
hdulist=pyfits.HDUList([pyfits.PrimaryHDU(),hdu,hduqso])
if args.dla:
hdulist.append(hdu_dla)
hdulist.append(hdu_dla)
if args.balprob:
hdulist.append(hdu_bal)
hdulist.append(hdu_bal)
hdulist.writeto(truth_filename, overwrite=True)
hdulist.close()

Expand All @@ -649,20 +652,20 @@ def simulate_one_healpix(ifilename,args,model,obsconditions,decam_and_wise_filte
('DELTACHI2', 'f8'),
('BRICKNAME', (str,8))]
zbest = Table(np.zeros(nqso, dtype=columns))
zbest["CHI2"][:] = 0.
zbest["Z"] = metadata['Z']
zbest["ZERR"][:] = 0.
zbest["CHI2"][:] = 0.
zbest["Z"][:] = metadata['Z']
zbest["ZERR"][:] = 0.

if args.gamma_kms_zfit:
log.info("Added zfit error with sigma {} to zbest".format(args.gamma_kms_zfit))
dz_fit=mod_cauchy(loc=0,scale=args.gamma_kms_zfit,size=nqso,cut=3000)*(1.+metadata['Z'])/c
zbest["Z"]+=dz_fit

zbest["ZERR"][:] = 0
zbest["ZWARN"][:] = 0
zbest["SPECTYPE"][:] = "QSO"
zbest["SUBTYPE"][:] = ""
zbest["TARGETID"] = fibermap["TARGETID"]
zbest["ZERR"][:] = 0
zbest["ZWARN"][:] = 0
zbest["SPECTYPE"][:] = 'QSO'
zbest["SUBTYPE"][:] = ''
zbest["TARGETID"][:] = metadata['MOCKID']
zbest["DELTACHI2"][:] = 25.
hzbest = pyfits.convenience.table_to_hdu(zbest); hzbest.name="ZBEST"
hfmap = pyfits.convenience.table_to_hdu(fibermap); hfmap.name="FIBERMAP"
Expand Down