How to create normal distribution from given data in C++?

I am trying to create a class, that is able to take numbers, calculate stddev, mean etc. for normal distribution and then, on method call with argument of confidence interval, return pair of values, which represent the borders of that confidence interval.

int main() {
    double arr[] = { 1, 2, 2, 3, 3, 3, 4, 4, 5 };
    int size = 9;
    std::vector<double> data;

    data.reserve( size );
    for ( int i = 0 ; i < size ; ++ i ) {
        data.push_back( arr[i] );
    }

    // 95 is percentage
    // return value should be a pair of values representing borders of that 95 % stripe in normal distribution
    auto pair = NormalDistribution( data ).borders_for_confidence_level( 95 );

    std::cout << "bottom: " << pair.first << std::endl
              << "top: " << pair.second << std::endl;

    return 0;
}

Currently, I have this:

class NormalDistribution {
public:
    explicit NormalDistribution( const std::vector<double> &data )
            : m_data_size( data.size() ), m_mean( mean( data ) ) {
        m_stddev = standard_deviation( data, m_mean );

        std::cout << "mean: " << m_mean << std::endl
                  << "stddev: " << m_stddev << std::endl;
    }

    std::pair<double, double> borders_for_confidence_level( double confidence_percentage ) const {
        double stddev_coefficient = stddev_coefficient_for_confidence( confidence_percentage );
        double delta = stddev_coefficient * m_stddev / sqrt( m_data_size );

        std::cout << "coeff: " << stddev_coefficient << std::endl;

        return std::make_pair( m_mean - delta, m_mean + delta );
    }

private:
    double m_stddev, m_mean;
    int m_data_size;

    static double standard_deviation( const std::vector<double> &data, double mu ) {
        return sqrt( variance( data, mu ) );
    }

    static double stddev_coefficient_for_confidence( double percentage ) {
        return 1.96;
    }

    static double mean( const std::vector<double> &data ) {
        return sum( data ) / data.size();
    }

    static double variance( const std::vector<double> &data, double mu ) {
        std::vector<double> squared_diff;

        squared_diff.reserve( data.size() );
        for ( auto n : data ) {
            squared_diff.push_back( pow( mu - n, 2 ) );
        }

        return mean( squared_diff );
    }

    static double sum( const std::vector<double> &data ) {
        double total = .0;

        for ( const auto &n : data ) {
            total += n;
        }

        return total;
    }
};

Current output from this is

mean: 3
stddev: 1.1547
coeff: 1.96
bottom: 2.2456
top: 3.7544

which looks correct.

What I need to do now is to implement

static double stddev_coefficient_for_confidence( double percentage ) {
    return 1.96; // this is hard-coded coefficient for 95 %
}

to convert confidence into the stddev coefficient.

I tried to use erfc, but it converts stddev to percentage. I need the opposite behaviour.

Any idea how to solve that?

Thank you.