; Messin' with H-HAS-MC-3.010
; Messin' with H-HAS-MC-3.010
;
;
; This is file: ~dd/idl/xrcf/acis_mc_cc_play.pro
;
; This file contains commands that could be executed at the IDL
; command line to play with the ACIS MC data.
;
; Starting from scratch it can be run by doing:
; IDL> astrolib
; IDL> @acis_mc_cc_play
;
; Note: The procedures pre_print_sqr(etc.) and post_print are used to create
; .ps files instead of screen display. They are in ~dd/idl/useful/
; The pre_print, post_print lines can be commented out and STOP command
; added after desired plot for interactive use.
;Reading-in and calibrating X,Y, and E
; Where is the FITS data file?
; At SAO this file is in: /data/pipes/p/113/H-HAS-MC-3.010/i001/acil1/out_R1U9.4/
; Thanks to Anastasia for processing it!
file_name = '/spectra/d6/970423/acis_hhasmc3010i001n000_evt.fits'
; read it in with ASTROLIB routine
; This line can be commented out for repeat playings with same data
; as all variables here are still at command line.
ev = mrdfits(file_name, 2)
; Show all the TAGS of ev:
print, 'TAGs of ev (FITS columns) : ', tag_names(ev)
; The output TAGs are:
; TIME EXPNO CCD_ID CHIPX CHIPY TDETX TDETY DETX DETY X Y PHA PI GRADE FLTGRADE
; Form an ACIS PHA histogram of all events
; Looks like pi is in units of eV
pha_hist = histogram(ev.pi,MIN=0,MAX=15000,BINSIZE=10)
pha_es = 10.0*(1.+indgen(n_elements(pha_hist))) ; eV
pre_print_sqr
plot_oo, pha_es/1000.0, pha_hist, PSYM=10, $
XRANGE=[0.1, 10.], YRANGE=[1.,1.E4], $
TITLE = 'ACIS-S spectrum of all events, H-HAS-MC-3.010', $
XTITLE = 'ACIS Energy (keV) [~pi/1000]', $
YTITLE='counts per 10 eV'
; Label the Al, and W peaks
xyouts, 1.48, 2000.0, 'Al-K'
xyouts, 1.78, 1000.0, 'W-M'
post_print, 'acis_mc_pha.ps'
; Turn the X,Y,and E values into "real" mm,mm,keV units
; with (X,Y) = (0.,0.) at zero-order
; E : looks like 1 eV per unit
acis_e = ev.pi/1000.0
; X,Y : put in the offsets to get zero-order at zero
; These offsets were derrived iteratively by looking at the
; following plots and playing with the data
acis_x = ev.X - 4094.0
acis_y = ev.Y - 4677.0
; Blur them by a "pixel" width
; Discrete (integer) values of X,Y will cause all sorts of problems
; with spectral binning and event plotting
acis_x = acis_x + RANDOMU(SEED,n_elements(acis_x)) - 0.5
acis_y = acis_y + RANDOMU(SEED,n_elements(acis_y)) - 0.5
; Plot the events in physical space
; Zero in on zero-order... with cross hairs...
; Adjust the acis_y correction above to get events on the line.
pre_print_sqr
plot, acis_x, acis_y, PSYM=3, XRANGE=[-500.,500.], YRANGE=[-500.,500.], $
TITLE='Close-up of Zero-order Region w/crosshairs', $
XTITLE='ACIS_X (ASC pixels)', YTITLE='ACIS_Y (ASC pixels)'
oplot, [-500.,500.],[0.,0.]
oplot, [0.,0.],[-500.,500.]
post_print, 'acis_mc_zeroin.ps'
; Calculate the floating time
; The TIME is a big number (1.E8 say) but is
; managable if the min value is subtracted off...
ftime = ev.time
ftime = ftime - MIN(ftime)
; Plot the acis_x event location vs ftime:
; shows dropped frames! of S3 chip after ftime ~140.
; Can also plot ACIS_X vs EXPNO: plot, ev.expno, acis_x, PSYM=3
pre_print_sqr
plot, ftime, acis_x, PSYM=3, TITLE='Events plotted as ACIS_X vs. Time', $
XTITLE='Time ( TIME - MIN(TIME) )', YTITLE='ACIS_X w.r.t 0-order (ASC pixels)'
post_print, 'acis_mc_XvsTime.ps'
; Look at Plus and Minus all-order spectra
; For this continous-clocking data use these plots to find zero-order in X
; by getting the Al-K lines to overlap (by adjusting the
; acis_x offset constant above.)
pre_print_sqr
!p.multi=[0,1,2] ; put two plots on one page
x_range = [0.,max(ABS(acis_x))]
;;x_range = [1450.,1490.] <-- uncomment to zero-in to Al-K lines
; Plus side...
x_hist = histogram(acis_x,MIN=0.,MAX=3000.,BINSIZE=1.)
y_range = [1., max(x_hist(x_range(0):(x_range(1) < n_elements(x_hist)-1)))]
plot_io, x_hist, XRANGE=x_range,YRANGE = y_range, PSYM=10, $
TITLE='Plus all-order events', XTITLE='ACIS_X w.r.t. zero-order', $
YTITLE='events per ASC pixel'
; indicate the Al-K line
xyouts, 1465., 300., 'Al-K'
; Negative side...
x_hist = histogram(-1.*acis_x,MIN=0.,MAX=3000.,BINSIZE=1.)
plot_io, x_hist, XRANGE=x_range, YRANGE=y_range, PSYM=10, LINESTYLE=2, $
TITLE='Minus all-order events', XTITLE='ACIS_X w.r.t. zero-order', $
YTITLE='events per ASC pixel'
xyouts, 1465., 300., 'Al-K'
!p.multi=0
post_print, 'acis_mc_pmspec.ps'
; Zooming in with x_range above yields a central channel of HEG Al line at 1465.
; The equivalent mm per pixel is defined to have this be 1.4867 keV.
; The COS(5 degrees) is thus included in this pixel value...
mmperpix = 8782.8*((12.3985/1.4867)/2000.95)/1465.0 ; dispersed um / pixel
print, 'mm per "dispersed pixel" is :', mmperpix
; which comes out to 25.0 - 0.0013 um. FYI the Cu MC data gave a
; value of mmperpix = 25.0047 um per ASC pixel
; Now put the acis event locations into dispersed mm s :
acis_x = mmperpix * acis_x
acis_y = mmperpix * acis_y
;Extracting the dispersed spectrum
; Now we're all set to extract the spectrum
; Extraction parameters
width = 1.0 ; mm about dispersion axis
offset = 0. ; mm shifts mask up or down (useful to get off pile up
; in non-cont colocking data)
e_width = 1.0 ; fraction of E that dispersed_E and acis_E have to agree to
lam_bin = 0.05 ; bin size in A to bin spectrum
ftime_max = 135. ; use only the beginning data that has no S3 drop outs
; HEG-specific parameters
angle = 0.0 ; 0.0 because it is continous clocking!
period = 2000.95 ; kind of by definition of mm above...
; Calculate the dispersion based energy assuming 1st order
disp_E = 12.3985/(period*(ABS(acis_X) > 0.1)/8782.8)
; Selection criteria for first-order photons
disp_select = where( acis_y LT offset+acis_x*SIN(!DTOR*angle)+width/2.0 AND $
acis_y GT offset+acis_x*SIN(!DTOR*angle)-width/2.0 AND $
acis_E LT disp_E*(1.+e_width/2.0) AND $
acis_E GT disp_E/(1.+e_width/2.0) AND $
ftime LT ftime_max)
; For the selected events plot Y vs. X
; and plot E vs. X
pre_print_sqr
plot, acis_x(disp_select), acis_y(disp_select), PSYM=3, $
TITLE='Y vs. X plot of selected First-order events', $
XTITLE='ACIS_X (mm)', YTITLE='ACIS_Y (mm)'
post_print, 'acis_mc_YvsX.ps'
pre_print_sqr
plot, acis_x(disp_select), acis_E(disp_select), $
PSYM=3, YRANGE=[0.,10.], $
TITLE='E vs. X plot of selected First-order events', $
XTITLE='ACIS_X (mm)', YTITLE='ACIS_E (keV)'
post_print, 'acis_mc_EvsX.ps'
; OK, now make the MC spectrum
; a histogram of these events' dispersed wavelengths
; (because wavelength bins are basically spatial bins so have const
; resolution...)
disp_lam_hist = histogram(12.3985/disp_E(disp_select), MIN=0., MAX=40.0, BINSIZE=lam_bin)
disp_lams = lam_bin*(1.+indgen(n_elements(disp_lam_hist)))
pre_print_sqr
plot_oo, 12.3985/disp_lams, disp_lam_hist, YRANGE = [10.,1.E4], PSYM=10, $
XRANGE = [1.0,10.], TITLE = 'HEG MC C Spectrum, H-HAS-MC-3.010', $
XTITLE = 'Dispersed Energy (keV)', YTITLE='counts/bin (bins ='+$
STRING(lam_bin,FORMAT='(F6.3)')+' A)'
post_print, 'acis_mc_spectrum.ps'
; Convert it to a Spectrum with units of counts/keV
;
; counts/keV = ( lambda^2 / hc ) * counts/A
;
spec_per_kev = ((disp_lams^2)/12.3985) * (disp_lam_hist/lam_bin)
spec_es = 12.3985/disp_lams
pre_print_sqr
plot_oo, spec_es, spec_per_kev, YRANGE = [100.,1.E5], PSYM=10, $
XRANGE = [1.0,10.], TITLE = 'HEG MC C Spectrum, H-HAS-MC-3.010', $
XTITLE = 'Dispersed Energy (keV)', YTITLE='counts/keV'
post_print, 'acis_mc_speckev.ps'
;Adding Predicted Curve
; OK, now make the MC spectrum with model
; curve using:
; - Dave H's HEG effective area that includes FI and BI, real gratings,
; and mirror (counts cm^2 / photon)
; - my Source simulation s/w for source spectrum (photons/cm^2/s/keV)
; read in the rdb version of Dave's HEG-ACIS-S Effective Area
heg_ea = rdb_read(!DDLOOKUPS+'/EA_HEG-AS.rdb')
; Read in the ACIS-phase CMDB and simulate the spectrum
cmdb_file_path = '~dd/idl/cmdb/Phase2/'
cmdb_file = 'req_run_2honly.cmdb'
cmdb_load
@cmdb_trw_integers
; Find where this measurement is
this_mc = where(cmdb_fields(*,c_id) EQ 'H-HAS-MC-3.010',nfound)
; simulate the spectrum
xss_sim, this_mc(0)+1
; Interpolate the smooth spectrum to Dave's grid
source_spec = INTERPOL(lx_tubespec, lx_es, heg_ea.energy)
; Multiply them and a time to get predicted counts / keV
predicted = heg_ea.area * source_spec * 135.0
pre_print_sqr
plot_oo, spec_es, spec_per_kev, YRANGE = [100.,1.E5], PSYM=10, $
XRANGE = [1.0,10.], $
TITLE = 'HEG-ACIS-S "Molecular Contamination Test" Spectrum and Prediction', $
XTITLE = 'Dispersed Energy (keV)', YTITLE='counts/keV'
oplot, heg_ea.energy, predicted, THICK=2
; Add some labels etc.
xyouts, 2.5, 5.E4, 'Solid curve - prediction, full HEG, t=135s'
xyouts, 2.5, 3.5E4,' Histogram - XRCF data, H-HAS-MC-3.010'
xyouts, 2.7, 1900.0, 'chip gap'
xyouts, 1.4, 7.E4, 'Al-K'
post_print, 'acis_mc_specwmod.ps'
;