;+ ; NAME: ; mdrp_extract_object ; ; PURPOSE: ; Performs all object extraction tasks ; 0) Locate bright fibers, and test image background ; 1) 4 Step Optimal extraction ; ; CALLING SEQUENCE: ; mdrp_extract_object, outname, objhdr, image, invvar, slitmap, wset, $ ; xarc, lambda, xtrace, fflat, fibermask, proftype=, color=, $ ; [ widthset=, dispset=, skylinefile=, plottitle=, superflatset=, $ ; /splitsky, ccdmask= , slithdr=, opfibers=] ; ; INPUTS: ; outname - Name of outputs FITS file ; objhdr - Header cards from object image ; image - Object image [nx,ny] ; invvar - Inverse Variance of object image [nx,ny] ; slitmap - Slitmap ; wset - Wavelength solution from arc line spectra ; xarc - centroids measured in arc line spectra ; lambda - air wavelengths corresponding to xarc ; xtrace - spatial traces from flat field ; fflat - 1d flat field vectors ; fibermask - Fiber status bits, set nonzero for bad status [NFIBER] ; proftype - Which type of profile should we use, (default=1 Gaussian) ; superflatset- If present, then divide by median superflat! ??? ; slithdr - Yanny header from the slitmap ; color - ??? ; widthset - ??? ; dispset - ??? ; opfibers - opFibers structure trimmed to this spectrograph ; ; REQUIRED KEYWORDS: ; color - camera color (red or blue) ; ; OPTIONAL KEYWORDS: ; plottitle - Prefix for titles in QA plots. ; ccdmask - If set, then use this to set some pixel values in pixmask ; ; OUTPUTS: ; A fits file is output in outname, which contains ; FLOAT flux [NX,NTRACE] ; FLOAT flux_invvar [NX,NTRACE] ; LONG manga pixelmask [NX,NTRACE] ; STRUCT vacuum wavelengths ; STRUCT wavelength dispersion ; STRUCT slitmap [NTRACE] ; ; OPTIONAL OUTPUTS: ; ; COMMENTS: ; ; EXAMPLES: ; ; BUGS: ; ; PROCEDURES CALLED: ; mdrp_divideflat ; djs_median() ; djs_oplot ; djs_plot ; extract_boxcar() ; extract_image ; drp2pixmask_bits() ; fitsn() ; fitvacset() ; get_tai ; heliocentric() ; mwrfits ; drp2pixmask_bits() ; qaplot_scatlight ; qaplot_skydev ; spadd_guiderinfo ; splog ; sxaddpar ; sxpar() ; traceset2xy ; find_whopping() ; ; INTERNAL SUPPORT ROUTINES: ; ; REVISION HISTORY: ; 24-Jan-2000 Written by S. Burles, Chicago ; 26-Jul-2001 Also pass proftype and superflatset ; 29-Mar-2011 Switching to pure IDL bundle-wise extraction (A. Bolton, U. Utah) ; 15-Aug-2011 Enabling split sky model across halves of CCD. (A. Bolton, U. Utah) ; Nov-2012 Modified for MaNGA survey (B. Cherinka, Toronto) ; 15-Jul-2013 Separated out from Boss (renamed mdrp_), removed all sky subtraction stuff ; 15-May-2014 Added quality control bitmask logic (David Law) ; 11-Sep-2014 Migrated drpqual to drp2qual (D. Law) ; 14-Jan-2015 Retool to use opfibers (D. Law) ; 31-Aug-2015 Use SKYDISPSET for science frames (D. Law) ;- ;------------------------------------------------------------------------------ pro mdrp_extract_object, outname, objhdr, image, invvar, slitmap, wset, $ xarc, lambda, xtrace, fflat, fibermask, color=color, proftype=proftype, skylinefile=skylinefile, $ width=width, dispset=dispset, plotname=plotname, plottitle=plottitle, superflatset=superflatset, $ ccdmask=ccdmask, slithdr=slithdr, VISUAL=VISUAL, SURVEY=survey, opfibers=opfibers, writefiles=writefiles on_error,0 compile_opt idl2 if (keyword_set(visual)) then begin set_plot,'x' device,decomposed=0 loadct,0,/silent !p.multi=[0,1,2] !p.thick=0 & !x.thick=0 & !y.thick=0 & !p.charthick=0 endif else if (keyword_set(plotname)) then begin set_plot,'ps' origp=!p origx=!x origy=!y !p.multi=[0,1,2] !p.thick=0 & !x.thick=0 & !y.thick=0 & !p.charthick=0 dfpsplot, plotname, /color endif configuration=obj_new('configuration', sxpar(objhdr, 'MJD')) objname = strtrim(sxpar(objhdr,'OBJFILE'),2) flavor = strtrim(sxpar(objhdr,'FLAVOR'),2) camera = strtrim(sxpar(objhdr,'CAMERAS'),2) spectrographid=fix(strmid(camera,1,1)) partialslit=slitmap[where(slitmap.spectrographid eq spectrographid)] ; Work out some bundle and fiber info for Brian's old QA routines ; and checks totfiber=n_elements(opfibers) uniqbundles=opfibers[uniq(opfibers.bundleid)].bundleid nbundle=n_elements(uniqbundles) nfiber=intarr(nbundle) bundleid=intarr(nbundle) radius=fltarr(nbundle) for i=0,nbundle-1 do begin index=where(opfibers.bundleid eq uniqbundles[i]) nfiber[i]=opfibers[index[0]].blocksize bundleid[i]=opfibers[index[0]].bundleid radius[i]=2. endfor ;------------------ ; Identify very bright objects ; Do a boxcar extraction, and look for fibers where the median counts are 10000 ADU per row. fextract = ml_extract_boxcar(image*(invvar GT 0), xtrace, nfiber=nfiber, radius=radius) scrunch = djs_median(fextract, 1) ; Find median counts/row in each fiber whopping = find_whopping(scrunch, 10000.0, whopct) scrunch_sort = sort(scrunch) i5 = n_elements(scrunch)/20 i95 = i5 * 19 splog, 'Whopping fibers: ', whopping splog, 'Median counts in all fibers = ', djs_median(scrunch) splog, 'Number of bright fibers = ', whopct if (whopct GT 20) then begin splog, 'WARNING: Disable whopping terms ' + objname whopping = -1 whopct = 0 endif ;------------------------------------------------------------ ; Check for bad pixels within 3 pixels of trace badcheck = ml_extract_boxcar((invvar LE 0), xtrace, radius=2.5, nfiber=nfiber) badplace = where(badcheck GT 0) nx = (size(fextract,/dim))[0] ny = (size(fextract,/dim))[1] pixelmask = lonarr(nx,ny) badcolumns = where(total(badcheck GT 0,1) GT 0.45 * nx) if (badplace[0] NE -1) then pixelmask[badplace] = pixelmask[badplace] OR drp2pixmask_bits('NEARBADPIXEL') if (badcolumns[0] NE -1) then fibermask[badcolumns] = fibermask[badcolumns] OR drp2pixmask_bits('MANYBADCOLUMNS') if (whopping[0] NE -1) then begin ; Set the mask bit for whopping fibers themselves fibermask[whopping] = fibermask[whopping] OR drp2pixmask_bits('WHOPPER') ; Set the mask bit for fibers near whopping fibers, excluding the whopping fibers themselves. Note that a fiber could still have both ; WHOPPER and NEARWHOPPER set if it is both bright and near another bright fiber. wp = [whopping - 2 , whopping -1, whopping+1 , whopping+2] wp = wp[ where(wp GE 0 AND wp LT ny) ] fibermask[wp] = fibermask[wp] OR drp2pixmask_bits('NEARWHOPPER') endif ;----- ; Inherit any mask bits from the ccdmask, by setting the pixmask bits for anything that would be hit in a boxcar extraction if (keyword_set(ccdmask)) then begin for ibit=0, 31 do begin thischeck = ml_extract_boxcar((ccdmask AND 2L^ibit) NE 0, xtrace, radius=2.5, nfiber=nfiber) pixelmask = pixelmask OR (2L^ibit * (thischeck GT 0)) endfor endif ;----------------------------------------------------------------------- ; This is a kludge to fix first and last column ??? ;----------------------------------------------------------------------- if (configuration->extract_object_fixcolumns()) then begin image[0,*] = image[0,*]*0.7 image[nx-1,*] = image[nx-1,*]*0.7 endif ; ; First we should attempt to shift trace to object flexure xnow = mdrp_match_trace(image, invvar, xtrace,radius=radius, nfiber=nfiber) bestlag = median(xnow-xtrace) splog, 'Shifting traces by match_trace ', bestlag if (abs(bestlag) GT 1.0) then begin splog, 'WARNING: pixel shift is large!' endif highrej = 10 ; just for first extraction steps lowrej = 10 ; just for first extraction steps ; We need to check npoly with new scattered light backgrounds npoly = 16 ; maybe more structure, lots of structure nrow = (size(image))[2] yrow = lindgen(nrow) nfirst = n_elements(yrow) splog, 'Extracting frame '+objname+' with 4 step process' sigma2=width ntrace = (size(sigma2,/dimens))[1] wfixed = [1,1] ; Fit gaussian height + width (fixed center position) nterms = n_elements(wfixed) ;----------------------------------------------------------------------- ; Now, subtract halo image and do final extraction with all rows ;----------------------------------------------------------------------- ; Scattered light correction if we think it is needed scat_warning=0 ; Define a mask of the 150 pixel wide boundary of the detector boundmask=intarr(size(image,/dimension)) boundmask[150:(size(image,/dimension))[0]-150,150:(size(image,/dimension))[1]-150]=1 theedge=where(boundmask eq 0) scat_warning=0 if ((median(image) lt 300)and(median(image[theedge]) gt 30)) then begin splog,'WARNING- bad scattered light, applying secondary scattered light model' ; Scat_warning=1 if this routine runs successfully, =2 if it fails image=mdrp_scattered(image,invvar,xnow,scatmodel=scatmodel,nbin=5,nord=4,nreq=400,maskrad=5,upper=5.,lower=5.,badthresh=1.5,xstep1=50,xstep2=500,errflag=scat_warning) endif ; (6) Final extraction splog, 'Step 6: Final Object extraction' ; DRL set these parameters May 22 2014 ; These give truly flat 30-sec flatfields?, and also good ; flat sky-subtracted data for the dark all-sky plate ; 7341-56693. ;highrej = 8 ;lowrej = 5 ; RY (Mar 1 2015) set these higher for blue and red camera to deal with bright standard stars if strmatch(color,'blue',/fold_case) then highrej=15 else highrej = 12 if strmatch(color,'blue',/fold_case) then lowrej=15 else lowrej = 12 wfixed = [1] ; Fit gaussian amplitude only (not width or centroid) nterms = n_elements(wfixed) reject=[0.1, 0.6, 0.6] npoly = 1L extract_image, image, invvar, xnow, sigma2, flux, fluxivar, proftype=proftype, wfixed=wfixed, ansimage=ansimage3, highrej=highrej, lowrej=lowrej, npoly=npoly, $ chisq=chisq, ymodel=ymodel, pixelmask=pixelmask, reject=reject, /relative, visual=visual, survey=survey ;write out post-extracted flux (uncalibrated) and ymodel from science extraction if keyword_set(writefiles) then begin sciextname = ml_strreplace(outname,'mgFrame','mgFrame-flux') mwrfits, flux, sciextname, /create scimodname=ml_strreplace(outname,'mgFrame','mgFrame-model') mwrfits, ymodel, scimodname, /create spawn, ['gzip', '-f', sciextname], /noshell spawn, ['gzip', '-f', scimodname], /noshell endif ;---------------------------------------------------------------------- ; Can we find cosmic rays by looking for outlandish ansimage ratios??? ; a = where(ansimage[lindgen(ntrace)*nterms, *] LT (-2*ansimage[lindgen(ntrace)*nterms+1, *]) ;------------------ ; QA chisq plot for fit calculated in extract image (make QAPLOT ???) xaxis = lindgen(n_elements(chisq)) + 1 ymax = 2.*median(chisq) djs_plot, xaxis, chisq, xrange=[0,N_elements(chisq)], xstyle=1, yrange=[0,ymax], ystyle=1, xtitle='Row number', ytitle = '\chi^2', title=plottitle+'Extraction chi^2 for '+objname djs_oplot, !x.crange, [1,1] ;x djs_oplot, yrow, secondchisq[yrow], color='blue' ;x djs_oplot, yrow, firstchisq[yrow], color='green' xyouts, 100, 0.05*!y.crange[0]+0.95*!y.crange[1], 'BLACK = Final chisq extraction' ;x xyouts, 100, 0.08*!y.crange[0]+0.92*!y.crange[1], 'BLUE = Initial-scattered chisq extraction' ;x xyouts, 100, 0.08*!y.crange[0]+0.89*!y.crange[1], 'GREEN = Initial chisq extraction' ;------------------ ; Flat-field the extracted object fibers with the global flat mdrp_divideflat, flux, fflat, invvar=fluxivar, /quiet pixelmask = pixelmask OR ((fflat LT 0.5) * drp2pixmask_bits('LOWFLAT')) ;------------------ ; Tweak up the wavelength solution to agree with the sky lines. ; DRL- Set reasonable wavelength ranges empirically for what ; lines to use for MaNGA if (color eq 'blue') then maxwave=6315. if (color eq 'red') then minwave=6000. locateskylines, skylinefile, flux, fluxivar, wset, xarc, arcshift=arcshift, xsky=xsky, skywaves=skywaves, skyshift=skyshift, radius=radius, nfiber=nfiber, fibermask=fibermask,maxwave=maxwave,minwave=minwave if keyword_set(VISUAL) then getwindow,/open qaplot_skyshift, wset, xsky, skywaves, skyshift, title=plottitle+'Sky Line Deviations for '+objname if ~keyword_set(arcshift) then splog, 'WARNING: Cannot shift to sky lines' ;------------------ ; Fit for the widths of the sky lines (relative to those of the arcs) ; The values in XSKY are noisy measurements, so replace them with ; the predicted positions based upon the arc solution. xsky = transpose(traceset2pix(wset, alog10(skywaves))) ; We should also apply the arcshift here xsky += skyshift ; Derive skydispset skydispset = skyline_dispersion(flux, fluxivar, xsky, dispset, skywaves, camera, everyn=150, quadavg=quadavg,fibermask=fibermask) splog, 'Applying skyline adjusted line widths' ;------------------ ; Apply heliocentric correction ; Correct LAMBDA, which is used to shift to vacuum wavelengths. helio=0.0 ra = sxpar(objhdr, 'RA', count=ct_ra) dec = sxpar(objhdr, 'DEC', count=ct_dec) if (ct_ra NE 1 OR ct_dec NE 1) then splog, 'WARNING: Missing RA and/or DEC from header' ;-------------------------------------------------------- ; Call standard proc to determine time-stamps get_tai, objhdr, tai_beg, tai_mid, tai_end ; Set TAI equal to the time half-way through the exposure ; If all these keywords are present in the header, they will be either type FLOAT or DOUBLE. Note that SDSS will put NaN in the header for these values if they are unknown. if ( size(ra, /tname) NE 'INT' AND size(dec, /tname) NE 'INT' AND size(tai_mid, /tname) NE 'INT' AND finite(ra) AND finite(dec) AND finite(tai_mid) ) then begin helio = heliocentric(ra, dec, tai=tai_mid) splog, 'Heliocentric correction = ', helio, ' km/s' sxaddpar, objhdr, 'HELIO_RV', helio, ' Heliocentric correction (added to velocities)' endif else begin splog, 'WARNING: Header info not present to compute heliocentric correction' endelse if (size(tai_mid, /tname) EQ 'INT' OR finite(tai_mid) EQ 0) then begin splog, 'WARNING: Header info not present to compute airmass correction to sky level' tai_mid = 0 endif ;------------------ ; Shift to skylines and fit to vacuum wavelengths vacset = fitvacset(xarc, lambda, wset, arcshift, fibermask=fibermask, helio=helio, airset=airset) ; No longer make the following QA plot ??? ;qaplot_skydev, flux, fluxivar, vacset, plugsort, color, title=plottitle+objname sxaddpar, objhdr, 'VACUUM', 'T', ' Wavelengths are in vacuum' ;------------------ ; If present, reconstruct superflat and normalize sxaddpar, objhdr, 'SFLATTEN', 'F', ' Superflat has not been applied' superfit = 0 if keyword_set(superflatset) AND keyword_set(airset) then begin superfit = float(smooth_superflat(superflatset, airset, plottitle=plottitle+'Smooth superflat for '+objname)) if keyword_set(superfit) then begin mdrp_divideflat, flux, superfit, invvar=fluxivar, /quiet sxaddpar, objhdr, 'SFLATTEN', 'T', ' Superflat has been applied' endif endif ; Save current pixelmask for later MaNGA use mangamask=pixelmask ;---------- ; Add keywords to object header objhdr=mdrp_makeheader(head=objhdr,/drp2d) if (keyword_set(osigma)) then sxaddpar, objhdr, 'OSIGMA', sigma, ' Original guess at spatial sigma in pix' sxaddpar, objhdr, 'PREJECT', reject[1], ' Profile area rejection threshold' sxaddpar, objhdr, 'LOWREJ', lowrej, ' Extraction: low rejection' sxaddpar, objhdr, 'HIGHREJ', highrej, ' Extraction: high rejection' sxaddpar, objhdr, 'SCATPOLY', npoly, ' Extraction: Order of scattered light polynomial' sxaddpar, objhdr, 'PROFTYPE', proftype, ' Extraction profile: 1=Gaussian' sxaddpar, objhdr, 'NFITPOLY', nterms, ' Extraction: Number of parameters in each profile' sxaddpar, objhdr, 'XCHI2', mean(chisq), ' Extraction: Mean chi^2' sxaddpar, objhdr, 'DESIGNID', yanny_par(slithdr,'designid'), ' Design ID' sxaddpar, objhdr, 'NAXIS1', (size(flux))[1], ' X Size of flux ext' sxaddpar, objhdr, 'NAXIS2', (size(flux))[2], ' Y Size of flux ext' sxaddpar, objhdr, 'CENRA', yanny_par(slithdr,'racen'), ' Plate center RA' sxaddpar, objhdr, 'CENDEC', yanny_par(slithdr,'deccen'), ' Plate center DEC' ; Add informative comment on the AZ keyword sxaddpar, objhdr, 'AZ', sxpar(objhdr,'AZ'), 'Azimuth axis pos. (approx, deg, 0=S, 90=E)' sxaddpar,objhdr,'EQUINOX',2000.0,after='DEC' sxaddpar,objhdr,'RADESYSa', 'FK5', after='EQUINOX' sxaddpar,objhdr,'AIRMASS', float(tai2airmass(ra, dec, tai=tai_mid)) * (ct_ra EQ 1) * (ct_dec EQ 1) * (tai_mid GT 0), after='ALT' ; Compute beginning, midpoint, and end airmass, lst, hour angle airmsbeg= float(tai2airmass(ra, dec, tai=tai_beg, ha=hrangbeg,lstdeg=lstbeg)) * (ct_ra EQ 1) * (ct_dec EQ 1) * (tai_beg GT 0) airmsmid= float(tai2airmass(ra, dec, tai=tai_mid, ha=hrangmid,lstdeg=lstmid)) * (ct_ra EQ 1) * (ct_dec EQ 1) * (tai_mid GT 0) airmsend= float(tai2airmass(ra, dec, tai=tai_end, ha=hrangend,lstdeg=lstend)) * (ct_ra EQ 1) * (ct_dec EQ 1) * (tai_end GT 0) sxaddpar,objhdr,'AIRMSBEG',airmsbeg,' Plate airmass at beginning of exposure' sxaddpar,objhdr,'AIRMSMID',airmsmid,' Plate airmass at midpoint of exposure' sxaddpar,objhdr,'AIRMSEND',airmsend,' Plate airmass at end of exposure' sxaddpar,objhdr,'LSTBEG',lstbeg/15.,' LST (hours) at beginning of exposure' sxaddpar,objhdr,'LSTMID',lstmid/15.,' LST (hours) at midpoint of exposure' sxaddpar,objhdr,'LSTEND',lstend/15.,' LST (hours) at end of exposure' sxaddpar,objhdr,'HRANGBEG',hrangbeg/15.,' Plate hour angle (hours) at beginning of exposure' sxaddpar,objhdr,'HRANGMID',hrangmid/15.,' Plate hour angle (hours) at midpoint of exposure' sxaddpar,objhdr,'HRANGEND',hrangend/15.,' Plate hour angle (hours) at end of exposure' ; Add information about mask name to header sxaddpar,objhdr,'MASKNAME','MANGA_DRP2PIXMASK','Bits in sdssMaskbits.par used by mask extension' sxaddpar, objhdr, 'EXTEND', 'T', after='NAXIS2' ; Set ivar to zero where NODATA, FULLREJECT, or LOWFLAT maskbits set index=where((pixelmask and sdss_flagval('MANGA_DRP2PIXMASK','NODATA')) $ or (pixelmask and sdss_flagval('MANGA_DRP2PIXMASK','FULLREJECT')) $ or (pixelmask and sdss_flagval('MANGA_DRP2PIXMASK','LOWFLAT')) ne 0,nindex) if (nindex ne 0) then fluxivar[index]=0. ; Interpolate over masked pixels for cosmetic purposes flux=djs_maskinterp(flux,fluxivar EQ 0, /const, iaxis=0) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Quality control ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Quality control: get DRP2QUAL bitmask from header drp2qual=long(fxpar(objhdr,'DRP2QUAL')) ; Set file is valid bit drp2qual=drp2qual OR sdss_flagval('MANGA_DRP2QUAL','VALIDFILE') ; Look for flexure warning flag (> 0.5 pixel ; based on year 1 data) if (quadavg gt 0.5) then drp2qual=drp2qual OR sdss_flagval('MANGA_DRP2QUAL','BADFLEXURE') sxaddpar,objhdr, 'FQUADAVG', quadavg, ' Average flexure quadrature term' ; Look for scattered light warning flags if (scat_warning ne 0) then drp2qual=drp2qual OR sdss_flagval('MANGA_DRP2QUAL','HIGHSCAT') if (scat_warning eq 2) then drp2qual=drp2qual OR sdss_flagval('MANGA_DRP2QUAL','SCATFAIL') ; Look for cases where less than 30% of the pixels have good value index=where(fluxivar ne 0,nindex) if (nindex/float(n_elements(fluxivar)) lt 0.3) then $ drp2qual=drp2qual OR sdss_flagval('MANGA_DRP2QUAL','EXTRACTBAD') ; Flag case where extracted fluxes are huge ; (e.g., perhaps this is a twilight flat) if (nindex gt 0) then begin if (median(flux[index]) gt 1e4) then $ drp2qual=drp2qual OR sdss_flagval('MANGA_DRP2QUAL','EXTRACTBRIGHT') endif ; Flag cases where the exposure time was abnormally short ; (less than 5 minutes) or long (over 20 minutes) if (float(sxpar(objhdr,'EXPTIME')) lt 300.) then $ drp2qual=drp2qual OR sdss_flagval('MANGA_DRP2QUAL','LOWEXPTIME') if (float(sxpar(objhdr,'EXPTIME')) gt 1200.) then $ drp2qual=drp2qual OR sdss_flagval('MANGA_DRP2QUAL','LOWEXPTIME') ; Loop over harnesses and look for entire missing IFUs ; (e.g., that might have fallen out) and IFUs with really ; bad throughput variations badifu='' harnames=partialslit[uniq(partialslit.harname)].harname nhar=n_elements(harnames) for i=0,nhar-1 do begin ; Identify the IFU fibers in this block thesefibers=where((partialslit.harname eq harnames[i])and(partialslit.fibertype eq 'IFU'),nthese) ; Sum the fluxes (silly way to do it, but it works) theflux=total(flux[*,thesefibers]) if (theflux eq 0.) then begin splog, strcompress('WARNING! Missing all IFU fibers in harness '+harnames[i]) if (badifu eq '') then badifu=harnames[i] $ else badifu+='-'+harnames[i] drp2qual=drp2qual OR sdss_flagval('MANGA_DRP2QUAL','BADIFU') endif else begin ; If the entire IFU wasn't unplugged, also look at ; throughput variations from the flatfield ; Do a sigma-clipped mean of the throughput values temp=ml_meanclip(fflat[2000,thesefibers],themean,thesig) ; If mean throughput is < 70%, or sigma> 15% then call the IFU bad if ((themean lt 0.7)or(thesig gt 0.15)) then begin splog, strcompress('WARNING! Bad IFU fibers in harness '+harnames[i]) if (badifu eq '') then badifu=harnames[i] $ else badifu+='-'+harnames[i] drp2qual=drp2qual OR sdss_flagval('MANGA_DRP2QUAL','BADIFU') endif endelse endfor ; Check the dither position information ; If MGDPOS doesn't match MGDRA and MGDDEC there's a problem mgdpos=strtrim(sxpar(objhdr,'MGDPOS'),2) mgdra=float(sxpar(objhdr,'MGDRA')) mgddec=float(sxpar(objhdr,'MGDDEC')) result=mdrp_checkdither(mgdpos,mgdra,mgddec,platescale=60.) if (result eq 1) then drp2qual=drp2qual OR sdss_flagval('MANGA_DRP2QUAL','BADDITHER') ; Add DRP2QUAL and MISSIFU to the FITS header sxaddpar, objhdr, 'DRP2QUAL', drp2qual, 'DRP-2d quality bitmask' sxaddpar, objhdr, 'BADIFU', badifu, 'Harness names of missing/bad ifus' ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Write extracted, lambda-calibrated spectra to disk following the MaNGA data model ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; obj0hdr=objhdr sxaddpar, objhdr, 'PCOUNT', 0, after='NAXIS2' sxaddpar, objhdr, 'GCOUNT', 1, after='PCOUNT' ml_mwrfits,dummyext,outname,hdr=obj0hdr, /create ; Blank ext 0 with full header ml_mwrfits, flux, outname, hdr=objhdr, extname='FLUX' ;extracted flux ml_mwrfits, fluxivar, outname, extname='IVAR' ; extracted inverse variance ml_mwrfits, mangamask, outname, extname='MASK' ;pixel mask ml_mwrfits, vacset, outname, extname='WSET' ;trace-set for wavelength sol ml_mwrfits, skydispset, outname, extname='DISPSET' ;trace-set for sky dispersion sol ml_mwrfits, slitmap, outname, extname='SLITMAP' ;slitmap ml_mwrfits, xnow, outname, extname='XPOS' ;x pos of traces on CCD ml_mwrfits, superfit, outname, extname='SUPERFLAT' ;superflat vector from quartz lamps spawn, ['gzip', '-f', outname], /noshell if ~keyword_set(visual) then begin dfpsclose ; Convert to PDF and remove old PS file spawn, strcompress('ps2pdf '+plotname+' '+ml_strreplace(plotname,'.ps','.pdf')) spawn, strcompress('rm -f '+plotname) !p.multi=0 !p=origp !x=origx !y=origy endif heap_gc obj_destroy,configuration return end