00001
00016 #ifdef _MSC_VER
00017
00018 #include "msdevstudio/MSconfig.h"
00019 #endif
00020
00021 #include "Bins2DHist.h"
00022
00023 #include "BinnerAxis.h"
00024
00025 #include "datasrcs/DataPointTuple.h"
00026 #include "datasrcs/NTuple.h"
00027
00028 #include <cmath>
00029
00030 #include <cassert>
00031
00032 #ifdef ITERATOR_MEMBER_DEFECT
00033 using namespace std;
00034 #else
00035 using std::fill;
00036 using std::list;
00037 using std::sqrt;
00038 using std::vector;
00039 #endif
00040
00041 using namespace hippodraw;
00042
00043 Bins2DHist::Bins2DHist ( )
00044 : Bins2DBase ( "Bins2DHist" )
00045 {
00046 }
00047
00048 Bins2DHist::Bins2DHist ( const Bins2DHist & binner )
00049 : Bins2DBase( binner ),
00050 m_variance( binner.m_variance )
00051 {
00052 m_num_bins = binner.m_num_bins;
00053
00054 for ( int i = 0; i < 3; i++ ) {
00055 m_x_moments[i] = binner.m_x_moments[i];
00056 m_y_moments[i] = binner.m_y_moments[i];
00057 }
00058 }
00059
00060 BinsBase *
00061 Bins2DHist::clone () const
00062 {
00063 return new Bins2DHist ( *this );
00064 }
00065
00066 Bins2DHist::~Bins2DHist ()
00067 {
00068 }
00069
00070 void
00071 Bins2DHist::
00072 setNumberOfBins ( hippodraw::Axes::Type axis, int nb )
00073 {
00074 assert ( axis == Axes::X || axis == Axes::Y );
00075
00076 if ( axis == Axes::X ) {
00077 m_data.resize ( nb + 2 );
00078 m_variance.resize ( nb );
00079 Bins2DBase::setNumberOfBins ( axis, nb );
00080 }
00081 else {
00082 Bins2DBase::setNumberOfBins ( axis, nb );
00083
00084 int number_y = numberOfBins ( Axes::Y );
00085
00086 unsigned int i = 0;
00087 for( ; i < m_data.size(); i++ ) {
00088 m_data[i].resize( number_y + 2 );
00089 }
00090
00091 for( i = 0; i < m_variance.size(); i++ ) {
00092 m_variance[i].resize( number_y );
00093 }
00094 }
00095
00096 reset ();
00097 }
00098
00099 void Bins2DHist::reset ()
00100 {
00101 unsigned int i = 0;
00102 for ( ; i < m_data.size(); i++ ) {
00103 fill ( m_data[i].begin(), m_data[i].end(), 0.0 );
00104 }
00105
00106 for ( i = 0; i < m_variance.size(); i++ ) {
00107 fill ( m_variance[i].begin(), m_variance[i].end(), 0.0 );
00108 }
00109
00110 fill ( m_x_moments, m_x_moments + 3, 0.0 );
00111 fill ( m_y_moments, m_y_moments + 3, 0.0 );
00112
00113 m_empty = true;
00114 }
00115
00116 void Bins2DHist::accumulate( double x, double y, double wt, double )
00117 {
00118 int i = binNumberX ( x );
00119 int j = binNumberY ( y );
00120
00121 if( i > 0 && i <= numberOfBins ( Axes::X ) &&
00122 j > 0 && j <= numberOfBins ( Axes::Y ) ) {
00123 m_variance[i-1][j-1] += wt * wt;
00124 m_x_moments[0] += wt;
00125 m_x_moments[1] += wt * x;
00126 m_x_moments[2] += wt * x * x;
00127 m_y_moments[0] += wt;
00128 m_y_moments[1] += wt * y;
00129 m_y_moments[2] += wt * y * y;
00130 }
00131 m_data[i][j] += wt;
00132
00133 m_empty = false;
00134 }
00135
00136 double Bins2DHist::getZValue ( double x, double y ) const
00137 {
00138 int i = binNumberX( x );
00139 int j = binNumberY( y );
00140 double widthX = binWidthX ( i-1 );
00141 double widthY = binWidthY ( j-1 );
00142 double v = m_data[i][j] / ( widthX * widthY );
00143 return v;
00144 }
00145
00146 NTuple *
00147 Bins2DHist::
00148 createNTuple () const
00149 {
00150 unsigned int size = numberOfBins ( Axes::X ) * numberOfBins ( Axes::Y );
00151 NTuple * ntuple = prepareNTuple ( size );
00152
00153 fillDataSource ( ntuple );
00154
00155 return ntuple;
00156 }
00157
00158 namespace dp = hippodraw::DataPoint3DTuple;
00159
00160 void
00161 Bins2DHist::
00162 fillDataSource ( DataSource * ntuple ) const
00163 {
00164 ntuple -> clear ();
00165
00166 double total = getNumberOfEntries ();
00167 double factor = m_is_scaling ? m_scale_factor / total : 1.0;
00168
00169 vector < double > row ( dp::SIZE );
00170
00171 double x = getLow ( Axes::X );
00172
00173 unsigned int num_x = numberOfBins ( Axes::X );
00174 unsigned int num_y = numberOfBins ( Axes::Y );
00175
00176 for ( unsigned int i = 0; i < num_x; i++ ) {
00177
00178 double widthX = binWidthX ( i );
00179 double half_widthX = 0.5 * widthX;
00180 x += half_widthX;
00181
00182 double y = getLow ( Axes::Y );
00183
00184 for ( unsigned int j = 0; j < num_y; j++ ) {
00185 double widthY = binWidthY ( j );
00186 double half_widthY = 0.5 * widthY;
00187 y += half_widthY;
00188
00189 double v = factor * ( m_data[i+1][j+1] / ( widthX * widthY ) );
00190 double verr = factor * ( sqrt ( m_variance[i][j] ) );
00191
00192 row[dp::X] = x;
00193 row[dp::Y] = y;
00194 row[dp::Z] = v;
00195 row [dp::XERR] = half_widthX;
00196 row [dp::YERR] = half_widthY;
00197 row [dp::ZERR] = verr;
00198
00199 ntuple -> addRow ( row );
00200
00201 y += half_widthY;
00202 }
00203 x += half_widthX;
00204 }
00205
00206 vector < unsigned int > shape ( 3 );
00207 shape[0] = num_x;
00208 shape[1] = num_y;
00209 shape[2] = dp::SIZE;
00210
00211 ntuple -> setShape ( shape );
00212 }
00213
00214 void
00215 Bins2DHist::
00216 setBinContents ( const DataSource * ntuple )
00217 {
00218 unsigned int num_x = numberOfBins ( Axes::X );
00219 unsigned int num_y = numberOfBins ( Axes::Y );
00220 unsigned int r = 0;
00221 for ( unsigned int i = 0; i < num_x; i++ ) {
00222 double widthX = binWidthX ( i );
00223 for ( unsigned int j = 0; j < num_y; j++ ) {
00224 double widthY = binWidthY ( j );
00225 const vector < double > & row = ntuple -> getRow ( r++ );
00226 double value = row [ dp::Z ];
00227 double verr = row [ dp::ZERR ];
00228 m_data[i+1][j+1] = value * ( widthX * widthY );
00229 m_variance[i][j] = verr * verr;
00230 }
00231 }
00232 }