davidt3/xga

X-ray: Generate and Analyse is a module designed to make the analysis of XMM observations simple and efficient. It provides an interface with SAS for the creation of XMM data products, as well as a way to easily perform fits (scalable for multiple observations) and retrieve information about an object, all within a Python package.

xga's Issues

sas.py command stacking system could be improved

Although the stacking functionality remains, I don't think it really needs to. If it were removed it would clean up sas.py slightly, as well as some of the BaseSource class.

Initial coordinates quite clearly inside an observation aren't being matched to that observation

Was exploring off axis clusters using the XCS-DES sample Reese provided, found a very bright off axis cluster that is ideal for my needs, his table reports it having XMM coordinates of:

RA: 313.98311792357896
DEC: -54.9302122469135

Which quite clearly lines up with a cluster and is solidly within this observation (0675010901).

However, when I define an XGA source it cannot find a matching source for some reason.

Some way to view PSFGrid components would be good

It would be nice to have some view() method that showed the PSF components in-situ somehow

HPC users won't be able to switch easily between submitting jobs and running locally

Could be very inconvenient for people that use interactive sessions to test things (for instance), then submit jobs for big samples.

Documentation is none existent

What it says in the title. I haven't written any documentation for this, and as it stands I'm the only one who knows how to use it...

As I want this to be used by others in XCS (undergraduates primarily), everything will need to be properly explained.

I have created a readthedocs project for this, but have been unable to get Sphinx working to generate some documentation from my docstrings.

GalaxyCluster that isn't detected in region files and that isn't declared with a custom region will not be allowed to continue

This seems to be an oversight on my part. When the ExtendedSource init is called with super from GalaxyCluster init, it checks to see if there is a custom region, if there is no custom region and no matching region files then it can't continue because it doesn't know about the overdensity radii required by a GalaxyCluster object.

Simple Match is very limited

Feel like I should attempt to calculate the amount of a given source that lies on a detector, raise flags if its too low, if the user coordinates lie in a chip gap etc.

evselect_spectrum should be able to generate core-excised spectra

Many global X-ray measurements of clusters are performed on 'core-excised' spectra, so that should be an option in XGA as well. A core-excised flag should be included in the file name - fitting and results storage mechanisms will need to be adjusted to be able to store fits to core-excised and non-core excised spectra of the same source

Suggested rewrite of file reading

XGA/xga/products.py

Line 259 in 30b6d1e

with fits.open(self.path, mode="readonly", memmap=False) as temp_im:

with fits.open(self.path, mode="readonly", memmap=False) as temp_im:
    self._data = temp_im[0].data.astype('float32')
    if self._data.min() < 0:
        warnings.warn(...)
        np.maximum(self._data, 0, out=self._data)

temp_im.close(verbose=True)
del temp_im
gc.collect()    

I'm sure this is probably what you had to start off with, but it achieves what you want with the creation of as few extra variables as possible (so we don't end up with unwanted references).

It's so frustrating that you can't reliably recreate this error...

Naming a source in the standard astronomy style results in spectrum generation issues

For instance, defining a source with name="XMMXCSJ142601.0+374937.0", results in a SAS error:

xga.exceptions.SASGenerationError: datasetCouldNotBeRead raised by rmfgen - The dataset with name '0112230201_mos1_XMMXCSJ142601.0' could not be read; either it does not exist or it is in an unrecognised format

The dataset name reported in the error is not the full spectrum path, its cut off just before a "+" symbol, suggesting some SAS procedure has an issue with that particular symbol.

The error does not seem to occur when the "+" is replaced by a "-"

If a source is not detected in region files, source declaration will just die

Averaging Lx values from different spectra is crude

Certainly not statistically valid to just average them and their uncertainties.

Not sure explicit copying is needed

I'm pretty sure you're copying the data above when you byteswap, so this line isn't needed :)

XGA/xga/products.py

Line 271 in 30b6d1e

new_data = np.empty_like(data)

The GalaxyCluster object needs to require overdensity radii passed in

For the GalaxyCluster source object to be useful, it needs to require the passing of overdensity radii. One of the simplest ways we characterise the extent of a cluster.

It then needs to setup regions with those radii, and centre them either on the user-passed coordinate, or the X-ray centroid.

Deal with missing overdensity radii

Should I try and deal with missing overdensity radii by using scaling relations and iterating?

view_brightness_profile is plotting out to the complete wrong radius in kpc, but only for some clusters

For instance, having set an r500 of 1000kpc for a really quite off axis cluster, the profile goes out to ~1700kpc.

