[Rivet-svn] r4300 - contrib/nlohisto

[blackhole at projects.hepforge.org]

Author: fsiegert
Date: Tue May 21 15:55:04 2013
New Revision: 4300

Log:
Add NLOHistogram1D as workaround for Rivet 1.x where correlated events are not taken into account automatically when calculating histogram errors.

Added:
   contrib/nlohisto/
   contrib/nlohisto/NLOHistogram1D.README
   contrib/nlohisto/NLOHistogram1D.cc

Added: contrib/nlohisto/NLOHistogram1D.README
==============================================================================
--- /dev/null	00:00:00 1970	(empty, because file is newly added)
+++ contrib/nlohisto/NLOHistogram1D.README	Tue May 21 15:55:04 2013	(r4300)
@@ -0,0 +1,49 @@
+This file adds a hack-ish way of allowing to take event-correlations into
+account when filling a histogram. It assumes, that correlated events are
+filled into Rivet subsequently and with the same (HepMC) event number.
+
+You can add the NLOHistogram1D functionality to your Rivet 1.x analysis like:
+
+---------------------------------------------------
+#include "Rivet/Analysis.hh"
+#include "Rivet/RivetAIDA.hh"
+#include "LWH/Histogram1D.h"
+
+namespace Rivet {
+
+  class MY_ANALYSIS : public Analysis {
+
+#include "NLOHistogram1D.cc"
+
+  public:
+
+    MY_ANALYSIS() : Analysis("MY_ANALYSIS")
+    {
+    }
+
+    void init() {
+      _h_pT_jet = bookNLOHistogram1D("pT_jet", logBinEdges(100, 20.0, 200.0));
+      _h_eta_jet = bookNLOHistogram1D("eta_jet", 45, -4.5, 4.5);
+    }
+
+    void analyze(const Event & event) {
+      Jet jet;
+      _h_pT_jet->fill(jet.momentum().pT(), event);
+      _h_eta_jet[i]->fill(jet.momentum().eta(), event);
+    }
+
+    void finalize() {
+      double norm=crossSection()/sumOfWeights();
+      _h_pT_jet->finalize(this, norm);
+      _h_eta_jet->finalize(this, norm);
+    }
+
+  private:
+
+    NLOHistogram1D * _h_pT_jet;
+    NLOHistogram1D * _h_eta_jet;
+  };
+
+}
+---------------------------------------------------
+

Added: contrib/nlohisto/NLOHistogram1D.cc
==============================================================================
--- /dev/null	00:00:00 1970	(empty, because file is newly added)
+++ contrib/nlohisto/NLOHistogram1D.cc	Tue May 21 15:55:04 2013	(r4300)
@@ -0,0 +1,125 @@
+class NLOHistogram1D : public AIDA::IHistogram1D {
+
+  AIDA::IHistogramFactory& _factory;
+  AIDA::IHistogram1D* _hist;
+  LWH::Histogram1D* _tmphist;
+  std::string _path;
+
+  int _current_event_number;
+
+  void _syncHists() {
+    const IAxis & ax = _tmphist->axis();
+    // the following loop is mightily inefficient, but to make it more
+    // efficient one needs to modify Rivet internals
+    for (int i=0; i<ax.bins(); ++i) {
+      _hist->fill(ax.binLowerEdge(i)+ax.binWidth(i)/2.0,
+                  _tmphist->binHeight(i));
+    }
+    _tmphist->reset();
+  }
+
+public:
+
+  NLOHistogram1D(AIDA::IHistogramFactory& factory, const string& path,
+                 size_t nbins, double lower, double upper) :
+    _factory(factory), _path(path), _current_event_number(-1) {
+    _hist = factory.createHistogram1D(_path, "", nbins, lower, upper);
+    _tmphist = dynamic_cast<LWH::Histogram1D*>
+      (factory.createHistogram1D(_path+"_tmp", "", nbins, lower, upper));
+    if (!_tmphist) throw Error("Cast failed.");
+  }
+  
+  NLOHistogram1D(AIDA::IHistogramFactory& factory, const string& path,
+                 const vector<double>& binedges) :
+    _factory(factory), _path(path), _current_event_number(-1) {
+    _hist = factory.createHistogram1D(_path, "", binedges);
+    _tmphist = dynamic_cast<LWH::Histogram1D*>
+      (factory.createHistogram1D(_path+"_tmp", "", binedges));
+    if (!_tmphist) throw Error("Cast failed.");
+  }
+
+  ~NLOHistogram1D() {}
+    
+  bool fill(double x, double weight = 1.) {
+    std::cout<<x<<" "<<weight<<std::endl;
+    throw Error("Called fill function without event.");
+    return false;
+  }
+
+  bool fill(double x, const Event& event)
+  {
+    if (_current_event_number==-1)
+      _current_event_number = event.genEvent().event_number();
+
+    if (event.genEvent().event_number()!=_current_event_number) {
+      _syncHists();
+      _current_event_number = event.genEvent().event_number();
+    }
+
+    return _tmphist->fill(x, event.weight());
+  }
+
+  void finalize(Analysis* ana, const double& xs_per_event)
+  {
+    _syncHists();
+    ana->scale(_hist, xs_per_event);
+    _factory.destroy(_tmphist);
+  }
+
+
+  double binMean(int index) const { return _hist->binMean(index); }
+  int binEntries(int index) const { return _hist->binEntries(index); }
+  double binHeight(int index) const { return _hist->binHeight(index); }
+  double binError(int index) const { return _hist->binError(index); }
+  double mean() const { return _hist->mean(); }
+  double rms() const { return _hist->rms(); }
+  const IAxis & axis() const { return _hist->axis(); }
+  int coordToIndex(double coord) const { return _hist->coordToIndex(coord); }
+  int allEntries() const { return _hist->allEntries(); }
+  int extraEntries() const { return _hist->extraEntries(); }
+  double equivalentBinEntries() const { return _hist->equivalentBinEntries(); }
+  double sumBinHeights() const { return _hist->sumBinHeights(); }
+  double sumAllBinHeights() const { return _hist->sumAllBinHeights(); }
+  double sumExtraBinHeights() const { return _hist->sumExtraBinHeights(); }
+  double minBinHeight() const { return _hist->minBinHeight(); }
+  double maxBinHeight() const { return _hist->maxBinHeight(); }
+  int dimension() const { return _hist->dimension(); }
+  bool reset() { return _hist->reset(); }
+  int entries() const { return _hist->entries(); }
+    
+  bool add(const IHistogram1D & hist) {
+    std::cout<<&hist<<std::endl;
+    throw Error("Called add function in NLO histo.");
+    return false;
+  }
+  bool scale(double scaleFactor) {
+    std::cout<<scaleFactor<<std::endl;
+    throw Error("Called scale function in NLO histo.");
+    return false;
+  }
+
+};
+
+
+NLOHistogram1D* bookNLOHistogram1D(const string& hname,
+                                   size_t nbins, double lower, double upper)
+{
+  myMakeHistoDir();
+  return new NLOHistogram1D(histogramFactory(), histoPath(hname), nbins, lower, upper);
+}
+
+NLOHistogram1D* bookNLOHistogram1D(const string& hname, const vector<double>& binedges)
+{
+  myMakeHistoDir();
+  return new NLOHistogram1D(histogramFactory(), histoPath(hname), binedges);
+}
+
+
+mutable bool _myMadeHistoDir;
+void myMakeHistoDir() {
+  if (!_myMadeHistoDir) {
+    if (! name().empty()) tree().mkdirs(histoDir());
+    _myMadeHistoDir = true;
+  }
+}
+