Ken Glotfelty # # Setup for CIAO (whatever your ailas happens to be) # source /soft/ciao/bin/ciao.csh -d /soft/ciao # # Local setup for easy cut-n-paste # set root = /data/ciao_demo/kjg ----------------- Tri Color Image ---------------- # # This demo shows some introductory varmm reading and s-lang accesing # of the data combined with some simple DM filtering and the new tool # dmimg2jpg. # cd $root/TriColorImage # # Display Cas-A image # ds9 acisf00214N002_evt2.fits # # or just ds9 img.fits (faster) # # # Want to talk about spectrum, see if there is a spatial dependence # on the energy of the detected events # # # Extract the spectrum. Note the use of short hand "clob" for # "clobber" and "+" for "yes" # dmextract acisf00214N002_evt2.fits"[bin pi=1:1024:1]" spectrum.fits clobb+ # # Now want to display the spectrum, let's goto chips ... # chips d=readfile("spectrum.fits") # # 'd' is now a variable structure that contains the data read from # the spectrum file. All forms of math functions are availalbe. C # programmers will recognize the foo.goo syntax for referencing # elements of a structure. Underneath lib is S-Lang (Thanks John Davis) # print(d) # # Plot up the data # clear plot x d.channels y d.counts c 1 symbol none c 1 step limits x 0 400 xlabel channels ylabel counts # # Okay, now, looking at the plot we see lots of structure in the # spectrum ... is there a spatial dependence? First, lets split the # spectrum into 3 parts ... low, med, and high energies and assign # them to red, green, and blue colors. # # By eye let's pick # # rd = where(d.channels < 120) # # "where" is an intrinsic S-Lang command that returns the index of # the array 'where' whatever the expression is is true. # Repeat for green and blue. gn = where(d.channels >= 120 & d.channels < 150) bu = where( d.channels >= 150) # # Okay, now lets add some color to the plot to see what we have # # redraw off # # Need to define temp variable to plot with. # xx1 = d.channels[rd] yy1 = d.counts[rd] # # We use the indexes in the 'rd' variable to select the data # Similiarly for the green and blue below. # xx2 = d.channels[gn] yy2 = d.counts[gn] xx3 = d.channels[bu] yy3 = d.counts[bu] # # Just breeze thru all this to save time # plot x xx1 y yy1 plot x xx2 y yy2 plot x xx3 y yy3 c 2 symbol none c 2 step c 2 red c 3 symbol none c 3 step c 3 green c 4 symbol none c 4 step c 4 blue limits x 0 400 D 2 LABEL 300.0 19000.0 "x < 120" D 2 LABEL red D 3 LABEL 300.0 18000.0 "120 < x < 150" D 3 LABEL green D 4 LABEL 300.0 17000.0 "150 < x" D 4 LABEL blue redraw # # Now can see the three color ranges == three energy ranges. Let's # see if there is a spatial dependence. Let's get out of chips. # quit # # Now, let's make 3 images using the same "channel" range. Note that # "channles" is the same as "pi" in the event files. # # # You can run the following, but they take a minute or so each, so you # might want to use the "canned" version. Basically explain the dm # binning and filtering syntax. You can try to switch the order of # binning and filtering...that's supposed to work now. # #dmcopy "acisf00214N002_evt2.fits[pi=0:119][bin x=3150:4503:2,y=3300:4650:2]" red.fits #dmcopy "acisf00214N002_evt2.fits[pi=120:149][bin x=3150:4503:2,y=3300:4650:2]" grn.fits #dmcopy "acisf00214N002_evt2.fits[pi=150:1024][bin x=3150:4503:2,y=3300:4650:2]" blu.fits # # Now let's compile all three images into a single 'true' color image. # We'll make a JPEG file # dmimg2jpg red.fits grn.fits blu.fits cas_a.jpg showgrid=no showaim- clobber+ # # And finially display. # xv cas_a.jpg ----------------------- Plotting orders ----------------------------- # # This shows a bit more advanced s-land/varmm usage and shows some # grating data # cd $root/Orders chips # # Let's read some grating data. This might take 30-45sec. # We are only reading in SOME of the data, and using DM filtering to # cut down the size (for speed) # evts=readfile("acism00459N000_evt2.fits[expno<1000][cols tg_m,tg_part,tg_lam,tg_r,tg_d,energy]") # # What do we have? # print(evts) # # And let's plot the data # clear redraw off plot x evts.tg_r y evts.energy xlabel "tg\_r" ylabel "energy" LABEL 0.25 19000 "all data" redraw # # Now, let's only select the non-background data. # # tg_part tells whether it is zero-th order, MEG, HEG, LEG, or # background # tg_m tells the grating order # energy is the event energy # # This command selects what we'll call "good" events. Note the use # of "abs(evts.tg_m)" to select +/- orders with the single command # along with the complicated expression for the "where" command. # good=where(evts.tg_part<5 & evts.energy<1e4 & abs(evts.tg_m) < 5 ) # # Select the data and plot only the good data # d1 = evts.tg_r[good] d2 = evts.energy[good] redraw off clear plot x d1 y d2 c 1 symbol point xlabel "tg\_r" ylabel "energy" LABEL 0.25 9500 "no background" redraw # # Okay, now lets see where the first, second, third and zero-th order # photons are. We'll color each differently from all the events. # clear redraw off plot x evts.tg_r y evts.energy c 1 symbol point c 1 symbol white limits y 0 1e4 xlabel "tg\_r" ylabel "energy" redraw # # Now, let's select +/- 1st (tg_m) order HEG (tg_part) photons # redraw off fo = where((evts.tg_part ==1) & (evts.energy<1e4) &( abs(evts.tg_m) == 1) ) fox = evts.tg_r[fo] foy = evts.energy[fo] plot x fox y foy c 2 symbol point c 2 symbol red redraw # # Repeat for second order, third order, and zeroth-order (do rapidly) # so = where((evts.tg_part==1) & (evts.energy<1e4) & (abs(evts.tg_m) == 2) ) sox = evts.tg_r[so] soy = evts.energy[so] plot x sox y soy c 3 symbol point c 3 symbol green to = where((evts.tg_part==1) & (evts.energy<1e4) & ( abs(evts.tg_m) == 3) ) tox = evts.tg_r[to] toy = evts.energy[to] plot x tox y toy c 4 symbol point c 4 symbol blue zo = where((evts.tg_part==0) & ( evts.energy<1e4) & (abs(evts.tg_m) == 0 )) zox = evts.tg_r[zo] zoy = evts.energy[zo] plot x zox y zoy c 5 symbol point c 5 symbol yellow D 5 LABEL 0.25 10000 "0th order" D 5 LABEL yellow D 2 LABEL 0.25 9500 "1st order HEG" D 2 LABEL red D 3 LABEL 0.25 9000 "2nd order HEG" D 3 LABEL green D 4 LABEL 0.25 8500 "3rd order HEG" D 4 LABEL blue redraw # #Done # quit ------------------dmgti/chips/s-lang/...----------------------------------------------- # # This shows how to make a GTI and see what limits were actually set. # cd $root/GTI prism acisf01843_000N001_mtl1.fits # # Select time and count_rate columns (first and last columns) # # Select Visualize and QuickLook Plot. # # # Now, let's say we want to make a good time interval GTI where the # count_rate is < 20 cnts/sec. We can run dmgti from the prism GUI # # # Goto "Analysis -> Tools -> Utilities -> dmgti" # # # Note that the file name has already been filed in! # # # In the "outfile" parameter, type "my_gti.fits" # In the "userlimit" parameter, type "count_rate<20" # # Click on "Expand/Contract" # Click on "Show hidden params" # Click on "Apply" # # In the "verbose" parameter, type "2" # # Click on "Run" # # Wait for the Status in the tool agent window to say "The # application has exited" and then click on "DISMISS" # # # You can view the file in prism if you wish by going to File -> Open # -> ... # # # Okay, now let's visualize what the limit means. A GTI is a special # datamodel filter, so we can use it to filter the mission time line # (MTL) file # # # Click on "Analysis -> Modeling/Fitting # This will bring up a new term w/ sherpa running. # # Note that all power of chips and s-lang are available to sherpa due # to "fall thru" nature of their parsers. # # # Okay, let's load the data, first the unfiltered data # orig=readfile("acisf01843_000N001_mtl1.fits[cols time,count_rate]") clear plot x orig.time y orig.count_rate xlabel "time" ylabel "count\_rate" # # Now, what was the "good" data, use dm syntax to apply the GTI # filt=readfile("acisf01843_000N001_mtl1.fits[@my_gti.fits][cols time,count_rate]") clear plt x filt.time y filt.count_rate xlabel "time" ylabel "count\_rate" # # Now, what about the "bad" data, can we see that too? Yup. We can # use the "exclude" command to negate the filter. You can also try # "!" to see if that works. There are subtle difference between the # two. # bad=readfile("acisf01843_000N001_mtl1.fits[exclude @my_gti.fits][cols time,count_rate]") clear plt x bad.time y bad.count_rate xlabel "time" ylabel "count\_rate" # # Now, let's play with colors and color code the good and bad # redraw off clear plt x filt.time y filt.count_rate c 1 symbol green plt x bad.time y bad.count_rate c 2 symbol red xlabel "time" ylabel "count\_rate" D 1 LABEL 84279000 55 "y > 20" D 1 LABEL red D 2 LABEL 84279000 52 "y < 20" D 2 LABEL green redraw # # It's now easy to tell which data was good and which were bad. # quit # # Note that much more complex expression are allowed by dmgti... can # play with some. ahelp "dmtcalc_expressions" to learn more # # # Exit from prism going to "File -> Exit" # ---------------dmcontour/ region-filtering-------------------------------- # # Now for some "powerful" region filtering # # cd $root/Contour # # Let's make a binned image (just use 'canned' file...quicker) # # #dmcopy "acisf00214N002_evt2.fits[bin x=3150:4503:4,y=3300:4650:4]" img_4.fits # ds9 img_4.fits # # Let's say we want to smooth the image, we can # # # #csmooth infile='img_4.fits' sclmap='' outfile='smooth.fits' outsigfile='sig.fits' \ #outsclfile='scl.fits' conmeth='fft' conkerneltype='gauss' sigmin='4' \ #sigmax='5' sclmin='INDEF' sclmax='INDEF' stepzero='0.01' sclmode='compute' \ #bkgmode='local' bkgmap='' bkgerr='' clobber='yes' verbose='0' kernel='default' mode='ql' # # # This takes about 2 hrs ... don't do it. Use the 'canned' results. # ds9 smooth.fits & # # Now, let's say that we trust that pixels with more than 3 counts # are "source" events ... and we want to make a filter. It'd be # really complicated, but we have a tool that makes it easy # dmcontour infile='smooth.fits' levels='3' outfile='reg.fits' \ verbose='0' clobber='yes' kernel='fits' mode='ql' # # What do we have ... # dmlist /tmp/reg.fits data,array | less # # This shows we have a 765 sided polygon ... pretty neat, huh? # # # You can view the region file in ds9 directly # ds9 img_4.fits # # Then goto Region -> region -> load # Select "reg.fits" # # # # Fine, so lets use the file as a filter for the image -- # dmcopy "smooth.fits[(x,y)=region(reg.fits)]" filtdata.fits clobber+ # # Note here the use of "(x,y)" to indicate SKY coordinate. If these # were omitted, it would filter in logical coords and the image would # be blank. # ds9 filtdata.fits # # We can use the filter on other images too. How about the image # with the smoothing scale (another output of csmooth) # dmcopy "scl.fits[(x,y)=region(reg.fits)]" filtdata2.fits clobber+ ds9 filtdata2.fits # # Finially, what about using on the original event list? Let's make # a spectrum from the BACKGROUND region ... we'll use the nifty # "exclude" syntax again ... # # #dmextract "acisf00214N002_evt2.fits[exclude sky=region(/tmp/reg.fits)][bin pi=1:1023:1]" \ #spectra2.fits clobber+ # # (takes a long time, use "canned" file) # # # And display results # prism spectra2.fits ------------------------------------------------