Page 1 of 1

CAMB sources

Posted: July 01 2016
by daan meerburg
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]

Re: CAMB sources

Posted: July 01 2016
by Antony Lewis
Can you use the newer scalCovCls output file?

CAMB sources

Posted: July 05 2016
by daan meerburg
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:

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
and here is the Fisher part (most of this was already there). Just added a loop and some options.

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 &#40;SampleL%l&#40;lmaxcuti&#41; < 425&#41; cycle
             ellmaxl = lmaxcuti !SampleL%l0
             !ellmaxl = SampleL%l0
             ellminl = 1
             !SampleL%l0 = SampleL%l0 - 10
             lmax= SampleL%l&#40;ellmaxl&#41; 
          endif
          !BispectrumParams%fskyE&#40;2&#41; = BispectrumParams%fskyE&#40;2&#41; + 0.01
          !BispectrumParams%fskyT&#40;2&#41; = BispectrumParams%fskyT&#40;2&#41; + 0.01 
          Noise&#40;1&#58;2&#41; = BispectrumParams%FisherNoise&#40;1&#58;2&#41;/ &#40;COBE_CMBTemp*1e6&#41;**2  !Planckish, dimensionless units  
          NoiseP&#40;1&#58;2&#41; = BispectrumParams%FisherNoisePol&#40;1&#58;2&#41;/ &#40;COBE_CMBTemp*1e6&#41;**2
          
          Noise&#40;1&#58;2&#41; = BispectrumParams%fskyT&#40;1&#58;2&#41;*Noise&#40;1&#58;2&#41;
          NoiseP&#40;1&#58;2&#41; = BispectrumParams%fskyE&#40;1&#58;2&#41;*NoiseP&#40;1&#58;2&#41;
          do i=lmin,lmax
             if &#40;CP%DoLensing&#41; then
              cl&#40;&#58;,i&#41; = CL_lensed&#40;i,1,CT_Temp&#58;CT_Cross&#41;
             else
              cl&#40;1,i&#41; = CL_Scalar&#40;i,1,C_Temp&#41;
              cl&#40;2,i&#41; = CL_Scalar&#40;i,1,C_E&#41;
              cl&#40;4,i&#41; = CL_Scalar&#40;i,1,C_Cross&#41;
              cl&#40;3,i&#41; = 0              
             end if
             if &#40;CP%WantTensors .and. i<= CP%Max_l_tensor .and. i>=2&#41; then
               cl&#40;&#58;,i&#41; = cl&#40;&#58;,i&#41; + Cl_tensor&#40;i,1,CT_Temp&#58;CT_Cross&#41;
             end if
          end do
          if &#40;.false.&#41; then
              call OpenTxtFile&#40;'CAMBdefault_lensedCls.dat',3&#41;
              do i=lmin,lmax
               !Assume T,E,B,X ordering
               read&#40;3,*&#41; j, cl&#40;1&#58;4,i&#41;
               if &#40;j<lmin&#41; read&#40;3,*&#41; j, cl&#40;1&#58;4,i&#41;
               cl&#40;&#58;,i&#41;=cl&#40;&#58;,i&#41;/&#40;COBE_CMBTemp*1e6&#41;**2 
              end do
              close&#40;3&#41;
          end if

          if &#40;Noise&#40;1&#41; >0&#41; then
           file_tag = concat&#40;file_tag,'_noise'&#41;
          end if
          xlc= 180*sqrt&#40;8.*log&#40;2.&#41;&#41;/3.14159
          sigma2&#40;1&#58;2&#41; = &#40;BispectrumParams%FisherNoiseFwhmArcmin&#40;1&#58;2&#41;/60/xlc&#41;**2
          write&#40;*,*&#41; "sigma2&#58;", sigma2
          allocate&#40;InvC&#40;lmax&#41;&#41;
          do l1= lmin, lmax
             tmp = l1*&#40;l1+1&#41;/&#40;2*pi&#41;
             !PDM 2016&#58;
             !allow to combine 2 different experiments &#40;e.g. Planck and CMB-S4&#41;
             if&#40;l1 .le. BispectrumParams%EllTrans&#40;1&#41;&#41; then
                Cl&#40;1,l1&#41; = Cl&#40;1,l1&#41;/tmp + Noise&#40;1&#41;*exp&#40;l1*&#40;l1+1&#41;*sigma2&#40;1&#41;&#41;
             else
                !PDM 2016&#58; factor 2.5 is to make Noise = 1 muK for fsky = 0.4
                Cl&#40;1,l1&#41; = Cl&#40;1,l1&#41;/tmp + 2.5*Noise&#40;2&#41;*exp&#40;l1*&#40;l1+1&#41;*sigma2&#40;2&#41;&#41;
             endif

             if&#40;l1 .le. BispectrumParams%EllTrans&#40;1&#41;&#41; then
                Cl&#40;2&#58;3,l1&#41; = Cl&#40;2&#58;3,l1&#41;/tmp + NoiseP&#40;1&#41;*exp&#40;l1*&#40;l1+1&#41;*sigma2&#40;1&#41;&#41;
	     else
                Cl&#40;2&#58;3,l1&#41; = Cl&#40;2&#58;3,l1&#41;/tmp + 2.5*NoiseP&#40;2&#41;*exp&#40;l1*&#40;l1+1&#41;*sigma2&#40;2&#41;&#41;
             endif
             Cl&#40;4,l1&#41; = Cl&#40;4,l1&#41;/tmp  
             allocate&#40;InvC&#40;l1&#41;%C&#40;nfields,nfields&#41;&#41;
             if &#40;nfields > 2&#41; stop 'Not implemented nfields>2 in detail'
             if &#40;nfields==1&#41; then
              InvC&#40;l1&#41;%C&#40;1,1&#41;=&#40;2*l1+1&#41;/cl&#40;1,l1&#41;/InternalScale
             else 
              InvC&#40;l1&#41;%C&#40;1,1&#41;=cl&#40;2,l1&#41;
              InvC&#40;l1&#41;%C&#40;1,2&#41;=-cl&#40;4,l1&#41;
              InvC&#40;l1&#41;%C&#40;2,1&#41;=-cl&#40;4,l1&#41;
              InvC&#40;l1&#41;%C&#40;2,2&#41;=cl&#40;1,l1&#41;              
              InvC&#40;l1&#41;%C= InvC&#40;l1&#41;%C * &#40;2*l1+1&#41;/&#40;cl&#40;1,l1&#41;*cl&#40;2,l1&#41;-cl&#40;4,l1&#41;**2&#41;/InternalScale
             end if
          end do

          print *,'Getting Fisher for lmax = ', lmax   

          if &#40;debugMsgs&#41; starttime=GetTestTime&#40;&#41;
          
          !mysum = 0       
          ifish_contribs=0
          tmp = 0
          tmp1 = 0
          tmp2 = 0
          tmpf = 0
          !lmax = SampleL%l&#40;59&#41; 
         !&#36;OMP PARAllEl DO DEFAUlT&#40;SHARED&#41;,SCHEDULE&#40;STATIC,3&#41; &
         !&#36;OMP PRIVATE&#40;il1,l1,l2,l3,fish_l1,bi_ix,min_l,max_l,a3j_00,a3j&#41;, &
         !&#36;OMP PRIVATE&#40;Bispectrum,Bispectrum2,minl2,bix,tmp,tmp1,tmp2,tmpf&#41;, &
         !&#36;OMP PRIVATE&#40;field1,field2,field3,f1,f2,f3,bispectrum_type,bispectrum_type2&#41; 
         !!&#36;OMP REDUCTION&#40;+&#58;fish_l1&#41;
          do il1= ellminl, ellmaxl
          !do il1= 58,  59 
            !write&#40;*,*&#41; lmaxcuti 
            allocate&#40;fish_l1&#40;nbispectra,nbispectra,nfields,nfields&#41;&#41; !last indices are field1,f1
            l1 = SampleL%l&#40;il1&#41;
            if &#40;l1< lstart&#41; cycle
            fish_l1=0
            bi_ix=0
            do l2 = l1,lmax
              if &#40;l2< lstart&#41; cycle
              min_l = max&#40;lstart,max&#40;abs&#40;l1-l2&#41;,l2&#41;&#41;
              if &#40;mod&#40;l1+l2+min_l,2&#41;/=0&#41; then
                 min_l = min_l+1
              end if 
              max_l = min&#40;lmax,l1+l2&#41; 
              call GetThreeJs&#40;a3j&#40;abs&#40;l2-l1&#41;&#41;,l1,l2,0,0&#41;
              do l3=min_l,max_l ,2    
                a3j_00&#40;l3&#41;=a3j&#40;l3&#41;**2
              end do
              
              tmp1= 1.d0/&#40;4*pi&#41;  !&#40;2l+1&#41; factors included in InvC
              minl2=min_l
              bix=bi_ix
               do field1=1,nfields
                do f1=1,nfields
                 tmpf&#40;1&#41;= InvC&#40;l1&#41;%C&#40;field1,f1&#41;*tmp1
                 do field2=1,nfields
                   do f2=1,nfields
                    tmpf&#40;2&#41;= InvC&#40;l2&#41;%C&#40;field2,f2&#41;*tmpf&#40;1&#41;
                     do field3=1,nfields
                      do bispectrum_type=1,nbispectra
                      if &#40;bispectrum_type==lens_bispectrum_ix&#41; then
                        Bispectrum=>SqueezedLensingKernel&#40;field2,field3&#41;
                      else
                        Bispectrum=>bispectra&#40;field1,field2,field3,bispectrum_type&#41;
                      end if
                      do f3=1,nfields
                        
                        do bispectrum_type2=bispectrum_type,nbispectra
                        if &#40;bispectrum_type2==lens_bispectrum_ix&#41; then
                          Bispectrum2=>SqueezedLensingKernel&#40;f2,f3&#41;
                        else
                          Bispectrum2=>Bispectra&#40;f1,f2,f3,bispectrum_type2&#41;
                        end if
   
          
                        min_l=minl2          
                        bi_ix=bix
                        if &#40;min_l==l2&#41; then
                            !Symmetry factors
                             bi_ix=bi_ix+1
                             l3=l2
                             if &#40;l2==l1&#41; then
                              !l1=l2=l3
                              tmp = Bispectrum%b&#40;bi_ix,il1&#41;*tmpf&#40;2&#41;*Bispectrum2%b&#40;bi_ix,il1&#41; &
                                      *InvC&#40;l3&#41;%C&#40;field3,f3&#41;*a3j_00&#40;l3&#41;/6  
                             else
                              !l3=l2 &#40;l3=l1<>l2 can't happen because l1<=l2<=l3&#41;
                              tmp = Bispectrum%b&#40;bi_ix,il1&#41;*tmpf&#40;2&#41;*Bispectrum2%b&#40;bi_ix,il1&#41; &
                                      * InvC&#40;l3&#41;%C&#40;field3,f3&#41;*a3j_00&#40;l3&#41;/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&#40;bi_ix,il1&#41;*Bispectrum2%b&#40;bi_ix,il1&#41; &
                                   * InvC&#40;l3&#41;%C&#40;field3,f3&#41;*a3j_00&#40;l3&#41;
                        end do  
                        if &#40;l2==l1&#41; then
                           tmp2=tmp2*tmpf&#40;2&#41;/2
                          else
                           tmp2=tmp2*tmpf&#40;2&#41;     
                        end if     
                        fish_l1&#40;bispectrum_type,bispectrum_type2,field1,f1&#41;= &
                         fish_l1&#40;bispectrum_type,bispectrum_type2,field1,f1&#41;+&#40;tmp+tmp2&#41; 
                        
                        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&#40;il1,&#58;,&#58;,&#58;,&#58;&#41;=fish_l1&#40;1&#58;nbispectra,1&#58;nbispectra,&#58;,&#58;&#41; /InternalScale
            
            deallocate&#40;fish_L1&#41; 
         
          end do
         !&#36;OMP END PARAllEl DO
Thanks!
Daan

Re: CAMB sources

Posted: July 07 2016
by Antony Lewis
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.