How can I plot mean and standard deviation error bars stripplot or swarmplot?

I created the following plot with the code and data posted at the end of this question:

The black dot represents the mean of the R2 Score over all retailers, and the black lines represent the corresponding standard deviation.

I want to achieve to display the mean and standard deviation in the typical way, as seen below:

I guess this must be possible with matplotlib errorbar or seaborn pointplot. But I'm working on this for ages and can not find a solution.

This answer with pointplot does not fulfill my needs, as I want one error bar over multiple categories, not one error bar per category.
I have a similar problem with this answer, working with swarmplot and pointplot.

The following is the corresponding code:

import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd

test = pd.read_csv('test.csv')

# Calculate mean and standard deviation
mean_data = test.groupby('featureset')['r2_score'].mean().values
std_data = test.groupby('featureset')['r2_score'].std().values
featuresets = ["c", "fc", "f", "s", "sc", "w"]

p = sns.stripplot(x="featureset",
                  y="r2_score",
                  hue="retailer",
                  data=test,
                  marker="^",
                  size=8)

# Plot stripplot with mean and standard deviation
sns.pointplot(x=featuresets,
              y=mean_data,
              join=False,
              color='black',
              markers='o',
              scale=2)
sns.pointplot(x=featuresets,
              y=mean_data - std_data,
              join=False,
              color='black',
              markers='_',
              scale=4)
sns.pointplot(x=featuresets,
              y=mean_data + std_data,
              join=False,
              color='black',
              markers='_',
              scale=4)

plt.legend(title='Retailer')
sns.move_legend(p, loc="upper left", bbox_to_anchor=(1, 1))

p.set(xlabel='Featureset', ylabel='R2 Score')

plt.savefig("plot.png", format="png", bbox_inches='tight')

For complete reproducibility, here add the used dataset that I named test.csv in this question:

r2_score,featureset,retailer
0.7055950484,c,S
0.942584686,c,K
0.8651950609,c,B
0.9051873402,c,H
0.5877088336,c,P
0.7944303127,c,O
0.6370605237,fc,S
0.9755270173,fc,K
0.9065356558,fc,B
0.921142567,fc,H
0.5798048892,fc,P
0.6580349995,fc,O
0.7217345443,f,S
0.9755270173,f,K
0.8839177116,f,B
0.921142567,f,H
0.5070612616,f,P
0.6580349995,f,O
0.5678318495,s,S
0.9637899061,s,K
0.9369641498,s,B
0.9297479733,s,H
0.5029283363,s,P
0.6580349995,s,O
0.5678318495,sc,S
0.9729308458,sc,K
0.8471079755,sc,B
0.9297479733,sc,H
0.497615548,sc,P
0.6580349995,sc,O
0.6624239947,w,S
0.889206858,w,K
0.7810312601,w,B
0.8562172874,w,H
0.4446346851,w,P
0.6580349995,w,O

EDIT: I updated my code to a point that fulfilles my needs better than before with the help of the answers, receiving plots in the manner of the attached example.

Please find the corresponding code below:

import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd

def plot(data, x_axis, hue, target, ordered_list=['S', 'K', 'B', 'H', 'P', 'O']):
    
    data = pd.read_csv(data)

    data = data[["r2_score", x_axis, hue]]

    # Calculate mean and standard deviation
    mean_data = data.groupby(x_axis, sort=False)['r2_score'].mean()
    std_data = data.groupby(x_axis, sort=False)['r2_score'].std()
    x = std_data.index.tolist()
    
    data_sorted = data.sort_values(hue, key=lambda x: x.map({v:k for k, v in enumerate(ordered_list)}))
    
    colorlist = ['yellowgreen', 'seagreen', 'lightseagreen', 'steelblue', 'royalblue', 'slateblue']
    
    for i in range(len(x)):
        plt.errorbar(x=i,
                     y=mean_data[i],
                     yerr=std_data[i],
                     color='grey',
                     fmt='_',
                     capsize=5,
                     elinewidth=1,
                     capthick=1)

    for i in range(len(ordered_list)):    
        p = sns.stripplot(x=x_axis,
                          y="r2_score",
                          hue=hue,
                          data=data.loc[data[hue] == ordered_list[i]],
                          marker='$' + ordered_list[i] + '$',
                          size=10,
                          palette=[colorlist[i]])

    plt.xlabel(x_axis.title(), size='xx-large')
    plt.ylabel("R2 Score", size='xx-large')
 
    p.get_legend().remove()

plot("test.csv", "featureset", "retailer", "focusproduct")

I still want to change one thing: I want that increase readability by prohibiting elements of the plot to overlap (e.g. the markers and the errorbar, or the markers among themselves). I can not find a way to do so.

You had the right idea. Errorbar works. You need the yerr and capsize arguments too.

for i, feature in enumerate(featuresets):
    plt.errorbar(x=feature, y=mean_data[i], yerr=std_data[i], color='black', fmt='_', capsize=3)

Output: