!+ !PURPOSE: ! The two following routines allow to account for the CMB lensing observables ! (C_l^Td and C_l^dd) in a CosmoMC parameter estimation. !AUTHOR: ! They are an adaptation of the original ones written by A. Lewis performed ! by L. Perotto and J. Lesgourgues !VERSION: ! Nov 2009: LP transfers the updates of the CosmoMC::ReadAllExact within the FUTURCMB::ReadAllExact ! May 2006: J. Lesgourgues redefines f_sky as the fraction of the observed sky ! (instead of its square root) ! March 2006: written by LP !- subroutine ReadAllExact(Ini,aset) Type(TIniFile) :: Ini Type (CMBdataset) :: aset character(LEN=Ini_max_string_len) :: fname integer l, idum, ncls, ncol integer file_unit !In this case we have data for every l, and use exact full-sky likelihood expression !with some fudge factor fsky_eff^2 to reduce the degrees of freedom: fsky_eff*(2l+1) ! (JL note: we change notations wrt Lewis since "his" fsky_eff was the square root ! of the actual f_sky) aset%num_points = 0 aset%all_l_lmax = Ini_Read_Int_File(Ini,'all_l_lmax') if (aset%all_l_lmax > lmax) stop 'cmbdata.f90::ReadAllExact: all_l_lmax > lmax' allocate(aset%all_l_obs(2:aset%all_l_lmax,num_cls)) if (aset%has_pol) then if (num_cls >= 5) then allocate(aset%all_l_noise(2:aset%all_l_lmax,3)) else allocate(aset%all_l_noise(2:aset%all_l_lmax,2)) endif else if (num_cls == 3) then allocate(aset%all_l_noise(2:aset%all_l_lmax,2)) else allocate(aset%all_l_noise(2:aset%all_l_lmax,1)) endif end if allocate(aset%all_l_fsky(2:aset%all_l_lmax)) fname = trim(Ini_Read_String_File(Ini,'all_l_file')) ncol = TxtFileColumns(fname) file_unit = new_file_unit() call OpenTxtFile(fname,file_unit) !File format: !l C_TT (C_TE C_EE [C_BB] C_dd C_Td) N_T (N_P N_d) fsky_eff !Must have num_cls set to correct value for file do l = 2, aset%all_l_lmax read (file_unit, *, end=100, err=100) idum, aset%all_l_obs(l,:), aset%all_l_noise(l,:), aset%all_l_fsky(l) if (idum /= l) stop 'Error reading all_l_file' end do call CloseFile(file_unit) return 100 stop 'Error reading all_l_file file' end subroutine ReadAllExact function ChiSqExact(cl, aset) !Compute -ln(Likelihood) !real cl(lmax,num_cls) ! on ajoute un test pour verifier l'effet de la normalisation choisie par AL real, dimension(:,:), intent(in) :: cl Type(CMBdataset) :: aset integer l,nell, ncls, nlcls, err, lge real ChiSqExact, determ, determobs, tcmb double precision chisq, CT, CE, CB, Cd, Ctd real dof integer :: num=0 logical :: do_lens num=num+1 lge=100 nlcls = size(cl,2) nell = size(cl,1) if (aset%has_pol) then if (nlcls >= 5 ) then do_lens = .true. else do_lens = .false. endif else if (nlcls == 3 ) then do_lens = .true. else do_lens = .false. endif endif if (do_lens) then ncls = nlcls - 2 else ncls = nlcls endif if (nlcls /= num_cls) then print*,"pas le bon nb de Cls CMB dans le fichier CAMB..." stop endif tcmb = 2.726e6 chisq=0 do l=2, aset%all_l_lmax dof = aset%all_l_fsky(l)*(2*l+1) !Ignoring l correlations but using f_sky_eff fudge factor may be a good approx !JL note: for us fsky is the portion of the sky (e.g. 0.65 for Planck) !while in the old Lewis implementation this fraction was noted as f_sky^2=0.65 !for nearly full sky observations CT = cl(l,1) + aset%all_l_noise(l,1) if (aset%has_pol) then CE = cl(l,3) + aset%all_l_noise(l,2) if (do_lens == .false.) then ! polarisation without lensing determ = CT*CE - cl(l,2)**2 determobs = aset%all_l_obs(l,1)*aset%all_l_obs(l,3) - & & aset%all_l_obs(l,2)**2 chisq = chisq + dof*( & ( CT*aset%all_l_obs(l,3) + CE*aset%all_l_obs(l,1) - & & 2*cl(l,2)*aset%all_l_obs(l,2))/determ & & + log( determ / determobs ) -2) else ! polarisation and lensing Ctd = cl(l,nlcls)/tcmb/twopi*(1.*(l+1)/(l*1.))**1.5 Cd = cl(l,nlcls-1)/tcmb**2/twopi*(1.*(l+1)/(l*1.))**2 + aset%all_l_noise(l,3) determ = CT*CE*Cd - cl(l,2)**2*Cd - Ctd**2*CE determobs = aset%all_l_obs(l,1)* & & aset%all_l_obs(l,3)* & & aset%all_l_obs(l,nlcls-1) - & & aset%all_l_obs(l,2)**2*aset%all_l_obs(l,nlcls-1) - & & aset%all_l_obs(l,nlcls)**2*aset%all_l_obs(l,3) chisq = chisq + dof*( & ( aset%all_l_obs(l,1)*CE*Cd + & & CT*aset%all_l_obs(l,3)*Cd + & & CT*CE*aset%all_l_obs(l,nlcls-1) - & & cl(l,2)**2*aset%all_l_obs(l,nlcls-1) - & & 2*cl(l,2)*aset%all_l_obs(l,2)*Cd - & & Ctd**2*aset%all_l_obs(l,3) - & & 2*Ctd*aset%all_l_obs(l,nlcls)*CE)/determ & & + log( determ / determobs ) - 3) endif if (ncls>3) then !add in BB CB = cl(l,ncls) + aset%all_l_noise(l,2) chisq = chisq + dof * (aset%all_l_obs(l,ncls)/CB & +log(CB/aset%all_l_obs(l,ncls)) - 1) end if else if (do_lens == .false.) then ! temperature without lensing chisq = chisq + dof * (aset%all_l_obs(l,1)/CT & +log(CT/aset%all_l_obs(l,1)) - 1) else ! temperature and lensing Cd = cl(l,2) + aset%all_l_noise(l,2) determ = CT*Cd - cl(l,3)**2 determobs = aset%all_l_obs(l,1)*aset%all_l_obs(l,2) - & & aset%all_l_obs(l,3)**2 chisq = chisq + dof*( & &(CT*aset%all_l_obs(l,2) + Cd*aset%all_l_obs(l,1) & & - 2 *cl(l,3)*aset%all_l_obs(l,3))/determ & & + log( determ / determobs ) -2) endif end if end do ChiSqExact = chisq end function ChiSqExact