caracal-pipeline / caracal

Containerized Automated Radio Astronomy Calibration (CARACal) pipeline
GNU General Public License v2.0
29 stars 6 forks source link

Unexpectedly small restoring beam for continuum image #1174

Open healytwin1 opened 4 years ago

healytwin1 commented 4 years ago

This is more a question than an issue (I think). The restoring beam in my continuum image (1300-1430MHz) is 3.8"x7.6" which seems very small. The output from wsclean has the theoretical beam at 5.62", so I am a little confused as to how it gets to ~3.8" in the final image. Any ideas?

I am running caracal 1.0.3, using wsclean 2.9.

KshitijT commented 4 years ago

It's very elliptical too ! How much flagging does that band have, btw?

healytwin1 commented 4 years ago

The ellipticity is expected - the field is pretty far north. Only about 10% of the field is flagged

paoloserra commented 4 years ago

Could you have a look at the dirty PSF? Is the FWHM of the restoring beam consistent with the dirty PSF? Or is this maybe a case in which the dirty PSF was poorly fit when determining the most appropriate restoring beam?

healytwin1 commented 4 years ago

I've had a look at the PSF, it is consistent with the restoring beam.

healytwin1 commented 4 years ago

So I thought this might be a wsclean 2.9 thing, I tried with wsclean 2.6 and 2.8, with the same result - a theoretical beam of 5.79'' and a fitted beam: major=7.67'', minor=3.79''. Having looked at the psf images, they look okay and match the fitted beam.

Not sure whether to accept this, or whether it is a bug

paoloserra commented 4 years ago

Maybe you could try to see whether based on the maximum uvdist a beam minor axis of 3.8'' is reasonable.

bennahugo commented 4 years ago

That does not sound too unreasonable to me if the field does not have a high peak elevation and/or poor hour angle coverage.

Here is a fitted beam at -0.6 robust with DDFacet for a short observation of one of my science fields: In [3]: 0.002291683546348191*3600.00 Out[3]: 8.250060766853489

In [4]: 0.001246866961825903*3600.0 Out[4]: 4.4887210625732505

On Thu, Jun 18, 2020 at 3:08 PM paoloserra notifications@github.com wrote:

Maybe you could try to see whether based on the maximum uvdist a beam minor axis of 3.8'' is reasonable.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/caracal-pipeline/caracal/issues/1174#issuecomment-646005260, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB4RE6UCNEROCHNSMOKFLEDRXIGTFANCNFSM4NZVUHDQ .

--

Benjamin Hugo

PhD. student, Centre for Radio Astronomy Techniques and Technologies Department of Physics and Electronics Rhodes University

Junior software developer Radio Astronomy Research Group South African Radio Astronomy Observatory Black River Business Park Observatory Cape Town

healytwin1 commented 4 years ago

Maybe you could try to see whether based on the maximum uvdist a beam minor axis of 3.8'' is reasonable.

The max uvdist is 6812m, so I would expect the minor axis to be a factor of 2 higher. I have been using robust of -1 and the central frequency is 1350MHz

paoloserra commented 4 years ago

You're using multi-frequency-synthesis imaging so what matters is the highest frequency in your data.

healytwin1 commented 4 years ago

Even with 1430MHz which is the highest frequency, I would expect a larger minor axis

paoloserra commented 4 years ago

Indeed, if I'm doing this right lambda/baseline is 6.4 arcsec at 1430 MHz for a baseline of 6.8 km. Not quite a factor of 2 compared to your 3.8 arcsec, but nearly. That's too much of a difference to. be caused by an incorrect fit.

With MeerKAT you should get in the best case a 7.7 km baseline at a frequency of 1750 MHz, which gives a lambda/baseline of 4.6 arcsec.

This seems to be consistent with previous images we've made in continuum (see, e.g. Filippo's paper on Fornax A).

So I'm a bit puzzled.

@bennahugo can you share more details on your case? What's the highest frequency and the longest baseline of the data for which you're getting the 4.5 arcsec beam?

bennahugo commented 4 years ago

I believe it is the sigma of the gaussian that is quoted, not the FWHM. In that case a value around 3 is quite sound.

On Thu, Jun 18, 2020 at 6:45 PM Julia H notifications@github.com wrote:

