diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 5d05249..78f81a0 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,26 +1,28 @@ set(example_dir "${CMAKE_CURRENT_SOURCE_DIR}") file(GLOB source_files ${example_dir}/*.cc) if(NOT ROOT_FOUND) - list(REMOVE_ITEM source_files "${example_dir}/AnalysisROOT.cc") + list(REMOVE_ITEM source_files "${example_dir}/AnalysisROOT.cc" "${example_dir}/WRapidityAnalysisROOT.cc") endif() foreach(source ${source_files}) get_filename_component(name ${source} NAME_WE) add_library(${name}_lib SHARED ${source}) set_target_properties(${name}_lib PROPERTIES OUTPUT_NAME "${name}") if (CMAKE_CXX_COMPILER_ID STREQUAL "Clang") # Clang warns about `std::unique_ptr` in C linked (`extern "C"`) function - # `make_analysis`. However this call is save from C++, and we don't care + # `make_analysis`. However this call is safe from C++, and we don't care # about pure C. We can ignore the warning. target_compile_options(${name}_lib PRIVATE "-Wno-return-type-c-linkage") endif() target_compile_options(${name}_lib PRIVATE "-fvisibility=hidden") ## link library target_link_libraries(${name}_lib HEJ) endforeach(source) if(ROOT_FOUND) - target_link_libraries(AnalysisROOT_lib ROOT::Core ROOT::Hist ROOT::RIO) + set(ROOT_LIBS ROOT::Core ROOT::Hist ROOT::RIO) + target_link_libraries(AnalysisROOT_lib $(ROOT_LIBS)) + target_link_libraries(WRapidityAnalysisROOT_lib $(ROOT_LIBS)) endif() diff --git a/examples/README.md b/examples/README.md new file mode 100644 index 0000000..f6a24f7 --- /dev/null +++ b/examples/README.md @@ -0,0 +1,26 @@ +# Example Analyses and Custom Scales for HEJ + +This directory contains a number of examples for HEJ analyses and +custom scale settings. Their usage is explained on +https://hej.hepforge.org/doc/current/user/analyses.html and +https://hej.hepforge.org/doc/current/user/scales.html +See the individual files for details. + +## Compilation + +The examples are compiled automatically together with the rest of HEJ, +see https://hej.hepforge.org/doc/current/user/installation.html for +instructions. The compiled libraries are saved in the `examples` +subdirectory of the build directory. + +To include examples requiring ROOT, you have to add +`-DEXCLUDE_ROOT=False` to the `cmake` command used to compile HEJ. You +can specify the ROOT installation directory with +`-DROOT_DIR=/directory/containing/ROOTConfig.cmake/`. + +It may also be necessary to preload the ROOT core library when running +HEJ or HEJFOG with a ROOT-based analysis: + +```sh +LD_PRELOAD=/path/to/libCore.so HEJ config.yml eventfile +``` diff --git a/examples/WRapidityAnalysisROOT.cc b/examples/WRapidityAnalysisROOT.cc new file mode 100644 index 0000000..b57f245 --- /dev/null +++ b/examples/WRapidityAnalysisROOT.cc @@ -0,0 +1,131 @@ +//! Simple HEJ analysis to plot y_W with ROOT +//! +//! See README.md for installation and running instructions. +//! +#include +#include +#include +#include // for std::unique_ptr +#include // for std::cout + +#include "HEJ/PDG_codes.hh" +#include "TFile.h" +#include "TH1.h" + +#include "HEJ/Analysis.hh" +#include "HEJ/Event.hh" +#include "HEJ/Version.hh" + +namespace YAML { + class Node; +} + +namespace LHEF { + class HEPRUP; +} +namespace { + void Fill2DHistogram(TH1D & h,double x,double wt) + { + const int binnumber=h.GetXaxis()->FindBin(x); + const double binwidth=h.GetXaxis()->GetBinWidth(binnumber); + h.Fill(x,wt/binwidth); + } + + class WRapidityAnalysisROOT: public HEJ::Analysis { + private: + TFile hfile_; + TH1D htotal_,hydif_,hyf_,hyb_,hyW_,hyWsampling_; + public: + WRapidityAnalysisROOT( + YAML::Node const & /* config */, LHEF::HEPRUP const & /* heprup */ + ):hfile_{"HEJ.root","RECREATE", + ("Output from "+HEJ::Version::package_name_full()).c_str()}, + htotal_{"total","total",1,-0.5,.5}, + hydif_{"ydif","ydif",40,0.,10.}, + hyf_{"yf","yf",40,-5.,5.}, + hyb_{"yb","yb",40,-5.,5.}, + hyW_{"yW","yW",40,-5.,5.}, + hyWsampling_{"yWsampling","yWsampling",40,-5.,5.} + { + // Setup histograms + htotal_.Sumw2(); + hydif_.Sumw2(); + hyf_.Sumw2(); + hyb_.Sumw2(); + hyW_.Sumw2(); + } + + void fill( + HEJ::Event const & event, + HEJ::Event const & /* FO_event */ + ) override { + const double weight = event.central().weight; + htotal_.Fill(0.0,weight); + auto const & jets = event.jets(); // jets are sorted in rapidity + // find the W momentum + auto const & outgoing = event.outgoing(); + auto the_boson = std::find_if( + outgoing.begin(), outgoing.end(), + [](auto const & particle) { return std::abs(particle.pid) == HEJ::pid::Wp; } + ); + // we already checked in `pass_cuts` that there is a boson, + // so the following assert cannot fail + assert(the_boson != outgoing.end()); + Fill2DHistogram(hydif_, + std::abs(jets.back().rapidity()-jets.front().rapidity()), weight); + Fill2DHistogram(hyf_, jets.front().rapidity(), weight); + Fill2DHistogram(hyb_, jets.back().rapidity(), weight); + Fill2DHistogram(hyW_, the_boson->rapidity(), weight); + Fill2DHistogram(hyWsampling_, the_boson->rapidity(),1); + } + + bool pass_cuts( + HEJ::Event const & event, + HEJ::Event const & /* FO_event */ + ) override { + if(event.jets().size() < 2) + return false; + auto const & outgoing = event.outgoing(); + + // Find the W boson + auto the_W_boson = std::find_if( + outgoing.begin(), outgoing.end(), + [](auto const & particle) { return std::abs(particle.pid) == HEJ::pid::Wp; } + ); + if(the_W_boson == outgoing.end()) + return false; + const auto W_pos = std::distance(outgoing.begin(), the_W_boson); + + // Find the decay products + const auto the_decay = event.decays().find(W_pos); + if(the_decay == event.decays().end()) + return false; + + // Find the charged lepton + const auto the_charged_lepton = std::find_if( + the_decay->second.begin(), the_decay->second.end(), + [](HEJ::Particle const & p) { + return HEJ::is_anylepton(p) && !HEJ::is_anyneutrino(p); + } + ); + if(the_charged_lepton == the_decay->second.end()) + return false; + return std::abs(the_charged_lepton->rapidity()) > 2.5; + } + + void finalise() override { + hfile_.Write(); + hfile_.Close(); + std::cout << "WRapidityAnalysisROOT bids you farewell " << std::endl; // be polite + } + + }; +} + +extern "C" +__attribute__((visibility("default"))) +std::unique_ptr make_analysis( + YAML::Node const & config, LHEF::HEPRUP const & heprup +){ + return std::make_unique(config, heprup); +}