CosmoCoffee Forum Index CosmoCoffee

 
 FAQFAQ   SearchSearch  MemberlistSmartFeed   MemberlistMemberlist    RegisterRegister 
   ProfileProfile   Log inLog in 
Arxiv New Filter | Bookmarks & clubs | Arxiv ref/author:

CAMB sources
 
Post new topic   Reply to topic    CosmoCoffee Forum Index -> Computers and software
View previous topic :: View next topic  
Author Message
daan meerburg



Joined: 20 Nov 2007
Posts: 21
Affiliation: CITA

PostPosted: July 01 2016  Reply with quote

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 CTT CEE CTE CPhi CPhiT Cwin_1 Cwin_2 ... CWin_T1 CWin_T2 .... Cwin1_win2 ...
Here all Cl are l(l+1)Cl/2pi except for CPhi and CPhiT which are CPhi= l4 Cl_phi where Cl_phi is the (CMB) lensing potential power spectrum and CPhiT = l3 Cl_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]
Back to top
View user's profile  
Antony Lewis



Joined: 23 Sep 2004
Posts: 1276
Affiliation: University of Sussex

PostPosted: July 01 2016  Reply with quote

Can you use the newer scalCovCls output file?
Back to top
View user's profile [ Hidden ] Visit poster's website
daan meerburg



Joined: 20 Nov 2007
Posts: 21
Affiliation: CITA

PostPosted: July 05 2016  Reply with quote

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 max and I have introduced (a partly exitsing) loop to run over max. However, this has a very persistent bug. When computing the 1 sigma error as a function of max, for very large max, the error bars increase! (so for low 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 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 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:
           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:

        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


Thanks!
Daan
Back to top
View user's profile  
Antony Lewis



Joined: 23 Sep 2004
Posts: 1276
Affiliation: University of Sussex

PostPosted: July 07 2016  Reply with quote

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.
Back to top
View user's profile [ Hidden ] Visit poster's website
Display posts from previous:   
Post new topic   Reply to topic    CosmoCoffee Forum Index -> Computers and software All times are GMT + 5 Hours
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © 2001, 2005 phpBB Group. Sponsored by WordWeb online dictionary and dictionary software.