;+ ; NAME: ; mdrp_fluxcal ; ; PURPOSE: ; Flux calibrate sky-subtracted MaNGA frames (mgSFrame) on a ; per-camera basis. This routine needs to be run once ; per spectrograph; innames should be ; a 2-element (N_camera=2) array of names pointing to the ; corresponding mgSFrame files. ; ; b1calib is essentially like extension 0 of the ; BOSS files spFluxcalib, which is calibimg = float( bspline_valu(loglam1, thisset, x2=x2) ) ; CALLING SEQUENCE: ; ... ; ; INPUTS: ; innames - filenames for one or many exposures of a plate in one of the ; spectrographs. ; This routine needs to be run per spectrograph. innames should be ; a N_camera (2) element array of names pointing to the ; corresponding mgSFrame files. N_camera=2 ; ; OPTIONAL INPUTS: ; plotname - File for the output plot (overrides default) ; /visual - Plots go to X-display instead of file ; ; OPTIONAL OUTPUTS: ; outnames - Array of names of the output mgFFrame files ; bcalib - ; rcalib - ; ; KEYWORDS ; nowritefile - when specified, it does not output the stellar model, mratio, ; mrativar, etc. ; for each star in a file named fluxcal-sp[1,2]-expno.fits ; ; PROCEDURES CALLED: ; ... ; ; REVISION HISTORY: ; 25-Mar-2014 First written based on BOSS modules (Renbin Yan) ; 06-May-2014 Tweaked to write out mgFFrame files (David Law) ; 15-May-2014 Added quality control bitmask logic (David Law) ; 08-Sep-2014 Nexp != 1 will break QA routines, make explicit that input ; must be on a per-exposure basis (D. Law) ; 11-Sep-2014 Migrated drpqual to drp2qual (D. Law) ; 04-Apr-2016 Added option for Munari models by Yanping Chen ; 04-Sep-2017 Changed smoothing scale in bspline, updated Telluric and ; stellar feature definition, updated meanthrupt curve and ; how it is used. (Renbin Yan) ; 12-Jan-2018 Add calculation of normalized throughputs (R. Yan, D. Law) ; 03-Jul-2018 Pass through super-sky model (D. Law, STScI) ; 19-Dec-2019 Revise handling of LSF pass-through (D. Law) ;- ;------------------------------------------------------------------------- function spflux_masklines, loglam, hwidth=hwidth, stellar=stellar, $ telluric=telluric if (NOT keyword_set(hwidth)) then $ hwidth = 5.7e-4 ; Default is to mask +/- 5.7 pix = 400 km/sec mask = bytarr(size(loglam,/dimens)) if (keyword_set(stellar)) then begin starwave = [ $ ; 3692.6 , $ ; H-16 ; 3698.2 , $ ; H-15 3704.9 , $ ; H-14 ; 3707.1 , $ ; Ca-II 3713.0 , $ ; H-13 3723.0 , $ ; H-12 3735.4 , $ ; H-11 ; 3738.0 , $ ; Ca-II 3751.2 , $ ; H-10 3771.7 , $ ; H-9 3799.0 , $ ; H-8 3836.5 , $ ; H-7 3890.2 , $ ; H-6 3934.8 , $ ; Ca_k 3969.6 , $ ; Ca_H 3971.2 , $ ; H-5 4102.9 , $ ; H-delta 4300. , $ ; G-band 4305. , $ ; G-band 4310. , $ ; more G-band 4341.7 , $ ; H-gamma 4862.7 , $ ; H-beta ; 4687.1 , $ ; He II 5168.8 , $ ; Mg I 5174.1 , $ ; Mg I 5185.0 , $ ; Mg I 5891.6 , $ ; Na I 5897.6 , $ ; Na I 6564.6 , $ ; H-alpha 8500.4 , $ ; Ca II 8544.4 , $ ; Ca II 8664.5 , $ ; Ca II 8752.9 , $ ; H I 8865.3 , $ ; H I ; RY: commented out the following three lines as these are in telluric absorption regions. ; If we mask them, we get bad bspline interpolation in these regions. ; They are not useful for stellar fitting anyway as they are in telluric regions. ; 9017.8 , $ ; H I ; 9232.2 , $ ; H I Pa-6 ; 9548.8 , $ ; H I Pa-5 10052.6 ] ; H I (Pa-delta) ; airtovac, starwave ; commented out because wavelengths of features have already been set in vacuum.RY Jul 13, 2015 for i=0L, n_elements(starwave)-1 do begin mask = mask OR (loglam GT alog10(starwave[i])-hwidth $ AND loglam LT alog10(starwave[i])+hwidth) endfor endif if (keyword_set(telluric)) then begin tellwave1 = [6842., 7146., 7588., 8105., 8910.] tellwave2 = [6980., 7390., 7730., 8440., 9880.] for i=0L, n_elements(tellwave1)-1 do begin mask = mask OR (loglam GT alog10(tellwave1[i]) $ AND loglam LT alog10(tellwave2[i])) endfor endif return, mask end ;------------------------------------------------------------------------------ ; Divide the spectrum by a median-filtered spectrum. ; The median-filtered version is computed ignoring stellar absorp. features. function spflux_medianfilt, loglam, objflux, objivar, width=width, $ newivar=newivar, _EXTRA=KeywordsForMedian ndim = size(objflux, /n_dimen) dims = size(objflux, /dimens) npix = dims[0] if (ndim EQ 1) then nspec = 1 $ else nspec = dims[1] ;---------- ; Loop over each spectrum medflux = 0 * objflux if (arg_present(objivar)) then newivar = 0 * objivar for ispec=0L, nspec-1 do begin ; For the median-filter, ignore points near stellar absorp. features, ; but keep points near telluric bands. qgood = 1 - spflux_masklines(loglam[*,ispec], /stellar,hwidth=8.e-4) ; Median-filter, but skipping masked points igood = where(qgood, ngood) thisback = fltarr(dims[0]) if (ngood GT 1) then begin thisback[igood] = djs_median(objflux[igood,ispec], width=width, $ _EXTRA=KeywordsForMedian) endif thisback = djs_maskinterp(thisback, (qgood EQ 0), /const) ; Force the ends of the background to be the same as the spectrum, ; which will force the ratio of the two to be unity. hwidth = ceil((width-1)/2.) thisback[0:hwidth] = objflux[0:hwidth,ispec] thisback[npix-1-hwidth:npix-1] = objflux[npix-1-hwidth:npix-1,ispec] czero2 = where(thisback eq 0., count2) if count2 gt 0 then thisback[czero2] = 1. medflux[*,ispec] = objflux[*,ispec] / thisback if (arg_present(objivar)) then $ newivar[*,ispec] = objivar[*,ispec] * thisback^2 endfor return, medflux end ;------------------------------------------------------------------------------ function spflux_bestmodel, loglam, objflux, objivar, dispimg, kindx=kindx1, $ plottitle=plottitle,template=template compile_opt defint32 filtsz = 99 ; the size of the window used in median-filter the spectra. cspeed = 2.99792458e5 ndim = size(objflux, /n_dimen) dims = size(objflux, /dimens) npix = dims[0] if (ndim EQ 1) then nspec = 1 $ else nspec = dims[1] ;---------- ; Median-filter the object fluxes medflux = spflux_medianfilt(loglam, objflux, objivar, $ width=filtsz, /reflect, newivar=medivar) sqivar = sqrt(medivar) ;---------- ; Mask out the telluric bands sqivar = sqivar * (1 - spflux_masklines(loglam, /telluric)) ;---------- ; Load the Kurucz models into memory ;junk = spflux_read_kurucz(kindx=kindx) CASE template OF 'kurucz': junk = spflux_read_kurucz(kindx=kindx,dslgpsize=dslgpsize) ;;Yanping test 'munari': junk = spflux_read_munari(kindx=kindx,dslgpsize=dslgpsize) ;;Yanping added 'BOSZ' : junk = spflux_read_bosz(kindx=kindx,dslgpsize=dslgpsize) ;;Yanping added else : message,"Flux calibration templates has to be specified and be one of the three: 'kurucz','munari', 'BOSZ'." ENDCASE nmodel = n_elements(kindx) ;---------- ; Fit the redshift just by using a canonical model ifud = where(kindx.teff EQ 6000 AND kindx.g EQ 4 AND kindx.feh EQ -1.5) if (ifud[0] EQ -1) then $ message, 'Could not find fiducial model!' nshift = ceil(1000./cspeed/alog(10.)/dslgpsize/2)*2 ; set this to cover +/-500 km/s ;logshift = (-nshift/2. + findgen(nshift)) * 1.d-4 logshift = (-nshift/2. + findgen(nshift)) * dslgpsize ;;Yanping test chivec = fltarr(nshift) for ishift=0L, nshift-1 do begin ;modflux = spflux_read_kurucz(loglam-logshift[ishift], $ ; dispimg, iselect=ifud) CASE template OF 'kurucz': modflux = spflux_read_kurucz(loglam-logshift[ishift], dispimg, iselect=ifud) ;;Yanping edited 'munari': modflux = spflux_read_munari(loglam-logshift[ishift], dispimg, iselect=ifud) ;;Yanping edited 'BOSZ' : modflux = spflux_read_bosz(loglam-logshift[ishift], dispimg, iselect=ifud) else : message,"Flux calibration templates has to be specified and be one of the three: 'kurucz','munari', 'BOSZ'." ENDCASE ; Median-filter this model medmodel = spflux_medianfilt(loglam, modflux, $ width=filtsz, /reflect) for ispec=0L, nspec-1 do begin chivec[ishift] = chivec[ishift] + computechi2(medflux[*,ispec], $ sqivar[*,ispec], medmodel[*,ispec]) endfor endfor zshift = (10.d^logshift - 1) ; Convert log-lambda shift to redshift zpeak = find_nminima(chivec, zshift, errcode=errcode) splog, 'Best-fit velocity for std star = ', zpeak * cspeed, ' km/s' if (errcode NE 0) then $ splog, 'Warning: Error code ', errcode, ' fitting std star' ; Warning messages if total(finite(chivec) eq 0) gt 0 then begin if total(medivar lt 0) gt 0 then begin message,'There are negative ivar values causing chi-square to be NaN or the likes.' endif else message,'chi-square are NaN or the likes, but not caused by negative ivar.' endif ;---------- ; Generate the Kurucz models at the specified wavelengths + dispersions, ; using the best-fit redshift ;modflux = spflux_read_kurucz(loglam-alog10(1.+zpeak), dispimg) CASE template OF 'kurucz': modflux = spflux_read_kurucz(loglam-alog10(1.+zpeak), dispimg) ;Yanping edited 'munari': modflux = spflux_read_munari(loglam-alog10(1.+zpeak), dispimg) ;Yanping edited 'BOSZ' : modflux = spflux_read_bosz(loglam-alog10(1.+zpeak), dispimg) ;Yanping edited else : message,"Flux calibration templates has to be specified and be one of the three: 'kurucz','munari', 'BOSZ'." ENDCASE ; Need to redo median-filter for the data with the correct redshift medflux = spflux_medianfilt(loglam-alog10(1.+zpeak), objflux, objivar, $ width=filtsz, /reflect, newivar=medivar) sqivar = sqrt(medivar) ;---------- ; Mask out the telluric bands sqivar = sqivar * (1 - spflux_masklines(loglam, /telluric)) ;---------- ; Loop through each model, computing the best chi^2 ; as the sum of the best-fit chi^2 to each of the several spectra ; for this same object. Counting only the regions around stellar absorption features. ; We do this after a median-filtering of both the spectra + the models. chiarr = fltarr(nmodel,nspec) chivec = fltarr(nmodel) medmodelarr = modflux*0.0 mlines = spflux_masklines(loglam-alog10(1.+zpeak), hwidth=12e-4, /stellar) linesqivar = sqivar * mlines linechiarr = fltarr(nmodel,nspec) linechivec = fltarr(nmodel) for imodel=0L, nmodel-1 do begin ; Median-filter this model medmodelarr[*,*,imodel] = spflux_medianfilt(loglam-alog10(1.+zpeak), modflux[*,*,imodel], $ width=filtsz, /reflect) for ispec=0L, nspec-1 do begin chiarr[imodel,ispec] = total((medflux[*,ispec]-medmodelarr[*,ispec,imodel])^2*sqivar[*,ispec]*sqivar[*,ispec]) linechiarr[imodel,ispec] = total((medflux[*,ispec]-medmodelarr[*,ispec,imodel])^2*linesqivar[*,ispec]*linesqivar[*,ispec]) endfor chivec[imodel] = total(chiarr[imodel,*]) linechivec[imodel] = total(linechiarr[imodel,*]) endfor ;---------- ; Return the best-fit model ; Computed both full spectra chi^2 and line chi^2, but use the line chi^2 to select the best model. linechi2 = min(linechivec, ibest) linedof = total(linesqivar NE 0) ; splog, 'Best-fit total chi2/DOF = ', minchi2/(dof>1) splog, 'Best-fit line chi2/DOF = ', linechi2/(linedof>1) bestflux = modflux[*,*,ibest] medmodel = spflux_medianfilt(loglam-alog10(1.+zpeak), bestflux, $ width=filtsz, /reflect) ;---------- ; Compute the median S/N for all the spectra of this object, ; and for those data just near the absorp. lines indx = where(objivar gt 0,ct) if (ct GT 1) then sn_median = median(objflux[indx] * sqrt(objivar[indx])) else sn_median=0 indx = where(mlines, ct) if (ct GT 1) then $ linesn_median = median(objflux[indx] * sqrt(objivar[indx])) $ else $ linesn_median = 0. splog, 'Full median S/N = ', sn_median splog, 'Line median S/N = ', linesn_median kindx1 = create_struct(kindx[ibest], $ 'IMODEL', ibest, $ 'Z', zpeak, $ 'SN_MEDIAN', float(sn_median), $ ; 'CHI2', minchi2, $ ; 'DOF', dof, $ 'LINESN_MEDIAN', linesn_median, $ 'LINECHI2', linechi2, $ 'LINEDOF', linedof) ;---------- ; Plot the filtered object spectrum, overplotting the best-fit Kurucz/Munari model ;;Yanping edited ; Select the observation to plot that has the highest S/N, ; and one that goes blueward of 4000 Ang. snvec = total(objflux * sqrt(objivar), 1) $ * (10.^loglam[0,*] LT 4000 OR 10.^loglam[npix-1,*] LT 4000) junk = max(snvec, iplot) ; Best blue exposure snvec = total(objflux * sqrt(objivar), 1) $ * (10.^loglam[0,*] GT 8600 OR 10.^loglam[npix-1,*] GT 8600) junk = max(snvec, jplot) ; Best red exposure csize = 0.85 djs_plot, [3840., 4120.], [0.0, 1.4], /xstyle, /ystyle, /nodata, $ xtitle='Wavelength [Ang]', ytitle='Normalized Flux', $ title=plottitle if (iplot[0] NE -1) then begin djs_oplot, 10^loglam[*,iplot], medflux[*,iplot] djs_oplot, 10^loglam[*,iplot], medmodel[*,iplot], color='red' endif xyouts, 3860, 1.25, kindx1.model, charsize=csize ; djs_xyouts, 4000, 0.3, charsize=csize, $ ; string(minchi2/(dof>1), format='("Total \chi^2/DOF=",f5.2)') djs_xyouts, 4000, 0.2, charsize=csize, $ string(linechi2/(linedof>1), format='("Lines \chi^2/DOF=",f5.2)') djs_xyouts, 3860, 0.1, string(kindx1.feh, kindx1.teff, kindx1.g, $ zpeak*cspeed, $ format='("Fe/H=", f4.1, " T_{eff}=", f6.0, " g=", f3.1, " cz=",f5.0)'), $ charsize=csize djs_plot, [8440., 9160.], [0.0, 1.4], /xstyle, /ystyle, /nodata, $ xtitle='Wavelength [Ang]', ytitle='Normalized Flux' if (jplot[0] NE -1) then begin djs_oplot, 10^loglam[*,jplot], medflux[*,jplot] djs_oplot, 10^loglam[*,jplot], medmodel[*,jplot], color='red' endif ; outmodflux = spflux_read_kurucz(outloglam-alog10(1.+zpeak), outdispimg,iselect=ibest) return,bestflux end ;------------------------------------------------------------------------------ function spflux_goodfiber, pixmask qgood = ((pixmask AND drp2pixmask_bits('NOPLUG')) EQ 0) $ AND ((pixmask AND drp2pixmask_bits('BADTRACE')) EQ 0) $ AND ((pixmask AND drp2pixmask_bits('BADFLAT')) EQ 0) $ AND ((pixmask AND drp2pixmask_bits('BADARC')) EQ 0) $ AND ((pixmask AND drp2pixmask_bits('MANYBADCOLUMNS')) EQ 0) $ AND ((pixmask AND drp2pixmask_bits('NEARWHOPPER')) EQ 0) $ AND ((pixmask AND drp2pixmask_bits('MANYREJECTED')) EQ 0) return, qgood end ;------------------------------------------------------------------------------ function spflux_bspline, loglam, mratio, mrativar, inmask=inmask, outmask=outmask, $ everyn=everyn, disp=disp, hwidth=hwidth, stellar=stellar, telluric=telluric isort = bsort(loglam) nord = 3 if NOT keyword_set(hwidth) then hwidth=12.e-4 if n_elements(inmask) eq 0 then begin ; Choose the break points using the EVERYN option, but masking ; out more pixels near stellar features and/or telluric just when selecting them. mask1 = 1 - spflux_masklines(loglam, hwidth=hwidth, stellar=stellar,telluric=telluric) endif else mask1 = 1 - inmask ii = where(mrativar[isort] GT 0 AND mask1[isort] EQ 1) bkpt = 0 fullbkpt = bspline_bkpts(loglam[isort[ii]], everyn=everyn, $ bkpt=bkpt, nord=nord) outmask1 = 0 if (keyword_set(disp)) then begin x2 = disp[isort] endif sset = bspline_iterfit(loglam[isort], mratio[isort], $ invvar=mrativar[isort], lower=3, upper=3, fullbkpt=fullbkpt, $ maxrej=ceil(0.05*n_elements(indx)), outmask=outmask1, nord=nord, $ x2=x2, npoly=2*keyword_set(disp), requiren=(everyn-1)>1) if (max(sset.coeff) EQ 0) then $ message, 'B-spline fit failed!!' if (arg_present(outmask)) then begin outmask = bytarr(size(loglam,/dimens)) outmask[isort] = outmask1 endif return, sset end ;------------------------------------------------------------------------------ function typingmodule, objflux,loglam,objivar,dispimg,sfd_ebv,psfmag,kindx=kindx,modflux=modflux,plottitle=plottitle,$ template=template, thekfile=thekfile, targetflag=targetflag npix = (size(objflux,/dimen))[0] ;---------- ; For each star, find the best-fit model. fm_unred,10.^loglam,loglam*0+1,1.0,unreddenfactor ; !p.multi = [0,1,2] modflux = 0 * objflux ; Find the best-fit model -- evaluated for each exposure [NPIX,NEXP] thismodel = spflux_bestmodel(loglam, objflux*unreddenfactor^sfd_ebv, objivar/unreddenfactor^(2*sfd_ebv),dispimg, $ kindx=kindx,plottitle=plottitle,template=template) ; CASE template OF ; 'kurucz': thismodel = spflux_bestmodel(loglam, objflux, $ ; objivar, dispimg, kindx=kindx, plottitle=plottitle,/kurucz) ; if keyword_set(munari) then thismodel = spflux_bestmodel(loglam, objflux, $ ; objivar, dispimg, kindx=kindx, plottitle=plottitle,/munari) ; Also evaluate this model over a big wavelength range [3000,11000] Ang. tmploglam = 3.4771d0 + lindgen(5644) * 1.d-4 tmpdispimg = 0 * tmploglam + 1.0 ; initializing this resolution vector bluedispimg = tmpdispimg & reddispimg = tmpdispimg bside = where(tmploglam lt alog10(6300.)) rside = where(tmploglam gt alog10(5900.)) middle = where(tmploglam lt alog10(6300.) and tmploglam gt alog10(5900.)) bluedispimg[bside] = interpol(dispimg[*,0],loglam[*,0],tmploglam[bside]) reddispimg[rside] = interpol(dispimg[*,1],loglam[*,1],tmploglam[rside]) tmpdispimg[bside] = bluedispimg[bside] tmpdispimg[rside] = reddispimg[rside] tmpdispimg[middle] = (bluedispimg[middle]+reddispimg[middle])/2. ;tmpflux = spflux_read_kurucz(tmploglam-alog10(1+kindx.z), tmpdispimg, $ ; iselect=kindx.imodel) CASE template OF 'kurucz': tmpflux = spflux_read_kurucz(tmploglam-alog10(1+kindx.z), tmpdispimg, $ iselect=kindx.imodel, thekfile=thekfile) 'munari': tmpflux = spflux_read_munari(tmploglam-alog10(1+kindx.z), tmpdispimg, $ iselect=kindx.imodel, thekfile=thekfile) 'BOSZ': tmpflux = spflux_read_bosz(tmploglam-alog10(1+kindx.z), tmpdispimg, $ iselect=kindx.imodel, thekfile=thekfile) else : message,'Template is not specified correctly.' ENDCASE ; The returned models are redshifted, but not fluxed or ; reddened. Do that now... we compare data vs. model reddened. ; extcurve1 = ext_odonnell(10.^loglam, 3.1) ; extinct,10.^loglam,extcurve1,/ccm,Rv=3.1 ; extcurve1 is A_lambda for Av=1 ; thismodel = thismodel * 10.^(-extcurve1 * 3.1 * sfd_ebv / 2.5) fm_unred,10.^loglam,loglam*0+1,-1.0,reddenfactor ; fluxreddened contain the reddening vector for E(B-V)=1.0 thismodel = thismodel*reddenfactor^sfd_ebv ; extcurve2 = ext_odonnell(10.^tmploglam, 3.1) ; extinct,10.^tmploglam,extcurve2,/ccm,Rv=3.1 ; tmpflux = tmpflux * 10.^(-extcurve2 * 3.1 * sfd_ebv / 2.5) fm_unred,10.^tmploglam,tmploglam*0+1,-1.0,reddenfactor2 tmpflux = tmpflux * reddenfactor2^sfd_ebv ; Now integrate the apparent magnitude for this spectrum, ; The units of FTHRU are such that m = -2.5*alog10(FTHRU) + (48.6-2.5*17) ; Note that these computed magnitudes, THISMAG, should be equivalent ; to THISINDX.MAG in the case of no reddening. wavevec = 10.d0^tmploglam flambda2fnu = wavevec^2 / 2.99792e18 photometry = 'sdss' ; Both APASS and SDSS are using sdss photometry system if (targetflag and sdss_flagval('MANGA_TARGET2',['STELLIB_PS1','STD_PS1_COM'])) ne 0 then photometry='ps1' if (targetflag and sdss_flagval('MANGA_TARGET2','STELLIB_GAIA')) ne 0 then photometry='gaiadr1' if (targetflag and sdss_flagval('MANGA_TARGET2','STELLIB_GAIADR2')) ne 0 then photometry='gaiadr2' CASE photometry OF 'ps1' : begin fthru = filter_thru(tmpflux * flambda2fnu, waveimg=wavevec, /toair, filternames=['ps1_g.txt','ps1_r.txt','ps1_i.txt','ps1_z.txt']) thismag = -2.5*alog10(fthru) - (48.6-2.5*17) thismag = [[-999.],[thismag]] end 'gaiadr1': begin ; Using the passbands from Gaia DR2 to compute the mag, which is not quite appropriate for DR1. Better to avoid using gaiadr1 fthru = filter_thru(tmpflux * flambda2fnu, waveimg=wavevec, /toair,filternames=['gaia_G.dat']) thismag = -2.5 * alog10(fthru) - (48.6-2.5*17) end 'gaiadr2': begin fthru = filter_thru(tmpflux * flambda2fnu, waveimg=wavevec, /toair, filternames=['gaia_G.dat','gaia_BP.dat','gaia_RP.dat']) thismag = -2.5*alog10(fthru) - (48.6-2.5*17) ; Convert from AB system using revised passband to Vega system in the Gaia DR2 as-released system. The zeropoints are from Evans et al. (2018) thismag = thismag + [[25.6884-25.7916],[25.3514-25.3862],[24.7619-25.1162]] thismag = [[-999.],[thismag],[-999.]] end 'sdss' : begin ; this applies to both SDSS and APASS standards. fthru = filter_thru(tmpflux * flambda2fnu, waveimg=wavevec, /toair) thismag = -2.5 * alog10(fthru) - (48.6-2.5*17) end else : ENDCASE ; !!!!!! IMPORTANT !!!!!!!!! ; !!!!!! THIS IS NOT THE ONLY PLACE WHERE THE MODEL SPECTRUM IS NORMALIZED.!!!! ; !!!!!! For MaStar plates, we renormalize them again after putting in ; individual extinction.!!!!! ; Compute SCALEFAC = (plugmap flux) / (flux of the model spectrum) CASE photometry OF 'gaiadr1': begin scalefac = 10.^((thismag-psfmag[1])/2.5) kindx.mag[1] = psfmag[1] end 'gaiadr2': begin scalefac = 10.^((thismag[1]-psfmag[1])/2.5) kindx.mag[1:3] = reform(thismag[1:3])+psfmag[1]-thismag[1] end 'ps1' : begin scalefac = 10.^((thismag[2]-psfmag[2])/2.5) kindx.mag[1:4] = reform(thismag[1:4])+psfmag[2]-thismag[2] end 'sdss' : begin scalefac = 10.^((thismag[2]-psfmag[2])/2.5) kindx.mag = reform(thismag)+psfmag[2]-thismag[2] end ENDCASE thismodel = thismodel * scalefac modflux = thismodel splog, prelog='' !p.multi = 0 return,0 end ;------------------------------------------------------------------------------ function matchbundleratio2, bundleratio, ratioivar, reffiber, distarr, cenwave=cenwave,dwave=dwave,predictratio=predictratio,usemask=usemask,weight=weight,ivarweight=ivarweight,lambda=lambda ;common com_fluxbundle, cube,xarr,yarr,wavearr common com_fluxbundle2, psfplane,xarr,wavearr if n_elements(cenwave) ne n_elements(dwave) then message,'Error: cenwave and dwave must have the same number of elements.' if n_elements(weight) ne n_elements(lambda) then message,'Error: lambda and weight must have the same number of elements.' if (size(weight,/dimen))[0] ne (size(ivarweight,/dimen))[0] then message,"Error: ivarweight's first dimension must match weight." ; stop xpos = interpol(findgen(n_elements(xarr)),xarr,distarr) nfiber = (size(distarr,/dimen))[0] ndistwave = (size(distarr,/dimen))[1] if nfiber ne (size(ivarweight,/dimen))[2] then message,"Error: ivarweight's third dimension must match nfiber." nwave = n_elements(cenwave) n_wavearr = n_elements(wavearr) if ndistwave ne n_wavearr then message,'Error: the 2nd dimension of distarrr does not match the dimension of wavearr.' siz = size(bundleratio,/dimen) nfile = siz[1] if siz[2] ne nfiber then message,'Error: the 3rd dimension of bundleratio does not match nfiber.' if siz[0] ne nwave then message,'Error: the 1st dimension of bundleratio does not match nwave.' covfn=fltarr(nwave,nfile,nfiber) for i=0,nfiber-1 do begin covfn_wavearr = interpolate(psfplane,xpos[i,*],indgen(n_wavearr)) covfn_allwave = interpol(covfn_wavearr,wavearr,lambda[*,*,i]) weightedflux = covfn_allwave*weight[*,*,i]*ivarweight[*,*,i] for ifile=0,nfile-1 do begin for j=0,nwave-1 do begin lam1=cenwave[j]-dwave[j]/2. lam2=cenwave[j]+dwave[j]/2. ind = where(lambda[*,ifile,i] ge lam1 and lambda[*,ifile,i] lt lam2 and ivarweight[*,ifile,i] gt 0, ngood) if ngood eq 0 then continue covfn[j,ifile,i] = total(weightedflux[ind,ifile])/total(ivarweight[ind,ifile,i]) endfor endfor endfor refpredict = rebin(covfn[*,*,reffiber],nwave,nfile,nfiber) predictratio=covfn/(refpredict + (refpredict eq 0)) useid = where(indgen(nfiber) ne reffiber and usemask,nuse) diff = (bundleratio[*,*,useid]-predictratio[*,*,useid])*sqrt(ratioivar[*,*,useid] < (1/0.003)^2) ind_valid = where(ratioivar[*,*,useid] gt 0,nvalid) if nvalid gt 0 then chi2 = total(diff[ind_valid]*diff[ind_valid])/(nvalid-1) else chi2 = 1.e38 return,chi2 end function mcmcbundleratio,fluxratio,fluxratioivar,offsetx0=offsetx0,offsety0=offsety0,fiberoffx=fiberoffx,fiberoffy=fiberoffy,nfiber=nfiber,n_wavearr=n_wavearr,reffiber=reffiber,centerwave=centerwave,dwave=dwave,innermask=innermask,weight=weight,ivarweight=ivarweight,lambda=lambda,scale=scale,rotation=rotation,xbest=xbest,ybest=ybest,nstep=nstep,debug=debug,errorreport=errorreport,fixlist=fixlist,seed=seed if n_elements(offsetx0) ne n_elements(offsety0) then message,'Error: Offsetx0 and Offsety0 do NOT have the same number of elements.' if n_elements(offsetx0) ne n_wavearr then message,'Error: Number of elements in offsetx0 (offsety0) does not equal to n_wavearr.' if NOT keyword_set(fixlist) then fixlist = [0,0,0,0] stepsize = [0.06,2,0.08] ; [positional move radius, rotation in degree, scale of ADR] if NOT keyword_set(nstep) then nstep =1000L ;splog,'TEST: ',randomu(seed,1) random=randomu(seed,nstep*4) rr = sqrt(random[0:nstep-1]) theta = random[nstep:2*nstep-1]*!pi*2 xstep = rr*cos(theta) ystep = rr*sin(theta) if fixlist[2] ne 1 then anglestep = random[2*nstep:3*nstep-1]-0.5 else anglestep = 0.0 ;*stepsize[1]*!dtor if fixlist[3] ne 1 then scalestep = random[3*nstep:4*nstep-1]-0.5 else scalestep = 0.0 ;*stepsize[2]) offsetx = offsetx0 offsety = offsety0 dxarr = fiberoffx#(fltarr(n_wavearr)+1)-(fltarr(nfiber)+1)#offsetx dyarr = fiberoffy#(fltarr(n_wavearr)+1)-(fltarr(nfiber)+1)#offsety distarr = sqrt(dxarr*dxarr+dyarr*dyarr) oldchi2 = matchbundleratio2(fluxratio,fluxratioivar,reffiber,distarr,cenwave=centerwave,dwave=dwave,usemask=innermask,weight=weight,ivarweight=ivarweight,lambda=lambda) xposarr = fltarr(nstep+1) yposarr = fltarr(nstep+1) scalearr = fltarr(nstep+1) rotationarr = fltarr(nstep+1) chi2arr = fltarr(nstep+1) stepnum = intarr(nstep+1) chi2arr[0] = oldchi2 newscale = 1.0 newrotation = 0.0 newx = fiberoffx newy = fiberoffy npoint=0 badstep_ct=0 bestcounter=0 bestchi2 = oldchi2 for i=0,nstep-1 do begin ; if i mod 1000 eq 0 then print,i if fixlist[3] ne 1 then begin tmpscale = (newscale*10^(scalestep[i]*stepsize[2]) < 1.2) > 0.8 offsetx=offsetx0*tmpscale offsety=offsety0*tmpscale endif else tmpscale = 1.0 if fixlist[2] ne 1 then begin tmprotation = (newrotation + anglestep[i]*stepsize[1]*!dtor < 0.17 ) > (-0.17) rotmatrix = [[cos(tmprotation),-sin(tmprotation)],[sin(tmprotation),cos(tmprotation)]] xy2= [[offsetx],[offsety]]#rotmatrix offsetx = xy2[*,0] offsety = xy2[*,1] endif else tmprotation=0.0 if fixlist[0] ne 1 then dx = xstep[i]*stepsize[0] else dx=0.0 if fixlist[1] ne 1 then dy = ystep[i]*stepsize[0] else dy=0.0 dxarr = (newx+dx)#(fltarr(n_wavearr)+1)-(fltarr(nfiber)+1)#offsetx dyarr = (newy+dy)#(fltarr(n_wavearr)+1)-(fltarr(nfiber)+1)#offsety distarr = sqrt(dxarr*dxarr+dyarr*dyarr) newchi2 = matchbundleratio2(fluxratio,fluxratioivar,reffiber,distarr,cenwave=centerwave,dwave=dwave,usemask=innermask,weight=weight,ivarweight=ivarweight,lambda=lambda) if newchi2 lt oldchi2 or randomu(seed,1) lt exp(oldchi2-newchi2) then begin npoint += 1 xposarr[npoint] = xposarr[npoint-1]+dx yposarr[npoint] = yposarr[npoint-1]+dy scalearr[npoint] = tmpscale rotationarr[npoint] = tmprotation chi2arr[npoint] = newchi2 oldchi2 = newchi2 newx = newx+dx newy = newy+dy newscale = tmpscale newrotation = tmprotation if newchi2 lt bestchi2*0.95 then begin oldbestchi2 = bestchi2 bestchi2 = newchi2 bestcounter = 0 endif else bestcounter += 1 stepnum[npoint]=i endif else begin badstep_ct +=1 if badstep_ct gt 300 then begin stepsize = stepsize*0.6 if stepsize[0] lt 0.003 and stepsize[1] lt 0.25 and stepsize[2] lt 0.004 then break badstep_ct=0 endif endelse if bestcounter gt 300 then break endfor ; print,'Iteration stopped on step ',i ; print,'Number of bad steps:',badstep_ct chi2arr = chi2arr[0:npoint] xposarr = xposarr[0:npoint] yposarr = yposarr[0:npoint] scalearr = scalearr[0:npoint] rotationarr = rotationarr[0:npoint] tmp = min(chi2arr,ind) xbest = xposarr[ind] ybest = yposarr[ind] scale = scalearr[ind] rotation = rotationarr[ind] fchi2 = chi2arr[ind] if keyword_set(errorreport) then begin q=where(chi2arr lt 2*min(chi2arr),nq) if nq gt 2 then begin xerr = stddev(xposarr[q]) yerr = stddev(yposarr[q]) print,'MCMC X uncertainty:',xerr print,'MCMC Y uncertainty:',yerr endif endif ; erase ; multiplot,[0,1,5,0,0] ; plot,stepnum,xposarr,ytitle='X POS' ; multiplot ; plot,stepnum,yposarr,ytitle='Y POS' ; multiplot ; plot,stepnum,scalearr,ytitle='Scale' ; multiplot ; plot,stepnum,rotationarr,ytitle='Rotation' ; multiplot ; plot,stepnum,chi2arr,ytitle='Chi Square' ; multiplot,[1,1],/default if keyword_set(debug) then stop return,fchi2 end ; innames should be an nexp x ncamera array of names ; pointing to the corresponding mgSFrame files. pro mdrp_fluxcal, innames, outnames=outnames, bcalib=bcalib,rcalib=rcalib,docams=docams, visual=visual, plotname=plotname, centerwave=centerwave, dwave=dwave,staroffx=staroffx,staroffy=staroffy, combinedir=combinedir,nowritefile=nowritefile,writename=writename, debugpsf=debugpsf,template=template,time=time common com_fluxbundle2, psfplane,xarr,wavearr common plotcolors, black, red, green, blue, cyan, magenta, yellow, white,grey time0 = systime(1) cc = 2.997925e10 ; in cm planck = 6.626d-27 ; planck's constant in cgs mirrorarea = 3.5814e4 ; in cm^2 ; diameter 2.5 meter with 1.3 meter obstruction (secondary+baffle) manga_redux_dir = getenv('MANGA_SPECTRO_REDUX') drpver = ~keyword_set(drpver) ? mangadrp_version(/simple) : strtrim(drpver,2) if NOT keyword_set(template) then template = 'BOSZ' if template ne 'munari' and template ne 'kurucz' and template ne 'BOSZ' then message,'Template keyword need to be either "munari", "kurucz", or "BOSZ".' ;;Yanping change 'or' to 'and' ; Set the output names outnames=strarr(n_elements(innames)) for i=0,n_elements(innames)-1 do begin outnames[i]=ml_strreplace(innames[i],'mgSFrame','mgFFrame') ; Make sure outname name doesn't end in '.gz' outnames[i]=ml_ensurenogz(outnames[i]) endfor nfile = n_elements(innames) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Quality control checks. If fails return from function ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Require that nfile=2; must be provided a blue and red frame ; for a single spectrograph if (nfile ne 2) then begin splog,'Incorrect number of files provided to mdrp_fluxcal, skipping!' return endif ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; End quality control checks ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Set up display if specified, otherwise set plot file 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 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 if (keyword_set(plotname)) then plotname=plotname $ else begin ; Set default plotname plotname=fileandpath(ml_psfilename(outnames[0]),path=path) path=path+'qa/' plotname=path+plotname if file_test(path,/directory) eq 0 then spawn, '\mkdir -p '+path plotname=ml_strreplace(plotname,'mgFFrame-b','mgFFrame-spec') endelse splog,'Printing to ',plotname dfpsplot, plotname, /color endelse lat = 32+46/60.+49/3600.d lon = -(105.+49/60.+13/3600.d) plateid = lonarr(nfile) mjd = lonarr(nfile) camname = strarr(nfile) camid = strarr(nfile) expnum = lonarr(nfile) spectroid = lonarr(nfile) npixarr = lonarr(nfile) dcnra = fltarr(nfile) dcndec = fltarr(nfile) jd = dblarr(nfile) ha = fltarr(nfile) tptemp = fltarr(nfile) taibeg = dblarr(nfile) exptimes = fltarr(nfile) for ifile=0, nfile-1 do begin infile=lookforgzip(innames[ifile]) ml_mgframeread, infile, hdr=hdr plateid[ifile] = strtrim(sxpar(hdr, 'PLATEID'),2) mjd[ifile] = strtrim(sxpar(hdr, 'MJD'),2) camname[ifile] = strtrim(sxpar(hdr, 'CAMERAS'),2) spectroid[ifile] = strmid(camname[ifile],1,1) camid[ifile] = strmid(camname[ifile],0,1) expnum[ifile] = sxpar(hdr, 'EXPOSURE') ; hdr1 = headfits(innames[ifile],exten='FLUX') npixarr[ifile] = sxpar(hdr, 'NAXIS1') dcnra[ifile] = sxpar(hdr,'MGDRA') dcndec[ifile] = sxpar(hdr,'MGDDEC') taibeg[ifile] = sxpar(hdr,'TAI-BEG') exptimes[ifile] = sxpar(hdr,'EXPTIME') jd[ifile] = (sxpar(hdr,'TAI-BEG')+sxpar(hdr,'EXPTIME')/2.)/86400.+2400000.5 tptemp[ifile] = sxpar(hdr,'TRANSPAR') endfor maxmjd = max(mjd) ; Initialize seed to something predictable seed=expnum[0] junk=randomu(seed,1) ;print,'SEEDS: ',expnum[0],junk uexp = uniq(expnum, bsort(expnum)) exposures = expnum[uexp] nexp = n_elements(uexp) ; Vectors to contain seeing and transparency seeing=fltarr(nexp) transparency=tptemp[uexp] for iexp=0,nexp-1 do begin ind = where(expnum eq exposures[iexp]) platestr = string(plateid[0],format='(i0.0)') mjdstr = string(mjd[0],format='(i0.0)') cogimg_dir = djs_filepath('', root_dir=ml_getenv('MANGA_SPECTRO_REDUX'), subdir=[drpver,platestr,mjdstr]) cogimg_file= file_search(cogimg_dir,'cogimg-'+mjdstr+'-'+string(exposures[iexp],format='(i8.8)')+'.fits*',count=ct) if ct eq 0 then $ coaddgimg,plateid[ind[0]],mjd[ind[0]],exposures[iexp] psf=psfparams(plateid[ind[0]],mjd[ind[0]],exposures[iexp],transparency=transpar) seeing[iexp]=psf.fwhm endfor goodexp = where(transparency gt 0.05,ngoodexp) if ngoodexp lt nexp then begin splog,'There are '+string(nexp-ngoodexp,format='(i0.0)')+' exposures with near-zero transparency. So, I cannot do flux calibration.' return endif ml_mgframeread,innames[0],slitmap=slitmap,hdr=thishdr ; Determine whether this is a MaNGA Dither or APOGEE lead exposure ; if *any* exposures are APOGEE led, this will set the method of ; flux calibration. thismode=sxpar(thishdr,'SRVYMODE') ; If srvymode doesn't contain MaNGA, treat as APOGEE-led if (strmatch(thismode,'*MaNGA*',/fold) eq 0) then srvymode='APOGEE' else srvymode = 'MaNGA' plugged=where(slitmap.plugstatus eq 0) eq2hor,slitmap[plugged[0]].cenra*(dblarr(nfile)+1),slitmap[plugged[0]].cendec*(dblarr(nfile)+1),jd,alt,az,ha,lat=lat,lon=lon,altitude=2788. usp = uniq(spectroid,bsort(spectroid)) if n_elements(usp) ne 1 then message, 'Run mdrp_fluxcal per spectrograph' ;Truncate the slitmap to one of the spectrograph slitmap = slitmap[where(slitmap.spectrographid eq spectroid[0])] ; iphoto = where(strtrim(slitmap.targettype,2) eq 'standard' and slitmap.plugstatus eq 0,nphoto) iphoto = where((slitmap.manga_target2 and $ 2L^20+2L^21+2L^22+2L^23+2L^24+2L^25+2L^26+2L^27) ne 0 $ and slitmap.plugstatus eq 0 $ and ((slitmap.psfmag[2] gt 0 and slitmap.psfmag[2] lt 30) or $ (slitmap.psfmag[1] gt 0 and slitmap.psfmag[1] lt 30 $ and (slitmap.manga_target2 and $ sdss_flagval('MANGA_TARGET2','STELLIB_GAIA')) ne 0) ), nphoto) isky = where(strmatch(slitmap.fibertype,'*SKY*',/fold_case) and strmatch(slitmap.targettype,'SKY*',/fold_case),nsky) if (nphoto EQ 0) then begin splog, 'WARNING: NO standard stars in this spectrograph for flux calibration' return endif npix = max(npixarr) loglam = dblarr(npix, nfile, nphoto) objflux = fltarr(npix, nfile, nphoto) objivar = fltarr(npix, nfile, nphoto) dispimg = fltarr(npix, nfile, nphoto) objmask = fltarr(npix, nfile, nphoto) ; airmass = fltarr(npix, nfile, nphoto1+nphoto2) for ifile=0,nfile-1 do begin ml_mgframeread,innames[ifile],iphoto,objflux=objflux1,objivar=objivar1,mask=mask1,loglam=loglam1,$ wset=wset1, lsfpost=lsfpost, lsfpre=lsfpre, superflat=superflat,hdr=hdr ; ml_mgframeread,innames[ifile],isky,objflux=skyflux1,objivar=skyivar1,mask=skymask1,loglam=skyloglam1,wset=skywset1, dispset=skydispset1,dispimg=skydispimg1,superflat=superflatsky ; DRL 12/19/19 ; Note that to avoid rewriting many parts of flux cal, we will ; continue to call the pre-pixel LSF 'disp' in this code dispimg1=lsfpre exptime = sxpar(hdr,'EXPTIME') ; convert flux and ivar from flat-field e- to actual e-/sec. objflux1 = objflux1*superflat/exptime objivar1 = objivar1/(superflat*superflat+(superflat eq 0))*(exptime*exptime) ; skyflux1 = skyflux1*superflatsky/exptime ; skyivar1 = skyivar1/(superflatsky*superflatsky+(superflatsky eq 0))*(exptime*exptime) ; Make a map of the size of each pixel in delta-(log10-Angstroms). ; Re-normalize the flux to ADU/(dloglam). ; Re-normalize the dispersion from /(raw pixel) to /(new pixel). correct_dlam,objflux1,objivar1,wset1,dlam=dloglam correct_dlam,dispimg1,0,wset1,dlam=dloglam,/inverse ; correct_dlam,skyflux1,skyivar1,skywset1,dlam=dloglam ; correct_dlam,skydispimg1,0,skywset1,dlam=dloglam,/inverse objivar1 = objivar1 * spflux_goodfiber(mask1) ; skyivar1 = skyivar1 * spflux_goodfiber(skymask1) if 0 then begin ; temporarily disable all sky subtraction in this code. ; Do a global sky subtraction to subtract residual skies. tmpskyarr = skyflux1*0.0 for i=0,nsky-1 do tmpskyarr[*,i] = interpol(skyflux1[*,i],skyloglam1[*,i],skyloglam1[*,0]) badmask = skyivar1 eq 0.0 tmpmaskarr = badmask*0 for i=0,nsky-1 do tmpmaskarr[*,i] = interpol(badmask[*,i],skyloglam1[*,i],skyloglam1[*,0]) newbadmask = tmpmaskarr gt 0 mediansky = dblarr(npixarr[ifile]) medianskyivar = dblarr(npixarr[ifile]) for k=0,npixarr[ifile]-1 do begin good = where(newbadmask[k,*] eq 0,ngood) if ngood lt 5 then continue mediansky[k] = median(tmpskyarr[k,good],/even) medianskyivar[k] = 1/ml_medianerr(tmpskyarr[k,good])^2 endfor for i=0,nphoto-1 do begin individualsky = interpol(mediansky,skyloglam1[*,0],loglam1[*,i]) individualskyivar = interpol(medianskyivar,skyloglam1[*,0],loglam1[*,i]) objflux1[*,i] = objflux1[*,i]-individualsky objivar1[*,i] = objivar1[*,i]*individualskyivar/(objivar1[*,i]+individualskyivar + ((objivar1[*,i]+individualskyivar) eq 0.0)) endfor endif loglam[0:npixarr[ifile]-1,ifile,*] = loglam1 ;it wont do having a tail of zeros in the wavelength so add some dummy values if (npix GT npixarr[ifile]) then begin dllam=loglam1[npixarr[ifile]-1,*]-loglam1[npixarr[ifile]-2,*] for j=0, nphoto-1 do $ loglam[npixarr[ifile]:*,ifile,j] = loglam1[npixarr[ifile]-1,j]+dllam[0,j]*(1+findgen(npix-npixarr[ifile])) endif objflux[0:npixarr[ifile]-1,ifile,*] = objflux1 ;hopefully the inverse variance of 0 of non-filled objects will indicate the uselessness ; of the extra objivar[0:npixarr[ifile]-1,ifile,*] = skymask(objivar1, mask1, mask1) dispimg[0:npixarr[ifile]-1,ifile,*] = dispimg1 objmask[0:npixarr[ifile]-1,ifile,*] = mask1 endfor lambda = 10^loglam energy = cc*planck/(lambda* 1.d-8)*1.d17 ; put energy also in units of 10^(-17) erg. ; dlambda = lambda-shift(lambda,[1,0]) ; because flux was in e- per pixel of 1e-4 in log(lambda) dlambda = 1.e-4*alog(10.)*lambda ; this is the pixel width in Angstrom unitfactor = energy/dlambda/mirrorarea objflux = objflux*unitfactor ; observed engery per sec per Angstrom per cm^2 * 10^17 objivar = objivar/unitfactor/unitfactor thru=mrdfits(getenv('MANGADRP_DIR')+'/etc/meanthrupt.fits',1) mjdselect= where(thru.mjd lt maxmjd , n_goodmjd) if n_goodmjd eq 0 then message,'No default throughput curve is available for this MJD :'+string(maxmjd,format='(i0.0)') tsel = mjdselect[n_goodmjd-1] if spectroid[0] eq 1 then init_thru = thru[tsel].sp1thrupt*thru[tsel].sp1corr if spectroid[0] eq 2 then init_thru = thru[tsel].sp2thrupt*thru[tsel].sp2corr thrupt = objflux*0.0 blues = where(camid eq 'b',nblue) if nblue gt 0 then thrupt[*,blues,*] = interpol(init_thru[*,0],thru[0].loglam,loglam[*,blues,*]) > 0.0 reds = where(camid eq 'r',nred) if nred gt 0 then thrupt[*,reds,*] = interpol(init_thru[*,1],thru[0].loglam,loglam[*,reds,*]) > 0.0 objflux = objflux/(thrupt+(thrupt eq 0)) objivar = objivar*thrupt*thrupt frlplug = slitmap[iphoto].frlplug uuplug = uniq(frlplug,bsort(frlplug)) starplug = frlplug[uuplug] nstar = n_elements(starplug) euler,slitmap[iphoto[uuplug]].ra,slitmap[iphoto[uuplug]].dec,l,b,1 sfd_ebv = dust_getval(l,b,/interp) starttime=systime(/sec) ; For MaStar exposures with exposure time shorter than 180, first we derive the velocity offset between blue and red cameras for short exposures. ; waveshiftflag = 0 totshift = 0.0 ; loglam0 = loglam if exptime le 180 then begin sn = objflux*sqrt(objivar) for istar=0,nstar-1 do begin bind = where(slitmap[iphoto].frlplug eq starplug[istar],nfiber) ss = bsort(slitmap[iphoto[bind]].fnum) bind = bind[ss] fibersn = total(total(sn[*,*,bind],1),1) tmp = max(fibersn,reffiber) ; temporary, reffiber reset later if max(fibersn) lt 1.e-6 then begin print,'No valid data for this ferrule in this exposure.' continue endif if n_elements(allid) eq 0 then allid=bind[reffiber] else allid=[allid,bind[reffiber]] endfor n_id = n_elements(allid) redvel = fltarr(n_id) redverr = fltarr(n_id) rederrcode = intarr(n_id) bluevel = fltarr(n_id) blueverr = fltarr(n_id) blueerrcode = intarr(n_id) bluecut = alog10(3622.) redcut = alog10(10354.) for iter=0,1 do begin for ii=0,n_id-1 do begin ind =where(loglam[*,0,allid[ii]] gt bluecut and loglam[*,0,allid[ii]] lt redcut and objivar[*,0,allid[ii]] gt 0) bluevel[ii] = spflux_velocity(loglam[ind,0,allid[ii]],objflux[ind,0,allid[ii]],objivar[ind,0,allid[ii]],dispimg[ind,0,allid[ii]],template=template,kindx=kindx,errcode=errcode,verr=verr) blueerrcode[ii]=errcode blueverr[ii]=verr ind =where(loglam[*,1,allid[ii]] gt bluecut and loglam[*,1,allid[ii]] lt redcut and objivar[*,1,allid[ii]] gt 0) redvel[ii] = spflux_velocity(loglam[ind,1,allid[ii]],objflux[ind,1,allid[ii]],objivar[ind,1,allid[ii]],dispimg[ind,1,allid[ii]],template=template,kindx=kindx,errcode=errcode,verr=verr) rederrcode[ii]=errcode redverr[ii]=verr endfor gvel = where(blueverr gt 0 and redverr gt 0,n_gvel) if n_gvel ge 3 then begin brdiff = bluevel[gvel]-redvel[gvel] brdifferr = sqrt(blueverr[gvel]^2+redverr[gvel]^2) brdiffmean = total(brdiff/brdifferr^2)/total(1/brdifferr^2) brdiffmeanerr = sqrt(total((brdiff-brdiffmean)^2/brdifferr^4))/total(1/brdifferr^2) keep = where(abs(brdiff-brdiffmean) lt brdifferr*3,nkeep) if nkeep ge 3 then begin brdiffmean = total(brdiff[keep]/brdifferr[keep]^2)/total(1/brdifferr[keep]^2) brdiffmeanerr = sqrt(total((brdiff[keep]-brdiffmean)^2/brdifferr[keep]^4))/total(1/brdifferr[keep]^2) endif ; Adjust the wavelength solution only if it is significantly different from zero. print,'Iteration ',iter,': Error-weighted Mean Velocity Offset:',brdiffmean,'+/-',brdiffmeanerr,'km/s',format='(a,i1,a,f7.1,a,f6.1,a)' if abs(brdiffmean) gt 2*brdiffmeanerr then begin ; If brdiffmean is positive, it means we need to reduce the wavelength at a given pixel location. waveshiftflag = 1 ; velshift = brdiffmean ; loglam[*,0,*] = loglam[*,0,*]-alog10(1+velshift*1.e5/cc) ; shifting blue camera to match with red camera. This is inaccurate as the shift in more constant in pixel unit rather than velocity unit. dlam = 10^loglam[1:npix-1,0,*]-10^loglam[0:npix-2,0,*] dlam = [dlam[0,*,*],dlam] pixelshift = 4800.*(brdiffmean*1.e5/cc)/1.10 dAngstrom = pixelshift*dlam ; The accuracies of the hardcoded numbers (4800 & 1.10) here do not matter much as we iterate once to get the exact shift in unit of pixels. The 4800*(velshift*1.e5/cc)/1.10 just converts the effective velocity shift to the shift in unit of pixels. We do not know what effective wavelength to adopt as the velocity is a compromised result between the two ends of the spectrum. Therefore, we need to iterate once to get to a more accurate result. In the 2nd iteration, the error introduced by the inaccurate effective wavelength becomes much smaller. loglam[*,0,*] = alog10(10^loglam[*,0,*]-dAngstrom) totshift = totshift+pixelshift endif else break endif endfor endif ; dAngstromTot = reform(10^loglam0[*,0,*]-10^loglam[*,0,*]) print,'Total shift (pixels):',totshift for istar=0,nstar-1 do begin starttime_istar=systime(/sec) print,'Processing standard star in Ferrule '+string(starplug[istar],format='(i0.0)') bind = where(slitmap[iphoto].frlplug eq starplug[istar],nfiber) ss = bsort(slitmap[iphoto[bind]].fnum) bind = bind[ss] ; slitmap[iphoto[bind]] gives the fibers in the bundle sorted in order of fnum ; bind also correspond to the row numbers in objflux and objivar. ; as those correspond to part of the original frame as read in with indices iphoto. platescale = 3.62730/60. ; (mm/arcsec) as given in Gunn et al. 2006 fiberx = slitmap[iphoto[bind]].xpmm/platescale fibery = slitmap[iphoto[bind]].ypmm/platescale ; Now identify which fibers belong to the central 7 fiber group in the bundle. This is for cases when we have standard star in a science IFU bundle. ; xmedian = median(slitmap[iphoto[bind]].xpmm)/platescale ; ymedian = median(slitmap[iphoto[bind]].ypmm)/platescale ; dist = sqrt((fiberx-xmedian)^2 + (fibery-ymedian)^2) ; tmp = min(dist,center) ; innermask = (sqrt((fiberx-fiberx[center])^2+(fibery-fibery[center])^2) lt 3) and (slitmap[iphoto[bind]].gbu eq 0) ; innerfiber = where(innermask,ninner) ; bind = bind[innerfiber] ; nfiber = ninner lam0 = 3500. if NOT (keyword_set(centerwave) and keyword_set(dwave)) then begin wave1 = [3500.,4000.,4500.,5000.,5500.,6500.,7500.,9000.] wave2 = [4000.,4500.,5000.,5500.,6500.,7500.,9000.,10500.] ; wave1 = [5000.] ; wave2 = [5500.] nwave = n_elements(wave1) centerwave = (wave1+wave2)/2. dwave = wave2-wave1 endif else begin if n_elements(centerwave) ne n_elements(dwave) then message,'centerwave and dwave must have the same number of elements.' nwave = n_elements(centerwave) wave1 = centerwave-dwave/2. wave2 = centerwave+dwave/2. endelse ; The following part of code do not work for multiple exposures yet. If we input multiple exposures, then we would need to select multiple sets of inner fibers. sn = objflux*sqrt(objivar) fibersn = total(total(sn[*,*,bind],1),1) tmp = max(fibersn,reffiber) ; temporary, reffiber reset later if max(fibersn) lt 1.e-6 then begin print,'No valid data for this ferrule in this exposure.' continue endif dist = sqrt((fiberx-fiberx[reffiber])^2+(fibery-fibery[reffiber])^2) ssdist = bsort(dist) innerfiber=ssdist[0:6] innerfiber = innerfiber[where(slitmap[iphoto[bind[innerfiber]]].gbu eq 0,ninner)] bind = bind[innerfiber] nfiber = ninner fibersn = fibersn[innerfiber] tmp = max(fibersn,reffiber) ; temporary, reffiber reset later-- though it won't change in principle. thisbundle = starplug[istar] thisfiber = iphoto[bind[reffiber]] plottitle = 'PLATE=' + string(plateid[0], format='(i0.0)') $ + ' MJD=' + string(maxmjd, format='(i0.0)') $ + ' Spectro-Photo Star' $ + ' Bundle '+ strtrim(thisbundle,2) $ + ' Fiber ' + strtrim(thisfiber,2) if (slitmap[iphoto[bind[reffiber]]].manga_target2 and sdss_flagval('MANGA_TARGET2','STELLIB_GAIA')) eq 0 then magid = 2 else magid=1 if slitmap[iphoto[bind[reffiber]]].psfmag[magid] lt 0.01 then message,'Error: This standard star has invalid psfmag. Cannot do flux calibration!!!' obsinfo = {plateid: 0L,ifudesign: 0, mangaid:' ',expnum: 0L,taibeg:0.0d, exptime:0.0d, cenra:0.d,cendec:0.0d, mjd: 0L,mgdra:0.0,mgddec:0.0,xfocal:0.0,yfocal:0.0,drp2qual:0L} obsinfo = replicate(obsinfo,nexp) obsinfo.plateid = plateid[0] obsinfo.ifudesign = slitmap[iphoto[bind[0]]].ifudesign obsinfo.mangaid = slitmap[iphoto[bind[0]]].mangaid obsinfo.cenra = slitmap[iphoto[bind[0]]].cenra obsinfo.cendec = slitmap[iphoto[bind[0]]].cendec obsinfo.xfocal = slitmap[iphoto[bind[0]]].xfocal obsinfo.yfocal = slitmap[iphoto[bind[0]]].yfocal obsinfo.expnum = exposures for iexp=0,nexp-1 do begin ind = where(expnum eq obsinfo[iexp].expnum) obsinfo[iexp].mgdra = dcnra[ind[0]] obsinfo[iexp].mgddec = dcndec[ind[0]] obsinfo[iexp].taibeg = taibeg[ind[0]] obsinfo[iexp].exptime = exptimes[ind[0]] obsinfo[iexp].mjd = mjd[ind[0]] endfor ; Construct flux_arr, ivar_arr, mask_arr, disp_arr, predisp_arr comloglam = alog10(ml_setwcalib()) nlampix = n_elements(comloglam) smscale =11 smflux = ml_smoothivar(objflux[*,*,bind],objivar[*,*,bind],[smscale,1,1]) smivar = (smooth(objivar[*,*,bind],[smscale,1,1],/edge_truncate,/nan) > 0)*smscale flux_arr = fltarr(nlampix,nfiber) ivar_arr = fltarr(nlampix,nfiber) mask_arr = lonarr(nlampix,nfiber) disp_arr = fltarr(nlampix,nfiber) predisp_arr = fltarr(nlampix,nfiber) for ifiber=0,nfiber-1 do begin bluesmflux = interpol(smflux[*,0,ifiber],loglam[*,0,bind[ifiber]],comloglam) bluesmivar = interpol(smivar[*,0,ifiber],loglam[*,0,bind[ifiber]],comloglam) redsmflux = interpol(smflux[*,1,ifiber],loglam[*,1,bind[ifiber]],comloglam) redsmivar = interpol(smivar[*,1,ifiber],loglam[*,1,bind[ifiber]],comloglam) bluedisp = interpol(dispimg[*,0,bind[ifiber]],loglam[*,0,bind[ifiber]],comloglam) reddisp = interpol(dispimg[*,1,bind[ifiber]],loglam[*,1,bind[ifiber]],comloglam) ivar_arr[*,ifiber] = bluesmivar+redsmivar flux_arr[*,ifiber] = (bluesmflux*bluesmivar+redsmflux*redsmivar)/(ivar_arr[*,ifiber]+ (ivar_arr[*,ifiber] eq 0)) disp_arr[*,ifiber] = (bluedisp*bluesmivar+reddisp*redsmivar)/(ivar_arr[*,ifiber] + (ivar_arr[*,ifiber] eq 0)) mask_arr[*,ifiber] = ivar_arr[*,ifiber] eq 0 endfor predisp_arr = disp_arr ; Just a placeholder. For determining the basic covering fraction. The predisp_arr do not matter. We will reset these later. ;stop ; starflux = mdrp_extractstar(flux_arr,ivar_arr,mask_arr,disp_arr,predisp_arr,obsinfo,slitmap[iphoto[bind]],starivar=starivar,starmask=starmask,stardisp=stardisp,starpredisp=starpredisp, reffiberindex=reffiberindex,allexpfiberid=allexpfiberid,allexpflux=allexpflux, allexpivar=allexpivar, allexpmask=allexpmask, allexpdisp=allexpdisp, allexppredisp=allexppredisp, mjdindex=mjdindex, allexpgrpfiberid=allexpgrpfiberid,allcovfn=allcovfn,scalebest=scalebest,rotationbest=rotationbest,mjdqual=mjdqual,allexpqual=allexpqual,allchi2=allchi2,expresult=expresult) ; Need to propagate seeing, transparency, factor, xbest, ybest, scalebest, rotationbest, and fluxratio_chi2 for each exposure, ; Here the fiber order in allcovfn is not necessarily the same as the order in flux_arr. Because the orders are sorted inside mdrp_extractstar. The way to figure out which index in allcovfn corresponds to the reffiberindex is to compare allexpgrpfiberid with allexpfiberid. Allexpgrpfiberid gives the fiberid in the same order as done in mdrp_extractstar. Therefore, it corresponds to the order given by allcovfn. refindex = where(allexpgrpfiberid[*,0] eq allexpfiberid[0]) ; Here the zero is for the first exposure. (preparation for implementing multiple exposures in the future.) refcovfn = interpol(allcovfn[*,refindex,0],comloglam,loglam[*,*,bind[reffiber]]) allexpflux = objflux[*,*,bind[reffiberindex]]/refcovfn allexpivar = objivar[*,*,bind[reffiberindex]]*refcovfn^2 ; stop if srvymode eq 'MaNGA' then ebv = sfd_ebv[istar] else ebv = 0.0 tmp = typingmodule(allexpflux,loglam[*,*,bind[reffiber]],allexpivar,dispimg[*,*,bind[reffiber]],ebv,slitmap[iphoto[bind[reffiber]]].psfmag,kindx=thisindx,modflux=modflux,plottitle=plottitle,template=template,thekfile=thekfile,targetflag=slitmap[iphoto[bind[reffiber]]].manga_target2) ; At this point, the modflux is the theoretical model. For MaNGA, it would contain reddening, for MaStar, it would not. modelspec = modflux*refcovfn ; includes the covering fraction. ; We will generate a correction curve for each camera/spectrograph/exposure. ; The two cameras from the same spectrograph and the same exposure should share a common positional offset and a common DAR model, but they will provide separate data points and get separate correction vector. ; !!! IMPORTANT !!! ; The flux in the following structure should be the original flux from ; mgSFrame. It should NOT be corrected for the covering fraction. The ; modflux does not include covering fraction while modelspec does. ; The ratio of the two would give the covfn. These are ; for the convenience of other diagnostic codes that are dependent ; on this assumption. result={plate:0L,mjd:0L,exposure:0L,dcnra:0.0,dcndec:0.0,frlplug:0L,ra:0.0d,dec:0.0d,$ fiberid:slitmap[iphoto[bind]].fiberid,fnum:slitmap[iphoto[bind]].fnum,$ psfmag:fltarr(5),manga_target2:slitmap[iphoto[bind[0]]].manga_target2, $ factor:0.0,xbest:0.0,ybest:0.0,$ ; fiberoffx:fiberoffx,fiberoffy:fiberoffy,$ ; offsetx:offsetx,offsety:offsety,$ scalebest:scalebest,rotationbest:rotationbest,fluxratio_chi2:allchi2[0],modelspec:modelspec,$ modflux:modflux, flux:objflux[*,*,bind[reffiberindex]],$ ivar:objivar[*,*,bind[reffiberindex]], disp:dispimg[*,*,bind[reffiber]],$ mratio:fltarr(npix,2),mrativar:fltarr(npix,2),$ loglam:loglam[*,*,bind[reffiber]],$ ;fluxratio:fluxratio,fluxratioerr:fluxratioerr,predictratio:predict,$ av:0.0,used:0} result = create_struct(result,thisindx) if n_elements(totresult) eq 0 then begin totresult = replicate(result,nexp*nstar) endif ; The following line is a reminder that we need to properly set dispimg if we indeed use multiple exposures together. Right now, it still only works for single exposure. for iexp=0,nexp-1 do result[iexp].disp = dispimg[*,*,bind[reffiber]] result.plate = plateid[uexp] result.mjd = mjd[uexp] result.exposure=expnum[uexp] result.dcnra = dcnra[uexp] result.dcndec = dcndec[uexp] result.frlplug=starplug[istar] result.ra = slitmap[iphoto[bind[0]]].ra result.dec = slitmap[iphoto[bind[0]]].dec result.psfmag = slitmap[iphoto[bind[0]]].psfmag result.factor= expresult.factor result.xbest= expresult.xbest result.ybest= expresult.ybest ; result.scalebest=scalebest ; result.rotationbest=rotationbest ; result.fluxratio_chi2 = expresult.fluxratio_chi2 ; totresult[iexp*nstar+istar]=result totresult[istar*nexp:(istar+1)*nexp-1] = result endfor print,'Total time used to process all exposures of all stars:',systime(/sec)-starttime if srvymode ne 'MaNGA' then begin iexp = 0 ; Assuming a single exposure is input. Need to modify if we were to use ; multiple exposures. ind = where(expnum eq exposures[iexp]) eq2hor,totresult.ra,totresult.dec,jd[ind[0]],alt,az,ha,lat=lat,lon=lon,altitude=2788. airmass = 1/sin(alt*!dtor) ratio = totresult.flux/(totresult.modelspec + (totresult.modelspec eq 0)) ratioivar = totresult.ivar*totresult.modelspec^2 sizes = size(ratio,/dimen) smscale = 11 smratio = ml_smoothivar(ratio,ratioivar,[smscale,1,1]) smivar = (smooth(ratioivar,[smscale,1,1],/edge_truncate,/nan) > 0)*smscale ; this sums the ivar of the pixels in the smoothing window. smivar = smivar/smscale ; this makes the variance larger to include the covariance. commonloglam = alog10(ml_setwcalib()) ncomloglam = n_elements(commonloglam) interpratio = fltarr(ncomloglam,sizes[1],sizes[2]) interpivar = fltarr(ncomloglam,sizes[1],sizes[2]) for i=0,sizes[2]-1 do begin &$ for j=0,sizes[1]-1 do begin &$ interpratio[*,j,i] = interpol(smratio[*,j,i],totresult[i].loglam[*,j],commonloglam) &$ interpivar[*,j,i] = interpol(smivar[*,j,i],totresult[i].loglam[*,j],commonloglam) > 0 &$ endfor &$ endfor galextcorr = ratio*0.0 loginterpratio = -2.5*alog10(interpratio) ; loginterpivar = 2.5/alog(10.)*interpivar/(interpratio + (interpratio eq 0.0)) ; this become negative, for which the loginterpratio would be NaN. loginterpivar = (alog(10.)/2.5*interpratio)^2 *interpivar sensibleratio = total(total(finite(loginterpratio),1),1) gt 1000 and total(total(finite(loginterpivar),1),1) gt 1000 photometry_arr = strarr(n_elements(totresult)) colorgood = intarr(n_elements(totresult)) c1b = 1.1517 & c1r = 0.6104 c2b = -0.0871 & c2r = -0.0170 c3b = -0.0333 & c3r = -0.0026 c4b = 0.0173 & c4r = -0.0017 c5b = -0.0230 & c5r = -0.0078 c6b = 0.0006 & c6r = 0.00005 c7b = 0.0043 & c7r = 0.0006 ebvarr = fltarr(n_elements(totresult)) for tt=0,n_elements(totresult)-1 do begin photometry_arr[tt] = 'sdss' ; Both APASS and SDSS are using sdss photometry system if (totresult[tt].manga_target2 and sdss_flagval('MANGA_TARGET2',['STELLIB_PS1','STD_PS1_COM'])) ne 0 then begin photometry_arr[tt]='ps1' colorgood[tt] = totresult[tt].psfmag[1] gt 0 and totresult[tt].psfmag[3] gt 0 if colorgood[tt] then begin modcolor = totresult[tt].mag[1]-totresult[tt].mag[3] psfcolor = totresult[tt].psfmag[1]-totresult[tt].psfmag[3] colordiff = psfcolor - modcolor ebvarr[tt] = colordiff/(3.172-1.682) endif endif if (totresult[tt].manga_target2 and sdss_flagval('MANGA_TARGET2','STELLIB_GAIA')) ne 0 then begin photometry_arr[tt]='gaiadr1' colorgood[tt] = 0 endif if (totresult[tt].manga_target2 and sdss_flagval('MANGA_TARGET2','STELLIB_GAIADR2')) ne 0 then begin photometry_arr[tt]='gaiadr2' colorgood[tt] = totresult[tt].psfmag[2] gt 0 and totresult[tt].psfmag[3] gt 0 if colorgood[tt] then begin modcolor = totresult[tt].mag[2]-totresult[tt].mag[3] psfcolor = totresult[tt].psfmag[2]-totresult[tt].psfmag[3] aa = c5b-c5r + (c7b-c7r)*modcolor bb = c1b-c1r + (c2b-c2r)*modcolor + (c3b-c3r)*modcolor*modcolor $ +(c4b-c4r)*modcolor*modcolor*modcolor dd = (modcolor - psfcolor) < 0 ebvarr[tt] = -(bb-sqrt(bb^2-4*aa*dd))/(2*aa)/3.1 endif endif if photometry_arr[tt] eq 'sdss' then begin colorgood[tt] = totresult[tt].psfmag[1] gt 0 and totresult[tt].psfmag[3] gt 0 if colorgood[tt] then begin modcolor = totresult[tt].mag[1]-totresult[tt].mag[3] psfcolor = totresult[tt].psfmag[1]-totresult[tt].psfmag[3] colordiff = psfcolor - modcolor ebvarr[tt] = colordiff/(3.303-1.698) endif endif endfor ; colorstar = where((totresult.manga_target2 and sdss_flagval('MANGA_TARGET2',['STELLIB_PS1','STELLIB_APASS','STELLIB_SDSS'])) ne 0 and totresult.psfmag[1] gt 0 and totresult.psfmag[3] gt 0 and sensibleratio, n_colorstar) colorstar = where(colorgood and sensibleratio, n_colorstar) if n_colorstar eq 0 then colorstar = indgen(n_elements(totresult)) n_nonzeroivar = total(loginterpivar[*,*,colorstar] gt 0,3) aveinterpratio = total(loginterpratio[*,*,colorstar]*(loginterpivar[*,*,colorstar] gt 0),3)/n_nonzeroivar aveinterpivar = 1/(total(1/loginterpivar[*,*,colorstar],3)/n_nonzeroivar) good = where(n_nonzeroivar eq max(n_nonzeroivar) and finite(aveinterpivar),n_good) if n_good eq 0 then message,'Error: No good pixels for an average correction curve.' mask = spflux_masklines(commonloglam,/stellar,/telluric) good = where(n_nonzeroivar eq max(n_nonzeroivar) and finite(aveinterpratio) and aveinterpratio ne 0.0 and [mask,mask] eq 0) ; extinct,10^commonloglam,aratio,/ccm, Rv=3.1 fm_unred,10^commonloglam,commonloglam*0+1,-1.0/3.1,freddened aratio = -2.5*alog10(freddened) ; xx = 1.e4/(10^commonloglam) ; aratio1 = ax(xx,/ccm) ; aratio2 = bx(xx,/ccm) if n_colorstar ne 0 then begin aveav = mean(ebvarr[colorstar])*3.1 ; modcolor_gi = totresult.mag[1]-totresult.mag[3] ; psfcolor_gi = totresult.psfmag[1]-totresult.psfmag[3] ; colordiff = psfcolor_gi - modcolor_gi ; ps1 = where((totresult[colorstar].manga_target2 and sdss_flagval('MANGA_TARGET2','STELLIB_PS1')) ne 0,n_ps1) ; sdssapass=where((totresult[colorstar].manga_target2 and sdss_flagval('MANGA_TARGET2',['STELLIB_SDSS','STELLIB_APASS'])) ne 0,n_sa) ; if n_ps1 gt 0 then ebvarr[ps1] = colordiff[ps1]/(3.172-1.682) ; if n_sa gt 0 then ebvarr[sdssapass] = colordiff[sdssapass]/(3.303-1.698) ; aveebv = median(ebvarr[colorstar]) ; aveav = aveebv*3.1 endif else begin refcoeff = linfit(([aratio,aratio])[good],aveinterpratio[good],measure_errors=1/sqrt(aveinterpivar[good]),chisq=chisq,sigma=sigma,yfit=yfit) aveav = refcoeff[1] endelse alambda = ([aratio,aratio])*aveav diff = loginterpratio-rebin(aveinterpratio,ncomloglam,sizes[1],sizes[2]) differr = sqrt(1/loginterpivar+1/rebin(aveinterpivar,ncomloglam,sizes[1],sizes[2])) coeff = fltarr(2,sizes[2]) sig = fltarr(2,sizes[2]) coeff3 = fltarr(3,sizes[2]) sig3 = fltarr(3,sizes[2]) residual = loginterpratio*0.0 residual3 = loginterpratio*0.0 chi2 = fltarr(sizes[2]) chi3 = fltarr(sizes[2]) for istar=0,nstar-1 do begin &$ ; extinct,10^totresult[istar].loglam,aratio_native,/ccm,Rv=3.1 &$ fm_unred,10^totresult[istar].loglam,totresult.loglam*0+1,-1/3.1,fred_native &$ aratio_native = -2.5*alog10(fred_native) &$ good = where(finite(diff[*,*,istar]) and diff[*,*,istar] ne 0.0 and finite(differr[*,*,istar]) and [mask,mask] eq 0) &$ coeff[*,istar] = linfit(([aratio,aratio])[good],(diff[*,*,istar])[good],measure_errors=(differr[*,*,istar])[good],chisq=chisq,sigma=sigma,yfit=yfit) &$ sig[*,istar] = sigma &$ chi2[istar] = chisq &$ ; xxx = transpose([[([aratio1,aratio1])[good]],[([aratio2,aratio2])[good]]]) &$ ; regressres = regress(xxx,(diff[*,*,istar])[good],measure_errors=(differr[*,*,istar])[good],chisq=chisq3,sigma=sigma3,yfit=yfit3,const=const) &$ ; coeff3[*,istar]= [const,regressres[0],regressres[1]] &$ ; sig3[1:2,istar] = sigma3 &$ ; chi3[istar] = chisq3 &$ residual[*,*,istar] = diff[*,*,istar]-coeff[0,istar]-coeff[1,istar]*[aratio,aratio] &$ ; residual3[*,*,istar] = diff[*,*,istar]-const-coeff3[1,istar]*[aratio1,aratio1]-coeff3[2,istar]*[aratio2,aratio2] &$ ; build the correction vector for the native wavelength grid. galextcorr[*,*,istar] = coeff[0,istar]+(coeff[1,istar]+aveav)*aratio_native &$ endfor totresult.av = reform(coeff[1,*]+aveav) ; Renormalize the extinction-corrected model. for istar=0,nstar-1 do begin ; covfn = totresult[istar].modelspec/totresult[istar].modflux ; allexpflux = totresult[istar].flux/covfn ; allexpivar = totresult[istar].ivar*covfn*covfn ; tmp = typingmodule(allexpflux,totresult[istar].loglam,allexpivar,totresult[istar].disp,totresult[istar].av/3.1,totresult[istar].psfmag,kindx=thisindx,modflux=modflux,plottitle=plottitle,template=template,thekfile=thekfile,targetflag=totresult[istar].manga_target2) ; Now include reddening in modflux and modelspec ; Renormalize the model to the observed magnitudes as we now have reddening. ; Compute SCALEFAC = (plugmap flux) / (flux of the model spectrum) thismodel = totresult[istar].modflux*10^(-0.4*galextcorr[*,*,istar]) modelspec = totresult[istar].modelspec*10^(-0.4*galextcorr[*,*,istar]) tmploglam = 3.4771d0 + lindgen(5644) * 1.d-4 s0 = bsort(totresult[istar].loglam) model = interpol(thismodel[s0],(totresult[istar].loglam)[s0],tmploglam) wavevec = 10.d0^tmploglam flambda2fnu = wavevec^2 / 2.99792e18 CASE photometry_arr[istar] OF 'ps1' : BEGIN fthru = filter_thru(model* flambda2fnu, waveimg=wavevec, /toair, filternames=['ps1_g.txt','ps1_r.txt','ps1_i.txt','ps1_z.txt']) thismag = -2.5*alog10(fthru) - (48.6-2.5*17) thismag = [[-999.],[thismag]] scalefac = 10.^((thismag[2]-totresult[istar].psfmag[2])/2.5) totresult[istar].mag = reform(thismag)+totresult[istar].psfmag[2]-thismag[2] END 'sdss': BEGIN fthru = filter_thru(model * flambda2fnu, waveimg=wavevec, /toair) thismag = -2.5 * alog10(fthru) - (48.6-2.5*17) scalefac = 10.^((thismag[2]-totresult[istar].psfmag[2])/2.5) totresult[istar].mag = reform(thismag)+totresult[istar].psfmag[2]-thismag[2] END 'gaiadr2': BEGIN fthru = filter_thru(model * flambda2fnu, waveimg=wavevec, /toair, filternames=['gaia_G.dat','gaia_BP.dat','gaia_RP.dat']) thismag = -2.5*alog10(fthru) - (48.6-2.5*17) ; convert from AB to Vega mag to compare with the Gaia DR2 mags thismag = thismag + [[25.6884-25.7916],[25.3514-25.3862],[24.7619-25.1162]] thismag = [[-999.],[thismag],[-999.]] scalefac = 10.^((thismag[1]-totresult[istar].psfmag[1])/2.5) totresult[istar].mag = reform(thismag)+totresult[istar].psfmag[1]-thismag[1] END 'gaiadr1': BEGIN fthru = filter_thru(model * flambda2fnu, waveimg=wavevec, /toair, filternames=['gaia_G.dat']) thismag = -2.5*alog10(fthru) - (48.6-2.5*17) thismag = [-999.,thismag,-999.,-999.,-999.] scalefac = 10.^((thismag[1]-totresult[istar].psfmag[1])/2.5) totresult[istar].mag = reform(thismag)+totresult[istar].psfmag[1]-thismag[1] END else: splog,'Error: standard stars have uncertain photometry source.' ENDCASE thismodel = thismodel * scalefac modelspec = modelspec * scalefac totresult[istar].modflux = thismodel totresult[istar].modelspec = modelspec endfor endif else totresult.av = sfd_ebv*3.1 ratiocorr = totresult.flux/(totresult.modelspec+(totresult.modelspec eq 0)) ratioivarcorr = totresult.ivar*totresult.modelspec^2 restloglam = totresult.loglam-transpose(rebin([alog10(1+totresult.z)],nstar,2,npix)) stellarmask = (spflux_masklines(restloglam,/stellar,hwidth=4.e-4) and 10^restloglam lt 3760) $ or (spflux_masklines(restloglam,/stellar,hwidth=8.e-4) and 10^restloglam gt 3760) ratioivarcorr = ratioivarcorr * (1-stellarmask) totresult.mratio = ratiocorr totresult.mrativar = ratioivarcorr ; Final loop over input frames ; Read in the frames, apply the calibrations, and write ; out the calibrated mgFFrame files. Note that these will ; have units of flux/pixel rather than flux per some uniform ; wavelength range. psffac = fltarr(nexp) for iexp=0,nexp-1 do begin ; looping through different exposures indexp = where(expnum eq exposures[iexp],nexpfile) inds = where(totresult.exposure eq exposures[iexp],n_inds) ; select only good stars --- rejecting those with S/N less than 1/3 of the median and those with Chi2 greater than 3 times the median if n_inds eq 0 then begin splog,'Error: No valid stars for exposure '+string(exposures[iexp],format='(i8.8)') return endif if n_inds gt 1 then begin blueratio = fltarr(n_inds) redratio = fltarr(n_inds) for jj=0,n_inds-1 do begin bluemid = where(totresult[inds[jj]].loglam[*,0] gt alog10(5300.) and totresult[inds[jj]].loglam[*,0] lt alog10(5350.)) redmid = where(totresult[inds[jj]].loglam[*,1] gt alog10(7800.) and totresult[inds[jj]].loglam[*,1] lt alog10(8000.)) blueratio[jj] = median(totresult[inds[jj]].mratio[bluemid,0],/even) redratio[jj] = median(totresult[inds[jj]].mratio[redmid,1],/even) endfor bluediff = abs(blueratio-median(blueratio,/even))/median(blueratio,/even) reddiff = abs(redratio-median(redratio,/even))/median(redratio,/even) bluescatter = (median(bluediff,/even)/0.6745) > 0.033; sigma scatter derived using median absolute deviation; set a minimum scatter of 3.3% so that we don't reject anything within +/-10% of median. Otherwise, in some cases the median absolute deviation could be super small (when there are three nearly identical values.) redscatter = (median(reddiff,/even)/0.6745) > 0.033 gstar = where(totresult[inds].sn_median gt $ median(totresult[inds].sn_median,/even)/3. and $ totresult[inds].linechi2/totresult[inds].linedof lt $ median(totresult[inds].linechi2/totresult[inds].linedof,/even)*3 and $ totresult[inds].fluxratio_chi2 lt (100 > median(totresult[inds].fluxratio_chi2,/even)) and $ (bluediff lt 3*bluescatter and reddiff lt 3*redscatter) and $ max(totresult[inds].mrativar[*,0],dimen=1) gt 0 and max(totresult[inds].mrativar[*,1],dimen=1) gt 0, ngstar) ; take out a lower-order polynomial relative to the mean correction vector ; This will hopefully average out the systematics introduced in imperfect ; model match. if ngstar gt 1 then begin flatarrblue = spflux_mratio_flatten(totresult[inds[gstar]].loglam[*,0],$ totresult[inds[gstar]].mratio[*,0],totresult[inds[gstar]].mrativar[*,0],pres=pres_b) flatarrred = spflux_mratio_flatten(totresult[inds[gstar]].loglam[*,1],$ totresult[inds[gstar]].mratio[*,1],totresult[inds[gstar]].mrativar[*,1],pres=pres_r) psffac[iexp] = median(totresult[inds[gstar]].factor,/even) endif else begin if ngstar eq 1 then begin flatarrblue = fltarr(npix)+1.0 flatarrred = fltarr(npix)+1.0 psffac[iexp] = totresult[inds[gstar]].factor endif else message,'There are stars but none of them gives good calibration.' endelse endif else begin ; this means inds has only 1 elements, as we have excluded the 0 element possibility above. gstar = 0 ngstar = 1 flatarrblue = fltarr(npix)+1.0 flatarrred = fltarr(npix)+1.0 psffac[iexp] = totresult[inds].factor endelse fmratio = fltarr(npix, 2, ngstar) fmrativar = fltarr(npix, 2, ngstar) fmratio[*,0,*] = totresult[inds[gstar]].mratio[*,0]/flatarrblue fmrativar[*,0,*] = totresult[inds[gstar]].mrativar[*,0]*flatarrblue^2 fmratio[*,1,*] = totresult[inds[gstar]].mratio[*,1]/flatarrred fmrativar[*,1,*] = totresult[inds[gstar]].mrativar[*,1]*flatarrred^2 ; In the following step, we do not need to specify the /stellar keyword because the masking has been already been applied to fmrativar. And it was done with the right redshift for each star. sset_b = spflux_bspline(totresult[inds[gstar]].loglam[*,0],$ fmratio[*,0,*],fmrativar[*,0,*],everyn=160*ngstar) mask = spflux_masklines(totresult[inds[gstar]].loglam[*,1],/telluric) tlc = where(mask,complement=smo) fmrativarcopy = fmrativar[*,1,*] fmrativarcopy[tlc] = 0.0 sset_rt = spflux_bspline(totresult[inds[gstar]].loglam[*,1],$ fmratio[*,1,*],fmrativar[*,1,*],$ everyn=1.0*ngstar);,disp=(totresult[inds[gstar]].disp[*,1])) ; without masking out stellar features for the telluric region.) sset_r = spflux_bspline(totresult[inds[gstar]].loglam[*,1],$ fmratio[*,1,*],fmrativarcopy,everyn=160*ngstar,/telluric) ; mask out telluric regions here. Stellar features have already been masked out in fmrativar. ; maskall = spflux_masklines() ; tlcpart = bspline_valu(totresult.loglam[*,1],sset_rt,x2=totresult.disp[*,1]) ; smopart = bspline_valu(totresult.loglam[*,1],sset_r) ; smopart[tlc] = tlcpart[tlc] !p.multi=[0,1,2] yrange = [0,2] djs_plot,10^totresult[inds[gstar]].loglam[*,0],fmratio[*,0,*],ps=1,symsize=0.1,xran=[3500,6500],xst=1,xtitle='Wavelength (A)',ytitle='Correction Vector',yran=yrange,yst=1 djs_oplot,10^totresult[0].loglam[*,0],bspline_valu(totresult[0].loglam[*,0],sset_b),color='red',ps=1,symsize=0.1 tlcversion = bspline_valu(totresult[inds[gstar]].loglam[*,1],sset_rt);,x2=totresult[inds[gstar]].disp[*,1]) smoversion = bspline_valu(totresult[inds[gstar]].loglam[*,1],sset_r) combinever = smoversion*0 combinever[tlc] = tlcversion[tlc] combinever[smo] = smoversion[smo] djs_plot,10^totresult[inds[gstar]].loglam[*,1],fmratio[*,1,*],ps=1,symsize=0.1,xran=[5800,10500],xst=1,xtitle='Wavelength (A)',ytitle='Correction Vector',yran=yrange,yst=1 djs_oplot,10^totresult[inds[gstar]].loglam[*,1],combinever,color='red',ps=1,symsize=0.1 !p.multi=0 ; thrupt_blue = bspline_valu(thru.loglam,sset_b)*init_thru[*,0] ; smo_standardlam = bspline_valu(thru.loglam,sset_r) ; tlc_standardlam = bspline_valu(thru.loglam,sset_rt) ; mask_standard = spflux_masklines(thru.loglam,/telluric) ; tlc_id = where(mask_standard,complement=smo_id) ; comb_standardlam= smo_standardlam*0 ; comb_standardlam[smo_id] = smo_standardlam[smo_id] ; comb_standardlam[tlc_id] = tlc_standardlam[tlc_id] ; thrupt_red = com_standardlam*init_thru[*,1] thruptref_blue = bspline_valu(thru[0].loglam[1000],sset_b)*init_thru[1000,0] thruptnorm_blue = bspline_valu(thru[0].loglam[2535],sset_b)*init_thru[2535,0] thruptref_red = bspline_valu(thru[0].loglam[4000],sset_r)*init_thru[4000,1] thruptnorm_red = bspline_valu(thru[0].loglam[2535],sset_r)*init_thru[2535,1] bluenormalized = thruptref_blue/thruptnorm_blue*0.08 rednormalized = thruptref_red/thruptnorm_red*0.18 mjdstr=string(mjd[indexp[0]],format='(i0.0)') if NOT keyword_set(writename) then writename='' if NOT keyword_set(nowritefile) then begin tmpname=fileandpath(ml_psfilename(outnames[0]),path=path) filename=path+'fluxcal-sp'+string(spectroid[0],format='(i0.0)')+'-'+string(exposures[iexp],format='(i8.8)')+writename+'.fits' totresult[inds[gstar]].used = 1 mwrfits,totresult[inds],filename,/create spawn, ['gzip', '-f', filename], /noshell endif for ifile=0,nexpfile-1 do begin ; looping through different cameras of the same exposure ml_mgframeread,innames[indexp[ifile]],loglam=loglam1, wset=wset1, objflux=tempflux, $ objivar=tempivar, mask=tempmask, lsfpost=lsfpost, lsfpre=lsfpre, slitmap=slitfull, ximg=tempximg, $ slithdr=slithdr, superflat=superflat1, airset=airset, $ skyflux=tempsky, supersky=supersky, hdr=thishdr, adderr=adderr if camid[indexp[ifile]] eq 'b' and waveshiftflag then begin ; loglam1 = loglam1-alog10(1+velshift*1.e5/cc) ; loglam1 = alog10(10^loglam1-dAngstromTot) ; Modify wset1 ; wset1.coeff[0,*] = wset1.coeff[0,*]-alog10(1+velshift*1.e5/cc) ; traceset2xy,wset1,xdummy,loglam_redundant ncoeff = (size(wset1.coeff,/dimen))[0] nrow = (size(wset1.coeff,/dimen))[1] xdummy = dindgen(wset1.xmax+1)#(dblarr(nrow)+1) dlam = 10^loglam1-shift(10^loglam1,[1,0]) dlam[0,*] = dlam[1,*] loglam1 = alog10(10^loglam1-totshift*dlam) xy2traceset,xdummy,loglam1,wsetnew,ncoeff=ncoeff,func='legendre' traceset2xy,wsetnew,dummy,newloglam print,'Maximum offset introduced by the re-fitting of wset:',max(abs(10^loglam1-10^newloglam)),' Angstrom',format='(a,f,a)' if max(abs(10^loglam1-10^newloglam)) lt 0.05 then begin wset1=wsetnew endif else begin ; if refitting of wset introduced too much error, then use the simple shift to the constant term in the existing wset. v_over_c = totshift*1.10/4800. wset1.coeff[0,*] = wset1.coeff[0,*]-alog10(1+v_over_c) endelse endif ; DRL 12/19/19 ; Note that to avoid rewriting many parts of flux cal, we will ; continue to call the pre-pixel LSF 'disp' in this code tempdisp=lsfpre ; convert flux and ivar from flat-field e- to actual e-/sec. exptime = sxpar(thishdr,'EXPTIME') tempflux = tempflux*superflat1/exptime tempivar = tempivar/(superflat1*superflat1+(superflat1 eq 0))*(exptime*exptime) tempsky = tempsky*superflat1/exptime ; Convert the super-sky model to actual e-/sec ; Use the superflat of the first fiber (which is ok even if ; fiber missing 'cuz the pipeline is smart and has filled in ; reasonable values) airx=0; Reset the x values so the loop doesn't save them traceset2xy,airset,airx,airloglam tempsupsky=supersky.flux*interpol(superflat1[*,0],airloglam[*,0],alog10(supersky.airwave))/exptime ; Make a map of the size of each pixel in 1e-4 (log10-Angstroms), ; and re-normalize the flux, ivar, sky to electrons/(dloglam) correct_dlam, tempflux, tempivar, wset1, appliedval=appliedval correct_dlam, tempsky, 0, wset1 correct_dlam, tempdisp, 0, wset1, /inverse ; Convert super-sky model by interpolating from first fiber again mdrp_divideflat,tempsupsky,interpol(appliedval[*,0],airloglam[*,0],alog10(supersky.airwave)),minval=0 lambda = 10^loglam1 energy = cc*planck/(lambda* 1.d-8)*1.d17 ; put energy also in units of 10^(-17) erg. ; because flux was in e- per pixel of 1e-4 in log(lambda) dlambda = 1.e-4*alog(10.)*lambda ; this is the pixel width in Angstrom unitfactor = energy/dlambda/mirrorarea tempflux = tempflux*unitfactor ; observed engery per sec per Angstrom per cm^2 * 10^17 tempivar = tempivar/unitfactor/unitfactor tempsky = tempsky*unitfactor ; Convert super-sky model by interpolating from first fiber factor tempsupsky = tempsupsky*interpol(unitfactor[*,0],airloglam[*,0],alog10(supersky.airwave)) splog,'Applying flux calibration vectors to mgSFrame inputs.' CASE camid[indexp[ifile]] OF 'b' : BEGIN thrupt = interpol(init_thru[*,0],thru[0].loglam,loglam1) > 0.0 bcalib = float(bspline_valu(loglam1,sset_b))*thrupt ; bluemask = spflux_masklines(loglam1,hwidth=12.e-4,/stellar) ; bcalib = djs_maskinterp(bcalib,bluemask,loglam1,iaxis=0,/const) minval = 0.05*median(bcalib) ; Apply calibration vector to flux, inverse variance, and sky extensions ; This will put it in units of flux/Angstrom mdrp_divideflat, tempflux, invvar=tempivar, bcalib, minval=minval mdrp_divideflat, tempsky, bcalib, minval=minval ; Convert super-sky model mdrp_divideflat, tempsupsky, interpol(bcalib[*,0],airloglam[*,0],alog10(supersky.airwave)),minval=minval ; Flag pixels with bad flux calibration factor tempmask = tempmask OR ((bcalib LE minval) * drp2pixmask_bits('BADFLUXFACTOR')) ; Add normalized throughput information sxaddpar, thishdr, 'NMTHRUPT', bluenormalized, 'Normalized throughput in the blue at 4354.1160A' END 'r' : BEGIN thrupt = interpol(init_thru[*,1],thru[0].loglam,loglam1) > 0.0 smocalib = bspline_valu(loglam1,sset_r) tlccalib = bspline_valu(loglam1,sset_rt);,x2=tempdisp) tlc = where(spflux_masklines(loglam1,/telluric),complement=smo) rcalib = smocalib*0 rcalib[tlc] = tlccalib[tlc] rcalib[smo] = smocalib[smo] rcalib *= thrupt ; redmask = spflux_masklines(loglam1,hwidth=12.e-4,/stellar) ; rcalib = djs_maskinterp(rcalib,redmask,loglam1,iaxis=0,/const) minval = 0.05*median(rcalib) ; Apply calibration vector to flux, inverse variance, and sky extensions ; This will put it in units of flux/Angstrom mdrp_divideflat, tempflux, invvar=tempivar, rcalib, minval=minval mdrp_divideflat, tempsky, rcalib, minval=minval ; Convert super-sky model mdrp_divideflat, tempsupsky, interpol(rcalib[*,0],airloglam[*,0],alog10(supersky.airwave)),minval=minval ; Flag pixels with bad flux calibration factor tempmask = tempmask OR ((rcalib LE minval) * drp2pixmask_bits('BADFLUXFACTOR')) ; Add normalized throughput information sxaddpar, thishdr, 'NMTHRUPT', rednormalized, 'Normalized throughput in the red at 8687.6037A' END ENDCASE ; Add flux information keywords sxaddpar, thishdr, 'BSCALE', 1., ' Intensity unit scaling' sxaddpar, thishdr, 'BZERO', 0., ' Intensity zeropoint' sxaddpar, thishdr, 'BUNIT', '1E-17 erg/s/cm^2/Ang/fiber', ' Specific intensity (per fiber-area)' ; Add PSF factor, seeing sxaddpar, thishdr, 'PSFFAC', psffac[iexp], 'Best-fit PSF size relative to guider measurements' sxaddpar, thishdr, 'SEEING', seeing[iexp], ' Best guider seeing FWHM in Arcsec' ; Add information about the template file used sxaddpar, thishdr, 'TPLDATA', fileandpath(thekfile),'Flux calibration template file used' ; Add quality control keywords drp2qual=long(fxpar(thishdr,'DRP2QUAL')) ; Add new flagging cases here ; None yet! ; Write back to header sxaddpar, thishdr, 'DRP2QUAL', drp2qual, 'DRP-2d quality bitmask' ; Remove old name sxdelpar, thishdr, 'EXTNAME' ; Add flux calibrated super-sky spectra back into the ; supersky table supersky.flux=tempsupsky ; Write sky-subtracted spectra to disk following the MaNGA data model ml_mwrfits, dummyext, outnames[indexp[ifile]],hdr=thishdr, /create ; Blank ext 0 with full header ml_mwrfits, tempflux, outnames[indexp[ifile]], hdr=thishdr, extname='FLUX' ;sky subtracted flux ml_mwrfits, tempivar, outnames[indexp[ifile]], extname='IVAR' ; sky subtracted inverse variance ml_mwrfits, tempmask, outnames[indexp[ifile]], extname='MASK' ; final pixel mask ml_mwrfits, wset1, outnames[indexp[ifile]], extname='WSET' ;trace-set for wavelength sol; wset ml_mwrfits, lsfpost, outnames[indexp[ifile]], extname='LSFPOST' ;post-pixel LSF ml_mwrfits, lsfpre, outnames[indexp[ifile]], extname='LSFPRE' ;pre-pixel LSF ml_mwrfits, slitfull, outnames[indexp[ifile]], extname='SLITMAP' ;slitmap used ml_mwrfits, tempximg, outnames[indexp[ifile]], extname='XPOS' ;x pos on CCD ml_mwrfits, superflat1, outnames[indexp[ifile]], extname='SUPERFLAT' ;superflat vector from quartz lamps ml_mwrfits, airset, outnames[indexp[ifile]], extname='AIRSET' ; Observed-frame wavelength traceset ml_mwrfits, tempsky, outnames[indexp[ifile]], extname='SKY' ;sky flux ml_mwrfits, supersky, outnames[indexp[ifile]], extname='SUPERSKY'; Supersky table ; gzip output file spawn, ['gzip', '-f', outnames[indexp[ifile]]], /noshell endfor ; finish looping through different cameras of the same exposure endfor ; finish looping through different exposures. ; Close out plots if ~keyword_set(visual) then begin dfpsclose ; Convert to PDF and remove old PS file ;stop spawn, strcompress('ps2pdf '+plotname+' '+ml_strreplace(plotname,'.ps','.pdf')) spawn, strcompress('rm -f '+plotname) !p.multi=0 !p=origp !x=origx !y=origy endif ; Take out the trash heap_gc if keyword_set(time) then print,'total lasted time= ',systime(1)-time0,' secs' end