Centre and Angles via Numerical Nous Using Bayesian Inference
CANNUBI.py
In real life, Cannubi is one of the most renown fields where the Barolo wine is produced.
Cannubi.py is a python script that uses 3D-Barolo to produce mock galaxies with the aim to derive the geometry of the disc.
Aim
​
Cannubi is a python script to estimate geometrical parameters of a rotating disc using a 3D modelling of the emission and a Markov chain Monte Carlo (MCMC) approach.
The kinematics of a rotating disc (of a galaxy) is derived by fitting a tilted-ring model. Given the large number of free parameters, it is customary to fix the geometrical parameters before fitting the kinematic ones (rotation velocities and velocity dispersions) or at least to have good first estimates.
The geometrical parameters that Cannubi can handle are:
1) coordinates of the centre of the rings;
2) inclination angle (orientation along the line of sight, 90 degrees = edge-on);
3) position angle (anticlockwise angle between the North and major axis of the receding side);
4) thickness of the disc;
5) radial extent of the disc.
Cannubi determined the geometrical parameters of a rotating disc, in particular the centre and the inclination, using the emission in the total flux map or the emission in the whole datacube (beta version). The code uses 3D-Barolo to create mock data cubes and total flux maps obtained from model with different parameters and minimises the residuals between models and data. Before producing the mock maps the model is convolved with the observational PSF.
Output of Cannubi: comparison between iso-surface brightness contours in the total flux map of the data (in this case an artificial galaxy) shown in black and the same contours from the best-fit model.
​
In this case all the parameters of the original artificial galaxy are recovered within the errors.
Download and Getting started
You can download the last version of Cannubi from here: cannubi.py
In order to run Cannubi you need to have:
- 3D-Barolo (1.4 or later)
- python 3 (should also work with 2.7)
- numpy, scipy, astropy, emcee, corner.
Examples
Three mock galaxies with the respective example parameter files can be found here:
From the command line run Cannubi as: python cannubi.py [parameter file.par]
Parameters​
Cannubi uses an input parameter file which is, for the first part very similar to 3D-Barolo .par file
For these parameter please refer to the 3D-Barolo's documentation
The cannubi-specific parameters are:
​
runMCMC [True]:
If Yes/True performs MCMC fit of the parameters. The default parameters that are going to be sampled are inclination, x0 and y0. The user can however specify other parameters to be included, see below.
If No/False performs a grid exploration of a 2D parameter space, of which the first parameter is always the inclination while the second can be specified, see below.
​
fitPA [False]:
If Yes/True the position angle is added to the parameters to be fitted with the MCMC.
​
fitZ0 [False]:
If Yes/True the thickness of the gaseous disc (assumed constant with radius) is added to the parameters to be fitted with the MCMC.
​
fitSep [False]:
If Yes/True the separation between the annuli (RADSEP parameter in 3D-Barolo) is added to the parameters to be fitted with the MCMC.
If this parameter is not fitted then it is highly recommended to fix is by exploring the 2D parameter space with (incl, radsep) and determining the location of the minimum in this space. To derive the correct value of the inclination of a disc it is important to have a good value for RADSEP as the outer parts of the disc are critical.
​
useAllCube [False]:
If Yes/True the residuals between model and data are minimised in the whole datacube.
If No/False only the total-flux map is used.
The option with the total-flux map is much faster (typically 1-2 orders of magnitude) and currently better tested.
​
fitkin [False]:
If Yes/True the kinematics is fitted (by running 3D-Barolo with 3DFIT=True) at each iteration.
If No/False the kinematics is fitted only at the first iteration and kept fixed in the others.
When one uses the total flux map to estimate the inclination (parameter useAllCube), given that this does not depend on the kinematics, this parameter should be False, which dramatically increases the speed of the calculation.
​
trustGuess [False]:
If Yes/True the guesses for inclination, x0 and y0 provided by the user are used (fixed) in the first fit.
If No/False the estimates of these three parameters are made by 3D-Barolo, see its documentation.
normaliseMod [True]:
If Yes/True Cannubi normalises the model using the brightest (parameters fromTop) pixels in the total flux map.
​
useMaskTot [True]:
If Yes/True the residuals in the total flux map are calculated within a mask derived from the mask that 3D-Barolo generates at the beginning (using smooth&search and the user-specified input parameters).
If No/False the residuals are calculated within a contour provided by the user (parameter lastContTot)
​
useMaskMod [Union]:
This keyword handles the type of mask to be built for the model, options are as follows.
Union: every pixel where the data are masked but there is significant emission in the model are used to calculate the residuals. The significance is taken at a threshold equal to the lowest value in the data inside the mask or at 0.5*noise per channel, respectively for useAllCube no and yes.
Same: the model is mask with the same mask used for the data.
Threshold: the model is masked by setting to zero all values below a contour specified by the user.
When useAllCube is False this is the last contour displayed in the total-flux map.
When useAllCube is true this is set by another parameter:
ThresholdAll []: in units of the r.m.s. noise
​
residualType [abs]:
The way the residuals are calculated.
Abs: sum of the absolute values of the differences data - model.
Chisq: sum of the squares of the differences.
​
errorType [marasco]:
The way uncertainties are calculated (the sigma value in the likelihood).
Marasco: the total flux map is rotated 180 degrees and subtracted from the original, the averaged absolute difference is used as uncertainty. This should take into account that the model cannot reproduce the data at a level better than the asymmetries present in the data. This is analogous to the approach in Marasco et al. 2019, 631, 50.
Dof: degrees of freedom, sometime to be preferred if residualType = chisq.
​
cutMask []:
Minimum number of channels that must be present in the mask to build a final totmask (for the total-flux map), typical value = 2.
​
fromTop []:
Percentage of pixel from the brightest one in the data to normalise the total flux map (between model and data), typical value = 30, meaning 30% from the top.
​
lastContTot []:
Lower contour in the total map, if known. This parameter is used only if useMaskMod=none.
​
noiseTimes []:
This is only used to plot the total flux map.
Cannubi tries to estimate the noise in this total map automatically and then plots the lowest contour at a noise level equal to this parameter, ideally between 2 and 3.
​
deltaCentre []:
Excursion of the flat prior of the MCMC in pixels around the initial value of the centre, typically = 5, i.e. pixels, for a well sampled datacube.
​
inclMin []:
Lowest value in the flat prior of the inclination (5-30)
inclMax []:
Highest value in the flat prior of the inclination (80-85)
​
z0min []:
Lowest value in the flat prior of the thickness (in arcsec)
inclMax []:
Lowest value in the flat prior of the thickness (in arcsec)
​
radsepMin []:
Multiplicative factor for the input value of RADSEP that defines the lowest value in the flat prior (0.8)
radsepMax []:
Multiplicative factor for the input value of RADSEP that defines the highest value in the flat prior (1.3)
​
nincl []:
Number of inclination bins if runMCMC = no (grid exploration of the 2D parameter space)
secondPar []:
Second parameter in the parameter space exploration, it can be: radsep, x0 or y0.
nSecondPar []:
Number of bins in the second parameter if runMCMC = no (grid exploration of the 2D parameter space).
​
nwalkers []:
If rumMCMC = yes, number of walkers.
runs []:
If rumMCMC = yes, number of iteration.
The user has the option to set this to "convergence", which is a way to try to make sure that the chains are converging. This is based on the calculation of the autocorrelation time, see here.
​