;+ ; NAME: ; evilmap ; ; PURPOSE: ; Correlate the output from a plate scane (from "evilscan") with object ; positions in a plugmap file. ; ; CALLING SEQUENCE: ; evilmap, [ platenum, scanfile=scanfile, plugmap=plugmap, $ ; maxiter=maxiter, scandir=scandir, plugdir=plugdir, outdir=outdir, $ ; wtime=wtime, copydir=copydir ] ; ; INPUTS: ; ; OPTIONAL INPUTS: ; platenum - Plate number for scan file; look for a the last file when ; sorting on names "fiberScan-XXXX-YYYYY-ZZ.par" where XXXX ; is the PLATENUM, YYYYY is the MJD, and ZZ is a scan index. ; scanfile - Specify the scan file explicity, rather than deriving its ; name from PLATENUM. Must specify either PLATENUM or SCANFILE. ; plugmap - Input plug map file. If not specified, then look for the ; first file matching "plPlugMapP-XXXX*.par" where XXXX is ; the plate number. ; maxiter - Maximum number of iterations for finding CCD to PLUGMAP ; matches; default to 4. ; scandir - Directory in which to search for SCANFILE; default to ; '$HOME/scan/*/' if it exists, otherwise to './'. ; plugdir - Directory in which to search for PLUGMAP; default to ; '/home/plates/*/' if it exists, otherwise to './'. ; outdir - Output directory for all output files; default to the ; same as SCANDIR. ; wtime - Time to wait after each plot; default to 2.0 sec. Set equal ; to 0 to disable all plotting to the screen. ; copydir - If set, then copy the output files to this directory using ; "scp1" copy. Normally, one want to set this to ; 'sdsshost:/data/spectro/plugMapM'. ; ; OUTPUT: ; ; OPTIONAL OUTPUTS: ; ; COMMENTS: ; The following files are input: ; plPlugMapP-*.fits - Input plugmap file ; ; The following files are output: ; plPlugMapM-*.par - Modified version of input plugmap file ; Unplugged-*.ps - Figure with unplugged fibers indicated by squares ; ; The throughput is normalized by a quadratic function as a function of ; distance from the plate center, since the mapping station sees less ; flux at larger distances. ; ; Note that for the sake of speed, SCANFILE is not read as a Yanny param ; file, but it is read as a fixed-format text file. This means that if one ; edits that file, one must be certain to maintain nearly the exact format. ; ; Presently, we search only a small parameter space for the scale of the ; CCD in terms of mm/pixel. This parameter space is specified as the ; loop over SCALX in BEST_OFFSET_SCALE. Should the camera be moved relative ; to the plate surface, or the focus of the camera dramatically change, ; then this loop will need to be changed. ; ; EXAMPLES: ; ; BUGS: ; If only 1 missing fiber, should we assign it? ; Copy output files to sdsshost? ; Has Matt changed the gain on the video camera?? Saturating now?? ; Write more of the plots to disk? Names of those plot files? ; ; PROCEDURES CALLED: ; djs_mean() ; djs_oplot ; readfmt ; yanny_read ; yanny_write ; ; INTERNAL SUPPORT ROUTINES: ; def_fib_struc() ; def_event_struc() ; histogram_time_offset() ; localize_fibers ; qplot1 ; qplot4 ; warp_apply ; plot_unplugged ; best_offset_scale ; warp_initial_guess ; add_hdr_line ; ; REVISION HISTORY: ; 17-Aug-1999 Written by Doug Finkbeiner & David Schlegel, Apache Point. ; 30-Aug-1999 Re-wrote step 1 to allow for multiple events at the same ; pixel (DJS). This has slowed the code down quite a bit. ;- ;------------------------------------------------------------------------------ ; Define structure for an event. An "event" is a pixel getting bright then ; faint again. It occurs in a single pixel but spans multiple time steps. ; function def_event_struc event1 = {event_struc, $ tmean: 0.0, $ tsig: 0.0, $ flux: 0.0, $ spectroID: 0B, $ x: 0, $ y: 0 } return, event1 end ;------------------------------------------------------------------------------ ; Find time offset between fibers by auto-correlating the time history. ; Find an accurate centroid to this peak by gaussian-fitting. ; Note that in real units, the time offset is this return value times TBIN. function histogram_time_offset, thist nstep = 50 corr = fltarr(nstep) ; Do not compute the first few points, since the zero-lag autocorrelation ; is always largest. for i=5, nstep-1 do begin corr[i] = mean(thist * shift(thist,i)) endfor yfit = gaussfit(findgen(nstep), corr, coeff) offset = coeff[1] return, offset end ;------------------------------------------------------------------------------ ; Fill in the FIBS.X, FIBS.Y positions by localizing the positions in ; time and space. ; The PLOT keyword is for debuggin purposes, and makes a little movie ; of the events (points) and fiber positions (square). pro localize_fibers, tmean, deltat, fibs, pixarr, rejdist=rejdist, plot=plot if (NOT keyword_set(rejdist)) then rejdist = 1.6 ; in CCD pixels for ifib=0, N_elements(fibs)-1 do begin w = where(abs(fibs[ifib].tstamp - tmean) LT deltat, ct) if (ct GT 0) then begin xw = pixarr[w].x yw = pixarr[w].y xcen = median([xw]) ycen = median([yw]) igood = where(abs(xw - xcen) LT rejdist $ AND abs(yw - ycen) LT rejdist, ngood) xcen = djs_mean( xw[igood] ) ycen = djs_mean( yw[igood] ) fibs[ifib].x = xcen fibs[ifib].y = ycen if (keyword_set(plot)) then begin plot, [xw], [yw], psym=3, xrange=[0,640], yrange=[0,600], $ xstyle=1, ystyle=1 oplot, [xcen], [ycen], psym=6, symsize=2 wait, 0.25 endif endif else begin print, ' Failed to localize any points near timestamp ' + $ string(fibs[ifib].tstamp) endelse endfor return end ;------------------------------------------------------------------------------ pro qplot1, timevec, thist, fibs, xrange=xrange plot, timevec, thist, xrange=xrange, /xstyle, $ xtitle='Time stamp', ytitle='Events' usersym, cos(findgen(11)*!pi/5), sin(findgen(11)*!pi/5) ; open circ oplot, [fibs.tstamp], [0*fibs.tstamp+2.5], psym=8 ibad = where(fibs.good EQ 0) if (ibad[0] NE -1) then $ djs_oplot, [fibs[ibad].tstamp], [0*fibs[ibad].tstamp+2.75], $ psym=2, symsize=1.5, color='red' return end ;------------------------------------------------------------------------------ pro qplot4, thist, fibs, tbin common evilmap_com, waittime !p.multi=[0,1,4] timevec = findgen(N_elements(thist)) * tbin tmax = max(timevec) qplot1, timevec, thist, fibs, xrange=[0.00*tmax,0.25*tmax] qplot1, timevec, thist, fibs, xrange=[0.25*tmax,0.50*tmax] qplot1, timevec, thist, fibs, xrange=[0.50*tmax,0.75*tmax] qplot1, timevec, thist, fibs, xrange=[0.75*tmax,1.00*tmax] !p.multi=[0,1,1] wait, waittime return end ;------------------------------------------------------------------------------ ; Compute the CCD positions after "warping" to the plugmap ; coordinate system pro warp_apply, xraw, yraw, k_x, k_y, xwarp, ywarp dims = size(k_x, /dimens) degree = dims[0]-1 ; Assume that k_x, k_y are [DEGREE+1,DEGREE+1] arrays xwarp = 0 * xraw ywarp = 0 * yraw for i=0, degree do begin for j=0, degree do begin xwarp = xwarp + k_x[i,j] * xraw^j * yraw^i ywarp = ywarp + k_y[i,j] * xraw^j * yraw^i endfor endfor end ;------------------------------------------------------------------------------ pro plot_unplugged, plugdat, plotfile=plotfile, title=title iobject = where(plugdat.holetype EQ 'OBJECT', nobject) iguide = where(plugdat.holetype EQ 'GUIDE' $ OR plugdat.holetype EQ 'COHERENT_SKY', nguide) iquality = where(plugdat.holetype EQ 'QUALITY', nquality) itrap = where(plugdat.holetype EQ 'LIGHT_TRAP', ntrap) pobj = plugdat[iobject] ; Determine which pixels are missing from plugmap list misslist = where(pobj.fiberid LE 0) if (keyword_set(plotfile)) then $ dfpsplot, plotfile, /square, /color plot, pobj.xfocal, pobj.yfocal, psym=1, symsize=0.5, $ xrange=[-350,350], yrange=[-350,350], xstyle=1, ystyle=1, $ title=title, xtitle='xFocal [mm]', ytitle='yFocal [mm]' oplot, [-320], [320], psym=1, symsize=0.5 xyouts, [-320], [320], ' Plugged fiber' if (misslist[0] NE -1) then begin djs_oplot, [pobj[misslist].xfocal], [pobj[misslist].yfocal], $ psym=2, symsize=2, color='red' endif djs_oplot, [-320], [300], psym=2, symsize=2, color='red' xyouts, [-320], [300], ' Unplugged fiber' if (nguide GT 0) then begin usersym, cos(findgen(21)*!pi/10), sin(findgen(21)*!pi/10) ; open circ djs_oplot, [plugdat[iguide].xfocal], [plugdat[iguide].yfocal], $ psym=8, symsize=3 xyouts, [plugdat[iguide].xfocal], [plugdat[iguide].yfocal-4], $ strtrim(string(plugdat[iguide].fiberid),2), alignment=0.5, charsize=1 endif oplot, [-320], [280], psym=8, symsize=3 xyouts, [-320], [280], ' GUIDE' if (ntrap GT 0) then begin djs_oplot, [plugdat[itrap].xfocal], [plugdat[itrap].yfocal], $ psym=6, symsize=1.5 endif oplot, [-320], [260], psym=6, symsize=1.5 xyouts, [-320], [260], ' LIGHT TRAP' if (nquality GT 0) then begin usersym, cos(findgen(21)*!pi/10), sin(findgen(21)*!pi/10), /fill ; circ djs_oplot, [plugdat[iquality].xfocal], [plugdat[iquality].yfocal], $ psym=8, symsize=1 endif oplot, [-320], [240], psym=8, symsize=1 xyouts, [-320], [240], ' QUALITY' if (keyword_set(plotfile)) then dfpsclose return end ;------------------------------------------------------------------------------ ; We pass this procedure X1 as the XFOCAL position from the plugmap file, ; and X2 as the observed fiber positions. We try rescaling X2 to XR, which ; maps the plugmap positions in units of mm. pro best_offset_scale, x1, x2, bestxoff, bestscalx common evilmap_com, waittime binsize = 0.5 ; Histogram bin size in millimeters (plugmap positions). ; Note that min fiber separation is 55-arcsec = 3.0 mm. nsmooth = 21 ; This smoothing scale is about 10 mm. bestcorr = -1.0 bestxoff = 0.0 bestscalx = 1.0 bestxaxis = [0] ; Saved only for plotting purposes besthist1 = [0] ; Saved only for plotting purposes besthist2 = [0] ; Saved only for plotting purposes ; With the current setup in Sep-1999, the scale is always ; about scalex=-1.475, scaly=-1.470 mm/pix. for scalx=-1.51, -1.44, 0.005 do begin xr = x2 * scalx xmin = min([x1,xr]) - (nsmooth+1)/2 xmax = max([x1,xr]) + (nsmooth+1)/2 ; Declare the following as type FLOAT to be sure HISTOGRAM creates ; exactly the same binning. hist1 = histogram(float(x1), min=xmin, max=xmax, binsize=binsize) hist2 = histogram(float(xr), min=xmin, max=xmax, binsize=binsize) hist1 = smooth(hist1+0.0, nsmooth) hist2 = smooth(hist2+0.0, nsmooth) ; Pad the histograms npad = N_elements(hist1) hist1 = [hist1, fltarr(npad)] ; Change the following padding for HIST2 to be absolutely certain ; that the number of bins is exactly the same as HIST1. ; hist2 = [hist2, fltarr(npad)] hist2 = [hist2, fltarr(n_elements(hist1)-n_elements(hist2))] ; LAGS indicates the amount to shift HIST1, so -LAGS is amount ; to shift HIST2 (e.g., amount to shift the observed fiber positions ; to match the plugmap positions in mm). ; The following two lines will allow for the absolutely maximal ; possible shifts between X1 and XR ; lagmin = fix( (min(xr)-max(x1)-5) / binsize) ; Align right sides of hist ; lagmax = fix( (max(xr)-min(x1)+5) / binsize) ; Align left sides of hist ; The following two lines will allow for lags that put the median XR ; position from the minimum value of X1 (left-most plug position) to ; the maximum value of X1 (right-most plug position) lagmin = fix( (median(xr)-max(x1)-5) / binsize) ; Align right side lagmax = fix( (median(xr)-min(x1)+5) / binsize) ; Align left side lagmax = lagmax > lagmin lags = indgen(lagmax-lagmin+1) + lagmin corr = c_correlate(hist1, hist2, lags) maxcorr = max(corr, i) if (maxcorr GT bestcorr) then begin bestxoff = lags[i] * binsize bestscalx = scalx bestcorr = maxcorr bestxaxis = xmin + findgen(N_elements(hist1)) * binsize besthist1 = hist1 besthist2 = hist2 endif endfor if (keyword_set(waittime)) then begin plot, bestxaxis, besthist1, psym=0, xrange=[-350,350], xstyle=1, $ xtitle='x/yFocal [mm]', ytitle='Histogram' djs_oplot, bestxaxis-bestxoff, besthist2, color='red', psym=0 wait, 1.5 endif print, ' Offset = ', bestxoff print, ' Scale = ', bestscalx print, ' Correlation = ', bestcorr return end ;------------------------------------------------------------------------------ pro warp_initial_guess, pobj, gfib, k_x, k_y print print, ' Initial guess in X:' x1 = pobj.xfocal x2 = gfib.x best_offset_scale, x1, x2, xoffset, scalx print print, ' Initial guess in Y:' y1 = pobj.yfocal y2 = gfib.y best_offset_scale, y1, y2, yoffset, scaly ; x,y values of fibers mapped to plug file focal plane (x,y) coords degree = 2 k_x = fltarr(degree+1,degree+1) k_x[0,0] = -xoffset k_x[0,1] = scalx k_y = fltarr(degree+1,degree+1) k_y[0,0] = -yoffset k_y[1,0] = scaly return end ;------------------------------------------------------------------------------ ; Add a comment line to HDR. Overwrite an existing line if the keyword ; is found as the first word of some line. Otherwise, add it before any ; occurence of 'typedef'. Or just append this line to the end. ; pro add_hdr_line, hdr, keyword, value addline = string(keyword,format='(A,T14)') + strtrim(string(value),2) ; If this keyword already exists, then overwrite its value i = (where(strmid(hdr,0,strlen(keyword)-1) EQ keyword))[0] if (i NE -1) then begin hdr[i] = addline endif else begin ; Add this field before the first typedef i = (where(strmid(hdr,0,7) EQ 'typedef'))[0] if (i NE -1) then begin hdr = [hdr[0:i-1], addline, hdr[i:N_elements(hdr)-1]] endif else begin hdr = [hdr, addline] endelse endelse return end ;------------------------------------------------------------------------------ ;------------------------------------------------------------------------------ pro evilmap, platenum, scanfile=scanfile, plugmap=plugmap, $ maxiter=maxiter, scandir=scandir, plugdir=plugdir, outdir=outdir, $ wtime=wtime, copydir=copydir fmapVersion = idlmapper_version() ; Set version for this program common evilmap_com, waittime if (NOT keyword_set(maxiter)) then maxiter=4 if (N_elements(scandir) EQ 0) then begin scandir = '$HOME/scan/*/' junk = findfile(scandir) if (junk[0] EQ '') then scandir = './' endif if (N_elements(plugdir) EQ 0) then begin plugdir = '/home/plates/*/' junk = findfile(plugdir) if (junk[0] EQ '') then plugdir = './' endif if (N_elements(wtime) EQ 0) then waittime = 2.0 $ else waittime = wtime ; Find the scanfile (but don't read it in yet) if (keyword_set(scanfile)) then begin ; Determine the plate number from the scan file name platenum = fix( strmid(scanfile,strpos(scanfile,'-')+1,4) ) endif else if (N_params() GT 0) then begin files = 'fiberScan-' + string(platenum,format='(i04.4)') + '-*.par' files = findfile(scandir+files, count=ct) if (ct GT 0) then begin ; Sort the scan files, taking the last one (which should be the ; most recent scan on the most recent MJD) ii = (reverse( sort( strmid(files,11) ) ))[0] j = rstrpos(files[ii], '/') if (j EQ -1) then begin scandir = './' scanfile = files[ii] endif else begin scandir = strmid(files[ii],0,j+1) scanfile = strmid(files[ii],j+1,strlen(files[ii])-j) endelse endif else begin message, 'No scan files found for plate number ' + string(platenum) endelse endif else begin print, 'Syntax - evilmap, [ platenum, scanfile=scanfile, plugmap=plugmap, wtime=wtime ]' return endelse if (N_elements(outdir) EQ 0) then outdir=scandir ; Find the plug map file (but don't read it in yet) platestr = string(platenum,format='(i04.4)') if (NOT keyword_set(plugmap)) then begin files = 'plPlugMapP-' + platestr + '*.par' files = findfile(plugdir+files, count=ct) print, files if (ct GT 1) then begin print, 'Several plug map files found for plate number ' $ + string(platenum) print, 'Using file ', files[0] files = files[0] ct = 1 endif if (ct EQ 1) then begin ; Take the only plugmap file j = rstrpos(files[0], '/') if (j EQ -1) then begin plugdir = './' plugmap = files[0] endif else begin plugdir = strmid(files[0],0,j+1) plugmap = strmid(files[0],j+1,strlen(files[0])-j) endelse endif if (ct EQ 0) then begin message, 'No plug map files found for plate number ' + string(platenum) endif endif print print, 'Setting directory for scan files to ', scandir print, 'Setting directory for plugmap files to ', plugdir print, 'Setting directory for output files to ', outdir print ; Construct the root name for the output files i1 = strpos(scanfile, '-') + 1 i2 = rstrpos(scanfile, '.') - 1 rootname = strmid(scanfile, i1, i2-i1+1) ; Read the ASCII scan file print, 'READING ', scandir+scanfile ; First set defaults for header values in case they aren't present fscanVersion = 'v?_?' pluggers = '' fscanMJD = 0L fscanId = 0 fscanMode = 'slow' fscanSpeed = 0 fscanBias = 0 motorId1 = 0 motorId2 = 0 ; Read the header get_lun, ilun openr, ilun, scandir+scanfile sline = '' stemp = '' readf, ilun, sline ww = str_sep(strtrim(sline,2), ' ') nskip = 0 while (strupcase(ww[0]) NE 'SCANPIX' $ AND NOT eof(ilun)) do begin if (nskip EQ 0) then header = sline $ else header = [header, sline] readf, ilun, sline ; Look for header cards that we want for the output plPlugMapM file ww = str_sep(strtrim(sline,2), ' ') if (ww[0] EQ 'fscanVersion') then $ reads, sline, stemp, fscanVersion, format='(A13,A)' if (ww[0] EQ 'pluggers') then $ reads, sline, stemp, pluggers, format='(A13,A)' if (ww[0] EQ 'fscanMJD') then $ reads, sline, stemp, fscanMJD, format='(A13,A)' if (ww[0] EQ 'fscanId') then $ reads, sline, stemp, fscanId, format='(A13,I)' if (ww[0] EQ 'fscanMode') then $ reads, sline, stemp, fscanMode, format='(A13,A)' if (ww[0] EQ 'fscanSpeed') then $ reads, sline, stemp, fscanSpeed, format='(A13,A)' if (ww[0] EQ 'motorId1') then $ reads, sline, stemp, motorId1, format='(A13,I)' if (ww[0] EQ 'motorId2') then $ reads, sline, stemp, motorId2, format='(A13,I)' if (ww[0] EQ 'fscanBias') then $ reads, sline, stemp, fscanBias, format='(A13,F)' nskip = nskip + 1 endwhile close, ilun free_lun, ilun ; Now read in the data lines readfmt, scandir+scanfile, 'A7,I3,I6,I6,F10.4,I4,I4,I4', $ junk, motornum, framenum, motorpos, tstamp, ypix, xpix, flux, skip=nskip junk = 0 ; Free memory if (keyword_set(waittime)) then begin i = where(motornum EQ 1) tt = framenum[i] + 0.0D mm = motorpos[i] + 0.0D res = linfit(tt, mm) plot, tt, mm - res[1]*tt - res[0], yrange=[-20,20], $ xtitle='Frame number', ytitle='Motor Position Residual' i = where(motornum EQ 2) tt = framenum[i] + 0.0D mm = motorpos[i] + 0.0D res = linfit(tt, mm) djs_oplot, tt, mm - res[1]*tt - res[0], color='blue' wait, waittime endif ; Subtract the bias value from all points, but retain the fluxes as integers print, 'Subtracting bias value of ', fscanBias flux = flux - fix(fscanBias) ; Choose which measure to use as fiber position. One could use the motor ; position (MOTORPOS), time stamp (TSTAMP), or frame number (FRAMENUM). ; Actually, it should always be the case that motor position is fairly ; accurate. ; Also set TBIN, which is used as the histogram binning size for ; rigorously solving for the fiber time difference. The code will ; basically work with TBIN set arbitrarily small, except that ; histogram_time_offset() would need to be modified. ; Also Set THWIDTH, which is the time half-width for grouping ; pixel data into an event. ; Note that the separation between fibers is about 150 motor steps. if (fscanMode EQ 'extreme') then begin print, 'Extremely fast scan: Using motor position for fiber position' tstamp = motorpos tbin = 5 ; in motor steps endif else if (fscanMode EQ 'fast') then begin print, 'Fast scan: Using motor position for fiber position' tstamp = motorpos tbin = 5 ; in motor steps endif else begin print, 'Slow scan: Using motor position for fiber position' tstamp = motorpos tbin = 5 ; in motor steps endelse thwidth = 8.0 * tbin framenum = 0 ; Free memory motorpos = 0 ; Free memory ; Apply software threshhold and trim all the data. ; Note that this should not be necessary, but it does result in fewer ; points to deal with. thresh = 35 ; Threshold after bias-subtraction wthresh = where(flux GE thresh) print, 'Trimming to ', N_elements(wthresh), $ ' data points above threshold of ', thresh motornum = motornum[wthresh] tstamp = tstamp[wthresh] xpix = xpix[wthresh] ypix = ypix[wthresh] flux = flux[wthresh] ; Read the cartridge ID and spectograph ID if ((motorId1 MOD 32) NE (motorId2 MOD 32) $ OR motorId1 EQ motorId2) then begin print, 'Failed sanity check on magnet ID values!' print, 'Check magnet IDs with "idtest"' print, 'Or manually edit the "motorId1" and "motorId2" values in' print, scandir+scanfile print, 'to read the cartridge number and the cartridge number plus 32.' message, 'Stop.' endif cartridgeId = motorId1 MOD 32 if (motorId1 LT 32) then begin ; This is the typical case, with motor#1 on spectro-2 ; Set spectroID=1 for motornum=2, and spectroID=2 for motornum=1 spectroID = 3 - motornum endif else begin ; This is the case where the motors are switched from their usual ; positions spectroID = motornum endelse ;--------------------------------------------------------------------------- ; STEP 1: Create pixel array in PIXARR structure ; ; Note that if a pixel is lit at very different time steps, then the ; following algorithm will split it into multiple events. ;--------------------------------------------------------------------------- print print, 'STEP 1: Creating event structure' ; Create the event structure. Make this structure as large as it could ; possibly be, then trim it down later npixtot = N_elements(xpix) pixarr = replicate(def_event_struc(), npixtot) ; Loop over each pixel, filling in PIXARR structure qpixdone = lonarr(npixtot) ; Set equal to 1 for each pixel that has been ; used already in the PIXARR structure nevent = 0 while (total(qpixdone) NE npixtot) do begin ; Select the next pixel location to consider ip = (where(qpixdone EQ 0))[0] ; Find all other data at that exact pixel location that are "unused" indx = where(xpix EQ xpix[ip] AND ypix EQ ypix[ip] $ AND qpixdone EQ 0, np) ; If this pixel is lit more than 20 timesteps, then assume that this ; is a hot pixel. Subtract away its median value, and trim to only ; values above a new threshold (set to 1.5 * median). ; Do the trimming by marking these pixels as "used" in the QPIXDONE ; array. if (np GT 20) then begin print, ' Hot pixel at ', xpix[ip], ypix[ip] medval = median( flux[indx] ) flux[indx] = flux[indx] - medval ibad = where(flux[indx] LT 1.5 * medval) qpixdone[indx[ibad]] = 1 ; Re-select "unused" pixels at this location. ; Note that there may now be none. indx = where(xpix EQ xpix[ip] AND ypix EQ ypix[ip] $ AND qpixdone EQ 0, np) endif if (np GT 0) then begin ; Localize to only those pixels within THWIDTH time of when this pixel ; is the brightest. junk = max(flux[indx], i0) tpeak = tstamp[indx[i0]] indx = indx[ where(tstamp[indx] GE tpeak-thwidth $ AND tstamp[indx] LE tpeak+thwidth) ] ; Take position from any one of the pixels (ie, the first one) i0 = indx[0] pixarr[nevent].spectroID = byte(spectroID[i0]) pixarr[nevent].x = xpix[i0] pixarr[nevent].y = ypix[i0] ; Compute the total flux, flux-weighted timestamp, and time dispersion pixarr[nevent].flux = total(flux[indx]) pixarr[nevent].tmean = total(flux[indx]*tstamp[indx]) $ / pixarr[nevent].flux pixarr[nevent].tsig = $ sqrt(total(flux[indx]*(tstamp[indx]-pixarr[nevent].tmean)^2) $ / pixarr[nevent].flux ) ; Add this event to the PIXARR only if it includes more than one ; timestamp, which will make TSIG greater than zero. ; If we do not update NEVENT, then this entry in PIXARR is ; over-written by the next event. if (pixarr[nevent].tsig GT 0.0) then nevent = nevent + 1 ; Mark these pixels as "used" qpixdone[indx] = 1 endif endwhile ; Trim the PIXARR structure to only non-zero elements if (nevent EQ 0) then message, 'No events found!' print, ' Number of events = ', nevent pixarr = pixarr[0:nevent-1] ;--------------------------------------------------------------------------- ; STEP 2: Find the timestamp sequence for fibers ; ; One side at a time ;--------------------------------------------------------------------------- print print, 'STEP 2: Find the timestamp sequence for fibers' side1 = where(pixarr.spectroID EQ 1) side2 = where(pixarr.spectroID EQ 2) ; For each side, shift their time-stamps to start at zero pixarr[side1].tmean = pixarr[side1].tmean - min(pixarr[side1].tmean) pixarr[side2].tmean = pixarr[side2].tmean - min(pixarr[side2].tmean) ; SIDE 1 tmean1 = pixarr[side1].tmean nbin = long( 1.05 * max(tmean1) / float(tbin) ) ; Add 5% padding thist1 = fltarr(nbin) for i=0L, N_elements(tmean1)-1 do $ thist1[tmean1[i]/tbin] = thist1[tmean1[i]/tbin] + 1 delt1 = histogram_time_offset(thist1) * tbin slit_position, thist1, delt1, 1, tbin, fibs1 if (keyword_set(waittime)) then begin qplot4, thist1, fibs1, tbin wait, waittime endif ; SIDE 2 tmean2 = pixarr[side2].tmean nbin = long( 1.05 * max(tmean2) / float(tbin) ) ; Add 5% padding thist2 = fltarr(nbin) for i=0L, N_elements(tmean2)-1 do $ thist2[tmean2[i]/tbin] = thist2[tmean2[i]/tbin] + 1 delt2 = histogram_time_offset(thist2) * tbin slit_position, thist2, delt2, 2, tbin, fibs2 if (keyword_set(waittime)) then begin qplot4, thist2, fibs2, tbin wait, waittime endif ;--------------------------------------------------------------------------- ; STEP 3: Localize X+Y fiber positions in time and space ; ; In time, we insist that each pixel event have a centroid within 0.25 ; times the event spacing. ;--------------------------------------------------------------------------- print print, 'STEP 3: Localize X+Y fiber positions in time and space' localize_fibers, tmean1, 0.25*delt1, fibs1, pixarr[side1] localize_fibers, tmean2, 0.25*delt2, fibs2, pixarr[side2] fibs = [fibs1, fibs2] gfib = fibs[where(fibs.good)] nfib = N_elements(gfib) ;--------------------------------------------------------------------------- ; Compute throughput of fibers by totalling the flux for all "events" ; within 2.0 pixels, and with a timestamp within DELTAT. rmax = 2.0 for j=0, nfib-1 do begin i = where(abs(pixarr.x - gfib[j].x) LT rmax $ AND abs(pixarr.y - gfib[j].y) LT rmax $ AND abs(pixarr.tmean - gfib[j].tstamp) LT 0.25*delt1) if (i[0] NE -1) then gfib[j].flux = total(pixarr[i].flux) $ else gfib[j].flux = 0 endfor ;--------------------------------------------------------------------------- ; STEP 4: Correlate fiber positions with plug map ;--------------------------------------------------------------------------- print print, 'STEP 4: Correlate fiber positions with plug map' ; Read plug map file and trim to only OBJECT entries. print, ' Reading plugmap file ', plugdir+plugmap yanny_read, plugdir+plugmap, pstruct, hdr=hdr, enums=enums, structs=structs plugdat = *pstruct[0] ; Assume the 0-th structure is the PLUGMAP data stname = tag_names(plugdat, /structure_name) if (stname NE 'PLUGMAPOBJ') then $ message, 'Invalid plugmap file' iobject = where(plugdat.holetype EQ 'OBJECT', nobject) pobj = plugdat[iobject] warp_initial_guess, pobj, gfib, k_x, k_y degree = (size(k_x))[1] - 1 ; IPLUG[i] will contain which plugmap-object number (0-639) corresponds to ; the i-th element of GFIB. iplug = lonarr(nfib) rmin = 4.5 ; Warp CCD positions to the plugmap coordinate system warp_apply, gfib.x, gfib.y, k_x, k_y, xwarp, ywarp if (keyword_set(waittime)) then begin plot, pobj.xfocal, pobj.yfocal, psym=6, symsize=1.5, $ xrange=[-350,350], yrange=[-350,350], xstyle=1, ystyle=1, $ title = 'INITIAL GUESS FOR PLATE SCALE + OFFSET', $ xtitle='xFocal [mm]', ytitle='yFocal [mm]', charsize=1.5 oplot, xwarp, ywarp, psym=1, symsize=1 oplot, [-310], [320], psym=6, symsize=1.5 xyouts, -310, 320, ' Plugmap position', charsize=1.5 oplot, [-310], [290], psym=1, symsize=1 xyouts, -310, 290, ' Scan event', charsize=1.5 wait, waittime endif iiter = 1 nmatch = 0 while (iiter LE maxiter AND nmatch LT nfib) do begin print print, format='("ITERATION ", i5, ":")', iiter ; Choose closest plugmap hole for each fiber for i=0, nfib-1 do begin dist = sqrt((xwarp[i] - pobj.xfocal)^2 + (ywarp[i] - pobj.yfocal)^2) if (min(dist, indmin) LT rmin) then begin iplug[i] = indmin endif else begin iplug[i] = -1 endelse endfor ; Determine which fibers are assigned to the same plug map positions. ; Discard all such fibers from the list of matches. for i=0, nfib-1 do begin if (iplug[i] NE -1) then begin indx = where(iplug EQ iplug[i], ct) if (ct GT 1) then begin print, ct, ' fibers point to the same plug position ', iplug[i] iplug[indx] = -1 endif endif endfor im = where(iplug NE -1, nmatch) print, ' Matched ', nmatch, ' of ', nfib, ' fibers' if (nmatch LT 3) then begin message, 'Catastrophic failure - too few matches!' end ; Find polynomial solution to matched points polywarp, pobj[iplug[im]].xfocal, pobj[iplug[im]].yfocal, $ gfib[im].x, gfib[im].y, degree, k_x, k_y ; Warp CCD positions to the plugmap coordinate system warp_apply, gfib.x, gfib.y, k_x, k_y, xwarp, ywarp iiter = iiter + 1 endwhile dx = xwarp[im] - pobj[iplug[im]].xfocal ; residuals dy = ywarp[im] - pobj[iplug[im]].yfocal dsig = mean( sqrt(dx^2+dy^2) ) if (keyword_set(waittime)) then begin plot, dx, dy, psym=1, symsize=0.5, $ title='POSITION RESIDUALS '+scanfile, $ xtitle='Delta X [mm]', ytitle='Delta Y [mm]' xyouts, 0.9 * !x.crange[0] + 0.1 * !x.crange[1], $ 0.9 * !y.crange[0] + 0.1 * !y.crange[1], $ 'RMS = '+string(dsig, format='(f6.3)')+' [mm]' endif ;--------------------------------------------------------------------------- ; Fill in the appropriate entries in the plugmap file ;--------------------------------------------------------------------------- ; Fill in the FIBS structure with the good fiber data fibs[where(fibs.good)] = gfib ; Fill in the PLUGDAT structure with SPECTROGRAPHID, FIBERID, THROUGHPUT ; PLUGDAT[IOBJECT] is the 640-element structure array of OBJECT plug data. ; GFIB contains the data for all the "good" fibers. plugdat[iobject].spectrographid = -1 ; Value for "bad" fibers plugdat[iobject].fiberid = -1 ; Value for "bad" fibers plugdat[iobject].throughput = 0 ; Value for "bad" fibers indx = where(iplug GE 0) ; Indices for assigned fibers in GFIB jndx = iobject[iplug[indx]] ; Corresponding indices in PLUGDAT plugdat[jndx].spectrographid = gfib[indx].spectroID plugdat[jndx].fiberid = $ (gfib[indx].spectroID - 1) * 320 + gfib[indx].num plugdat[jndx].throughput = gfib[indx].flux pobj = plugdat[iobject] ; Normalize the throughput measurements by fitting a quadratic function ; with radius from the plate center, then dividing by that function. ; Then rescale by the value of that function at radius=0. radius = sqrt(plugdat[jndx].xfocal^2 + plugdat[jndx].yfocal^2) thru = plugdat[jndx].throughput coeff = poly_fit(radius, thru, 2, thrufit) plugdat[jndx].throughput = plugdat[jndx].throughput * coeff[0] / thrufit if (keyword_set(waittime)) then begin !p.multi=[0,1,2] plot, radius, thru, psym=1, xrange=[0,350], xstyle=1, $ xtitle='Radius [mm]', ytitle='Throughput', $ title='THROUGHPUT AS FUNCTION OF RADIUS' xfit = findgen(350) yfit = coeff[0] + coeff[1] * xfit + coeff[2] * xfit^2 djs_oplot, xfit, yfit, color='red' plot, radius, plugdat[jndx].throughput, psym=1, $ xrange=[0,350], xstyle=1, $ xtitle='Radius [mm]', ytitle='Throughput', $ title='THROUGHPUT SCALED BY RADIAL CORRECTION' djs_oplot, !x.crange, [coeff[0],coeff[0]], color='red' !p.multi=[0,1,1] wait, waittime endif if (keyword_set(waittime)) then begin plot_unplugged, plugdat, $ title='UNPLUGGED HOLES '+scanfile wait, waittime endif plotfile1 = outdir+'/Unplugged-'+rootname+'.ps' plot_unplugged, plugdat, plotfile=plotfile1, $ title='UNPLUGGED HOLES '+scanfile add_hdr_line, hdr, 'fscanVersion', fscanVersion add_hdr_line, hdr, 'fmapVersion', fmapVersion add_hdr_line, hdr, 'fscanId', fscanId add_hdr_line, hdr, 'fscanMJD', fscanMJD add_hdr_line, hdr, 'fscanMode', fscanMode add_hdr_line, hdr, 'fscanSpeed', fscanSpeed add_hdr_line, hdr, 'pluggers', pluggers add_hdr_line, hdr, 'cartridgeId', cartridgeId add_hdr_line, hdr, ' ', '' outfile1 = outdir+'/plPlugMapM-'+rootname+'.par' yanny_write, outfile1, ptr_new(plugdat), $ hdr=hdr, enums=enums, structs=structs ; mwrfits, plugdat, 'plPlugMapM-'+rootname+'.fits' ;--------------------------------------------------------------------------- ; STEP 5: Print summary ;--------------------------------------------------------------------------- igood = where(pobj.fiberid GT 0, ngood) ibad = where(pobj.fiberid LE 0, nbad) print print, 'EVILMAP SUMMARY ' + systime() print print, ' Fibers (total): ', N_elements(pobj) print, ' Fibers (matched): ', ngood print print, ' Bad fiber list: ' if (nbad EQ 0) then begin print, ' ' endif else begin for i=1, N_elements(pobj)-1 do $ if ((where(pobj.fiberid EQ i))[0] EQ -1) then $ print, ' Fiber', string(i, format='(I4)') endelse print print, ' RMS agreement (mm): ', dsig ; Spawn Unix tasks to display and print the final plot. ; Copy that plot and ; the plPlugMapM file to the /data/spectro/plugMapM ; directory on sdsshost. print print, 'FOR THE PLATE DATABASE:' print, ' Plate =', platenum print, ' MJD =', fscanMJD print, ' mapId =', fscanId print, ' Cartridge =', cartridgeId if (keyword_set(copydir)) then begin print print, 'Copying files to ', copydir, '... (hangs if sdsshost is down)' spawn, 'gs ' + plotfile1 + ' &' spawn, 'lpr ' + plotfile1 + ' &' spawn, 'scp1 ' + plotfile1 + ' ' + copydir spawn, 'scp1 ' + outfile1 + ' ' + copydir spawn, 'scp1 ' + scandir+scanfile + ' ' + copydir ; spawn, 'cp ' + plotfile1 + ' /data/spectro/plugMapM' ; spawn, 'cp ' + outfile1 + ' /data/spectro/plugMapM' ; spawn, 'cp ' + scandir+scanfile + ' /data/spectro/plugMapM' print, 'Done.' endif return end ;------------------------------------------------------------------------------