diff --git a/v2.6/toolkit/src/fnlo-tk-statunc.cc b/v2.6/toolkit/src/fnlo-tk-statunc.cc
index 751fb5a71904271df0a2e44cac0d6dbac54eccda..5c8ace158eed905fe1733d207ce35a118ab3ab8f 100644
--- a/v2.6/toolkit/src/fnlo-tk-statunc.cc
+++ b/v2.6/toolkit/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
 
@@ -495,28 +504,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);
    }
@@ -524,29 +530,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/v2.6/toolkit/src/fnlo-tk-yodaout.cc b/v2.6/toolkit/src/fnlo-tk-yodaout.cc
index 23b6d571c36a979609e2ae66883afe2f00d92b31..de4381573dadc06395b42d9dd169a1ddc0b2df54 100644
--- a/v2.6/toolkit/src/fnlo-tk-yodaout.cc
+++ b/v2.6/toolkit/src/fnlo-tk-yodaout.cc
@@ -29,7 +29,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
 
@@ -598,26 +607,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*xrescale);
-         explus.push_back((bins[iobs].second - bins[iobs].first)/2.0*xrescale);
-         exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0*xrescale);
-         y.push_back(XsUnc.xs[iobs]);
-         eyplus.push_back(XsUnc.xs[iobs]*XsUnc.dxsu[iobs]*dxsrescale);
-         eyminus.push_back(std::abs(XsUnc.xs[iobs]*XsUnc.dxsl[iobs])*dxsrescale);
+         points.emplace_back(
+            (bins[iobs].second + bins[iobs].first)/2.0*xrescale,
+            XsUnc.xs[iobs],
+            (bins[iobs].second - bins[iobs].first)/2.0*xrescale,
+            (bins[iobs].second - bins[iobs].first)/2.0*xrescale,
+            std::abs(XsUnc.xs[iobs]*XsUnc.dxsl[iobs])*dxsrescale,
+            XsUnc.xs[iobs]*XsUnc.dxsu[iobs]*dxsrescale
+         );
          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
@@ -627,23 +633,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*xrescale);
-            explus.push_back((bins[iobs].second - bins[iobs].first)/2.0*xrescale);
-            exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0*xrescale);
-            y.push_back(XsUnc.xs[iobs]);
-            eyplus.push_back(XsUnc.xs[iobs]*XsUnc.dxsu[iobs]*dxsrescale);
-            eyminus.push_back(std::abs(XsUnc.xs[iobs]*XsUnc.dxsl[iobs])*dxsrescale);
+            points.emplace_back(
+               (bins[iobs].second + bins[iobs].first)/2.0*xrescale,
+               XsUnc.xs[iobs],
+               (bins[iobs].second - bins[iobs].first)/2.0*xrescale,
+               (bins[iobs].second - bins[iobs].first)/2.0*xrescale,
+               std::abs(XsUnc.xs[iobs]*XsUnc.dxsl[iobs])*dxsrescale,
+               XsUnc.xs[iobs]*XsUnc.dxsu[iobs]*dxsrescale
+            );
             iobs++;
          }
          /// Derive histogram counter
@@ -660,7 +663,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
@@ -674,23 +677,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*xrescale);
-               explus.push_back((bins[iobs].second - bins[iobs].first)/2.0*xrescale);
-               exminus.push_back((bins[iobs].second - bins[iobs].first)/2.0*xrescale);
-               y.push_back(XsUnc.xs[iobs]);
-               eyplus.push_back(XsUnc.xs[iobs]*XsUnc.dxsu[iobs]*dxsrescale);
-               eyminus.push_back(std::abs(XsUnc.xs[iobs]*XsUnc.dxsl[iobs])*dxsrescale);
+               points.emplace_back(
+                 (bins[iobs].second + bins[iobs].first)/2.0*xrescale,
+                 XsUnc.xs[iobs],
+                 (bins[iobs].second - bins[iobs].first)/2.0*xrescale,
+                 (bins[iobs].second - bins[iobs].first)/2.0*xrescale,
+                 std::abs(XsUnc.xs[iobs]*XsUnc.dxsl[iobs])*dxsrescale,
+                 XsUnc.xs[iobs]*XsUnc.dxsu[iobs]*dxsrescale
+               );
                iobs++;
             }
             /// Derive histogram counter
@@ -707,7 +707,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