NewtonFitter.C

Go to the documentation of this file.

// Class NewtonFitter
//
// Author: Benno List
// Last update: $Date: 2004/06/01 16:40:10 $
//          by: $Author: blist $
// 
// Description: kinematic fit using Newton method
//               

#undef NDEBUG

#include "jbltools/kinfit/NewtonFitter.h" 

#include<iostream>
#include<cmath>
#include<cassert>

#include "jbltools/kinfit/BaseFitObject.h"
#include "jbltools/kinfit/BaseConstraint.h"
#include "jbltools/kinfit/ftypes.h"
#include "jbltools/kinfit/cernlib.h"

using std::cout;
using std::cerr;
using std::endl;

// constructor
NewtonFitter::NewtonFitter() 
: npar (0), ncon (0) 
{}

// destructor
NewtonFitter::~NewtonFitter() {};

double NewtonFitter::fit() {

  // cout statements
  bool debug = false;

  // order parameters etc
  initialize();
  
  // initialize eta, etasv, y 
  int idim = npar+ncon;
  FDouble *x  = new FDouble [idim];
  FDouble *dx = new FDouble [idim];
  FDouble *y  = new FDouble [idim];
  FDouble *M  = new FDouble [idim*idim];
  
  for (int i = 0; i < idim) x[i] = y[i] = 0;
  
  // Get starting values
  for (FitObjectIterator i = fitobjects.begin(); i != fitobjects.end(); ++i) {
    BaseFitObject *fo = *i;
    assert (fo);
    for (int ilocal = 0; ilocal < fo.getNPar(); ++ilocal) {
      int iglobal = fo->getGlobalParNum (ilocal);
      assert (iglobal >= 0 && iglobal < npar);
      x[iglobal] = fo->getParam (ilocal);
    }
  }
  
  // Set initial lambda values all to 1
  for (int i = npar; i < idim; i++) x[i] = 1;
  
  // Set M to 0
  for (int i = 0; i < idim*idim) M[i] = 0;
  // Compose M: First, all terms d^2 chi^2/dx1 dx2
  for (FitObjectIterator i = fitobjects.begin(); i != fitobjects.end(); ++i) {
    BaseFitObject *fo = *i;
    assert (fo);
    fo->addToGlobalDerMatrix (idim, M);
  }
  
  
  
  
  
  // Update values in Fitobjects
  for (FitObjectIterator i = fitobjects.begin(); i != fitobjects.end(); ++i) {
    BaseFitObject *fo = *i;
    assert (fo);
    for (int ilocal = 0; ilocal < fo.getNPar(); ++ilocal) {
      int iglobal = fo->getGlobalParNum (ilocal);
      assert (iglobal >= 0 && iglobal < npar);
      fo->getParam (ilocal) = x[iglobal];
    }
  }
  

  
// *-- End of iterations - calculate errors.

//    ====> HERE GOES ERROR CALCULATION <======

  delete[] x;
  delete[] dx;
  delete[] y;
  delete[] M;
  

// *-- Turn chisq into probability.
  FReal chi = FReal(chinew);

  return prob(chi,ncon-nunm);
    
};

bool NewtonFitter::initialize() {
  // tell fitobjects the global ordering of their parameters:
  int iglobal = 0;
  // measured parameters first!
  for (int ifitobj = 0; ifitobj < fitobjects.size(); ++ifitobj) {
    for (int ilocal = 0; ilocal < fitobjects[ifitobj]->getNPar(); ++ilocal) {
      fitobjects[ifitobj]->setGlobalParNum (ilocal, iglobal);
      ++iglobal;
    }
  }
  npar = iglobal;
  int nunm = 0;
  // noe unmeasured parameters!
  for (int ifitobj = 0; ifitobj < fitobjects.size(); ++ifitobj) {
    for (int ilocal = 0; ilocal < fitobjects[ifitobj]->getNPar(); ++ilocal) {
      if (!fitobjects[ifitobj]->getMeasured(ilocal)) ++nunm;
    }
  }
  
  // set number of constraints
  ncon = constraints.size();
  
  if (nunm > ncon) {
    cerr << "NewtonFitter::initialize: nunm=" << nunm << " > ncon=" << ncon << endl;
  }
  
  return true;

};
  

---