Writing custom analyses

HEJ 2 and the HEJ fixed-order generator can generate HepMC files, so you can always run a Rivet analysis on these. However if you compiled HEJ 2 with Rivet you can use the native Rivet interface. For example

analyses:
   - rivet: [MC_XS, MC_JETS]
     output: HEJ

would call the generic MC_XS and MC_JETS analysis and write the result into HEJ[.Scalename].yoda. HEJ 2 will then run Rivet over all different scales seperatly and write out each into a different yoda file. Alternatively instead of using Rivet, you can provide a custom analysis inside a C++ library.

An analysis is a class that derives from the abstract Analysis base class provided by HEJ 2. It has to implement four public functions:

• The pass_cuts member function return true if and only if the given event (first argument) passes the analysis cuts

• The fill member function adds an event to the analysis, which for example can be used to fill histograms. HEJ 2 will only pass events for which pass_cuts has returned true.

• The set_xs_scale member function updates the ratio between the total cross section and the sum of event weights.

• The finalise member function is called after all events have been processed. It can be used, for example, to print out or save the analysis results.

The pass_cuts and fill functions take two arguments: the resummation event generated by HEJ 2 and the original fixed-order input event. Usually, the second argument can be ignored. It can be used, for example, for implementing cuts that depend on the ratio of the weights between the fixed-order and the resummation event.

In addition to the three member functions, there has to be a global make_analysis function that takes the analysis parameters in the form of a YAML Node and returns a std::unique_ptr to the Analysis.

The following code creates the simplest conceivable analysis.

//! simple HEJ analysis that doesn't do anything

#include <memory>     // for std::unique_ptr
#include <iostream>

#include "HEJ/Analysis.hh"

namespace YAML {
  class Node;
}

namespace LHEF {
  class HEPRUP;
}

namespace {
  class AnalysisTemplate: public HEJ::Analysis {
    public:
      AnalysisTemplate(
        YAML::Node const & /* config */, LHEF::HEPRUP const & /* heprup */
      ) {}

      void fill(
        HEJ::Event const & /* event */,
        HEJ::Event const & /* FO_event */
      ) override {
      }

      bool pass_cuts(
        HEJ::Event const & /* event */,
        HEJ::Event const & /* FO_event */
      ) override {
        return true;
      }

      void set_xs_scale(double /* xsscale */) override {
        // typically store this information for later
      }

      void finalise() override {
        std::cout << "bye" << std::endl; // be polite
      }
    };
}

extern "C"
__attribute__((visibility("default")))
std::unique_ptr<HEJ::Analysis> make_analysis(
    YAML::Node const & config, LHEF::HEPRUP const & heprup
){
  return std::make_unique<AnalysisTemplate>(config, heprup);
}

You can save this code to a file, for example myanalysis.cc, and compile it into a shared library. Using the g++ compiler, the library can be built with

g++ $(HEJ-config --cxxflags) -fPIC -shared -O2 -fvisibility=hidden -Wl,-soname,libmyanalysis.so -o libmyanalysis.so myanalysis.cc

It is also good practice to add __attribute__((visibility("default"))) after extern "C" in the above code snippet and then compile with the additional flag -fvisibility=hidden to prevent name clashes.

You can use the analysis in HEJ 2 or the HEJ fixed-order generator by adding

analyses:
   - plugin: /path/to/libmyanalysis.so

to the .yml configuration file.

As a more interesting example, here is the code for an analysis that sums up the total cross section and prints the result to both standard output and a file specified in the .yml config with

analyses:
   - plugin: analysis/build/directory/src/libmy_analysis.so
     output: outfile

To access the configuration at run time, HEJ 2 uses the yaml-cpp library; for more details see the yaml-cpp tutorial. The analysis code itself is

//! HEJ analysis to output the cross section to a file

#include <cmath>
#include <fstream>
#include <iostream>
#include <memory>
#include <string>

#include "HEJ/Analysis.hh"
#include "HEJ/Event.hh"

#include "yaml-cpp/yaml.h"

#include "LHEF/LHEF.h"

namespace {
  class AnalysisPrint: public HEJ::Analysis {
  public:
    AnalysisPrint(YAML::Node const & config, LHEF::HEPRUP const & heprup):
      wt_sum_{0.}, wt2_sum_{0.},
      outfile_{config["output"].as<std::string>()},
      generators_{heprup.generators}
    {}

    void fill(
      HEJ::Event const & event,
      HEJ::Event const & /* FO_event */
    ) override {
      const double wt = event.central().weight;
      wt_sum_ += wt;
      wt2_sum_ += wt*wt; // this error estimate is too small
    }

    bool pass_cuts(
      HEJ::Event const & /* event */,
      HEJ::Event const & /* FO_event */
    ) override {
      return true;
    }

    void set_xs_scale(double xsscale) override {
      xs_scale_ = xsscale;
    }

    void finalise() override {
      // compute cross section from sum of weights and scaling factor
      const double xsection = wt_sum_ * xs_scale_;
      const double xsection_error = std::sqrt(wt2_sum_) * xs_scale_;

      // print to screen
      std::cout << "Generated with:\n";
      for(auto const & generator: generators_)
        std::cout << generator.name << " " << generator.version << "\n";
      std::cout << "cross section: " << xsection << " +- "
        << xsection_error << std::endl;

      // print to file
      std::ofstream fout{outfile_};
      fout << "Generated with\n";
      for(auto const & generator: generators_)
        fout << generator.name << " " << generator.version << "\n";
      fout << "cross section: " << xsection << " +- "
        << xsection_error << std::endl;
    }

    private:
    double wt_sum_, wt2_sum_;
    double xs_scale_{1.0};
    std::string outfile_;
    std::vector<LHEF::Generator> generators_;
  };
}

extern "C"
__attribute__((visibility("default")))
std::unique_ptr<HEJ::Analysis> make_analysis(
    YAML::Node const & config, LHEF::HEPRUP const & heprup
){
  return std::make_unique<AnalysisPrint>(config, heprup);
}