1 | Quick and Dirty ATNF Spectra line Analysis Package Cookbook |
---|
2 | |
---|
3 | |
---|
4 | The following "cookbook" is intended to help asap testers get up and |
---|
5 | running. It is not intended to be a comprehensive user guide. If |
---|
6 | anything is not clear please get in contact with |
---|
7 | Chris.Phillips@csiro.au or Malte.Marquarding@csiro.au. |
---|
8 | |
---|
9 | *********** |
---|
10 | * General * |
---|
11 | *********** |
---|
12 | |
---|
13 | The asap user interface is the python scripting language (ipython to |
---|
14 | be exact). It currently uses an object oriented approach for the user |
---|
15 | interface. One job of the testers is to give feed back to how they |
---|
16 | like this or what they would prefer as the command line interface. |
---|
17 | |
---|
18 | The main object in use is a scantable which initially contains the |
---|
19 | scans from a set of observations as individual correlator |
---|
20 | cycles. There are currently some limitations on how this data can be |
---|
21 | averaged. For practical purposes at the moment each scan has to be |
---|
22 | copied from the initial scan table before averaging or other |
---|
23 | reduction. |
---|
24 | |
---|
25 | Asap comes with built in help (help functionname), see the examples |
---|
26 | further down in this document. ipython allows you to run shell |
---|
27 | built-ins (ls, pwd, cd, etc). It also allows tab completion of function |
---|
28 | names etc, so if you cannot remember full command name you can just |
---|
29 | type the first few characters and type tab (eg pol<tab>). |
---|
30 | |
---|
31 | Asap unses "zero-based" indexing. This means the first scan in a |
---|
32 | scantable is scan 0. The first IF, polarisation and beam would have |
---|
33 | index [0,0,0]. |
---|
34 | |
---|
35 | **************** |
---|
36 | * Running asap * |
---|
37 | **************** |
---|
38 | |
---|
39 | Currently in the testing phase, asap only runs on Linux at the ATNF. |
---|
40 | |
---|
41 | Ssh onto brage or vulcan at Epping. |
---|
42 | |
---|
43 | > source /nfs/aips++/weekly/aipsinit.csh # Tempory measure |
---|
44 | > asap |
---|
45 | |
---|
46 | This will start up the asap ipython interpretor. |
---|
47 | |
---|
48 | **************** |
---|
49 | * Reading data * |
---|
50 | **************** |
---|
51 | |
---|
52 | To read data from disk, a reader object needs to be created then the |
---|
53 | scans copied into a scan table. This scan table will contain the |
---|
54 | averaged (!) correlator cycles for every scan in the rpfits file. |
---|
55 | |
---|
56 | ASAP> scans = scantable('2003-03-16_082048_t0002.rpf') |
---|
57 | |
---|
58 | The scantable contains all the data from the rpfits file. To reduce |
---|
59 | the data you need to separate the individual scans. Before you do this |
---|
60 | is is advisable to set the restfrequency etc on all scans. |
---|
61 | |
---|
62 | ASAP> scans.set_restfreqs([1665.4018], 'MHz') |
---|
63 | ASAP> scans.set_freqframe('LSRK') |
---|
64 | |
---|
65 | This will set the rest frequency and tell asap to use LSR velocity |
---|
66 | reference. |
---|
67 | |
---|
68 | |
---|
69 | ******************************************** |
---|
70 | * Extract scans and make a quotient spectra * |
---|
71 | ******************************************** |
---|
72 | |
---|
73 | ASAP> q = scans.get_scan(0) |
---|
74 | ASAP> s = scans.get_scan(1) |
---|
75 | |
---|
76 | This copies the data for the first and second scan. |
---|
77 | |
---|
78 | ASAP> scan = quotient(s, q) |
---|
79 | |
---|
80 | "scan" now contains the quotient spectrum, which you might want to look at. |
---|
81 | The default display is to have polarisation colour stacked and scans across panels. |
---|
82 | |
---|
83 | ASAP> plotter.plot(scan) |
---|
84 | |
---|
85 | ******************************* |
---|
86 | * Time average separate scans * |
---|
87 | ******************************* |
---|
88 | |
---|
89 | If you have observed the source with multiple source/reference cycles you |
---|
90 | will want to average the quotient spectra together. |
---|
91 | **Note** - You have to have time_averaged the individual scans before you |
---|
92 | average different scans. But you had to time average the individual |
---|
93 | scans before you do the quotient spectra, so you should be OK. |
---|
94 | - When time averaging multiple scans, asap uses rms weighting. |
---|
95 | Generally you have to set a mask so it can calculate the rms. |
---|
96 | - The default setting for a mask it to specify the channel |
---|
97 | range you want to *include* in the fit. The "invert=True" option |
---|
98 | below lets to choose the range you want to exclude (rather than |
---|
99 | include). |
---|
100 | |
---|
101 | ASAP> msk = s1.create_mask([0,100],[450,550],[900,1024]) # Include these |
---|
102 | ASAP> msk = s1.create_mask([250,450],[550,750],invert=True) # or exclude these |
---|
103 | |
---|
104 | ASAP> av = average_time(s1, s2, s3, mask=msk) |
---|
105 | |
---|
106 | ******************* |
---|
107 | * Baseline fitting * |
---|
108 | ******************** |
---|
109 | |
---|
110 | To make a baseline fit, you must first create a mask of channels to |
---|
111 | use in the baseline fit. Remember that the default is the range to |
---|
112 | exclude, not include. |
---|
113 | |
---|
114 | ASAP> msk = scan.create_mask([1600,1700],[1900,2000]) |
---|
115 | ASAP> scan2 = poly_baseline(scan, msk, 1) |
---|
116 | |
---|
117 | This will fit and subtract a first order polynomial to the masked region |
---|
118 | |
---|
119 | ***************************** |
---|
120 | * Average the polarisations * |
---|
121 | ***************************** |
---|
122 | |
---|
123 | If you are just interested in highest SNR for total intensity you |
---|
124 | want to average the parallel polarisations together. |
---|
125 | |
---|
126 | ASAP> iav = average_pol(s) |
---|
127 | |
---|
128 | |
---|
129 | ********************* |
---|
130 | * Plotting the data * |
---|
131 | ********************* |
---|
132 | |
---|
133 | You can plot the data at any time. |
---|
134 | ASAP> plotter.plot(s) |
---|
135 | |
---|
136 | This will bring up a plot window (a separate one for each scantable at |
---|
137 | the moment). |
---|
138 | |
---|
139 | The plot window has 6 buttons on the bottom left hand side. |
---|
140 | |
---|
141 | - Home: This will unzoom the window to the original zoom factor. A current |
---|
142 | bug means that if you change the axis units (say from channel |
---|
143 | to velocity) strange things happen |
---|
144 | |
---|
145 | - Plot history (left ad right arrow). The left arrow sets the plot zoom |
---|
146 | to the previous value. The right arrow returns it. This allows you |
---|
147 | to zoom in on one feature then return the plot to how it was |
---|
148 | previously. |
---|
149 | |
---|
150 | - Pan: The Cross sets the cursor to pan, or scroll mode allowing you to shift |
---|
151 | the plot within the window. Useful when zoomed in on a feature. |
---|
152 | |
---|
153 | - Zoom (the letter with the magnifying glass) lets you draw a rectangle around |
---|
154 | a region of interest then zooms in on that region. Use the plot |
---|
155 | history to unzoom again. |
---|
156 | plotter.set_range(100,500) |
---|
157 | sets a permanent viewing window on the abcissa between in the range 100 to 500 in the current unit. |
---|
158 | plotter.set_range() |
---|
159 | resets this |
---|
160 | |
---|
161 | - Save (floppy disk). Save the plot as a postscript or .png file. |
---|
162 | plotter.save('myfile.ps') doers the same from the command line |
---|
163 | |
---|
164 | |
---|
165 | *********** |
---|
166 | * Fitting * |
---|
167 | *********** |
---|
168 | |
---|
169 | Currently multicomponent Gaussian function is available. This is done |
---|
170 | by creating a fitting object, setting up the fit and actually fitting |
---|
171 | the data. |
---|
172 | |
---|
173 | ASAP> f = fitter() |
---|
174 | ASAP> f.set_function(gauss=1) |
---|
175 | ASAP> f.set_scan(s) |
---|
176 | ASAP> f.fit() |
---|
177 | ASAP> f.plot() |
---|
178 | ASAP> f.get_parameters() |
---|
179 | |
---|
180 | If you want to confine the fitting to a smaller range (e.g. to avoid |
---|
181 | band edge effects or RFI |
---|
182 | |
---|
183 | ASAP> f = fitter() |
---|
184 | ASAP> f.set_function(gauss=2) |
---|
185 | ASAP> msk = s.create_mask([2000,2100]) |
---|
186 | ASAP> f.set_scan(s,msk) |
---|
187 | ASAP> f.fit() |
---|
188 | ASAP> f.plot() |
---|
189 | ASAP> f.get_parameters() |
---|
190 | |
---|
191 | If you want to set the initial paramter guesses and fix specific parameters: |
---|
192 | |
---|
193 | ASAP> f = fitter() |
---|
194 | ASAP> f.set_function(gauss=2) |
---|
195 | ASAP> f.set_scan(s,msk) |
---|
196 | ASAP> f.fit() |
---|
197 | ASAP> p = [1.2,160.,50.,0.7,280.,100] |
---|
198 | # order is: [peak0,centre0,FWHM0,peak1,centre1,FWHM1] |
---|
199 | ASAP> fxd = [0,1,0,0,1,0] # fix the centres |
---|
200 | ASAP> f.set_parameters(p,fxd) # set the parameters |
---|
201 | ASAP> f.fit() |
---|
202 | |
---|
203 | ********* |
---|
204 | * Units * |
---|
205 | ********* |
---|
206 | |
---|
207 | Restfrequency, preferred velocity frame and active units need to be set |
---|
208 | for plotting etc to give sensible velocity units. |
---|
209 | |
---|
210 | * Set the rest frequency of the transition (in MHz, GHz etc). |
---|
211 | |
---|
212 | ASAP> s.set_restfreqs([1666.4018], 'MHz') |
---|
213 | |
---|
214 | * Set the preferred velocity reference frame ('REST', 'TOPO', 'LSRD', |
---|
215 | 'LSRK', 'BARY', 'GEO', 'GALACTO', 'LGROUP', 'CMB') |
---|
216 | |
---|
217 | ASAP> s.set_freqframe('LSRK') |
---|
218 | |
---|
219 | * Set the units for plotting, mask selection etc ('km/s', 'channel', 'MHz') |
---|
220 | |
---|
221 | ASAP> s.set_unit('km/s') |
---|
222 | |
---|
223 | **NOTE** When the unit has been set to, say, km/s then this is the unit used |
---|
224 | for the create_mask function etc, ie give the mask range in km/s |
---|
225 | and not channels. |
---|
226 | |
---|
227 | *************************** |
---|
228 | * Miscellaneous reduction * |
---|
229 | *************************** |
---|
230 | |
---|
231 | * Scale a spectra by a fixed amount |
---|
232 | |
---|
233 | ASAP> ss = scale(s, 1.2) |
---|
234 | |
---|
235 | * Hanning smooth the spectrum |
---|
236 | |
---|
237 | ASAP> ss = hanning(s) |
---|
238 | |
---|
239 | * Bin average the spectrum and reduce the number of spectral points |
---|
240 | |
---|
241 | ASAP> ss = bin(s, 4) |
---|
242 | |
---|
243 | * Save a scan |
---|
244 | |
---|
245 | ASAP> s.save('scan.asap') |
---|
246 | |
---|
247 | * Read a saved scan from disk |
---|
248 | |
---|
249 | ASAP> s = scantable('scan.asap') |
---|
250 | |
---|
251 | * Execute a python script named 'myscript' |
---|
252 | |
---|
253 | ASAP> execfile('myscript') |
---|
254 | |
---|
255 | * To see a summary of a individual or set of scans |
---|
256 | |
---|
257 | ASAP> print s |
---|
258 | |
---|
259 | or alternatively: |
---|
260 | |
---|
261 | ASAP> s.summary() |
---|
262 | |
---|
263 | * Get all (reference) scans which match the pattern '*_R': |
---|
264 | |
---|
265 | ASAP> s = scans.get_scan('*_R') |
---|
266 | |
---|
267 | * Or just and individual: |
---|
268 | |
---|
269 | ASAP> s = scans.get_scan('232p621_S') |
---|
270 | |
---|
271 | |
---|
272 | ******************* |
---|
273 | * Object clean up * |
---|
274 | ******************* |
---|
275 | |
---|
276 | A asap session has the potential to create a lot of scantable |
---|
277 | objects. To list the current objects and remove unneeded ones use the |
---|
278 | list_scan and del functions. Overwriting deletes the old one automatically |
---|
279 | |
---|
280 | ASAP> list_scan |
---|
281 | ----> list_scans() |
---|
282 | The user created scantables are: |
---|
283 | ['av', 'q1', 's2', 'scans', 'q2', 's', 'ss', 'r1', 'r2', 's1'] |
---|
284 | |
---|
285 | ASAP> del s1 |
---|
286 | ASAP> del s2 |
---|
287 | ASAP> del q1 |
---|
288 | ASAP> del q2 |
---|
289 | |
---|
290 | **************** |
---|
291 | * Builtin help * |
---|
292 | **************** |
---|
293 | |
---|
294 | All asap functions have built-in documentation. To get a summary of |
---|
295 | all functions use |
---|
296 | |
---|
297 | ASAP> commands |
---|
298 | |
---|
299 | To get help on individual tasks use the help command. For the Maths |
---|
300 | functions do this as: |
---|
301 | |
---|
302 | ASAP> help quotient |
---|
303 | |
---|
304 | For the scantable or fitter functions use the syntax |
---|
305 | |
---|
306 | ASAP> help scantable.set_freqframe |
---|
307 | ASAP> help fitter.set_scan |
---|
308 | |
---|
309 | **NOTE** If you type "help" only, you get ipython help. To exit this mode |
---|
310 | you need to type ^-d (control-d). |
---|
311 | |
---|
312 | * Notes |
---|
313 | |
---|
314 | - Only does one IF/Pol/Beam. Need to change cursor |
---|
315 | |
---|
316 | |
---|
317 | ****************** |
---|
318 | * Worked example * |
---|
319 | ****************** |
---|
320 | |
---|
321 | commands |
---|
322 | scans = scantable('2001-10-07_0337.rpf') |
---|
323 | # data is auto-averaged |
---|
324 | print scans |
---|
325 | scans.set_restfreqs([90.663543], 'GHz') |
---|
326 | scans.set_freqframe('LSRK') |
---|
327 | r1 = scans.get_scan(0) |
---|
328 | s1 = scans.get_scan(1) |
---|
329 | r2 = scans.get_scan(2) |
---|
330 | s2 = scans.get_scan(3) |
---|
331 | q1 = quotient(s1,r1) |
---|
332 | q2 = quotient(s2,r2) |
---|
333 | del s1,s2 |
---|
334 | del r1,r2 |
---|
335 | plotter.plot(q1) |
---|
336 | msk = q1.create_mask([100,500],[620,1024]) |
---|
337 | help scantable.average_time |
---|
338 | av = average_time(q1, q2, mask=msk) |
---|
339 | plotter.plot(av) |
---|
340 | # reuse the old mask as data range hasn't changed |
---|
341 | s = poly_baseline(av,msk,1) |
---|
342 | plotter.plot(s) |
---|
343 | ss = average_pol(s) |
---|
344 | plotter.plot(ss) |
---|
345 | ss.save('asap.test') |
---|
346 | |
---|
347 | f = fitter() |
---|
348 | f.set_function(gauss=1) |
---|
349 | f.set_scan(ss) |
---|
350 | f.fit() |
---|
351 | f.plot() |
---|
352 | f.get_parameters() |
---|
353 | ^d |
---|