Histogram2.h

// -*- C++ -*-

// (C) 2007-2009 Simon Plaetzer -- sp@particle.uni-karlsruhe.de

#ifndef Analysis2_Histogram2_H
#define Analysis2_Histogram2_H
//
// This is the declaration of the Histogram2 class.
//

#include <cassert>

#include "ThePEG/Interface/Interfaced.h"
#include "Histogram2.fh"

#include "ThePEG/Utilities/StringUtils.h"

namespace Analysis2 {

using namespace ThePEG;

class HistogramChannel {

public:


  inline HistogramChannel ();

  inline explicit HistogramChannel (unsigned int bins, bool counting = true);

  inline explicit HistogramChannel (const vector<pair<double,double> >&,
                                    double underflow = 0., double overflow = 0.);

public:

  inline void book (unsigned int bin, double weight = 1.);

  inline void bookUnderflow (double weight = 1.);

  inline void bookOverflow (double weight = 1.);

  inline void nanEvent ();

public:

  HistogramChannel& operator += (const HistogramChannel&);

  HistogramChannel& operator -= (const HistogramChannel&);

  HistogramChannel& operator *= (const HistogramChannel&);

  HistogramChannel& operator *= (double);

  HistogramChannel& operator += (double);

  HistogramChannel& operator *= (pair<double,double>);

  HistogramChannel& operator /= (const HistogramChannel&);

  inline HistogramChannel& operator /= (double);

  HistogramChannel& operator /= (pair<double,double>);

public:

  inline bool isCountingChannel () const;

  inline vector<pair<double,double> > bins () const;

  inline pair<double,double> bin (unsigned int) const;

  inline vector<unsigned long> binEntries () const;

  inline unsigned long binEntries (unsigned int) const;

  inline void bin (unsigned int, pair<double,double>, unsigned long en = 0);

  inline pair<double,double> outOfRange () const;

  inline unsigned long visible () const;

  inline unsigned long total () const;

  inline vector<unsigned long> nanWeights () const;

  unsigned long nanWeightEvents () const;

  inline unsigned long nanEvents () const;

public:


  inline void finish ();

  inline double binVariance (unsigned int) const;

  inline double weightMean (unsigned int) const;

  inline double weightVariance (unsigned int) const;

  void differential (const vector<pair<double,double> >&);

  pair<double,double> binSum () const;

  inline pair<double,double> binAverage () const;

  pair<double,double> integrate (const vector<pair<double,double> >&) const;

  pair<double,double> average (const vector<pair<double,double> >&) const;

  HistogramChannel delta (const HistogramChannel& channel) const;

  HistogramChannel chi2 (const HistogramChannel& channel, double minfrac = .0) const;

  HistogramChannel profile () const;


public:


  void write (ostream&, const string&);

  string read (istream&);


public:

  void persistentOutput(PersistentOStream & os) const;

  void persistentInput(PersistentIStream & is);

private:

  bool _isCountingChannel;

  vector<pair<double,double> > _bins;

  vector<unsigned long> _binEntries;

  pair<double,double> _outOfRange;

  unsigned long _visible;

  unsigned long _total;

  unsigned long _nanEvents;

  vector<unsigned long> _nanWeights;

  bool _finished;

};

inline PersistentOStream& operator << (PersistentOStream& os, const HistogramChannel&);

inline PersistentIStream& operator >> (PersistentIStream& is, HistogramChannel&);

inline HistogramChannel operator + (const HistogramChannel& a, const HistogramChannel& b);

inline HistogramChannel operator - (const HistogramChannel& a, const HistogramChannel& b);

inline HistogramChannel operator * (const HistogramChannel& a, const HistogramChannel& b);

inline HistogramChannel operator / (const HistogramChannel& a, const HistogramChannel& b);

inline string getNextTag (istream&);

template<class T>
inline void fromString (const string& str, T& result) {
  istringstream buffer (str);
  buffer >> result;
}

namespace ChannelOutput {

  const unsigned int Default     = 0;

  const unsigned int Bincenters  = 1;

  const unsigned int NoErrorbars = 1 << 1;

  const unsigned int NanEvents   = 1 << 2;

  const unsigned int Statistics  = 1 << 3;

}

class Histogram2: public Interfaced {

public:


  inline Histogram2 ();

  virtual ~Histogram2();

public:

  explicit Histogram2 (double low, double high,
                       unsigned int bins,
                       const string& name = "mc");

  explicit Histogram2 (const vector<pair<double,double> >&, const string& name);

  explicit Histogram2 (const string& dataFile, const string& dataName);

  inline void setName (const string&);


public:

  void book (const string& name, double event, double weight = 1.);

  inline bool haveChannel (const string&) const;

  inline HistogramChannel& channel (const string&);

  inline HistogramChannel channel (const string&) const;

  inline void insertChannel (const string&, const HistogramChannel&);

  inline void insertChannel (const string&);

  inline HistogramChannel removeChannel (const string&);

public:

  inline const vector<pair<double,double> >& binning () const;

  inline const map<double,unsigned int>& binhash () const;

  inline pair<double,double> range () const;

public:


  inline void finish (const string&);

  inline void differential (const string&);

  inline pair<double,double> integrate (const string&) const;

  inline void normalise (const string&);

  inline void normalise (const string&, const string&);

  inline void rescale (const string&, double);

  inline void rescale (double);

  inline void normaliseToCrossSection(const string& name = "mc");

  inline pair<double,double> chi2perDOF (const string&, const string&, double minfrac = .0) const;

public:

  vector<string> channels () const;

  void output (ostream&, const string& name,
               unsigned int flags = 0,
               char comment = '#') const;

public:

  /*@name Support for parallel runs */

  void store (const string& name);

  bool load (const string& name);

  Histogram2Ptr loadToHistogram (const string& name) const;

  void combine (const string& prefix, const string& name,
                unsigned int numRuns, const string& dataChannel = "",
                const string& mcChannel = "MC");

  inline CrossSection xSec () const;

  inline void xSec (CrossSection);


public:


  explicit Histogram2 (vector<double> limits);

  explicit Histogram2 (vector<double> limits,
                       vector<double> data,
                       vector<double> dataerror);

  inline void operator+=(double);

  inline void addWeighted(double data, double weight);

  inline unsigned int numberOfBins() const;

  inline void normaliseToData();

  inline void chiSquared(double & chisq, 
                         unsigned int & ndegrees,
                         double minfrac=0.) const;

  Histogram2 ratioWith(const Histogram2& h2) const;


public:

  void persistentOutput(PersistentOStream & os) const;

  void persistentInput(PersistentIStream & is, int version);

  static void Init();

protected:

  inline virtual IBPtr clone() const;

  inline virtual IBPtr fullclone() const;

// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).

private:

  vector<pair<double,double> > _binning;

  pair<double,double> _range;

  map<double, unsigned int> _binhash;

  map<string,HistogramChannel> _channels;

  CrossSection _xSec;

  static ClassDescription<Histogram2> initHistogram2;

  Histogram2 & operator=(const Histogram2 &);

};

}

#include "ThePEG/Utilities/ClassTraits.h"

namespace ThePEG {

template <>
struct BaseClassTrait<Analysis2::Histogram2,1> {
  typedef Interfaced NthBase;
};

template <>
struct ClassTraits<Analysis2::Histogram2>
  : public ClassTraitsBase<Analysis2::Histogram2> {
  static string className() { return "Analysis2::Histogram2"; }
  static string library() { return "Analysis2.so"; }
};

}

#include "Histogram2.icc"
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
// #include "Histogram2.tcc"
#endif

#endif /* Analysis2_Histogram2_H */