However, in the same script I set the same radius for an on axis cluster and it plots correctly.

Suspect its to do with how I do radius conversions in ann_radii.

evselect_spectrum and single_temp_apec should save positions and radii in filesnames, not just region type

Currently, a custom region will just save with a 'custom' flag in the spectrum filename, but it perhaps it should depend on where the custom region is so that multiple custom regions can be saved and analysed?

Account for PSF effects in images

Be able to de-convolve a PSF model with ratemaps.

Should be useful for multiple things:

• Should probably be applied before peak finding - though you'd hope it wouldn't have much of an effect
• Must be accounted for so that we can construct accurate brightness profiles - and by extension density profiles
• Hopefully useful for point source analysis.

PSF files aren't re-read in on source declaration

What it says in the title. When I was writing it I wasn't really bothered because they're really quick to generate.

Good to throw a custom warning if appropriate, people can catch them if it's annoying

i.e. create an "UnphysicalDataWarning", and warn whenever there is data that doesn't make sense

people can set warnings.simplefilter('ignore', UnphysicalDataWarning) if they don't want it

XGA/xga/products.py

Line 267 in 30b6d1e

# warnings.warn("You are loading an {} with elements that are < 0, "

Need a good demonstration notebook for the alpha release - in lieu of documentation

I should make a new folder with demonstration notebooks in, showing off everything I've managed to implement so far.

get_masks successfully returns an r500 mask when also supplied with obs_id and instrument

I don't currently generate overdensity radii masks/ custom masks for individual observations and instruments, so get_masks shouldn't successfully return anything in this case.

TypeError: unsupported operand type(s) for +: 'float' and 'NoneType' when constructing masks with all interlopers

Don't know why this is happening. Presented on a cluster that is quite off axis on one observation.

Initial thoughts are there is a None is making its way into the interloper source list.

Generate radial brightness profiles

From the ratemaps and masks construct radial brightness profiles.

Ideally I'd like to be able to do this with a dependence on angle as well, and scrap some of the assumption of spherical symmetry, but initially I think just radial profiles will be fine.

Should generate overdensity/custom radius masks for individual observations, not just combined

I think it would be a useful thing to do, and it wouldn't be hard. People might want to look at individual observations.

Store the point clusters used to select the X-ray peak

The hierarchical clustering peak finding method seems to work very well, but I should return the point clusters from the method and store them as attributes of the source class.

That way they can be used to make plots to check the peak finding has worked well.

Maybe double check the byte ordering before byteswapping?

Might be good to double-check if byteswapping is needed before executing.

This could also be done in-place to save on RAM

XGA/xga/products.py

Line 262 in 30b6d1e

data = data.newbyteorder().byteswap()

Paul gave me feedback on the general run through notebook

Bullet points of his main thoughts:

• Add a table of contents at the top, and state what the point and end goal of the notebook is more clearly.
• I noticed that the order didn't make too much sense in places, shuffle it around a bit.
• Make it really clear that the events lists going into XGA should be cleaned.
• Explain the configuration file a bit, link to an example.
• Highlight more that the sources are the most important part of the whole module.
• Explain regions better, with better distinction between custom, overdensity, and region file regions.
• Explain that the different source objects match to region files differently.
• Move the mask retrieval to the cell where its used.
• The product section is messy, can be split up a bit.
• Remove extra information about eexpmap generating images if necessary, doesn't matter to the user.
• Explain better what products are.

Add SFTP guide to README

Add advice and instructions on setting up an SFTP connection and using XGA like that. Possibly most useful for XCS members

Setup.py requirements seem to cause problems sometimes

I tried to install XGA v0.1.0 using pip on the workstation and I think the version requirements were causing some issues. Either way I should simplify them to greater than equal statements were possible. Regions has to remain at v0.4.

This will also need to change in the requirements.txt.

Also, as something also in setup.py that isn't worth of its own issue, I need to add a long description to setup.py so it appears on PyPi on the next release.

XSPEC interface failure running in Jupyter Notebook

When an XSPEC fit is run from a Jupyter Notebook, something funny goes on with the creation of the summary fits file. The actual XSPEC run and fit complete fine, but the script fails when it comes to write out the results.

Errors like this get thrown:
Unable to redirect prompts to the /dev/tty (at headas_stdio.c:152)

ERROR: No such device or address

Task ftcreate 0.0 terminating with status 6

The same fits work fine when run outside of Jupyter Notebook

Extended Source Peak finding fails if no redshift is passed

I made the design choice to use a search aperture of radius 500kpc for the iterative peak finding algorithm, but of course that proper radius can't be converted to degrees and then to pixels without redshift information.

