The code has been rewritten since December 6th 2024. To tidy things up, the git history of the public repository was deleted. For last version of the 'old' code, please see the Legacy release.
This package can be obtained using pip:
pip install halo-fdca.
See also: https://pypi.org/project/halo-fdca/1.0/
This software is created to automate flux density (and power) estimations of radio halos in galaxy clusters. This is done by fitting the surface brightness profile of halos to a mathematical model using Bayesian Inference. From the resulting fit, the flux density profile can be calculatied analytically. The text below provides a step-by-step example and walk trough of the algorithm as well as a summary of installation and machine requirements. A full text on the methodology can be found under 'Citation'.
This software is open source and still in development. Suggestions, remarks, bugfixes, questions etc. are more than welcome and can be sent to j.boxelaar@ira.inaf.it.
If you make use of this code in its original form or portions of it, please cite:
Boxelaar J.M. et al., A robust model for flux density calculations of radio halos in galaxy clusters: Halo-FDCA, Astronomy and Computing, Volume 35, 2021, 100464, ISSN 2213-1337, https://doi.org/10.1016/j.ascom.2021.100464.
Also available on the arXiv: https://arxiv.org/abs/2103.08554
This software is written and tested in the Python 3.11(!!) and is meant for application in astronomical research. Note: Multiprocessing does not work for python versions below 3.11, in case you use an older vresion of Python, always set num_cpus=1 in the fdca.Fit() function.
Key software requirements to run the pipeline include:
Astropy v4.0 and up
https://www.astropy.org
Scipy v1.4 and up
Numpy v1.19 and up
Matplotlib v3.2 and up
https://www.scipy.org/docs.html
emcee v3.0 and up
https://emcee.readthedocs.io/en/stable/
corner v2.2
https://corner.readthedocs.io/en/latest/index.html
skimage v0.17
https://scikit-image.org
tqdm
https://tqdm.github.io/ (this is not a strict requirement, but it is useful if you want to follow the progress of the run)
This software will run a multithreaded process by default using the multiprocessing module. This requires significant computing power and it is advised to use a dedicated external computing mashine. Multithreading can be turned off but this will increase the run time up to a factor 10 (depending on the settings).
The code is tested on Unix machines using both PYTHON 3.6 and PYTHON 3.8. Please inform us when you run into any issues related toincompatible software versions.
This flowchart gives a general overview of the pipeline. On code level, the pipeline works roughly like this:
First, a Radio_Halo object must be initiated. When initiating the object, all relevant properties are processed. This class is documented in HaloObject.py.
This class also handles the very first fit such that the image coordinates can be related to sky coordinates without using the header. The Radio_Halo class handles steps up to 'Second preliminary fit'. 'Rotating and regridding' is not performed by the class.
The Markov Chain Monte Carlo (MCMC) algorithm is performed in the blue part of the flowchart. This takes place in the fitting class in markov_chain_monte_carlo.py (see class documentation there and below). This class takes a Radio Halo object as input and from there starts the profile fitting based on the extra settings given as input. 'Chains' that are the result of MCMC are saved in new FITS files (./outputpath/Samples/), settings of the speciffic run are saved in the header.
The processing class in mcmc_fitting_multicomponent.py takes a halo object as input and processes the MCMC results by generating figures, statistical analysis and final flux density and parameter estimations. these results are found in the log files.
The code requires very specific input to be able to work properly.
-
.FITS file:
This includes a FITS image with a clean, preferably point source subtracted radio image of a galaxy cluster containing diffuse emission. -
.reg file:
A DS9 region file with regions drawn around contaminating sources/artifacts in the image, saved in physical coordinates. The shape of the region can be of any kind. A region file is not mandatory to run the code.
Currently, cluster characteristics such as location and redshift can be retrieved from catalogues on VIZIER. This automatic retrieval of cluster information is availible for MCXC, PSZ2, Abell and WHL clusters. A catalogue search is not always succesful, it is adviced to give the essential cluster information, redshift and sky location, as input.
In v2.0, Halo-FDCA is used as a package. See develop.py as example of how to run Halo-FDCA. we are currently working on making the package availible through pip. For now, if one wants to use the code as a package, run export PYTHONPATH=/path_to_code/halo_fdca:$PYTHONPATH in your terminal.
Creating a RadioHalo object passing the minimum amount of information:
import halo_fdca as fdca
Halo = fdca.RadioHalo('A2744','Example/Data_dir/A2744_JVLA.image.fits')Creating a Fitting object and start a run:
import halo_fdca as fdca
Halo = fdca.RadioHalo('A2744','Example/Data_dir/A2744_JVLA.image.fits')
# in Fit(), it is possible to specify number of steps and walkers
fit = fdca.Fit(Halo, model='circle')
fit.run()To prevent rerunning the MCMC, one can save a run with fit.save(path) where the path to the saved samples is optional. Load a previous run with
fit = fdca.Fit(Halo, model='circle').load()
# OR
fit = fdca.load('samples.json')where the path points to the json generated by fit.save().
Analyzing results after a fit is run or loaded is done by
fit = fdca.Fit(Halo, model='circle').load()
fit.results.plot()
chi2 = fit.results.get_chi2()
flux, flux_uncertainty = fit.results.get_flux()
power, power_uncertainty = fit.results.get_power()
samples = fit.get_samples()
parameter_names = fit.get_param_names()Samples file in json format containing the walker strings and run information (found in outputpath/Samples). This file is used to process the routine and retrieve flux denisty and parameter values. All values and their one sigma uncertainties will be printed in the log file after running the pipeline. The figures will be in outputpath/Plots. There, all relevant figures will be saved. This includes the corner plot, walker plot and radio image with model overlay for the original and regridded image.
The .log file outputs the final fit information and it looks something like this:
Run information for object A1033:
RMS noise: 0.0002863803820218891 Jy / beam
Model: circle
Walkers: 10
Steps: 100
Burntime: 12
Mask: True
Rebin: True
K_exponent: True
Offset: False
Fit results:
Flux density at 143.7 MHz: 336.81546 mJy +/- 4.79004 mJy
Reduced chi-squared: 4.859490839643051
I0: 25.543 +/- 0.2336 ($\mu$Jy arcsec$^{-2}$)
x0: 157.928 +/- 0.0002 (deg)
y0: 35.040 +/- 0.0002 (deg)
r1: 100.547 +/- 1.1243 (kpc)
Uncertainties (lower, upper):
[2.38880646e-01 1.90126966e-04 1.85981806e-04 1.27793972e+00
5.43166949e-20]
[2.28331040e-01 2.48196633e-04 2.07639059e-04 9.70618808e-01
8.45588410e-20]
The same information can be printed to the terminal by printing the resuts: print(fit.results).
Fitting multiple exponential shapes to an image in now possible, but untested. The ability to add multiple components adds to possibility to link the location of the shapes (usefull for mega halos) or to change the profile of the exponential (e.g. to a Gaussian istead of a exponential., usefull for fitting the contribution from an AGN). As a example, given a RadioHalo, we want to fit a circle model with exponential shape and have a circle model with a Gaussian shape at the same location accounting for the mega halo shape:
fit = fdca.Fit(Halo, model=['circle','circle'], link_loc=[True, True], profiles=["default", "gaussian"])where "default" is the default exponential shape. Note: it is not yet possible to get results from the fits, to analyse the results, use the outputted samples directly:
samples = fit.get_samples()
param_names = fit.get_param_names()