Prof Richard Massey
Dark Matter cartographer

Module to include COSMOS 3D weak lensing data in CosmoMC

Description Download Installing the new modules in your version of CosmoMC Further modifications for simultaneously using weak lensing and Lyman-α data Referring to this code

Description

This page describes the installation of the modules WeakLen (file 'WeakLen.f90') and Angular (file 'AngularCorr.f90'), written by Julien Lesgourgues (LAPTH, Annecy), to be used with the public MCMC code CosmoMC by Sarah Bridle and Antony Lewis.

The modules WeakLen and Angular extend CosmoMC to incorporate weak lensing (cosmic shear) data into the analysis. The modules were written for the paper "A combined analysis of Lyman-α forest, 3D Weak Lensing and WMAP year three data" by Lesgourgues et al., JCAP 11 (2007) 8. They are designed and optimized for the HST COSMOS data released in "COSMOS: 3D weak lensing and the growth of structure" by Massey et al., ApJS 172 (2007) 239. The COSMOS data is included. However, the modules could also be easily adapted to other weak lensing data.

If it is a while since the papers were published, it might be worth checking with the authors for any subsequent updates to code or lensing data.

Download

Download cosmos.tar.gz. If you expand it inside your cosmomc/ directory via "tar -zxvf cosmos.tar.gz" you will get the files

 source/WeakLen.f90
 source/AngularCorr.f90
 data/cosmos_ctheta.txt
 data/cosmos_eta.txt
 data/cosmos_invcov.txt

You may be asked whether you want to overwrite the existing file source/WeakLen.f90. In the current distribution of CosmoMC (October 2006), this existing file is essentially empty: you can safely overwrite it. We have not included documentation for the fortan code, but the files contain embedded explanatory comments.

Installing the new modules in your version of CosmoMC

If you have CosmoMC installed (we refer here to the version of October 2006), you have to make a few simple modifications in order to make it compatible with the two new source files source/WeakLen.f90 and source/AngularCorr.f90:

• The code must be aware of the presence of three new nuisance parameters A, B, C. These describe the systematic errors in COSMOS weak lensing data — see JCAP 11 (2007) 8. To enable this you should:
  • increase the parameter num_norm in source/settings.f90 by 3. In the public code, the default value is 3 and should be increased to 6:
    integer, parameter :: num_norm = 6
  • give the names norm_A, norm_B, norm_C to the index of the three new parameters. In source/cmbtypes.f90, just after the line:
    integer, parameter :: norm_As=1, norm_amp_ratio=2
    add a new line:
    integer, parameter :: norm_A=4, norm_B=5, norm_C=6
  • define the new parameters and their range in params.ini: at the very end of this file, add:
# nuisance param A (COSMOS systematics)
param14 = 1. 0. 2. 0.06 0.06
# nuisance param B (COSMOS systematics)
param15 = 1. 1. 2. 0.05 0.05
# nuisance param C (COSMOS systematics)
param16 = 1. 1. 2. 0.10 0.10
    Note that B and C must be greater than one by definition, hence the lower limit 1 in the previous lines.
  • if you use a covariance matrix, don't forget that its dimension will increase by 3 due to the new parameters. Also in distparams.ini the number of columns and parameters for the covariance matrix should be increase.
• The code should know that there are new modules being used. This requires minimal modifications since the public version of CosmoMC already includes a sketch of module called Weaklen. The few additions remaining to be done are:
  • in source/CMB_Cls_simple.f90, just after the line:
    logical :: use_nonlinear = .false.
    add a new line:
    logical :: Use_WeakLen = .false.
  • in source/Makefile, in the definition of OBJFILES, add AngularCorr.o betwwen cmbdata.o and WeakLen.o
  • in params.ini, together with the lines:
    use_CMB = T
    ...
    use_min_zre = 0
    add a new line:
    use_WeakLen = T
    set it to T/F when you want/don't want to include this data set.
    IMPORTANT: if use_WeakLen = T don't forget to include non-linear corrections by setting also nonlinear_pk = T, otherwise your results will be incorrect.
• CosmoMC should ask CAMB to compute the matter power spectrum P(k,z) at various redshifts between 0 and 3, until higher values of k than in the default code. For this you should:
  • open source/cmbtypes.f90 and change the value of three parameters at the beginning. the recommended values are:
    integer, parameter :: num_matter_power = 98 ! default was 74; 98 is fine for COSMOS
