Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/AnalysisHandler.hh"
#include "Rivet/Projections/DirectFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/VetoedFinalState.hh"

namespace Rivet {


  /// Analysis of the generated event-weight distributions
  class MC_WEIGHTS : public Analysis {
  public:

    RIVET_DEFAULT_ANALYSIS_CTOR(MC_WEIGHTS);


    /// @name Analysis methods
    /// @{

    /// Book histograms and initialise projections before the run
    void init() {

      DirectFinalState leps(Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
      VetoedFinalState vfs;
      vfs.vetoFinalState(leps);
      declare(FastJets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::DECAY, JetInvisibles::DECAY), "jets");

      book(_n["weight_10+7"], "weight_10+7", 200, -1e+03, 1e+03);
      book(_n["weight_10+6"], "weight_10+6", 200, -1e+02, 1e+02);
      book(_n["weight_10+5"], "weight_10+5", 200, -1e+01, 1e+01);
      book(_n["weight_10+4"], "weight_10+4", 200, -1e+04, 1e+04);
      book(_n["weight_10+3"], "weight_10+3", 200, -1e+03, 1e+03);
      book(_n["weight_10+2"], "weight_10+2", 200, -1e+02, 1e+02);
      book(_n["weight_10"],   "weight_10",   200, -1e+01, 1e+01);
      book(_n["weight_1"],    "weight_1",    200, -1e+00, 1e+00);

      book(_n["logweight_pos"], "logweight_pos", logspace(100, 1e-5, 1e5));
      book(_n["logweight_neg"], "logweight_neg", logspace(100, 1e-5, 1e5));

      book(_h["xs_10+7"], "xs_10+7", 200, -1e+03, 1e+03);
      book(_h["xs_10+6"], "xs_10+6", 200, -1e+02, 1e+02);
      book(_h["xs_10+5"], "xs_10+5", 200, -1e+01, 1e+01);
      book(_h["xs_10+4"], "xs_10+4", 200, -1e+04, 1e+04);
      book(_h["xs_10+3"], "xs_10+3", 200, -1e+03, 1e+03);
      book(_h["xs_10+2"], "xs_10+2", 200, -1e+02, 1e+02);
      book(_h["xs_10"],   "xs_10",   200, -1e+01, 1e+01);
      book(_h["xs_1"],    "xs_1",    200, -1e+00, 1e+00);

      tripleBook("Njets",    7, -0.5,   6.5);
      tripleBook("Nbjets",   5, -0.5,   4.5);
      tripleBook("HT",      50,   0., 1500.);
      tripleBook("jet1pT", 100,  20., 1000.);
      tripleBook("bjet1pT", 50,  20.,  500.);

      book(_e_xsfraction_neg, "xsfraction_neg");
      book(_e_fracneff, "neff_frac");
      book(_e_neff, "neff");
    }

    void tripleBook(const string& name, size_t nbins, double xmin, double xmax) {
      book(_h[name+"_neg"], name+"_neg", nbins, xmin, xmax);
      book(_h[name+"_pos"], name+"_pos", nbins, xmin, xmax);
      book(_e[name], name, nbins, xmin, xmax);
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {

      Jets jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 20*GeV);
      Jets bjets = select(jets, hasBTag(Cuts::pT > 5*GeV));

      const size_t numWeights = event.weights().size();
      for (size_t m = 0; m < numWeights; ++m) {
        const double weight = event.weights()[m];
        _n["weight_10+7"].get()->persistent(m)->fill(weight/1e4, 1.0);
        _n["weight_10+6"].get()->persistent(m)->fill(weight/1e4, 1.0);
        _n["weight_10+5"].get()->persistent(m)->fill(weight/1e4, 1.0);
        _n["weight_10+4"].get()->persistent(m)->fill(weight, 1.0);
        _n["weight_10+3"].get()->persistent(m)->fill(weight, 1.0);
        _n["weight_10+2"].get()->persistent(m)->fill(weight, 1.0);
        _n["weight_10"].get()->persistent(m)->fill(weight, 1.0);
        _n["weight_1"].get()->persistent(m)->fill(weight, 1.0);

        _h["xs_10+7"].get()->persistent(m)->fill(weight/1e4, fabs(weight));
        _h["xs_10+6"].get()->persistent(m)->fill(weight/1e4, fabs(weight));
        _h["xs_10+5"].get()->persistent(m)->fill(weight/1e4, fabs(weight));
        _h["xs_10+4"].get()->persistent(m)->fill(weight, fabs(weight));
        _h["xs_10+3"].get()->persistent(m)->fill(weight, fabs(weight));
        _h["xs_10+2"].get()->persistent(m)->fill(weight, fabs(weight));
        _h["xs_10"].get()->persistent(m)->fill(weight, fabs(weight));
        _h["xs_1"].get()->persistent(m)->fill(weight, fabs(weight));

        if (weight < 0.) {
          _n["logweight_neg"].get()->persistent(m)->fill(fabs(weight), 1.0);
          _h["Njets_neg"].get()->persistent(m)->fill(jets.size(), fabs(weight));
          _h["Nbjets_neg"].get()->persistent(m)->fill(bjets.size(), fabs(weight));
          _h["HT_neg"].get()->persistent(m)->fill(sum(jets, Kin::pT, 0.0)/GeV, fabs(weight));
          if (jets.size()) {
            _h["jet1pT_neg"].get()->persistent(m)->fill(jets[0].pT()/GeV, fabs(weight));
          }
          if (bjets.size()) {
            _h["bjet1pT_neg"].get()->persistent(m)->fill(bjets[0].pT()/GeV, fabs(weight));
          }
        } else {
          _n["logweight_pos"].get()->persistent(m)->fill(weight, 1.0);
          _h["Njets_pos"].get()->persistent(m)->fill(jets.size(), weight);
          _h["Nbjets_pos"].get()->persistent(m)->fill(bjets.size(), weight);
          _h["HT_pos"].get()->persistent(m)->fill(sum(jets, Kin::pT, 0.0)/GeV, weight);
          if (jets.size()) {
            _h["jet1pT_pos"].get()->persistent(m)->fill(jets[0].pT()/GeV, weight);
          }
          if (bjets.size()) {
            _h["bjet1pT_pos"].get()->persistent(m)->fill(bjets[0].pT()/GeV, weight);
          }
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {

      scale(_n, 1.0 / numEvents());
      scale(_h, crossSection() / picobarn / sumOfWeights());

      const double totalSumW  = _n["logweight_neg"]->sumW()  + _n["logweight_pos"]->sumW();
      const double totalSumW2 = _n["logweight_neg"]->sumW2() + _n["logweight_pos"]->sumW2();
      const double negFrac    = _n["logweight_neg"]->sumW() / totalSumW;
      const double negFracErr = negFrac * totalSumW / sqrt(totalSumW2);

      for (auto& item : _e) {
        const string& aop = item.second->path();
        *item.second = *_h[item.first+"_neg"] / (*_h[item.first+"_neg"] + *_h[item.first+"_pos"]);
        item.second->setPath(aop);
      }

      _e_xsfraction_neg->set(negFrac, negFracErr);
      _e_neff->setVal(effNumEvents());
      _e_fracneff->setVal(effNumEvents()/numEvents());
    }

    /// @}


    /// @name Histograms
    /// @{
    map<string,Histo1DPtr> _h, _n;
    map<string,Estimate1DPtr> _e;
    Estimate0DPtr _e_xsfraction_neg, _e_neff, _e_fracneff;
    /// @}

  };



  RIVET_DECLARE_PLUGIN(MC_WEIGHTS);

}