editeodoro / Bbarolo

Bbarolo is a 3D fitting tool to derive the kinematics of galaxies from emission-line observations.
http://editeodoro.github.io/Bbarolo/
GNU General Public License v2.0
29 stars 8 forks source link

Question about the moment map task #32

Closed aparnavenkat27 closed 5 months ago

aparnavenkat27 commented 9 months ago

Hello Enrico!

Thank you for the updated version of Barolo - I am keen on trying it out soon! I have questions about the S/N maps constructed in the moment map task. You once pointed out that it is a non-trivial operation to determine the contours on a moment 0 map based on the S/N cuts we apply to a datacube. I read your references, but unfortunately, I kept going in circles due to my specific case.

I obtained my datacube using the CASA task tclean using a robust weighting parameter of +0.5, i.e., favouring resolution over sensitivity. Is using "UNIFORM" in the TAPER parameter correct in my case? I did not smooth the datacube in the analysis, so the other two keywords seem inappropriate. Furthermore, would it be possible to know the sigma/noise value at the end of the task so I can use it to apply contours to my maps directly?

Cheers, and thanks for all the work you do!

editeodoro commented 9 months ago

Hello Aparna! If you made your cube in CASA without any tapering of the data (as I think), just use TAPER=UNIFORM. The SNMAP task will produce a RMS map, with the values of the RMS as a function of spatial pixel (which will depend on the masking options), and a S/N map, i.e. the variation of S/N as a function of spatial pixel (this is just moment0 / RMSmap) . This should be all you need, and you can use these maps to draw the S/N contours on the 0th moment.

Hope this is clear! Enrico

aparnavenkat27 commented 9 months ago

Thanks for the quick response, as always!

I have a slightly bigger problem now - BBarolo v1.7 isn't working on my machine. Until this past Friday, I was able to run v1.6 with no issues on MacOS 14.0. I updated to v1.7 using homebrew. How do I find out what is going on?

To answer your earlier assumption - I do not apply any tapering via CASA, but I do apply the Briggs weighting scheme. I hope that once this new issue is fixed, I will be able to make the moment maps ASAP!

editeodoro commented 9 months ago

Well, I need some more information. What problem are you having?

aparnavenkat27 commented 9 months ago

Screen Shot 2023-10-23 at 19 19 38 PM

I attached a screenshot - this is all it shows. The only difference is that I updated to v1.7 a few hours ago.

Edit - I modified the ruby script slightly to reinstall v1.6, which also does not work for codes that worked three days ago.

editeodoro commented 9 months ago

Hi Aparna, it looks like there is a problem with BBarolo when compiled with the latest Apple default compiler (Clang). I have now updated the Homebrew recipe to use GNU GCC instead. Please try again to download the file bbarolo.rb from the repo and re-run it through Homebrew. The resulting binary file should now work properly!

Let me know if this does not fix the issue. Enrico

aparnavenkat27 commented 9 months ago

Thank you for the quick fix! I re-ran all my old codes and the new moment map tasks with no issues. My only question is now about the noise level estimated by the code. Is there a way to supply the sigma value for the datacube, the way we do for the search task, and then determine the S/N map for the Moment 0 map? Or is this what happens when the moment map task uses the search task to establish a mask?

editeodoro commented 9 months ago

I am not sure to understand your question....the whole point here is that, when you make a moment0 with a mask, you're summing up a different number of channels in each spatial pixel, therefore the noise level in the moment0 map varies as a function of spatial position (= the RMSmap produced by the code). This is independent from how you make the mask; whatever the mask chosen (source finder, sigma clipping, smooth, etc...), the code will calculate the RMS map based on the mask.

For the masking, you can set all the parameters supported by the corresponding task that makes the mask. For example, if you want to use the SEARCH algorithm for masking, you can provide your own value of noise level with the THRESHOLD parameter.

Not sure I answered though...

aparnavenkat27 commented 8 months ago

Right, now it makes sense! I was getting twisted over the RMS value to report and got confused in the process. Thank you again!

aparnavenkat27 commented 6 months ago

Hello! I am reopening this because I have a few quick questions about the results produced by Barolo in the moment map tasks - I see two RMS maps in the directory after running the program with RMSMAP=True. Why are there two, and what do they signify? At first glance, the file with "_0th_RMS" appears to have a Jy/beam.km/s unit, and the file with "_RMS" map has Jy/beam. I presume it is the former that is used to generate S/N maps as you stated above due to the units? In this case, why is it that when my input cube contains values less than 0, the "_0th_RMS" file appears as a binary-valued (0's and 1's) fits file?

Next, if I wanted to plot the datacube RMS on the residual maps generated in 3DFIT to find if they have emissions greater than 1sigma, how do I go about it? Is there a relationship I can use to say if 0.5mJy is the noise in a datacube with x channels and y channel width, then z is the corresponding noise in the moment 0 map? Is there a possibility to have a goodness-of-fit plot to quantitatively explain how well the model fits the data?

Cheers, and a happy new year to you! :)

editeodoro commented 6 months ago

Hello Aparna, the "_0th_RMS" file appears when you set TOTALMAP=true together with SNMAP=true. In this case, the code produces the RMS map of the total intensity map, which is in Jy/beam*km/s and depends on the spatial position because you're using a mask, i.e. because in each spatial pixel you're summing over different number of channels. If you do not use a mask (MASK=NONE), you'll see that the "_0th_RMS.fits" file is a constant across the entire field.

Instead the "_RMS" file is produced when RMSMAP=true. In this case, the code, for each spaxel computes the RMS in emission-free channels, so the final map is in the same units as the datacube (Jy/beam). This does not concern the moment maps, but it can be useful to see how the noise varies with position in a datacube (for example, if you've done a OTF map with a single dish).

I am not sure to understand this: "why is it that when my input cube contains values less than 0, the "_0th_RMS" file appears as a binary-valued (0's and 1's) fits file?". This sounds wrong, can you send me a working example (FITS file and parameter file) to reproduce this?

To evaluate the goodness of the fit, I would not rely on moment maps. Why don't you just calculate a residual cube (data-model) and see in each channel if you have regions where you have residual emission (or negative residuals!) at a n*sigma level, where sigma is the RMS of the datacube?

Enrico

aparnavenkat27 commented 6 months ago

Thanks for explaining the two files! I did not realise that the "0th_RMS" file had a constant value across, I just saw the histogram flat and wondered what was going on. That is my bad!

I have at the moment performed a "data-model" in the code I am plotting with, and I make a moment 0,1 and 2 map out of this residual cube. If I were to use the "0th_RMS" file to divide the residual moment 0 map that I build, and then plot a 3/5 sigma contour in matplotlib, does this give me the correct result I need? I need to plot the residuals so that visually there is an understanding of how well the model fits the data.

Lastly, I use the SEARCH algorithm to find a source in a given sub-cube, provide (growth-)thresholds, and perform a fitting. When I am preparing the above residual maps, I assume that the SEARCH algorithm should be used again with the moment maps task to maintain consistency between the contours and data?

editeodoro commented 6 months ago

I think it does not make much sense to calculate a residual cube and then make moment maps with it, this would not give you much information on the goodness of the fit. Especially the residual 0th moment does not tell anything on the kinematics, so what's the point??!!

But in general, the 0th_RMS file is calculated for a given mask, e.g. for a mask that at each spatial position makes a sum over a number of channels. So, you need to use the same exact mask, not only the SEARCH algorithm with the same parameters (because it could give a different mask if applied to either the original cube or to the residuals).