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

analysis:
   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 three 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 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 two 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.

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

#include "HEJ/Analysis.hh"

class MyAnalysis: public HEJ::Analysis {
public:
  MyAnalysis(YAML::Node const & /* config */) {}

  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 finalise() override {
  }

};

extern "C"
std::unique_ptr<HEJ::Analysis> make_analysis(
    YAML::Node const & config
){
  return std::make_unique<MyAnalysis>(config);
}

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 -Wl,-soname,libmyanalysis.so -o libmyanalysis.so myanalysis.cc

With g++ 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

analysis:
   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

analysis:
   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

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

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

#include "yaml-cpp/yaml.h"

class MyAnalysis: public HEJ::Analysis {
public:
  MyAnalysis(YAML::Node const & config):
  xsection_{0.}, xsection_error_{0.},
  outfile_{config["output"].as<std::string>()}
  {}

  void fill(
    HEJ::Event const & event,
    HEJ::Event const & /* FO_event */
  ) override {
    const double wt = event.central().weight;
    xsection_ += wt;
    xsection_error_ += wt*wt;
  }

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

  void finalise() override {
    std::cout << "cross section: " << xsection_ << " +- "
      << std::sqrt(xsection_error_) << "\n";
    std::ofstream fout{outfile_};
    fout << "cross section: " << xsection_ << " +- "
      << std::sqrt(xsection_error_) << "\n";
  }

  private:
  double xsection_, xsection_error_;
  std::string outfile_;

};

extern "C"
std::unique_ptr<HEJ::Analysis> make_analysis(
    YAML::Node const & config
){
  return std::make_unique<MyAnalysis>(config);
}