Even with 1430MHz which is the highest frequency, I would expect a larger minor axis

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/caracal-pipeline/caracal/issues/1174#issuecomment-646160327, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB4RE6UTUJIGBHR7KL4CMMDRXJADJANCNFSM4NZVUHDQ .

--

Benjamin Hugo

PhD. student, Centre for Radio Astronomy Techniques and Technologies Department of Physics and Electronics Rhodes University

Junior software developer Radio Astronomy Research Group South African Radio Astronomy Observatory Black River Business Park Observatory Cape Town

bennahugo commented 4 years ago

This dataset has the longest spacings in it, but each image is only 2 hrs of integration so the beam is very elliptical as expected. I synthesize with the full bandwidth.

On Thu, Jun 18, 2020 at 7:02 PM Benna Hugo bennahugo@gmail.com wrote:

I believe it is the sigma of the gaussian that is quoted, not the FWHM. In that case a value around 3 is quite sound.

On Thu, Jun 18, 2020 at 6:45 PM Julia H notifications@github.com wrote:

Even with 1430MHz which is the highest frequency, I would expect a larger minor axis

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/caracal-pipeline/caracal/issues/1174#issuecomment-646160327, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB4RE6UTUJIGBHR7KL4CMMDRXJADJANCNFSM4NZVUHDQ .

--

Benjamin Hugo

PhD. student, Centre for Radio Astronomy Techniques and Technologies Department of Physics and Electronics Rhodes University

Junior software developer Radio Astronomy Research Group South African Radio Astronomy Observatory Black River Business Park Observatory Cape Town

--

Benjamin Hugo

PhD. student, Centre for Radio Astronomy Techniques and Technologies Department of Physics and Electronics Rhodes University

Junior software developer Radio Astronomy Research Group South African Radio Astronomy Observatory Black River Business Park Observatory Cape Town

paoloserra commented 4 years ago

I believe it is the sigma of the gaussian that is quoted, not the FWHM. In that case a value around 3 is quite sound.

Mmm, really, PSF size in sigma? That'd be a first, but yes, let's check whether that is the case. Maybe @o-smirnov knows?

bennahugo commented 4 years ago

Well it is a Gaussian we are talking about... That is what is fitted for in DDF. I haven't perused that part of the source code of WSClean but maybe oleg knows.

Any case the value does not surprise me.

Cheers,

On Thu, Jun 18, 2020 at 8:16 PM paoloserra notifications@github.com wrote:

I believe it is the sigma of the gaussian that is quoted, not the FWHM. In that case a value around 3 is quite sound.

Mmm, really, PSF size in sigma? That'd be a first, but yes, let's check whether that is the case. Maybe @o-smirnov https://github.com/o-smirnov knows?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/caracal-pipeline/caracal/issues/1174#issuecomment-646227464, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB4RE6Q67POW4T5OLUX5LTDRXJKXFANCNFSM4NZVUHDQ .

--

Benjamin Hugo

PhD. student, Centre for Radio Astronomy Techniques and Technologies Department of Physics and Electronics Rhodes University

Junior software developer Radio Astronomy Research Group South African Radio Astronomy Observatory Black River Business Park Observatory Cape Town

healytwin1 commented 4 years ago

I compared the psf image created with an earlier set that I created in December to experiment with different robust weightings. See the attached plot, the green/purple are cross-sections through the major/minor axis of the December images, and the orange is from my most recent image -- there is a very clear difference. The green is robust of -1.5 and the purple is robust of 1.5 image

I believe it is the sigma of the gaussian that is quoted, not the FWHM. In that case a value around 3 is quite sound.

I do not believe the numbers quoted are the sigma are the gaussian, but rather than the FWHM as the numbers match what I get when I fit a gaussian to the psf image and work out what the FWHM is.

o-smirnov commented 4 years ago

Yeah I seem to recall it's the FWHM actually, but what do I know.

Have the short spacings gotten flagged into submission somehow?

healytwin1 commented 4 years ago

Have the short spacings gotten flagged into submission somehow?

We thought this may be a possibility, but only about 10% of the field is flagged.

paoloserra commented 4 years ago

So I got a small .MS from 1407 to 1414 MHz. When imaging with robust = -1.5 WSClean says the following:

