diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc
index db2d620fa0185668b31f202bdddf99fe8537724e..42deea421e610eedda4756e0c2de27b86ad26b0e 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc
@@ -352,6 +352,10 @@ bool fastNLOCoeffAddFix::IsEquivalent(const fastNLOCoeffBase& other, double rtol
       return false;
    }
 
+   if (NPDFDim != op->GetNPDFDim()) {
+      debug["IsEquivalent"] << "NPDFDim is different: " << NPDFDim << " <-> " << op->GetNPDFDim() << endl;
+      return false;
+   }
    if (!fastNLOTools::SameTails(GetAllXNodes1(), op->GetAllXNodes1(), rtol)) {
       debug["IsEquivalent"] << "XNode1 not equivalent, see above." << endl;
       return false;
@@ -366,11 +370,13 @@ bool fastNLOCoeffAddFix::IsEquivalent(const fastNLOCoeffBase& other, double rtol
       return false;
    }
    for (unsigned int obsBin = 0; obsBin < SigmaTilde.size(); obsBin++) {
-      unsigned int tOffsetXN1;
-      unsigned int oOffsetXN1;
+      unsigned int tOffsetXN1, oOffsetXN1, tOffsetXN2, oOffsetXN2;
       int tnb1 = GetNxtot1(obsBin);
       int onb1 = op->GetNxtot1(obsBin);
-      int xMax = std::min(tnb1, onb1);
+      int tnb2 = GetNxtot2(obsBin);
+      int onb2 = op->GetNxtot2(obsBin);
+      int x1Max = std::min(tnb1, onb1);
+      int x2Max = std::min(tnb2, onb2);
       if (tnb1 > onb1) {
          tOffsetXN1 = tnb1 - onb1;
          oOffsetXN1 = 0;
@@ -378,6 +384,18 @@ bool fastNLOCoeffAddFix::IsEquivalent(const fastNLOCoeffBase& other, double rtol
          tOffsetXN1 = 0;
          oOffsetXN1 = onb1 - tnb1;
       }
+      if (NPDFDim > 1) {
+         if (tnb2 > onb2) {
+            tOffsetXN2 = tnb2 - onb2;
+            oOffsetXN2 = 0;
+         } else {
+            tOffsetXN2 = 0;
+            oOffsetXN2 = onb2 - tnb2;
+         }
+      } else {
+         tOffsetXN2 = tOffsetXN1;
+         oOffsetXN2 = oOffsetXN1;
+      }
 
       fastNLO::v4d tst4 = SigmaTilde[obsBin];
       fastNLO::v4d ost4 = ost5[obsBin];
@@ -393,10 +411,12 @@ bool fastNLOCoeffAddFix::IsEquivalent(const fastNLOCoeffBase& other, double rtol
          for (unsigned int scaleNode = 0; scaleNode < tst3.size(); scaleNode++) {
             fastNLO::v2d tst2 = tst3[scaleNode];
             fastNLO::v2d ost2 = ost3[scaleNode];
-            for (int xHigh = 0; xHigh < xMax; xHigh++) {
-               for (int xLow = 0; xLow <= xHigh; xLow++) {
-                  fastNLO::v1d tst1 = tst2[GetXIndex(obsBin, xLow + tOffsetXN1, xHigh + tOffsetXN1)];
-                  fastNLO::v1d ost1 = ost2[op->GetXIndex(obsBin, xLow + oOffsetXN1, xHigh + oOffsetXN1)];
+            int xIt2Max = NPDFDim > 1 ? x2Max : x1Max;
+            for (int xIt2 = 0; xIt2 < xIt2Max; xIt2++) {
+               int xIt1Max = NPDFDim > 1 ? x1Max : xIt2 + 1;
+               for (int xIt1 = 0; xIt1 < xIt1Max; xIt1++) {
+                  fastNLO::v1d tst1 = tst2[GetXIndex(obsBin, xIt1 + tOffsetXN1, xIt2 + tOffsetXN2)];
+                  fastNLO::v1d ost1 = ost2[op->GetXIndex(obsBin, xIt1 + oOffsetXN1, xIt2 + oOffsetXN2)];
                   if (tst1.size() != ost1.size()) {
                      return false;
                   }
@@ -404,7 +424,7 @@ bool fastNLOCoeffAddFix::IsEquivalent(const fastNLOCoeffBase& other, double rtol
                      double rdiff = std::abs((tst1[subProcess] - ost1[subProcess]) / tst1[subProcess]);
                      if (rdiff > rtol) {
                         debug["IsEquivalent"] << "SigmaTilde rdiff too high: obsBin=" << obsBin << " scaleVar=" << scaleVar
-                           << " scaleNode=" << scaleNode << " xHigh=" << (xHigh - xMax) << " xLow=" << (xLow - xMax)
+                           << " scaleNode=" << scaleNode << " x1=" << (xIt1 - x1Max) << " x2=" << (xIt2 - x2Max)
                            << " subProcess=" << subProcess << " t=" << tst1[subProcess] << " o=" << ost1[subProcess]
                            << " rdiff=" << rdiff << endl;
                         return false;
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFlex.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFlex.cc
index 47637be9e0ffd341658bad0e6cbd769721486b49..7c4717d57656d832d204a6cc78789a3ce3fe52d9 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFlex.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFlex.cc
@@ -612,17 +612,21 @@ bool fastNLOCoeffAddFlex::IsEquivalent(const fastNLOCoeffBase& other, double rto
 //________________________________________________________________________________________________________________ //
 bool fastNLOCoeffAddFlex::IsSigmaTildeEquivalent(const fastNLOCoeffAddFlex* op, const fastNLO::v5d *tst5, const fastNLO::v5d *ost5, double rtol, std::string name) const {
    if (tst5->size() != ost5->size()) {
-      debug["IsEquivalent"] << name << ": Number of observable bins is different." << endl;
+      debug["IsSigmaTildeEquivalent"] << name << ": Number of observable bins is different." << endl;
+      return false;
+   }
+   if (NPDFDim != op->GetNPDFDim()) {
+      debug["IsSigmaTildeEquivalent"] << name << ": NPDFDim is different: " << NPDFDim << " <-> " << op->GetNPDFDim() << endl;
       return false;
    }
    for (unsigned int obsBin = 0; obsBin < tst5->size(); obsBin++) {
-      fastNLO::v4d tst4 = tst5->at(obsBin);
-      fastNLO::v4d ost4 = ost5->at(obsBin);
-      unsigned int tOffsetXN1;
-      unsigned int oOffsetXN1;
+      unsigned int tOffsetXN1, oOffsetXN1, tOffsetXN2, oOffsetXN2;
       int tnb1 = GetNxtot1(obsBin);
       int onb1 = op->GetNxtot1(obsBin);
-      int xMax = std::min(tnb1, onb1);
+      int tnb2 = GetNxtot2(obsBin);
+      int onb2 = op->GetNxtot2(obsBin);
+      int x1Max = std::min(tnb1, onb1);
+      int x2Max = std::min(tnb2, onb2);
       if (tnb1 > onb1) {
          tOffsetXN1 = tnb1 - onb1;
          oOffsetXN1 = 0;
@@ -630,10 +634,28 @@ bool fastNLOCoeffAddFlex::IsSigmaTildeEquivalent(const fastNLOCoeffAddFlex* op,
          tOffsetXN1 = 0;
          oOffsetXN1 = onb1 - tnb1;
       }
-      for (int xHigh = 0; xHigh < xMax; xHigh++) {
-         for (int xLow = 0; xLow <= xHigh; xLow++) {
-            fastNLO::v3d tst3 = tst4[GetXIndex(obsBin, xLow + tOffsetXN1, xHigh + tOffsetXN1)];
-            fastNLO::v3d ost3 = ost4[op->GetXIndex(obsBin, xLow + oOffsetXN1, xHigh + oOffsetXN1)];
+      if (NPDFDim > 1) {
+         if (tnb2 > onb2) {
+            tOffsetXN2 = tnb2 - onb2;
+            oOffsetXN2 = 0;
+         } else {
+            tOffsetXN2 = 0;
+            oOffsetXN2 = onb2 - tnb2;
+         }
+      } else {
+         tOffsetXN2 = tOffsetXN1;
+         oOffsetXN2 = oOffsetXN1;
+      }
+
+      fastNLO::v4d tst4 = tst5->at(obsBin);
+      fastNLO::v4d ost4 = ost5->at(obsBin);
+
+      int xIt2Max = NPDFDim > 1 ? x2Max : x1Max;
+      for (int xIt2 = 0; xIt2 < xIt2Max; xIt2++) {
+         int xIt1Max = NPDFDim > 1 ? x1Max : xIt2 + 1;
+         for (int xIt1 = 0; xIt1 < xIt1Max; xIt1++) {
+            fastNLO::v3d tst3 = tst4[GetXIndex(obsBin, xIt1 + tOffsetXN1, xIt2 + tOffsetXN2)];
+            fastNLO::v3d ost3 = ost4[op->GetXIndex(obsBin, xIt1 + oOffsetXN1, xIt2 + oOffsetXN2)];
             if (tst3.size() != ost3.size()) {
                return false;
             }
@@ -652,8 +674,8 @@ bool fastNLOCoeffAddFlex::IsSigmaTildeEquivalent(const fastNLOCoeffAddFlex* op,
                   for (unsigned int subProcess = 0; subProcess < tst1.size(); subProcess++) {
                      double rdiff = std::abs((tst1[subProcess] - ost1[subProcess]) / tst1[subProcess]);
                      if (rdiff > rtol) {
-                        debug["IsEquivalent"] << "SigmaTilde rdiff too high: obsBin=" << obsBin << " scaleNode1=" << scaleNode1
-                           << " scaleNode2=" << scaleNode2 << " xHigh=" << (xHigh - xMax) << " xLow=" << (xLow - xMax)
+                        debug["IsSigmaTildeEquivalent"] << name << " rdiff too high: obsBin=" << obsBin << " scaleNode1="
+                           << scaleNode1 << " scaleNode2=" << scaleNode2 << " x1=" << (xIt1 - x1Max) << " x2=" << (xIt2 - x2Max)
                            << " subProcess=" << subProcess << " t=" << tst1[subProcess] << " o=" << ost1[subProcess]
                            << " rdiff=" << rdiff << endl;
                         return false;
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
index 82d490bb8db8f07bb02a98c175ab10b605af7a6d..5fc33787a0e0f9a650a9774633a3fad0565d4459 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
@@ -2475,7 +2475,7 @@ void fastNLOCreate::FillAllSubprocesses(const vector<vector<fnloEvent> >& events
       vector<pair<int,double> > nxlo = fKernX1[ObsBin]->GetNodeValues(xmin);
       vector<pair<int,double> > nxup;
       unsigned int oldNxtot1 = c->XNode1[ObsBin].size();
-      unsigned int oldNxtot2, newNxtot2;
+      unsigned int oldNxtot2 = 0, newNxtot2 = 0;
       if (NPDFDim > 1) {
          nxup = fKernX2[ObsBin]->GetNodeValues(xmax);
          fastNLOTools::ExtendHead(c->XNode1[ObsBin], fKernX1[ObsBin]->fgrid);
@@ -3015,11 +3015,12 @@ void fastNLOCreate::FillContributionFlexHHC(fastNLOCoeffAddFlex* c, int ObsBin,
    }
 
    if (fEvent._w == 0 && fEvent._wf==0 && fEvent._wr==0 && fEvent._wrr==0 && fEvent._wff==0 && fEvent._wrf==0) return;   // nothing todo.
+   int NPDFDim = c->GetNPDFDim();
 
    // do interpolation
    //cout<<"try to interpol. ObsBin="<<ObsBin<<" ,x1="<<fEvent._x1<<", x2="<<fEvent._x2<<", mu1="<<Scenario._m1<<", mu2="<<Scenario._m2<<endl;
-   double xmin = GetTheCoeffTable()->GetNPDFDim() == 1 ? std::min(fEvent._x1,fEvent._x2) : fEvent._x1;
-   double xmax = GetTheCoeffTable()->GetNPDFDim() == 1 ? std::max(fEvent._x1,fEvent._x2) : fEvent._x2;
+   double xmin = NPDFDim == 1 ? std::min(fEvent._x1,fEvent._x2) : fEvent._x1;
+   double xmax = NPDFDim == 1 ? std::max(fEvent._x1,fEvent._x2) : fEvent._x2;
    vector<pair<int,double> > nxup;
    vector<pair<int,double> > nxlo = fKernX1[ObsBin]->GetNodeValues(xmin);
    unsigned int oldNxtot1 = c->GetNxtot1(ObsBin);
@@ -3032,6 +3033,8 @@ void fastNLOCreate::FillContributionFlexHHC(fastNLOCoeffAddFlex* c, int ObsBin,
       nxup = fKernX1[ObsBin]->GetNodeValues(xmax);
       fastNLOTools::ExtendHead(c->XNode1[ObsBin], fKernX1[ObsBin]->fgrid);
    }
+   ExtendAllFlexSigmaTilde(c, oldNxtot1, oldNxtot2, ObsBin);
+
    const vector<pair<int,double> >& nmu1 = fKernMu1[ObsBin]->GetNodeValues(fScenario._m1);
    const vector<pair<int,double> >& nmu2 = fKernMu2[ObsBin]->GetNodeValues(fScenario._m2);
 
@@ -3057,7 +3060,6 @@ void fastNLOCreate::FillContributionFlexHHC(fastNLOCoeffAddFlex* c, int ObsBin,
          HalfMatrixCheck(fEvent._x1,fEvent._x2,xminbin,xmaxbin,p);
          //HalfMatrixCheck(xminbin,xmaxbin,p);
          int ixHM = GetXIndex(ObsBin,xminbin,xmaxbin);
-         ExtendAllFlexSigmaTilde(c, oldNxtot1, oldNxtot2, ObsBin);
 
          for (unsigned int m1 = 0 ; m1<nmu1.size() ; m1++) {
             for (unsigned int mu2 = 0 ; mu2<nmu2.size() ; mu2++) {
@@ -3263,26 +3265,28 @@ void fastNLOCreate::ExtendAllFlexSigmaTilde(
    fastNLO::v3d insertValue = fastNLO::v3d(
       c->GetNScaleNode1(ObsBin),
       fastNLO::v2d(
-         c->GetNScaleNode1(ObsBin),
+         c->GetNScaleNode2(ObsBin),
          fastNLO::v1d(c->GetNSubproc())));
+   unsigned int newNxtot1 = c->GetNxtot1(ObsBin);
+   unsigned int newNxtot2 = c->GetNxtot2(ObsBin);
    fastNLOTools::ExtendSigmaTildeX(
-      c->SigmaTildeMuIndep[ObsBin], OldNxtot1, c->GetNxtot1(ObsBin),
-      OldNxtot2, c->GetNxtot2(ObsBin), c->NPDFDim, insertValue);
+      c->SigmaTildeMuIndep[ObsBin], OldNxtot1, newNxtot1,
+      OldNxtot2, newNxtot2, c->NPDFDim, insertValue);
    fastNLOTools::ExtendSigmaTildeX(
-      c->SigmaTildeMuRDep[ObsBin], OldNxtot1, c->GetNxtot1(ObsBin),
-      OldNxtot2, c->GetNxtot2(ObsBin), c->NPDFDim, insertValue);
+      c->SigmaTildeMuRDep[ObsBin], OldNxtot1, newNxtot1,
+      OldNxtot2, newNxtot2, c->NPDFDim, insertValue);
    fastNLOTools::ExtendSigmaTildeX(
-      c->SigmaTildeMuFDep[ObsBin], OldNxtot1, c->GetNxtot1(ObsBin),
-      OldNxtot2, c->GetNxtot2(ObsBin), c->NPDFDim, insertValue);
+      c->SigmaTildeMuFDep[ObsBin], OldNxtot1, newNxtot1,
+      OldNxtot2, newNxtot2, c->NPDFDim, insertValue);
    fastNLOTools::ExtendSigmaTildeX(
-      c->SigmaTildeMuRRDep[ObsBin], OldNxtot1, c->GetNxtot1(ObsBin),
-      OldNxtot2, c->GetNxtot2(ObsBin), c->NPDFDim, insertValue);
+      c->SigmaTildeMuRRDep[ObsBin], OldNxtot1, newNxtot1,
+      OldNxtot2, newNxtot2, c->NPDFDim, insertValue);
    fastNLOTools::ExtendSigmaTildeX(
-      c->SigmaTildeMuFFDep[ObsBin], OldNxtot1, c->GetNxtot1(ObsBin),
-      OldNxtot2, c->GetNxtot2(ObsBin), c->NPDFDim, insertValue);
+      c->SigmaTildeMuFFDep[ObsBin], OldNxtot1, newNxtot1,
+      OldNxtot2, newNxtot2, c->NPDFDim, insertValue);
    fastNLOTools::ExtendSigmaTildeX(
-      c->SigmaTildeMuRFDep[ObsBin], OldNxtot1, c->GetNxtot1(ObsBin),
-      OldNxtot2, c->GetNxtot2(ObsBin), c->NPDFDim, insertValue);
+      c->SigmaTildeMuRFDep[ObsBin], OldNxtot1, newNxtot1,
+      OldNxtot2, newNxtot2, c->NPDFDim, insertValue);
 }