Chisq_Constrainer.h

Go to the documentation of this file.

//
// $Id: Chisq_Constrainer.h,v 1.2 2011/10/13 09:49:48 snaumann Exp $
//
// File: hitfit/Chisq_Constrainer.h
// Purpose: Minimize a chisq subject to a set of constraints.
//          Based on the SQUAW algorithm.
// Created: Jul, 2000, sss, based on run 1 mass analysis code.
//
// For full details on the algorithm, see
//
//    @phdthesis{sssthesis,
//      author =       "Scott Snyder",
//      school =       "State University of New York at Stony Brook",
//      month =        may,
//      year =         "1995 (unpublished)"}
//    @comment{  note =         "available from {\tt http://www-d0.fnal.gov/publications\_talks/thesis/ snyder/thesis-ps.html}"
//    }
//
// CMSSW File      : interface/Chisq_Constrainer.h
// Original Author : Scott Stuart Snyder <snyder@bnl.gov> for D0
// Imported to CMSSW by Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>

#ifndef HITFIT_CHISQ_CONSTRAINER_H
#define HITFIT_CHISQ_CONSTRAINER_H

#include "TopQuarkAnalysis/TopHitFit/interface/Base_Constrainer.h"
#include "TopQuarkAnalysis/TopHitFit/interface/matutil.h"
#include <iosfwd>


namespace hitfit {


class Defaults;


class Chisq_Constrainer_Args
//
// Purpose: Hold on to parameters for the Chisq_Constrainer class.
//
// Parameters controlling the operation of the fitter:
//   bool printfit      - If true, print a trace of the fit to cout.
//   bool use_G         - If true, check the chisq formula by computing
//                        chisq directly from G.  This requires that G_i
//                        be invertable.
//
// Parameters affecting the fit:
//   float constraint_sum_eps - Convergence threshold for sum of constraints.
//   float chisq_diff_eps - onvergence threshold for change in chisq.
//   int maxit          - Maximum number of iterations permitted.
//   int max_cut        - Maximum number of cut steps permitted.
//   float cutsize      - Fraction by which to cut steps.
//   float min_tot_cutsize - Smallest fractional cut step permitted.
//
// Parameters affecting testing modes:
//   float chisq_test_eps - When use_G is true, the maximum relative
//                          difference permitted between the two chisq
//                          calculations.
//
{
public:
  // Constructor.  Initialize from a Defaults object.

  Chisq_Constrainer_Args (const Defaults& defs);

  // Retrieve parameter values.

  bool printfit () const;

  bool use_G () const;

  double constraint_sum_eps () const;

  double chisq_diff_eps () const;

  unsigned  maxit () const;

  unsigned  max_cut () const;

  double cutsize () const;

  double min_tot_cutsize () const;

  double chisq_test_eps () const;

  // Arguments for subobjects.

  const Base_Constrainer_Args& base_constrainer_args () const;


private:
  // Hold on to parameter values.

  bool _printfit;

  bool _use_G;

  double _constraint_sum_eps;

  double _chisq_diff_eps;

  int  _maxit;

  int  _max_cut;

  double _cutsize;

  double _min_tot_cutsize;

  double _chisq_test_eps;


  const Base_Constrainer_Args _base_constrainer_args;
};


//*************************************************************************


class Chisq_Constrainer
//
// Purpose: Minimize a chisq subject to a set of constraints.
//          Based on the SQUAW algorithm.
//
  : public Base_Constrainer
{
public:
  // Constructor, destructor.
  // ARGS holds the parameter settings for this instance.

  Chisq_Constrainer (const Chisq_Constrainer_Args& args);

  virtual ~Chisq_Constrainer () {}

  // Do the fit.
  // Call the number of well-measured variables Nw, the number of
  // poorly-measured variables Np, and the number of constraints Nc.
  // Inputs:
  //   CONSTRAINT_CALCULATOR is the object that will be used to evaluate
  //     the constraints.
  //   XM(Nw) and YM(Np) are the measured values of the well- and
  //     poorly-measured variables, respectively.
  //   X(Nw) and Y(Np) are the starting values for the fit.
  //   G_I(Nw,Nw) is the error matrix for the well-measured variables.
  //   Y(Np,Np) is the inverse error matrix for the poorly-measured variables.
  //
  // Outputs:
  //   X(Nw) and Y(Np) is the point at the minimum.
  //   PULLX(Nw) and PULLY(Np) are the pull quantities.
  //   Q(Nw,Nw), R(Np,Np), and S(Nw,Np) are the final error matrices
  //     between all the variables.
  //
  // The return value is the final chisq.  Returns a value < 0 if the
  // fit failed to converge.

  virtual double fit (Constraint_Calculator& constraint_calculator,
                      const Column_Vector& xm,
                      Column_Vector& x,
                      const Column_Vector& ym,
                      Column_Vector& y,
                      const Matrix& G_i,
                      const Diagonal_Matrix& Y,
                      Column_Vector& pullx,
                      Column_Vector& pully,
                      Matrix& Q,
                      Matrix& R,
                      Matrix& S);

  // Print out any internal state to S.
  virtual std::ostream& print (std::ostream& s) const;


private:
  // Parameter settings.
  const Chisq_Constrainer_Args _args;
};

} // namespace hitfit


#endif // not HITFIT_CHISQ_CONSTRAINER_H