# +-+-J- 0  I   0  0  0  0  1407-1414 (7)
# 
# Reordering /stimela_mount/msdir/1588316462_sdp_l0-NGC1365_cont.ms into 1 x 1 parts.
# Reordering: 0%....10%....20%....30%....40%....50%....60%....70%....80%....90%....100%
# Initializing model visibilities: 0%....10%....20%....30%....40%....50%....60%....70%....80%....90%....100%
# Detected 440.8 GB of system memory, usage not limited.
#  == Constructing PSF ==
# Opening reordered part 0 spw 0 for /stimela_mount/msdir/1588316462_sdp_l0-NGC1365_cont.ms
# Precalculating weights for Briggs'(-1.5) weighting... DONE
# Determining min and max w & theoretical beam size... DONE (w=[5.01233e-05:30491.7] lambdas, maxuvw=36247.4 lambda)
# Theoretic beam = 5.69''

This makes sense to me. For maxuvw = 36247.4 lambda at 1414 MHz (max_baseline = 7684 m) I obtain lambda/max_baseline = 5.69''.

So WSClean reports FWHM, not sigma, as customary.

Then WSClean continues:

# Rendering sources to restored image (beam=4.06''-4.98'', PA=137.61 deg)... DONE

That's a fair bit smaller than the theoretical beam.

Note that I had imaged with a ridiculously small pixel size of 0.2 arcsec in order to inspect the PSF visually without being limited by its sampling. This is what I see:

Screen Shot 2020-06-19 at 11 06 35

The background image shows the dirty PSF. The black contour is drawn at half-maximum. The weird colour-coding is meant to show that the black contour is a reliable description of the half-maximum level of the dirty PSF. The major and minor axis measured (not fitted) from the black contour are 5.7'' and 5.4'', respectively.

The green contour shows the restoring beam taken from the .FITS header (which is 4.06'' x 4.98'' as reported by WSClean). The restoring beam is indeed sgnificantly smaller than the dirty PSF, as in @healytwin1's case.

I've repeated this experiment with the maximum pixel size allowed by the theoretical PSF, i.e., 5.7''/3 = 1.9''. In this case WSClean states the same theoretical PSF FWHM while it restores as follows:

# Rendering sources to restored image (beam=5.48''-5.87'', PA=133.64 deg)... DONE

So using a larger pixel size delivered a restoring beam FWHM closer to the theoretical PSF FWHM.

For some details on this see https://sourceforge.net/p/wsclean/wiki/DeterminingPSFSize/ .

I realise that I'm not offering a solution to @healytwin1's problem yet, but hopefully this clarifies what's going on.

PeterKamphuis commented 4 years ago

Could this be driven by the beam fitting area (Maybe due to the first negative sidelobe)?

The wsclean wiki does say:

In WSClean 2.6, the option -beam-fitting-size was added. It sets the box size to be used during the fit. If the beam that is fitted is not what you can expect, you can try changing this value. It seems that in particular lowering the value to 1-3 can sometimes solve certain issues.

I found it helped me with the highly elliptical beams caused by non-circular sidelobes but I'm not sure of this issue.

paoloserra commented 4 years ago

Yes, the fitting is surely going funny and it might be worth trying a smaller fitting box size than default. Maybe something for @healytwin1 to try.

healytwin1 commented 4 years ago

I'm not sure how to run wsclean outside of caracal, so am not sure how to change the beam-fitting-size parameter.

paoloserra commented 4 years ago

That option is not in Stimela, so it's not available in CARACal either.

But @healytwin1, I'm puzzled by your beam profiles. In your case it does seem like the latest PSF (orange profile) is very different from the old ones. Is that the same data?

healytwin1 commented 4 years ago

But @healytwin1, I'm puzzled by your beam profiles. In your case it does seem like the latest PSF (orange profile) is very different from the old ones. Is that the same data?

Yes, although the purple/green profiles come from when I created the HI cube in Jan/Feb and the orange one from this week which while created from the same data, has had the calibration redone using the newer more stable caracal.

paoloserra commented 4 years ago

The beam depends on the UV coverage. The calibration is irrelevant. At what frequency were the old (HI) beams made? For an HI cube you get a beam per channel. Also, your old beams never go negative before the first sidelobe (left panel). That's not what I expect for MeerKAT.

gigjozsa commented 4 years ago

