lucabaldini / ixpeobssim

Simulation and analysis framework for the Imaging X-ray Polarimetry Explorer
GNU General Public License v3.0
10 stars 13 forks source link

Trouble changing the binning in an event file before producing a PMAP cube with xpbin #729

Open app46 opened 3 months ago

app46 commented 3 months ago

Hi! I am using a separation method to extract parts of an original event file. I then intend to use the xpbin routine to create PMAP cubes and study the polarization of the extracted emission.

However, the separation method needs good statistics, so I had to rebin the original event file before applying it. For that reason, the resulting event file I want to apply does not have the same binning as the original. I tried to change the WCS information in the header of the fits file to reflect this change (in particular the TCDLT7 and TCDLT8 keywords, describing the X and Y pixel size in degree) but after numerous attempts, I still can't make it work.

Either I try to write the new WCS information in the header while saving in Python, but somehow the relevant keywords don't even appear in the header anymore, or I try to add them manually using fv. If I do this, I get this message while trying to run xpbin:

`>>> Reading WCS information from the EVENTS HDU... WARNING: FITSFixedWarning: CRPIX7 = 50.5 / Pixel coordinate of reference point image-header keyword CRPIXja in binary table. [astropy.wcs.wcs] WARNING: FITSFixedWarning: CRPIX8 = 50.5 / Pixel coordinate of reference point image-header keyword CRPIXja in binary table. [astropy.wcs.wcs] WARNING: FITSFixedWarning: CDELT7 = -0.0029166666666667 / [deg] Coordinate increment at reference point image-header keyword CDELTia in binary table. [astropy.wcs.wcs] WARNING: FITSFixedWarning: CDELT8 = 0.0029166666666667 / [deg] Coordinate increment at reference point image-header keyword CDELTia in binary table. [astropy.wcs.wcs] WARNING: FITSFixedWarning: CUNIT7 = 'deg ' / Units of coordinate increment and value image-header keyword CUNITia in binary table. [astropy.wcs.wcs] WARNING: FITSFixedWarning: CUNIT8 = 'deg ' / Units of coordinate increment and value image-header keyword CUNITia in binary table. [astropy.wcs.wcs] WARNING: FITSFixedWarning: CTYPE7 = 'RA---TAN' / Right ascension, gnomonic projection image-header keyword CTYPEia in binary table. [astropy.wcs.wcs] WARNING: FITSFixedWarning: CTYPE8 = 'DEC--TAN' / Declination, gnomonic projection image-header keyword CTYPEia in binary table. [astropy.wcs.wcs] WARNING: FITSFixedWarning: CRVAL7 = 83.63304 / [deg] Coordinate value at reference point image-header keyword CRVALia in binary table. [astropy.wcs.wcs] WARNING: FITSFixedWarning: CRVAL8 = 22.01449 / [deg] Coordinate value at reference point image-header keyword CRVALia in binary table. [astropy.wcs.wcs] WARNING: FITSFixedWarning: LONPOLE = 180.0 / [deg] Native longitude of celestial pole image-header keyword LONPOLEa in binary table. [astropy.wcs.wcs] WARNING: FITSFixedWarning: LATPOLE = 22.01449 / [deg] Native latitude of celestial pole image-header keyword LATPOLEa in binary table. [astropy.wcs.wcs] WARNING: FITSFixedWarning: MJDREF = 0.0 / [d] MJD of fiducial time image-header keyword MJDREF in binary table. [astropy.wcs.wcs] WARNING: FITSFixedWarning: RADESYS = 'ICRS' / Equatorial coordinate system image-header keyword RADESYSa in binary table. [astropy.wcs.wcs]

If no header is provided, keysel may not be provided either. Cannot parse the WCS information. Creating a standard WCS on the fly based on RA_OBJ and DEC_OBJ... This might get you going, but you should fix the input file. Building wcs object... 600 pixel(s) @ 0.00072 deg (0.433 deg image size) WCS Keywords

Number of WCS axes: 2 CTYPE : 'RA---TAN' 'DEC--TAN' CRVAL : nan nan CRPIX : 300.5 300.5 PC1_1 PC1_2 : 1.0 0.0 PC2_1 PC2_2 : 0.0 1.0 CDELT : -0.0007222222222222223 0.0007222222222222223 NAXIS : 600 600`

So for some reason the values I entered are not taken into account, the routine discards them and uses the original parameters, which of course doesn't work in the end. I don't really understand the warnings and I don't see what I could do differently, has anybody ever run into an issue like this? Would you know how I could fix it?

lucabaldini commented 3 months ago

Hi there, and my apologies for the late reply.

I am not sure I understand the question, particularly when you refer to "rebinning the original event file", since an event file is not binned by definition. (And I am sure I am missing something, but I am assuming by event file you mean a Level 2 file, right?)

Could you provide the exact sequence of commands you are issuing? That will help understand what the problem is.

Luca

app46 commented 3 months ago

Hi! The event file appears to be coded as a 600 per 600 grid, and I would like to reduce that to 100 per 100. I was actually able to change the header information by messing around with different astropy functions, but now xpbin does not seem to take the new information into account, and reads the new WCS as such:

CTYPE : 'RA---TAN' 'DEC--TAN' CRVAL : 83.63304 22.01449 CRPIX : 50.5 50.5 PC1_1 PC1_2 : 1.0 0.0 PC2_1 PC2_2 : 0.0 1.0 CDELT : -0.0029166666666667 0.0029166666666667 NAXIS : 600 600

where CRPIX and CDELT correspond to the new information I wrote in the header, but NAXIS remains 600 while it should be 100. As "600" does not appear in any header keyword anymore, I think the event file size might be hard coded in the xpbin routine?

And by the way, thank you for your taking the time to answer!

lucabaldini commented 3 months ago

Ok we don't seem to be understanding each other :-)

The WCS information in the event files serves the sole purpose to convert from pixel to sky coordinates the event positions (columns X and Y in the EVENTS binary table), and you should not change that by hand. I maintain that, assuming that by "event file" you mean a "Level 2 file", that is not "binned".

If you want to produce a 100 x 100 binned map with xpbin, you should use the --npix command-line switch (and, optionally, --pixsize). So, assuming that you have a level-2 file and you want to create a 100 x 100 pixel polarization map cube, the command would be something along the lines of

xpbin(.py) --npix 100 path_to_the_input_l2_file.fits
app46 commented 3 months ago

Isn't the data always necessarily binned, in that this is computing and continuum doesn't exist? In the header from the event file there is TLMIN and TLMAX keywords for X, Y, and E, defining the maximum resolution of the file: 600 x 600 x 375. It is way finer than the actual physical resolution of IXPE, so it may be considered unbinned, but this is precisely this resolution I am trying to change. But as what I am trying to do does not seem to be very common, I suspect this "maximum resolution" info to be hard coded, so impossible to modify prior to using xpbin...