Bins1DProfile.cxx

Go to the documentation of this file.

#ifdef HAVE_CONFIG_H
#include "config.h"
#else
#ifdef _MSC_VER
#include "msdevstudio/MSconfig.h"
#endif
#endif

#include "binners/Bins1DProfile.h"

#include "datasrcs/DataPointTuple.h"
#include "datasrcs/NTuple.h"

#include <algorithm>
#include <numeric>

#include <cmath>

#include <cassert>

using namespace hippodraw;

#ifdef ITERATOR_MEMBER_DEFECT
using namespace std;
#else
using std::fill;
using std::list;
using std::sqrt;
using std::string;
using std::vector;
using std::min_element;
using std::max_element;
#endif

Bins1DProfile::Bins1DProfile ( )
  : Bins1DBase ( "Bins1DProfile" )
{
}

Bins1DProfile::Bins1DProfile ( const Bins1DProfile & binner )
  : Bins1DBase ( binner ),
    m_sum( binner.m_sum ), 
    m_sumsq( binner.m_sumsq ),
    m_num( binner.m_num )
{
}

BinsBase *
Bins1DProfile::clone () const
{
  return new Bins1DProfile(*this);
}

Bins1DProfile::~Bins1DProfile ()
{
}

double Bins1DProfile::minBin()
{
  return *min_element ( m_sum.begin() + 1, m_sum.end ( ) - 1 );
}

double Bins1DProfile::maxBin()
{
  return *max_element ( m_sum.begin () + 1, m_sum.end () - 1 );
}

void Bins1DProfile::resize ( int number )
{
  m_sum.resize( number + 2 );
  m_num.resize( number + 2 );
  m_sumsq.resize( number + 2 );

  reset();

  m_values_dirty = true;
}

void Bins1DProfile::reset()
{
  fill( m_sum.begin(), m_sum.end(), 0.0 );
  fill( m_num.begin(), m_num.end(), 0 ); // This is an int array
  fill( m_sumsq.begin(), m_sumsq.end(), 0.0 );

  m_empty = true;
}

int Bins1DProfile::getNumberOfEntries () const
{
  return static_cast < int > ( std::accumulate ( m_num.begin() + 1, 
                                                 m_num.end() - 2, 0.0 ) );
}

int Bins1DProfile::getNumberOfEntries ( int i ) const
{
  return *( m_num.begin() + i + 1 );
}

void Bins1DProfile::accumulate( double x, double y, double, double )
{
  int i = binNumber ( x );

  m_sum[i] += y;
  m_sumsq[i] += y * y;
  m_num[i] ++;

  m_empty = false;
}


NTuple *
Bins1DProfile::
createNTuple () const
{
  unsigned int size = m_sum.size ();
  NTuple * ntuple = prepareNTuple ( size );

  fillProjectedValues ( ntuple );

  return ntuple;
}

namespace dp = hippodraw::DataPoint2DTuple;

void
Bins1DProfile::
fillProjectedValues ( DataSource * ntuple ) const
{
  ntuple -> clear();
  vector < double > row ( dp::SIZE );

//   double total = getNumberOfEntries ();
//   double factor = m_is_scaling ? m_scale_factor / total : 1.0;

  double width;
  double x = getLow ( Axes::X );
  int x_number = numberOfBins ( Axes::X );
  for ( int i = 0; i < x_number; i++ ) {
    width = binWidth ( i );
    x += 0.5 * width;
    double y = 0;
    double num = m_num [ i+1 ];
    if ( num <= 0 ) {
      x += 0.5 * width;
      continue;
    }
    if ( m_sum[i+1] != 0 ) {
      y = m_sum[i+1] / m_num[i+1];
    }
    else y = 0.;

    double xerr = width;
    double yerr = 0.0;

    if ( m_num[i+1] > 1. ) {
      double num = m_num[i+1];
//       yerr = sqrt ( ( m_sumsq[i+1] / ( num - 1. )  - y * y ) 
//                  / ( factor * num ) ) ;
// Restore profile error bars as the RMS of the distribution.
      yerr = sqrt ( ( m_sumsq[i+1] / ( num - 1. )  - y * y ) );
    }
    else {
      yerr = 0.0;
    }

    row[dp::X] = x;
    row[dp::Y] = y;
    row[dp::XERR] = xerr;
    row[dp::YERR] = yerr;

    ntuple -> addRow ( row );

    x += 0.5 * width;
  }
}

void
Bins1DProfile::
fillDataSource ( DataSource * ntuple ) const
{
  ntuple -> clear();
  vector < double > row ( dp::SIZE );

//   double total = getNumberOfEntries ();
//   double factor = m_is_scaling ? m_scale_factor / total : 1.0;

  double width;
  double x = getLow ( Axes::X );
  int x_number = numberOfBins ( Axes::X );
  for ( int i = 0; i < x_number; i++ ) {
    width = binWidth ( i );
    double half_width = 0.5 * width;
//     x += 0.5 * width;
    x += half_width;
    double y = 0;
    double num = m_num [ i+1 ];
    if ( num <= 0 ) {
//       x += 0.5 * width;
      x += half_width;
      continue;
    }
    if ( m_sum[i+1] != 0 ) {
      y = m_sum[i+1] / m_num[i+1];
    }
    else y = 0.;

//     double xerr = width;
    double yerr = 0.0;

    if ( m_num[i+1] > 1. ) {
      double num = m_num[i+1];
//       yerr = sqrt ( ( m_sumsq[i+1] / ( num - 1. )  - y * y ) 
//                  / ( factor * num ) ) ;
// Restore profile error bars as the RMS of the distribution.
      yerr = sqrt ( ( m_sumsq[i+1] / ( num - 1. )  - y * y ) );
    }
    else {
      yerr = 0.0;
    }

    row[dp::X] = x;
    row[dp::Y] = y;
//     row[dp::XERR] = xerr;
    row[dp::XERR] = half_width;
    row[dp::YERR] = yerr;

    ntuple -> addRow ( row );

//     x += 0.5 * width;
    x += half_width;
  }
}

void
Bins1DProfile::
setBinContents ( const DataSource * ntuple )
{
}

---