!Module storing observed datasets, their points and window functions !and routines for computing the likelihood !Cls passed to these routines are in in MicroK^2, no 2pi or l(l+1) !WMAP data is handled separately as a special case !If windows_are_bandpowers=T Windows functions read in from files are assumed to be W_l/l. !(unless windows_are_bare=T, in which case we assume they are W_l) !The band powers are obtained from !W_l by \sum_l (l+1/2) W_l C_l/2pi where W_l is normalized so \sum_l (l+1/2)/l(l+1) W_l = 1. !If windows_are_normalized=T the windows are assumed to be already normalized; !this must be the case for polarization in which the window files contain ! l TT TE EE BB !contributions to each band power. Usually all but two columns will be zeros. !If windows_are_bandpowers=F then window functions are assumed raw, and directly related !to the bandpowers by B_i=sum_l W_{il} C_l !This code is very indebted to Sarah's cosmarg !Analytic marginalization over calibration and beam from astro-ph/0112114 !x factors are described in astro-ph/9808264. Essentially they are a first correction !for the bandpower errors being non-Gaussian. If no x-factors we assume they are gaussian. !Using x-factors we transform to a variable Z = log (bandpower + x), and assume Z has Gaussian errors. !When x-factors are used the inverse covariance matrix is assumed to be that for Z. !See the RADPACK page for data http://bubba.ucdavis.edu/~knox/radpack.html !This version August 2006 !Mar 04: added first_band parameter to .dataset files, added format for doing exact likelihoods !Jul 05: Readdataset_bcp changes for BOOM/CBI data, allowing for band cuts !Mar 06: changed to WMAP 3-year likelihood !Aug 06: added cal**2 scaling of x-factors !Oct 06: edited ReadAllExact to auto-account for number of cls (e.g. missing BB) !Oct 07: added Planck-like CMBLikes format !March 08: switched to support WMAP5 !May 09: modified ReadDataset_bcp for QUaD, allowing beam errors to be explicitly provided for each band module cmbdata use settings use cmbtypes use MatrixUtils use CMBLikes implicit none !if CMBdataset%has_xfactors then obs, var and N_inv are for the offset-lognormal variables Z Type CMBdatapoint integer win_min, win_max !Ranges in which window is non-zero real, pointer, dimension(:,:) :: window !Normalized window function in l real obs, err_minus, err_plus, sigma, var !Observed value of l(l+1) C_l/2pi bandpowers in MicroK^2, with errors real beam_err !fractional beam error (file is value in MicroK^2) logical inc_pol end Type CMBdatapoint Type CMBdataset logical :: use_set logical :: has_pol, windows_are_bandpowers,windows_are_normalized logical :: has_sz_template real :: calib_uncertainty logical :: beam_uncertain, has_corr_errors, has_xfactors integer :: num_points, file_points character(LEN=80) :: name Type(CMBdatapoint), pointer, dimension(:) :: points real, pointer, dimension(:,:) :: N_inv real, pointer, dimension(:) :: xfactors logical, pointer, dimension(:) :: has_xfactor !whether each bin has one logical :: all_l_exact logical :: CMBLike !New format integer :: all_l_lmax integer :: nuisance_parameters real, pointer, dimension(:,:) :: all_l_obs, all_l_noise real, pointer, dimension(:) :: all_l_fsky real, pointer, dimension(:) :: sz_template Type (TCMBLikes), pointer :: CMBLikes end Type CMBdataset integer :: num_datasets = 0 Type(CMBdataset) datasets(10) logical :: init_MAP = .true. integer :: cl_bin_width =1 integer, parameter :: halfsteps = 5 !do 2*halfsteps+1 steps in numerical marginalization real margeweights(-halfsteps:halfsteps) real :: margenorm = 0 private halfsteps, margeweights, margenorm contains subroutine ReadWindow(AP, aname, are_bare, aset) Type(CMBdatapoint) AP character(LEN=*), intent(IN) :: aname logical, intent(IN) :: are_bare Type (CMBdataset) :: aset integer l, ncls real wpol(1:num_cls-1),ll, IW,w character(LEN=200) tmp if (Feedback > 1) write (*,*) 'reading window: '//trim(aname) if (aset%has_pol) then ncls = num_cls else ncls = 1 endif allocate(AP%window(ncls,lmax)) AP%window = 0 call OpenTxtFile(aname, tmp_file_unit) do if (aset%has_pol) then read(tmp_file_unit,'(a)',end=1) tmp read(tmp,*, end=1) ll, w, wpol else read(tmp_file_unit,*, end=1) ll, w end if l=nint(ll) if (abs(l-ll) > 1e-4) stop 'non-integer l in window file' if (l>=2 .and. l<=lmax) then AP%window(1,l) = w if(aset%has_pol) then AP%window(2:num_cls,l) = wpol end if if (.not. are_bare) AP%window(:,l) = AP%window(:,l)*l else if (l>lmax .and. w /= 0) then write (*,*) 'Warning: Window function non-zero outside l_max: ',trim(aname) write (*,*) 'assuming contribution is negligible' exit end if end if end do 1 close(tmp_file_unit) do l=2, lmax if (any(AP%window(:,l)/=0)) then AP%win_min = l exit end if end do do l=lmax,2,-1 if (any(AP%window(:,l)/=0)) then AP%win_max = l exit end if end do AP%inc_pol = .false. if (aset%windows_are_bandpowers) then IW = 0 do l = AP%win_min, AP%win_max IW = IW + ((l+0.5)*AP%window(1,l))/(l*(l+1)) AP%window(:,l) = AP%window(:,l)*(l+0.5)/(2*pi) end do if (.not. aset%windows_are_normalized) then if (aset%has_pol) stop 'Polarization windows must be normalized' AP%window(1,AP%win_min:AP%win_max) = & AP%window(1,AP%win_min:AP%win_max)/IW end if end if !If has_pol we are assuming windows are normalized if (aset%has_pol) AP%inc_pol = any(AP%window(2:num_cls,:)/=0) end subroutine ReadWindow subroutine ReadAllExact(Ini,aset) Type(TIniFile) :: Ini Type (CMBdataset) :: aset character(LEN=Ini_max_string_len) :: fname integer l, idum, ncls, ncol real inobs(4) 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) 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' if (aset%has_pol) then allocate(aset%all_l_obs(2:aset%all_l_lmax,num_cls)) allocate(aset%all_l_noise(2:aset%all_l_lmax,2)) else allocate(aset%all_l_obs(2:aset%all_l_lmax,1)) allocate(aset%all_l_noise(2:aset%all_l_lmax,1)) 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) if (ncol==7) then ncls = 3 elseif (ncol==8) then ncls = 4 elseif (ncol==4) then ncls=1 if (aset%has_pol) stop 'cmbdata.f90::ReadAllExact: has_pol wrong' else stop 'cmbdata.f90::ReadAllExact: wrong number of columns' end if file_unit = new_file_unit() call OpenTxtFile(fname,file_unit) !File format: !l C_TT (C_TE C_EE [C_BB]) N_T (N_P) 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,inobs(1:ncls), aset%all_l_noise(l,:), aset%all_l_fsky(l) if (idum /= l) stop 'Error reading all_l_file' !set BB to pure noise if not in file if (aset%has_pol .and. ncls < num_cls) inobs(num_cls) = aset%all_l_noise(l,2) aset%all_l_obs(l,:) = inobs(1:num_cls) 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) Type(CMBdataset) :: aset integer l real ChiSqExact, chisq, term, CT, CE, CB, CC real CThat, CChat, CEhat real dof integer i chisq=0 do l=2, 30 dof = aset%all_l_fsky(l)**2*(2*l+1) !Ignoring l correlations but using f_sky^2_eff fudge factor may be a good approx !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) term = CT*CE- cl(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))/term & + log( term/ (aset%all_l_obs(l,1)*aset%all_l_obs(l,3) - aset%all_l_obs(l,2)**2)) -2) if (num_cls>3) then !add in BB CB = cl(l,num_cls) + aset%all_l_noise(l,2) chisq = chisq + dof * (aset%all_l_obs(l,num_cls)/CB & +log(CB/aset%all_l_obs(l,num_cls)) - 1) end if else chisq = chisq + dof * (aset%all_l_obs(l,1)/CT & +log(CT/aset%all_l_obs(l,1)) - 1) end if end do do l=31 , aset%all_l_lmax, cl_bin_width dof = 0 CT=0 CE=0 CC=0 CThat=0 CEhat=0 CChat=0 do i=l,l+cl_bin_width-1 dof = dof + aset%all_l_fsky(i)**2*(2*i+1) CT = CT + (cl(l,1) + aset%all_l_noise(l,1))*(2*i+1) CThat = CThat + aset%all_l_obs(l,1)*(2*i+1) if (aset%has_pol) then CE = CE + (cl(l,3) + aset%all_l_noise(l,2) ) *(2*i+1) CC = CC + cl(l,2)*(2*i+1) CEhat = CEhat + aset%all_l_obs(l,3)*(2*i+1) CChat = CChat + aset%all_l_obs(l,2)*(2*i+1) end if end do if (aset%has_pol) then term = CT*CE - CC**2 chisq = chisq + dof*( & (CT*CEHat + CE*CThat - 2 *CC*CCHat)/term & + log( term/ (CTHat*CEHat - CCHat**2)) -2) if (num_cls>3) then !add in BB CB = cl(l,num_cls) + aset%all_l_noise(l,2) chisq = chisq + dof * (aset%all_l_obs(l,num_cls)/CB & +log(CB/aset%all_l_obs(l,num_cls)) - 1) end if else chisq = chisq + dof * (CTHat/CT +log(CT/CTHat) - 1) end if end do ChiSqExact = ChiSq end function ChiSqExact subroutine ReadDataset(aname) use CMBLikes character(LEN=*), intent(IN) :: aname character(LEN=Ini_max_string_len) :: InLine, window_dir, Ninv_file, xfact_file, band_file logical bad, windows_are_bare Type (CMBdataset) :: aset integer i, first_band, use_i real, pointer, dimension(:,:) :: tmp_mat real, pointer, dimension(:) :: tmp_arr character(LEN=Ini_max_string_len) :: data_format integer file_unit Type(TIniFile) :: Ini num_datasets = num_datasets + 1 if (num_datasets > 10) stop 'too many datasets' aset%has_sz_template = .false. aset%nuisance_parameters = 0 !Special cases if (aname == 'MAP' .or. aname == 'WMAP') then aset%name = 'WMAP' datasets(num_datasets) = aset return elseif( aname(LEN_TRIM(aname)-5:LEN_TRIM(aname)) == 'newdat') then !Carlo format for polarized Boomerang et al. if (Feedback > 0) write(*,*) 'Reading BCP data set: ' // TRIM(aname) call ReadDataset_bcp(aname) return end if file_unit = new_file_unit() call Ini_Open_File(Ini,aname, file_unit, bad, .false.) if (bad) then write (*,*) 'Error opening dataset file '//trim(aname) stop end if Ini_fail_on_not_found = .true. aset%name = Ini_Read_String_File(Ini,'name') aset%use_set =.true. aset%num_points = 0 if (Feedback > 0) write (*,*) 'reading: '//trim(aset%name) Ini_fail_on_not_found = .false. data_format = Ini_Read_String_File(Ini,'dataset_format') aset%CMBlike = data_format == 'CMBLike' aset%all_l_exact = (data_format =='all_l_exact') & .or. Ini_Read_Logical_File(Ini,'all_l_exact',.false.) if (aset%CMBLike) then allocate(aset%CMBLikes) call CMBLikes_ReadData(aset%CMBLikes, Ini, ExtractFilePath(aname)) aset%nuisance_parameters = aset%CMBLikes%num_nuisance_parameters else if (aset%all_l_exact) then aset%has_pol = Ini_Read_Logical_File(Ini,'has_pol',.false.) call ReadAllExact(Ini,aset) else if (data_format/='') then write(*,*) 'Error in '//trim(aname) write(*,*) 'Unknown data_format: '//trim(data_format) stop else !Otherwise do usual guassian/offset lognormal stuff aset%has_pol = Ini_Read_Logical_File(Ini,'has_pol',.false.) aset%num_points = Ini_Read_Int_File(Ini,'num_points') aset%calib_uncertainty = Ini_Read_Double_File(Ini,'calib_uncertainty') aset%beam_uncertain = Ini_Read_logical_File(Ini,'beam_uncertainty') window_dir = Ini_Read_String_File(Ini,'window_dir') windows_are_bare = Ini_Read_Logical_File(Ini,'windows_are_bare',.false.) aset%windows_are_bandpowers = Ini_Read_Logical_File(Ini,'windows_are_bandpowers',.true.) aset%windows_are_normalized = Ini_Read_Logical_File(Ini,'windows_are_normalized',.false.) aset%file_points = Ini_read_Int_File(Ini,'file_points',aset%num_points) first_band = Ini_read_Int_File(Ini,'first_band',1) if (first_band + aset%num_points > aset%file_points+1) then write (*,*) 'Error with dataset file '//trim(aname) write (*,*) 'first_band + num_points > file_points' stop end if !Read in the observed values, errors and beam uncertainties allocate(aset%points(aset%num_points)) band_file = Ini_Read_String_File(Ini,'bandpowers') if (band_file /= '') call OpenTxtFile(band_file, 51) Ini_fail_on_not_found = .true. do i=1, aset%num_points + first_band -1 if (band_file /= '') then read(51,'(a)') InLine else InLine = Ini_Read_String_File(Ini,numcat('data',i)) end if if (i < first_band) cycle use_i = i - first_band + 1 if (aset%beam_uncertain) then read(InLine,*) aset%points(use_i)%obs, aset%points(use_i)%err_minus, aset%points(use_i)%err_plus, & aset%points(use_i)%beam_err aset%points(use_i)%beam_err = aset%points(use_i)%beam_err/aset%points(use_i)%obs else read(InLine,*) aset%points(use_i)%obs, aset%points(use_i)%err_minus, aset%points(use_i)%err_plus aset%points(use_i)%beam_err = 0 end if aset%points(use_i)%sigma = (aset%points(use_i)%err_minus + aset%points(use_i)%err_plus)/2 aset%points(use_i)%var = aset%points(use_i)%sigma**2 call ReadWindow(aset%points(use_i),trim(window_dir)//'/'//trim(numcat(aset%name,i)),windows_are_bare,aset) end do if (band_file /= '') Close(51) !See if the inverse covariance matrix is given (otherwise assume diagonal) Ini_fail_on_not_found = .false. Ninv_file = Ini_Read_String_File(Ini,'N_inv') aset%has_corr_errors = Ninv_file /= '' if (aset%has_corr_errors) then allocate(tmp_mat(aset%file_points,aset%file_points)) allocate(aset%N_inv(aset%num_points,aset%num_points)) call ReadMatrix(Ninv_file, tmp_mat,aset%file_points,aset%file_points) if (aset%num_points /= aset%file_points) then !!Changed to truncation of covariance matrix, AL: May 2003 call Matrix_Inverse(tmp_mat) aset%N_inv = tmp_mat(first_band:first_band+aset%num_points-1,& first_band:first_band + aset%num_points -1) call Matrix_Inverse(aset%N_inv) else aset%N_inv = tmp_mat(1:aset%num_points,1:aset%num_points) end if deallocate(tmp_mat) end if if (Ini_Read_Logical_File(Ini,'use_hyperparameter',.false.)) stop 'Hyperparameters deprecated' !See if xfactors are given xfact_file = Ini_Read_String_File(Ini,'xfactors') aset%has_xfactors = xfact_file /= '' if (aset%has_xfactors) then allocate(tmp_arr(aset%num_points + first_band -1)) call ReadVector(xfact_file, tmp_arr,aset%num_points+ first_band -1) allocate(aset%xfactors(aset%num_points)) allocate(aset%has_xfactor(aset%num_points)) aset%has_xfactor = .true. aset%xfactors = tmp_arr(first_band:first_band+aset%num_points-1) deallocate(tmp_arr) aset%points(:)%var = aset%points(:)%var/(aset%points(:)%obs +aset%xfactors)**2 aset%points(:)%obs = log(aset%points(:)%obs +aset%xfactors) end if end if !not all_l_exact or cut sky unbinned call Ini_Close_File(Ini) call ClearFileUnit(file_unit) datasets(num_datasets) = aset end subroutine ReadDataset subroutine ReadSZTemplate(aset, aname, ascale) Type (CMBdataset) :: aset real, intent(in) :: ascale character(LEN=*), intent(IN) :: aname integer l, unit real sz allocate(aset%sz_template(2:lmax)) aset%sz_template = 0 aset%has_sz_template = .true. call OpenTxtFile(aname, tmp_file_unit) do read(tmp_file_unit,*,end=2) l, sz if (l>=2 .and. l<=lmax) aset%sz_template(l) = ascale * sz/(l*(l+1)/twopi) end do 2 Close(tmp_file_unit) end subroutine ReadSZTemplate function GetWinBandPower(AP, cl) real GetWinBandPower real cl(lmax,num_cls) Type(CMBdatapoint) AP integer l real bandpower bandpower = sum(cl(AP%win_min:AP%win_max,1)*AP%window(1,AP%win_min:AP%win_max)) if (AP%inc_pol) then do l= AP%win_min, AP%win_max bandpower = bandpower + sum(cl(l,2:num_cls)*AP%window(2:num_cls,l)) end do endif GetWinBandPower = bandpower end function GetWinBandPower subroutine InitNumericalMarge integer i do i= -halfsteps, halfsteps margeweights(i) = exp(-(i*3/real(halfsteps))**2/2) end do margenorm = sum(margeweights) end subroutine InitNumericalMarge function GetCalibMargexChisq(bandpowers, aset) !Numerically integrate over the calibration uncertainty !Assume Gaussian prior, as for analytic calculation without x-factors !Could also Monte-Carlo real GetCalibMargexChisq Type(CMBdataset) aset real bandpowers(aset%num_points),beambandpowers(aset%num_points),diffs(aset%num_points) real calib, chisq(-halfsteps:halfsteps),chisqcalib(-halfsteps:halfsteps) real minchisq integer i,j, ibeam if (margenorm == 0) call InitNumericalMarge do ibeam= -halfsteps, halfsteps if (aset%beam_uncertain) then beambandpowers = bandpowers*(1 + aset%points(:)%beam_err*ibeam*3/real(halfsteps))!Go out to 3 sigma else beambandpowers = bandpowers end if do i=-halfsteps,halfsteps calib = 1 + aset%calib_uncertainty*i*3./halfsteps !Go out to 3 sigma if (aset%has_xfactors) then do j=1, aset%num_points if (aset%has_xfactor(j)) then diffs(j) = aset%points(j)%obs- log(calib*beambandpowers(j) + aset%xfactors(j)) else diffs(j) = aset%points(j)%obs - calib*beambandpowers(j) endif end do else diffs = aset%points(:)%obs - calib*beambandpowers end if if (aset%has_corr_errors) then chisq(i) = SUM(diffs*MATMUL(aset%N_inv,diffs)) else chisq(i) = SUM(diffs**2/aset%points(:)%var) end if end do minchisq = minval(chisq) chisqcalib(ibeam) = -2*log(sum(margeweights*exp(max(-30.,-(chisq-minchisq)/2)))/margenorm) + minchisq if (.not. aset%beam_uncertain) then GetCalibMargexChisq = chisqcalib(ibeam) return end if end do minchisq = minval(chisqcalib) GetCalibMargexChisq = -2*log(sum(margeweights*exp(max(-30.,-(chisqcalib-minchisq)/2)))/margenorm) + minchisq end function GetCalibMargexChisq !Routine by Carlo Contaldi to read .newdat file format (Boomerang et al) !Modified to account for offset lognormal toggle per band !AL July 2005: modified to allow band selection !MLB May 09: modified to allow provision of per-band beam errors SUBROUTINE ReadDataset_bcp(aname) !File Format: !name !n_bands_TT n_EE, n_BB, n_EB, n_TE, n_TB !has_calib_uncertain calib(amplitude) calib_err(power) !has_beam_uncertain beam beam_err !beam control flag: if has_beam_uncertain = 0 => no beam errors, ! 1 => use FWHM and FWHM error as provided, 2 => read per-band errors from tables !ilike (0: Gaussian, 1: all x-factor, 2: specified have x-factor) !loop over { ! band-types ! band info: num obs + - x_factor l_min l_max use_x ! correlation-matrix (ignored) ! } ! covariance matrix CHARACTER(LEN=*), INTENT(IN) :: aname CHARACTER(LEN=100) :: instr CHARACTER(LEN=3), DIMENSION(1:6) :: ch_types LOGICAL windows_are_bare TYPE (CMBdataset) :: aset INTEGER i, j, k,use_i, n_types, xin REAL, POINTER, DIMENSION(:) :: tmp_x REAL, POINTER, DIMENSION(:,:) :: lb real, allocatable, dimension(:,:) :: tmp_mat integer, allocatable, dimension(:) :: used_bands INTEGER :: npol(6), minmax(2,6) INTEGER :: file_i,ijunk, ilike, ibeam REAL :: cal, beam_width, beam_sigma, l_mid integer, parameter :: like_xfactorall=1, like_xfactorsome = 2 !to be compatible with some older CITA output files LOGICAL :: FISHER_T_CMB CALL OpenTxtFile(aname, 62) READ(62,'(a)') instr FISHER_T_CMB = .FALSE. IF(instr == 'FISHER_T_CMB') THEN FISHER_T_CMB = .TRUE. READ(62,'(a)') instr WRITE(*,'(a)') 'FISHER_T_CMB is set for :'//TRIM(ADJUSTL(instr)) ENDIF aset%name = TRIM(ADJUSTL(instr)) WRITE(*,*) 'Reading: '//TRIM(ADJUSTL(aset%name)) aset%use_set =.TRUE. READ(62,*) npol(1:6) aset%has_pol = any(npol(2:6) /=0) aset%all_l_exact = .FALSE. aset%file_points = SUM(npol) aset%num_points = SUM(npol) n_types = count(npol /= 0) READ(62,'(a)') instr IF(instr == 'BAND_SELECTION') THEN !list of 'first_band last_band' for each pol type !if first_band=0 then ignore that pol type aset%num_points = 0 aset%has_pol = .false. if(feedback>0) WRITE(*,*) 'Using selected band ranges' do i=1,6 READ(62,*) minmax(1:2,i) if (minmax(1,i)/=0) then aset%num_points = aset%num_points + minmax(2,i) - minmax(1,i) + 1 if (i>1) aset%has_pol = .true. else minmax(2,i) = 0 end if end do READ(62,'(a)') instr ELSE !use all bands in file do i=1,6 minmax(1,i)=1 minmax(2,i)=npol(i) end do ENDIF READ(instr,*) ijunk, cal, aset%calib_uncertainty IF(ijunk == 0) aset%calib_uncertainty = 0.e0 READ(62,*) ibeam, beam_width, beam_sigma aset%beam_uncertain = ibeam /= 0 !this agrees with latest windows coming out of MPIlikely windows_are_bare = .FALSE. aset%windows_are_bandpowers = .TRUE. aset%windows_are_normalized = .TRUE. ALLOCATE(aset%points(aset%num_points)) allocate(used_bands(aset%num_points)) READ(62,*) ilike aset%has_xfactors = ilike ==like_xfactorsome .or. ilike==like_xfactorall !1 : all bands are offset lognormal !2 : only bands specified have offset lognormal IF(aset%has_xfactors) then ALLOCATE(aset%has_xfactor(1:aset%num_points)) aset%has_xfactor = .true. end if ALLOCATE(tmp_x(1:aset%num_points)) ALLOCATE(lb(1:aset%num_points,1:2)) k = 0 use_i = 0 file_i = 0 DO j=1,n_types READ(62,'(a2)') ch_types(j) k = k + 1 DO i=1,20 IF(npol(k) == 0 ) THEN k = k + 1 ELSE EXIT ENDIF ENDDO if(feedback>1) write(*,*) j, ch_types(j), minmax(1,k), minmax(2,k) DO i=1,npol(k) file_i = file_i+1 if (i>=minmax(1,k) .and. i<=minmax(2,k)) then use_i = use_i + 1 used_bands(use_i) = file_i IF(ibeam .le. 1) THEN IF(ilike /= like_xfactorsome) THEN READ(62,'(a)') instr READ(instr,*) ijunk, aset%points(use_i)%obs,aset%points(use_i)%err_minus, & aset%points(use_i)%err_plus,tmp_x(use_i),lb(use_i,1),lb(use_i,2) ELSE !like_xfactorsome !read also offset switch per band READ(62,*) ijunk, aset%points(use_i)%obs,aset%points(use_i)%err_minus, & aset%points(use_i)%err_plus,tmp_x(use_i),lb(use_i,1),lb(use_i,2),xin aset%has_xfactor(use_i) = xin/=0 ENDIF !beam error in bandpower l_mid = (lb(use_i,2)-lb(use_i,1))/2.d0 + lb(use_i,1) aset%points(use_i)%beam_err = exp(-l_mid*(l_mid+1.d0)*1.526e-8*2.d0*beam_sigma*beam_width)-1.d0 aset%points(use_i)%beam_err = abs(aset%points(use_i)%beam_err) ELSE IF(ilike /= like_xfactorsome) THEN READ(62,'(a)') instr READ(instr,*) ijunk, aset%points(use_i)%obs,aset%points(use_i)%err_minus, & aset%points(use_i)%err_plus,tmp_x(use_i),lb(use_i,1),lb(use_i,2),aset%points(use_i)%beam_err ELSE !like_xfactorsome !read also offset switch per band READ(62,*) ijunk, aset%points(use_i)%obs,aset%points(use_i)%err_minus, & aset%points(use_i)%err_plus,tmp_x(use_i),lb(use_i,1),lb(use_i,2),xin,aset%points(use_i)%beam_err aset%has_xfactor(use_i) = xin/=0 ENDIF l_mid = (lb(use_i,2)-lb(use_i,1))/2.d0 + lb(use_i,1) ENDIF aset%points(use_i)%sigma = (aset%points(use_i)%err_minus + aset%points(use_i)%err_plus)/2 if (Feedback>1 ) print*, aset%beam_uncertain, l_mid, aset%points(use_i)%beam_err !recalibrate aset%points(use_i)%obs = cal**2 * aset%points(use_i)%obs aset%points(use_i)%sigma = cal**2 * aset%points(use_i)%sigma aset%points(use_i)%var = aset%points(use_i)%sigma**2 CALL ReadWindow(aset%points(use_i),'data/windows/'//& TRIM(numcat(aset%name,file_i)),windows_are_bare,aset) else !discard band READ(62,'(a)') instr end if ENDDO !discard correlation submatrix READ(62,'(a)') (instr,i=1,npol(k)) ENDDO !assume always have the matrix aset%has_corr_errors = .TRUE. allocate(tmp_mat(aset%file_points,aset%file_points)) ALLOCATE(aset%N_inv(aset%num_points,aset%num_points)) READ(62,*) (tmp_mat(1:aset%file_points,i),i=1,aset%file_points) aset%N_inv = tmp_mat(used_bands,used_bands) deallocate(tmp_mat) deallocate(used_bands) !READ(62,*,err=101,end=101) instr !stop 'ReadDataset_bcp:Should be at end of file after reading matrix' 101 CLOSE(62) !recalibrate and change units as required !some older output had final fisher matrix in !in units of T_CMB whitle bandpowers are in units !of \microK^2 IF(FISHER_T_CMB) THEN aset%N_inv = cal**4 * aset%N_inv *2.725**4 * 1.e24 ELSE aset%N_inv = cal**4 * aset%N_inv ENDIF IF(aset%has_pol) WRITE(*,*) 'has pols: ', ADJUSTR(ch_types(1:n_types)) !transform into Z_B = ln(C_B+x_B) for bandpowers with xfactors IF(aset%has_xfactors) THEN ALLOCATE(aset%xfactors(1:aset%num_points)) aset%xfactors(1:aset%num_points) = cal**2*tmp_x(1:aset%num_points) DO i=1,aset%num_points DO j=1,aset%num_points if(aset%has_xfactor(i)) aset%N_inv(i,j) = aset%N_inv(i,j)/(aset%points(i)%obs + aset%xfactors(i)) if(aset%has_xfactor(j)) aset%N_inv(i,j) = aset%N_inv(i,j)/(aset%points(j)%obs + aset%xfactors(j)) ENDDO ENDDO DO i=1,aset%num_points IF(aset%has_xfactor(i)) THEN aset%points(i)%var = aset%points(i)%var/(aset%points(i)%obs +aset%xfactors(i))**2 aset%points(i)%obs = LOG(aset%points(i)%obs + aset%xfactors(i)) ENDIF ENDDO CALL Matrix_Inverse(aset%N_inv) ELSE CALL Matrix_Inverse(aset%N_inv) ENDIF !add to full dataset datasets(num_datasets) = aset DEALLOCATE(tmp_x) END SUBROUTINE ReadDataset_bcp function CalcLnLike(cl, aset,nuisance_params) !Compute -ln(Likelihood) real cl(lmax,num_cls), CalcLnLike real, intent(in) :: nuisance_params(:) Type(CMBdataset) aset integer i real chisq real chi2op, chi2pp, wpp, wdd real chi2dd,chi2pd,chi2od real bandpowers(aset%num_points), diffs(aset%num_points), tmp(aset%num_points), beam(aset%num_points) real denom if (.not. aset%use_set) then CalcLnLike = 0 return end if if (aset%CMBLike) then chisq = CMBLikes_CMBLike(aset%CMBLikes, cl, nuisance_params) else if (aset%all_l_exact) then chisq = ChiSqExact(cl,aset) else denom = 1 !Assume Prob \propto exp(-chisq/2)/sqrt(denom) do i=1, aset%num_points bandpowers(i) = GetWinBandPower(aset%points(i), cl) end do if (aset%has_xfactors .and. (aset%calib_uncertainty > 1e-4 .or. aset%beam_uncertain)) then chisq = GetCalibMargexChisq(bandpowers,aset) else if (aset%has_xfactors) then do i=1, aset%num_points if (aset%has_xfactor(i)) then !obs in this case is Z = log(observed + x) diffs(i) = aset%points(i)%obs- log(bandpowers(i) + aset%xfactors(i)) else diffs(i) = aset%points(i)%obs - bandpowers(i) end if end do else diffs = aset%points(:)%obs - bandpowers end if if (aset%has_corr_errors) then chisq = SUM(diffs*MATMUL(aset%N_inv,diffs)) else chisq = SUM(diffs**2/aset%points(:)%var) end if if (aset%calib_uncertainty > 1e-4 .or. aset%beam_uncertain) then if (aset%has_corr_errors) then tmp = MATMUL(aset%N_inv,bandpowers) else tmp = bandpowers/aset%points(:)%var end if chi2op = SUM(diffs*tmp) chi2pp = SUM(bandpowers*tmp) if (aset%beam_uncertain) then beam = aset%points(:)%beam_err*bandpowers if (aset%has_corr_errors) then tmp = MATMUL(aset%N_inv,beam) else tmp = beam/aset%points(:)%var end if chi2dd = SUM(beam*tmp) chi2pd = SUM(bandpowers*tmp) chi2od = SUM(diffs*tmp) end if if (aset%calib_uncertainty > 1e-4) then !analytic marginalization over calibration uncertainty wpp = 1/(chi2pp+1/aset%calib_uncertainty**2) chisq = chisq - wpp*chi2op**2 denom = denom/wpp*aset%calib_uncertainty**2 else wpp = 0 end if if (aset%beam_uncertain) then !analytic marginalization over beam uncertainty wdd=1/(chi2dd-wpp*chi2pd**2+1) chisq = chisq - wdd*(chi2od-wpp*chi2op*chi2pd)**2 denom = denom/wdd end if end if end if if (denom /= 1) chisq = chisq + log(denom) end if CalcLnLike = chisq/2 end function CalcLnLike function CMBLnLike(cl, sznorm, nuisance_params) real, intent(in) :: cl(lmax,num_cls), sznorm real CMBLnLike real,intent(in) :: nuisance_params(:) real szcl(lmax,num_cls) integer i integer nuisance real tot(num_datasets) nuisance =1 do i=1, num_datasets szcl = cl if (datasets(i)%has_sz_template) then szcl(2:lmax,1) = szcl(2:lmax,1) + sznorm*datasets(i)%sz_template(2:lmax) end if if (datasets(i)%name == 'WMAP') then tot(i) = MAPLnLike(szcl) else tot(i) = CalcLnLike(szcl,datasets(i), nuisance_params(nuisance:)) if (datasets(i)%CMBLike) nuisance = nuisance + datasets(i)%CMBLikes%num_nuisance_parameters end if end do CMBLnLike = SUM(tot) end function function MAPLnLike(cl) #ifndef NOWMAP use wmap_likelihood_5yr use WMAP_OPTIONS use WMAP_UTIL #endif real cl(lmax,num_cls), MAPLnLike #ifndef NOWMAP real(8), dimension(2:ttmax) :: cl_tt,cl_te,cl_ee,cl_bb real(8) :: like(num_WMAP),like_tot integer l if (Init_MAP) then #ifdef WMAPNOHIGHLTT use_TT = .false. use_TT_beam_ptsrc = .false. #endif if (lmax0) write(*,*) 'reading WMAP5 data' Init_MAP = .false. end if do l = 2, ttmax cl_tt(l) = cl(l,1)*l*(l+1)/twopi cl_te(l) = cl(l,2)*l*(l+1)/twopi cl_ee(l) = cl(l,3)*l*(l+1)/twopi if(num_cls == 4) then cl_bb(l) = cl(l,num_cls)*l*(l+1)/twopi else cl_bb(l) = 0.0d0 end if end do like=0.0d0 call wmap_likelihood_compute(cl_tt,cl_te,cl_ee,cl_bb,like) !call wmap_likelihood_error_report if (wmap_likelihood_ok) then MAPLnLike = sum(like) else MAPLnLike = LogZero endif #else MAPLnLike=cl(2,1) !just stop unused symbol warnings stop 'Compiled without WMAP' #endif end function !WMAP1 ! function MAPLnLike(cl) ! use WMAP ! real cl(lmax,num_cls), MAPLnLike ! integer l ! real(WMAP_precision) clTT(WMAP_lmax_TT), clTE(WMAP_lmax_TT), ClEE(WMAP_lmax_TT) ! integer stat ! character(LEN=20) :: ttFile = 'WMAP/tt_diag.dat' ! character(LEN=20) :: ttOffDiag ='WMAP/tt_offdiag.dat' ! character(LEN=20) :: teFile = 'WMAP/te_diag.dat' ! character(LEN=20) :: teOffDiag ='WMAP/te_offdiag.dat' ! ! if (Init_MAP) then ! if (lmax0) write(*,*) 'reading WMAP data' ! Call WMAP_init(ttFile, ttOffDiag, teFile, teOffDiag, stat)! ! ! if (stat /=0) stop 'Error reading WMAP files' ! if (Feedback>0) write(*,*) 'WMAP read' ! Init_MAP = .false. ! end if ! do l = 2, WMAP_lmax_TT ! clTT(l) = cl(l,1)*l*(l+1)/twopi ! clTE(l) = cl(l,2)*l*(l+1)/twopi ! clEE(l) = cl(l,3)*l*(l+1)/twopi ! end do ! MAPLnLike = -(WMAP_LnLike_TT(clTT) + WMAP_LnLike_TE(clTT, clTE, clEE)) ! end function end module cmbdata