[375] | 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 |
---|