@healytwin1 , maybe you could quickly FFT your beam and look at the resulting amplitudes. That would provide you with the applied sampling function (otherwise looking at the uv coverage suffices). If that is looking fishy, we have a problem.

gigjozsa commented 4 years ago

If this is too complicated using Python (I would need to find out myself), Miriad is your friend to do it very quickly using the task fits and fft

KshitijT commented 4 years ago

I'm not sure how to run wsclean outside of caracal, so am not sure how to change the beam-fitting-size parameter.

On ILIFU, you could use singularity shell /idia/software/containers/STIMELA_IMAGES_1.6.1/stimela_wsclean_1.6.0.sif or singularity exec for scripts?

I'll add the parameter to stimela - there's already an issue open on it.

"The calibration is irrelevant."

Great tagline for the CARACal logo?

healytwin1 commented 4 years ago

The HI beams are around 1312MHz, so not quite the same as the new one, but I still don't expect such a huge discrepancy.

Also, your old beams never go negative before the first sidelobe (left panel). That's not what I expect for MeerKAT.

Yeah, that was interesting to me too when I first saw this back in Feb, but there were no comments from the others that looked at the images at the time, so I didn't think too much of it. The pixel scale is 4'' for the HI beams, so it is possible that the negative sidelobes get a little washed out because of that.

@gigjozsa I haven't used miriad before, I can try have a look in python over the weekend but I haven't ever tried an FFT with python before

gigjozsa commented 4 years ago

@healytwin1 if you can safe one image of the beam (1 plane) here or somewhere private I can do that quickly for you and send back the result and a manual (which will be very brief).

paoloserra commented 4 years ago

The pixel scale is 4'' for the HI beams, so it is possible that the negative sidelobes get a little washed out because of that.

4'' pixels with no uv tapering?

paoloserra commented 4 years ago

also, just noticed, why don't you plots peak at 1?

PeterKamphuis commented 4 years ago

In the old profiles it does look like the minor axis FWHM is quite undersampled.

PeterKamphuis commented 4 years ago

And, I'm a bit out of it but, do we have command logging already back in Caracal? In that case you can just run wsclean by copy pasting the command, modify the file names and add the -beam-fitting-size parameter.

healytwin1 commented 4 years ago

@healytwin1 if you can safe one image of the beam (1 plane) here or somewhere private I can do that quickly for you and send back the result and a manual (which will be very brief).

I tried doing it in python, but I'm not sure I did it probably. I can't upload the psf here, but have uploaded it to a google drive (https://drive.google.com/file/d/1H7EIlY97Khzkib1jwIxJF-ixRb85OKN0/view?usp=sharing)

The pixel scale is 4'' for the HI beams, so it is possible that the negative sidelobes get a little washed out because of that.

4'' pixels with no uv tapering?

Yes, that is correct.

also, just noticed, why don't you plots peak at 1?

That is my bad, I wasn't going through the centre of the image - see below: image

paoloserra commented 4 years ago

The pixel scale is 4'' for the HI beams, so it is possible that the negative sidelobes get a little washed out because of that.

4'' pixels with no uv tapering?

Yes, that is correct.

I think that needs to be corrected. You need to sample the longest baseline resolution (FWHM) with 3 pixels, so based on the numbers above you need 2'' pixels.

gigjozsa commented 4 years ago

If you want to be pedantic and you know the baseline lengths you can also divide 1 by the longest baseline in units of lambda, then divide it again by (2 x sqrt(2)) (that's approximately 3) to get the largest allowable pixel size in radians.

healytwin1 commented 4 years ago

Pixel sizes aside, any ideas on what could be causing the smaller than expected restoring beam? Most of my confusion stems from the fact that older versions of caracal have yielded the expected beamsize, but the more recent versions seem to be under estimating the size. I admit that I don't know enough about how all the software pieces work or fit together to identify what the source of the problem is. Regardless of pixel size and robust value, I expect the restoring beam around 1430 MHz to be (at its smallest) around 6arcsec, but that is clearly not the case here. I have tried pixel sizes ranging from 1.5arcsecs to 2.5arcsec (which is the pixel scale for our final HI cubes), see figure below (all psfs are from roughly the same frequency - 1430MHz). The blue profile comes from the final HI cube psf created in Jan/Feb while the green & orange come from my most recent continuum images. The 4 arcsec pixels were for a test to see how the beam varied with robust value. image