function out = fitSpMuZ(X,s,p,m,mask,plotFig)

% Computes misfit between spallation + muon profile and 
% measurements for GISP2 data. Zero erosion. 
% We allow multipliers (e.g., apparent ages) to float for Be-10 and Al-26 separately.
% Depth is also allowed to float, i.e. we want to know what depth makes
% the measurements look like the production profile. 
%
% syntax is out = fitSpMuZ(X,s,p,m,mask,plotFig)
% 
% So X has [t10 t26 Z] (units: ka ka g/cm2) 
% s is sample data structure
% p is production rate data structure. 
% m is muon data structure. 
% mask is vector of data numbers to consider, normally [1:7]
% plotFig is output of plotGISPdata.m
%
% Greg Balco -- Berkeley Geochronology Center -- May 2016



for a = 1:length(s);
    this_td = s{a}.tdg + X(3); 
    this_bd = s{a}.bdg + X(3);
    this_mq = s{a}.mq;
    
    if s{a}.nuclide == 10; 
        pfun = @(z,t) PofZ(z,m,p,10); % Simple production
        texp = X(1).*1e3;
    else
        pfun = @(z,t) PofZ(z,m,p,26); % Simple production
        texp = X(2).*1e3;
    end;

    % Integrate
    
    TempTotN = zeros(1,length(this_td));
    for b = 1:length(this_td);
        % Integral of production rate in interval divided by interval
        % thickness
        TempTotN(b) = integral2(pfun,this_td(b),this_bd(b),0,texp,'RelTol',1e-3,'AbsTol',1e-3) ./ (this_bd(b) - this_td(b));
    end;
    
    % Now account for quartz yields
    NpAvg(a) = (sum(TempTotN.*this_mq))./sum(this_mq);
    
    % Compute miss
    miss(a) = ((NpAvg(a) - s{a}.N)./s{a}.dN).^2;
    
end;

if ~isempty(plotFig);
    axes(plotFig(2));
    for a = 1:length(s);
        xx = [NpAvg(a) NpAvg(a)];
        for b = 1:length(s{a}.td);
            yy = [s{a}.td(b) s{a}.bd(b)];
            if isempty(findobj('tag',['predN' int2str(a) int2str(b)])); 
                plot(xx,yy,'g','tag',['predN' int2str(a) int2str(b)]);
                set(findobj('tag',['predN' int2str(a) int2str(b)]),'LineWidth',2);
            else
                set(findobj('tag',['predN' int2str(a) int2str(b)]),'xdata',xx,'ydata',yy);
            end;
        end;
    end;
   
    % Also plot continuous p(z)
    pz = 0:0.05:2; pzm = pz.*100.*2.7 + X(3);
    pN10 = X(1).*1e3.*PofZ(pzm,m,p,10);
    pN26 = X(2).*1e3.*PofZ(pzm,m,p,26);
    if isempty(findobj('tag','predN10'));
        plot(pN10,pz,'g','tag','predN10');
    else
        set(findobj('tag','predN10'),'xdata',pN10);
    end;
    if isempty(findobj('tag','predN26'));
        plot(pN26,pz,'g','tag','predN26');
    else
        set(findobj('tag','predN26'),'xdata',pN26);
    end;
    
    drawnow;

end;

out = sum(miss(mask));