| 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 | |
|---|
