diff --git a/configure.ac b/configure.ac
--- a/configure.ac
+++ b/configure.ac
@@ -1,315 +1,315 @@
## Process this file with autoconf to produce a configure script.
AC_PREREQ(2.59)
AC_INIT([Rivet],[2.6.0-pre],[rivet@projects.hepforge.org],[Rivet])
## Check and block installation into the src/build dir
if test "$prefix" = "$PWD"; then
AC_MSG_ERROR([Installation into the build directory is not supported: use a different --prefix argument])
fi
## Force default prefix to have a path value rather than NONE
if test "$prefix" = "NONE"; then
prefix=/usr/local
fi
AC_CONFIG_SRCDIR([src/Core/Analysis.cc])
AC_CONFIG_HEADERS([include/Rivet/Config/DummyConfig.hh include/Rivet/Config/RivetConfig.hh include/Rivet/Config/BuildOptions.hh])
AM_INIT_AUTOMAKE([dist-bzip2 -Wall])
m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
m4_ifdef([AM_PROG_AR], [AM_PROG_AR])
AC_CONFIG_MACRO_DIR([m4])
AC_SUBST(LT_OBJDIR)
## Package-specific #defines
AC_DEFINE_UNQUOTED(RIVET_VERSION, "$PACKAGE_VERSION", "Rivet version string")
AC_DEFINE_UNQUOTED(RIVET_NAME, "$PACKAGE_NAME", "Rivet name string")
AC_DEFINE_UNQUOTED(RIVET_STRING, "$PACKAGE_STRING", "Rivet name and version string")
AC_DEFINE_UNQUOTED(RIVET_TARNAME, "$PACKAGE_TARNAME", "Rivet short name string")
AC_DEFINE_UNQUOTED(RIVET_BUGREPORT, "$PACKAGE_BUGREPORT", "Rivet contact email address")
## OS X
AC_CEDAR_OSX
## Work out the LCG platform tag
AC_LCG_TAG
## Set default compiler flags
if test "x$CXXFLAGS" == "x"; then CXXFLAGS="-O2"; fi
## Compiler setup
AC_LANG(C++)
AC_PROG_CXX
AX_CXX_COMPILE_STDCXX([11], [noext], [mandatory])
## Store and propagate the compiler identity and flags
RIVETCXX="$CXX"
AC_SUBST(RIVETCXX)
RIVETCXXFLAGS="$CXXFLAGS"
AC_SUBST(RIVETCXXFLAGS)
## Checks for programs.
AC_PROG_INSTALL
AC_PROG_LN_S
AC_DISABLE_STATIC
AC_LIBTOOL_DLOPEN
AC_PROG_LIBTOOL
-
+AX_EXECINFO
## YODA histogramming library
AC_CEDAR_LIBRARYANDHEADERS([YODA], , , [AC_MSG_ERROR([YODA is required])])
YODABINPATH=$YODALIBPATH/../bin
AC_SUBST(YODABINPATH)
AC_PATH_PROG(YODACONFIG, yoda-config, [], [$YODALIBPATH/../bin:$PATH])
YODA_PYTHONPATH=""
if test -f "$YODACONFIG"; then
AC_MSG_CHECKING([YODA version using yoda-config])
YODA_VERSION=`$YODACONFIG --version`
AC_MSG_RESULT([$YODA_VERSION])
YODA_VERSION1=[`echo $YODA_VERSION | cut -d. -f1 | sed -e 's/\([0-9]*\).*/\1/g'`]
YODA_VERSION2=[`echo $YODA_VERSION | cut -d. -f2 | sed -e 's/\([0-9]*\).*/\1/g'`]
YODA_VERSION3=[`echo $YODA_VERSION | cut -d. -f3 | sed -e 's/\([0-9]*\).*/\1/g'`]
let YODA_VERSION_INT=YODA_VERSION1*10000+YODA_VERSION2*100+YODA_VERSION3
if test $YODA_VERSION_INT -lt 10500; then
AC_MSG_ERROR([YODA version isn't sufficient: at least version 1.5.0 required])
fi
AC_MSG_CHECKING([YODA Python path using yoda-config])
YODA_PYTHONPATH=`$YODACONFIG --pythonpath`
AC_MSG_RESULT([$YODA_PYTHONPATH])
fi
AC_SUBST(YODA_PYTHONPATH)
## HepMC event record library
AC_CEDAR_LIBRARYANDHEADERS([HepMC], , , [AC_MSG_ERROR([HepMC is required])])
oldCPPFLAGS=$CPPFLAGS
CPPFLAGS="$CPPFLAGS -I$HEPMCINCPATH"
if test -e "$HEPMCINCPATH/HepMC/HepMCDefs.h"; then
AC_LANG_CONFTEST([AC_LANG_SOURCE([#include <iostream>
#include "HepMC/HepMCDefs.h"
int main() { std::cout << HEPMC_VERSION << std::endl; return 0; }])])
else
AC_LANG_CONFTEST([AC_LANG_SOURCE([#include <iostream>
#include "HepMC/defs.h"
int main() { std::cout << VERSION << std::endl; return 0; }])])
fi
if test -f conftest.cc; then
$CXX $CPPFLAGS conftest.cc -o conftest 2>&1 1>&5
elif test -f conftest.C; then
$CXX $CPPFLAGS conftest.C -o conftest 2>&1 1>&5
else
$CXX $CPPFLAGS conftest.cpp -o conftest 2>&1 1>&5
fi
hepmc_version=`./conftest`
if test x$hepmc_version != x; then
let hepmc_major=`echo "$hepmc_version" | cut -d. -f1`
let hepmc_minor=`echo "$hepmc_version" | cut -d. -f2`
fi
rm -f conftest conftest.cpp conftest.cc conftest.C
HEPMC_VERSION=$hepmc_major$hepmc_minor
AC_MSG_NOTICE([HepMC version is $hepmc_version -> $HEPMC_VERSION])
AC_SUBST(HEPMC_VERSION)
CPPFLAGS=$oldCPPFLAGS
## FastJet clustering library
AC_CEDAR_LIBRARYANDHEADERS([fastjet], , , [AC_MSG_ERROR([FastJet is required])])
AC_PATH_PROG(FJCONFIG, fastjet-config, [], $FASTJETPATH/bin:$PATH)
if test -f "$FJCONFIG"; then
AC_MSG_CHECKING([FastJet version using fastjet-config])
fjversion=`$FJCONFIG --version`
AC_MSG_RESULT([$fjversion])
fjmajor=$(echo $fjversion | cut -f1 -d.)
fjminor=$(echo $fjversion | cut -f2 -d.)
fjmicro=$(echo $fjversion | cut -f3 -d.)
if test "$fjmajor" -lt 3; then
AC_MSG_ERROR([FastJet version 3.0.0 or later is required])
fi
FASTJETCONFIGLIBADD="$($FJCONFIG --plugins --shared --libs)"
else
FASTJETCONFIGLIBADD="-L$FASTJETLIBPATH -l$FASTJETLIBNAME"
FASTJETCONFIGLIBADD="$FASTJETCONFIGLIBADD -lSISConePlugin -lsiscone -lsiscone_spherical"
FASTJETCONFIGLIBADD="$FASTJETCONFIGLIBADD -lCDFConesPlugin -lD0RunIIConePlugin -lNestedDefsPlugin"
FASTJETCONFIGLIBADD="$FASTJETCONFIGLIBADD -lTrackJetPlugin -lATLASConePlugin -lCMSIterativeConePlugin"
FASTJETCONFIGLIBADD="$FASTJETCONFIGLIBADD -lEECambridgePlugin -lJadePlugin"
fi;
AC_MSG_NOTICE([FastJet LIBADD = $FASTJETCONFIGLIBADD])
AC_SUBST(FASTJETCONFIGLIBADD)
# Check for FastJet headers that require the --enable-all(cxx)plugins option
FASTJET_ERRMSG="Required FastJet plugin headers were not found: did you build FastJet with the --enable-allcxxplugins option?"
oldCPPFLAGS=$CPPFLAGS
CPPFLAGS="$CPPFLAGS -I$FASTJETINCPATH"
AC_CHECK_HEADER([fastjet/D0RunIIConePlugin.hh], [], [AC_MSG_ERROR([$FASTJET_ERRMSG])])
AC_CHECK_HEADER([fastjet/TrackJetPlugin.hh], [], [AC_MSG_ERROR([$FASTJET_ERRMSG])])
CPPFLAGS=$oldCPPFLAGS
# ## GNU Scientific Library
# AC_SEARCH_GSL
# AC_CEDAR_HEADERS([gsl], , , [AC_MSG_ERROR([GSL (GNU Scientific Library) is required])])
# oldCPPFLAGS=$CPPFLAGS
# CPPFLAGS="$CPPFLAGS -I$GSLINCPATH"
# AC_CHECK_HEADER([gsl/gsl_vector.h], [], [AC_MSG_ERROR([GSL vectors not found.])])
# CPPFLAGS=$oldCPPFLAGS
## Disable build/install of standard analyses
AC_ARG_ENABLE([analyses],
[AC_HELP_STRING(--disable-analyses, [don't try to build or install standard analyses])],
[], [enable_analyses=yes])
if test x$enable_analyses != xyes; then
AC_MSG_WARN([Not building standard Rivet analyses, by request])
fi
AM_CONDITIONAL(ENABLE_ANALYSES, [test x$enable_analyses = xyes])
## Build LaTeX docs if possible...
AC_PATH_PROG(PDFLATEX, pdflatex)
AM_CONDITIONAL(WITH_PDFLATEX, [test x$PDFLATEX != x])
## ... unless told otherwise!
AC_ARG_ENABLE([pdfmanual],
[AC_HELP_STRING(--enable-pdfmanual, [build and install the manual])],
[], [enable_pdfmanual=no])
if test x$enable_pdfmanual = xyes; then
AC_MSG_WARN([Building Rivet manual, by request])
fi
AM_CONDITIONAL(ENABLE_PDFMANUAL, [test x$enable_pdfmanual = xyes])
## Build Doxygen documentation if possible
AC_ARG_ENABLE([doxygen],
[AC_HELP_STRING(--disable-doxygen, [don't try to make Doxygen documentation])],
[], [enable_doxygen=yes])
if test x$enable_doxygen = xyes; then
AC_PATH_PROG(DOXYGEN, doxygen)
fi
AM_CONDITIONAL(WITH_DOXYGEN, [test x$DOXYGEN != x])
## Build asciidoc docs if possible
AC_PATH_PROG(ASCIIDOC, asciidoc)
AM_CONDITIONAL(WITH_ASCIIDOC, [test x$ASCIIDOC != x])
## Python extension
AC_ARG_ENABLE(pyext, [AC_HELP_STRING(--disable-pyext,
[don't build Python module (default=build)])],
[], [enable_pyext=yes])
## Basic Python checks
if test x$enable_pyext == xyes; then
AX_PYTHON_DEVEL([>= '2.6'])
AC_SUBST(PYTHON_VERSION)
RIVET_PYTHONPATH=`$PYTHON -c "import distutils.sysconfig; print distutils.sysconfig.get_python_lib(prefix='$prefix', plat_specific=True);"`
AC_SUBST(RIVET_PYTHONPATH)
if test -z "$PYTHON"; then
AC_MSG_ERROR([Can't build Python extension since python can't be found])
enable_pyext=no
fi
if test -z "$PYTHON_CPPFLAGS"; then
AC_MSG_ERROR([Can't build Python extension since Python.h header file cannot be found])
enable_pyext=no
fi
fi
AM_CONDITIONAL(ENABLE_PYEXT, [test x$enable_pyext == xyes])
## Cython checks
if test x$enable_pyext == xyes; then
AM_CHECK_CYTHON([0.24.0], [:], [:])
if test x$CYTHON_FOUND = xyes; then
AC_MSG_NOTICE([Cython >= 0.24 found: Python extension source can be rebuilt (for developers)])
fi
AC_CHECK_FILE([pyext/rivet/core.cpp],
[],
[if test "x$CYTHON_FOUND" != "xyes"; then
AC_MSG_ERROR([Cython is required for --enable-pyext, no pre-built core.cpp was found.])
fi])
cython_compiler=$CXX
## Set extra Python extension build flags (to cope with Cython output code oddities)
PYEXT_CXXFLAGS="$CXXFLAGS"
AC_CEDAR_CHECKCXXFLAG([-Wno-unused-but-set-variable], [PYEXT_CXXFLAGS="$PYEXT_CXXFLAGS -Wno-unused-but-set-variable"])
AC_CEDAR_CHECKCXXFLAG([-Wno-sign-compare], [PYEXT_CXXFLAGS="$PYEXT_CXXFLAGS -Wno-sign-compare"])
AC_SUBST(PYEXT_CXXFLAGS)
AC_MSG_NOTICE([All Python build checks successful: 'rivet' Python extension will be built])
fi
AM_CONDITIONAL(WITH_CYTHON, [test x$CYTHON_FOUND = xyes])
## Set default build flags
AM_CPPFLAGS="-I\$(top_srcdir)/include -I\$(top_builddir)/include"
#AM_CPPFLAGS="$AM_CPPFLAGS -I\$(top_srcdir)/include/eigen3"
#AM_CPPFLAGS="$AM_CPPFLAGS \$(GSL_CPPFLAGS)"
dnl AM_CPPFLAGS="$AM_CPPFLAGS \$(BOOST_CPPFLAGS)"
AM_CPPFLAGS="$AM_CPPFLAGS -I\$(YODAINCPATH)"
AM_CPPFLAGS="$AM_CPPFLAGS -I\$(HEPMCINCPATH)"
AM_CPPFLAGS="$AM_CPPFLAGS -I\$(FASTJETINCPATH)"
AC_CEDAR_CHECKCXXFLAG([-pedantic], [AM_CXXFLAGS="$AM_CXXFLAGS -pedantic"])
AC_CEDAR_CHECKCXXFLAG([-Wall], [AM_CXXFLAGS="$AM_CXXFLAGS -Wall"])
AC_CEDAR_CHECKCXXFLAG([-Wno-long-long], [AM_CXXFLAGS="$AM_CXXFLAGS -Wno-long-long"])
AC_CEDAR_CHECKCXXFLAG([-Wno-format], [AM_CXXFLAGS="$AM_CXXFLAGS -Wno-format"])
dnl AC_CEDAR_CHECKCXXFLAG([-Wno-unused-variable], [AM_CXXFLAGS="$AM_CXXFLAGS -Wno-unused-variable"])
AC_CEDAR_CHECKCXXFLAG([-Werror=uninitialized], [AM_CXXFLAGS="$AM_CXXFLAGS -Werror=uninitialized"])
AC_CEDAR_CHECKCXXFLAG([-Werror=delete-non-virtual-dtor], [AM_CXXFLAGS="$AM_CXXFLAGS -Werror=delete-non-virtual-dtor"])
## Add OpenMP-enabling flags if possible
AX_OPENMP([AM_CXXFLAGS="$AM_CXXFLAGS $OPENMP_CXXFLAGS"])
## Debug flag (default=-DNDEBUG, enabled=-g)
AC_ARG_ENABLE([debug], [AC_HELP_STRING(--enable-debug,
[build with debugging symbols @<:@default=no@:>@])], [], [enable_debug=no])
if test x$enable_debug == xyes; then
AM_CXXFLAGS="$AM_CXXFLAGS -g"
fi
## Extra warnings flag (default=none)
AC_ARG_ENABLE([extra-warnings], [AC_HELP_STRING(--enable-extra-warnings,
[build with extra compiler warnings (recommended for developers) @<:@default=no@:>@])], [], [enable_extra_warnings=no])
if test x$enable_extra_warnings == xyes; then
AC_CEDAR_CHECKCXXFLAG([-Wextra], [AM_CXXFLAGS="$AM_CXXFLAGS -Wextra "])
fi
AC_SUBST(AM_CPPFLAGS)
AC_SUBST(AM_CXXFLAGS)
AC_EMPTY_SUBST
AC_CONFIG_FILES(Makefile Doxyfile)
AC_CONFIG_FILES(include/Makefile include/Rivet/Makefile)
AC_CONFIG_FILES(src/Makefile)
AC_CONFIG_FILES(src/Core/Makefile src/Core/yamlcpp/Makefile)
AC_CONFIG_FILES(src/Tools/Makefile)
AC_CONFIG_FILES(src/Projections/Makefile)
AC_CONFIG_FILES(src/AnalysisTools/Makefile)
AC_CONFIG_FILES(analyses/Makefile)
AC_CONFIG_FILES(test/Makefile)
AC_CONFIG_FILES(pyext/Makefile pyext/rivet/Makefile pyext/setup.py)
AC_CONFIG_FILES(data/Makefile data/texmf/Makefile)
AC_CONFIG_FILES(doc/Makefile)
AC_CONFIG_FILES(doc/rivetversion.sty doc/diffanas)
AC_CONFIG_FILES(bin/Makefile bin/rivet-config bin/rivet-buildplugin)
AC_CONFIG_FILES(rivetenv.sh rivetenv.csh rivet.pc)
AC_OUTPUT
if test x$enable_pyrivet == xyes; then
cat <<EOF
************************************************************
RIVET CONFIGURED!
Now build and install (to the $prefix tree) with e.g.
make -j2 && make -j2 install
To use Rivet, we recommend reading HepMC files from a file
or pipe (the latter may be made with mkfifo) using the
'rivet' executable.
For a more pleasant command line experience, you can include
the data/rivet-completion file into your .bashrc file,
or your bash_completion.d directory if you have one.
The rivetenv.*sh files will not be installed, but can help you
to set up a Rivet runtime environment in future.
************************************************************
EOF
fi
diff --git a/include/Rivet/Analyses/MC_JetAnalysis.hh b/include/Rivet/Analyses/MC_JetAnalysis.hh
--- a/include/Rivet/Analyses/MC_JetAnalysis.hh
+++ b/include/Rivet/Analyses/MC_JetAnalysis.hh
@@ -1,66 +1,67 @@
// -*- C++ -*-
#ifndef RIVET_MC_JetAnalysis_HH
#define RIVET_MC_JetAnalysis_HH
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
namespace Rivet {
/// @brief Base class providing common functionality for MC jet validation analyses
/// @todo Could reduce duplication by inheriting this from MC_ParticleAnalysis, with minor tweaking
class MC_JetAnalysis : public Analysis {
public:
/// Default constructor.
MC_JetAnalysis(const string& name,
size_t njet,
const string& jetpro_name,
double jetptcut=20*GeV);
/// @name Analysis methods
//@{
virtual void init();
virtual void analyze(const Event& event);
virtual void finalize();
//@}
protected:
/// The number of jets for which histograms are to be initialised
size_t _njet;
/// The name of the jet projection to be used for this analysis
/// (this projection has to be registered by the derived analysis!)
const std::string _jetpro_name;
/// Jet pT cutoff
double _jetptcut;
/// @todo Add jet masses and d(rap)
/// @name Histograms
//@{
std::vector<Histo1DPtr> _h_pT_jet;
std::vector<Histo1DPtr> _h_eta_jet;
std::vector<Histo1DPtr> _h_eta_jet_plus, _h_eta_jet_minus;
std::vector<Histo1DPtr> _h_rap_jet;
std::vector<Histo1DPtr> _h_rap_jet_plus, _h_rap_jet_minus;
std::vector<Histo1DPtr> _h_mass_jet;
+ std::vector<Scatter2DPtr> tmpeta, tmprap;
std::map<std::pair<size_t, size_t>, Histo1DPtr> _h_deta_jets;
std::map<std::pair<size_t, size_t>, Histo1DPtr> _h_dphi_jets;
std::map<std::pair<size_t, size_t>, Histo1DPtr> _h_dR_jets;
Histo1DPtr _h_jet_multi_exclusive;
Histo1DPtr _h_jet_multi_inclusive;
Scatter2DPtr _h_jet_multi_ratio;
Histo1DPtr _h_jet_HT, _h_mjj_jets;
//@}
};
}
#endif
diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh
--- a/include/Rivet/AnalysisHandler.hh
+++ b/include/Rivet/AnalysisHandler.hh
@@ -1,263 +1,271 @@
// -*- C++ -*-
#ifndef RIVET_RivetHandler_HH
#define RIVET_RivetHandler_HH
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Particle.hh"
#include "Rivet/AnalysisLoader.hh"
#include "Rivet/Tools/RivetYODA.hh"
namespace Rivet {
// Forward declaration and smart pointer for Analysis
class Analysis;
typedef std::shared_ptr<Analysis> AnaHandle;
// Needed to make smart pointers compare equivalent in the STL set
struct CmpAnaHandle {
bool operator() (const AnaHandle& a, const AnaHandle& b) {
return a.get() < b.get();
}
};
/// A class which handles a number of analysis objects to be applied to
/// generated events. An {@link Analysis}' AnalysisHandler is also responsible
/// for handling the final writing-out of histograms.
class AnalysisHandler {
public:
/// @name Constructors and destructors. */
//@{
/// Preferred constructor, with optional run name.
AnalysisHandler(const string& runname="");
/// @brief Destructor
/// The destructor is not virtual, as this class should not be inherited from.
~AnalysisHandler();
//@}
private:
/// Get a logger object.
Log& getLog() const;
public:
/// @name Run properties
//@{
/// Get the name of this run.
string runName() const;
/// Get the number of events seen. Should only really be used by external
/// steering code or analyses in the finalize phase.
size_t numEvents() const;
/// Get the sum of the event weights seen - the weighted equivalent of the
/// number of events. Should only really be used by external steering code
/// or analyses in the finalize phase.
double sumOfWeights() const {
return _eventCounter->sumW();
}
size_t numWeights() const {
return _weightNames.size();
}
/// Set the active weight.
void setActiveWeight(unsigned int iWeight);
/// Is cross-section information required by at least one child analysis?
bool needCrossSection() const;
/// Set the cross-section for the process being generated.
AnalysisHandler& setCrossSection(double xs);
/// Get the cross-section known to the handler.
double crossSection() const {
return _xs;
}
/// Whether the handler knows about a cross-section.
bool hasCrossSection() const;
/// Set the beam particles for this run
AnalysisHandler& setRunBeams(const ParticlePair& beams) {
_beams = beams;
MSG_DEBUG("Setting run beams = " << beams << " @ " << sqrtS()/GeV << " GeV");
return *this;
}
/// Get the beam particles for this run, usually determined from the first event.
const ParticlePair& beams() const { return _beams; }
/// Get beam IDs for this run, usually determined from the first event.
/// @deprecated Use standalone beamIds(ah.beams()), to clean AH interface
PdgIdPair beamIds() const;
/// Get energy for this run, usually determined from the first event.
/// @deprecated Use standalone sqrtS(ah.beams()), to clean AH interface
double sqrtS() const;
/// Setter for _ignoreBeams
void setIgnoreBeams(bool ignore=true);
//@}
/// @name Handle analyses
//@{
/// Get a list of the currently registered analyses' names.
std::vector<std::string> analysisNames() const;
/// Get the collection of currently registered analyses.
const std::set<AnaHandle, CmpAnaHandle>& analyses() const {
return _analyses;
}
/// Get a registered analysis by name.
const AnaHandle analysis(const std::string& analysisname) const;
/// Add an analysis to the run list by object
AnalysisHandler& addAnalysis(Analysis* analysis);
/// @brief Add an analysis to the run list using its name.
///
/// The actual Analysis to be used will be obtained via
/// AnalysisLoader::getAnalysis(string). If no matching analysis is found,
/// no analysis is added (i.e. the null pointer is checked and discarded.
AnalysisHandler& addAnalysis(const std::string& analysisname);
/// @brief Add analyses to the run list using their names.
///
/// The actual {@link Analysis}' to be used will be obtained via
/// AnalysisHandler::addAnalysis(string), which in turn uses
/// AnalysisLoader::getAnalysis(string). If no matching analysis is found
/// for a given name, no analysis is added, but also no error is thrown.
AnalysisHandler& addAnalyses(const std::vector<std::string>& analysisnames);
/// Remove an analysis from the run list using its name.
AnalysisHandler& removeAnalysis(const std::string& analysisname);
/// Remove analyses from the run list using their names.
AnalysisHandler& removeAnalyses(const std::vector<std::string>& analysisnames);
//@}
/// @name Main init/execute/finalise
//@{
/// Initialize a run, with the run beams taken from the example event.
void init(const GenEvent& event);
/// @brief Analyze the given \a event by reference.
///
/// This function will call the AnalysisBase::analyze() function of all
/// included analysis objects.
void analyze(const GenEvent& event);
/// @brief Analyze the given \a event by pointer.
///
/// This function will call the AnalysisBase::analyze() function of all
/// included analysis objects, after checking the event pointer validity.
void analyze(const GenEvent* event);
/// Finalize a run. This function calls the AnalysisBase::finalize()
/// functions of all included analysis objects.
void finalize();
//@}
/// @name Histogram / data object access
//@{
/// Add a vector of analysis objects to the current state.
void addData(const std::vector<YODA::AnalysisObjectPtr>& aos);
/// Read analysis plots into the histo collection (via addData) from the named file.
void readData(const std::string& filename);
/// Get all analyses' plots as a vector of analysis objects.
std::vector<YODA::AnalysisObjectPtr> getData() const;
std::vector<reference_wrapper<MultiweightAOPtr> > getRivetAOs() const;
std::vector<YODA::AnalysisObjectPtr> getYodaAOs() const;
/// Get all analyses' plots as a vector of analysis objects.
void setWeightNames(const GenEvent& ge);
/// Do we have named weights?
bool haveNamedWeights();
/// Write all analyses' plots (via getData) to the named file.
void writeData(const std::string& filename) const;
//@}
+ /// Indicate which Rivet stage we're in.
+ /// At the moment, only INIT is used to enable booking.
+ enum class Stage { OTHER, INIT };
+
+ /// Which stage are we in?
+ Stage stage() const { return _stage; }
private:
+ /// Current handler stage
+ Stage _stage = Stage::OTHER;
/// The collection of Analysis objects to be used.
set<AnaHandle, CmpAnaHandle> _analyses;
/// @name Run properties
//@{
/// Weight names
std::vector<std::string> _weightNames;
std::vector<std::valarray<double> > _subEventWeights;
size_t _numWeightTypes; // always == WeightVector.size()
/// Run name
std::string _runname;
mutable CounterPtr _eventCounter;
/// Cross-section known to AH
double _xs, _xserr;
/// Beams used by this run.
ParticlePair _beams;
/// Flag to check if init has been called
bool _initialised;
/// Flag whether input event beams should be ignored in compatibility check
bool _ignoreBeams;
/// Current event number
int _eventNumber;
//@}
private:
/// The assignment operator is private and must never be called.
/// In fact, it should not even be implemented.
AnalysisHandler& operator=(const AnalysisHandler&);
/// The copy constructor is private and must never be called. In
/// fact, it should not even be implemented.
AnalysisHandler(const AnalysisHandler&);
};
}
#endif
diff --git a/include/Rivet/Tools/RivetYODA.hh b/include/Rivet/Tools/RivetYODA.hh
--- a/include/Rivet/Tools/RivetYODA.hh
+++ b/include/Rivet/Tools/RivetYODA.hh
@@ -1,419 +1,438 @@
#ifndef RIVET_RIVETYODA_HH
#define RIVET_RIVETYODA_HH
/// @author Andy Buckley
/// @date 2009-01-30
/// @author David Grellscheid
/// @date 2011-07-18
/// @author David Grellscheid
/// @date 2016-09-27
#include "Rivet/Config/RivetCommon.hh"
#include "YODA/AnalysisObject.h"
#include "YODA/Counter.h"
#include "YODA/Histo1D.h"
#include "YODA/Histo2D.h"
#include "YODA/Profile1D.h"
#include "YODA/Profile2D.h"
#include "YODA/Scatter1D.h"
#include "YODA/Scatter2D.h"
#include "YODA/Scatter3D.h"
#include <map>
namespace YODA {
typedef std::shared_ptr<YODA::AnalysisObject> AnalysisObjectPtr;
// typedef std::shared_ptr<YODA::Scatter1D> Scatter1DPtr;
// typedef std::shared_ptr<YODA::Scatter2D> Scatter2DPtr;
// typedef std::shared_ptr<YODA::Scatter3D> Scatter3DPtr;
}
namespace Rivet {
class AnalysisObjectPtr {
public:
virtual ~AnalysisObjectPtr() {}
virtual YODA::AnalysisObject* operator->() = 0;
virtual YODA::AnalysisObject* operator->() const = 0;
virtual const YODA::AnalysisObject & operator*() const = 0;
/// @todo
/// rename to setActive(Idx)?
virtual void setActiveWeightIdx(unsigned int iWeight) = 0;
+ virtual void blockDestructor(bool) = 0;
+
bool operator ==(const AnalysisObjectPtr& p) { return (this == &p); }
protected:
/// @todo do we need this?
// virtual void reset() = 0;
};
/// @todo
/// implement scatter1dptr and scatter2dptr here
/// these need to be multi-weighted eventually.
/*
class Scatter1DPtr : public AnalysisObjectPtr {
public:
Scatter1DPtr() : _persistent() { }
Scatter1DPtr(size_t len_of_weightvec, const YODA::Scatter1D& p) {
for (size_t m = 0; m < len_of_weightvec; ++m)
_persistent.push_back(make_shared<YODA::Scatter1D>(p));
}
bool operator!() const { return !_persistent; }
operator bool() const { return bool(_persistent); }
YODA::Scatter1D* operator->() { return _persistent.get(); }
YODA::Scatter1D* operator->() const { return _persistent.get(); }
YODA::Scatter1D & operator*() { return *_persistent; }
const YODA::Scatter1D & operator*() const { return *_persistent; }
protected:
vector<YODA::Scatter1DPtr> _persistent;
};
class Scatter2DPtr : public AnalysisObjectPtr {
public:
Scatter2DPtr(size_t len_of_weightvec, const YODA::Scatter2D& p) {
for (size_t m = 0; m < len_of_weightvec; ++m)
_persistent.push_back(make_shared<YODA::Scatter2D>(p));
}
Scatter2DPtr() : _persistent() { }
bool operator!() { return !_persistent; }
operator bool() { return bool(_persistent); }
YODA::Scatter2D* operator->() { return _persistent.get(); }
YODA::Scatter2D* operator->() const { return _persistent.get(); }
YODA::Scatter2D & operator*() { return *_persistent; }
const YODA::Scatter2D & operator*() const { return *_persistent; }
protected:
vector<YODA::Scatter2DPtr> _persistent;
};
class Scatter3DPtr : public AnalysisObjectPtr {
public:
Scatter3DPtr(size_t len_of_weightvec, const YODA::Scatter3D& p) {
for (size_t m = 0; m < len_of_weightvec; ++m)
_persistent.push_back(make_shared<YODA::Scatter3D>(p));
}
Scatter3DPtr() : _persistent() { }
bool operator!() { return !_persistent; }
operator bool() { return bool(_persistent); }
YODA::Scatter3D* operator->() { return _persistent.get(); }
YODA::Scatter3D* operator->() const { return _persistent.get(); }
YODA::Scatter3D & operator*() { return *_persistent; }
const YODA::Scatter3D & operator*() const { return *_persistent; }
protected:
vector<YODA::Scatter3DPtr> _persistent;
};
*/
class MultiweightAOPtr : public AnalysisObjectPtr {
public:
virtual void newSubEvent() = 0;
virtual void pushToPersistent(const vector<valarray<double> >& weight) = 0;
virtual YODA::AnalysisObjectPtr activeYODAPtr() const = 0;
};
using Weight = double;
template <class T>
using Fill = pair<typename T::FillType, Weight>;
template <class T>
using Fills = multiset<Fill<T>>;
// TODO TODO TODO
// need to override the old fill method too!
// otherwise we can't intercept existing fill calls in analysis code
// TODO TODO TODO
template <class T>
class TupleWrapper;
template<>
class TupleWrapper<YODA::Counter> : public YODA::Counter {
public:
typedef shared_ptr<TupleWrapper<YODA::Counter>> Ptr;
TupleWrapper(const YODA::Counter & h) : YODA::Counter(h) {}
// todo: do we need to deal with users using fractions directly?
void fill( double weight=1.0, double fraction=1.0 ) {
fills_.insert( {YODA::Counter::FillType(),weight} );
}
void reset() { fills_.clear(); }
const Fills<YODA::Counter> & fills() const { return fills_; }
private:
// x / weight pairs
Fills<YODA::Counter> fills_;
};
template<>
class TupleWrapper<YODA::Histo1D> : public YODA::Histo1D {
public:
typedef shared_ptr<TupleWrapper<YODA::Histo1D>> Ptr;
TupleWrapper(const YODA::Histo1D & h) : YODA::Histo1D(h) {}
// todo: do we need to deal with users using fractions directly?
void fill( double x, double weight=1.0, double fraction=1.0 ) {
if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
fills_.insert( { x , weight } );
}
void reset() { fills_.clear(); }
const Fills<YODA::Histo1D> & fills() const { return fills_; }
private:
// x / weight pairs
Fills<YODA::Histo1D> fills_;
};
template<>
class TupleWrapper<YODA::Profile1D> : public YODA::Profile1D {
public:
typedef shared_ptr<TupleWrapper<YODA::Profile1D>> Ptr;
TupleWrapper(const YODA::Profile1D & h) : YODA::Profile1D(h) {}
// todo: do we need to deal with users using fractions directly?
void fill( double x, double y, double weight=1.0, double fraction=1.0 ) {
if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
fills_.insert( { YODA::Profile1D::FillType{x,y}, weight } );
}
void reset() { fills_.clear(); }
const Fills<YODA::Profile1D> & fills() const { return fills_; }
private:
// x / weight pairs
Fills<YODA::Profile1D> fills_;
};
template<>
class TupleWrapper<YODA::Histo2D> : public YODA::Histo2D {
public:
typedef shared_ptr<TupleWrapper<YODA::Histo2D>> Ptr;
TupleWrapper(const YODA::Histo2D & h) : YODA::Histo2D(h) {}
// todo: do we need to deal with users using fractions directly?
void fill( double x, double y, double weight=1.0, double fraction=1.0 ) {
if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
fills_.insert( { YODA::Histo2D::FillType{x,y}, weight } );
}
void reset() { fills_.clear(); }
const Fills<YODA::Histo2D> & fills() const { return fills_; }
private:
// x / weight pairs
Fills<YODA::Histo2D> fills_;
};
template<>
class TupleWrapper<YODA::Profile2D> : public YODA::Profile2D {
public:
typedef shared_ptr<TupleWrapper<YODA::Profile2D>> Ptr;
TupleWrapper(const YODA::Profile2D & h) : YODA::Profile2D(h) {}
// todo: do we need to deal with users using fractions directly?
void fill( double x, double y, double z, double weight=1.0, double fraction=1.0 ) {
if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
if ( std::isnan(z) ) throw YODA::RangeError("Z is NaN");
fills_.insert( { YODA::Profile2D::FillType{x,y,z}, weight } );
}
void reset() { fills_.clear(); }
const Fills<YODA::Profile2D> & fills() const { return fills_; }
private:
// x / weight pairs
Fills<YODA::Profile2D> fills_;
};
template<>
class TupleWrapper<YODA::Scatter1D> : public YODA::Scatter1D {
public:
typedef shared_ptr<TupleWrapper<YODA::Scatter1D>> Ptr;
TupleWrapper(const YODA::Scatter1D & h) : YODA::Scatter1D(h) {}
};
template<>
class TupleWrapper<YODA::Scatter2D> : public YODA::Scatter2D {
public:
typedef shared_ptr<TupleWrapper<YODA::Scatter2D>> Ptr;
TupleWrapper(const YODA::Scatter2D & h) : YODA::Scatter2D(h) {}
};
template<>
class TupleWrapper<YODA::Scatter3D> : public YODA::Scatter3D {
public:
typedef shared_ptr<TupleWrapper<YODA::Scatter3D>> Ptr;
TupleWrapper(const YODA::Scatter3D & h) : YODA::Scatter3D(h) {}
};
template <class T>
class Wrapper : public MultiweightAOPtr {
+ friend class Analysis;
+
public:
/* @todo
* some things are not really well-defined here
* for instance: fill() in the finalize() method and integral() in
* the analyze() method.
*/
Wrapper() = default;
Wrapper(size_t len_of_weightvec, const T & p);
+ ~Wrapper();
+
typename T::Ptr active() const;
/* @todo this probably need to loop over all? */
bool operator!() const { return !_active; } // Don't use active() here, assert will catch
operator bool() const { return static_cast<bool>(_active); } // Don't use active() here, assert will catch
T * operator->() { return active().get(); }
T * operator->() const { return active().get(); }
T & operator*() { return *active(); }
const T & operator*() const { return *active(); }
/* @todo
* these need to be re-thought out.
void reset() { active()->reset(); }
*/
/* @todo
* these probably need to loop over all?
* do we even want to provide equality?
*/
/* @todo
* how about no.
friend bool operator==(Wrapper a, Wrapper b){
if (a._persistent.size() != b._persistent.size())
return false;
for (size_t i = 0; i < a._persistent.size(); i++) {
if (a._persistent.at(i) != b._persistent.at(i)) {
return false;
}
}
return true;
}
friend bool operator!=(Wrapper a, Wrapper b){
return !(a == b);
}
friend bool operator<(Wrapper a, Wrapper b){
if (a._persistent.size() >= b._persistent.size())
return false;
for (size_t i = 0; i < a._persistent.size(); i++) {
if (*(a._persistent.at(i)) >= *(b._persistent.at(i))) {
return false;
}
}
return true;
}
*/
private:
void setActiveWeightIdx(unsigned int iWeight) {
_active = _persistent.at(iWeight);
}
/* this is for dev only---we shouldn't need this in real runs. */
void unsetActiveWeight() { _active.reset(); }
void newSubEvent();
virtual YODA::AnalysisObjectPtr activeYODAPtr() const { return _active; }
const vector<typename T::Ptr> & persistent() const { return _persistent; }
/* to be implemented for each type */
void pushToPersistent(const vector<valarray<double> >& weight);
+ /// Set destructor blocking flag.
+ /// Once booked in an analysis, we need to block the destructor
+ /// from being called until we're done with finalize.
+ /// This gives a clearer indication to the users that they're booking with
+ /// temporary objects.
+ void blockDestructor(bool b) { _blockDestructor = b; }
+
/* M of these, one for each weight */
vector<typename T::Ptr> _persistent;
/* N of these, one for each event in evgroup */
vector<typename TupleWrapper<T>::Ptr> _evgroup;
typename T::Ptr _active;
+ /// Destructor blocking flag.
+ /// Once booked in an analysis, we need to block the destructor
+ /// from being called until we're done with finalize.
+ ///
+ bool _blockDestructor = false;
+
// do we need implicit cast?
// operator typename T::Ptr () {
// return _active;
// }
friend class AnalysisHandler;
};
// every object listed here needs a virtual fill method in YODA,
// otherwise the Tuple fakery won't work.
- using Histo1DPtr = Wrapper<YODA::Histo1D>;
- using Histo2DPtr = Wrapper<YODA::Histo2D>;
+ using Histo1DPtr = Wrapper<YODA::Histo1D>;
+ using Histo2DPtr = Wrapper<YODA::Histo2D>;
using Profile1DPtr = Wrapper<YODA::Profile1D>;
using Profile2DPtr = Wrapper<YODA::Profile2D>;
- using CounterPtr = Wrapper<YODA::Counter>;
+ using CounterPtr = Wrapper<YODA::Counter>;
using Scatter1DPtr = Wrapper<YODA::Scatter1D>;
using Scatter2DPtr = Wrapper<YODA::Scatter2D>;
using Scatter3DPtr = Wrapper<YODA::Scatter3D>;
using YODA::Counter;
using YODA::Histo1D;
using YODA::HistoBin1D;
using YODA::Histo2D;
using YODA::HistoBin2D;
using YODA::Profile1D;
using YODA::ProfileBin1D;
using YODA::Profile2D;
using YODA::ProfileBin2D;
using YODA::Scatter1D;
using YODA::Point1D;
using YODA::Scatter2D;
using YODA::Point2D;
using YODA::Scatter3D;
using YODA::Point3D;
/// Function to get a map of all the refdata in a paper with the
/// given @a papername.
map<string, YODA::AnalysisObjectPtr> getRefData(const string& papername);
/// @todo Also provide a Scatter3D getRefData() version?
/// Get the file system path to the reference file for this paper.
string getDatafilePath(const string& papername);
}
#endif
diff --git a/m4/ax_execinfo.m4 b/m4/ax_execinfo.m4
new file mode 100644
--- /dev/null
+++ b/m4/ax_execinfo.m4
@@ -0,0 +1,67 @@
+# ===========================================================================
+# https://www.gnu.org/software/autoconf-archive/ax_execinfo.html
+# ===========================================================================
+#
+# SYNOPSIS
+#
+# AX_EXECINFO([ACTION-IF-EXECINFO-H-IS-FOUND], [ACTION-IF-EXECINFO-H-IS-NOT-FOUND], [ADDITIONAL-TYPES-LIST])
+#
+# DESCRIPTION
+#
+# Checks for execinfo.h header and if the len parameter/return type can be
+# found from a list, also define backtrace_size_t to that type.
+#
+# By default the list of types to try contains int and size_t, but should
+# some yet undiscovered system use e.g. unsigned, the 3rd argument can be
+# used for extensions. I'd like to hear of further suggestions.
+#
+# Executes ACTION-IF-EXECINFO-H-IS-FOUND when present and the execinfo.h
+# header is found or ACTION-IF-EXECINFO-H-IS-NOT-FOUND in case the header
+# seems unavailable.
+#
+# Also adds -lexecinfo to LIBS on BSD if needed.
+#
+# LICENSE
+#
+# Copyright (c) 2014 Thomas Jahns <jahns@dkrz.de>
+#
+# Copying and distribution of this file, with or without modification, are
+# permitted in any medium without royalty provided the copyright notice
+# and this notice are preserved. This file is offered as-is, without any
+# warranty.
+
+#serial 2
+
+AC_DEFUN([AX_EXECINFO],
+ [AC_CHECK_HEADERS([execinfo.h])
+ AS_IF([test x"$ac_cv_header_execinfo_h" = xyes],
+ [AC_CACHE_CHECK([size parameter type for backtrace()],
+ [ax_cv_proto_backtrace_type],
+ [AC_LANG_PUSH([C])
+ for ax_cv_proto_backtrace_type in size_t int m4_ifnblank([$3],[$3 ])none; do
+ AS_IF([test "${ax_cv_proto_backtrace_type}" = none],
+ [ax_cv_proto_backtrace_type= ; break])
+ AC_COMPILE_IFELSE([AC_LANG_PROGRAM([
+#include <execinfo.h>
+extern
+${ax_cv_proto_backtrace_type} backtrace(void **addrlist, ${ax_cv_proto_backtrace_type} len);
+char **backtrace_symbols(void *const *buffer, ${ax_cv_proto_backtrace_type} size);
+])],
+ [break])
+ done
+ AC_LANG_POP([C])])])
+ AS_IF([test x${ax_cv_proto_backtrace_type} != x],
+ [AC_DEFINE_UNQUOTED([backtrace_size_t], [$ax_cv_proto_backtrace_type],
+ [Defined to return type of backtrace().])])
+ AC_SEARCH_LIBS([backtrace],[execinfo])
+ AS_IF([test x"${ax_cv_proto_backtrace_type}" != x -a x"$ac_cv_header_execinfo_h" = xyes -a x"$ac_cv_search_backtrace" != xno],
+ [AC_DEFINE([HAVE_BACKTRACE],[1],
+ [Defined if backtrace() could be fully identified.])
+ ]m4_ifnblank([$1],[$1
+]),m4_ifnblank([$2],[$2
+]))])
+dnl
+dnl Local Variables:
+dnl mode: autoconf
+dnl End:
+dnl
diff --git a/src/AnalysisTools/MC_JetAnalysis.cc b/src/AnalysisTools/MC_JetAnalysis.cc
--- a/src/AnalysisTools/MC_JetAnalysis.cc
+++ b/src/AnalysisTools/MC_JetAnalysis.cc
@@ -1,190 +1,185 @@
// -*- C++ -*-
#include "Rivet/Analyses/MC_JetAnalysis.hh"
#include "Rivet/Projections/FastJets.hh"
namespace Rivet {
MC_JetAnalysis::MC_JetAnalysis(const string& name,
size_t njet,
const string& jetpro_name,
double jetptcut)
: Analysis(name), _njet(njet), _jetpro_name(jetpro_name), _jetptcut(jetptcut),
_h_pT_jet(njet),
_h_eta_jet(njet), _h_eta_jet_plus(njet), _h_eta_jet_minus(njet),
_h_rap_jet(njet), _h_rap_jet_plus(njet), _h_rap_jet_minus(njet),
- _h_mass_jet(njet)
+ _h_mass_jet(njet), tmpeta(njet), tmprap(njet)
{
setNeedsCrossSection(true); // legitimate use, since a base class has no .info file!
}
// Book histograms
void MC_JetAnalysis::init() {
const double sqrts = sqrtS() ? sqrtS() : 14000.*GeV;
for (size_t i = 0; i < _njet; ++i) {
const string pTname = "jet_pT_" + to_str(i+1);
const double pTmax = 1.0/(double(i)+2.0) * sqrts/GeV/2.0;
const int nbins_pT = 100/(i+1);
if (pTmax > 10) { // Protection aginst logspace exception, needed for LEP
book(_h_pT_jet[i] ,pTname, logspace(nbins_pT, 10.0, pTmax));
}
const string massname = "jet_mass_" + to_str(i+1);
const double mmax = 100.0;
const int nbins_m = 100/(i+1);
book(_h_mass_jet[i] ,massname, logspace(nbins_m, 1.0, mmax));
const string etaname = "jet_eta_" + to_str(i+1);
book(_h_eta_jet[i] ,etaname, i > 1 ? 25 : 50, -5.0, 5.0);
book(_h_eta_jet_plus[i], "_" + etaname + "_plus", i > 1 ? 15 : 25, 0, 5);
book(_h_eta_jet_minus[i], "_" + etaname + "_minus", i > 1 ? 15 : 25, 0, 5);
const string rapname = "jet_y_" + to_str(i+1);
book(_h_rap_jet[i] ,rapname, i>1 ? 25 : 50, -5.0, 5.0);
book(_h_rap_jet_plus[i], "_" + rapname + "_plus", i > 1 ? 15 : 25, 0, 5);
book(_h_rap_jet_minus[i], "_" + rapname + "_minus", i > 1 ? 15 : 25, 0, 5);
+ book(tmpeta[i], "jet_eta_pmratio_" + to_str(i+1));
+ book(tmprap[i], "jet_y_pmratio_" + to_str(i+1));
+
for (size_t j = i+1; j < min(size_t(3), _njet); ++j) {
const std::pair<size_t, size_t> ij = std::make_pair(i, j);
+ const string ijstr = to_str(i+1) + to_str(j+1);
- string detaname = "jets_deta_" + to_str(i+1) + to_str(j+1);
- Histo1DPtr tmpeta;
- book(tmpeta, detaname, 25, -5.0, 5.0);
- _h_deta_jets.insert(make_pair(ij, tmpeta));
+ string detaname = "jets_deta_" + ijstr;
+ book(_h_deta_jets[ij], detaname, 25, -5.0, 5.0);
- string dphiname = "jets_dphi_" + to_str(i+1) + to_str(j+1);
- Histo1DPtr tmpphi;
- book(tmpphi, dphiname, 25, 0.0, M_PI);
- _h_dphi_jets.insert(make_pair(ij, tmpphi));
+ string dphiname = "jets_dphi_" + ijstr;
+ book(_h_dphi_jets[ij], dphiname, 25, 0.0, M_PI);
- string dRname = "jets_dR_" + to_str(i+1) + to_str(j+1);
- Histo1DPtr tmpR;
- book(tmpR, dRname, 25, 0.0, 5.0);
- _h_dR_jets.insert(make_pair(ij, tmpR));
+ string dRname = "jets_dR_" + ijstr;
+ book(_h_dR_jets[ij], dRname, 25, 0.0, 5.0);
}
}
book(_h_jet_multi_exclusive ,"jet_multi_exclusive", _njet+3, -0.5, _njet+3-0.5);
book(_h_jet_multi_inclusive ,"jet_multi_inclusive", _njet+3, -0.5, _njet+3-0.5);
book(_h_jet_multi_ratio, "jet_multi_ratio");
book(_h_jet_HT ,"jet_HT", logspace(50, _jetptcut, sqrts/GeV/2.0));
book(_h_mjj_jets, "jets_mjj", 40, 0.0, sqrts/GeV/2.0);
}
// Do the analysis
void MC_JetAnalysis::analyze(const Event & e) {
const Jets& jets = apply<FastJets>(e, _jetpro_name).jetsByPt(_jetptcut);
for (size_t i = 0; i < _njet; ++i) {
if (jets.size() < i+1) continue;
_h_pT_jet[i]->fill(jets[i].pT()/GeV);
// Check for numerical precision issues with jet masses
double m2_i = jets[i].mass2();
if (m2_i < 0) {
if (m2_i < -1e-4) {
MSG_WARNING("Jet mass2 is negative: " << m2_i << " GeV^2. "
<< "Truncating to 0.0, assuming numerical precision is to blame.");
}
m2_i = 0.0;
}
// Jet mass
_h_mass_jet[i]->fill(sqrt(m2_i)/GeV);
// Jet eta
const double eta_i = jets[i].eta();
_h_eta_jet[i]->fill(eta_i);
(eta_i > 0.0 ? _h_eta_jet_plus : _h_eta_jet_minus)[i]->fill(fabs(eta_i));
// Jet rapidity
const double rap_i = jets[i].rapidity();
_h_rap_jet[i]->fill(rap_i);
(rap_i > 0.0 ? _h_rap_jet_plus : _h_rap_jet_minus)[i]->fill(fabs(rap_i));
// Inter-jet properties
for (size_t j = i+1; j < min(size_t(3),_njet); ++j) {
if (jets.size() < j+1) continue;
std::pair<size_t, size_t> ij = std::make_pair(i, j);
double deta = jets[i].eta()-jets[j].eta();
double dphi = deltaPhi(jets[i].momentum(),jets[j].momentum());
double dR = deltaR(jets[i].momentum(), jets[j].momentum());
_h_deta_jets[ij]->fill(deta);
_h_dphi_jets[ij]->fill(dphi);
_h_dR_jets[ij]->fill(dR);
}
}
// Multiplicities
_h_jet_multi_exclusive->fill(jets.size());
for (size_t i = 0; i < _njet+2; ++i) {
if (jets.size() >= i) {
_h_jet_multi_inclusive->fill(i);
}
}
// HT
double HT = 0.0;
foreach (const Jet& jet, jets) {
HT += jet.pT();
}
_h_jet_HT->fill(HT);
// mjj
if (jets.size() > 1) {
double mjj = (jets[0].momentum() + jets[1].momentum()).mass();
_h_mjj_jets->fill(mjj);
}
}
// Finalize
void MC_JetAnalysis::finalize() {
const double scaling = crossSection()/sumOfWeights();
for (size_t i = 0; i < _njet; ++i) {
scale(_h_pT_jet[i], scaling);
scale(_h_mass_jet[i], scaling);
scale(_h_eta_jet[i], scaling);
scale(_h_rap_jet[i], scaling);
// Create eta/rapidity ratio plots
- Scatter2DPtr tmpeta, tmprap;
- book(tmpeta, "jet_eta_pmratio_" + to_str(i+1));
- book(tmprap, "jet_y_pmratio_" + to_str(i+1));
- divide(*_h_eta_jet_plus[i], *_h_eta_jet_minus[i], tmpeta);
- divide(*_h_rap_jet_plus[i], *_h_rap_jet_minus[i], tmprap);
+ divide(*_h_eta_jet_plus[i], *_h_eta_jet_minus[i], tmpeta[i]);
+ divide(*_h_rap_jet_plus[i], *_h_rap_jet_minus[i], tmprap[i]);
}
// Scale the d{eta,phi,R} histograms
typedef map<pair<size_t, size_t>, Histo1DPtr> HistMap;
for (HistMap::value_type& it : _h_deta_jets) scale(it.second, scaling);
for (HistMap::value_type& it : _h_dphi_jets) scale(it.second, scaling);
for (HistMap::value_type& it : _h_dR_jets) scale(it.second, scaling);
// Fill inclusive jet multi ratio
int Nbins = _h_jet_multi_inclusive->numBins();
for (int i = 0; i < Nbins-1; ++i) {
_h_jet_multi_ratio->addPoint(i+1, 0, 0.5, 0);
if (_h_jet_multi_inclusive->bin(i).sumW() > 0.0) {
const double ratio = _h_jet_multi_inclusive->bin(i+1).sumW()/_h_jet_multi_inclusive->bin(i).sumW();
const double relerr_i = _h_jet_multi_inclusive->bin(i).relErr();
const double relerr_j = _h_jet_multi_inclusive->bin(i+1).relErr();
const double err = ratio * (relerr_i + relerr_j);
_h_jet_multi_ratio->point(i).setY(ratio, err);
}
}
scale(_h_jet_multi_exclusive, scaling);
scale(_h_jet_multi_inclusive, scaling);
scale(_h_jet_HT, scaling);
scale(_h_mjj_jets, scaling);
}
}
diff --git a/src/Core/Analysis.cc b/src/Core/Analysis.cc
--- a/src/Core/Analysis.cc
+++ b/src/Core/Analysis.cc
@@ -1,885 +1,926 @@
// -*- C++ -*-
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Analysis.hh"
#include "Rivet/AnalysisHandler.hh"
#include "Rivet/AnalysisInfo.hh"
#include "Rivet/Tools/BeamConstraint.hh"
+#include "DummyConfig.hh"
+#ifdef HAVE_EXECINFO_H
+#include <execinfo.h>
+#endif
+
namespace Rivet {
Analysis::Analysis(const string& name)
: _crossSection(-1.0),
_gotCrossSection(false),
_analysishandler(NULL)
{
ProjectionApplier::_allowProjReg = false;
_defaultname = name;
unique_ptr<AnalysisInfo> ai = AnalysisInfo::make(name);
assert(ai);
_info = move(ai);
assert(_info);
}
double Analysis::sqrtS() const {
return handler().sqrtS();
}
const ParticlePair& Analysis::beams() const {
return handler().beams();
}
const PdgIdPair Analysis::beamIds() const {
return handler().beamIds();
}
const string Analysis::histoDir() const {
/// @todo Cache in a member variable
string _histoDir;
if (_histoDir.empty()) {
_histoDir = "/" + name();
if (handler().runName().length() > 0) {
_histoDir = "/" + handler().runName() + _histoDir;
}
replace_all(_histoDir, "//", "/"); //< iterates until none
}
return _histoDir;
}
const string Analysis::histoPath(const string& hname) const {
const string path = histoDir() + "/" + hname;
return path;
}
const string Analysis::histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
return histoDir() + "/" + makeAxisCode(datasetId, xAxisId, yAxisId);
}
const string Analysis::makeAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
stringstream axisCode;
axisCode << "d";
if (datasetId < 10) axisCode << 0;
axisCode << datasetId;
axisCode << "-x";
if (xAxisId < 10) axisCode << 0;
axisCode << xAxisId;
axisCode << "-y";
if (yAxisId < 10) axisCode << 0;
axisCode << yAxisId;
return axisCode.str();
}
Log& Analysis::getLog() const {
string logname = "Rivet.Analysis." + name();
return Log::getLog(logname);
}
size_t Analysis::numEvents() const {
return handler().numEvents();
}
double Analysis::sumOfWeights() const {
return handler().sumOfWeights();
}
///////////////////////////////////////////
bool Analysis::isCompatible(const ParticlePair& beams) const {
return isCompatible(beams.first.pid(), beams.second.pid(),
beams.first.energy(), beams.second.energy());
}
bool Analysis::isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const {
PdgIdPair beams(beam1, beam2);
pair<double,double> energies(e1, e2);
return isCompatible(beams, energies);
}
bool Analysis::isCompatible(const PdgIdPair& beams, const pair<double,double>& energies) const {
// First check the beam IDs
bool beamIdsOk = false;
foreach (const PdgIdPair& bp, requiredBeams()) {
if (compatible(beams, bp)) {
beamIdsOk = true;
break;
}
}
if (!beamIdsOk) return false;
// Next check that the energies are compatible (within 1% or 1 GeV, whichever is larger, for a bit of UI forgiveness)
/// @todo Use some sort of standard ordering to improve comparisons, esp. when the two beams are different particles
bool beamEnergiesOk = requiredEnergies().size() > 0 ? false : true;
typedef pair<double,double> DoublePair;
foreach (const DoublePair& ep, requiredEnergies()) {
if ((fuzzyEquals(ep.first, energies.first, 0.01) && fuzzyEquals(ep.second, energies.second, 0.01)) ||
(fuzzyEquals(ep.first, energies.second, 0.01) && fuzzyEquals(ep.second, energies.first, 0.01)) ||
(abs(ep.first - energies.first) < 1*GeV && abs(ep.second - energies.second) < 1*GeV) ||
(abs(ep.first - energies.second) < 1*GeV && abs(ep.second - energies.first) < 1*GeV)) {
beamEnergiesOk = true;
break;
}
}
return beamEnergiesOk;
/// @todo Need to also check internal consistency of the analysis'
/// beam requirements with those of the projections it uses.
}
///////////////////////////////////////////
Analysis& Analysis::setCrossSection(double xs) {
_crossSection = xs;
_gotCrossSection = true;
return *this;
}
double Analysis::crossSection() const {
if (!_gotCrossSection || std::isnan(_crossSection)) {
string errMsg = "You did not set the cross section for the analysis " + name();
throw Error(errMsg);
}
return _crossSection;
}
double Analysis::crossSectionPerEvent() const {
return _crossSection/sumOfWeights();
}
////////////////////////////////////////////////////////////
// Histogramming
void Analysis::_cacheRefData() const {
if (_refdata.empty()) {
MSG_TRACE("Getting refdata cache for paper " << name());
_refdata = getRefData(name());
}
}
CounterPtr & Analysis::book(CounterPtr & ctr,
const string& cname,
const string& title) {
const string path = histoPath(cname);
ctr = CounterPtr(handler().numWeights(), Counter(path, title));
addAnalysisObject(ctr);
MSG_TRACE("Made counter " << cname << " for " << name());
return ctr;
}
CounterPtr & Analysis::book(CounterPtr & ctr, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
const string& title) {
// const string& xtitle,
// const string& ytitle) {
const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return book(ctr, axisCode, title);
}
Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
size_t nbins, double lower, double upper,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
Histo1D hist = Histo1D(nbins, lower, upper, path, title);
hist.setAnnotation("XLabel", xtitle);
hist.setAnnotation("YLabel", ytitle);
histo = Histo1DPtr(handler().numWeights(), hist);
addAnalysisObject(histo);
MSG_TRACE("Made histogram " << hname << " for " << name());
return histo;
// Histo1DPtr hist;
// try { // try to bind to pre-existing
// // AnalysisObjectPtr ao = getAnalysisObject(path);
// // hist = dynamic_pointer_cast<Histo1D>(ao);
// hist = getHisto1D(hname);
// /// @todo Test that cast worked
// /// @todo Also test that binning is as expected?
// MSG_TRACE("Bound pre-existing histogram " << hname << " for " << name());
// } catch (...) { // binding failed; make it from scratch
// hist = make_shared<Histo1D>(nbins, lower, upper, histoPath(hname), title);
// addAnalysisObject(hist);
// MSG_TRACE("Made histogram " << hname << " for " << name());
// }
}
Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
const initializer_list<double>& binedges,
const string& title,
const string& xtitle,
const string& ytitle) {
return book(histo, hname, vector<double>{binedges}, title, xtitle, ytitle);
}
Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
const vector<double>& binedges,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
// Histo1DPtr hist;
// try { // try to bind to pre-existing
// // AnalysisObjectPtr ao = getAnalysisObject(path);
// // hist = dynamic_pointer_cast<Histo1D>(ao);
// hist = getHisto1D(hname);
// /// @todo Test that cast worked
// /// @todo Also test that binning is as expected?
// MSG_TRACE("Bound pre-existing histogram " << hname << " for " << name());
// } catch (...) { // binding failed; make it from scratch
// hist = make_shared<Histo1D>(binedges, histoPath(hname), title);
// addAnalysisObject(hist);
// MSG_TRACE("Made histogram " << hname << " for " << name());
// }
Histo1D hist = Histo1D(binedges, path, title);
hist.setAnnotation("XLabel", xtitle);
hist.setAnnotation("YLabel", ytitle);
histo = Histo1DPtr(handler().numWeights(), hist);
addAnalysisObject(histo);
MSG_TRACE("Made histogram " << hname << " for " << name());
return histo;
}
Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
const Scatter2D& refscatter,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
// Histo1DPtr hist;
// try { // try to bind to pre-existing
// // AnalysisObjectPtr ao = getAnalysisObject(path);
// // hist = dynamic_pointer_cast<Histo1D>(ao);
// hist = getHisto1D(hname);
// /// @todo Test that cast worked
// /// @todo Also test that binning is as expected?
// MSG_TRACE("Bound pre-existing histogram " << hname << " for " << name());
// } catch (...) { // binding failed; make it from scratch
// hist = make_shared<Histo1D>(refscatter, histoPath(hname));
// if (hist->hasAnnotation("IsRef")) hist->rmAnnotation("IsRef");
// addAnalysisObject(hist);
// MSG_TRACE("Made histogram " << hname << " for " << name());
// }
Histo1D hist = Histo1D(refscatter, path);
hist.setTitle(title);
hist.setAnnotation("XLabel", xtitle);
hist.setAnnotation("YLabel", ytitle);
histo = Histo1DPtr(handler().numWeights(), hist);
addAnalysisObject(histo);
MSG_TRACE("Made histogram " << hname << " for " << name());
return histo;
}
Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
const string& title,
const string& xtitle,
const string& ytitle) {
const Scatter2D& refdata = refData(hname);
return book(histo, hname, refdata, title, xtitle, ytitle);
}
Histo1DPtr & Analysis::book(Histo1DPtr & histo, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
const string& title,
const string& xtitle,
const string& ytitle) {
const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return book(histo, axisCode, title, xtitle, ytitle);
}
/// @todo Add booking methods which take a path, titles and *a reference Scatter from which to book*
/////////////////
Histo2DPtr & Analysis::book(Histo2DPtr & h2d,const string& hname,
size_t nxbins, double xlower, double xupper,
size_t nybins, double ylower, double yupper,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
const string path = histoPath(hname);
Histo2D hist(nxbins, xlower, xupper, nybins, ylower, yupper, path, title);
hist.setAnnotation("XLabel", xtitle);
hist.setAnnotation("YLabel", ytitle);
hist.setAnnotation("ZLabel", ztitle);
h2d = Histo2DPtr(handler().numWeights(), hist);
addAnalysisObject(h2d);
MSG_TRACE("Made 2D histogram " << hname << " for " << name());
return h2d;
}
Histo2DPtr & Analysis::book(Histo2DPtr & h2d,const string& hname,
const initializer_list<double>& xbinedges,
const initializer_list<double>& ybinedges,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
return book(h2d, hname, vector<double>{xbinedges}, vector<double>{ybinedges}, title, xtitle, ytitle, ztitle);
}
Histo2DPtr & Analysis::book(Histo2DPtr & h2d,const string& hname,
const vector<double>& xbinedges,
const vector<double>& ybinedges,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
const string path = histoPath(hname);
Histo2D hist(xbinedges, ybinedges, path, title);
hist.setAnnotation("XLabel", xtitle);
hist.setAnnotation("YLabel", ytitle);
hist.setAnnotation("ZLabel", ztitle);
h2d = Histo2DPtr(handler().numWeights(), hist);
addAnalysisObject(h2d);
MSG_TRACE("Made 2D histogram " << hname << " for " << name());
return h2d;
}
Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname,
size_t nbins, double lower, double upper,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
Profile1D prof(nbins, lower, upper, path, title);
prof.setAnnotation("XLabel", xtitle);
prof.setAnnotation("YLabel", ytitle);
p1d = Profile1DPtr(handler().numWeights(), prof);
addAnalysisObject(p1d);
MSG_TRACE("Made profile histogram " << hname << " for " << name());
return p1d;
}
Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname,
const initializer_list<double>& binedges,
const string& title,
const string& xtitle,
const string& ytitle) {
return book(p1d, hname, vector<double>{binedges}, title, xtitle, ytitle);
}
Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname,
const vector<double>& binedges,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
Profile1D prof(binedges, path, title);
prof.setAnnotation("XLabel", xtitle);
prof.setAnnotation("YLabel", ytitle);
p1d = Profile1DPtr(handler().numWeights(), prof);
addAnalysisObject(p1d);
MSG_TRACE("Made profile histogram " << hname << " for " << name());
return p1d;
}
Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname,
const Scatter2D& refscatter,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
Profile1D prof(refscatter, path);
prof.setTitle(title);
prof.setAnnotation("XLabel", xtitle);
prof.setAnnotation("YLabel", ytitle);
p1d = Profile1DPtr(handler().numWeights(), prof);
addAnalysisObject(p1d);
MSG_TRACE("Made profile histogram " << hname << " for " << name());
return p1d;
// if (prof.hasAnnotation("IsRef")) prof.rmAnnotation("IsRef");
}
Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname,
const string& title,
const string& xtitle,
const string& ytitle) {
const Scatter2D& refdata = refData(hname);
book(p1d, hname, refdata, title, xtitle, ytitle);
return p1d;
}
Profile1DPtr & Analysis::book(Profile1DPtr & p1d,unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
const string& title,
const string& xtitle,
const string& ytitle) {
const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return book(p1d, axisCode, title, xtitle, ytitle);
}
Profile2DPtr & Analysis::book(Profile2DPtr & p2d, const string& hname,
size_t nxbins, double xlower, double xupper,
size_t nybins, double ylower, double yupper,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
const string path = histoPath(hname);
Profile2D prof(nxbins, xlower, xupper, nybins, ylower, yupper, path, title);
prof.setAnnotation("XLabel", xtitle);
prof.setAnnotation("YLabel", ytitle);
prof.setAnnotation("ZLabel", ztitle);
p2d = Profile2DPtr(handler().numWeights(), prof);
addAnalysisObject(p2d);
MSG_TRACE("Made 2D profile histogram " << hname << " for " << name());
return p2d;
}
Profile2DPtr & Analysis::book(Profile2DPtr & p2d, const string& hname,
const initializer_list<double>& xbinedges,
const initializer_list<double>& ybinedges,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
return book(p2d, hname, vector<double>{xbinedges}, vector<double>{ybinedges}, title, xtitle, ytitle, ztitle);
}
Profile2DPtr & Analysis::book(Profile2DPtr & p2d, const string& hname,
const vector<double>& xbinedges,
const vector<double>& ybinedges,
const string& title,
const string& xtitle,
const string& ytitle,
const string& ztitle)
{
const string path = histoPath(hname);
Profile2D prof(xbinedges, ybinedges, path, title);
prof.setAnnotation("XLabel", xtitle);
prof.setAnnotation("YLabel", ytitle);
prof.setAnnotation("ZLabel", ztitle);
p2d = Profile2DPtr(handler().numWeights(), prof);
addAnalysisObject(p2d);
MSG_TRACE("Made 2D profile histogram " << hname << " for " << name());
return p2d;
}
Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
bool copy_pts,
const string& title,
const string& xtitle,
const string& ytitle) {
const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
return book(s2d, axisCode, copy_pts, title, xtitle, ytitle);
}
Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, const string& hname,
bool copy_pts,
const string& title,
const string& xtitle,
const string& ytitle) {
Scatter2D scat;
const string path = histoPath(hname);
if (copy_pts) {
const Scatter2D& refdata = refData(hname);
scat = Scatter2D(refdata, path);
for (Point2D& p : scat.points()) p.setY(0, 0);
} else {
scat = Scatter2D(path);
}
scat.setTitle(title);
scat.setAnnotation("XLabel", xtitle);
scat.setAnnotation("YLabel", ytitle);
s2d = Scatter2DPtr(handler().numWeights(), scat);
addAnalysisObject(s2d);
MSG_TRACE("Made scatter " << hname << " for " << name());
// if (scat.hasAnnotation("IsRef")) scat.rmAnnotation("IsRef");
return s2d;
}
Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, const string& hname,
size_t npts, double lower, double upper,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
Scatter2D scat;
const double binwidth = (upper-lower)/npts;
for (size_t pt = 0; pt < npts; ++pt) {
const double bincentre = lower + (pt + 0.5) * binwidth;
scat.addPoint(bincentre, 0, binwidth/2.0, 0);
}
scat.setTitle(title);
scat.setAnnotation("XLabel", xtitle);
scat.setAnnotation("YLabel", ytitle);
s2d = Scatter2DPtr(handler().numWeights(), scat);
addAnalysisObject(s2d);
MSG_TRACE("Made scatter " << hname << " for " << name());
return s2d;
}
Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, const string& hname,
const vector<double>& binedges,
const string& title,
const string& xtitle,
const string& ytitle) {
const string path = histoPath(hname);
Scatter2D scat;
for (size_t pt = 0; pt < binedges.size()-1; ++pt) {
const double bincentre = (binedges[pt] + binedges[pt+1]) / 2.0;
const double binwidth = binedges[pt+1] - binedges[pt];
scat.addPoint(bincentre, 0, binwidth/2.0, 0);
}
scat.setTitle(title);
scat.setAnnotation("XLabel", xtitle);
scat.setAnnotation("YLabel", ytitle);
s2d = Scatter2DPtr(handler().numWeights(), scat);
addAnalysisObject(s2d);
MSG_TRACE("Made scatter " << hname << " for " << name());
return s2d;
}
void Analysis::divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const {
const string path = s->path();
*s = *c1 / *c2;
s->setPath(path);
}
void Analysis::divide(const Counter& c1, const Counter& c2, Scatter1DPtr s) const {
const string path = s->path();
*s = c1 / c2;
s->setPath(path);
}
void Analysis::divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const {
const string path = s->path();
*s = *h1 / *h2;
s->setPath(path);
}
void Analysis::divide(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const {
const string path = s->path();
*s = h1 / h2;
s->setPath(path);
}
void Analysis::divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const {
const string path = s->path();
*s = *p1 / *p2;
s->setPath(path);
}
void Analysis::divide(const Profile1D& p1, const Profile1D& p2, Scatter2DPtr s) const {
const string path = s->path();
*s = p1 / p2;
s->setPath(path);
}
void Analysis::divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const {
const string path = s->path();
*s = *h1 / *h2;
s->setPath(path);
}
void Analysis::divide(const Histo2D& h1, const Histo2D& h2, Scatter3DPtr s) const {
const string path = s->path();
*s = h1 / h2;
s->setPath(path);
}
void Analysis::divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const {
const string path = s->path();
*s = *p1 / *p2;
s->setPath(path);
}
void Analysis::divide(const Profile2D& p1, const Profile2D& p2, Scatter3DPtr s) const {
const string path = s->path();
*s = p1 / p2;
s->setPath(path);
}
/// @todo Counter and Histo2D efficiencies and asymms
void Analysis::efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const {
const string path = s->path();
*s = YODA::efficiency(*h1, *h2);
s->setPath(path);
}
void Analysis::efficiency(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const {
const string path = s->path();
*s = YODA::efficiency(h1, h2);
s->setPath(path);
}
void Analysis::asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const {
const string path = s->path();
*s = YODA::asymm(*h1, *h2);
s->setPath(path);
}
void Analysis::asymm(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const {
const string path = s->path();
*s = YODA::asymm(h1, h2);
s->setPath(path);
}
void Analysis::scale(CounterPtr cnt, double factor) {
if (!cnt) {
MSG_WARNING("Failed to scale counter=NULL in analysis " << name() << " (scale=" << factor << ")");
return;
}
if (std::isnan(factor) || std::isinf(factor)) {
MSG_WARNING("Failed to scale counter=" << cnt->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")");
factor = 0;
}
MSG_TRACE("Scaling counter " << cnt->path() << " by factor " << factor);
try {
cnt->scaleW(factor);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not scale counter " << cnt->path());
return;
}
}
void Analysis::normalize(Histo1DPtr histo, double norm, bool includeoverflows) {
if (!histo) {
MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")");
return;
}
MSG_TRACE("Normalizing histo " << histo->path() << " to " << norm);
try {
histo->normalize(norm, includeoverflows);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not normalize histo " << histo->path());
return;
}
}
void Analysis::scale(Histo1DPtr histo, double factor) {
if (!histo) {
MSG_WARNING("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")");
return;
}
if (std::isnan(factor) || std::isinf(factor)) {
MSG_WARNING("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")");
factor = 0;
}
MSG_TRACE("Scaling histo " << histo->path() << " by factor " << factor);
try {
histo->scaleW(factor);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not scale histo " << histo->path());
return;
}
}
void Analysis::normalize(Histo2DPtr histo, double norm, bool includeoverflows) {
if (!histo) {
MSG_ERROR("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")");
return;
}
MSG_TRACE("Normalizing histo " << histo->path() << " to " << norm);
try {
histo->normalize(norm, includeoverflows);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not normalize histo " << histo->path());
return;
}
}
void Analysis::scale(Histo2DPtr histo, double factor) {
if (!histo) {
MSG_ERROR("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")");
return;
}
if (std::isnan(factor) || std::isinf(factor)) {
MSG_ERROR("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")");
factor = 0;
}
MSG_TRACE("Scaling histo " << histo->path() << " by factor " << factor);
try {
histo->scaleW(factor);
} catch (YODA::Exception& we) {
MSG_WARNING("Could not scale histo " << histo->path());
return;
}
}
void Analysis::integrate(Histo1DPtr h, Scatter2DPtr s) const {
// preserve the path info
const string path = s->path();
*s = toIntegralHisto(*h);
s->setPath(path);
}
void Analysis::integrate(const Histo1D& h, Scatter2DPtr s) const {
// preserve the path info
const string path = s->path();
*s = toIntegralHisto(h);
s->setPath(path);
}
-
+}
/// @todo 2D versions of integrate... defined how, exactly?!?
//////////////////////////////////
+namespace {
+ void errormsg(std::string name) {
+#ifdef HAVE_BACKTRACE
+ void * buffer[4];
+ backtrace(buffer, 4);
+ backtrace_symbols_fd(buffer, 4 , 1);
+#endif
+ std::cerr << name << ": Can't book objects outside of init().\n";
+ assert(false);
+ }
+}
+namespace Rivet {
// @todo
// special handling for scatters
void Analysis::addAnalysisObject(const shared_ptr<Scatter1DPtr>& ao) {
- _scatters.push_back(ao);
+ if (handler().stage() == AnalysisHandler::Stage::INIT) {
+ _scatters.push_back(ao);
+ ao.get()->blockDestructor(true);
+ }
+ else {
+ errormsg(name());
+ }
}
void Analysis::addAnalysisObject(const shared_ptr<Scatter2DPtr>& ao) {
- _scatters.push_back(ao);
+ if (handler().stage() == AnalysisHandler::Stage::INIT) {
+ _scatters.push_back(ao);
+ ao.get()->blockDestructor(true);
+ }
+ else {
+ errormsg(name());
+ }
}
void Analysis::addAnalysisObject(const shared_ptr<Scatter3DPtr>& ao) {
- _scatters.push_back(ao);
+ if (handler().stage() == AnalysisHandler::Stage::INIT) {
+ _scatters.push_back(ao);
+ ao.get()->blockDestructor(true);
+ }
+ else {
+ errormsg(name());
+ }
}
void Analysis::addAnalysisObject(MultiweightAOPtr & ao) {
- _analysisobjects.push_back(ao);
+ if (handler().stage() == AnalysisHandler::Stage::INIT) {
+ _analysisobjects.push_back(ao);
+ ao.blockDestructor(true);
+ }
+ else {
+ errormsg(name());
+ }
}
void Analysis::removeAnalysisObject(const string& path) {
for (auto it = _analysisobjects.begin();
it != _analysisobjects.end(); ++it) {
if ((*it).get()->path() == path) {
_analysisobjects.erase(it);
break;
}
}
for (auto it = _scatters.begin(); it != _scatters.end(); ++it) {
if ((**it)->path() == path) {
_scatters.erase(it);
break;
}
}
}
/// @todo can we really remove (multiweighted) analysis objects by == operator??
void Analysis::removeAnalysisObject(const Scatter1DPtr& ao) {
for (auto it = _scatters.begin(); it != _scatters.end(); ++it) {
if (**it == ao) {
_scatters.erase(it);
break;
}
}
}
void Analysis::removeAnalysisObject(const Scatter2DPtr& ao) {
for (auto it = _scatters.begin(); it != _scatters.end(); ++it) {
if (**it == ao) {
_scatters.erase(it);
break;
}
}
}
void Analysis::removeAnalysisObject(const Scatter3DPtr& ao) {
for (auto it = _scatters.begin(); it != _scatters.end(); ++it) {
if (**it == ao) {
_scatters.erase(it);
break;
}
}
}
void Analysis::removeAnalysisObject(const MultiweightAOPtr& ao) {
for (auto it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) {
if ((*it).get() == ao) {
_analysisobjects.erase(it);
break;
}
}
}
}
diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc
--- a/src/Core/AnalysisHandler.cc
+++ b/src/Core/AnalysisHandler.cc
@@ -1,519 +1,524 @@
// -*- C++ -*-
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/AnalysisHandler.hh"
#include "Rivet/Analysis.hh"
#include "Rivet/Tools/ParticleName.hh"
#include "Rivet/Tools/BeamConstraint.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Projections/Beam.hh"
#include "YODA/ReaderYODA.h"
#include "YODA/WriterYODA.h"
#include <regex>
namespace {
inline std::vector<std::string> split(const std::string& input, const std::string& regex) {
// passing -1 as the submatch index parameter performs splitting
std::regex re(regex);
std::sregex_token_iterator
first{input.begin(), input.end(), re, -1},
last;
return {first, last};
}
}
namespace Rivet {
AnalysisHandler::AnalysisHandler(const string& runname)
: _runname(runname),
_eventCounter(0, Counter()), _xs(NAN),
_initialised(false), _ignoreBeams(false)
{}
AnalysisHandler::~AnalysisHandler()
{}
Log& AnalysisHandler::getLog() const {
return Log::getLog("Rivet.Analysis.Handler");
}
/// http://stackoverflow.com/questions/4654636/how-to-determine-if-a-string-is-a-number-with-c
bool is_number(const std::string& s)
{
std::string::const_iterator it = s.begin();
while (it != s.end() && std::isdigit(*it)) ++it;
return !s.empty() && it == s.end();
}
/// Check if any of the weightnames is not a number
bool AnalysisHandler::haveNamedWeights() {
bool dec=false;
for (unsigned int i=0;i<_weightNames.size();++i) {
string s = _weightNames[i];
if (!is_number(s)) {
dec=true;
break;
}
}
return dec;
}
void AnalysisHandler::init(const GenEvent& ge) {
if (_initialised)
throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!");
setRunBeams(Rivet::beams(ge));
MSG_DEBUG("Initialising the analysis handler");
_eventNumber = ge.event_number();
setWeightNames(ge);
if (haveNamedWeights())
MSG_INFO("Using named weights");
else
MSG_INFO("NOT using named weights. Using first weight as nominal weight");
_numWeightTypes = _weightNames.size();
_eventCounter = CounterPtr(_numWeightTypes, Counter("_EVTCOUNT"));
// Check that analyses are beam-compatible, and remove those that aren't
const size_t num_anas_requested = analysisNames().size();
vector<string> anamestodelete;
for (const AnaHandle a : _analyses) {
if (!_ignoreBeams && !a->isCompatible(beams())) {
//MSG_DEBUG(a->name() << " requires beams " << a->requiredBeams() << " @ " << a->requiredEnergies() << " GeV");
anamestodelete.push_back(a->name());
}
}
for (const string& aname : anamestodelete) {
MSG_WARNING("Analysis '" << aname << "' is incompatible with the provided beams: removing");
removeAnalysis(aname);
}
if (num_anas_requested > 0 && analysisNames().empty()) {
cerr << "All analyses were incompatible with the first event's beams\n"
<< "Exiting, since this probably wasn't intentional!" << endl;
exit(1);
}
// Warn if any analysis' status is not unblemished
for (const AnaHandle a : analyses()) {
if (toUpper(a->status()) == "PRELIMINARY") {
MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!");
} else if (toUpper(a->status()) == "OBSOLETE") {
MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!");
} else if (toUpper(a->status()).find("UNVALIDATED") != string::npos) {
MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!");
}
}
// Initialize the remaining analyses
+ _stage = Stage::INIT;
for (AnaHandle a : _analyses) {
MSG_DEBUG("Initialising analysis: " << a->name());
try {
// Allow projection registration in the init phase onwards
a->_allowProjReg = true;
a->init();
//MSG_DEBUG("Checking consistency of analysis: " << a->name());
//a->checkConsistency();
} catch (const Error& err) {
cerr << "Error in " << a->name() << "::init method: " << err.what() << endl;
exit(1);
}
MSG_DEBUG("Done initialising analysis: " << a->name());
}
+ _stage = Stage::OTHER;
_initialised = true;
MSG_DEBUG("Analysis handler initialised");
}
void AnalysisHandler::setWeightNames(const GenEvent& ge) {
/// reroute the print output to a stringstream and process
/// The iteration is done over a map in hepmc2 so this is safe
ostringstream stream;
ge.weights().print(stream); // Super lame, I know
string str = stream.str();
std::regex re("(([^()]+))"); // Regex for stuff enclosed by parentheses ()
for(std::sregex_iterator i = std::sregex_iterator(str.begin(), str.end(), re);
i != std::sregex_iterator(); ++i ) {
std::smatch m = *i;
vector<string> temp = ::split(m.str(), "[,]");
if (temp.size() ==2) {
MSG_DEBUG("Name of weight #" << _weightNames.size() << ": " << temp[0]);
_weightNames.push_back(temp[0]);
}
}
}
void AnalysisHandler::analyze(const GenEvent& ge) {
// Call init with event as template if not already initialised
if (!_initialised) init(ge);
assert(_initialised);
// Ensure that beam details match those from the first event (if we're checking beams)
if ( !_ignoreBeams ) {
const PdgIdPair beams = Rivet::beamIds(ge);
const double sqrts = Rivet::sqrtS(ge);
if (!compatible(beams, _beams) || !fuzzyEquals(sqrts, sqrtS())) {
cerr << "Event beams mismatch: "
<< PID::toBeamsString(beams) << " @ " << sqrts/GeV << " GeV" << " vs. first beams "
<< this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl;
exit(1);
}
}
// Create the Rivet event wrapper
/// @todo Filter/normalize the event here
Event event(ge);
// won't happen for first event because _eventNumber is set in
// init()
if (_eventNumber != ge.event_number()) {
/// @todo
/// can we get away with not passing a matrix?
MSG_TRACE("AnalysisHandler::analyze(): Pushing _eventCounter to persistent.");
_eventCounter.pushToPersistent(_subEventWeights);
// if this is indeed a new event, push the temporary
// histograms and reset
for (const AnaHandle& a : _analyses) {
for (MultiweightAOPtr & ao : a->analysisObjects()) {
MSG_TRACE("AnalysisHandler::analyze(): Pushing " << a->name() << "'s " << ao->name() << " to persistent.");
ao.pushToPersistent(_subEventWeights);
}
MSG_TRACE("AnalysisHandler::analyze(): finished pushing " << a->name() << "'s objects to persistent.");
}
_eventNumber = ge.event_number();
MSG_DEBUG("nominal event # " << _eventCounter._persistent[0]->numEntries());
MSG_DEBUG("nominal sum of weights: " << _eventCounter._persistent[0]->sumW());
MSG_DEBUG("Event has " << _subEventWeights.size() << " sub events.");
_subEventWeights.clear();
}
_eventCounter.newSubEvent();
for (const AnaHandle& a : _analyses) {
for (MultiweightAOPtr & ao : a->analysisObjects()) {
ao.newSubEvent();
}
}
_subEventWeights.push_back(event.weights());
MSG_DEBUG("Analyzing subevent #" << _subEventWeights.size() - 1 << ".");
// Cross-section
#ifdef HEPMC_HAS_CROSS_SECTION
if (ge.cross_section()) {
_xs = ge.cross_section()->cross_section();
_xserr = ge.cross_section()->cross_section_error();
}
#endif
_eventCounter->fill();
// Run the analyses
for (AnaHandle a : _analyses) {
MSG_TRACE("About to run analysis " << a->name());
try {
a->analyze(event);
} catch (const Error& err) {
cerr << "Error in " << a->name() << "::analyze method: " << err.what() << endl;
exit(1);
}
MSG_TRACE("Finished running analysis " << a->name());
}
}
void AnalysisHandler::analyze(const GenEvent* ge) {
if (ge == NULL) {
MSG_ERROR("AnalysisHandler received null pointer to GenEvent");
//throw Error("AnalysisHandler received null pointer to GenEvent");
}
analyze(*ge);
}
void AnalysisHandler::finalize() {
if (!_initialised) return;
MSG_INFO("Finalising analyses");
MSG_TRACE("AnalysisHandler::finalize(): Pushing analysis objects to persistent.");
_eventCounter.pushToPersistent(_subEventWeights);
for (const AnaHandle& a : _analyses) {
for (MultiweightAOPtr & ao : a->analysisObjects())
ao.pushToPersistent(_subEventWeights);
}
for (const AnaHandle& a : _analyses) {
a->setCrossSection(_xs);
for (size_t iW = 0; iW < numWeights(); iW++) {
_eventCounter.setActiveWeightIdx(iW);
for (MultiweightAOPtr & ao : a->analysisObjects())
ao.setActiveWeightIdx(iW);
MSG_TRACE("running " << a->name() << "::finalize() for weight " << iW << ".");
try {
a->finalize();
} catch (const Error& err) {
cerr << "Error in " << a->name() << "::finalize method: " << err.what() << endl;
exit(1);
}
}
+ // allow AO destruction again
+ for (MultiweightAOPtr & ao : a->analysisObjects())
+ ao.blockDestructor(false);
}
// Print out number of events processed
MSG_INFO("Processed " << numEvents() << " event" << (numEvents() == 1 ? "" : "s"));
// // Delete analyses
// MSG_DEBUG("Deleting analyses");
// _analyses.clear();
// Print out MCnet boilerplate
cout << endl;
cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << endl;
cout << "Please acknowledge plots made with Rivet analyses, and cite arXiv:1003.0694 (http://arxiv.org/abs/1003.0694)" << endl;
}
AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname) {
// Check for a duplicate analysis
/// @todo Might we want to be able to run an analysis twice, with different params?
/// Requires avoiding histo tree clashes, i.e. storing the histos on the analysis objects.
for (const AnaHandle& a : _analyses) {
if (a->name() == analysisname) {
MSG_WARNING("Analysis '" << analysisname << "' already registered: skipping duplicate");
return *this;
}
}
AnaHandle analysis( AnalysisLoader::getAnalysis(analysisname) );
if (analysis.get() != 0) { // < Check for null analysis.
MSG_DEBUG("Adding analysis '" << analysisname << "'");
analysis->_analysishandler = this;
_analyses.insert(analysis);
} else {
MSG_WARNING("Analysis '" << analysisname << "' not found.");
}
// MSG_WARNING(_analyses.size());
// for (const AnaHandle& a : _analyses) MSG_WARNING(a->name());
return *this;
}
AnalysisHandler& AnalysisHandler::removeAnalysis(const string& analysisname) {
std::shared_ptr<Analysis> toremove;
for (const AnaHandle a : _analyses) {
if (a->name() == analysisname) {
toremove = a;
break;
}
}
if (toremove.get() != 0) {
MSG_DEBUG("Removing analysis '" << analysisname << "'");
_analyses.erase(toremove);
}
return *this;
}
void AnalysisHandler::addData(const std::vector<YODA::AnalysisObjectPtr>& aos) {
for (const YODA::AnalysisObjectPtr ao : aos) {
const string path = ao->path();
if (path.size() > 1) { // path > "/"
try {
const string ananame = ::split(path, "/")[0];
AnaHandle a = analysis(ananame);
//MultiweightAOPtr mao = ????; /// @todo generate right Multiweight object from ao
//a->addAnalysisObject(mao); /// @todo Need to statistically merge...
} catch (const Error& e) {
MSG_WARNING(e.what());
}
}
}
}
void AnalysisHandler::readData(const string& filename) {
vector<YODA::AnalysisObjectPtr> aos;
try {
/// @todo Use new YODA SFINAE to fill the smart ptr vector directly
vector<YODA::AnalysisObject*> aos_raw;
YODA::ReaderYODA::read(filename, aos_raw);
for (YODA::AnalysisObject* aor : aos_raw) aos.push_back(YODA::AnalysisObjectPtr(aor));
} catch (const YODA::ReadError & e) {
throw UserError("Unexpected error in reading file: " + filename);
}
if (!aos.empty()) addData(aos);
}
vector<reference_wrapper<MultiweightAOPtr> > AnalysisHandler::getRivetAOs() const {
vector<reference_wrapper<MultiweightAOPtr> > rtn;
for (AnaHandle a : _analyses) {
for (const auto & ao : a->analysisObjects()) {
rtn.push_back(ao);
}
}
// Should event counter be included here?
rtn.push_back(_eventCounter);
return rtn;
}
vector<YODA::AnalysisObjectPtr> AnalysisHandler::getYodaAOs() const {
vector<YODA::AnalysisObjectPtr> rtn;
for (MultiweightAOPtr & rao : getRivetAOs()) {
// need to set the index
// before we can search the PATH
rao.setActiveWeightIdx(0);
if (rao->path().find("/TMP/") != string::npos)
continue;
for (size_t iW = 0; iW < numWeights(); iW++) {
rao.setActiveWeightIdx(iW);
// add the weight name in brackets unless we recognize a
// nominal weight
if (_weightNames[iW] != "Weight"
&& _weightNames[iW] != "0"
&& _weightNames[iW] != "Default")
rao->setPath(rao->path() + "[" + _weightNames[iW] + "]");
rtn.push_back(rao.activeYODAPtr());
}
}
YODA::Scatter1D::Points pts; pts.insert(YODA::Point1D(_xs, _xserr));
rtn.push_back( make_shared<Scatter1D>(pts, "/_XSEC") );
sort(rtn.begin(), rtn.end(),
[](YODA::AnalysisObjectPtr a, YODA::AnalysisObjectPtr b) {
return a->path() < b->path();
}
);
return rtn;
}
vector<YODA::AnalysisObjectPtr> AnalysisHandler::getData() const {
return getYodaAOs();
}
void AnalysisHandler::writeData(const string& filename) const {
const vector<YODA::AnalysisObjectPtr> aos = getData();
try {
YODA::WriterYODA::write(filename, aos.begin(), aos.end());
} catch ( YODA::WriteError ) {
throw UserError("Unexpected error in writing file: " + filename);
}
}
string AnalysisHandler::runName() const { return _runname; }
size_t AnalysisHandler::numEvents() const { return _eventCounter->numEntries(); }
/*
* why is this here?
void AnalysisHandler::setSumOfWeights(const double& sum) {
sumOfWeights() = sum;
}
*/
std::vector<std::string> AnalysisHandler::analysisNames() const {
std::vector<std::string> rtn;
for (AnaHandle a : _analyses) {
rtn.push_back(a->name());
}
return rtn;
}
const AnaHandle AnalysisHandler::analysis(const std::string& analysisname) const {
for (const AnaHandle a : analyses())
if (a->name() == analysisname) return a;
throw Error("No analysis named '" + analysisname + "' registered in AnalysisHandler");
}
AnalysisHandler& AnalysisHandler::addAnalyses(const std::vector<std::string>& analysisnames) {
for (const string& aname : analysisnames) {
//MSG_DEBUG("Adding analysis '" << aname << "'");
addAnalysis(aname);
}
return *this;
}
AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector<std::string>& analysisnames) {
for (const string& aname : analysisnames) {
removeAnalysis(aname);
}
return *this;
}
bool AnalysisHandler::needCrossSection() const {
bool rtn = false;
for (const AnaHandle a : _analyses) {
if (!rtn) rtn = a->needsCrossSection();
if (rtn) break;
}
return rtn;
}
AnalysisHandler& AnalysisHandler::setCrossSection(double xs) {
_xs = xs;
return *this;
}
bool AnalysisHandler::hasCrossSection() const {
return (!std::isnan(crossSection()));
}
AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) {
analysis->_analysishandler = this;
_analyses.insert(AnaHandle(analysis));
return *this;
}
PdgIdPair AnalysisHandler::beamIds() const {
return Rivet::beamIds(beams());
}
double AnalysisHandler::sqrtS() const {
return Rivet::sqrtS(beams());
}
void AnalysisHandler::setIgnoreBeams(bool ignore) {
_ignoreBeams=ignore;
}
}
diff --git a/src/Tools/RivetYODA.cc b/src/Tools/RivetYODA.cc
--- a/src/Tools/RivetYODA.cc
+++ b/src/Tools/RivetYODA.cc
@@ -1,368 +1,392 @@
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Tools/RivetYODA.hh"
#include "Rivet/Tools/RivetPaths.hh"
#include "YODA/ReaderYODA.h"
#include "YODA/ReaderAIDA.h"
+#include "DummyConfig.hh"
+#ifdef HAVE_EXECINFO_H
+#include <execinfo.h>
+#endif
+
using namespace std;
namespace Rivet {
+template <class T>
+Wrapper<T>::~Wrapper() {
+ // checking the use_count, because we only want to prevent destruction of the last one,
+ // which seems to have a count of 2 for some reason
+ if ( _blockDestructor // we're a registered AO
+ && !_persistent.empty() // with some entries
+ && _persistent[0] // that are non-null
+ && _persistent[0].use_count() <= 2 ) { // and have no refs left => can't destruct
+#ifdef HAVE_BACKTRACE
+ void * buffer[4];
+ backtrace(buffer, 4);
+ backtrace_symbols_fd(buffer, 4 , 1);
+#endif
+ cerr << "***\n"
+ << "* Cannot destruct temporary AO before finalize.\n"
+ << "* All booked AOs must be class members.\n"
+ << "***\n";
+ assert(false);
+ }
+}
template <class T>
-Wrapper<T>::Wrapper(size_t len_of_weightvec, const T & p) {
+Wrapper<T>::Wrapper(size_t len_of_weightvec, const T & p)
+{
for (size_t m = 0; m < len_of_weightvec; ++m)
_persistent.push_back(make_shared<T>(p));
}
template <class T>
typename T::Ptr Wrapper<T>::active() const {
- assert(_active);
+ if ( !_active ) {
+ assert(false && "No active pointer set. Was this object booked in init()?");
+ }
return _active;
}
template <class T>
void Wrapper<T>::newSubEvent() {
typename TupleWrapper<T>::Ptr tmp
= make_shared<TupleWrapper<T>>(_persistent[0]->clone());
tmp->reset();
_evgroup.push_back( tmp );
_active = _evgroup.back();
assert(_active);
}
string getDatafilePath(const string& papername) {
/// Try to find YODA otherwise fall back to try AIDA
const string path1 = findAnalysisRefFile(papername + ".yoda");
if (!path1.empty()) return path1;
const string path2 = findAnalysisRefFile(papername + ".aida");
if (!path2.empty()) return path2;
throw Rivet::Error("Couldn't find ref data file '" + papername + ".yoda" +
" in data path, '" + getRivetDataPath() + "', or '.'");
}
map<string, YODA::AnalysisObjectPtr> getRefData(const string& papername) {
const string datafile = getDatafilePath(papername);
// Make an appropriate data file reader and read the data objects
/// @todo Remove AIDA support some day...
YODA::Reader& reader = (datafile.find(".yoda") != string::npos) ? \
YODA::ReaderYODA::create() : YODA::ReaderAIDA::create();
vector<YODA::AnalysisObject *> aovec;
reader.read(datafile, aovec);
// Return value, to be populated
map<string, YODA::AnalysisObjectPtr> rtn;
foreach ( YODA::AnalysisObject* ao, aovec ) {
YODA::AnalysisObjectPtr refdata(ao);
if (!refdata) continue;
const string plotpath = refdata->path();
// Split path at "/" and only return the last field, i.e. the histogram ID
const size_t slashpos = plotpath.rfind("/");
const string plotname = (slashpos+1 < plotpath.size()) ? plotpath.substr(slashpos+1) : "";
rtn[plotname] = refdata;
}
return rtn;
}
}
namespace {
using Rivet::Fill;
using Rivet::Fills;
using Rivet::TupleWrapper;
template <class T>
double get_window_size(const typename T::Ptr & histo,
typename T::BinType x) {
// the bin index we fall in
const auto binidx = histo->binIndexAt(x);
// gaps, overflow, underflow don't contribute
if ( binidx == -1 )
return 0;
const auto & b = histo->bin(binidx);
// if we don't have a valid neighbouring bin,
// we use infinite width
typename T::Bin b1(-1.0/0.0, 1.0/0.0);
// points in the top half compare to the upper neighbour
if ( x > b.xMid() ) {
size_t nextidx = binidx + 1;
if ( nextidx < histo->bins().size() )
b1 = histo->bin(nextidx);
}
else { // compare to the lower neighbour
int nextidx = binidx - 1;
if ( nextidx >= 0 )
b1 = histo->bin(nextidx);
}
// the factor 2 is arbitrary, could poss. be smaller
return min( b.width(), b1.width() ) / 2.0;
}
template <class T>
typename T::BinType
fillT2binT(typename T::FillType a) {
return a;
}
template <>
YODA::Profile1D::BinType
fillT2binT<YODA::Profile1D>(YODA::Profile1D::FillType a) {
return get<0>(a);
}
template <>
YODA::Profile2D::BinType
fillT2binT<YODA::Profile2D>(YODA::Profile2D::FillType a) {
return YODA::Profile2D::BinType{ get<0>(a), get<1>(a) };
}
template <class T>
void commit(vector<typename T::Ptr> & persistent,
const vector< vector<Fill<T>> > & tuple,
const vector<valarray<double>> & weights ) {
// TODO check if all the xs are in the same bin anyway!
// Then no windowing needed
assert(persistent.size() == weights[0].size());
for ( const auto & x : tuple ) {
double maxwindow = 0.0;
for ( const auto & xi : x ) {
// TODO check for NOFILL here
// persistent[0] has the same binning as all the other persistent objects
double window = get_window_size<T>(persistent[0], fillT2binT<T>(xi.first));
if ( window > maxwindow )
maxwindow = window;
}
const double wsize = maxwindow;
// all windows have same size
set<double> edgeset;
// bin edges need to be in here!
for ( const auto & xi : x ) {
edgeset.insert(fillT2binT<T>(xi.first) - wsize);
edgeset.insert(fillT2binT<T>(xi.first) + wsize);
}
vector< std::tuple<double,valarray<double>,double> > hfill;
double sumf = 0.0;
auto edgit = edgeset.begin();
double ehi = *edgit;
while ( ++edgit != edgeset.end() ) {
double elo = ehi;
ehi = *edgit;
valarray<double> sumw(0.0, persistent.size()); // need m copies of this
bool gap = true; // Check for gaps between the sub-windows.
for ( size_t i = 0; i < x.size(); ++i ) {
// check equals comparisons here!
if ( fillT2binT<T>(x[i].first) + wsize >= ehi
&&
fillT2binT<T>(x[i].first) - wsize <= elo ) {
sumw += x[i].second * weights[i];
gap = false;
}
}
if ( gap ) continue;
hfill.push_back( make_tuple( (ehi + elo)/2.0, sumw, (ehi - elo) ) );
sumf += ehi - elo;
}
for ( auto f : hfill )
for ( size_t m = 0; m < persistent.size(); ++m )
persistent[m]->fill( get<0>(f), get<1>(f)[m], get<2>(f)/sumf );
// Note the scaling to one single fill
}
}
template<>
void commit<YODA::Histo2D>(vector<YODA::Histo2D::Ptr> & persistent,
const vector< vector<Fill<YODA::Histo2D>> > & tuple,
const vector<valarray<double>> & weights)
{}
template<>
void commit<YODA::Profile2D>(vector<YODA::Profile2D::Ptr> & persistent,
const vector< vector<Fill<YODA::Profile2D>> > & tuple,
const vector<valarray<double>> & weights)
{}
template <class T>
double distance(T a, T b) {
return abs(a - b);
}
template <>
double distance<tuple<double,double> >(tuple<double,double> a, tuple<double,double> b) {
return Rivet::sqr(get<0>(a) - get<0>(b)) + Rivet::sqr(get<1>(a) - get<1>(b));
}
}
/// fills is a vector of sub-event with an ordered set of x-values of
/// the fills in each sub-event. NOFILL should be an "impossible"
/// value for this histogram. Returns a vector of sub-events with
/// an ordered vector of fills (including NOFILLs) for each sub-event.
template <class T>
vector< vector<Fill<T> > >
match_fills(const vector<typename TupleWrapper<T>::Ptr> & evgroup, const Fill<T> & NOFILL)
{
vector< vector<Fill<T> > > matched;
// First just copy subevents into vectors and find the longest vector.
unsigned int maxfill = 0; // length of biggest vector
int imax = 0; // index position of biggest vector
for ( const auto & it : evgroup ) {
const auto & subev = it->fills();
if ( subev.size() > maxfill ) {
maxfill = subev.size();
imax = matched.size();
}
matched.push_back(vector<Fill<T> >(subev.begin(), subev.end()));
}
// Now, go through all subevents with missing fills.
const vector<Fill<T>> & full = matched[imax]; // the longest one
for ( auto & subev : matched ) {
if ( subev.size() == maxfill ) continue;
// Add NOFILLs to the end;
while ( subev.size() < maxfill ) subev.push_back(NOFILL);
// Iterate from the back and shift all fill values backwards by
// swapping with NOFILLs so that they better match the full
// subevent.
for ( int i = maxfill - 1; i >= 0; --i ) {
if ( subev[i] == NOFILL ) continue;
size_t j = i;
while ( j + 1 < maxfill && subev[j + 1] == NOFILL &&
distance(fillT2binT<T>(subev[j].first),
fillT2binT<T>(full[j].first))
>
distance(fillT2binT<T>(subev[j].first),
fillT2binT<T>(full[j + 1].first)) )
{
swap(subev[j], subev[j + 1]);
++j;
}
}
}
// transpose
vector<vector<Fill<T>>> result(maxfill,vector<Fill<T>>(matched.size()));
for (size_t i = 0; i < matched.size(); ++i)
for (size_t j = 0; j < maxfill; ++j)
result.at(j).at(i) = matched.at(i).at(j);
return result;
}
namespace Rivet {
template <class T>
void Wrapper<T>::pushToPersistent(const vector<valarray<double> >& weight) {
-
assert( _evgroup.size() == weight.size() );
// have we had subevents at all?
const bool have_subevents = _evgroup.size() > 1;
if ( ! have_subevents ) {
// simple replay of all tuple entries
// each recorded fill is inserted into all persistent weightname histos
for ( size_t m = 0; m < _persistent.size(); ++m )
for ( const auto & f : _evgroup[0]->fills() )
_persistent[m]->fill( f.first, f.second * weight[0][m] );
} else {
// outer index is subevent, inner index is jets in the event
vector<vector<Fill<T>>> linedUpXs
= match_fills<T>(_evgroup, {typename T::FillType(), 0.0});
commit<T>( _persistent, linedUpXs, weight );
}
_evgroup.clear();
_active.reset();
}
template <>
void Wrapper<YODA::Counter>::pushToPersistent(const vector<valarray<double> >& weight) {
for ( size_t m = 0; m < _persistent.size(); ++m ) {
for ( size_t n = 0; n < _evgroup.size(); ++n ) {
for ( const auto & f : _evgroup[n]->fills() ) {
_persistent[m]->fill( f.second * weight[n][m] );
}
}
}
_evgroup.clear();
_active.reset();
}
template <>
void Wrapper<YODA::Scatter1D>::pushToPersistent(const vector<valarray<double> >& weight) {
cout << ("WARNING: filling scatters in the event loop is not a well-defined behavior!!") << endl;
-
_evgroup.clear();
_active.reset();
}
template <>
void Wrapper<YODA::Scatter2D>::pushToPersistent(const vector<valarray<double> >& weight) {
cout << ("WARNING: filling scatters in the event loop is not a well-defined behavior!!") << endl;
-
_evgroup.clear();
_active.reset();
}
template <>
void Wrapper<YODA::Scatter3D>::pushToPersistent(const vector<valarray<double> >& weight) {
cout << ("WARNING: filling scatters in the event loop is not a well-defined behavior!!") << endl;
-
_evgroup.clear();
_active.reset();
}
// explicitly instantiate all wrappers
template class Wrapper<YODA::Histo1D>;
template class Wrapper<YODA::Histo2D>;
template class Wrapper<YODA::Profile1D>;
template class Wrapper<YODA::Profile2D>;
template class Wrapper<YODA::Counter>;
template class Wrapper<YODA::Scatter1D>;
template class Wrapper<YODA::Scatter2D>;
template class Wrapper<YODA::Scatter3D>;
}
diff --git a/test/testApi.cc b/test/testApi.cc
--- a/test/testApi.cc
+++ b/test/testApi.cc
@@ -1,34 +1,34 @@
#include "Rivet/AnalysisHandler.hh"
#include "HepMC/GenEvent.h"
#include "HepMC/IO_GenEvent.h"
using namespace std;
int main() {
Rivet::AnalysisHandler ah;
- Rivet::Log::setLevel("Rivet", Rivet::Log::TRACE);
+ Rivet::Log::setLevel("Rivet", Rivet::Log::DEBUG);
// Specify the analyses to be used
ah.addAnalysis("EXAMPLE");
ah.addAnalyses({{ "MC_JETS", "EXAMPLE_CUTS", "EXAMPLE_SMEAR" }});
std::ifstream file("testApi.hepmc");
HepMC::IO_GenEvent hepmcio(file);
HepMC::GenEvent* evt = hepmcio.read_next_event();
while (evt) {
// Analyse current event
ah.analyze(*evt);
// Clean up and get next event
delete evt; evt = nullptr;
hepmcio >> evt;
}
file.close();
ah.setCrossSection(1.0);
ah.finalize();
ah.writeData("out.yoda");
return 0;
}