real, parameter :: matter_power_minkh = 0.999e-4 ! unchanged
real, parameter :: matter_power_dlnkh = 0.143911568 ! unchanged
real, parameter :: matter_power_maxz = 3. ! default was 0; 3 is fine for COSMOS
integer, parameter :: matter_power_lnzsteps = 10 ! default was 1; 10 is fine for COSMOS
  • open source/CMB_Cls_simple.f90 and change the value of one parameter:
    P%Transfer%kmax = 100. ! default value was 1.2; 100 is OK for COSMOS
• The code should store somewhere the values of comobile distance r and Hubble rate H=dz/dr for each redshift were P(k,z) is computed: the knowledge of r(z) and H(z) is necessary for computing the shear correlation functions. For this you should implement the following changes:
  • define a new array z_and_r(:,:) where these values will be stored. In source/cmbtypes.f90, inside the definition of the type Type CosmoTheory, insert new lines:
    Type CosmoTheory
...
real z_and_r(3,matter_power_lnzsteps)
end Type CosmoTheory
  • at the begining of source/CMB_Cls_simple.f90, in the line starting by use CAMB, only :, add another function DeltaTime to the list:
    use CAMB, only : CAMB_GetResults, CAMB_GetAge, CAMBParams, CAMB_SetDefParams,Transfer_GetMatterPower, &
AccuracyBoost, Cl_scalar, Cl_tensor, Cl_lensed, outNone, w_lam, &
CAMBParams_Set,Re_OpticalDepthAtZ, MT, CAMBdata, NonLinear_Pk,&
CAMB_GetTransfers,CAMB_FreeCAMBdata,CAMB_InitCAMBdata, CAMB_TransfersToPowers, &
initial_adiabatic,initial_vector,initial_iso_baryon,initial_iso_neutrino, initial_iso_neutrino_vel, &
DeltaTime
  • still in source/CMB_Cls_simple.f90, in the subroutine subroutine GetCls(CMB,Info, Cls, error), at the end of the declaration of variables, add:
    real*8 aaa
real*8 dtauda
external dtauda
  • still in the same subroutine but a few lines later, just after the lines:
    if (Feedback > 1) write (*,*) 'CAMB done'
    add the new group of lines:
    if (Use_WeakLen) then
  do zix = 1,matter_power_lnzsteps
   Info%Theory%z_and_r(1,zix)=P%Transfer%redshifts(zix)
   aaa=1.d0/(1.+P%Transfer%redshifts(zix))
   Info%Theory%z_and_r(2,zix)=DeltaTime(aaa,1.d0)
   Info%Theory%z_and_r(3,zix)=1./(aaa**2*dtauda(aaa))
  end do
end if

You are done! In principle you should be able to compile as usual ( cd camb, make clean, make all, cd ../source, make clean, make cosmomc, cd ..) and run with ./cosmomc params.ini. If you start in chatty mode (Feedback = 2) the new module will output some useful information if use_WeakLen = T.

Further modifications for simultaneously using weak lensing and Lyman-α data

When using weak lensing data (use_WeakLen = T) it is crucial to include non-linear corrections by setting also nonlinear_pk = T. However, some other data sets require nonlinear_pk to be set to F. This is the case for current Lyman-alpha data modules, which provide constraints on the linear power spectrum (non linear corrections are accounted for earlier during the analysis pipeline of Lyman-alpha raw data). This raises a contradiction which is not difficult to solve. Our solution is to duplicate various subroutines in camb/module.f90 and source/cmbtypes.f90in order to define two independent functions:
function MatterPowerAt(T,kh) ! non-linear matter power spectrum
function MatterPowerLinearAt(T,kh) ! linear matter power spectrum
Once this is done, the trick is simply to call MatterPowerLinearAt(T,kh) in the Lyman-alpha module(s), and MatterPowerAt(T,kh) in the other modules (weak lensing and galaxy power spectrum). In order to see the list of necessary changes, you can:

• Have a look at our camb/module.f90 and source/cmbtypes.f90 files. Modifications are always indicated by the comment ! LAPTH NEW. In addition:
• in source/CMB_Cls_simple.f90, you should adapt accordingly the following two lines:
call Transfer_GetMatterPower(Info%Transfers%MTrans,&
Info%Theory%matter_power(:,zix),&
Info%Theory%matter_power_linear(:,zix),& !LAPTH NEW: note the new argument
matter_power_lnzsteps-zix+1,&
1,matter_power_minkh, matter_power_dlnkh,num_matter_power)
  and:
call Transfer_GetMatterPower(MT,&
Theory%matter_power(:,zix),&
Theory%matter_power_linear(:,zix),& !LAPTH NEW: note the new argument
matter_power_lnzsteps-zix+1,1,matter_power_minkh,&
matter_power_dlnkh,num_matter_power)

Reference to this code — publishing

If you use the code and publish results, please refer to having used this code and cite the paper JCAP 11 (2007) 8.