After you have calculated the relevant parameters of your fit with the help of the Gauss snippet, you may want to calculate the quality of your fit with a Chi-squared test.

The program cl.cpp calculates the confidence level from the χ2 and the number of constraints. The version included here is the original FORTRAN routine copied from Lau Gatignon (1980) with some later modifications.

/***************************************************
 * Computation of the confidence level from Chi-squared (chi2)
 * and number of constraints (ncons).
 * 
 * Note :	
 * for even ncons ==> same result as the Cernlib function PROB
 * for odd	ncons ==> confidence level < result of the Cernlib function PROB
 * 
 * --- Original FORTRAN routine copied from Lau Gatignon 1980
 * --- Modified by NvE 29-sep-1980 K.U. Nijmegen
 * --- Modified by NvE 24-apr-1985 NIKHEF-H Amsterdam
 * --- Converted to C++ by Nve 06-nov-1998 UU-SAP Utrecht
 *
 *  - Compiling:   g++ cl.cpp -o cl
 *  - Running:     ./cl
 ***************************************************/

#include <math.h>
#include <cstdlib>
#include <iostream>

using namespace std;

double conlev(double chi2, int ncons) {
	const double a1 = 0.705230784e-1, a2 = 0.422820123e-1,
                 a3 = 0.92705272e-2,  a4 = 0.1520143e-3,
                 a5 = 0.2765672e-3,   a6 = 0.430638e-4;

	// Change to dummy variables
	int n     = ncons;
	double y0 = chi2;

	// Set CL to zero in case ncons <= 0
	if (n <= 0)
		return 0;

	if (y0 <= 0.) {
		if (y0 < 0.) 
			return 0;
		else 
			return 1;
	}

	if (n > 100) {
		double x     = sqrt(2.*y0) - sqrt(double(n + n - 1));
		double t     = fabs(x) / 1.1412135;
		double denom = 1.0 + t*(a1 + t*(a2 + t*(a3 + t*(a4 + t*(a5 + t*a6)))));
		double vfun  = 1.0 / denom;

		// Prevent underflow
		if (fabs(vfun) < 1.3e-5)
			vfun = 0.; 

		vfun     = pow(vfun, 16);
		double v = 0.5*vfun;

		if (x < 0.)
			v = 1.-v;

		if (v < 0.)
			v = 0.;

		return v;
	}

	else {
		double y1 = 0.5*y0;
		int m     = n / 2;

		if (2*m == n) {
		 double sum=1.;

		if (m > 1) {
			double term=1.;
			
			for (int i=2; i<=m; i++) {
				// Prevent overflow
				if (term > 1.e20)
					break;
 
				 double fi = double(i-1);
				 term      = term * y1 / fi;
				 sum      += term;
			}
		}

		double rnick = y1;

		// Prevent underflow
		if (rnick >= 1.75e2)
			rnick  = 1.75e2; 

		double v = sum*exp(-rnick);

		if (v < 0.)
			v = 0.;

			return v;
		}

		else {
			double x     = sqrt(y0);
			double t     = fabs(x) / 1.1412135;
			double denom = 1.0 + t*(a1 + t*(a2 + t*(a3 + t*(a4 + t*(a5 + t*a6)))));
			double vfun  = 1.0 / denom;

			// Prevent underflow
			if (fabs(vfun) < 1.3e-5)
				vfun=0.; 

			vfun     = pow(vfun,16);
			double v = vfun;

			if (n < 1)
				return 0;

			if (n == 1) {
				if (v < 0.)
					v=0.;
				return v;
			}

			double value = v;
			double sum   = 1.; 

			if (n >= 4) {
				double term = 1.;
				int k       = m - 1;

				for (int j=1; j<=k; j++) {
					// Prevent overflow
					if (term > 1.e20)
						break;

					double fj = double(j);
					term      = term * y0 / (fj + fj + 1.);
					sum      += term;
				}
			}

			// Prevent underflow
			if (y1 > 1.75e2)
				y1 = 1.75e2; 

			double vi = sum * 0.797846 * sqrt(y0) * exp(-y1);
			v         = vi + value;

			if (v < 0.)
				v = 0.;
	
			return v;
		}
	}
}


int main() {
	cout << conlev(1,2) << endl;
}