| [1636] | 1 | ==========================================================================
|
|---|
| 2 | Tutorial 3 - Spectral-line gaussian component fitting and rms noise levels
|
|---|
| 3 | ==========================================================================
|
|---|
| 4 |
|
|---|
| 5 | .. sectionauthor:: Shari Breen
|
|---|
| 6 |
|
|---|
| 7 | Files
|
|---|
| 8 | -----
|
|---|
| 9 |
|
|---|
| 10 | * mon.asap: Data file
|
|---|
| 11 |
|
|---|
| 12 | Data Log
|
|---|
| 13 | --------
|
|---|
| 14 |
|
|---|
| 15 | * ON-OFF Position switching mode with Hobart
|
|---|
| 16 |
|
|---|
| 17 | * 4 scans (2 OFF and 2 ON)
|
|---|
| 18 |
|
|---|
| 19 | * Dual pol, 1 IF
|
|---|
| 20 |
|
|---|
| 21 | Instructions
|
|---|
| 22 | ------------
|
|---|
| 23 |
|
|---|
| 24 | * Work through the list of commands given in the text file to calibrate the data
|
|---|
| 25 | and fit gaussian components to the emission. The commands should be typed
|
|---|
| 26 | into ASAP line-by-line until you get to the ADVANCED FITTING step which
|
|---|
| 27 | will require you to experiment with the parameters.
|
|---|
| 28 |
|
|---|
| 29 | * ADDITIONAL exercise: Determine the rms noise in the spectrum using the
|
|---|
| 30 | create mask function you learned in the first section and the stat command
|
|---|
| 31 | specified in the text file.
|
|---|
| 32 |
|
|---|
| 33 | Examples - Fitting gaussians
|
|---|
| 34 | ----------------------------
|
|---|
| 35 |
|
|---|
| 36 | If you only see a file called *mon.asap.tar*, untar it before getting started
|
|---|
| 37 |
|
|---|
| 38 | .. code-block: shell
|
|---|
| 39 |
|
|---|
| 40 | tar xf mon.asap.tar
|
|---|
| 41 | asap
|
|---|
| 42 |
|
|---|
| 43 | .. code-block:: python
|
|---|
| 44 |
|
|---|
| 45 | # Load the scantable for the source into ASAP
|
|---|
| 46 | mon = scantable("mon.asap")
|
|---|
| 47 | mon.summary()
|
|---|
| 48 | # Select the on-source scans
|
|---|
| 49 | source=mon.get_scan([1,3])
|
|---|
| 50 | # Select the reference scans
|
|---|
| 51 | ref=mon.get_scan([0,2])
|
|---|
| 52 | # Make the quotient spectra and plot it
|
|---|
| 53 | quot=quotient(source,ref)
|
|---|
| 54 | plotter.plot(quot)
|
|---|
| 55 | # Remove the baseline and plot
|
|---|
| 56 | quot.auto_poly_baseline()
|
|---|
| 57 | plotter.plot(quot)
|
|---|
| 58 | # Average the 2 scans together
|
|---|
| 59 | av=quot.average_time(align=True)
|
|---|
| 60 | plotter.plot(av)
|
|---|
| 61 | # Check the the rest frequency is set correctly and define unit as km/s
|
|---|
| 62 | print av.get_restfreqs()
|
|---|
| 63 | av.set_unit("km/s")
|
|---|
| 64 | plotter.plot(av)
|
|---|
| 65 | # Scale the two polarizations separately
|
|---|
| 66 | sel=selector()
|
|---|
| 67 | sel.set_polarisations(0)
|
|---|
| 68 | av.set_selection(sel)
|
|---|
| 69 | av.scale(50.6)
|
|---|
| 70 | sel.set_polarisations(1)
|
|---|
| 71 | av.set_selection(sel)
|
|---|
| 72 | av.scale(32.5)
|
|---|
| 73 | # unset the selection
|
|---|
| 74 | av.set_selection()
|
|---|
| 75 | # Average the polarizations together
|
|---|
| 76 | avpol=av.average_pol()
|
|---|
| 77 | # Plot the data
|
|---|
| 78 | plotter.plot(avpol)
|
|---|
| 79 | # Create a mask that contains the emission
|
|---|
| 80 | msk=avpol.create_mask([-10,30])
|
|---|
| 81 | # Set up some fitting parameters
|
|---|
| 82 | f=fitter()
|
|---|
| 83 | f.set_scan(avpol,msk)
|
|---|
| 84 | f.set_function(gauss=2)
|
|---|
| 85 | # Fit a gaussian to the emission, plot the emission and fit
|
|---|
| 86 | f.fit()
|
|---|
| 87 | f.plot()
|
|---|
| 88 | # Save the plot
|
|---|
| 89 | plotter.save("monr2.eps")
|
|---|
| 90 | # Get the fit parameters
|
|---|
| 91 | f.get_parameters()
|
|---|
| 92 | # Plot the fit residuals
|
|---|
| 93 | f.plot(residual=True)
|
|---|
| 94 | # Plot each of the fitted gaussians separately, with or without the fit parameters overlaid
|
|---|
| 95 | f.plot(components=0,plotparms=True)
|
|---|
| 96 |
|
|---|
| 97 | f.plot(components=1)
|
|---|
| 98 |
|
|---|
| 99 |
|
|---|
| 100 | ADVANCED FITTING
|
|---|
| 101 | ----------------
|
|---|
| 102 |
|
|---|
| 103 | Sometimes it is necessary to force some of the fitting pa-
|
|---|
| 104 | rameters such as peak or FWHM. For example, if you want to fix the peak of the
|
|---|
| 105 | first gaussian component to 39 Jy you can use the following command:
|
|---|
| 106 |
|
|---|
| 107 | .. code-block:: python
|
|---|
| 108 |
|
|---|
| 109 | f.set_gauss_parameters(39, 10.2, 1, peakfixed=1, component=0)
|
|---|
| 110 |
|
|---|
| 111 | You can now experiment with forcing the values of the different parameters. To
|
|---|
| 112 | look at the help file on this function type
|
|---|
| 113 | help fitter.set_gauss parameters
|
|---|
| 114 | Additional Exercise: Determine the rms noise in the spectrum.
|
|---|
| 115 |
|
|---|
| 116 | 1) Create a mask that excludes the emission (if you leave it in the rms noise will be
|
|---|
| 117 | significantly inflated), using the same command as before except this time you will
|
|---|
| 118 | have to specify two ranges of values either side of the emission.
|
|---|
| 119 |
|
|---|
| 120 | 2) Find the rms noise of the spectrum:
|
|---|
| 121 |
|
|---|
| 122 | .. code-block:: python
|
|---|
| 123 |
|
|---|
| 124 | avpol.stats(stat="rms",mask=mymask)
|
|---|
| 125 |
|
|---|
| 126 | The stats function can be used to extract many other statistics from the data
|
|---|
| 127 | such as the maximum and minimum values, the median value and many more. Use
|
|---|
| 128 | help to find out how to extract these statistics from the spectrum.
|
|---|
| 129 |
|
|---|
