Bins1DProfile.cxx
Go to the documentation of this file.
1 
12 #ifdef _MSC_VER
13 #include "msdevstudio/MSconfig.h"
14 #endif
15 
16 #include "binners/Bins1DProfile.h"
17 
19 #include "datasrcs/NTuple.h"
20 
21 #include <algorithm>
22 #include <numeric>
23 
24 #include <cmath>
25 
26 #include <cassert>
27 
28 
29 #ifdef ITERATOR_MEMBER_DEFECT
30 using namespace std;
31 #else
32 using std::fill;
33 using std::list;
34 using std::sqrt;
35 using std::string;
36 using std::vector;
37 using std::min_element;
38 using std::max_element;
39 #endif
40 
41 using namespace hippodraw;
42 
43 Bins1DProfile::Bins1DProfile ( )
44  : Bins1DBase ( "Bins1DProfile" )
45 {
46 }
47 
49  : Bins1DBase ( binner ),
50  m_sum( binner.m_sum ),
51  m_sumsq( binner.m_sumsq ),
52  m_num( binner.m_num ),
53  m_sumwtsq (binner.m_sumwtsq )
54 {
55 }
56 
57 BinsBase *
59 {
60  return new Bins1DProfile(*this);
61 }
62 
64 {
65 }
66 
68 {
69  return *min_element ( m_sum.begin() + 1, m_sum.end ( ) - 1 );
70 }
71 
73 {
74  return *max_element ( m_sum.begin () + 1, m_sum.end () - 1 );
75 }
76 
77 void Bins1DProfile::resize ( int number )
78 {
79  m_sum.resize( number + 2 );
80  m_num.resize( number + 2 );
81  m_sumsq.resize( number + 2 );
82  m_sumwtsq.resize( number + 2);
83 
84  reset();
85 
86  m_values_dirty = true;
87 }
88 
90 {
91  fill( m_sum.begin(), m_sum.end(), 0.0 );
92  fill( m_num.begin(), m_num.end(), 0 ); // This is an int array
93  fill( m_sumsq.begin(), m_sumsq.end(), 0.0 );
94  fill( m_sumwtsq.begin(), m_sumwtsq.end(), 0.0 );
95 
96  m_empty = true;
97 }
98 
100 {
101  return static_cast < int > ( std::accumulate ( m_num.begin() + 1,
102  m_num.end() - 2, 0.0 ) );
103 }
104 
106 {
107 // return *( m_num.begin() + i + 1 );
108  return static_cast < int > (m_num [ i+1] );
109 }
110 
112 {
113  return -1;
114 }
115 
117 {
118  return -1;
119 }
120 
121 
122 void Bins1DProfile::accumulate( double x, double y, double wt, double )
123 {
124  int i = binNumber ( x );
125 
126  m_sum[i] += y * wt; // sum(w*y)
127  m_sumsq[i] += y * y * wt; // sum(w*y*y)
128  m_num[i] += wt; // sum(w), actually m_sumwt
129 
130  m_sumwtsq[i] += wt * wt; // sum(w*w)
131 
132  m_empty = false;
133 }
134 
135 
136 NTuple *
138 createNTuple () const
139 {
140  unsigned int size = m_sum.size ();
141  NTuple * ntuple = prepareNTuple ( size );
142 
143  fillDataSource ( ntuple );
144 
145  return ntuple;
146 }
147 
148 namespace dp = hippodraw::DataPoint2DTuple;
149 
150 void
152 fillDataSource ( DataSource * ntuple ) const
153 {
154  ntuple -> clear();
155  vector < double > row ( dp::SIZE );
156 
157  double width;
158  double x = getLow ( Axes::X );
159  int x_number = numberOfBins ( Axes::X );
160  for ( int i = 0; i < x_number; i++ ) {
161  width = binWidth ( i );
162  double half_width = 0.5 * width;
163  x += half_width;
164  double y = 0;
165  double num = m_num [ i+1 ];
166  if ( num <= 0 ) {
167  x += half_width;
168  continue;
169  }
170  if ( m_sum[i+1] != 0 ) {
171  y = m_sum[i+1] / m_num[i+1];
172  }
173  else y = 0.;
174 
175  double yerr = 0.0;
176 
177  if ( m_num[i+1] > 1. ) {
178  double num = m_num[i+1];
179  /* Weighted function from
180  http://pygsl.sourceforge.net/reference/pygsl/node36.html
181 
182  Special case: when all wt==1;
183  yerr = sqrt ( m_sumsq[i+1] / (num - 1 ) - y * y * ( num/ (num-1) ) );
184 
185  When num >> 1, same as unweighed case. But what if num == 2?
186  */
187  yerr = sqrt ( (num/((num*num)-m_sumwtsq[i+1]))*
188  (m_sumsq[i+1]-2.0*y*m_sum[i+1]+y*y*num) );
189 
190  /* Unweight RMS function
191  yerr = sqrt ( ( m_sumsq[i+1] / ( num - 1. ) - y * y ) );
192  */
193  }
194  else {
195  yerr = 0.0;
196  }
197 
198  row[dp::X] = x;
199  row[dp::Y] = y;
200  row[dp::XERR] = half_width;
201  row[dp::YERR] = yerr;
202 
203  ntuple -> addRow ( row );
204 
205  x += half_width;
206  }
207 }
208 
212 void
214 setBinContents ( const DataSource * ) // ntuple )
215 {
216 }

Generated for HippoDraw Class Library by doxygen