Not difficult to fix, I'll probably just add an if statement that says if there is no redshift use a search aperture of 5 arcminutes or something.

evselect_spectrum doesn't support grouping spectra

Currently it is not possible to generate grouped spectra, as many people do. This should be added as an option in evselect_spectrum.

In the case where someone requests a grouped spectrum, I can just add a specgroup command on the end of the stack, and overwrite the spectrum file.

Fit Result Units

Would really like to read out XSPEC results with the right units if possible.

Will be harder for the user-defined fits, but very possible for ones I've set up (like single temperature apec).

User supplied images/exposure maps could be binned differently from how XGA creates them

Not really an immediate issue, but there should be some checking for resolutions other than the XCS/XGA standard of 512x512, otherwise things like PSFgen are going to get really confused. - If there are not XCS resolution then they should probably only be used to read in region files.

Use Union type hint for function parameters with multiple allowed input types

What it says in the title really. I didn't know Union typing existed in Python before, so I only stuck with one type

RateMap background subtraction

For the ratemaps to be of any real use in proper analysis (such as measuring accurate gas masses), I'm going to have to deal with the background. The simplest solution is simply to define a background annulus (which I have done), average it and subtract from the ratemap - this is probably what should be implemented first.

Otherwise, I know it is possible to generate some more sophisticated background subtracted images with eSAS, so maybe I'll look more into that.

Need to implement HPC support

This is a planned feature which is essential for making this module viable for large scale analyses. Doesn't matter on the short term, but certainly should be done for beta release.

What I suspect will happen is I'll have to look out for the correct log files to know when a job has finished - then I can read the log files in as stdout and stderr and carry on much like I do already.

Nuisance point sources aren't necessarily all removed from background region

If a nuisance point source is detected as interloping in a cluster, it won't be removed from a background region. So if a point source is right on the edge of a region and present in both the source and background regions, it won't be properly removed.

I've written the DocStrings a bit wrong

I think I'll have to do a pass over everything and write the docstrings better.

Definitely haven't written the return parts right for a lot of them.

Quantity unit checking in some methods fails for equivalent, but not identical, units

For instance, if a method wants an angle argument, and the calculation expects it to be degrees, some methods (like _setup_new_region), will fail if the argument is not in degrees, even though an argument in arcseconds would be totally fine.

GalaxyCluster object needs to take common measurements as arguments

The GC class needs to take common cluster observable measurements as arguments (such as richness, Y, WL mass etc.)

This will make the cluster source objects more useful when it comes to making scaling relations, as well as giving some multi-wavelength context to the X-ray observables.

Novel binning methods - Sanders et al.

Sanders et al papers from the early 2000s had some really nice binning methods for ratemaps in them.

Could be really good to implement them in XGA, because the regions could then be patched in to the spectrum generating functionality to measure more localised temperature maps.

Also binning like that could allow me to come up with a better deprojection method that takes into account localised temperature and density changes better.

Have to make the README more useful

Before I can expose this to other people, I have to make the README more than a collection of random thoughts.

Have to explain how to install it, what the non-python software requirements are, and how to setup the config file.

I should also explain from scratch how to setup an SFTP link for those not working on Apollo.

And have some disclaimers about what doesn't work yet.

Fitting radial profiles

Have to be able to fit a model to the radial profiles to then integrate and retrieve the gas mass.

Plenty of examples in literature, but we'll see if I can improve on them at all.

Read in PSF products if load_products is True

On declaration of a source, XGA can hunt round for products it generated previously. I need to expand it to support PSFs, and eventually the deconvolved image files I will make.

Even if use_peak is false in ExtendedSource declaration, still attempts to find peaks

The custom regions/overdensity regions may not use it as the centre, but the peak finder is still run.

Think this behaviour is wrong, need to correct it.

Create some 'how things work internally' notebooks

Not for everything, just for the key or 'new' bits. The motivation of this is mainly because this is a science module, and one of the major points of making this code public and installable is so it can be understood and appraised/criticised by anyone in the community.

For instance I should do one on the hierarchical clustering peak finding method I came up with.

Also for the PSF deconvolution stuff I'm planning (wrt the rebinning etc), assuming this works at all and makes it into the module.

Split up the big files

The structure of the module could be made a little better, in terms of the internal file structure. I've already split up the product classes into more sensible files and chunks, I should probably do the same for the rest of the project.

Sometimes image data is None

For some reason image data is coming out as None when I try and plot it, for some images.

However a ratemap made of the same image seems to work fine?