NTupleFCN.cxx

Go to the documentation of this file.

#ifdef _MSC_VER
#include "msdevstudio/MSconfig.h"
#endif

#include "NTupleFCN.h"

#include "datasrcs/DataPointTuple.h"
#include "datasrcs/DataSource.h"
#include "datasrcs/TupleCut.h"

#include "functions/FunctionBase.h"

#include <algorithm>
#include <functional>

using std::bind2nd;
using std::count;
using std::find_if;
using std::not_equal_to;
using std::vector;

using namespace hippodraw;

NTupleFCN::
NTupleFCN ( )
  :  m_fit_cut ( 0 ),
     m_fit_range ( false ),
     m_ntuple ( 0 )
{
}

NTupleFCN::
NTupleFCN ( const NTupleFCN & fcn )
  : StatedFCN ( fcn ),
    m_fit_cut ( 0 ),
    m_fit_range ( fcn.m_fit_range ),
    m_ntuple ( fcn.m_ntuple )
{
}

NTupleFCN::
NTupleFCN ( const DataSource * ntuple, FunctionBase * function )
  : StatedFCN ( function ),
    m_fit_cut ( 0 ),
    m_fit_range ( false ),
    m_ntuple ( ntuple )
{
}

namespace dp2 = hippodraw::DataPoint2DTuple;
namespace dp3 = hippodraw::DataPoint3DTuple;

void
NTupleFCN::
setDataSource ( const DataSource * ntuple )
{
  unsigned int size = dp2::SIZE;
  if ( ntuple -> columns () == dp3::SIZE ) size = dp3::SIZE;
  vector < int > indices ( size );

  if ( size == dp2::SIZE ) {
    indices [ dp2::X ] = dp2::X;
    indices [ dp2::Y ] = dp2::Y;
    indices [ dp2::XERR ] = dp2::XERR;
    indices [ dp2::YERR ] = dp2::YERR;
  }

  if ( size == dp3::SIZE ) {
    indices [ dp3::X ] = dp3::X;
    indices [ dp3::Y ] = dp3::Y;
    indices [ dp3::Z ] = dp3::Z;
    indices [ dp3::XERR ] = dp3::XERR;
    indices [ dp3::YERR ] = dp3::YERR;
    indices [ dp3::ZERR ] = dp3::ZERR;
  }
  int dimension = size == dp2::SIZE ? 1 : 2;

  setDataSource ( ntuple, dimension, indices );
}

void
NTupleFCN::
setDataSource ( const DataSource * ntuple, 
                int dimension,
                const std::vector < int > & indices )
{
  m_ntuple = ntuple;
  m_indices = indices;
  m_dimen = dimension;
}

int
NTupleFCN::
getErrorColumn () const
{
  int ie = -1;
  if ( m_dimen == 1 ) {
    ie = m_indices [ dp2::YERR ];
  }
  else {
    ie = m_indices [ dp3::ZERR ];
  }
  return ie;
}

bool
NTupleFCN::
hasErrors ( ) const
{
  int ie = getErrorColumn ();
  const vector < double > & errors = m_ntuple -> getColumn ( ie );

  if ( errors.empty() ) return false;

  vector < double >::const_iterator first
    = find_if ( errors.begin(), errors.end (),
                bind2nd ( not_equal_to < double > (), 0.0 ) );

  return first != errors.end ();
}

bool
NTupleFCN::
setUseErrors ( bool yes )
{
  bool didit = false;
  if ( yes ) {
    if ( hasErrors () ) {
      m_has_errors = true;
      didit = true;
    }
    else m_has_errors = false;
  }
  else {
    m_has_errors = false;
      didit = true;
  }
  return didit;
}

int
NTupleFCN::
degreesOfFreedom () const
{
  int ie = getErrorColumn ();

  const vector < double > & errors = m_ntuple -> getColumn ( ie );
  int number_points = errors.size();
  if ( m_has_errors ) {
    int zeros = count ( errors.begin(), errors.end(), 0.0 );
    number_points -= zeros;
  }

  vector< double > free_parms;
  return number_points - getNumberFreeParms ();
}

void 
NTupleFCN::
reset ( std::vector < std::vector < double > > & alpha,
        std::vector < double > & beta,
        unsigned int size )
{
  beta.clear ();
  beta.resize ( size, 0.0 );

  alpha.resize ( size );

  for ( unsigned int i = 0; i < alpha.size (); i++ ) {
    alpha[i].clear ();
    alpha[i].resize ( size, 0.0 );
  }
}

void
NTupleFCN::
calcAlphaBeta ( std::vector < std::vector < double > > & alpha,
                std::vector < double > & beta )
{
  int ix = m_indices [ dp2::X ];
  int iy = m_indices [ dp2::Y ];
  int ie = m_indices [ dp2::YERR ];
  unsigned int num_parms = getNumberFreeParms ();
  reset ( alpha, beta, num_parms );

  unsigned int rows = m_ntuple -> rows ();
  for ( unsigned int i = 0; i < rows; i++ ) {
    if ( acceptRow ( i ) ) {
      const vector < double > & row = m_ntuple -> getRow ( i );

      double err = ie < 0 ? 0. : row [ ie ];
      if ( err == 0.0 && m_has_errors ) continue;
      if ( m_has_errors == false ) err = 1.0;

      double x = row [ ix ];
      double y = row [ iy ];

      double y_diff = y - m_function -> operator () ( x );
      vector < double > derives;
      fillFreeDerivatives ( derives, x );

      for ( unsigned int j = 0; j < num_parms; j++ ) {
        double t = derives[j] / ( err * err );

        for ( unsigned int k = 0; k <= j; k++ ) {
          alpha[j][k] = alpha[j][k] + t * derives[k];
        }

        beta[j] += t * y_diff;
      }
    }
  }
}

void
NTupleFCN::
setFitCut ( TupleCut * cut )
{
  int ix = m_indices [ dp2::X ];
  cut -> setColumn ( ix );

  m_fit_cut = cut;
}

void
NTupleFCN::
setFitRange ( bool yes )
{
  m_fit_range = yes;
}

bool
NTupleFCN::
acceptRow ( unsigned int row ) const
{
    bool yes = true;
    if ( m_fit_cut != 0 &&
         m_fit_cut -> isEnabled () &&
         m_fit_range ) {
      yes = m_fit_cut -> acceptRow ( m_ntuple, row );
    }

    return yes;
}

---