diff --git a/analyses/pluginMC/MC_PDFS.cc b/analyses/pluginMC/MC_PDFS.cc --- a/analyses/pluginMC/MC_PDFS.cc +++ b/analyses/pluginMC/MC_PDFS.cc @@ -1,109 +1,109 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" // #include "Rivet/Projections/ChargedFinalState.hh" namespace Rivet { /// Generic analysis looking at various distributions of final state particles class MC_PDFS : public Analysis { public: /// Constructor MC_PDFS() : Analysis("MC_PDFS") { } public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Projections // declare(ChargedFinalState((Cuts::etaIn(-5.0, 5.0) && Cuts::pT >= 500*MeV)), "CFS"); // Histograms book(_histPdfX ,"PdfX", logspace(50, 0.000001, 1.0)); book(_histPdfXmin ,"PdfXmin", logspace(50, 0.000001, 1.0)); book(_histPdfXmax ,"PdfXmax", logspace(50, 0.000001, 1.0)); book(_histPdfQ ,"PdfQ", 50, 0.0, 30.0); book(_histPdfXQ,"PdfXQ", logspace(50, 0.000001, 1.0), linspace(50, 0.0, 30.0)); //book( _histPdfTrackptVsX ,"PdfTrackptVsX", logspace(50, 0.000001, 1.0)); //book( _histPdfTrackptVsQ ,"PdfTrackptVsQ", 50, 0.0, 30.0); } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = 1.0; // This analysis needs a valid HepMC PDF info object to do anything if (event.genEvent()->pdf_info() == 0) vetoEvent; PdfInfo pdfi = *(event.genEvent()->pdf_info()); -#ifdef ENABLE_HEPMC_3 +#ifdef RIVET_ENABLE_HEPMC_3 MSG_DEBUG("PDF Q = " << pdfi.scale<< " for (id, x) = " << "(" << pdfi.pdf_id[0] << ", " << pdfi.x[0] << ") " << "(" << pdfi.pdf_id[1] << ", " << pdfi.x[1] << ")"); _histPdfX->fill(pdfi.x[0], weight); _histPdfX->fill(pdfi.x[1], weight); _histPdfXmin->fill(std::min(pdfi.x[0], pdfi.x[1]), weight); _histPdfXmax->fill(std::max(pdfi.x[0], pdfi.x[1]), weight); _histPdfQ->fill(pdfi.scale, weight); // always in GeV? _histPdfXQ->fill(pdfi.x[0], pdfi.scale, weight); // always in GeV? _histPdfXQ->fill(pdfi.x[1], pdfi.scale, weight); // always in GeV? #else MSG_DEBUG("PDF Q = " << pdfi.scalePDF() << " for (id, x) = " << "(" << pdfi.id1() << ", " << pdfi.x1() << ") " << "(" << pdfi.id2() << ", " << pdfi.x2() << ")"); _histPdfX->fill(pdfi.x1(), weight); _histPdfX->fill(pdfi.x2(), weight); _histPdfXmin->fill(std::min(pdfi.x1(), pdfi.x2()), weight); _histPdfXmax->fill(std::max(pdfi.x1(), pdfi.x2()), weight); _histPdfQ->fill(pdfi.scalePDF(), weight); // always in GeV? _histPdfXQ->fill(pdfi.x1(), pdfi.scalePDF(), weight); // always in GeV? _histPdfXQ->fill(pdfi.x2(), pdfi.scalePDF(), weight); // always in GeV? #endif // const FinalState& cfs = apply(event, "CFS"); // for (const Particle& p : cfs.particles()) { // if (fabs(eta) < 2.5 && p.pT() > 10*GeV) { // _histPdfTrackptVsX->fill(pdfi.x1(), p.pT()/GeV, weight); // _histPdfTrackptVsX->fill(pdfi.x2(), p.pT()/GeV, weight); // _histPdfTrackptVsQ->fill(pdfi.scalePDF(), p.pT()/GeV, weight); // } // } } /// Finalize void finalize() { scale(_histPdfX, 1/sumOfWeights()); scale(_histPdfXmin, 1/sumOfWeights()); scale(_histPdfXmax, 1/sumOfWeights()); scale(_histPdfQ, 1/sumOfWeights()); } //@} private: /// @name Histograms //@{ Histo1DPtr _histPdfX, _histPdfXmin, _histPdfXmax, _histPdfQ; Histo2DPtr _histPdfXQ; // Profile1DPtr _histPdfTrackptVsX, _histPdfTrackptVsQ; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(MC_PDFS); } diff --git a/analyses/pluginMC/MC_PRINTEVENT.cc b/analyses/pluginMC/MC_PRINTEVENT.cc --- a/analyses/pluginMC/MC_PRINTEVENT.cc +++ b/analyses/pluginMC/MC_PRINTEVENT.cc @@ -1,267 +1,267 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" namespace Rivet { /// @author Andy Buckley class MC_PRINTEVENT : public Analysis { public: /// Constructor MC_PRINTEVENT() : Analysis("MC_PRINTEVENT") { } /// @name Analysis methods //@{ void init() { /// Set up particle name map // quarks & gluons _pnames[1] = "d"; _pnames[-1] = "d~"; _pnames[2] = "u"; _pnames[-2] = "u~"; _pnames[3] = "s"; _pnames[-3] = "s~"; _pnames[4] = "c"; _pnames[-4] = "c~"; _pnames[5] = "b"; _pnames[-5] = "b~"; _pnames[6] = "t"; _pnames[-6] = "t~"; // bosons _pnames[21] = "g"; _pnames[22] = "gamma"; _pnames[23] = "Z0"; _pnames[24] = "W+"; _pnames[-24] = "W-"; _pnames[25] = "h0"; _pnames[26] = "h0"; // leptons _pnames[11] = "e-"; _pnames[-11] = "e+"; _pnames[13] = "mu-"; _pnames[-13] = "mu+"; _pnames[15] = "tau-"; _pnames[-15] = "tau+"; _pnames[12] = "nu_e"; _pnames[-12] = "nu_e~"; _pnames[14] = "nu_mu"; _pnames[-14] = "nu_mu~"; _pnames[16] = "nu_tau"; _pnames[-16] = "nu_tau~"; // common hadrons _pnames[111] = "pi0"; _pnames[211] = "pi+"; _pnames[-211] = "pi-"; _pnames[221] = "eta"; _pnames[331] = "eta'"; _pnames[113] = "rho0"; _pnames[213] = "rho+"; _pnames[-213] = "rho-"; _pnames[223] = "omega"; _pnames[333] = "phi"; _pnames[130] = "K0L"; _pnames[310] = "K0S"; _pnames[311] = "K0"; _pnames[-311] = "K0"; _pnames[321] = "K+"; _pnames[-321] = "K-"; _pnames[313] = "K*0"; _pnames[-313] = "K*0~"; _pnames[323] = "K*+"; _pnames[-323] = "K*-"; _pnames[411] = "D+"; _pnames[-411] = "D-"; _pnames[421] = "D0"; _pnames[-421] = "D0~"; _pnames[413] = "D*+"; _pnames[-413] = "D*-"; _pnames[423] = "D*0"; _pnames[-423] = "D*0~"; _pnames[431] = "Ds+"; _pnames[-431] = "Ds-"; _pnames[433] = "Ds*+"; _pnames[-433] = "Ds*-"; _pnames[511] = "B0"; _pnames[-511] = "B0~"; _pnames[521] = "B+"; _pnames[-521] = "B-"; _pnames[513] = "B*0"; _pnames[-513] = "B*0~"; _pnames[523] = "B*+"; _pnames[-523] = "B*-"; _pnames[531] = "B0s"; _pnames[541] = "Bc+"; _pnames[-541] = "Bc-"; _pnames[441] = "eta_c(1S)"; _pnames[443] = "J/psi(1S)"; _pnames[551] = "eta_b(1S)"; _pnames[553] = "Upsilon(1S)"; _pnames[2212] = "p+"; _pnames[-2212] = "p-"; _pnames[2112] = "n"; _pnames[-2112] = "n~"; _pnames[2224] = "Delta++"; _pnames[2214] = "Delta+"; _pnames[2114] = "Delta0"; _pnames[1114] = "Delta-"; _pnames[3122] = "Lambda"; _pnames[-3122] = "Lambda~"; _pnames[3222] = "Sigma+"; _pnames[-3222] = "Sigma+~"; _pnames[3212] = "Sigma0"; _pnames[-3212] = "Sigma0~"; _pnames[3112] = "Sigma-"; _pnames[-3112] = "Sigma-~"; _pnames[4122] = "Lambda_c+"; _pnames[-4122] = "Lambda_c-"; _pnames[5122] = "Lambda_b"; // exotic _pnames[32] = "Z'"; _pnames[34] = "W'+"; _pnames[-34] = "W'-"; _pnames[35] = "H0"; _pnames[36] = "A0"; _pnames[37] = "H+"; _pnames[-37] = "H-"; // shower-specific _pnames[91] = "cluster"; _pnames[92] = "string"; _pnames[9922212] = "remn"; _pnames[1103] = "dd"; _pnames[2101] = "ud0"; _pnames[2103] = "ud1"; _pnames[2203] = "uu"; } /// Perform the per-event analysis void analyze(const Event& event) { /// @todo Wouldn't this be nice... if HepMC::IO_AsciiParticles was sane :-/ // printEvent(event.genEvent()); - #ifdef ENABLE_HEPMC_3 + #ifdef RIVET_ENABLE_HEPMC_3 /// @todo gonna try this instead of replicating everything below RivetHepMC::Print::content(*(event.genEvent())); #else const GenEvent* evt = event.genEvent(); cout << string(120, '=') << "\n" << endl; // Weights cout << "Weights(" << evt->weights().size() << ")="; /// @todo Re-enable // for (double w, evt->weights()) // cout << w << " "; cout << "\n" << "EventScale " << evt->event_scale() << " [energy] \t alphaQCD=" << evt->alphaQCD() << "\t alphaQED=" << evt->alphaQED() << endl; if (evt->pdf_info()) { cout << "PdfInfo: id1=" << evt->pdf_info()->id1() << " id2=" << evt->pdf_info()->id2() << " x1=" << evt->pdf_info()->x1() << " x2=" << evt->pdf_info()->x2() << " q=" << evt->pdf_info()->scalePDF() << " xpdf1=" << evt->pdf_info()->pdf1() << " xpdf2=" << evt->pdf_info()->pdf2() << endl; } else { cout << "PdfInfo: EMPTY"; } // Print a legend to describe the particle info char particle_legend[120]; sprintf( particle_legend," %9s %8s %-15s %4s %8s %8s (%9s,%9s,%9s,%9s,%9s)", "Barcode","PDG ID","Name","Stat","ProdVtx","DecayVtx","Px","Py","Pz","E ","m"); cout << endl; cout << " GenParticle Legend\n" << particle_legend << "\n"; // if (m_vertexinfo) { // sprintf( particle_legend," %60s (%9s,%9s,%9s,%9s)"," ","Vx","Vy","Vz","Vct "); // cout << particle_legend << endl; // } // cout << string(120, '_') << endl; // Print all particles // const HepPDT::ParticleDataTable* pdt = m_ppsvc->PDT(); for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p) { int p_bcode = (*p)->barcode(); int p_pdg_id = (*p)->pdg_id(); double p_px = (*p)->momentum().px(); double p_py = (*p)->momentum().py(); double p_pz = (*p)->momentum().pz(); double p_pe = (*p)->momentum().e(); int p_stat = (*p)->status(); int p_prodvtx = 0; if ((*p)->production_vertex() && (*p)->production_vertex()->barcode() != 0) { p_prodvtx = (*p)->production_vertex()->barcode(); } int p_endvtx = 0; if ((*p)->end_vertex() && (*p)->end_vertex()->barcode() != 0) { p_endvtx=(*p)->end_vertex()->barcode(); } // double v_x = 0; // double v_y = 0; // double v_z = 0; // double v_ct = 0; // if ((*p)->production_vertex()) { // v_x = (*p)->production_vertex()->position().x(); // v_y = (*p)->production_vertex()->position().y(); // v_z = (*p)->production_vertex()->position().z(); // v_ct = (*p)->production_vertex()->position().t(); // } // Mass (prefer generated mass if available) double p_mass = (*p)->generated_mass(); if (p_mass == 0 && !(p_stat == 1 && p_pdg_id == 22)) p_mass = (*p)->momentum().m(); // Particle names string sname = (_pnames.find(p_pdg_id) != _pnames.end()) ? _pnames[p_pdg_id] : ""; const char* p_name = sname.c_str() ; char particle_entries[120]; sprintf(particle_entries, " %9i %8i %-15s %4i %8i %8i (%+9.3g,%+9.3g,%+9.3g,%+9.3g,%9.3g)", p_bcode, p_pdg_id, p_name, p_stat, p_prodvtx, p_endvtx, p_px, p_py, p_pz, p_pe, p_mass); cout << particle_entries << "\n"; // if (m_vertexinfo) { // sprintf(particle_entries," %60s (%+9.3g,%+9.3g,%+9.3g,%+9.3g)"," ",v_x, v_y, v_z, v_ct); // cout << particle_entries << "\n"; // } } cout << "\n" << endl; #endif // VERSION_CODE >= 3000000 } /// Normalise histograms etc., after the run void finalize() { } //@} private: map _pnames; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(MC_PRINTEVENT); } diff --git a/analyses/pluginMC/MC_XS.cc b/analyses/pluginMC/MC_XS.cc --- a/analyses/pluginMC/MC_XS.cc +++ b/analyses/pluginMC/MC_XS.cc @@ -1,92 +1,92 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" -#ifndef ENABLE_HEPMC_3 +#ifndef RIVET_ENABLE_HEPMC_3 #include "HepMC/HepMCDefs.h" #endif namespace Rivet { /// @brief Analysis for the generated cross section class MC_XS : public Analysis { public: /// @name Constructors etc. //@{ /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(MC_XS); //@} public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { /// @todo Convert to Scatter1D or Counter book(_h_XS, "XS"); book(_h_N, "N", 1, 0.0, 1.0); book(_h_pmXS, "pmXS", 2, -1.0, 1.0); book(_h_pmN, "pmN", 2, -1.0, 1.0); _mc_xs = _mc_error = 0.; } /// Perform the per-event analysis void analyze(const Event& event) { - #if defined ENABLE_HEPMC_3 + #if defined RIVET_ENABLE_HEPMC_3 //@todo HepMC3::GenCrossSection methods aren't const accessible :( RivetHepMC::GenCrossSection gcs = *(event.genEvent()->cross_section()); _mc_xs = gcs.xsec(); _mc_error = gcs.xsec_err(); #elif defined HEPMC_HAS_CROSS_SECTION _mc_xs = event.genEvent()->cross_section()->cross_section(); _mc_error = event.genEvent()->cross_section()->cross_section_error(); #endif // VERSION_CODE >= 3000000 for (size_t m = 0; m < event.weights().size(); ++m) { const double weight = event.weights()[m]; _h_pmXS.get()->_getPersistent(m)->fill(0.5*(weight > 0 ? 1. : -1), abs(weight)); _h_pmN.get()->_getPersistent(m)->fill(0.5*(weight > 0 ? 1. : -1), 1.); _h_N.get()->_getPersistent(m)->fill(0.5, 1.0); } } /// Normalise histograms etc., after the run void finalize() { scale(_h_pmXS, crossSection()/sumOfWeights()); scale(_h_pmN, 1.0/numEvents()); scale(_h_N, 1.0/numEvents()); #ifndef HEPMC_HAS_CROSS_SECTION _mc_xs = crossSection(); _mc_error = 0.0; #endif _h_XS->addPoint(0, _mc_xs, 0.5, _mc_error); } //@} private: /// @name Histograms //@{ Scatter2DPtr _h_XS; Histo1DPtr _h_pmXS, _h_pmN, _h_N; double _mc_xs, _mc_error; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(MC_XS); } diff --git a/configure.ac b/configure.ac --- a/configure.ac +++ b/configure.ac @@ -1,441 +1,441 @@ ## Process this file with autoconf to produce a configure script. AC_PREREQ(2.59) AC_INIT([Rivet],[3.0.1],[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]) AM_INIT_AUTOMAKE([dist-bzip2 -Wall 1.10]) 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 AC_FUNC_STRERROR_R ## 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 10704; then AC_MSG_ERROR([YODA version isn't sufficient: at least version 1.7.4 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 if test "${with_hepmc+set}" = set; then if test "${with_hepmc3+set}" = set; then AC_MSG_ERROR([HepMC2 *OR* HepMC3 is required. Both were specified!]) fi fi if test "${with_hepmc3+set}" = set; then AC_CEDAR_LIBRARYANDHEADERS([HepMC3], , [use_hepmc3=yes], [use_hepmc3=no]) AM_CONDITIONAL(WITH_HEPMC,false) AM_CONDITIONAL(WITH_HEPMCINC,false) AM_CONDITIONAL(WITH_HEPMCLIB,false) AM_CONDITIONAL(WITHOUT_HEPMC,true) AM_CONDITIONAL(WITHOUT_HEPMCINC,true) AM_CONDITIONAL(WITHOUT_HEPMCLIB,true) else AC_CEDAR_LIBRARYANDHEADERS([HepMC], , [use_hepmc2=yes], [use_hepmc2=no]) AM_CONDITIONAL(WITH_HEPMC3,false) AM_CONDITIONAL(WITH_HEPMC3INC,false) AM_CONDITIONAL(WITH_HEPMC3LIB,false) AM_CONDITIONAL(WITHOUT_HEPMC3,true) AM_CONDITIONAL(WITHOUT_HEPMC3INC,true) AM_CONDITIONAL(WITHOUT_HEPMC3LIB,true) fi if test x$use_hepmc2 = xno && test x$use_hepmc3 = xno ; then AC_MSG_ERROR([HepMC2 or HepMC3 is required]) fi if test x$use_hepmc2 = xyes; then oldCPPFLAGS=$CPPFLAGS CPPFLAGS="$CPPFLAGS -I$HEPMCINCPATH" if test -e "$HEPMCINCPATH/HepMC/HepMCDefs.h"; then AC_LANG_CONFTEST([AC_LANG_SOURCE([#include #include "HepMC/HepMCDefs.h" int main() { std::cout << HEPMC_VERSION << std::endl; return 0; }])]) else AC_LANG_CONFTEST([AC_LANG_SOURCE([#include #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 else oldCPPFLAGS=$CPPFLAGS CPPFLAGS="$CPPFLAGS -I$HEPMC3INCPATH" AC_LANG_CONFTEST([AC_LANG_SOURCE([#include #include "HepMC3/Version.h" int main() { std::cout << HepMC3::version() << std::endl; return 0; }])]) 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` let hepmc_third=`echo "$hepmc_version" | cut -d. -f3` fi rm -f conftest conftest.cpp conftest.cc conftest.C HEPMC_VERSION=$hepmc_major$hepmc_minor$hepmc_third AC_MSG_NOTICE([HepMC version is $hepmc_version -> $HEPMC_VERSION]) if test $HEPMC_VERSION -le 310; then AC_MSG_ERROR([HepMC3 version 3.1.1 or later is required.]) fi AC_SUBST(HEPMC_VERSION) + AC_DEFINE([RIVET_ENABLE_HEPMC_3],[1], [Defined if version 3 of HepMC is used.]) CPPFLAGS=$oldCPPFLAGS fi AM_CONDITIONAL([ENABLE_HEPMC_3], [test x$hepmc_major = x3]) -AM_COND_IF([ENABLE_HEPMC_3],[CPPFLAGS="$CPPFLAGS -DENABLE_HEPMC_3=true"]) ## 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 FASTJETLIBADD="$($FJCONFIG --plugins --shared --libs)" else FASTJETLIBADD="-L$FASTJETLIBPATH -l$FASTJETLIBNAME" FASTJETLIBADD="$FASTJETLIBADD -lSISConePlugin -lsiscone -lsiscone_spherical" FASTJETLIBADD="$FASTJETLIBADD -lCDFConesPlugin -lD0RunIIConePlugin -lNestedDefsPlugin" FASTJETLIBADD="$FASTJETLIBADD -lTrackJetPlugin -lATLASConePlugin -lCMSIterativeConePlugin" FASTJETLIBADD="$FASTJETLIBADD -lEECambridgePlugin -lJadePlugin" fi; FASTJETLDFLAGS="-L$FASTJETLIBPATH" ## FastJet contrib (allow paths specified separately from FJ) AC_ARG_WITH([fjcontrib], [AC_HELP_STRING(--with-fjcontrib, path to FastJet Contrib @<:@$FASTJETPATH@:>@)], [], []) if test -n "$with_fjcontrib"; then AC_MSG_NOTICE([FastJet Contrib base path given: $with_fjcontrib]) FJCONTRIBINCPATH=$with_fjcontrib/include FJCONTRIBLIBPATH=$with_fjcontrib/lib FJCONTRIBCPPFLAGS=-I$FJCONTRIBINCPATH AC_SUBST(FJCONTRIBINCPATH) AC_SUBST(FJCONTRIBLIBPATH) AC_SUBST(FJCONTRIBCPPFLAGS) FASTJETLIBADD="$FASTJETLIBADD -L$FJCONTRIBLIBPATH" FASTJETLDFLAGS="$FASTJETLDFLAGS -L$FJCONTRIBLIBPATH" fi FASTJETLIBADD="$FASTJETLIBADD -lfastjetcontribfragile -lfastjettools" AC_SUBST(FASTJETLIBADD) AC_SUBST(FASTJETLDFLAGS) AC_MSG_NOTICE([FastJet lib flags: $FASTJETLIBADD]) AC_MSG_NOTICE([FastJet ld flags: $FASTJETLDFLAGS]) ## Test FastJet + Contrib ## Check for FastJet headers that require the --enable-all(cxx)plugins option oldCPPFLAGS=$CPPFLAGS CPPFLAGS="$CPPFLAGS $FASTJETCPPFLAGS $FJCONTRIBCPPFLAGS" oldLIBS=$LIBS LIBS="$LIBS $FASTJETLIBADD" FASTJET_ERRMSG="Required FastJet plugin headers were not found: did you build FastJet with the --enable-allcxxplugins option?" AC_CHECK_HEADER([fastjet/D0RunIIConePlugin.hh], [], [AC_MSG_ERROR([$FASTJET_ERRMSG])]) AC_CHECK_HEADER([fastjet/TrackJetPlugin.hh], [], [AC_MSG_ERROR([$FASTJET_ERRMSG])]) FASTJET_ERRMSG="Required FastJet Contrib headers were not found: did you install FastJet Contrib?" AC_CHECK_HEADER([fastjet/contrib/SoftDrop.hh], [], [AC_MSG_ERROR([$FASTJET_ERRMSG])]) AC_MSG_CHECKING([FastJet test program compilation]) FASTJET_ERRMSG="Couldn't build FastJet test program: did you install FastJet Contrib with 'make fragile-shared-install'?" AC_LINK_IFELSE([AC_LANG_SOURCE([#include #include "fastjet/PseudoJet.hh" #include "fastjet/ClusterSequence.hh" #include "fastjet/contrib/SoftDrop.hh" int main() { std::vector particles{{1,1,0,2}, {1,-1,0,2}}; fastjet::ClusterSequence cs(particles, fastjet::JetDefinition(fastjet::antikt_algorithm, 0.4)); fastjet::contrib::SoftDrop sd(0.5, 0.1); for (const fastjet::PseudoJet& j : cs.inclusive_jets()) { const fastjet::PseudoJet& jsd = sd(j); } return 0; } ])], [AC_MSG_RESULT([successful])], [AC_MSG_RESULT([failed]) AC_MSG_ERROR([$FASTJET_ERRMSG])]) CPPFLAGS=$oldCPPFLAGS LIBS=$oldLIBS # ## 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.7.3']) AC_SUBST(PYTHON_VERSION) RIVET_PYTHONPATH=`$PYTHON -c "from __future__ import print_function; 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]) dnl dnl setup.py puts its build artifacts into a labelled path dnl this helps the test scripts to find them locally instead of dnl having to install first dnl RIVET_SETUP_PY_PATH=$(${PYTHON} -c 'from __future__ import print_function; import distutils.util as u, sys; vi=sys.version_info; print("lib.%s-%s.%s" % (u.get_platform(),vi.major, vi.minor))') AC_SUBST(RIVET_SETUP_PY_PATH) ## 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)" AM_CPPFLAGS="$AM_CPPFLAGS -I\$(YODAINCPATH)" if test x$use_hepmc2 = xyes ; then AM_CPPFLAGS="$AM_CPPFLAGS -I\$(HEPMCINCPATH)" else AM_CPPFLAGS="$AM_CPPFLAGS -I\$(HEPMC3INCPATH)" fi AM_CPPFLAGS="$AM_CPPFLAGS \$(FASTJETCPPFLAGS) \$(FJCONTRIBCPPFLAGS)" 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"]) ## Optional zlib support for gzip-compressed data streams/files AX_CHECK_ZLIB ## 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) 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 < -#ifdef ENABLE_HEPMC_3 +#ifdef RIVET_ENABLE_HEPMC_3 #include "HepMC3/HepMC3.h" #include "HepMC3/Relatives.h" #include "HepMC3/Reader.h" namespace Rivet{ namespace RivetHepMC = HepMC3; using RivetHepMC::ConstGenParticlePtr; using RivetHepMC::ConstGenVertexPtr; using RivetHepMC::Relatives; using RivetHepMC::ConstGenHeavyIonPtr; using HepMC_IO_type = RivetHepMC::Reader; using PdfInfo = RivetHepMC::GenPdfInfo; } #else #include "HepMC/GenEvent.h" #include "HepMC/GenParticle.h" #include "HepMC/HeavyIon.h" #include "HepMC/GenVertex.h" #include "HepMC/Version.h" #include "HepMC/GenRanges.h" #include "HepMC/IO_GenEvent.h" namespace Rivet{ namespace RivetHepMC = HepMC; // HepMC 2.07 provides its own #defines typedef const HepMC::GenParticle* ConstGenParticlePtr; typedef const HepMC::GenVertex* ConstGenVertexPtr; typedef const HepMC::HeavyIon* ConstGenHeavyIonPtr; /// @brief Replicated the HepMC3 Relatives syntax using HepMC2 IteratorRanges /// This is necessary mainly because of capitalisation differences class Relatives{ public: constexpr Relatives(HepMC::IteratorRange relo): _internal(relo){} constexpr HepMC::IteratorRange operator()() const {return _internal;} operator HepMC::IteratorRange() const {return _internal;} const static Relatives PARENTS; const static Relatives CHILDREN; const static Relatives ANCESTORS; const static Relatives DESCENDANTS; private: const HepMC::IteratorRange _internal; }; using HepMC_IO_type = HepMC::IO_GenEvent; using PdfInfo = RivetHepMC::PdfInfo; } #endif #include "Rivet/Tools/RivetSTL.hh" #include "Rivet/Tools/Exceptions.hh" namespace Rivet { using RivetHepMC::GenEvent; using ConstGenEventPtr = std::shared_ptr; /// @todo Use mcutils? namespace HepMCUtils{ ConstGenParticlePtr getParticlePtr(const RivetHepMC::GenParticle & gp); std::vector particles(ConstGenEventPtr ge); std::vector particles(const GenEvent *ge); std::vector vertices(ConstGenEventPtr ge); std::vector vertices(const GenEvent *ge); std::vector particles(ConstGenVertexPtr gv, const Relatives &relo); std::vector particles(ConstGenParticlePtr gp, const Relatives &relo); int uniqueId(ConstGenParticlePtr gp); int particles_size(ConstGenEventPtr ge); int particles_size(const GenEvent *ge); std::pair beams(const GenEvent *ge); std::shared_ptr makeReader(std::istream &istr, std::string * errm = 0); bool readEvent(std::shared_ptr io, std::shared_ptr evt); void strip(GenEvent & ge, const set & stripid = {1, -1, 2, -2, 3,-3, 21}); vector weightNames(const GenEvent & ge); pair crossSection(const GenEvent & ge); std::valarray weights(const GenEvent & ge); } } #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,716 +1,721 @@ #ifndef RIVET_RIVETYODA_HH #define RIVET_RIVETYODA_HH #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 #include namespace YODA { typedef std::shared_ptr AnalysisObjectPtr; typedef std::shared_ptr CounterPtr; typedef std::shared_ptr Histo1DPtr; typedef std::shared_ptr Histo2DPtr; typedef std::shared_ptr Profile1DPtr; typedef std::shared_ptr Profile2DPtr; typedef std::shared_ptr Scatter1DPtr; typedef std::shared_ptr Scatter2DPtr; typedef std::shared_ptr Scatter3DPtr; } namespace Rivet { class AnalysisObjectWrapper { public: virtual ~AnalysisObjectWrapper() {} 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; /// @todo Set active object for finalize virtual void setActiveFinalWeightIdx(unsigned int iWeight) = 0; virtual void unsetActiveWeight() = 0; bool operator ==(const AnalysisObjectWrapper& p) { return (this == &p); } bool operator !=(const AnalysisObjectWrapper& 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(p)); } bool operator!() const { return !_persistent; } explicit 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 _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(p)); } Scatter2DPtr() : _persistent() { } bool operator!() { return !_persistent; } explicit 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 _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(p)); } Scatter3DPtr() : _persistent() { } bool operator!() { return !_persistent; } explicit 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 _persistent; }; */ class MultiweightAOWrapper : public AnalysisObjectWrapper { public: using Inner = YODA::AnalysisObject; virtual void newSubEvent() = 0; virtual void pushToPersistent(const vector >& weight) = 0; virtual void pushToFinal() = 0; virtual YODA::AnalysisObjectPtr activeYODAPtr() const = 0; virtual string basePath() const = 0; }; using Weight = double; template using Fill = pair; template using Fills = multiset>; // 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 /// Wrappers for analysis objects to store all fills unaggregated, until collapsed template class TupleWrapper; template<> class TupleWrapper : public YODA::Counter { public: typedef shared_ptr> 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 ) { + (void)fraction; fills_.insert( {YODA::Counter::FillType(),weight} ); } void reset() { fills_.clear(); } const Fills & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Histo1D { public: typedef shared_ptr> 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 ) { + (void)fraction; if ( std::isnan(x) ) throw YODA::RangeError("X is NaN"); fills_.insert( { x , weight } ); } void reset() { fills_.clear(); } const Fills & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Profile1D { public: typedef shared_ptr> 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 ) { + (void)fraction; 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 & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Histo2D { public: typedef shared_ptr> 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 ) { + (void)fraction; 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 & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Profile2D { public: typedef shared_ptr> 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 ) { + (void)fraction; 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 & fills() const { return fills_; } private: // x / weight pairs Fills fills_; }; template<> class TupleWrapper : public YODA::Scatter1D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Scatter1D & h) : YODA::Scatter1D(h) {} }; template<> class TupleWrapper : public YODA::Scatter2D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Scatter2D & h) : YODA::Scatter2D(h) {} }; template<> class TupleWrapper : public YODA::Scatter3D { public: typedef shared_ptr> Ptr; TupleWrapper(const YODA::Scatter3D & h) : YODA::Scatter3D(h) {} }; template class Wrapper : public MultiweightAOWrapper { friend class Analysis; public: using Inner = T; /* @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(const vector& weightnames, 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 explicit operator bool() const { return static_cast(_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(); } // can be useful for weight analysis (see e.g. MC_WEIGHTS for use) T * _getPersistent (unsigned int iWeight) { return _persistent.at(iWeight).get(); } /* @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); } void setActiveFinalWeightIdx(unsigned int iWeight) { _active = _final.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 & persistent() const { return _persistent; } const vector & final() const { return _final; } /* to be implemented for each type */ void pushToPersistent(const vector >& weight); void pushToFinal(); /* M of these, one for each weight */ vector _persistent; /* This is the copy of _persistent that will be passed to finalize(). */ vector _final; /* N of these, one for each event in evgroup */ vector::Ptr> _evgroup; typename T::Ptr _active; string basePath() const { return _basePath; } string _basePath; // do we need implicit cast? // operator typename T::Ptr () { // return _active; // } friend class AnalysisHandler; }; /// We need our own shared_ptr class, so we can dispatch -> and * /// all the way down to the inner YODA analysis objects /// /// TODO: provide remaining functionality that shared_ptr has (not needed right now) /// template class rivet_shared_ptr { public: typedef T value_type; rivet_shared_ptr() = default; rivet_shared_ptr(decltype(nullptr)) : _p(nullptr) {} /// Convenience constructor, pass through to the Wrapper constructor rivet_shared_ptr(const vector& weightNames, const typename T::Inner & p) : _p( make_shared(weightNames, p) ) {} template rivet_shared_ptr(const shared_ptr & p) : _p(p) {} template rivet_shared_ptr(const rivet_shared_ptr & p) : _p(p.get()) {} // Goes right through to the active Wrapper object's members T & operator->() { return *_p; } const T & operator->() const { return *_p; } // The active YODA object typename T::Inner & operator*() { return **_p; } const typename T::Inner & operator*() const { return **_p; } bool operator!() const { return !_p || !(*_p); } explicit operator bool() const { return _p && bool(*_p); } template bool operator==(const rivet_shared_ptr & other) const { return _p == other._p; } template bool operator!=(const rivet_shared_ptr & other) const { return _p != other._p; } template bool operator<(const rivet_shared_ptr & other) const { return _p < other._p; } template bool operator>(const rivet_shared_ptr & other) const { return _p > other._p; } template bool operator<=(const rivet_shared_ptr & other) const { return _p <= other._p; } template bool operator>=(const rivet_shared_ptr & other) const { return _p >= other._p; } shared_ptr get() const { return _p; } private: shared_ptr _p; }; // every object listed here needs a virtual fill method in YODA, // otherwise the Tuple fakery won't work. using MultiweightAOPtr = rivet_shared_ptr; using Histo1DPtr = rivet_shared_ptr>; using Histo2DPtr = rivet_shared_ptr>; using Profile1DPtr = rivet_shared_ptr>; using Profile2DPtr = rivet_shared_ptr>; using CounterPtr = rivet_shared_ptr>; using Scatter1DPtr = rivet_shared_ptr>; using Scatter2DPtr = rivet_shared_ptr>; using Scatter3DPtr = rivet_shared_ptr>; 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 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); /// Traits class to access the type of the AnalysisObject in the /// reference files. template struct ReferenceTraits {}; template<> struct ReferenceTraits { typedef Counter RefT; }; template<> struct ReferenceTraits { typedef Scatter1D RefT; }; template<> struct ReferenceTraits { typedef Scatter2D RefT; }; template<> struct ReferenceTraits { typedef Scatter2D RefT; }; template<> struct ReferenceTraits { typedef Scatter2D RefT; }; template<> struct ReferenceTraits { typedef Scatter3D RefT; }; template<> struct ReferenceTraits { typedef Scatter3D RefT; }; template<> struct ReferenceTraits { typedef Scatter3D RefT; }; /// If @a dst and @a src both are of same subclass T, copy the /// contents of @a src into @a dst and return true. Otherwise return /// false. template inline bool aocopy(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst) { shared_ptr tsrc = dynamic_pointer_cast(src); if ( !tsrc ) return false; shared_ptr tdst = dynamic_pointer_cast(dst); if ( !tdst ) return false; *tdst = *tsrc; return true; } /// If @a dst and @a src both are of same subclass T, add the /// contents of @a src into @a dst and return true. Otherwise return /// false. template inline bool aoadd(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale) { shared_ptr tsrc = dynamic_pointer_cast(src); if ( !tsrc ) return false; shared_ptr tdst = dynamic_pointer_cast(dst); if ( !tdst ) return false; tsrc->scaleW(scale); *tdst += *tsrc; return true; } /// If @a dst is the same subclass as @a src, copy the contents of @a /// src into @a dst and return true. Otherwise return false. bool copyao(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst); /// If @a dst is the same subclass as @a src, scale the contents of /// @a src with @a scale and add it to @a dst and return true. Otherwise /// return false. bool addaos(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale); /// Check if two analysis objects have the same binning or, if not /// binned, are in other ways compatible. template inline bool bookingCompatible(TPtr a, TPtr b) { return a->sameBinning(*b); } inline bool bookingCompatible(CounterPtr, CounterPtr) { return true; } inline bool bookingCompatible(Scatter1DPtr a, Scatter1DPtr b) { return a->numPoints() == b->numPoints(); } inline bool bookingCompatible(Scatter2DPtr a, Scatter2DPtr b) { return a->numPoints() == b->numPoints(); } inline bool bookingCompatible(Scatter3DPtr a, Scatter3DPtr b) { return a->numPoints() == b->numPoints(); } inline bool bookingCompatible(YODA::CounterPtr, YODA::CounterPtr) { return true; } inline bool bookingCompatible(YODA::Scatter1DPtr a, YODA::Scatter1DPtr b) { return a->numPoints() == b->numPoints(); } inline bool bookingCompatible(YODA::Scatter2DPtr a, YODA::Scatter2DPtr b) { return a->numPoints() == b->numPoints(); } inline bool bookingCompatible(YODA::Scatter3DPtr a, YODA::Scatter3DPtr b) { return a->numPoints() == b->numPoints(); } /// class representing a YODA path with all its components. class AOPath { public: /// Constructor AOPath(string fullpath) : _valid(false), _path(fullpath), _raw(false), _tmp(false), _ref(false) { _valid = init(fullpath); } /// The full path. string path() const { return _path; } /// The analysis name. string analysis() const { return _analysis; } /// The analysis name with options. string analysisWithOptions() const { return _analysis + _optionstring; } /// The base name of the analysis object. string name() const { return _name; } /// The weight name. string weight() const { return _weight; } /// Is This a RAW (filling) object? bool isRaw() const { return _raw; } // Is This a temporary (filling) object? bool isTmp() const { return _tmp; } /// Is This a reference object? bool isRef() const { return _ref; } /// The string describing the options passed to the analysis. string optionString() const { return _optionstring; } /// Are there options passed to the analysis? bool hasOptions() const { return !_options.empty(); } /// Don't pass This optionto the analysis void removeOption(string opt) { _options.erase(opt); fixOptionString(); } /// Pass this option to the analysis. void setOption(string opt, string val) { _options[opt] = val; fixOptionString();} /// Was This option passed to the analyisi. bool hasOption(string opt) const { return _options.find(opt) != _options.end(); } /// Get the value of this option. string getOption(string opt) const { auto it = _options.find(opt); if ( it != _options.end() ) return it->second; return ""; } /// Reset the option string after changes; void fixOptionString(); /// Creat a full path (and set) for this. string mkPath() const; string setPath() { return _path = mkPath(); } /// Print out information void debug() const; /// Make this class ordered. bool operator<(const AOPath & other) const { return _path < other._path; } /// Check if path is valid. bool valid() const { return _valid; }; bool operator!() const { return !valid(); } private: /// Internal functions for disassembling a path name bool init(string fullpath); bool chopweight(string & fullpath); bool chopoptions(string & anal); bool _valid; string _path; string _analysis; string _optionstring; string _name; string _weight; bool _raw; bool _tmp; bool _ref; map _options; }; } #endif diff --git a/src/Projections/Beam.cc b/src/Projections/Beam.cc --- a/src/Projections/Beam.cc +++ b/src/Projections/Beam.cc @@ -1,147 +1,147 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/Beam.hh" namespace Rivet { ParticlePair beams(const Event& e) { // First try the official way: ask the GenEvent for the beam pointers assert(HepMCUtils::particles_size(e.genEvent()) >= 2); std::pair thebeams = HepMCUtils::beams(e.originalGenEvent()); if ( thebeams.first && thebeams.second ) { return ParticlePair{thebeams.first, thebeams.second}; } // Ok, that failed: let's find the status = 4 particles by hand const vector pstat4s = e.allParticles([](const Particle& p){ return p.genParticle()->status() == 4; }); if (pstat4s.size() >= 2) { return ParticlePair{pstat4s[0], pstat4s[1]}; } /// There are no barcodes in HepMC3 /// @todo implement some other fallback rubric? - #ifndef ENABLE_HEPMC_3 + #ifndef RIVET_ENABLE_HEPMC_3 // Hmm, this sucks. Last guess is that barcodes 1 and 2 are the beams if (e.genEvent()->barcode_to_particle(1) && e.genEvent()->barcode_to_particle(2)) { return ParticlePair{e.genEvent()->barcode_to_particle(1), e.genEvent()->barcode_to_particle(2)}; } #endif // Give up: return null beams return ParticlePair{Particle(), Particle()}; } double sqrtS(const FourMomentum& pa, const FourMomentum& pb) { const double mom1 = pa.pz(); const double e1 = pa.E(); const double mom2 = pb.pz(); const double e2 = pb.E(); double sqrts = sqrt( sqr(e1+e2) - sqr(mom1+mom2) ); return sqrts; } double asqrtS(const FourMomentum& pa, const FourMomentum& pb) { const static double MNUCLEON = 939*MeV; //< nominal nucleon mass return sqrtS(pa/(pa.mass()/MNUCLEON), pb/(pb.mass()/MNUCLEON)); } double asqrtS(const ParticlePair& beams) { return sqrtS(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second)); } FourMomentum acmsBoostVec(const FourMomentum& pa, const FourMomentum& pb) { const static double MNUCLEON = 939*MeV; //< nominal nucleon mass const double na = pa.mass()/MNUCLEON, nb = pb.mass()/MNUCLEON; return cmsBoostVec(pa/na, pb/nb); } FourMomentum acmsBoostVec(const ParticlePair& beams) { return cmsBoostVec(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second)); } Vector3 cmsBetaVec(const FourMomentum& pa, const FourMomentum& pb) { // const Vector3 rtn = (pa.p3() + pb.p3()) / (pa.E() + pb.E()); const Vector3 rtn = (pa + pb).betaVec(); return rtn; } Vector3 acmsBetaVec(const FourMomentum& pa, const FourMomentum& pb) { const static double MNUCLEON = 939*MeV; //< nominal nucleon mass const Vector3 rtn = cmsBetaVec(pa/(pa.mass()/MNUCLEON), pb/(pb.mass()/MNUCLEON)); return rtn; } Vector3 acmsBetaVec(const ParticlePair& beams) { const Vector3 rtn = cmsBetaVec(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second)); return rtn; } Vector3 cmsGammaVec(const FourMomentum& pa, const FourMomentum& pb) { // const Vector3 rtn = (pa + pb).gammaVec(); const double gamma = (pa.E() + pb.E()) / sqrt( sqr(pa.mass()) + sqr(pb.mass()) + 2*(pa.E()*pb.E() - dot(pa.p3(), pb.p3())) ); const Vector3 rtn = gamma * (pa.p3() + pb.p3()).unit(); return rtn; } Vector3 acmsGammaVec(const FourMomentum& pa, const FourMomentum& pb) { const static double MNUCLEON = 939*MeV; //< nominal nucleon mass Vector3 rtn = cmsGammaVec(pa/(pa.mass()/MNUCLEON), pb/(pb.mass()/MNUCLEON)); return rtn; } Vector3 acmsGammaVec(const ParticlePair& beams) { Vector3 rtn = cmsGammaVec(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second)); return rtn; } LorentzTransform cmsTransform(const FourMomentum& pa, const FourMomentum& pb) { /// @todo Automatically choose to construct from beta or gamma according to which is more precise? return LorentzTransform::mkFrameTransformFromGamma(cmsGammaVec(pa, pb)); } LorentzTransform acmsTransform(const FourMomentum& pa, const FourMomentum& pb) { /// @todo Automatically choose to construct from beta or gamma according to which is more precise? return LorentzTransform::mkFrameTransformFromGamma(acmsGammaVec(pa, pb)); } LorentzTransform acmsTransform(const ParticlePair& beams) { return LorentzTransform::mkFrameTransformFromGamma(acmsGammaVec(beams)); } ///////////////////////////////////////////// void Beam::project(const Event& e) { _theBeams = Rivet::beams(e); MSG_DEBUG("Beam particles = " << _theBeams << " => sqrt(s) = " << sqrtS()/GeV << " GeV"); } FourVector Beam::pv() const { RivetHepMC::FourVector v1, v2; const ParticlePair bpair = beams(); if (bpair.first.genParticle() && bpair.first.genParticle()->end_vertex()) v1 = bpair.first.genParticle()->end_vertex()->position(); if (bpair.second.genParticle() && bpair.second.genParticle()->end_vertex()) v2 = bpair.second.genParticle()->end_vertex()->position(); const FourVector rtn = (v1 == v2) ? FourVector(v1.t(), v1.x(), v1.y(), v1.z()) : FourVector(); MSG_DEBUG("Beam PV 4-position = " << rtn); return rtn; } } diff --git a/src/Projections/HepMCHeavyIon.cc b/src/Projections/HepMCHeavyIon.cc --- a/src/Projections/HepMCHeavyIon.cc +++ b/src/Projections/HepMCHeavyIon.cc @@ -1,104 +1,104 @@ // -*- C++ -*- #include "Rivet/Projections/HepMCHeavyIon.hh" -#ifdef ENABLE_HEPMC_3 +#ifdef RIVET_ENABLE_HEPMC_3 #define IMPLEMENTATION(rettype, functionname, defret) \ rettype HepMCHeavyIon::functionname() const { \ return _hi? _hi->functionname: defret; \ } #define IMPLEMENTATION_ALT_HEPMC2(rettype, functionname, hepmc2name, defret) \ rettype HepMCHeavyIon::functionname() const { \ return _hi? _hi->functionname: defret; \ } #define IMPLEMENTATION_NO_HEPMC2(rettype, functionname, defret) \ rettype HepMCHeavyIon::functionname() const { \ return _hi? _hi->functionname: defret; \ } #else #define IMPLEMENTATION(rettype, functionname, defret) \ rettype HepMCHeavyIon::functionname() const { \ return _hi? _hi->functionname(): defret; \ } #define IMPLEMENTATION_ALT_HEPMC2(rettype, functionname, hepmc2name, defret) \ rettype HepMCHeavyIon::functionname() const { \ return _hi? _hi->, hepmc2name(): defret; \ } #define IMPLEMENTATION_NO_HEPMC2(rettype, functionname, defret) \ rettype HepMCHeavyIon::functionname() const { \ MSG_WARNING("HeavyIon::" #functionname " is only avialable in HepMC3"); \ return defret; \ } #endif namespace Rivet { HepMCHeavyIon::HepMCHeavyIon() { setName("HepMCHeavyIon"); } void HepMCHeavyIon::project(const Event& e) { _hi = e.genEvent()->heavy_ion(); if ( !_hi ) MSG_WARNING("Could not find the HepMC HeavyIon object"); } IMPLEMENTATION(int, Ncoll_hard, -1) IMPLEMENTATION(int, Npart_proj, -1) IMPLEMENTATION(int, Npart_targ, -1) IMPLEMENTATION(int, Ncoll, -1) IMPLEMENTATION(int, N_Nwounded_collisions, -1) IMPLEMENTATION(int, Nwounded_N_collisions, -1) IMPLEMENTATION(int, Nwounded_Nwounded_collisions, -1) IMPLEMENTATION(double, impact_parameter, -1.0) IMPLEMENTATION(double, event_plane_angle, -1.0) IMPLEMENTATION(double, sigma_inel_NN, -1.0) IMPLEMENTATION_NO_HEPMC2(double, centrality, -1.0) IMPLEMENTATION_NO_HEPMC2(double, user_cent_estimate, -1.0) IMPLEMENTATION_NO_HEPMC2(int, Nspec_proj_n, -1) IMPLEMENTATION_NO_HEPMC2(int, Nspec_targ_n, -1) IMPLEMENTATION_NO_HEPMC2(int, Nspec_proj_p, -1) IMPLEMENTATION_NO_HEPMC2(int, Nspec_targ_p, -1) map HepMCHeavyIon::participant_plane_angles() const { -#ifdef ENABLE_HEPMC_3 +#ifdef RIVET_ENABLE_HEPMC_3 return _hi? _hi->participant_plane_angles: map(); #else MSG_WARNING("HeavyIon::participant_plane_angles is only avialable in HepMC3"); return map(); #endif } map HepMCHeavyIon::eccentricities() const { -#ifdef ENABLE_HEPMC_3 +#ifdef RIVET_ENABLE_HEPMC_3 return _hi? _hi->eccentricities: map(); #else MSG_WARNING("HeavyIon::eccentricities is only avialable in HepMC3"); return map(); #endif } }