1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/AnalysisHandler.hh"
 4
 5#ifndef RIVET_ENABLE_HEPMC_3
 6#include "HepMC/HepMCDefs.h"
 7#endif
 8
 9namespace Rivet {
10
11
12  /// Analysis of the generated cross-section
13  class MC_XS : public Analysis {
14  public:
15
16    RIVET_DEFAULT_ANALYSIS_CTOR(MC_XS);
17
18
19    /// @name Analysis methods
20    /// @{
21
22    /// Book histograms and initialise projections before the run
23    void init() {
24      /// @todo Convert to Scatter1D or Counter
25      book(_h_XS, "XS", {-0.5, 0.5});
26      book(_h_N, "N", 1, 0.0, 1.0);
27      book(_h_pmXS, "pmXS", 2, -1.0, 1.0);
28      book(_h_pmN, "pmN", 2, -1.0, 1.0);
29    }
30
31
32    /// Perform the per-event analysis
33    void analyze(const Event& event) {
34
35      const size_t numWeights = event.weights().size();
36      const vector<pair<double,double>> xsecs = event.crossSections();
37      for (size_t m = 0; m < numWeights; ++m) {
38        #if defined RIVET_ENABLE_HEPMC_3 || defined HEPMC_HAS_CROSS_SECTION
39        size_t idx = (xsecs.size() == numWeights)? m : 0;
40        const double xs    = xsecs[idx].first;
41        const double xserr = xsecs[idx].second;
42        _h_XS.get()->_getPersistent(m)->point(0).setY(xs, xserr);
43        # endif
44        const double weight = event.weights()[m];
45        _h_pmXS.get()->_getPersistent(m)->fill(0.5*(weight > 0 ? 1. : -1), abs(weight));
46        _h_pmN.get()->_getPersistent(m)->fill(0.5*(weight > 0 ? 1. : -1), 1.);
47        _h_N.get()->_getPersistent(m)->fill(0.5, 1.0);
48      }
49
50    }
51
52
53    /// Normalise histograms etc., after the run
54    void finalize() {
55      scale(_h_pmXS, crossSection()/sumOfWeights());
56      #if !defined RIVET_ENABLE_HEPMC_3 && !defined HEPMC_HAS_CROSS_SECTION
57      _h_XS->point(0).setY(crossSection(), 0.);
58      #endif
59    }
60
61    /// @}
62
63
64    /// @name Histograms
65    /// @{
66    Scatter2DPtr _h_XS;
67    Histo1DPtr _h_pmXS, _h_pmN, _h_N;
68    /// @}
69
70  };
71
72
73
74  RIVET_DECLARE_PLUGIN(MC_XS);
75
76}