search

caracal-pipeline / caracal

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

Velocities after primary beam correction look strange #1470

Closed SpaceMeerkat closed 4 months ago

SpaceMeerkat commented 1 year ago

Perhaps related to closed issue: https://github.com/caracal-pipeline/caracal/issues/1230#issue-688516338

After running primary beam correction in the line worker..

    ...
  make_cube:
    enable: true
    npix: [1200]
    cell: 6
    robust: 0.5
    taper: 20
    niter: 1000000
    wscl_mgain: 1
    wscl_nrdeconvsubimg: 1
    wscl_sofia_niter: 3
  pb_cube:
    enable: true
    apply_pb: true
  remove_stokes_axis:
    enable: true
  freq_to_vel:
    enable: true 
  report: true

... the velocities in the resulting corrected cube are looking far too large. (Pre apply_pb: 32000 to -22000 km/s, after pb_corr: 299788.6497 to 299787.8834 km/s for my current MeerKAT cube).

Is this an artefact of the pb_cube vs freq_to_vel ordering still? And can it be easily undone by rerunning the line worker in a different way?

Many thanks!

paoloserra commented 1 year ago

Could you please copy here the header of the cube before and after freq_to_vel? Thanks!

SpaceMeerkat commented 1 year ago

I'm not sure how to get the cube prior to freq_to_vel without reurunning but here it is prior to pb_cube:

SIMPLE = T / Standard FITS
BITPIX = -32 / Floating point (32 bit)
NAXIS = 3
NAXIS1 = 1200
NAXIS2 = 1200
NAXIS3 = 1224
BSCALE = 1.000000000000E+00 / PHYSICAL = PIXEL*BSCALE + BZERO
BZERO = 0.000000000000E+00
BTYPE = Intensity
OBJECT = Abell4038
BUNIT = Jy/beam / Brightness (pixel) unit
EQUINOX = 2000
RADESYS = FK5
LONPOLE = 1.800000000000E+02
LATPOLE = -2.820277777778E+01
PC1_1 = 1.000000000000E+00
PC2_1 = 0.000000000000E+00
PC3_1 = 0.000000000000E+00
PC1_2 = 0.000000000000E+00
PC2_2 = 1.000000000000E+00
PC3_2 = 0.000000000000E+00
PC1_3 = 0.000000000000E+00
PC2_3 = 0.000000000000E+00
PC3_3 = 1.000000000000E+00
CTYPE1 = RA---SIN
CRVAL1 = -3.120416666667E+00
CDELT1 = -1.666666666667E-03
CRPIX1 = 601
CUNIT1 = deg
CTYPE2 = DEC--SIN
CRVAL2 = -2.820277777778E+01
CDELT2 = 1.666666666667E-03
CRPIX2 = 601
CUNIT2 = deg
CTYPE3 = FREQ
CRVAL3 = 1.270317544597E+09
CDELT3 = 2.090025580001E+05
CRPIX3 = 1
CUNIT3 = Hz
PV2_1 = 0.000000000000E+00
PV2_2 = 0.000000000000E+00
RESTFRQ = 1.420405751767E+09 / Rest Frequency (Hz)
TELESCOP = MeerKAT
OBSERVER = Sharmila
DATE-OBS = 2018-11-17T13:27:12.600000
TIMESYS = UTC
OBSGEO-X = 5.109217420824E+06
OBSGEO-Y = 2.006797675146E+06
OBSGEO-Z = -3.239133867335E+06
WSCDATAC = DATA
WSCVERSI = 3.0.1
WSCWEIGH = Briggs'(0.5)
WSCCHANE = 1.224000000000E+03
WSCCHANS = 0.000000000000E+00
WSCENVIS = 0.000000000000E+00
WSCFIELD = 0.000000000000E+00
WSCGAIN = 1.000000000000E-01
WSCGKRNL = 7.000000000000E+00
WSCIMGWG = 0.000000000000E+00
WSCMAJOR = 1.000000000000E+00
WSCMGAIN = 1.000000000000E+00
WSCMINOR = 0.000000000000E+00
WSCNEGCM = 1.000000000000E+00
WSCNEGST = 0.000000000000E+00
WSCNITER = 1.000000000000E+06
WSCNORMF = 0.000000000000E+00
WSCNVIS = 5.685000000000E+03
WSCNWLAY = 6.400000000000E+01
WSCTHRES = 0.000000000000E+00
WSCVWSUM = 0.000000000000E+00
SPECSYS3 = BARYCENT

and after:

