;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO rd_eval_fef, fn, f, p, npts, nsig f=mrdfits(fn,1,hf, /silent) p=eval_fef(f[fix(findgen(npts)/npts*n_elements(f))],hf,nsig) END ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FUNCTION eval_fef, ti, hdr, nsig ; quick-and-dirty fef evaluation (many ad hoc assumptions about the files!) ; ti : one or more rows of an fef table structure ; hdr : fef fits file header ; result is a structure containing tags x and y. ; if input is an array (more than one row), output is 2D array ; ; evaluate fef function, which is assumed to contain the structure: ; Troz> help,/st,f1 ; ** Structure <8257e4c>, 34 tags, length=136, refs=1: ; ENERGY FLOAT 0.0500000 ; B1_XLOW FLOAT 1.00000 ; B1_XHI FLOAT 4096.00 ; B1_AMPL FLOAT 1.00000 ; G1_FWHM FLOAT 0.00000 ; G1_POS FLOAT -106.610 ; G1_AMPL FLOAT 0.00000 ; G2_FWHM FLOAT 0.00000 ; G2_POS FLOAT -47.9500 ; G2_AMPL FLOAT 0.00000 ; G3_FWHM FLOAT 0.00000 ; G3_POS FLOAT 358.910 ; G3_AMPL FLOAT 0.00000 ; G4_FWHM FLOAT 0.00000 ; G4_POS FLOAT -48.4700 ; G4_AMPL FLOAT 0.00000 ; G5_FWHM FLOAT 0.00000 ; G5_POS FLOAT -3.97000 ; G5_AMPL FLOAT 0.00000 ; G6_FWHM FLOAT 0.00000 ; G6_POS FLOAT -6.45000 ; G6_AMPL FLOAT 0.00000 ; G7_FWHM FLOAT 22.0665 ; G7_POS FLOAT -2.25000 ; G7_AMPL FLOAT 901.850 ; G8_FWHM FLOAT 0.00000 ; G8_POS FLOAT 272.600 ; G8_AMPL FLOAT 0.00000 ; G9_FWHM FLOAT 11.5150 ; G9_POS FLOAT 6.49000 ; G9_AMPL FLOAT 14288.0 ; G10_FWHM FLOAT 27.5420 ; G10_POS FLOAT 16.8700 ; G10_AMPL FLOAT 74104.6 ; i.e., a clipping box and 10 gaussians, at one energy ; (the assumption could be relaxed, but the the header has to be parsed.) ; PROBLEM: the structure above is only for FI CCDs. ; BI has 2 gaussians, no box. Hacked below. np = n_params() IF np NE 3 THEN BEGIN message, /inf, 'Wrong number of parameters' message, /inf, 'Usage: result=eval_fef(fef_struct, fef_hdr, nsig )' return, -1 ENDIF pmin = sxpar(hdr, 'flmin2') ; assumes pha is in keyword index 2 pmax = sxpar(hdr, 'flmax2') nrows = n_elements(ti) xx = findgen(pmax-pmin) + pmin yy = fltarr(nrows, n_elements(xx)) ee = fltarr(nrows) ; energy vector, if multiple rows. frac = fltarr(nrows) elo = fltarr(nrows) ehi = fltarr(nrows) tags = tag_names(ti) nt = n_elements(tags) IF nt GT 7 THEN BEGIN ; guess FI ngauss = (nt-4) / 3 ENDIF ELSE BEGIN ; guess BI ngauss = (nt-1) / 3 ENDELSE ksig = 1./(2.*sqrt(2.*alog(2.))) FOR ir = 0, nrows-1 DO BEGIN ee[ir] = ti[ir].energy FOR ig = 0, ngauss-1 DO BEGIN ; loop over gaussians and sum them IF ngauss GT 2 THEN BEGIN ; guess FI gidx = ig*3 + 4 ; ad hoc; ; could regex parse tag_names(ti) to ; look for '_FWHM', '_AMPL', '_POS' ; and index them. Instead assume 0-3 ; are box, then gaussians come as ; FWHM, POS, AMPL,... ENDIF ELSE BEGIN ; guess BI gidx = ig*3 + 1 ENDELSE ifwhm = gidx ipos = gidx+1 iampl = gidx+2 amp = ti[ir].(iampl) pos = ti[ir].(ipos) sig = ti[ir].(ifwhm)*ksig yy[ir, *]= yy[ir, *] + amp * exp(-(xx-pos)^2/(2.*sig^2)) ; message, /inf, 'debug+++++++++++++++++++++++++' ; print, ir, ig, ifwhm, ipos, iampl, sig/ksig, pos, amp ; message, /inf, 'debug-------------------------' ENDFOR ; endfor ig (gaussians) IF ngauss GT 2 THEN BEGIN ; guess FI limits = where(outside(xx, ti[ir].b1_xlow, ti[ir].b1_xhi), nout) IF nout GT 0 THEN yy[ir, limits] = 0.0 ENDIF yy[ir, *] = yy[ir, *]/total(yy[ir, *]) ; normalize ; find peak; assumes g10_fwhm is main ; peak; IF ngauss GT 2 THEN BEGIN ; guess FI dpha = nsig * ti[ir].g10_fwhm * ksig ENDIF ELSE BEGIN ; guess BI dpha = nsig * ti[ir].g1_fwhm * ksig ENDELSE ip = (where(yy[ir, *] EQ max(yy[ir, *])))[0] ; find peak iplo = 0 > min(where(xx GE (xx[ip] - dpha) ) ) ; peak-n*sigma iphi = min(where(xx GE (xx[ip] + dpha) ) ) ; peak+n*sigma frac[ir] = total(yy[ir, iplo:iphi]) ; fractional response in +-n*sig ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;+++++++ ; PROVISIONAL - to test chip_id=4, cell x00,y00 g0 = 0.0 ; adu g1 = 0.26978 ; adu/eV; from FEF 00,00 chip_id=4 ; acis4a_x00_y00_FP-110_D1999-09-16fef_phaN0002.fits.gz ; agrees w/ gain table acisD1999-09-16gainN0004.fits.gz elo[ir] = (xx[ip] - dpha)/g1/1000. ehi[ir] = (xx[ip] + dpha)/g1/1000. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;-------- ENDFOR ; endfor irow (table rows) result = {e:ee, x:xx, y:yy, frac:frac, elo:elo, ehi:ehi} return, result END ; e.g. usage: ; nn=256 ; f=mrdfits('acis4a_x00_y00_FP-110_D1999-09-16fef_phaN0002.fits.gz',1,hf) ; p=eval_fef(f[fix(findgen(nn)/nn*1000)],hf,3) ; plot,p.x,p.y[nn-1,*],xra=[1,4000],yra=[1.e-7,4.e-1],/yty,/nod,/xst;,/xty ; write_gif,'tst.gif',tvrd(),/mult ; for ii=0,nn-1 do begin oplot,p.x,p.y[ii,*],col=(ii*(210./nn)+40) &$ ; xyouts, 0.1,0.9, string(p.e[ii],format='(f6.2)'),/norm &$ ; xyouts, 0.2,0.9, string(p.frac[ii],format='(f6.2)'),/norm &$ ; wait,0.04 &$ ; write_gif,'tst.gif',tvrd(),/mult&$ ; oplot,p.x,p.y[ii,*],col=0 &$ ; xyouts, 0.1,0.9, string(p.e[ii],format='(f6.2)'),/norm,col=0 &$ ; xyouts, 0.2,0.9, string(p.frac[ii],format='(f6.2)'),/norm,col=0 &$ ; endfor ; write_gif,'tst.gif',tvrd(),/close ; more example : compare osip to dph: ; ; f=mrdfits('acis4a_x00_y00_FP-110_D1999-09-16fef_phaN0002.fits.gz',1,hf) ; g=mrdfits('acisD1999-09-16gainN0004.fits.gz',1,hg) ; o=mrdfits('acisD1999-09-16osipN0001.fits.gz',1,ho) ; p=eval_fef(f[fix(findgen(nn)/nn*1000)],hf,3) ; pa=eval_fef(f[fix(findgen(nn)/nn*1000)],hf,2) ; plot,!hc/p.e,p.ehi/p.e,xra=[0,30],xtit='wavelength',ytit='order',tit='ccdid=4,cell 00,00; 3sigma',yra=[0.6,1.4] ; oplot,!hc/p.e,p.elo/p.e ; ll=where( (o.ccd_id eq 4) and (o.chipx_min lt 33) and (o.chipy_min lt 33)) ; oplot,!hc/o[ll].energy*1000,o[ll].energ_hi/o[ll].energy,psy=4,thick=2,col=199 ; oplot,!hc/o[ll].energy*1000,o[ll].energ_lo/o[ll].energy,psy=4,thick=2,col=199 ; oplot,!hc/pa.e,pa.elo/pa.e,col=55 ; oplot,!hc/pa.e,pa.ehi/pa.e,col=55 ; legend,['dph 3*sig','dph 2*sig','gea osip 2*sig'],psy=[0,0,4],col=[250,55,199] ; wgif,0,'Cmp_dph_osip_c4x00y00.gif' ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO rmf_movie_gif, p, fout nn = n_elements(p.y[*, 0]) plot,p.x,p.y[nn-1,*],xra=[1,3000],yra=[1.e-7,4.e-1],/yty,/nod,/xst, /yst do_gif = 0 IF n_elements(fout) NE 0 THEN do_gif = 1 plot,p.x,p.y[0,*],xra=[1,3000],yra=[1.e-7,4.e-1],/yty,/xst, /yst IF do_gif THEN write_gif, fout, tvrd(), /multiple FOR ii=1,nn-1 DO BEGIN plot,p.x,p.y[ii,*],xra=[1,3000],yra=[1.e-7,4.e-1],/yty,/xst, /yst xyouts, 0.1,0.9, string(p.e[ii],format='(f6.2)'),/norm xyouts, 0.2,0.9, string(p.frac[ii],format='(f6.2)'),/norm empty ; will this make display smoother? IF do_gif THEN write_gif, fout, tvrd(), /multiple ENDFOR IF do_gif THEN BEGIN write_gif, fout, /close ENDIF END