00001
00002
00003
00004
00005
00006
00007
00008
00009 #ifndef LWH_Histogram1D_H
00010 #define LWH_Histogram1D_H
00011
00012
00013
00014
00015 #include "AIHistogram1D.h"
00016 #include "ManagedObject.h"
00017 #include "Axis.h"
00018 #include "VariAxis.h"
00019 #include <vector>
00020 #include <stdexcept>
00021
00022 namespace LWH {
00023
00024 using namespace AIDA;
00025
00029 class Histogram1D: public IHistogram1D, public ManagedObject {
00030
00031 public:
00032
00034 friend class HistogramFactory;
00035
00036 public:
00037
00041 Histogram1D(int n, double lo, double up)
00042 : fax(new Axis(n, lo, up)), vax(0),
00043 sum(n + 2), sumw(n + 2), sumw2(n + 2), sumxw(n + 2), sumx2w(n + 2) {
00044 ax = fax;
00045 }
00046
00050 Histogram1D(const std::vector<double> & edges)
00051 : fax(0), vax(new VariAxis(edges)),
00052 sum(edges.size() + 1), sumw(edges.size() + 1), sumw2(edges.size() + 1),
00053 sumxw(edges.size() + 1), sumx2w(edges.size() + 1) {
00054 ax = vax;
00055 }
00056
00060 Histogram1D(const Histogram1D & h)
00061 : IBaseHistogram(h), IHistogram(h), IHistogram1D(h), ManagedObject(h),
00062 fax(0), vax(0), sum(h.sum), sumw(h.sumw), sumw2(h.sumw2),
00063 sumxw(h.sumxw), sumx2w(h.sumx2w) {
00064 const VariAxis * hvax = dynamic_cast<const VariAxis *>(h.ax);
00065 if ( vax ) ax = vax = new VariAxis(*hvax);
00066 else ax = fax = new Axis(dynamic_cast<const Axis &>(*h.ax));
00067 }
00068
00070 virtual ~Histogram1D() {
00071 delete ax;
00072 }
00073
00078 std::string title() const {
00079 return theTitle;
00080 }
00081
00086 std::string name() const {
00087 return theTitle;
00088 }
00089
00095 bool setTitle(const std::string & title) {
00096 theTitle = title;
00097 return true;
00098 }
00099
00103 IAnnotation & annotation() {
00104 throw std::runtime_error("LWH cannot handle annotations");
00105 }
00106
00110 const IAnnotation & annotation() const {
00111 throw std::runtime_error("LWH cannot handle annotations");
00112 }
00113
00118 int dimension() const {
00119 return 1;
00120 }
00121
00126 bool reset() {
00127 sum = std::vector<int>(ax->bins() + 2);
00128 sumw = std::vector<double>(ax->bins() + 2);
00129 sumxw = std::vector<double>(ax->bins() + 2);
00130 sumx2w = std::vector<double>(ax->bins() + 2);
00131 sumw2 = std::vector<double>(ax->bins() + 2);
00132 return true;
00133 }
00134
00140 int entries() const {
00141 int si = 0;
00142 for ( int i = 2; i < ax->bins() + 2; ++i ) si += sum[i];
00143 return si;
00144 }
00145
00153 int allEntries() const {
00154 return entries() + extraEntries();
00155 }
00156
00161 int extraEntries() const {
00162 return sum[0] + sum[1];
00163 }
00164
00170 double equivalentBinEntries() const {
00171 double sw = 0.0;
00172 double sw2 = 0.0;
00173 for ( int i = 2; i < ax->bins() + 2; ++i ) {
00174 sw += sumw[i];
00175 sw2 += sumw2[i];
00176 }
00177 return sw2/(sw*sw);
00178 }
00179
00186 double sumBinHeights() const {
00187 double sw = 0.0;
00188 for ( int i = 2; i < ax->bins() + 2; ++i ) sw += sumw[i];
00189 return sw;
00190 }
00191
00197 double sumAllBinHeights() const {
00198 return sumBinHeights() + sumExtraBinHeights();
00199 }
00200
00205 double sumExtraBinHeights() const {
00206 return sumw[0] + sumw[1];
00207 }
00208
00214 double minBinHeight() const {
00215 double minw = sumw[2];
00216 for ( int i = 3; i < ax->bins() + 2; ++i ) minw = std::min(minw, sumw[i]);
00217 return minw;
00218 }
00219
00225 double maxBinHeight() const{
00226 double maxw = sumw[2];
00227 for ( int i = 3; i < ax->bins() + 2; ++i ) maxw = std::max(maxw, sumw[i]);
00228 return maxw;
00229 }
00230
00238 bool fill(double x, double weight = 1.) {
00239 int i = ax->coordToIndex(x) + 2;
00240 ++sum[i];
00241 sumw[i] += weight;
00242 sumxw[i] += x*weight;
00243 sumx2w[i] += x*x*weight;
00244 sumw2[i] += weight*weight;
00245 return weight >= 0 && weight <= 1;
00246 }
00247
00253 double binMean(int index) const {
00254 int i = index + 2;
00255 return sumw[i] != 0.0? sumxw[i]/sumw[i]:
00256 ( vax? vax->binMidPoint(index): fax->binMidPoint(index) );
00257 };
00258
00264 double binRms(int index) const {
00265 int i = index + 2;
00266 return sumw[i] == 0.0 || sum[i] < 2? ax->binWidth(index):
00267 std::sqrt(std::max(sumw[i]*sumx2w[i] - sumxw[i]*sumxw[i], 0.0))/sumw[i];
00268 };
00269
00276 int binEntries(int index) const {
00277 return sum[index + 2];
00278 }
00279
00286 double binHeight(int index) const {
00287 return sumw[index + 2];
00288 }
00289
00296 double binError(int index) const {
00297 return std::sqrt(sumw2[index + 2]);
00298 }
00299
00304 double mean() const {
00305 double s = 0.0;
00306 double sx = 0.0;
00307 for ( int i = 2; i < ax->bins() + 2; ++i ) {
00308 s += sumw[i];
00309 sx += sumxw[i];
00310 }
00311 return s != 0.0? sx/s: 0.0;
00312 }
00313
00318 double rms() const {
00319 double s = 0.0;
00320 double sx = 0.0;
00321 double sx2 = 0.0;
00322 for ( int i = 2; i < ax->bins() + 2; ++i ) {
00323 s += sumw[i];
00324 sx += sumxw[i];
00325 sx2 += sumx2w[i];
00326 }
00327 return s != 0.0? std::sqrt(std::max(s*sx2 - sx*sx, 0.0))/s:
00328 ax->upperEdge() - ax->lowerEdge();
00329 }
00330
00335 const IAxis & axis() const {
00336 return *ax;
00337 }
00338
00346 int coordToIndex(double coord) const {
00347 return ax->coordToIndex(coord);
00348 }
00349
00355 bool add(const Histogram1D & h) {
00356 if ( ax->upperEdge() != h.ax->upperEdge() ||
00357 ax->lowerEdge() != h.ax->lowerEdge() ||
00358 ax->bins() != h.ax->bins() ) return false;
00359 for ( int i = 0; i < ax->bins() + 2; ++i ) {
00360 sum[i] += h.sum[i];
00361 sumw[i] += h.sumw[i];
00362 sumxw[i] += h.sumxw[i];
00363 sumx2w[i] += h.sumx2w[i];
00364 sumw2[i] += h.sumw2[i];
00365 }
00366 return true;
00367 }
00368
00374 bool add(const IHistogram1D & hist) {
00375 return add(dynamic_cast<const Histogram1D &>(hist));
00376 }
00377
00382 bool scale(double s) {
00383 for ( int i = 0; i < ax->bins() + 2; ++i ) {
00384 sumw[i] *= s;
00385 sumxw[i] *= s;
00386 sumx2w[i] *= s;
00387 sumw2[i] *= s*s;
00388 }
00389 return true;
00390 }
00391
00399 void normalize(double intg) {
00400 double oldintg = sumAllBinHeights();
00401 if ( oldintg == 0.0 ) return;
00402 for ( int i = 0; i < ax->bins() + 2; ++i ) {
00403 double fac = intg/oldintg;
00404 if ( i >= 2 ) fac /= (ax->binUpperEdge(i - 2) - ax->binLowerEdge(i - 2));
00405 sumw[i] *= fac;
00406 sumxw[i] *= fac;
00407 sumx2w[i] *= fac;
00408 sumw2[i] *= fac*fac;
00409 }
00410 }
00411
00416 double integral() const {
00417 double intg = sumw[0] + sumw[1];
00418 for ( int i = 2; i < ax->bins() + 2; ++i )
00419 intg += sumw[i]*(ax->binUpperEdge(i - 2) - ax->binLowerEdge(i - 2));
00420 return intg;
00421 }
00422
00427 void * cast(const std::string &) const {
00428 return 0;
00429 }
00430
00434 bool writeXML(std::ostream & os, std::string path, std::string name) {
00435 os << " <histogram1d name=\"" << name
00436 << "\"\n title=\"" << title()
00437 << "\" path=\"" << path
00438 << "\">\n <axis max=\"" << ax->upperEdge()
00439 << "\" numberOfBins=\"" << ax->bins()
00440 << "\" min=\"" << ax->lowerEdge()
00441 << "\" direction=\"x\"";
00442 if ( vax ) {
00443 os << ">\n";
00444 for ( int i = 0, N = ax->bins() - 1; i < N; ++i )
00445 os << " <binBorder value=\"" << ax->binUpperEdge(i) << "\"/>\n";
00446 os << " </axis>\n";
00447 } else {
00448 os << "/>\n";
00449 }
00450 os << " <statistics entries=\"" << entries()
00451 << "\">\n <statistic mean=\"" << mean()
00452 << "\" direction=\"x\"\n rms=\"" << rms()
00453 << "\"/>\n </statistics>\n <data1d>\n";
00454 for ( int i = 0; i < ax->bins() + 2; ++i ) if ( sum[i] ) {
00455 os << " <bin1d binNum=\"";
00456 if ( i == 0 ) os << "UNDERFLOW";
00457 else if ( i == 1 ) os << "OVERFLOW";
00458 else os << i - 2;
00459 os << "\" entries=\"" << sum[i]
00460 << "\" height=\"" << sumw[i]
00461 << "\"\n error=\"" << std::sqrt(sumw2[i])
00462 << "\" error2=\"" << sumw2[i]
00463 << "\"\n weightedMean=\"" << binMean(i - 2)
00464 << "\" weightedRms=\"" << binRms(i - 2)
00465 << "\"/>\n";
00466 }
00467 os << " </data1d>\n </histogram1d>" << std::endl;
00468 return true;
00469 }
00470
00471
00476 bool writeFLAT(std::ostream & os, std::string path, std::string name) {
00477 os << "# " << path << "/" << name << " " << ax->lowerEdge()
00478 << " " << ax->bins() << " " << ax->upperEdge()
00479 << " \"" << title() << " \"" << std::endl;
00480 for ( int i = 2; i < ax->bins() + 2; ++i )
00481 os << binMean(i - 2) << " "
00482 << sumw[i] << " " << sqrt(sumw2[i]) << " " << sum[i] << std::endl;
00483 os << std::endl;
00484 return true;
00485 }
00486
00487 private:
00488
00490 std::string theTitle;
00491
00493 IAxis * ax;
00494
00496 Axis * fax;
00497
00499 VariAxis * vax;
00500
00502 std::vector<int> sum;
00503
00505 std::vector<double> sumw;
00506
00508 std::vector<double> sumw2;
00509
00511 std::vector<double> sumxw;
00512
00514 std::vector<double> sumx2w;
00515
00517 IAnnotation * anno;
00518
00519 };
00520
00521 }
00522
00523 #endif
