This is minuitFCN.c in view mode; [Download] [Up]
/*
* These are the C routines which handle evaluating the Chi^2 for Minuit
*
* Paul Rensing, Jan 27, 1993
*
* $Id: minuitFCN.c,v 1.6 1993/03/23 22:50:18 rensing Exp $
*/
#include <math.h>
#include <float.h>
#include "minuitFCN.h"
const char minuitFCN_c_rcsid[] =
"$Id: minuitFCN.c,v 1.6 1993/03/23 22:50:18 rensing Exp $";
const char minuitFCN_h_rcsid[] = MINUITFCN_H_RCSID;
#define EPS 3.0e-11
#define ERRLIM 1e-6
#define ABSERRLIM 1e-20
static double xyChi2(FCNParms *fcnParms);
static double intHistoChi2(FCNParms *fcnParms, double lastParm);
static double likelyhood(FCNParms *fcnParms);
static double binInt(display dsp, double xlow, double xhigh );
static double qgaus(display dsp, double xlow, double xhigh, int order);
static void gauleg(double x1,double x2,double *x,double *w,int n);
void minuitfcn_(double *chi2, double *parm, int *npar, FCNParms *fcnParms)
{
/*
* First, get the InspectPFunc to update the parameter values.
*/
loadParams(parm);
switch (fcnParms->fitType)
{
case FIT_ONESAMPLE:
*chi2 = xyChi2(fcnParms);
break;
case FIT_INTEGRAL:
case FIT_INTFIXEDNORM:
case FIT_AVERAGE:
*chi2 = intHistoChi2(fcnParms, parm[*npar-1]);
break;
case FIT_MAXLIKE:
*chi2 = likelyhood(fcnParms);
break;
}
}
/*
* Standard Chi^2 calculation
*/
static double xyChi2(FCNParms *fcnParms)
{
double sum = 0.0;
double dy;
/* do this crud for optimization. */
display dsp = fcnParms->disp;
const float *x = fcnParms->x;
const float *y = fcnParms->y;
const float *yerr = fcnParms->yerr;
const float *last = &fcnParms->x[fcnParms->npts];
for (; x < last; x++, y++, yerr++) {
dy = h_funcSum(dsp, *x) - *y;
if (*yerr) {
dy /= (*yerr);
sum += dy * dy;
} else if (fcnParms->ebType == EB_ONE) {
sum += dy * dy;
}
}
return sum;
}
/*
* Chi^2 calculation: y value is integration over bin.
*/
static double intHistoChi2(FCNParms * fcnParms, double lastParm)
{
double sum = 0.0;
double dy, yint;
double ysum = 0.0;
double ydsum = 0.0;
/* do this crud for optimization. */
display dsp = fcnParms->disp;
const float *x = fcnParms->x;
const float *y = fcnParms->y;
const float *yerr = fcnParms->yerr;
const float *last = &fcnParms->x[fcnParms->npts];
const float *xerr = fcnParms->xerr;
float binWidth;
if (h_getDispType(dsp) == HISTOGRAM)
binWidth = h_getBinWidth(dsp, XAXIS);
else
binWidth = 1.0;
for (; x < last; x++, y++, yerr++, xerr++) {
/*
* do the integral. Divide by binWidth because the function
* multiplies by it to simulate the integral.
*/
yint = binInt(dsp, *x - *xerr, *x + *xerr) / binWidth;
if (fcnParms->fitType == FIT_AVERAGE)
yint /= (2.0 * (*xerr));
if (*yerr) {
dy = (yint - *y) / (*yerr);
sum += dy * dy;
} else if (fcnParms->ebType == EB_ONE) {
dy = yint - *y;
sum += dy * dy;
}
ysum += yint; /* keep sums needed for norm constraint */
ydsum += *y;
}
/*
* If the normalization is fixed to sum of bins, add constraint using
* Lagrange multiplier.
*/
if (fcnParms->fitType == FIT_INTFIXEDNORM)
sum += 2.0 * lastParm * (ysum - ydsum);
return sum;
}
static double binInt(display dsp, double xlow, double xhigh )
{
int order = 10;
double sum, oldsum;
oldsum = qgaus(dsp,xlow,xhigh,order++);
sum = qgaus(dsp,xlow,xhigh,order++);
while ( order < 40
&& (fabs(oldsum - sum) > ERRLIM*fabs(sum)
&& fabs(oldsum - sum) > ABSERRLIM )
)
{
oldsum = sum;
sum = qgaus(dsp,xlow,xhigh,order++);
}
return sum;
}
static double qgaus(display dsp, double xlow, double xhigh, int order)
{
double x[order], w[order];
int j;
double s = 0.0;
gauleg(xlow,xhigh,x,w,order);
for (j=0; j<order; j++)
s += w[j]*h_funcSum(dsp, x[j]);
return s;
}
/*
* Routine to calculate the Gauss-Legendre weights and abcissas for
* quadrature integration.
*
* from Numerical Recipes in C.
*/
static void gauleg(double x1,double x2,double *x,double *w,int n)
{
int m,j,i;
double z1,z,xm,xl,pp,p3,p2,p1;
m=(n+1)/2;
xm=0.5*(x2+x1);
xl=0.5*(x2-x1);
for (i=1;i<=m;i++) {
z=cos(3.141592654*(i-0.25)/(n+0.5));
do {
p1=1.0;
p2=0.0;
for (j=1;j<=n;j++) {
p3=p2;
p2=p1;
p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
}
pp=n*(z*p1-p2)/(z*z-1.0);
z1=z;
z=z1-p1/pp;
} while (fabs(z-z1) > EPS);
x[i-1]=xm-xl*z;
x[n-i]=xm+xl*z;
w[i-1]=2.0*xl/((1.0-z*z)*pp*pp);
w[n-i]=w[i-1];
}
}
/*
* Maximum likelihood function
*/
static double likelyhood(FCNParms *fcnParms)
{
return 0.0;
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.