function post_200709b(ds,ss,s,t)
%
% Used to make bpwf_test part of posting 200709_bkfinal_06blt
%
% check CMB-S4 mean of sims against expvals
% See post_200528.m
%
% ds=dataset='06b'
% ss=simset='sigl','sim'
% s=spectra, 3=EE, 4=BB
% t=plottype, 1=meanofsims/expv, 3=ratio

switch ds
 case '06b'
  % get the BK15 data
  %load('final/4325/real_06bp00_AL1p0_jack0_directbpwf_broken_bpwf.mat');
  load('final/4325/real_06bp00_AL1p0_jack0_directbpwf.mat');
  A_L=1;
  sp=1:10; %length(opt.finalopt.mapname);
 case '04'
  load('final/1699/real_04p00_AL1p0_jack0_directbpwf.mat');
  A_L=1.0;
  sp=[1:9,NaN];
end

sn=opt.finalopt.mapname;

%n=max(find(cellfun(@isempty,strfind(opt.finalopt.mapname,'x'))))
n=length(sp);

% The provided expv are for unlensed LCDM - and these sims are not passed forward
% in the final file. Need to get the model and re-make
inpmodsl=load_cmbfast(opt.finalopt.inpmodsl);
inpmodsl.Cs_l(:,4)=inpmodsl.Cs_l(:,4)*A_L;
r=calc_expvals(r,inpmodsl,bpwf);

% don't plot the junk bin
r(1).l(1)=NaN;

% get expv via multicomp code
switch s
 case 3
  likeopt.fields = {'E'}; yl=[0,25];
 case 4
  likeopt.fields = {'B'}; yl=[0,0.2];
end
likeopt.l=1:17;
ld = get_likedata(ds,likeopt);
switch ss
 case 'sim'
  % from cmbs4/make_s4_gcomps.m
  mod = [0,A_L,3.8,4.25,-3.1,-0.6,1.6,-0.4,2,2,0,19.6]
 case 'sigl'
  % same without sync/dust
  mod = [0,A_L,0.0,0.00,-3.1,-0.6,1.6,-0.4,2,2,0,19.6]
end
mlh=like_getexpvals(mod,ld);

setwinsize(gcf,800,700); clf
% plot the auto/cross spectra
for i=1:n
  for j=i:n
    if(~any(isnan([sp(i),sp(j)])))
      k=aps_getxind(r,sp(i),sp(j))
      l=likeopt_getxind(ld.opt,sp(i),sp(j));
      
      [i j k l]
      subplot_grid(n,n,(j-1)*n+i,1,1);

      m=mlh(:,l)';
      switch ss
       case 'sim'
        % r=0 in odd numbered rlz and non-zero in even - see make_s4_comb_XXp00
        % that is for zero based indecies - so here r=0 in even numbered
        sims=squeeze(r(k).(ss)(:,s,2:2:end))';
       case 'sigl'
        % sigl has no r so may as well use all the rlz
        sims=squeeze(r(k).(ss)(:,s,:))';
      end
      x=mean(sims);
      sd=std(sims);
      e=sd/sqrt(size(sims,1));
      
      switch t
       case 1
        % mean of sims with errorbars plotted with expvals
        plot(r(1).l,r(k).expv(:,s),'r');
        hold on
        plot(r(1).l,m,'g');
        errorbar2(r(1).l,x,e,'.'); %ylim(yl);
        xlim([0,330]);  
       case 2
        % fractional deviations of above
        plot(r(1).l,(x-m)./e,'.-'); ylim([-3,3]);
        xlim([0,600]);
       case 3
        % mean of sims divided by expvals
        errorbar2(r(1).l,x./m,e./m,'.'); ylim([0.9,1.1]);
        grid
        xlim([0,600]);
       case 4
        % Ben posting style
        errorbar2(r(1).l,(x-m)./sd,ones(size(r(1).l))/sqrt(size(sims,1)),'.');
       case 5
        % ratio of expvals from old code versus likelihood code -
        % should be the same
        plot(r(1).l,r(k).expv(:,s)./m'); ylim([0.9,1.1]);

      end
      
      if(i==j)
        title(sn{i});
      end
      subplot_grid2(n,n,(i-1)*n+j,1,1);
    end
  end
end

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
post_200709b('06b','sim',4,1);

% build the plots for posting
post_200709b('06b','sigl',3,1); mkpng('temp/06b_sigl_ee_spec.png',0);
post_200709b('06b','sigl',3,3); mkpng('temp/06b_sigl_ee_ratio.png',0);
post_200709b('06b','sigl',4,1); mkpng('temp/06b_sigl_bb_spec.png',0);
post_200709b('06b','sigl',4,3); mkpng('temp/06b_sigl_bb_ratio.png',0);

post_200709b('06b','sim',3,1); mkpng('temp/06b_sim_ee_spec.png',0);
post_200709b('06b','sim',3,3); mkpng('temp/06b_sim_ee_ratio.png',0);
post_200709b('06b','sim',4,1); mkpng('temp/06b_sim_bb_spec.png',0);
post_200709b('06b','sim',4,3); mkpng('temp/06b_sim_bb_ratio.png',0);

% 04 does not contain .sigl field
post_200709b('04','sim',3,1); mkpng('temp/04_sim_ee_spec.png',0);
post_200709b('04','sim',3,3); mkpng('temp/04_sim_ee_ratio.png',0);
post_200709b('04','sim',4,1); mkpng('temp/04_sim_bb_spec.png',0);
post_200709b('04','sim',4,3); mkpng('temp/04_sim_bb_ratio.png',0);

