Bins1DHist.cxx
Go to the documentation of this file.
1 
12 #ifdef _MSC_VER
13 #include "msdevstudio/MSconfig.h"
14 #endif
15 
16 
17 #include "Bins1DHist.h"
18 
20 #include "datasrcs/NTuple.h"
21 
22 #include <algorithm>
23 #include <numeric>
24 
25 #include <cmath>
26 
27 #include <cassert>
28 
29 using namespace hippodraw;
30 
31 using std::fill;
32 using std::list;
33 using std::max_element;
34 using std::min_element;
35 using std::string;
36 using std::vector;
37 
39  : Bins1DBase ( "Bins1DHist" )
40 {
41  m_min_entries=0;
42 }
43 
45  : Bins1DBase ( binner ),
46  m_data( binner.m_data ),
47  m_variance( binner.m_variance ),
48  m_min_entries ( binner.m_min_entries )
49 {
50  for ( int i = 0; i < 3; i++ ) m_moments[i] = binner.m_moments[i];
51 }
52 
53 BinsBase *
55 {
56  return new Bins1DHist ( *this );
57 }
58 
60 {
61 }
62 
63 void Bins1DHist::resize ( int number )
64 {
65  m_data.resize ( number + 2 );
66  m_variance.resize ( number );
67  reset();
68 
69  m_values_dirty = true;
70 }
71 
73 {
74  return *min_element ( m_data.begin() + 1, m_data.end ( ) - 1 );
75 }
76 
78 {
79  return *max_element ( m_data.begin () + 1, m_data.end () - 1 );
80 }
81 
83 {
84  fill( m_data.begin(), m_data.end(), 0.0 );
85  fill( m_variance.begin(), m_variance.end(), 0.0 );
86  fill( m_moments, m_moments + 3, 0.0 );
87 
88  m_values_dirty = true;
89  m_empty = true;
90 }
91 
93 {
94  return static_cast < int > ( std::accumulate ( m_data.begin()+1,
95  m_data.end()-1, 0.0 ) );
96 }
97 
98 int Bins1DHist::getNumberOfEntries ( int i ) const
99 {
100  int retval = static_cast < int > ( m_data[i+1] );
101 
102  return retval;
103 }
104 
106 {
107  return static_cast < int > ( m_data[0] );
108 }
109 
111 {
112  return static_cast < int > ( *(m_data.end()-1) );
113 }
114 
115 void Bins1DHist::accumulate ( double x, double w, double, double )
116 {
117  int i = binNumber ( x );
118  if ( i > 0 && i <= static_cast < int > (m_variance.size() ) ) {
119  m_variance[i-1] += w * w;
120  m_moments[0] += w;
121  m_moments[1] += w * x;
122  m_moments[2] += w * x * x;
123  }
124  m_data[i] += w;
125 
126  m_empty = false;
127 }
128 
129 NTuple *
131 createNTuple ( ) const
132 {
133  unsigned int size = m_data.size ();
134  NTuple * ntuple = prepareNTuple ( size );
135 
136  fillDataSource ( ntuple );
137 
138  return ntuple;
139 }
140 
141 namespace dp = hippodraw::DataPoint2DTuple;
142 
143 void
145 fillDataSource ( DataSource * ntuple ) const
146 {
147  ntuple -> clear();
148  vector < double > row ( dp::SIZE );
149 
150  double total = std::accumulate ( m_data.begin() + 1,
151  m_data.end () - 1, 0.0 );
152  double factor = m_is_scaling ? m_scale_factor / total : 1.0;
153 
154  std::vector<double>::const_iterator dit = m_data.begin() + 1;
155  std::vector<double>::const_iterator vit = m_variance.begin();
156  std::vector<double>::const_iterator first_non_zero = m_data.begin() + 1;
157  std::vector<double>::const_iterator last_non_zero = m_data.end() - 2;
158 
159  if ( total > 0 ) {
160  // Get the first and last non_zero elements.
161  while ( (*first_non_zero) == 0 ) {
162  first_non_zero++;
163  }
164  while ( (*last_non_zero) == 0 ) {
165  last_non_zero--;
166  }
167  }
168  else {
169  swap ( first_non_zero, last_non_zero );
170  }
171 
172  double x = getLow ( Axes::X );
173 
174  int i = 0;
175 
176  if ( last_non_zero < first_non_zero ) { // all non zero
177  for (;dit !=m_data.end()-1; dit ++ ) {
178  double width = binWidth ( i++ );
179  double half_width = 0.5 * width;
180  double y = factor * ( *dit / width );
181  double yerr = factor * ( sqrt( *vit++ ) / width );
182  x += half_width;
183  row[dp::X] = x;
184  row[dp::Y] = y;
185  row[dp::XERR] = half_width;
186  row[dp::YERR] = yerr;
187 
188  ntuple -> addRow ( row );
189 
190  x += half_width;
191  }
192  return;
193  }
194 
195  // Process the leading zeros.
196  for ( ; dit != first_non_zero; dit++ ) {
197  double width = binWidth ( i++ );
198  double half_width = 0.5 * width;
199  double y = factor * ( *dit / width );
200  double yerr = factor * ( sqrt( *vit++ ) / width );
201  x += half_width;
202  row[dp::X] = x;
203  row[dp::Y] = y;
204  row[dp::XERR] = half_width;
205  row[dp::YERR] = yerr;
206 
207  ntuple -> addRow ( row );
208 
209  x += half_width;
210  }
211 
212 
213  // Redo the binning.
214  for( ; dit != last_non_zero+1; dit++ ) {
215  double width = binWidth( i++ );
216  double entries = *dit;
217 
218  double var = sqrt(*vit);
219  while ( entries < m_min_entries ) {
220  if (dit == last_non_zero) break; // need to change last row
221  dit++; vit++;
222  width += binWidth( i++ );
223  entries += *dit;
224  var += sqrt(*vit);
225  }
226 
227  double half_width = 0.5 * width;
228  x += half_width;
229  double y = factor * ( entries / width );
230  double yerr = factor * ( var / width );
231  vit++;
232 
233  // Update last bin if the remain bins has too few entries
234  if (( entries < m_min_entries ) && ( dit==last_non_zero )) {
235  unsigned int numOfRows = ntuple->rows();
236  unsigned int lastIndex = numOfRows-1;
237  vector <double> lastRow = ntuple -> getRow( lastIndex );
238 
239  x = lastRow[dp::X] + half_width;
240  y = ( lastRow[dp::Y] * lastRow[dp::XERR] * 2 + entries )
241  / ( width + lastRow[dp::XERR] *2 );
242  yerr = ( lastRow[dp::YERR] * lastRow[dp::XERR] * 2 + var )
243  / ( width + lastRow[dp::XERR] * 2 ) ;
244  half_width += lastRow[dp::XERR];
245 
246  ntuple -> eraseRow ( lastIndex );
247  }
248 
249  row[dp::X] = x;
250  row[dp::Y] = y;
251  row[dp::XERR] = half_width;
252  row[dp::YERR] = yerr;
253 
254  ntuple -> addRow ( row );
255 
256  x += half_width;
257  }
258 
259  // Processing the ending zeros.
260  for ( ; dit != m_data.end()-1; dit++ ) {
261  double width = binWidth ( i++ );
262  double half_width = 0.5 * width;
263  double y = factor * ( *dit / width );
264  double yerr = factor * ( sqrt( *vit++ ) / width );
265  x += half_width;
266  row[dp::X] = x;
267  row[dp::Y] = y;
268  row[dp::XERR] = half_width;
269  row[dp::YERR] = yerr;
270 
271  ntuple -> addRow ( row );
272 
273  x += half_width;
274  }
275 
276 }
277 
281 void
283 setBinContents ( const DataSource * ntuple )
284 {
285  unsigned int size = ntuple -> rows ();
286  assert ( size == m_variance.size () );
287 
288  for ( unsigned int i = 0; i < size; i++ ) {
289  const vector < double > & row = ntuple -> getRow ( i );
290  m_data[i+1] = row[dp::Y];
291  double yerr = row[dp::YERR];
292  double width = row[dp::XERR];
293  m_variance[i] = yerr * yerr * width * width;
294  }
295 }
296 
297 void
299 setMinEntries ( int entries )
300 {
301  m_min_entries = entries;
302 }
303 
304 int
307 {
308  return m_min_entries;
309 }

Generated for HippoDraw Class Library by doxygen