function [h,p] = plotFamilyBarGraphComparisonODnorm(fam2plot, legend2plot, SIC1_idx, SIC2_idx, SIC3_idx, SIC4_idx, L_taxonomy, L_5_relative_abundance_cleaned_forexp, L_5_abs_abundance_ODnorm_forexp) 

idx_fam = find(~cellfun(@isempty,strfind(L_taxonomy{5},fam2plot),'un',1)); % find certain family

% 14 mb plus BtWT
for i=1:3
    % create vectors for relative abundance
    SIC1_rel_reps(i,:) = L_5_relative_abundance_cleaned_forexp(SIC1_idx(i):96:288, idx_fam)';
    SIC2_rel_reps(i,:) = L_5_relative_abundance_cleaned_forexp(SIC2_idx(i):96:288, idx_fam)';
    SIC3_rel_reps(i,:) = L_5_relative_abundance_cleaned_forexp(SIC3_idx(i):96:288, idx_fam)';
    SIC4_rel_reps(i,:) = L_5_relative_abundance_cleaned_forexp(SIC4_idx(i):96:288, idx_fam)';
    
    % create vectors for absolute abundance
    SIC1_abs_reps(i,:) = L_5_abs_abundance_ODnorm_forexp(SIC1_idx(i):96:288, idx_fam)';
    SIC2_abs_reps(i,:) = L_5_abs_abundance_ODnorm_forexp(SIC2_idx(i):96:288, idx_fam)';
    SIC3_abs_reps(i,:) = L_5_abs_abundance_ODnorm_forexp(SIC3_idx(i):96:288, idx_fam)';
    SIC4_abs_reps(i,:) = L_5_abs_abundance_ODnorm_forexp(SIC4_idx(i):96:288, idx_fam)';
end

% mean relative abundance 
SIC1_rel_mean = mean(SIC1_rel_reps);
SIC2_rel_mean = mean(SIC2_rel_reps);
SIC3_rel_mean = mean(SIC3_rel_reps);
SIC4_rel_mean = mean(SIC4_rel_reps);

% std relative abundance
SIC1_rel_std = std(SIC1_rel_reps);
SIC2_rel_std = std(SIC2_rel_reps);
SIC3_rel_std = std(SIC3_rel_reps);
SIC4_rel_std = std(SIC4_rel_reps);

% mean relative abundance 
SIC1_abs_mean = mean(SIC1_abs_reps);
SIC2_abs_mean = mean(SIC2_abs_reps);
SIC3_abs_mean = mean(SIC3_abs_reps);
SIC4_abs_mean = mean(SIC4_abs_reps);

% std relative abundance
SIC1_abs_std = std(SIC1_abs_reps);
SIC2_abs_std = std(SIC2_abs_reps);
SIC3_abs_std = std(SIC3_abs_reps);
SIC4_abs_std = std(SIC4_abs_reps);
    
% make data structures for bar plot
d_rel = []; s_rel = []; d_abs = []; s_abs = [];
for i = 1:3
    % relative abundance
    d_rel = [d_rel; [SIC1_rel_mean(i) SIC2_rel_mean(i) SIC3_rel_mean(i) SIC4_rel_mean(i)]];
    s_rel = [s_rel; [SIC1_rel_std(i) SIC2_rel_std(i) SIC3_rel_std(i) SIC4_rel_std(i)]];
    
    % absolute abundance
    d_abs = [d_abs; [SIC1_abs_mean(i) SIC2_abs_mean(i) SIC3_abs_mean(i) SIC4_abs_mean(i)]];
    s_abs = [s_abs; [SIC1_abs_std(i) SIC2_abs_std(i) SIC3_abs_std(i) SIC4_abs_std(i)]];
end

% plot relative abundance
figure
b = bar(d_rel, 'grouped'); hold all;

% Calculate the number of groups and number of bars in each group
[ngroups,nbars] = size(d_rel);
% Get the x coordinate of the bars
x = nan(nbars, ngroups);
for i = 1:nbars
    x(i,:) = b(i).XEndPoints;
end
% Plot the errorbars
errorbar(x',d_rel,s_rel,'k','linestyle','none'); hold all;
xticklabels([1, 2, 3])
xlabel('passages');
ylabel('relative abundance');
ylim([0 0.60])
xlim([0.3 3.5])
h = patch([0,1.5,1.5,0,0], [0,0,1,1,0], 'r');
set(h, 'FaceAlpha', 0.1); hold all;
legend(legend2plot, 'location', 'northeast');

% plot absolute abundance
figure
b = bar(d_abs, 'grouped'); hold all;

[ngroups,nbars] = size(d_abs); % Calculate the number of groups and number of bars in each group
x = nan(nbars, ngroups); % Get the x coordinate of the bars

for i = 1:nbars
    x(i,:) = b(i).XEndPoints;
end

% Plot the errorbars
errorbar(x',d_abs,s_abs,'k','linestyle','none'); hold all;
xticklabels([1, 2, 3])
xlabel('passages');
ylabel('abs abundance (norm to final OD)');
ylim([0 0.6])
xlim([0.3 3.5])

x = [0,1.5,1.5,0,0];
y = [0,0,40,40,0];
h = patch(x,y, 'r');
set(h, 'FaceAlpha', 0.1); hold all;
legend(legend2plot, 'location', 'northeast');

h = [];
for i = 1:3
    varname = strcat('pass', string(i), 'rel');
    [h(1),p.(varname)(1,2)] = ttest(SIC1_rel_reps(1:3,i), SIC2_rel_reps(1:3,i));
    [h(2),p.(varname)(1,3)] = ttest(SIC1_rel_reps(1:3,i), SIC3_rel_reps(1:3,i));
    [h(3),p.(varname)(1,4)] = ttest(SIC1_rel_reps(1:3,i), SIC4_rel_reps(1:3,i));
    [h(7),p.(varname)(2,3)] = ttest(SIC2_rel_reps(1:3,i), SIC3_rel_reps(1:3,i));
    [h(8),p.(varname)(2,4)] = ttest(SIC2_rel_reps(1:3,i), SIC4_rel_reps(1:3,i));
    [h(9),p.(varname)(3,4)] = ttest(SIC3_rel_reps(1:3,i), SIC4_rel_reps(1:3,i));
    
    varname = strcat('pass', string(i), 'abs');
    [h(4),p.(varname)(1,2)] = ttest(SIC1_abs_reps(1:3,i), SIC2_abs_reps(1:3,i));
    [h(5),p.(varname)(1,3)] = ttest(SIC1_abs_reps(1:3,i), SIC3_abs_reps(1:3,i));
    [h(6),p.(varname)(1,4)] = ttest(SIC1_abs_reps(1:3,i), SIC4_abs_reps(1:3,i));
    [h(10),p.(varname)(2,3)] = ttest(SIC2_abs_reps(1:3,i), SIC3_abs_reps(1:3,i));
    [h(11),p.(varname)(2,4)] = ttest(SIC2_abs_reps(1:3,i), SIC4_abs_reps(1:3,i));
    [h(12),p.(varname)(3,4)] = ttest(SIC3_abs_reps(1:3,i), SIC4_abs_reps(1:3,i));
end