Page MenuHomeHEPForge
Files F8725598

No OneTemporary
Actions

View Options

diff --git a/Utilities/Histogram.cc b/Utilities/Histogram.cc
--- a/Utilities/Histogram.cc
+++ b/Utilities/Histogram.cc
@@ -1,362 +1,362 @@
// -*- C++ -*-
//
// Histogram.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the Histogram class.
//
#include "Histogram.h"
#include "HerwigStrategy.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "ThePEG/Handlers/EventHandler.h"
using namespace Herwig;
NoPIOClassDescription<Histogram> Histogram::initHistogram;
// Definition of the static class description member.
void Histogram::Init() {
static ClassDocumentation<Histogram> documentation
("The Histogram class implements a simple histogram include data"
" points for comparision with experimental results.");
}
void Histogram::topdrawOutput(ostream & out,
unsigned int flags,
string colour,
string title, string titlecase,
string left, string leftcase,
string bottom, string bottomcase
) const {
using namespace HistogramOptions;
bool frame = ( flags & Frame ) == Frame;
bool errorbars = ( flags & Errorbars ) == Errorbars;
bool xlog = ( flags & Xlog ) == Xlog;
bool ylog = ( flags & Ylog ) == Ylog;
bool smooth = ( flags & Smooth ) == Smooth;
bool rawcount = ( flags & Rawcount ) == Rawcount;
// output the title info if needed
if(frame) {
out << "NEW FRAME\n";
if(_havedata) out << "SET WINDOW X 1.6 8 Y 3.5 9\n";
else out << "SET WINDOW X 1.6 8 Y 1.6 9\n";
out << "SET FONT DUPLEX\n";
out << "TITLE TOP \"" << title << "\"\n";
out << "CASE \"" << titlecase << "\"\n";
out << "TITLE LEFT \"" << left << "\"\n";
out << "CASE \"" << leftcase << "\"\n";
out << (errorbars ? "SET ORDER X Y DX DY \n" : "SET ORDER X Y DX\n");
if (HerwigStrategy::version != "") {
out << "TITLE RIGHT \"" << HerwigStrategy::version << "\"\n";
out << "CASE \"\"\n";
}
if(_havedata) out << "SET AXIS BOTTOM OFF\n";
else {
out << "TITLE BOTTOM \"" << bottom << "\"\n";
out << "CASE \"" << bottomcase << "\"\n";
}
}
// scales
if(xlog && frame) out << "SET SCALE X LOG " << endl;
if(ylog && frame) out << "SET SCALE Y LOG " << endl;
// set the x limits
const unsigned int lastDataBinIndx = _bins.size()-2;
if (frame) {
out << "SET LIMITS X " << _bins[1].limit << " "
<< _bins[lastDataBinIndx+1].limit << endl;
}
// work out the y points
vector<double> yout;
double ymax=-9.8765e34,ymin=9.8765e34;
- unsigned int numPoints = _globalStats.numberOfPoints();
+ double numPoints = _total;
if (numPoints == 0) ++numPoints;
for(unsigned int ix=1; ix<=lastDataBinIndx; ++ix) {
double delta = 0.5*(_bins[ix+1].limit-_bins[ix].limit);
double factor = rawcount ? _prefactor : 0.5 * _prefactor / (numPoints * delta);
double value = factor*_bins[ix].contents;
yout.push_back(value);
ymax=max(ymax, max(value, _bins[ix].data+_bins[ix].dataerror) );
if(yout.back()>0.) ymin=min(ymin,value);
if(_bins[ix].data>0) ymin=min(ymin,_bins[ix].data);
}
if (ymin > 1e34) ymin = 1e-34;
if (ymax < 1e-33) ymax = 1e-33;
if (ymax < 10*ymin) ymin = 0.1*ymax;
if (ylog && frame) {
out << "SET LIMITS Y " << ymin << " " << ymax << endl;
}
// the histogram from the event generator
for(unsigned int ix=1; ix<=lastDataBinIndx; ++ix) {
double delta = 0.5*(_bins[ix+1].limit-_bins[ix].limit);
double factor = rawcount ? _prefactor : 0.5 * _prefactor / (numPoints * delta);
out << _bins[ix].limit+delta << '\t' << yout[ix-1] << '\t' << delta;
if (errorbars) {
out << '\t' << factor*sqrt(_bins[ix].contentsSq);
}
out << '\n';
}
// N.B. in td smoothing only works for histograms with uniform binning.
if(!smooth) {
out << "HIST " << colour << endl;
} else {
out << "SMOOTH Y LEVEL 2 " << endl;
out << "JOIN " << colour << endl;
}
if (_havedata) {
// the real experimental data
for(unsigned int ix=1; ix<=lastDataBinIndx; ++ix) {
double delta = 0.5*(_bins[ix+1].limit-_bins[ix].limit);
out << _bins[ix].limit+delta << '\t' << _bins[ix].data << '\t' << delta;
if (errorbars) out << '\t' << _bins[ix].dataerror;
out << '\n';
}
out << "PLOT " << endl;
out << "SET WINDOW X 1.6 8 Y 2.5 3.5\n";
out << "SET LIMITS X " << _bins[1].limit << " "
<< _bins[lastDataBinIndx+1].limit << "\n";
double ymax=0.;
out << _bins[1].limit << "\t" << _bins[1].dataerror/_bins[1].data << "\n";
for(unsigned int ix=1; ix<=lastDataBinIndx; ++ix) {
double delta = 0.5*(_bins[ix+1].limit-_bins[ix].limit);
if(_bins[ix].data!=0.) {
if(_bins[ix].dataerror/_bins[ix].data>ymax)
ymax=_bins[ix].dataerror/_bins[ix].data;
out << _bins[ix].limit+delta << '\t'
<< _bins[ix].dataerror/_bins[ix].data << '\n';
}
else {
out << _bins[ix].limit+delta << '\t'
<< 1. << '\n';
}
}
if(_bins[lastDataBinIndx].data!=0.) {
out << _bins[lastDataBinIndx+1].limit << "\t"
<< _bins[lastDataBinIndx].dataerror/_bins[lastDataBinIndx].data << "\n";
out << _bins[lastDataBinIndx+1].limit << "\t"
<<-_bins[lastDataBinIndx].dataerror/_bins[lastDataBinIndx].data << "\n";
}
else {
out << _bins[lastDataBinIndx+1].limit << "\t" << 1. << "\n";
out << _bins[lastDataBinIndx+1].limit << "\t" << -1. << "\n";
}
for(unsigned int ix=lastDataBinIndx;ix>=1;--ix) {
double delta = 0.5*(_bins[ix+1].limit-_bins[ix].limit);
if(_bins[ix].data!=0.) {
out << _bins[ix].limit+delta << '\t'
<< -_bins[ix].dataerror/_bins[ix].data << '\n';
}
else {
out << _bins[ix].limit+delta << '\t'
<< -1. << '\n';
}
}
if(_bins[1].data!=0.) {
out << _bins[1].limit << "\t" << -_bins[1].dataerror/_bins[1].data << "\n";
}
else {
out << _bins[1].limit << "\t" << -1. << "\n";
}
out << "set scale y lin\n";
out << "set limits y " << -ymax << " " << ymax << "\n";
out << "set fill full\n";
out << "join yellow fill yellow\n";
for(unsigned int ix=1; ix<=lastDataBinIndx; ++ix) {
double delta = 0.5*(_bins[ix+1].limit-_bins[ix].limit);
if(_bins[ix].data!=0.) {
out << _bins[ix].limit+delta << "\t"
<< (yout[ix-1]-_bins[ix].data)/_bins[ix].data << "\n";
}
else if(_bins[ix].dataerror!=0.) {
out << _bins[ix].limit+delta << "\t"
<< (yout[ix-1]-_bins[ix].data)/_bins[ix].dataerror << "\n";
}
else {
out << _bins[ix].limit+delta << "\t" << 0. << "\n";
}
}
out << "join\n";
out << "SET WINDOW X 1.6 8 Y 1.6 2.5\n";
out << "SET LIMITS X " << _bins[1].limit << " "
<< _bins[lastDataBinIndx+1].limit << "\n";
out << "SET AXIS BOTTOM ON\n";
out << "TITLE BOTTOM \"" << bottom << "\"\n";
out << "CASE \"" << bottomcase << "\"\n";
ymax =0.;
double ymin=0.;
for(unsigned int ix=1; ix<=lastDataBinIndx; ++ix) {
double delta = 0.5*(_bins[ix+1].limit-_bins[ix].limit);
double error = sqrt(sqr(0.5*sqrt(_bins[ix].contentsSq)/(delta*numPoints))+
sqr(_bins[ix].dataerror));
double point=(yout[ix-1]-_bins[ix].data)/error;
if(point<ymin) ymin=point;
if(point>ymax) ymax=point;
out << _bins[ix].limit+delta << '\t'
<< point << '\n';
}
out << "set limits y " << ymin << " " << ymax << "\n";
out << "JOIN" << endl;
}
}
double Histogram::dataNorm() const {
double norm(0.0);
if (_havedata) {
const unsigned int lastDataBinIndx = _bins.size()-2;
for(unsigned int ix=1; ix<=lastDataBinIndx; ++ix) {
double delta = _bins[ix+1].limit-_bins[ix].limit;
double value = _bins[ix].data;
norm += delta*value;
}
} else {
norm = -1.0;
}
return norm;
}
unsigned int Histogram::visibleEntries() const {
unsigned int numPoints(0);
const unsigned int lastDataBinIndx = _bins.size()-2;
for(unsigned int ix=1; ix<=lastDataBinIndx; ++ix) {
numPoints += static_cast<unsigned int>( _bins[ix].contents );
}
return numPoints;
}
void Histogram::simpleOutput(ostream & out, bool errorbars,
bool normdata) {
// simple ascii output (eg for gnuplot)
// work out the y points
vector<double> yout;
// unsigned int numPoints = _globalStats.numberOfPoints();
unsigned int numPoints = visibleEntries();
if (numPoints == 0) ++numPoints;
double datanorm(1.0);
double chisq(0.0), minfrac(0.05);
unsigned int ndof(0);
if (_havedata) {
if (normdata) datanorm = dataNorm();
normaliseToData();
chiSquared(chisq, ndof, minfrac);
}
prefactor(1.0);
const unsigned int lastDataBinIndx = _bins.size()-2;
for(unsigned int ix=1; ix<=lastDataBinIndx; ++ix) {
double delta = 0.5*(_bins[ix+1].limit-_bins[ix].limit);
double value = 0.5*_prefactor*_bins[ix].contents / (delta*numPoints);
yout.push_back(value);
}
out << "# " << numPoints << " entries, mean +- sigma = "
<< _globalStats.mean() << " +- "
<< _globalStats.stdDev() << "\n";
if (_havedata) {
out << "# chi^2/dof = " << chisq << "/" << ndof << " = "
<< chisq/double(ndof) << " (min err = " << minfrac << ")\n";
if (datanorm) {
out << "# data normalised by factor " << datanorm << "\n";
}
}
out << "# xlo xhi ynorm "
<< (errorbars ? "ynorm_err " : "")
<< (_havedata ? "data " : "")
<< (_havedata && errorbars ? "dataerr " : "")
<< "y_entr\n";
// the histogram from the event generator
for(unsigned int ix=1; ix<=lastDataBinIndx; ++ix) {
double delta = 0.5*(_bins[ix+1].limit-_bins[ix].limit);
out << _bins[ix].limit << " "
<< _bins[ix+1].limit << " "
<< yout[ix-1];
if (errorbars) {
out << " " << 0.5*sqrt(_bins[ix].contentsSq)/(delta*numPoints);
}
if (_havedata) {
out << " " << _bins[ix].data/datanorm;
if (errorbars)
out << " " << _bins[ix].dataerror/datanorm;
}
out << " " << _bins[ix].contents << '\n';
}
}
vector<double> Histogram::dumpBins() const {
vector<double> bincontents(_bins.size());
for (size_t i=0; i < _bins.size(); ++i)
bincontents[i] = _bins[i].contents;
return bincontents;
}
Histogram Histogram::ratioWith(const Histogram & h2) const {
const size_t numBins = _bins.size();
assert( numBins > 2 && numBins == h2._bins.size());
Histogram ratio(*this);
for (size_t i=0; i < numBins; ++i) {
assert(_bins[i].limit == h2._bins[i].limit);
if (h2._bins[i].contents > 0.0)
ratio._bins[i].contents /= h2._bins[i].contents;
else
ratio._bins[i].contents = 0.0;
}
return ratio;
}
void Histogram::normaliseToData()
{
double numer(0.),denom(0.);
unsigned int numPoints = _globalStats.numberOfPoints();
for(unsigned int ix=1;ix<_bins.size()-1;++ix)
{
double delta = 0.5*(_bins[ix+1].limit-_bins[ix].limit);
double value = 0.5*_bins[ix].contents / (delta*numPoints);
if(_bins[ix].dataerror>0.)
{
double var = sqr(_bins[ix].dataerror);
numer += _bins[ix].data * value/var;
denom += sqr(value)/var;
}
}
_prefactor = numer/denom;
}
void Histogram::chiSquared(double & chisq,
unsigned int & ndegrees, double minfrac) const {
chisq =0.;
ndegrees=0;
unsigned int numPoints = _globalStats.numberOfPoints();
for(unsigned int ix=1;ix<_bins.size()-1;++ix) {
double delta = 0.5*(_bins[ix+1].limit-_bins[ix].limit);
double value = 0.5*_prefactor*_bins[ix].contents / (delta*numPoints);
double error = _bins[ix].dataerror;
if(error>0.) {
if(abs(error/_bins[ix].data) < minfrac) error = minfrac*_bins[ix].data;
double var=sqr(error)
+ _bins[ix].contentsSq * sqr(0.5*_prefactor / (delta*numPoints));
chisq += sqr(_bins[ix].data - value) / var;
++ndegrees;
}
}
}
void Histogram::normaliseToCrossSection() {
unsigned int numPoints = _globalStats.numberOfPoints();
if (numPoints == 0) ++numPoints;
_prefactor=CurrentGenerator::current().eventHandler()->histogramScale()*
numPoints/nanobarn;
}
diff --git a/Utilities/Histogram.h b/Utilities/Histogram.h
--- a/Utilities/Histogram.h
+++ b/Utilities/Histogram.h
@@ -1,349 +1,354 @@
// -*- C++ -*-
//
// Histogram.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_Histogram_H
#define HERWIG_Histogram_H
//
// This is the declaration of the Histogram class.
//
#include "Histogram.fh"
#include "ThePEG/Interface/Interfaced.h"
#include "Statistic.h"
#include <string>
// workaround for OS X bug where isnan() and isinf() are hidden
// when <iostream> is included
extern "C" int isnan(double) throw();
namespace Herwig {
using namespace ThePEG;
/**
* Options for histogram output.
* They can be combined using the '|' operator, e.g. 'Frame | Ylog'
*/
namespace HistogramOptions {
const unsigned int None = 0; /**< No options */
const unsigned int Frame = 1; /**< Plot on new frame */
const unsigned int Errorbars = 1 << 1; /**< Plot error bars */
const unsigned int Xlog = 1 << 2; /**< log scale for x-axis */
const unsigned int Ylog = 1 << 3; /**< log scale for y-axis */
const unsigned int Smooth = 1 << 4; /**< smooth the line */
const unsigned int Rawcount = 1 << 5; /**< don't normalize to unit area */
}
/**
* The Histogram class is a simple histogram for the Analysis handlers.
*
* @see \ref HistogramInterfaces "The interfaces"
* defined for Histogram.
*/
class Histogram: public Interfaced {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
* @param lower The lower limit of the histogram
* @param upper The upper limit of the histogram
* @param nbin Number of bins
*/
inline Histogram(double lower=0., double upper=0., unsigned int nbin=0);
/**
* Constructor for variable width bins
* @param limits The lower limits for the bins followed by the upper limit of the last bin
*/
inline Histogram(vector<double> limits);
/**
* Constructor with data included
* @param limits The lower limits for the bins followed by the upper limit of the last bin
* @param data The data
* @param dataerror The errors on the data
*/
inline Histogram(vector<double> limits, vector<double> data, vector<double> dataerror);
//@}
public:
/**
* Operator to add a point to the histogrma with unit weight
*/
inline void operator+=(double);
/**
* Function to add a weighted point to the histogram
*/
inline void addWeighted(double data, double weight);
/**
* Number of bins (not counting the overflow)
*/
inline unsigned int numberOfBins() const;
/**
* Get the prefactor
*/
inline double prefactor() const;
/**
* Set the prefactor
*/
inline void prefactor(double );
/**
* Access to the statistics on the total entry of the histogram
*/
inline const Statistic & globalStatistics() const;
/**
* Normalise the distributions to the data
*/
void normaliseToData();
/**
* Normalise the distributions to the total cross section
*/
void normaliseToCrossSection();
/**
* Return the chi squared
* @param chisq The chi squared
* @param ndegrees The number of points
* @param minfrac The minimum fractional error on the data point
*/
void chiSquared(double & chisq,
unsigned int & ndegrees, double minfrac=0.) const;
/**
* Output as a topdrawer file. The histogram is normalised to unit area
* @param out The output stream
* @param flags A bitmask of flags from HistogramOptions, e.g. Frame|Ylog
* @param colour The colour for the line
* @param title The title for the top of the plot
* @param titlecase topdraw format for the title
* @param left Left axis lable
* @param leftcase topdraw format for left axis label
* @param bottom Bottom axis lable
* @param bottomcase Bottom axis lable ofr topdraw
* N.B. in td smoothing only works for histograms with uniform binning.
*/
void topdrawOutput(ostream & out,
unsigned int flags = 0,
string colour = string("BLACK"),
string title = string(),
string titlecase = string(),
string left = string(),
string leftcase = string(),
string bottom = string(),
string bottomcase = string()
) const;
/**
* get the number of visible entries (all entries without those in the
* under- and overflow bins) in the histogram. This assumes integer
* entries, ie it gives wrong results for weighted histograms.
*/
unsigned int visibleEntries() const;
/**
* Compute the normalisation of the data.
*/
double dataNorm() const;
/**
* Output into a simple ascii file, easily readable by gnuplot.
*/
void simpleOutput(ostream & out, bool errorbars, bool normdata=false);
/**
* Dump bin data into a vector
*/
vector<double> dumpBins() const;
/**
* Returns a new histogram containing bin-by-bin ratios of two histograms
*/
Histogram ratioWith(const Histogram & h2) const;
public:
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
inline virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
inline virtual IBPtr fullclone() const;
//@}
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static NoPIOClassDescription<Histogram> initHistogram;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
Histogram & operator=(const Histogram &);
private:
/**
* Global statistics of all data that went into the histogram.
*/
Statistic _globalStats;
/**
* Set to true if there is experimental data available
*/
bool _havedata;
/**
* One bin of the histogram. limit is the _lower_ bound of the bin.
*/
struct Bin {
/**
* Default constructor
*/
Bin() : contents(0.0), contentsSq(0.0),
limit(0.0), data(0.0), dataerror(0.0) {}
/**
* Contents of the bin
*/
double contents;
/**
* Contents squared for the error
*/
double contentsSq;
/**
* The limit for the bin
*/
double limit;
/**
* The experimental value for the bin
*/
double data;
/**
* The error on the experimental value for the bin
*/
double dataerror;
};
/**
* The histogram bins. _bins[0] is the underflow, _bins.back() the overflow
*/
vector<Bin> _bins;
/**
* Prefactors to multiply the output by
*/
double _prefactor;
+
+ /**
+ * Total entry
+ */
+ double _total;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of Histogram. */
template <>
struct BaseClassTrait<Herwig::Histogram,1> {
/** Typedef of the first base class of Histogram. */
typedef Herwig::Statistic NthBase;
};
/** This template specialization informs ThePEG about the name of
* the Histogram class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::Histogram>
: public ClassTraitsBase<Herwig::Histogram> {
/** Return a platform-independent class name */
static string className() { return "Herwig::Histogram"; }
};
/** @endcond */
}
#include "Histogram.icc"
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
// #include "Histogram.tcc"
#endif
// void SampleHistogram::printMoments(char* name, double Nmax, double dN,
// double x0, double x1) {
// ofstream out(name);
// if (!out) {
// cerr << "SampleHistoGram::printMoments: ERROR! Can't open file" << endl;
// }
// time_t now_t;
// now_t = time(0);
// out << "# created " << ctime(&now_t)
// << "# by SampleHistogram::printMoments(..., "
// << Nmax << ", " << dN << ")" << endl
// << "# " << this->samples() << " entries, mean +/- sigma = "
// << this->mean() << " +/- " << this->stdDev() << endl;
// double x0N, x1N, delta, hi;
// for (double N=dN; N < Nmax; N += dN) {
// double fN = 0.0;
// for(int i = 0; i < howManyBuckets-1; i++) {
// x0N = pow(bucketLimit[i], N);
// x1N = pow(bucketLimit[i+1], N);
// delta = (bucketLimit[i+1] - bucketLimit[i]);
// if (delta > 0 && this->samples() > 0
// && bucketLimit[i] >= x0 && bucketLimit[i] <= x1
// && bucketLimit[i+1] >= x0 && bucketLimit[i+1] <= x1) {
// hi = double(bucketCount[i+1]/(delta*(this->samples())));
// fN += hi*(x1N-x0N)/N;
// }
// }
// out << N
// << " "
// << fN << endl;
// }
// out.close();
// }
#endif /* HERWIG_Histogram_H */
diff --git a/Utilities/Histogram.icc b/Utilities/Histogram.icc
--- a/Utilities/Histogram.icc
+++ b/Utilities/Histogram.icc
@@ -1,97 +1,98 @@
// -*- C++ -*-
//
// Histogram.icc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the inlined member functions of
// the Histogram class.
//
namespace Herwig {
inline Histogram::Histogram(double lower, double upper, unsigned int nbin)
- : _globalStats(), _havedata(false), _bins(nbin+2),_prefactor(1.)
+ : _globalStats(), _havedata(false), _bins(nbin+2),_prefactor(1.),_total(0.)
{
if (upper<lower)
swap(upper,lower);
_bins[0].limit=-1.e100;
double limit(lower);
double width((upper-lower)/nbin);
for(unsigned int ix=1; ix <= nbin; ++ix)
{
_bins[ix].limit=limit;
limit += width;
}
_bins.back().limit=limit;
}
inline Histogram::Histogram(vector<double> limits)
- : _globalStats(), _havedata(false), _bins(limits.size()+1), _prefactor(1.)
+ : _globalStats(), _havedata(false), _bins(limits.size()+1), _prefactor(1.),_total(0.)
{
_bins[0].limit=-1.e100;
for (size_t i=1; i<=limits.size(); ++i)
_bins[i].limit=limits[i-1];
}
inline Histogram::Histogram(vector<double> limits, vector<double> data,
vector<double> dataerror)
- : _globalStats(), _havedata(true), _bins(limits.size()+1), _prefactor(1.)
+ : _globalStats(), _havedata(true), _bins(limits.size()+1), _prefactor(1.),_total(0.)
{
_bins[0].limit=-1.e100;
for (size_t i=1; i<=limits.size(); ++i)
_bins[i].limit=limits[i-1];
// no data goes into _bins[0] or _bins.back()!
for (size_t i=1; i<=min(limits.size()-1,data.size()); ++i)
_bins[i].data=data[i-1];
for (size_t i=1; i<=min(limits.size()-1,dataerror.size()); ++i)
_bins[i].dataerror=dataerror[i-1];
}
inline IBPtr Histogram::clone() const {
return new_ptr(*this);
}
inline IBPtr Histogram::fullclone() const {
return new_ptr(*this);
}
inline void Histogram::operator+=(double input) {
addWeighted(input,1.0);
}
inline void Histogram::addWeighted(double input, double weight) {
if(isnan(input)) return;
unsigned int ibin;
for(ibin=1; ibin<_bins.size(); ++ibin) {
if(input<_bins[ibin].limit)
break;
}
_bins[ibin-1].contents += weight;
_bins[ibin-1].contentsSq += sqr(weight);
_globalStats += weight * input;
+ _total += weight;
}
inline unsigned int Histogram::numberOfBins() const {
return _bins.size()-2;
}
inline double Histogram::prefactor() const {
return _prefactor;
}
inline void Histogram::prefactor(double in) {
_prefactor=in;
}
inline const Statistic & Histogram::globalStatistics() const {
return _globalStats;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 2:01 AM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4243584
Default Alt Text
(25 KB)

Event Timeline