diff --git a/pkgs/applications/science/physics/sherpa/default.nix b/pkgs/applications/science/physics/sherpa/default.nix
index 2dfe3f7d7bb2c..b1fbbab8781fb 100644
--- a/pkgs/applications/science/physics/sherpa/default.nix
+++ b/pkgs/applications/science/physics/sherpa/default.nix
@@ -1,12 +1,12 @@
-{ lib, stdenv, fetchurl, autoconf, gfortran, hepmc2, fastjet, lhapdf, rivet, sqlite }:
+{ lib, stdenv, fetchurl, autoconf, gfortran, hepmc3, fastjet, lhapdf, rivet, sqlite }:
 
 stdenv.mkDerivation rec {
   pname = "sherpa";
-  version = "2.2.15";
+  version = "2.2.16";
 
   src = fetchurl {
     url = "https://www.hepforge.org/archive/sherpa/SHERPA-MC-${version}.tar.gz";
-    sha256 = "sha256-3zvLa1k/bm7uOWKUsTyQM39cPBXJJlF1OgPgznl1hks=";
+    sha256 = "sha256-AntSN5BhtJFuDBoOFvrzoCr/W4SnX5CeAXiTcz9MjUs=";
   };
 
   postPatch = lib.optionalString (stdenv.hostPlatform.libc == "glibc") ''
@@ -21,7 +21,7 @@ stdenv.mkDerivation rec {
 
   configureFlags = [
     "--with-sqlite3=${sqlite.dev}"
-    "--enable-hepmc2=${hepmc2}"
+    "--enable-hepmc3=${hepmc3}"
     "--enable-fastjet=${fastjet}"
     "--enable-lhapdf=${lhapdf}"
     "--enable-rivet=${rivet}"
diff --git a/pkgs/development/libraries/physics/fastnlo-toolkit/default.nix b/pkgs/development/libraries/physics/fastnlo-toolkit/default.nix
index 3c781689bbf83..820afcbd8f68b 100644
--- a/pkgs/development/libraries/physics/fastnlo-toolkit/default.nix
+++ b/pkgs/development/libraries/physics/fastnlo-toolkit/default.nix
@@ -22,6 +22,11 @@ stdenv.mkDerivation rec {
     sha256 = "sha256-7aIMYCOkHC/17CHYiEfrxvtSJxTDivrS7BQ32cGiEy0=";
   };
 
+  patches = [
+    # Compatibility with YODA 2.x
+    ./yoda2_support.patch
+  ];
+
   buildInputs = [
     boost
     gfortran
diff --git a/pkgs/development/libraries/physics/fastnlo-toolkit/yoda2_support.patch b/pkgs/development/libraries/physics/fastnlo-toolkit/yoda2_support.patch
new file mode 100644
index 0000000000000..d9f4aaa5159f8
--- /dev/null
+++ b/pkgs/development/libraries/physics/fastnlo-toolkit/yoda2_support.patch
@@ -0,0 +1,247 @@
+diff --git a/src/fnlo-tk-statunc.cc b/src/fnlo-tk-statunc.cc
+index 62d1eec..ad62cac 100644
+--- a/src/fnlo-tk-statunc.cc
++++ b/src/fnlo-tk-statunc.cc
+@@ -30,7 +30,16 @@
+ #include "fastnlotk/fastNLOLHAPDF.h"
+ #include "fastnlotk/speaker.h"
+ #ifdef WITH_YODA
++#if defined __has_include
++#if !__has_include("YODA/WriterAIDA.h")
++#define WITH_YODA2
++#endif
++#endif
++#ifdef WITH_YODA2
++#include "YODA/Scatter.h"
++#else
+ #include "YODA/Scatter2D.h"
++#endif
+ #include "YODA/WriterYODA.h"
+ #endif
+ 
+@@ -493,28 +502,25 @@ int main(int argc, char** argv) {
+ 
+    //! --- 1D
+    if (NDim == 1) {
+-      //! Vectors to fill 2D scatter plot
+-      vector < double > x;
+-      vector < double > y;
+-      vector < double > exminus;
+-      vector < double > explus;
+-      vector < double > eyminus;
+-      vector < double > eyplus;
++      //! Vector to fill 2D scatter plot
++      vector<YODA::Point2D> points;
+       //! Loop over bins in outer (1st) dimension
+       for (unsigned int k =0 ; k<NDimBins[0] ; k++) {
+-         x.push_back((bins[iobs].second + bins[iobs].first)/2.0);
+-         explus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
+-         exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
+-         y.push_back(xs[iobs]);
+-         eyplus.push_back(dxsu[iobs]);
+-         eyminus.push_back(std::abs(dxsl[iobs]));
++         points.emplace_back(
++            (bins[iobs].second + bins[iobs].first)/2.0,
++            xs[iobs],
++            (bins[iobs].second - bins[iobs].first)/2.0,
++            (bins[iobs].second - bins[iobs].first)/2.0,
++            std::abs(dxsl[iobs]),
++            dxsu[iobs]
++         );
+          iobs++;
+       }
+       stringstream plotno;                                                                         // To make i+1 from int
+       plotno << offset;                                                                            // to a string for the naming
+       //      RivetId.replace( capital_pos +3 - plotno.str().size(), plotno.str().size(), plotno.str());   // Next plot name
+       // Pointer in order not to be deleted after we exit the loop, so we can then save them into the yoda file
+-      YODA::Scatter2D * plot = new YODA::Scatter2D(x,y,exminus,explus,eyminus,eyplus,"/" + RivetId,LineName);
++      YODA::Scatter2D * plot = new YODA::Scatter2D(points, "/" + RivetId, LineName);
+       // Insert the plot pointer into the vector of analysis object pointers
+       aos.push_back(plot);
+    }
+@@ -522,29 +528,26 @@ int main(int argc, char** argv) {
+    else if (NDim == 2) {
+       //! Loop over bins in outer (1st) dimension
+       for (unsigned int j=0; j<NDimBins[0]; j++) {
+-         //! Vectors to fill 2D scatter plot
+-         vector < double > x;
+-         vector < double > y;
+-         vector < double > exminus;
+-         vector < double > explus;
+-         vector < double > eyminus;
+-         vector < double > eyplus;
++         //! Vector to fill 2D scatter plot
++         vector<YODA::Point2D> points;
+          //! Loop over bins in inner (2nd) dimension
+          NDimBins[1] = fnlo.GetNDim1Bins(j);
+          for (unsigned int k = 0; k<NDimBins[1]; k++) {
+-            x.push_back((bins[iobs].second + bins[iobs].first)/2.0);
+-            explus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
+-            exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
+-            y.push_back(xs[iobs]);
+-            eyplus.push_back(dxsu[iobs]);
+-            eyminus.push_back(std::abs(dxsl[iobs]));
++            points.emplace_back(
++               (bins[iobs].second + bins[iobs].first)/2.0,
++               xs[iobs],
++               (bins[iobs].second - bins[iobs].first)/2.0,
++               (bins[iobs].second - bins[iobs].first)/2.0,
++               std::abs(dxsl[iobs]),
++               dxsu[iobs]
++            );
+             iobs++;
+          }
+          stringstream plotno;                                                                         // To make i+1 from int
+          plotno << offset+j;                                                                          // to a string for the naming
+          RivetId.replace( capital_pos +3 - plotno.str().size(), plotno.str().size(), plotno.str());   // Next plot name
+          // Pointer in order not to be deleted after we exit the loop, so we can then save them into the yoda file
+-         YODA::Scatter2D * plot = new YODA::Scatter2D(x,y,exminus,explus,eyminus,eyplus,"/" + RivetId,LineName);
++         YODA::Scatter2D * plot = new YODA::Scatter2D(points, "/" + RivetId, LineName);
+          // Insert the plot pointer into the vector of analysis object pointers
+          aos.push_back(plot);
+       }
+diff --git a/src/fnlo-tk-yodaout.cc b/src/fnlo-tk-yodaout.cc
+index 9aa8ac3..dabe924 100644
+--- a/src/fnlo-tk-yodaout.cc
++++ b/src/fnlo-tk-yodaout.cc
+@@ -25,7 +25,16 @@
+ #include "fastnlotk/fastNLOTools.h"
+ #include "fastnlotk/speaker.h"
+ #ifdef WITH_YODA
++#if defined __has_include
++#if !__has_include("YODA/WriterAIDA.h")
++#define WITH_YODA2
++#endif
++#endif
++#ifdef WITH_YODA2
++#include "YODA/Scatter.h"
++#else
+ #include "YODA/Scatter2D.h"
++#endif
+ #include "YODA/WriterYODA.h"
+ #endif
+ 
+@@ -548,26 +557,23 @@ int main(int argc, char** argv) {
+ 
+    //! --- 1D
+    if (NDim == 1) {
+-      //! Vectors to fill 2D scatter plot
+-      vector < double > x;
+-      vector < double > y;
+-      vector < double > exminus;
+-      vector < double > explus;
+-      vector < double > eyminus;
+-      vector < double > eyplus;
++      //! Vector to fill 2D scatter plot
++      vector < YODA::Point2D > points;
+       //! Loop over bins in outer (1st) dimension
+       for (unsigned int k =0 ; k<NDimBins[0] ; k++) {
+-         x.push_back((bins[iobs].second + bins[iobs].first)/2.0);
+-         explus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
+-         exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
+-         y.push_back(xs[iobs]);
+-         eyplus.push_back(dxsu[iobs]);
+-         eyminus.push_back(std::abs(dxsl[iobs]));
++         points.emplace_back(
++           (bins[iobs].second + bins[iobs].first)/2.0,
++           xs[iobs],
++           (bins[iobs].second - bins[iobs].first)/2.0,
++           (bins[iobs].second - bins[iobs].first)/2.0,
++           std::abs(dxsl[iobs]),
++           dxsu[iobs]
++         );
+          iobs++;
+       }
+ #ifdef WITH_YODA
+       /// Pointer in order not to be deleted after we exit the loop, so we can then save them into the yoda file
+-      YODA::Scatter2D * plot = new YODA::Scatter2D(x,y,exminus,explus,eyminus,eyplus,"/" + RivetId,LineName);
++      YODA::Scatter2D * plot = new YODA::Scatter2D(points, "/" + RivetId, LineName);
+       /// Insert the plot pointer into the vector of analysis object pointers
+       aos.push_back(plot);
+ #endif
+@@ -577,23 +583,20 @@ int main(int argc, char** argv) {
+       //! Loop over bins in outer (1st) dimension
+       int nhist = 0;
+       for (unsigned int j=0; j<NDimBins[0]; j++) {
+-         //! Vectors to fill 2D scatter plot
++         //! Vector to fill 2D scatter plot
+          nhist++;
+-         vector < double > x;
+-         vector < double > y;
+-         vector < double > exminus;
+-         vector < double > explus;
+-         vector < double > eyminus;
+-         vector < double > eyplus;
++         vector < YODA::Point2D > points;
+          //! Loop over bins in inner (2nd) dimension
+          NDimBins[1] = fnlo->GetNDim1Bins(j);
+          for (unsigned int k = 0; k<NDimBins[1]; k++) {
+-            x.push_back((bins[iobs].second + bins[iobs].first)/2.0);
+-            explus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
+-            exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
+-            y.push_back(xs[iobs]);
+-            eyplus.push_back(dxsu[iobs]);
+-            eyminus.push_back(std::abs(dxsl[iobs]));
++            points.emplace_back(
++              (bins[iobs].second + bins[iobs].first)/2.0,
++              xs[iobs],
++              (bins[iobs].second - bins[iobs].first)/2.0,
++              (bins[iobs].second - bins[iobs].first)/2.0,
++              std::abs(dxsl[iobs]),
++              dxsu[iobs]
++            );
+             iobs++;
+          }
+          /// Derive histogram counter
+@@ -610,7 +613,7 @@ int main(int argc, char** argv) {
+          RivetId.replace(capital_pos +3 - histno.str().size(), histno.str().size(), histno.str());
+ #ifdef WITH_YODA
+          /// Pointer in order not to be deleted after we exit the loop, so we can then save the plots into the yoda file
+-         YODA::Scatter2D * plot = new YODA::Scatter2D(x,y,exminus,explus,eyminus,eyplus,"/" + RivetId,LineName);
++         YODA::Scatter2D * plot = new YODA::Scatter2D(points, "/" + RivetId, LineName);
+          /// Insert the plot pointer into the vector of analysis object pointers
+          aos.push_back(plot);
+ #endif
+@@ -624,23 +627,20 @@ int main(int argc, char** argv) {
+          //! Loop over bins in middle (2nd) dimension
+          NDimBins[1] = fnlo->GetNDim1Bins(j);
+          for (unsigned int k = 0; k<NDimBins[1]; k++) {
+-            //! Vectors to fill 2D scatter plot
++            //! Vector to fill 2D scatter plot
+             nhist++;
+-            vector < double > x;
+-            vector < double > y;
+-            vector < double > exminus;
+-            vector < double > explus;
+-            vector < double > eyminus;
+-            vector < double > eyplus;
++            vector < YODA::Point2D > points;
+             //! Loop over bins in inner (3rd) dimension
+             NDimBins[2] = fnlo->GetNDim2Bins(j,k);
+             for (unsigned int l = 0; l<NDimBins[2]; l++) {
+-               x.push_back((bins[iobs].second + bins[iobs].first)/2.0);
+-               explus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
+-               exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0);
+-               y.push_back(xs[iobs]);
+-               eyplus.push_back(dxsu[iobs]);
+-               eyminus.push_back(std::abs(dxsl[iobs]));
++               points.emplace_back(
++                  (bins[iobs].second + bins[iobs].first)/2.0,
++                  xs[iobs],
++                  (bins[iobs].second - bins[iobs].first)/2.0,
++                  (bins[iobs].second - bins[iobs].first)/2.0,
++                  std::abs(dxsl[iobs]),
++                  dxsu[iobs]
++               );
+                iobs++;
+             }
+             /// Derive histogram counter
+@@ -657,7 +657,7 @@ int main(int argc, char** argv) {
+             RivetId.replace(capital_pos +3 - histno.str().size(), histno.str().size(), histno.str());
+ #ifdef WITH_YODA
+             /// Pointer in order not to be deleted after we exit the loop, so we can then save the plots into the yoda file
+-            YODA::Scatter2D * plot = new YODA::Scatter2D(x,y,exminus,explus,eyminus,eyplus,"/" + RivetId,LineName);
++            YODA::Scatter2D * plot = new YODA::Scatter2D(points, "/" + RivetId, LineName);
+             /// Insert the plot pointer into the vector of analysis object pointers
+             aos.push_back(plot);
+ #endif
diff --git a/pkgs/development/libraries/physics/rivet/default.nix b/pkgs/development/libraries/physics/rivet/default.nix
index c8db0fb22d0ee..3246d9e06ec38 100644
--- a/pkgs/development/libraries/physics/rivet/default.nix
+++ b/pkgs/development/libraries/physics/rivet/default.nix
@@ -1,12 +1,12 @@
-{ lib, stdenv, fetchurl, fastjet, fastjet-contrib, ghostscript, hepmc, imagemagick, less, python3, rsync, texliveBasic, yoda, which, makeWrapper }:
+{ lib, stdenv, fetchurl, fastjet, fastjet-contrib, ghostscript, hdf5, hepmc3, highfive, imagemagick, less, pkg-config, python3, rsync, texliveBasic, yoda, which, makeWrapper }:
 
 stdenv.mkDerivation rec {
   pname = "rivet";
-  version = "3.1.10";
+  version = "4.0.1";
 
   src = fetchurl {
     url = "https://www.hepforge.org/archive/rivet/Rivet-${version}.tar.bz2";
-    hash = "sha256-RYuODfHec46ZctJLJg6qCH3xLJnU/p3uU3fUfqakmRk=";
+    hash = "sha256-ToaS1uilOWHHeYPra6SJPDdlzyP3BXieTYZb5Iku/3k=";
   };
 
   latex = texliveBasic.withPackages (ps: with ps; [
@@ -25,43 +25,43 @@ stdenv.mkDerivation rec {
     xstring
   ]);
 
-  nativeBuildInputs = [ rsync makeWrapper ];
-  buildInputs = [ hepmc imagemagick python3 latex python3.pkgs.yoda ];
-  propagatedBuildInputs = [ fastjet fastjet-contrib ];
+  nativeBuildInputs = [ rsync makeWrapper pkg-config ];
+  buildInputs = [ hepmc3 highfive imagemagick python3 latex python3.pkgs.yoda ];
+  propagatedBuildInputs = [ hdf5 fastjet fastjet-contrib ];
 
   preConfigure = ''
+    substituteInPlace configure \
+      --replace-fail 'if test $HEPMC_VERSION -le 310; then' 'if false; then'
     substituteInPlace bin/rivet-build.in \
-      --replace 'num_jobs=$(getconf _NPROCESSORS_ONLN)' 'num_jobs=''${NIX_BUILD_CORES:-$(getconf _NPROCESSORS_ONLN)}' \
-      --replace 'which' '"${which}/bin/which"' \
-      --replace 'mycxx=' 'mycxx=${stdenv.cc}/bin/${if stdenv.cc.isClang or false then "clang++" else "g++"}  #' \
-      --replace 'mycxxflags="' "mycxxflags=\"$NIX_CFLAGS_COMPILE $NIX_CXXSTDLIB_COMPILE $NIX_CFLAGS_LINK "
+      --replace-fail 'num_jobs=$(getconf _NPROCESSORS_ONLN)' 'num_jobs=''${NIX_BUILD_CORES:-$(getconf _NPROCESSORS_ONLN)}' \
+      --replace-fail 'which' '"${which}/bin/which"' \
+      --replace-fail 'mycxx=' 'mycxx=${stdenv.cc}/bin/${if stdenv.cc.isClang or false then "clang++" else "g++"}  #' \
+      --replace-fail 'mycxxflags="' "mycxxflags=\"$NIX_CFLAGS_COMPILE $NIX_CXXSTDLIB_COMPILE $NIX_CFLAGS_LINK "
   '';
 
   preInstall = ''
     substituteInPlace bin/make-plots \
-      --replace '"which"' '"${which}/bin/which"' \
-      --replace '"latex"' '"'$latex'/bin/latex"' \
-      --replace '"dvips"' '"'$latex'/bin/dvips"' \
-      --replace '"ps2pdf"' '"${ghostscript}/bin/ps2pdf"' \
-      --replace '"ps2eps"' '"${ghostscript}/bin/ps2eps"' \
-      --replace '"kpsewhich"' '"'$latex'/bin/kpsewhich"' \
-      --replace '"convert"' '"${imagemagick.out}/bin/convert"'
+      --replace-fail '"which"' '"${which}/bin/which"' \
+      --replace-fail '"latex"' '"'$latex'/bin/latex"' \
+      --replace-fail '"dvips"' '"'$latex'/bin/dvips"' \
+      --replace-fail '"ps2pdf"' '"${ghostscript}/bin/ps2pdf"' \
+      --replace-fail '"ps2eps"' '"${ghostscript}/bin/ps2eps"' \
+      --replace-fail '"kpsewhich"' '"'$latex'/bin/kpsewhich"' \
+      --replace-fail '"convert"' '"${imagemagick.out}/bin/convert"'
     substituteInPlace bin/rivet \
-      --replace '"less"' '"${less}/bin/less"'
-    substituteInPlace bin/rivet-mkhtml \
-      --replace '"make-plots"' \"$out/bin/make-plots\" \
-      --replace '"rivet-cmphistos"' \"$out/bin/rivet-cmphistos\" \
-      --replace 'ch_cmd = [sys.executable, os.path.join(os.path.dirname(__file__),' 'ch_cmd = [('
+      --replace-fail '"less"' '"${less}/bin/less"'
+    substituteInPlace bin/rivet-mkhtml-tex \
+      --replace-fail '"make-plots"' \"$out/bin/make-plots\" \
+      --replace-fail '"rivet-cmphistos"' \"$out/bin/rivet-cmphistos\" \
+      --replace-fail 'ch_cmd = [sys.executable, os.path.join(os.path.dirname(__file__),' 'ch_cmd = [('
   '';
 
   configureFlags = [
     "--with-fastjet=${fastjet}"
     "--with-yoda=${yoda}"
-  ] ++ (if lib.versions.major hepmc.version == "3" then [
-    "--with-hepmc3=${hepmc}"
-  ] else [
-    "--with-hepmc=${hepmc}"
-  ]);
+    "--with-hepmc3=${hepmc3}"
+    "--with-highfive=${highfive}"
+  ];
 
   enableParallelBuilding = true;
 
diff --git a/pkgs/development/libraries/physics/thepeg/default.nix b/pkgs/development/libraries/physics/thepeg/default.nix
index b55d6869e4dac..d4274b9dbed49 100644
--- a/pkgs/development/libraries/physics/thepeg/default.nix
+++ b/pkgs/development/libraries/physics/thepeg/default.nix
@@ -1,4 +1,4 @@
-{ lib, stdenv, fetchurl, autoreconfHook, boost, fastjet, gsl, hepmc2, lhapdf, rivet, zlib }:
+{ lib, stdenv, fetchurl, fetchpatch, autoreconfHook, boost, fastjet, gsl, hepmc3, lhapdf, rivet, zlib }:
 
 stdenv.mkDerivation rec {
   pname = "thepeg";
@@ -9,12 +9,21 @@ stdenv.mkDerivation rec {
     hash = "sha256-rDWXmuicKWCMqSwVakn/aKrOeloSoMkvCgGoM9LTRXI=";
   };
 
+  patches = [
+    # Rivet 4 support
+    (fetchpatch {
+      url = "https://github.com/hep-mirrors/thepeg/commit/d974704fe48876c13cb7f544cabcf6ef30c00f48.diff";
+      hash = "sha256-HzyNigbhWzSpjvvYw3+LZvnrSoV6Pmzghw/5VY5nlqk=";
+    })
+  ];
+
   nativeBuildInputs = [ autoreconfHook ];
 
-  buildInputs = [ boost fastjet gsl hepmc2 lhapdf rivet zlib ];
+  buildInputs = [ boost fastjet gsl hepmc3 lhapdf rivet zlib ];
 
   configureFlags = [
-    "--with-hepmc=${hepmc2}"
+    "--with-hepmc=${hepmc3}"
+    "--with-hepmcversion=3"
     "--with-rivet=${rivet}"
     "--without-javagui"
   ];
diff --git a/pkgs/development/libraries/physics/yoda/default.nix b/pkgs/development/libraries/physics/yoda/default.nix
index a6ad9cbd92744..1b478e368c8d5 100644
--- a/pkgs/development/libraries/physics/yoda/default.nix
+++ b/pkgs/development/libraries/physics/yoda/default.nix
@@ -1,6 +1,7 @@
 { lib
 , stdenv
-, fetchurl
+, fetchFromGitLab
+, autoreconfHook
 , python
 , root
 , makeWrapper
@@ -10,14 +11,17 @@
 
 stdenv.mkDerivation rec {
   pname = "yoda";
-  version = "1.9.10";
+  version = "2.0.1";
 
-  src = fetchurl {
-    url = "https://www.hepforge.org/archive/yoda/YODA-${version}.tar.bz2";
-    hash = "sha256-CnCO6dcElF0zh8xDexX/3fOCxw/lurOe0r2/g8LCjG8=";
+  src = fetchFromGitLab {
+    owner = "hepcedar";
+    repo = pname;
+    rev = "yoda-${version}";
+    hash = "sha256-WeMOd1dsy+qwiODQ7EWiLa1GNfUWXwpTtkS4EKUWF+c=";
   };
 
   nativeBuildInputs = with python.pkgs; [
+    autoreconfHook
     cython
     makeWrapper
   ];
diff --git a/pkgs/top-level/all-packages.nix b/pkgs/top-level/all-packages.nix
index f3a9536c70ec7..d5cb22997c571 100644
--- a/pkgs/top-level/all-packages.nix
+++ b/pkgs/top-level/all-packages.nix
@@ -36814,7 +36814,6 @@ with pkgs;
   };
 
   rivet = callPackage ../development/libraries/physics/rivet {
-    hepmc = hepmc2;
     imagemagick = graphicsmagick-imagemagick-compat;
   };