1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/AnalysisHandler.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// Analysis of the generated event-weight distributions
 9  class MC_WEIGHTS : public Analysis {
10  public:
11
12    RIVET_DEFAULT_ANALYSIS_CTOR(MC_WEIGHTS);
13
14
15    /// @name Analysis methods
16    /// @{
17
18    /// Book histograms and initialise projections before the run
19    void init() {
20      /// @todo Convert to Scatter1D or Counter
21      book(_h_weight_100, "weight_100", 200, -100.0, 100.0);
22      book(_h_weight_10,  "weight_10",  200,  -10.0,  10.0);
23      book(_h_logweight_pos, "logweight_pos", logspace(100, 0.1, 10000.0));
24      book(_h_logweight_neg, "logweight_neg", logspace(100, 0.1, 10000.0));
25
26      const vector<string> edges = { "Negative weight fraction" };
27      book(_h_xsfraction_neg, "xsfraction_neg", edges);
28    }
29
30
31    /// Perform the per-event analysis
32    void analyze(const Event& event) {
33
34      const size_t numWeights = event.weights().size();
35      for (size_t m = 0; m < numWeights; ++m) {
36        const double weight = event.weights()[m];
37        _h_weight_100.get()->persistent(m)->fill(weight, 1.0);
38        _h_weight_10.get()->persistent(m)->fill(weight, 1.0);
39        if (weight < 0.) {
40          _h_logweight_neg.get()->persistent(m)->fill(fabs(weight), 1.0);
41        } else {
42          _h_logweight_pos.get()->persistent(m)->fill(weight, 1.0);
43        }
44      }
45    }
46
47
48    /// Normalise histograms etc., after the run
49    void finalize() {
50      const double sf = 1.0 / numEvents();
51      scale(_h_weight_100, sf);
52      scale(_h_weight_10, sf);
53      scale(_h_logweight_pos, sf);
54      scale(_h_logweight_neg, sf);
55
56      const double totalSumW  = _h_logweight_neg->sumW() + _h_logweight_pos->sumW();
57      const double totalSumW2 = _h_logweight_neg->sumW2() + _h_logweight_pos->sumW2();
58      const double negFrac = _h_logweight_neg->sumW() / totalSumW;
59      const double negFracErr = negFrac * totalSumW / sqrt(totalSumW2);
60      _h_xsfraction_neg->bin(1).set(negFrac, negFracErr);
61    }
62
63    /// @}
64
65
66    /// @name Histograms
67    /// @{
68    BinnedEstimatePtr<string> _h_xsfraction_neg;
69    Histo1DPtr _h_weight_100, _h_weight_10, _h_logweight_pos, _h_logweight_neg;
70    /// @}
71
72  };
73
74
75
76  RIVET_DECLARE_PLUGIN(MC_WEIGHTS);
77
78}