87 lines
1.8 KiB
C
87 lines
1.8 KiB
C
#include "lib.h"
|
|
|
|
|
|
// Minimization wrapper.
|
|
//
|
|
double chi2(double p_max, void* params)
|
|
{
|
|
// Pars parameters.
|
|
struct parameters p = *((struct parameters*) params);
|
|
struct bin* histo = p.histo;
|
|
size_t n = p.n;
|
|
double step = p.step;
|
|
|
|
// Compute χ².
|
|
double chi = 0;
|
|
double expecto = 0;
|
|
for (size_t i = 0; i < n; i++)
|
|
{
|
|
expecto = expected((i + 0.5) * step, p_max);
|
|
chi += pow(histo[i].sum - expecto, 2)/expecto;
|
|
};
|
|
return chi;
|
|
}
|
|
|
|
|
|
// Expected function.
|
|
//
|
|
double expected (double x, double p_max)
|
|
{
|
|
// When p_max < x, the argument under sqrt is negative. Return great number
|
|
// to lead away the minimization.
|
|
if (p_max < x)
|
|
{
|
|
double num = 500;
|
|
return num;
|
|
}
|
|
return x * log(p_max/x)/atan(sqrt(fabs(pow(p_max, 2)/pow(x,2) - 1)));
|
|
}
|
|
|
|
|
|
// First derivative of the expected function.
|
|
//
|
|
double exp1d (double x, double p_max)
|
|
{
|
|
// When p_max < x, the argument under sqrt is negative. Return great number
|
|
// to lead away the minimization.
|
|
if (p_max < x)
|
|
{
|
|
double num = 500;
|
|
return num;
|
|
}
|
|
double a = p_max;
|
|
double A = sqrt(fabs(pow(a, 2)/pow(x,2) - 1));
|
|
double B = atan(A);
|
|
double C = log(a/x);
|
|
return x/(a*B) - x*C/(a*A*pow(B, 2));
|
|
}
|
|
|
|
//
|
|
// Second derivative of the expected function.
|
|
//
|
|
double exp2d (double x, double p_max)
|
|
{
|
|
// When p_max < x, the argument under sqrt is negative. Return great number
|
|
// to lead away the minimization.
|
|
if (p_max < x)
|
|
{
|
|
double num = 500;
|
|
return num;
|
|
}
|
|
double a = p_max;
|
|
double a2 = pow(a, 2);
|
|
|
|
double A = sqrt(fabs(a2/pow(x,2) - 1));
|
|
double B = atan(A);
|
|
double C = log(a/x);
|
|
double D = 1 - a2/pow(x,2);
|
|
|
|
double E = -2/(a2*D*B);
|
|
double F = 1/(a2*A);
|
|
double G = 1/(pow(x, 2)*pow(A, 3));
|
|
double H = 1/a2 + 2/(a2 * A * B);
|
|
|
|
return x*((E + F + G) * C/B - H)/B;
|
|
}
|
|
|