source: trunk/doc/cookbook.txt@ 474

Last change on this file since 474 was 375, checked in by mar637, 20 years ago

ASAP cookbook

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.8 KB
Line 
1Quick and Dirty ATNF Spectra line Analysis Package Cookbook
2
3
4The following "cookbook" is intended to help asap testers get up and
5running. It is not intended to be a comprehensive user guide. If
6anything is not clear please get in contact with
7Chris.Phillips@csiro.au or Malte.Marquarding@csiro.au.
8
9***********
10* General *
11***********
12
13The asap user interface is the python scripting language (ipython to
14be exact). It currently uses an object oriented approach for the user
15interface. One job of the testers is to give feed back to how they
16like this or what they would prefer as the command line interface.
17
18The main object in use is a scantable which initially contains the
19scans from a set of observations as individual correlator
20cycles. There are currently some limitations on how this data can be
21averaged. For practical purposes at the moment each scan has to be
22copied from the initial scan table before averaging or other
23reduction.
24
25Asap comes with built in help (help functionname), see the examples
26further down in this document. ipython allows you to run shell
27built-ins (ls, pwd, cd, etc). It also allows tab completion of function
28names etc, so if you cannot remember full command name you can just
29type the first few characters and type tab (eg pol<tab>).
30
31Asap unses "zero-based" indexing. This means the first scan in a
32scantable is scan 0. The first IF, polarisation and beam would have
33index [0,0,0].
34
35****************
36* Running asap *
37****************
38
39Currently in the testing phase, asap only runs on Linux at the ATNF.
40
41Ssh onto brage or vulcan at Epping.
42
43 > source /nfs/aips++/weekly/aipsinit.csh # Tempory measure
44 > asap
45
46This will start up the asap ipython interpretor.
47
48****************
49* Reading data *
50****************
51
52To read data from disk, a reader object needs to be created then the
53scans copied into a scan table. This scan table will contain the
54averaged (!) correlator cycles for every scan in the rpfits file.
55
56 ASAP> scans = scantable('2003-03-16_082048_t0002.rpf')
57
58The scantable contains all the data from the rpfits file. To reduce
59the data you need to separate the individual scans. Before you do this
60is 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
65This will set the rest frequency and tell asap to use LSR velocity
66reference.
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
89If you have observed the source with multiple source/reference cycles you
90will 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
110To make a baseline fit, you must first create a mask of channels to
111use in the baseline fit. Remember that the default is the range to
112exclude, not include.
113
114 ASAP> msk = scan.create_mask([1600,1700],[1900,2000])
115 ASAP> scan2 = poly_baseline(scan, msk, 1)
116
117This will fit and subtract a first order polynomial to the masked region
118
119*****************************
120* Average the polarisations *
121*****************************
122
123If you are just interested in highest SNR for total intensity you
124want to average the parallel polarisations together.
125
126 ASAP> iav = average_pol(s)
127
128
129*********************
130* Plotting the data *
131*********************
132
133You can plot the data at any time.
134 ASAP> plotter.plot(s)
135
136This will bring up a plot window (a separate one for each scantable at
137the moment).
138
139The 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
169Currently multicomponent Gaussian function is available. This is done
170by creating a fitting object, setting up the fit and actually fitting
171the 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
180If you want to confine the fitting to a smaller range (e.g. to avoid
181band 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
191If 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
207Restfrequency, preferred velocity frame and active units need to be set
208for 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
249ASAP> 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
276A asap session has the potential to create a lot of scantable
277objects. To list the current objects and remove unneeded ones use the
278list_scan and del functions. Overwriting deletes the old one automatically
279
280 ASAP> list_scan
281----> list_scans()
282The 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
294All asap functions have built-in documentation. To get a summary of
295all functions use
296
297 ASAP> commands
298
299To get help on individual tasks use the help command. For the Maths
300functions do this as:
301
302 ASAP> help quotient
303
304For 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
321commands
322scans = scantable('2001-10-07_0337.rpf')
323# data is auto-averaged
324print scans
325scans.set_restfreqs([90.663543], 'GHz')
326scans.set_freqframe('LSRK')
327r1 = scans.get_scan(0)
328s1 = scans.get_scan(1)
329r2 = scans.get_scan(2)
330s2 = scans.get_scan(3)
331q1 = quotient(s1,r1)
332q2 = quotient(s2,r2)
333del s1,s2
334del r1,r2
335plotter.plot(q1)
336msk = q1.create_mask([100,500],[620,1024])
337help scantable.average_time
338av = average_time(q1, q2, mask=msk)
339plotter.plot(av)
340# reuse the old mask as data range hasn't changed
341s = poly_baseline(av,msk,1)
342plotter.plot(s)
343ss = average_pol(s)
344plotter.plot(ss)
345ss.save('asap.test')
346
347f = fitter()
348f.set_function(gauss=1)
349f.set_scan(ss)
350f.fit()
351f.plot()
352f.get_parameters()
353^d
Note: See TracBrowser for help on using the repository browser.