function sigmaData = buildupFit8stats(sigmaData)

options = optimset('display', 'off');


fprintf('\nFitting the buildup intensities\n')

% Fit sigma to cross peak intensities

for i = 1:length(sigmaData)
    
       time                                 = sigmaData(i).buildups.mixTimes;
       exp_build_up_H1H2                    = sigmaData(i).buildups.DPraw;
       exp_rho_H1                           = sigmaData(i).buildups.rhoi;
       exp_rho_H2                           = sigmaData(i).buildups.rhoj;
       exp_delta_I_zero_H1                  = sigmaData(i).buildups.Izero;
       sigmaEff_scaling                     = sigmaData(i).numcode.sigmaEffScaling;
       exp_build_up_H1H2_norm               = exp_build_up_H1H2 / exp_delta_I_zero_H1;
              
% % %        start                = -( exp_rho_H1 + exp_rho_H2 )/2;
       start                                = -0.5; 
       [uncorrected_estimated_exp_sigma, sse_target_fun_sigma, exitflag_sigma] = fminsearch(@(x)fitsigma(x, exp_rho_H1, exp_rho_H2, time ,exp_build_up_H1H2_norm, sigmaEff_scaling), start, options);
       
       
       % statistics
              
       lambda_p                             = lambdaPlusScaledFun(exp_rho_H1 ,exp_rho_H2,sigmaEff_scaling,uncorrected_estimated_exp_sigma);
       lambda_m                             = lambdaMinusScaledFun(exp_rho_H1 ,exp_rho_H2,sigmaEff_scaling,uncorrected_estimated_exp_sigma);
 
       two_spins_build_up_H1H2              = two_spins_build_up_Fun(lambda_p,lambda_m,uncorrected_estimated_exp_sigma,time,sigmaEff_scaling);
       
       diff                                 = (two_spins_build_up_H1H2 - exp_build_up_H1H2_norm).*10^3;
       chisq                                = sum((diff).^ 2) / (length(time) - 1)*10^4;
       chisq05                              = (chisq)^0.5;
       maxDP                                = max(abs(exp_build_up_H1H2_norm));
       chiN                                 = chisq05 / maxDP;
       
       
% % %        sigmaData(i).buildups.DP                         = exp_build_up_H1H2_norm;
       sigmaData(i).buildups.exp_build_up_H1H2_norm 	= exp_build_up_H1H2_norm;
       sigmaData(i).buildups.Pointtheo                  = two_spins_build_up_H1H2;
       sigmaData(i).buildups.sigmaij                    = uncorrected_estimated_exp_sigma;
       sigmaData(i).buildups.LambdaPlus                 = lambda_p;
       sigmaData(i).buildups.LambdaMinus                = lambda_m;
       sigmaData(i).buildups.diff                       = diff;  
       sigmaData(i).buildups.chisq                      = chisq;
       sigmaData(i).buildups.chisq05                    = chisq05;
       sigmaData(i).buildups.maxDP                      = maxDP;
       sigmaData(i).buildups.chiN                       = chiN;
       
end

% Write results to the sigmaData structure
for i = 1:length(sigmaData)
    
    sigmaData(i).sigmaij                    = sigmaData(i).buildups.sigmaij;  
    sigmaData(i).rhoi                       = sigmaData(i).buildups.rhoi;
    sigmaData(i).rhoj                       = sigmaData(i).buildups.rhoj;
    sigmaData(i).Izeroij                    = sigmaData(i).buildups.Izero;
    sigmaData(i).chisq05                    = sigmaData(i).buildups.chisq05;
    sigmaData(i).chiN                       = sigmaData(i).buildups.chiN;
    
end