.. _`Writing custom analyses`:

Writing custom analyses
=======================

HEJ 2 and the HEJ fixed-order generator can generate HepMC files, so you
can always run a `Rivet <https://rivet.hepforge.org/>`_ analysis on
these. However if you compiled HEJ 2 with Rivet you can use the native
Rivet interface. For example

.. code-block:: YAML

  analysis:
     rivet: [MC_XS, MC_JETS]
     output: HEJ

would call the generic
`MC_XS <https://rivet.hepforge.org/analyses/MC_XS.html>`_ and
`MC_JETS <https://rivet.hepforge.org/analyses/MC_JETS.html>`_ analysis
and write the result into :code:`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 :code:`Analysis`
base class provided by HEJ 2. It has to implement three public
functions:

* The :code:`pass_cuts` member function return true if and only if the
  given event (first argument) passes the analysis cuts

* The :code:`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 :code:`pass_cuts` has returned true.

* The :code:`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 :code:`pass_cuts` and :code:`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 :code:`Node` and returns a :code:`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 :code:`myanalysis.cc`, and
compile it into a shared library. Using the :code:`g++` compiler, the
library can be built with

.. code-block:: sh

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

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

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

.. code-block:: YAML

  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

.. code-block:: YAML

  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
<https://github.com/jbeder/yaml-cpp/wiki/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);
 }