NMFit Namespace Reference

Functions
void	Fit (bool mwt, const std::vector< double > x, const std::vector< double > y, const std::vector< double > sig, double &a, double &b, double &siga, double &sigb, double &chi2, double &q)
double	gammq (double a, double x)
void	gser (double &gamser, double a, double x, double &gln)
void	gcf (double &gammcf, double a, double x, double &gln)
double	gammln (double xx)
double	FitPol (int pow, int n, const double *x, const double *y, const double *w, double *par, double *pErr)
	Fits points to a polynomial of a given degree, returns chi squared, and if desired, errors on parameters.

---

Function Documentation

void NMFit::Fit	(	bool	mwt,
		const std::vector< double >	x,
		const std::vector< double >	y,
		const std::vector< double >	sig,
		double &	a,
		double &	b,
		double &	siga,
		double &	sigb,
		double &	chi2,
		double &	q
	)

Definition at line 23 of file NMFit.cxx.

References gammq().

Referenced by bestChi2(), ComputeADCOffset(), main(), and RosieFCN().

                                 {

  double wt,t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss,sigdat;

  if (x.empty()) return;

  std::cout << "I am inside NMFit::Fit()" << std::endl;

  b=0.0;
  if (mwt) {
    ss=0.0;
    for (unsigned int i=0; i<x.size(); i++) {
      std::cout << "With weights: x=" << x[i] << " z=" << y[i]
        << " sig=" << sig[i] << std::endl;
      wt=1.0/(sig[i]*sig[i]);
      ss = ss + wt;
      sx = sx + x[i]*wt;
      sy = sy + (y[i]+4000.0)*wt;
    }
  }
  else {
    for (unsigned int i=0; i<x.size(); i++) {
      std::cout << "W/out weights: x=" << x[i] << " z=" << y[i]
        << " sig=" << sig[i] << std::endl;
      sx = sx + x[i];
      sy = sy + (y[i]+4000.0);
    }
    ss=(double)x.size();
  }
  sxoss=sx/ss;

  std::cout << "ss/sx/sy/sxoss=" 
        << ss << "/" << sx  << "/" << sy  << "/" << sxoss << std::endl;

  if (mwt) {
    for (unsigned int i=0; i<x.size(); i++) {
      t = (x[i]-sxoss)/sig[i];
      st2 = st2 + t*t;
      b = b + t*y[i]/sig[i];
    }
  }
  else {
    for (unsigned int i=0; i<x.size(); i++) {
      t = x[i] - sxoss;
      st2 = st2 + t*t;
      b = b + t*y[i];
    }
  }

  b = b/st2;
  a = (sy-sx*b)/ss;
  siga = sqrt((1.0+sx*sx/(ss*st2))/ss);
  sigb = sqrt(1.0/st2);
  chi2 = 0.0;
  if (!mwt) {
    for (unsigned int i=0; i<x.size(); i++)
      chi2 = chi2 + ((y[i]-a-b*x[i])*(y[i]-a-b*x[i]));
    q = 1.0;
    sigdat = sqrt((chi2)/(x.size()-2));
    siga = siga*sigdat;
    sigb = sigb*sigdat;

    // tmp output
    std::cout << "Finished NMFit (no weights): slope/intercept=" 
          << a << "/" << b 
          << "\n                   with errors= "
          << siga << "/" << sigb 
          << "\n                   with chi^2/q= "
          << chi2 << "/" << q
          << std::endl;
  }
  else {
    for (unsigned int i=0; i<x.size(); i++)
      chi2 = chi2 + ((y[i]-a-b*x[i])/sig[i])*((y[i]-a-b*x[i])/sig[i]);
    q=gammq(0.5*(x.size()-2),0.5*chi2);

    // tmp output
    std::cout << "Finished NMFit (w/weights): slope/intercept=" 
          << a << "/" << b 
          << "\n                   with errors= "
          << siga << "/" << sigb 
          << "\n                   with chi^2/q= "
          << chi2 << "/" << q
          << std::endl;
  }
}

double NMFit::gammq	(	double	a,
		double	x
	)

Definition at line 118 of file NMFit.cxx.

References gcf(), and gser().

Referenced by Fit().

                                      {

  double gamser,gammcf,gln;

  if (x < 0.0 || a <= 0.0) {
    std::cout << "Invalid arguments in routine GAMMQ" << std::endl;
  }
  if (x < (a+1.0)) {
    gser(gamser,a,x,gln);
    return 1.0-gamser;
  }
  else {
    gcf(gammcf,a,x,gln);
    return gammcf;
  }
}

void NMFit::gser	(	double &	gamser,
		double	a,
		double	x,
		double &	gln
	)

Definition at line 138 of file NMFit.cxx.

References EPS, gammln(), and ITMAX.

