CAMB sources
-
- Posts: 27
- Joined: November 20 2007
- Affiliation: Cambridge
CAMB sources
Hi,
I ran the latest versions of CAMB sources straight out the of the box and I and I am trying to understand where the 21 cm spectra are printed (and if they are printed at al?)
In the params_21cm file I set Scalars = T as it says:
#for broad window and other sources use get_scalar_cls=T
It does print the power spectra for a sharp window if you set transfer to T.
I am not sure about this verion, but i think it used to work (used it in 2013 and it worked then). According to the website:
scalCls file
Angular power spectra are also output to output_root_scalCls.dat. The columns are in the more confusing but less redundant form as used by earlier versions of CAMB sources:
l C_TT C_EE C_TE C_Phi C_PhiT C_win_1 C_win_2 ... C_Win_T1 C_Win_T2 .... C_win1_win2 ...
Here all C_l are l(l+1)C_l/2pi except for C_Phi and C_PhiT which are C_Phi= l^4 C_l_phi where C_l_phi is the (CMB) lensing potential power spectrum and C_PhiT = l^3 C_l_phiT. So for two window functions, columns 7 and 8 are their auto-spectra, and 11 is their cross-spectrum.[/quote]
which suggests column # 7 is the first 21cm spectrum. However there is no 7th column and when I add a window there is 8th either. So clearly the file is simply not printing the 21cm spectra.
Am I missing something?
Any help would be appreciated.
Cheers
Daan[/quote]
I ran the latest versions of CAMB sources straight out the of the box and I and I am trying to understand where the 21 cm spectra are printed (and if they are printed at al?)
In the params_21cm file I set Scalars = T as it says:
#for broad window and other sources use get_scalar_cls=T
It does print the power spectra for a sharp window if you set transfer to T.
I am not sure about this verion, but i think it used to work (used it in 2013 and it worked then). According to the website:
scalCls file
Angular power spectra are also output to output_root_scalCls.dat. The columns are in the more confusing but less redundant form as used by earlier versions of CAMB sources:
l C_TT C_EE C_TE C_Phi C_PhiT C_win_1 C_win_2 ... C_Win_T1 C_Win_T2 .... C_win1_win2 ...
Here all C_l are l(l+1)C_l/2pi except for C_Phi and C_PhiT which are C_Phi= l^4 C_l_phi where C_l_phi is the (CMB) lensing potential power spectrum and C_PhiT = l^3 C_l_phiT. So for two window functions, columns 7 and 8 are their auto-spectra, and 11 is their cross-spectrum.[/quote]
which suggests column # 7 is the first 21cm spectrum. However there is no 7th column and when I add a window there is 8th either. So clearly the file is simply not printing the 21cm spectra.
Am I missing something?
Any help would be appreciated.
Cheers
Daan[/quote]
-
- Posts: 1943
- Joined: September 23 2004
- Affiliation: University of Sussex
- Contact:
Re: CAMB sources
Can you use the newer scalCovCls output file?
-
- Posts: 27
- Joined: November 20 2007
- Affiliation: Cambridge
CAMB sources
Hi Antony
Thanks. That will probably work.
On a completely different note, I added modules to the SeperableBispectrum.f90 module to include equilateral, orthogonal and folded bispectra and updated the Fisher code to compute the error bars and allow for multiple experiments to make forecasts.
However, I have been struggling with one existing piece of code that does the Fisher forecasting. It computes the error bar as a function of \ell_max and I have introduced (a partly exitsing) loop to run over \ell_max. However, this has a very persistent bug. When computing the 1 sigma error as a function of \ell_max, for very large \ell_max, the error bars increase! (so for low \ell it behaves as expected, for high ell it will turn over at some point; this is true for all shapes, including local which was already there; the effect however is biggest for equilateral, where the error bars for example increases by 20 percent when going from \ell_max = 2500 to 3000). I have been trying to understand this, since one is dealing with a positive definite quantity inside the loop. So changing \ell_max to large values should aways improve the bounds.
I could only think of 2 possibilities:
1) there is a non-positive piece somewhere. To check, I absoluted all quantities. This made no difference.
2) I parallelized the loop and maybe I made a mistake; I checked carefully, but could not find any; I compiled the whole code without opnmp and the results is almost the same (not exactly, which is already worrying, only in the 4th digit though).
There seems to be something inherently off about the Fisher loop since again, it is just receiving a positive definite quantity; so even if I made a mistake in the computation of equilateral for example (where the effect is largest), I still do not see how this could ever happen.
Maybe you have an idea. Here is my bispectrum part:
and here is the Fisher part (most of this was already there). Just added a loop and some options.
Thanks!
Daan
Thanks. That will probably work.
On a completely different note, I added modules to the SeperableBispectrum.f90 module to include equilateral, orthogonal and folded bispectra and updated the Fisher code to compute the error bars and allow for multiple experiments to make forecasts.
However, I have been struggling with one existing piece of code that does the Fisher forecasting. It computes the error bar as a function of \ell_max and I have introduced (a partly exitsing) loop to run over \ell_max. However, this has a very persistent bug. When computing the 1 sigma error as a function of \ell_max, for very large \ell_max, the error bars increase! (so for low \ell it behaves as expected, for high ell it will turn over at some point; this is true for all shapes, including local which was already there; the effect however is biggest for equilateral, where the error bars for example increases by 20 percent when going from \ell_max = 2500 to 3000). I have been trying to understand this, since one is dealing with a positive definite quantity inside the loop. So changing \ell_max to large values should aways improve the bounds.
I could only think of 2 possibilities:
1) there is a non-positive piece somewhere. To check, I absoluted all quantities. This made no difference.
2) I parallelized the loop and maybe I made a mistake; I checked carefully, but could not find any; I compiled the whole code without opnmp and the results is almost the same (not exactly, which is already worrying, only in the 4th digit though).
There seems to be something inherently off about the Fisher loop since again, it is just receiving a positive definite quantity; so even if I made a mistake in the computation of equilateral for example (where the effect is largest), I still do not see how this could ever happen.
Maybe you have an idea. Here is my bispectrum part:
Code: Select all
do il1= 1, SampleL%l0
l1 = SampleL%l(il1)
bi_ix=0
do l2= max(lmin,l1), lmax
min_l = max(abs(l1-l2),l2)
if (mod(l1+l2+min_l,2)/=0) then
min_l = min_l+1
end if
max_l = min(lmax,l1+l2)
do field1=1,nfields
do field2=1,nfields
!PDM 2016: adding other shapes
if(shape == shape_local) then
tmp1 = 2*term*(res_l(l1,1,field1)*resP_l(l2,1,field2) + &
res_l(l2,1,field2)*resP_l(l1,1,field1))
tmp2 = 2*term*(resP_l(l1,1,field1)*resP_l(l2,1,field2))
elseif(shape == shape_equil) then
!extra 2 because of equileteral normalization
tmp1 = -6*term*(res_l(l1,1,field1)*res_l(l2,4,field2) + &
res_l(l2,1,field2)*res_l(l1,4,field1))
tmp2 = -6*term*(res_l(l1,4,field1)*res_l(l2,4,field2))
!multiplicity 2 (since it can be T or E). All terms are the same
!so no two mirrored terms as in local as above
!tmp3 already used
tmp4 = -12*term*(res_l(l1,3,field1)*res_l(l2,3,field2))
tmp5 = 6*term*(res_l(l1,2,field1)*res_l(l2,3,field2) + &
res_l(l2,2,field2)*res_l(l1,3,field1))
tmp6 = 6*term*(res_l(l1,3,field1)*res_l(l2,4,field2) + &
res_l(l2,3,field2)*res_l(l1,4,field1))
tmp7 = 6*term*(res_l(l1,2,field1)*res_l(l2,4,field2) + &
res_l(l2,2,field2)*res_l(l1,4,field1))
elseif(shape == shape_flat) then
tmp1 = 6*term*(res_l(l1,1,field1)*res_l(l2,4,field2) + &
res_l(l2,1,field2)*res_l(l1,4,field1))
tmp2 = 6*term*(res_l(l1,4,field1)*res_l(l2,4,field2))
tmp4 = 18*term*(res_l(l1,3,field1)*res_l(l2,3,field2))
tmp5 = -6*term*(res_l(l1,2,field1)*res_l(l2,3,field2) + &
res_l(l2,2,field2)*res_l(l1,3,field1))
tmp6 = -6*term*(res_l(l1,3,field1)*res_l(l2,4,field2) + &
res_l(l2,3,field2)*res_l(l1,4,field1))
tmp7 = -6*term*(res_l(l1,2,field1)*res_l(l2,4,field2) + &
res_l(l2,2,field2)*res_l(l1,4,field1))
elseif(shape == shape_ortho) then
tmp1 = -18*term*(res_l(l1,1,field1)*res_l(l2,4,field2) + &
res_l(l2,1,field2)*res_l(l1,4,field1))
tmp2 = -18*term*(res_l(l1,4,field1)*res_l(l2,4,field2))
tmp4 = -48*term*(res_l(l1,3,field1)*res_l(l2,3,field2))
tmp5 = 18*term*(res_l(l1,2,field1)*res_l(l2,3,field2) + &
res_l(l2,2,field2)*res_l(l1,3,field1))
tmp6 = 18*term*(res_l(l1,3,field1)*res_l(l2,4,field2) + &
res_l(l2,3,field2)*res_l(l1,4,field1))
tmp7 = 18*term*(res_l(l1,2,field1)*res_l(l2,4,field2) + &
res_l(l2,2,field2)*res_l(l1,4,field1))
endif
do field3=1,nfields
Bispectrum => Bispectra(field1,field2,field3,fnl_bispectrum_ix)
bix=bi_ix
do l3=min_l,max_l ,2
bix=bix+1
!PDM 2016: adding other shapes
if(shape == shape_local) then
Bispectrum%b(bix,il1) = Bispectrum%b(bix,il1) + &
(tmp1*resP_l(l3,1,field3) + tmp2*res_l(l3,1,field3))
else!equilateral, orthogonal and flat
Bispectrum%b(bix,il1) = Bispectrum%b(bix,il1) + &
(21)**2*32*pi*(tmp1*res_l(l3,4,field3) + tmp2*res_l(l3,1,field3) + &
tmp4*res_l(l3,3,field3) + tmp5*res_l(l3,4,field3) + &
tmp6*res_l(l3,2,field3) + tmp7*res_l(l3,3,field3))
endif
end do
end do
end do
end do
bi_ix=bix
end do !l2
end do !il1
Code: Select all
do lmaxcuti=1, SampleL%l0
if(Fisher_lmin) then
lstart = 2
!lstart = 2 + (lmaxcuti-1)*4
ellmaxl = SampleL%l0
ellminl = lmaxcuti
lmax= SampleL%l(SampleL%l0)
else
!can also change lstart not to include lowest
!multipole
lstart = 2
if (SampleL%l(lmaxcuti) < 425) cycle
ellmaxl = lmaxcuti !SampleL%l0
!ellmaxl = SampleL%l0
ellminl = 1
!SampleL%l0 = SampleL%l0 - 10
lmax= SampleL%l(ellmaxl)
endif
!BispectrumParams%fskyE(2) = BispectrumParams%fskyE(2) + 0.01
!BispectrumParams%fskyT(2) = BispectrumParams%fskyT(2) + 0.01
Noise(1:2) = BispectrumParams%FisherNoise(1:2)/ (COBE_CMBTemp*1e6)**2 !Planckish, dimensionless units
NoiseP(1:2) = BispectrumParams%FisherNoisePol(1:2)/ (COBE_CMBTemp*1e6)**2
Noise(1:2) = BispectrumParams%fskyT(1:2)*Noise(1:2)
NoiseP(1:2) = BispectrumParams%fskyE(1:2)*NoiseP(1:2)
do i=lmin,lmax
if (CP%DoLensing) then
cl(:,i) = CL_lensed(i,1,CT_Temp:CT_Cross)
else
cl(1,i) = CL_Scalar(i,1,C_Temp)
cl(2,i) = CL_Scalar(i,1,C_E)
cl(4,i) = CL_Scalar(i,1,C_Cross)
cl(3,i) = 0
end if
if (CP%WantTensors .and. i<= CP%Max_l_tensor .and. i>=2) then
cl(:,i) = cl(:,i) + Cl_tensor(i,1,CT_Temp:CT_Cross)
end if
end do
if (.false.) then
call OpenTxtFile('CAMBdefault_lensedCls.dat',3)
do i=lmin,lmax
!Assume T,E,B,X ordering
read(3,*) j, cl(1:4,i)
if (j<lmin) read(3,*) j, cl(1:4,i)
cl(:,i)=cl(:,i)/(COBE_CMBTemp*1e6)**2
end do
close(3)
end if
if (Noise(1) >0) then
file_tag = concat(file_tag,'_noise')
end if
xlc= 180*sqrt(8.*log(2.))/3.14159
sigma2(1:2) = (BispectrumParams%FisherNoiseFwhmArcmin(1:2)/60/xlc)**2
write(*,*) "sigma2:", sigma2
allocate(InvC(lmax))
do l1= lmin, lmax
tmp = l1*(l1+1)/(2*pi)
!PDM 2016:
!allow to combine 2 different experiments (e.g. Planck and CMB-S4)
if(l1 .le. BispectrumParams%EllTrans(1)) then
Cl(1,l1) = Cl(1,l1)/tmp + Noise(1)*exp(l1*(l1+1)*sigma2(1))
else
!PDM 2016: factor 2.5 is to make Noise = 1 muK for fsky = 0.4
Cl(1,l1) = Cl(1,l1)/tmp + 2.5*Noise(2)*exp(l1*(l1+1)*sigma2(2))
endif
if(l1 .le. BispectrumParams%EllTrans(1)) then
Cl(2:3,l1) = Cl(2:3,l1)/tmp + NoiseP(1)*exp(l1*(l1+1)*sigma2(1))
else
Cl(2:3,l1) = Cl(2:3,l1)/tmp + 2.5*NoiseP(2)*exp(l1*(l1+1)*sigma2(2))
endif
Cl(4,l1) = Cl(4,l1)/tmp
allocate(InvC(l1)%C(nfields,nfields))
if (nfields > 2) stop 'Not implemented nfields>2 in detail'
if (nfields==1) then
InvC(l1)%C(1,1)=(2*l1+1)/cl(1,l1)/InternalScale
else
InvC(l1)%C(1,1)=cl(2,l1)
InvC(l1)%C(1,2)=-cl(4,l1)
InvC(l1)%C(2,1)=-cl(4,l1)
InvC(l1)%C(2,2)=cl(1,l1)
InvC(l1)%C= InvC(l1)%C * (2*l1+1)/(cl(1,l1)*cl(2,l1)-cl(4,l1)**2)/InternalScale
end if
end do
print *,'Getting Fisher for lmax = ', lmax
if (debugMsgs) starttime=GetTestTime()
!mysum = 0
ifish_contribs=0
tmp = 0
tmp1 = 0
tmp2 = 0
tmpf = 0
!lmax = SampleL%l(59)
!$OMP PARAllEl DO DEFAUlT(SHARED),SCHEDULE(STATIC,3) &
!$OMP PRIVATE(il1,l1,l2,l3,fish_l1,bi_ix,min_l,max_l,a3j_00,a3j), &
!$OMP PRIVATE(Bispectrum,Bispectrum2,minl2,bix,tmp,tmp1,tmp2,tmpf), &
!$OMP PRIVATE(field1,field2,field3,f1,f2,f3,bispectrum_type,bispectrum_type2)
!!$OMP REDUCTION(+:fish_l1)
do il1= ellminl, ellmaxl
!do il1= 58, 59
!write(*,*) lmaxcuti
allocate(fish_l1(nbispectra,nbispectra,nfields,nfields)) !last indices are field1,f1
l1 = SampleL%l(il1)
if (l1< lstart) cycle
fish_l1=0
bi_ix=0
do l2 = l1,lmax
if (l2< lstart) cycle
min_l = max(lstart,max(abs(l1-l2),l2))
if (mod(l1+l2+min_l,2)/=0) then
min_l = min_l+1
end if
max_l = min(lmax,l1+l2)
call GetThreeJs(a3j(abs(l2-l1)),l1,l2,0,0)
do l3=min_l,max_l ,2
a3j_00(l3)=a3j(l3)**2
end do
tmp1= 1.d0/(4*pi) !(2l+1) factors included in InvC
minl2=min_l
bix=bi_ix
do field1=1,nfields
do f1=1,nfields
tmpf(1)= InvC(l1)%C(field1,f1)*tmp1
do field2=1,nfields
do f2=1,nfields
tmpf(2)= InvC(l2)%C(field2,f2)*tmpf(1)
do field3=1,nfields
do bispectrum_type=1,nbispectra
if (bispectrum_type==lens_bispectrum_ix) then
Bispectrum=>SqueezedLensingKernel(field2,field3)
else
Bispectrum=>bispectra(field1,field2,field3,bispectrum_type)
end if
do f3=1,nfields
do bispectrum_type2=bispectrum_type,nbispectra
if (bispectrum_type2==lens_bispectrum_ix) then
Bispectrum2=>SqueezedLensingKernel(f2,f3)
else
Bispectrum2=>Bispectra(f1,f2,f3,bispectrum_type2)
end if
min_l=minl2
bi_ix=bix
if (min_l==l2) then
!Symmetry factors
bi_ix=bi_ix+1
l3=l2
if (l2==l1) then
!l1=l2=l3
tmp = Bispectrum%b(bi_ix,il1)*tmpf(2)*Bispectrum2%b(bi_ix,il1) &
*InvC(l3)%C(field3,f3)*a3j_00(l3)/6
else
!l3=l2 (l3=l1<>l2 can't happen because l1<=l2<=l3)
tmp = Bispectrum%b(bi_ix,il1)*tmpf(2)*Bispectrum2%b(bi_ix,il1) &
* InvC(l3)%C(field3,f3)*a3j_00(l3)/2
end if
min_l = min_l+2
else
tmp=0
end if
tmp2=0
do l3=min_l,max_l ,2
bi_ix=bi_ix+1
tmp2 = tmp2 + Bispectrum%b(bi_ix,il1)*Bispectrum2%b(bi_ix,il1) &
* InvC(l3)%C(field3,f3)*a3j_00(l3)
end do
if (l2==l1) then
tmp2=tmp2*tmpf(2)/2
else
tmp2=tmp2*tmpf(2)
end if
fish_l1(bispectrum_type,bispectrum_type2,field1,f1)= &
fish_l1(bispectrum_type,bispectrum_type2,field1,f1)+(tmp+tmp2)
end do !bispectrum_type2
end do
end do !bispectrum_type
end do !field3
end do !f2
end do !field2
end do !f1
end do !field1
end do !l2
ifish_contribs(il1,:,:,:,:)=fish_l1(1:nbispectra,1:nbispectra,:,:) /InternalScale
deallocate(fish_L1)
end do
!$OMP END PARAllEl DO
Daan
-
- Posts: 1943
- Joined: September 23 2004
- Affiliation: University of Sussex
- Contact:
Re: CAMB sources
Hmm, don't know. At least if you have a clear contradiction it should be easy enough to pin download to which L/which term negative contributions are coming from.