00001
00012 #ifdef _MSC_VER
00013 #include "msdevstudio/MSconfig.h"
00014 #endif
00015
00016
00017 #include "Bins1DHist.h"
00018
00019 #include "datasrcs/DataPointTuple.h"
00020 #include "datasrcs/NTuple.h"
00021
00022 #include <algorithm>
00023 #include <numeric>
00024
00025 #include <cmath>
00026
00027 #include <cassert>
00028
00029 using namespace hippodraw;
00030
00031 using std::fill;
00032 using std::list;
00033 using std::max_element;
00034 using std::min_element;
00035 using std::string;
00036 using std::vector;
00037
00038 Bins1DHist::Bins1DHist ( )
00039 : Bins1DBase ( "Bins1DHist" )
00040 {
00041 m_min_entries=0;
00042 }
00043
00044 Bins1DHist::Bins1DHist ( const Bins1DHist & binner )
00045 : Bins1DBase ( binner ),
00046 m_data( binner.m_data ),
00047 m_variance( binner.m_variance ),
00048 m_min_entries ( binner.m_min_entries )
00049 {
00050 for ( int i = 0; i < 3; i++ ) m_moments[i] = binner.m_moments[i];
00051 }
00052
00053 BinsBase *
00054 Bins1DHist::clone () const
00055 {
00056 return new Bins1DHist ( *this );
00057 }
00058
00059 Bins1DHist::~Bins1DHist ()
00060 {
00061 }
00062
00063 void Bins1DHist::resize ( int number )
00064 {
00065 m_data.resize ( number + 2 );
00066 m_variance.resize ( number );
00067 reset();
00068
00069 m_values_dirty = true;
00070 }
00071
00072 double Bins1DHist::minBin()
00073 {
00074 return *min_element ( m_data.begin() + 1, m_data.end ( ) - 1 );
00075 }
00076
00077 double Bins1DHist::maxBin()
00078 {
00079 return *max_element ( m_data.begin () + 1, m_data.end () - 1 );
00080 }
00081
00082 void Bins1DHist::reset()
00083 {
00084 fill( m_data.begin(), m_data.end(), 0.0 );
00085 fill( m_variance.begin(), m_variance.end(), 0.0 );
00086 fill( m_moments, m_moments + 3, 0.0 );
00087
00088 m_values_dirty = true;
00089 m_empty = true;
00090 }
00091
00092 int Bins1DHist::getNumberOfEntries () const
00093 {
00094 return static_cast < int > ( std::accumulate ( m_data.begin()+1,
00095 m_data.end()-1, 0.0 ) );
00096 }
00097
00098 int Bins1DHist::getNumberOfEntries ( int i ) const
00099 {
00100 int retval = static_cast < int > ( m_data[i+1] );
00101
00102 return retval;
00103 }
00104
00105 int Bins1DHist::getUnderflow () const
00106 {
00107 return static_cast < int > ( m_data[0] );
00108 }
00109
00110 int Bins1DHist::getOverflow () const
00111 {
00112 return static_cast < int > ( *(m_data.end()-1) );
00113 }
00114
00115 void Bins1DHist::accumulate ( double x, double w, double, double )
00116 {
00117 int i = binNumber ( x );
00118 if ( i > 0 && i <= static_cast < int > (m_variance.size() ) ) {
00119 m_variance[i-1] += w * w;
00120 m_moments[0] += w;
00121 m_moments[1] += w * x;
00122 m_moments[2] += w * x * x;
00123 }
00124 m_data[i] += w;
00125
00126 m_empty = false;
00127 }
00128
00129 NTuple *
00130 Bins1DHist::
00131 createNTuple ( ) const
00132 {
00133 unsigned int size = m_data.size ();
00134 NTuple * ntuple = prepareNTuple ( size );
00135
00136 fillDataSource ( ntuple );
00137
00138 return ntuple;
00139 }
00140
00141 namespace dp = hippodraw::DataPoint2DTuple;
00142
00143 void
00144 Bins1DHist::
00145 fillDataSource ( DataSource * ntuple ) const
00146 {
00147 ntuple -> clear();
00148 vector < double > row ( dp::SIZE );
00149
00150 double total = std::accumulate ( m_data.begin() + 1,
00151 m_data.end () - 1, 0.0 );
00152 double factor = m_is_scaling ? m_scale_factor / total : 1.0;
00153
00154 std::vector<double>::const_iterator dit = m_data.begin() + 1;
00155 std::vector<double>::const_iterator vit = m_variance.begin();
00156 std::vector<double>::const_iterator first_non_zero = m_data.begin() + 1;
00157 std::vector<double>::const_iterator last_non_zero = m_data.end() - 2;
00158
00159 if ( total > 0 ) {
00160
00161 while ( (*first_non_zero) == 0 ) {
00162 first_non_zero++;
00163 }
00164 while ( (*last_non_zero) == 0 ) {
00165 last_non_zero--;
00166 }
00167 }
00168 else {
00169 swap ( first_non_zero, last_non_zero );
00170 }
00171
00172 double x = getLow ( Axes::X );
00173
00174 int i = 0;
00175
00176 if ( last_non_zero < first_non_zero ) {
00177 for (;dit !=m_data.end()-1; dit ++ ) {
00178 double width = binWidth ( i++ );
00179 double half_width = 0.5 * width;
00180 double y = factor * ( *dit / width );
00181 double yerr = factor * ( sqrt( *vit++ ) / width );
00182 x += half_width;
00183 row[dp::X] = x;
00184 row[dp::Y] = y;
00185 row[dp::XERR] = half_width;
00186 row[dp::YERR] = yerr;
00187
00188 ntuple -> addRow ( row );
00189
00190 x += half_width;
00191 }
00192 return;
00193 }
00194
00195
00196 for ( ; dit != first_non_zero; dit++ ) {
00197 double width = binWidth ( i++ );
00198 double half_width = 0.5 * width;
00199 double y = factor * ( *dit / width );
00200 double yerr = factor * ( sqrt( *vit++ ) / width );
00201 x += half_width;
00202 row[dp::X] = x;
00203 row[dp::Y] = y;
00204 row[dp::XERR] = half_width;
00205 row[dp::YERR] = yerr;
00206
00207 ntuple -> addRow ( row );
00208
00209 x += half_width;
00210 }
00211
00212
00213
00214 for( ; dit != last_non_zero+1; dit++ ) {
00215 double width = binWidth( i++ );
00216 double entries = *dit;
00217
00218 double var = sqrt(*vit);
00219 while ( entries < m_min_entries ) {
00220 if (dit == last_non_zero) break;
00221 dit++; vit++;
00222 width += binWidth( i++ );
00223 entries += *dit;
00224 var += sqrt(*vit);
00225 }
00226
00227 double half_width = 0.5 * width;
00228 x += half_width;
00229 double y = factor * ( entries / width );
00230 double yerr = factor * ( var / width );
00231 vit++;
00232
00233
00234 if (( entries < m_min_entries ) && ( dit==last_non_zero )) {
00235 unsigned int numOfRows = ntuple->rows();
00236 unsigned int lastIndex = numOfRows-1;
00237 vector <double> lastRow = ntuple -> getRow( lastIndex );
00238
00239 x = lastRow[dp::X] + half_width;
00240 y = ( lastRow[dp::Y] * lastRow[dp::XERR] * 2 + entries )
00241 / ( width + lastRow[dp::XERR] *2 );
00242 yerr = ( lastRow[dp::YERR] * lastRow[dp::XERR] * 2 + var )
00243 / ( width + lastRow[dp::XERR] * 2 ) ;
00244 half_width += lastRow[dp::XERR];
00245
00246 ntuple -> eraseRow ( lastIndex );
00247 }
00248
00249 row[dp::X] = x;
00250 row[dp::Y] = y;
00251 row[dp::XERR] = half_width;
00252 row[dp::YERR] = yerr;
00253
00254 ntuple -> addRow ( row );
00255
00256 x += half_width;
00257 }
00258
00259
00260 for ( ; dit != m_data.end()-1; dit++ ) {
00261 double width = binWidth ( i++ );
00262 double half_width = 0.5 * width;
00263 double y = factor * ( *dit / width );
00264 double yerr = factor * ( sqrt( *vit++ ) / width );
00265 x += half_width;
00266 row[dp::X] = x;
00267 row[dp::Y] = y;
00268 row[dp::XERR] = half_width;
00269 row[dp::YERR] = yerr;
00270
00271 ntuple -> addRow ( row );
00272
00273 x += half_width;
00274 }
00275
00276 }
00277
00281 void
00282 Bins1DHist::
00283 setBinContents ( const DataSource * ntuple )
00284 {
00285 unsigned int size = ntuple -> rows ();
00286 assert ( size == m_variance.size () );
00287
00288 for ( unsigned int i = 0; i < size; i++ ) {
00289 const vector < double > & row = ntuple -> getRow ( i );
00290 m_data[i+1] = row[dp::Y];
00291 double yerr = row[dp::YERR];
00292 double width = row[dp::XERR];
00293 m_variance[i] = yerr * yerr * width * width;
00294 }
00295 }
00296
00297 void
00298 Bins1DHist::
00299 setMinEntries ( int entries )
00300 {
00301 m_min_entries = entries;
00302 }
00303
00304 int
00305 Bins1DHist::
00306 getMinEntries ()
00307 {
00308 return m_min_entries;
00309 }