SIMPLE = T / Standard FITS
BITPIX = -32 / Floating point (32 bit)
NAXIS = 3
NAXIS1 = 1200
NAXIS2 = 1200
NAXIS3 = 1224
BSCALE = 1.000000000000E+00 / PHYSICAL = PIXEL*BSCALE + BZERO
BZERO = 0.000000000000E+00
BTYPE = Intensity
OBJECT = Abell4038
BUNIT = Jy/beam / Brightness (pixel) unit
EQUINOX = 2000
RADESYS = FK5
LONPOLE = 1.800000000000E+02
LATPOLE = -2.820277777778E+01
PC1_1 = 1.000000000000E+00
PC2_1 = 0.000000000000E+00
PC3_1 = 0.000000000000E+00
PC1_2 = 0.000000000000E+00
PC2_2 = 1.000000000000E+00
PC3_2 = 0.000000000000E+00
PC1_3 = 0.000000000000E+00
PC2_3 = 0.000000000000E+00
PC3_3 = 1.000000000000E+00
CTYPE1 = RA---SIN
CRVAL1 = -3.120416666667E+00
CDELT1 = -1.666666666667E-03
CRPIX1 = 601
CUNIT1 = deg
CTYPE2 = DEC--SIN
CRVAL2 = -2.820277777778E+01
CDELT2 = 1.666666666667E-03
CRPIX2 = 601
CUNIT2 = deg
CTYPE3 = FREQ
CRVAL3 = 1.270317544597E+09
CDELT3 = 2.090025580001E+05
CRPIX3 = 1
CUNIT3 = Hz
PV2_1 = 0.000000000000E+00
PV2_2 = 0.000000000000E+00
TELESCOP = MeerKAT
OBSERVER = Sharmila
DATE-OBS = 2018-11-17T13:27:12.600000
TIMESYS = UTC
OBSGEO-X = 5.109217420824E+06
OBSGEO-Y = 2.006797675146E+06
OBSGEO-Z = -3.239133867335E+06
WSCDATAC = DATA
WSCVERSI = 3.0.1
WSCWEIGH = Briggs'(0.5)
WSCCHANE = 1.224000000000E+03
WSCCHANS = 0.000000000000E+00
WSCENVIS = 0.000000000000E+00
WSCFIELD = 0.000000000000E+00
WSCGAIN = 1.000000000000E-01
WSCGKRNL = 7.000000000000E+00
WSCIMGWG = 0.000000000000E+00
WSCMAJOR = 1.000000000000E+00
WSCMGAIN = 1.000000000000E+00
WSCMINOR = 0.000000000000E+00
WSCNEGCM = 1.000000000000E+00
WSCNEGST = 0.000000000000E+00
WSCNITER = 1.000000000000E+06
WSCNORMF = 0.000000000000E+00
WSCNVIS = 5.685000000000E+03
WSCNWLAY = 6.400000000000E+01
WSCTHRES = 0.000000000000E+00
WSCVWSUM = 0.000000000000E+00
SPECSYS3 = BARYCENT

CDELT3 doesn't look crazily different and they start in roughly the same spot.. which is strange

paoloserra commented 1 year ago

but cunit3 is still Hz after freq_to_vel?

SpaceMeerkat commented 1 year ago

Yes.. this also confused me. If you open the cube (pre pb-correction) in CARTA then it all looks fine but the header definitely still says the 3rd axis is in Frequency units.

paoloserra commented 1 year ago

Can you share the CARACal log?

SpaceMeerkat commented 1 year ago

log-caracal.txt

No problem!

SpaceMeerkat commented 1 year ago

I've noticed the convert-spectral_header lines have (None), I wonder if that's where the problem lies

paoloserra commented 1 year ago

I'll look at the log in a moment, but judging from the headers you've sent no conversion to velocity was made. I wonder what CARTA does in those cases. Does it attempt a conversion to velocity? If, so based on what rest frequency? In your second header there is no rest frequency, so maybe that is the problem with the strange velocity values in CARTA. While, on the CARACal side, the problem is not a wrong conversion, but the fact that there is no conversion at all.

SpaceMeerkat commented 1 year ago

I think it must, because for the pre-pb-corrected cubes, both frequency and velocity axes are listed in CARTA, but not for post-pb correction. So you're probably right in that CARTA requires an anchor value for the rest frequency and automatically shows velocities when it can. In which case... why no conversion on the CARACal side, that's the question as you say!

Or correction, what am I doing wrong to cause CARACal to not perform the conversion is probably more likely..

paoloserra commented 1 year ago

I've noticed the convert-spectral_header lines have (None), I wonder if that's where the problem lies

I think that (None) is normally where the version of the Stimela cab goes. See the other instances of Adding cab in the log. In this case, since freq_to_vel is a cab (actually, a Python function) defined in the line worker itself there is no version.

SpaceMeerkat commented 1 year ago

Ahhh I see, that makes sense!

SpaceMeerkat commented 1 year ago

I've noticed that the line worker is performing checks for 'restfreq' in the cube headers for freq_to_vel (see https://github.com/caracal-pipeline/caracal/blob/4f716f7c7f16ff5b7face76af053fb459a06bd77/caracal/workers/line_worker.py#L54). I've tried rerunning the line worker without applying freq_to_vel to spot any differences in the cube headers, and it seems that the header term RESTFRQ is added in somewhere (note the mismatch of restfreq and restfrq). Could this be a clue as to the problem, or would the log file show warnings from not finding the header keyword?

SpaceMeerkat commented 1 year ago

My mistake; after checking the headers in python instead of CARTA, it looks like the initial cube is in fact in velocities so freq_to_vel is working, but CARTA is showing the wrong header information. However, the primary beam corrected cube is definitely still output in frequency rather than velocity. So the problem is isolated to freq_to_vel not operating at the pb correction stage after all.

paoloserra commented 1 year ago

Sorry, I'm getting a little confused now. So the only problem you think we have is that freq_to_vel does not run on the PB-corrected cube? Or on the PB cube itself?

SpaceMeerkat commented 1 year ago

Sorry, yep I made it a little confusing there! Correct, freq_to_vel doesn't run on the PB corrected cube, the PB cube itself is fine.

paoloserra commented 4 months ago

@SpaceMeerkat sorry this took so long. It should be fixed now.