Referenced by gammq().

                                                             {

  int n;
  double sum,del,ap;

  gln=gammln(a);
  if (x <= 0.0) {
    if (x < 0.0) std::cout << "x less than 0 in routine GSER" << std::endl;
    gamser=0.0;
    return;
  }
  else {
    ap=a;
    sum=1.0/a;
    del=1.0/a;
    for (n=1;n<=ITMAX;n++) {
      ap += 1.0;
      del *= x/ap;
      sum += del;
      if (fabs(del) < fabs(sum)*EPS) {
    gamser=sum*exp(-x+a*log(x)-gln);
    return;
      }
    }
    std::cout << "a too large, ITMAX too small in routine GSER" << std::endl;
    return;
  }
}

void NMFit::gcf	(	double &	gammcf,
		double	a,
		double	x,
		double &	gln
	)

Definition at line 170 of file NMFit.cxx.

References EPS, gammln(), and ITMAX.

Referenced by gammq().

                                                            {

  int n;
  double gold=0.0,g,fac=1.0,b1=1.0;
  double b0=0.0,anf,ana,an,a1,a0=1.0;

  gln=gammln(a);
  a1=x;
  for (n=1;n<=ITMAX;n++) {
    an=(double) n;
    ana=an-a;
    a0=(a1+a0*ana)*fac;
    b0=(b1+b0*ana)*fac;
    anf=an*fac;
    a1=x*a0+anf*a1;
    b1=x*b0+anf*b1;
    if (a1) {
      fac=1.0/a1;
      g=b1*fac;
      if (fabs((g-gold)/g) < EPS) {
    gammcf=exp(-x+a*log(x)-gln)*g;
    return;
      }
      gold=g;
    }
  }
  std::cout << "a too large, ITMAX too small in routine GCF" << std::endl;
}

double NMFit::gammln

(

double

xx

)

Definition at line 202 of file NMFit.cxx.

References x.

Referenced by gcf(), and gser().

                              {

  double x,tmp,ser;
  static double cof[6]={76.18009173,-86.50532033,24.01409822,
            -1.231739516,0.120858003e-2,-0.536382e-5};
  int j;

  x=xx-1.0;
  tmp=x+5.5;
  tmp -= (x+0.5)*log(tmp);
  ser=1.0;
  for (j=0;j<=5;j++) {
    x += 1.0;
    ser += cof[j]/x;
  }
  return -tmp+log(2.50662827465*ser);
}

double NMFit::FitPol	(	int	pow,
		int	n,
		const double *	x,
		const double *	y,
		const double *	w,
		double *	par,
		double *	pErr = 0
	)

Fits points to a polynomial of a given degree, returns chi squared, and if desired, errors on parameters.

If w is 0 (NULL pointer), then the errors are all set to 1

Definition at line 225 of file NMFit.cxx.

References ok, and p.

Referenced by AlignChamW0::AdjustBField(), TOFTempCorrector::FitGraph(), FitGraph(), VtxConFit::FitVertex(), main(), and RosieFCN().

{
  pow++;
  memset(par, 0, pow * sizeof(double));

  const int kMaxDegree = 16;

  if (pow > kMaxDegree) {
    cout << "NMFit::Pol() does not support " << pow - 1
     << "-degree polynomial" << endl;
    exit(0);
  }
  if (n < pow) {
    cout << "NMFit::Pol() cannot do a " << pow - 1 << "-degree polynomial fit"
     << " with " << n << " points" << endl;
    return 0;
  }

  double mat[kMaxDegree * kMaxDegree];
  double vec[kMaxDegree];

  memset(mat, 0, kMaxDegree * kMaxDegree * sizeof(double));
  memset(vec, 0, kMaxDegree * sizeof(double));

  double val1; // w * x^(2j)
  double val2; // w * y * x^j
  double val3; // w * x^j * x^k
  double chi2 = 0; // chi^2 = sum_i(w * y^2) - p * v

  // Compute elements of symmetric matrix and vector
  // m[j][k] = sum_i (w[i] * x[i]^j * x[i]^k)
  // v[j] = sum_i (w[i] * y[i] * x[i]^j)
  //
  for (int i = 0; i < n; ++i) {
    if (w) val1 = fabs(w[i]); // Make sure that weights are definite positive
    else val1 = 1.0;

    val2 = val1 * y[i];
    chi2 += val2 * y[i]; 

    for (int j = 0; j < pow; ++j) {
      vec[j] += val2;
      val3 = val1;
      for (int k = j; k < pow; ++k) {
    mat[j * pow + k] += val3;
    mat[k * pow + j] = mat[j * pow + k];
    val3 *= x[i];
      }

      val1 *= x[i] * x[i]; // We go from element [i][i] to [i+1][i+1], or 2 powerso f x
      val2 *= x[i]; 
    }

  }

  TMatrixDSym m(pow, mat);
  TVectorD v(pow, vec);

  TDecompLU lu(m);
  lu.Decompose();
  Bool_t ok;
  TVectorD p = lu.Solve(v, ok);

  // If errors on parameters are requested, we need matrix inverse
  if (pErr) {
    TMatrixD inv(pow, pow);
    lu.Invert(inv);
    for (int i = 0; i < pow; ++i) pErr[i] = sqrt(inv[i][i]);
  }

  // Copy parameters and finish computation of chi squared
  for (int i = 0; i < pow; ++i) {
    par[i] = p[i];
    chi2 -= par[i] * vec[i];
  }

  return chi2;
}

---