diff --git a/CMakeLists.txt b/CMakeLists.txt
index caf4d8b..0c6ef65 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -8,7 +8,7 @@ project (ThermalFIST)
 # The version number.
 set (ThermalFIST_VERSION_MAJOR 1)
 set (ThermalFIST_VERSION_MINOR 3)
-set (ThermalFIST_VERSION_DEVEL 3)
+set (ThermalFIST_VERSION_DEVEL 4)
 
 # configure a header file to pass some of the CMake settings
 # to the source code
diff --git a/src/gui/QtThermalFIST/correlationsdialog.cpp b/src/gui/QtThermalFIST/correlationsdialog.cpp
index 717d08c..95919e6 100644
--- a/src/gui/QtThermalFIST/correlationsdialog.cpp
+++ b/src/gui/QtThermalFIST/correlationsdialog.cpp
@@ -232,7 +232,12 @@ void CorrelationsDialog::recalculate()
             if (pdg1 == pdg2) {
               qDebug() << pdg1 << " " << corr << " " << N1 << " " << N2 << endl;
             }
-            tableCorr->setItem(i, j, new QTableWidgetItem(QString::number(corr / sqrt(N1 * N2))));
+            if (i != j) {
+              tableCorr->setItem(i, j, new QTableWidgetItem("N/A"));
+            }
+            else {
+              tableCorr->setItem(i, j, new QTableWidgetItem(QString::number(corr / N1)));
+            }
           }
           // Pearson
           else if (comboQuantity->currentIndex() == 3) {
@@ -303,7 +308,12 @@ void CorrelationsDialog::recalculate()
               double corr = model->Susc((ConservedCharge::Name)i, (ConservedCharge::Name)j) * model->Volume() * pow(model->Parameters().T, 3) * pow(xMath::GeVtoifm(), 3);
               // Scaled moment
               if (comboQuantity->currentIndex() == 2) {
-                tableCorr->setItem(i, j, new QTableWidgetItem(QString::number(corr / sqrt(N1 * N2))));
+                if (i != j) {
+                  tableCorr->setItem(i, j, new QTableWidgetItem("N/A"));
+                }
+                else {
+                  tableCorr->setItem(i, j, new QTableWidgetItem(QString::number(corr / N1)));
+                }
               }
               // Pearson
               else if (comboQuantity->currentIndex() == 3) {
@@ -311,7 +321,7 @@ void CorrelationsDialog::recalculate()
                 tableCorr->setItem(i, j, new QTableWidgetItem(QString::number(corr / sqrt(var1 * var2))));
               }
               // Delta
-              if (comboQuantity->currentIndex() == 4) {
+              else if (comboQuantity->currentIndex() == 4) {
                 double DeltaN1N2 = (N1 * wn2 - N2 * wn1) / (N2 - N1);
                 tableCorr->setItem(i1, i2, new QTableWidgetItem(QString::number(DeltaN1N2)));
               }
@@ -434,7 +444,7 @@ void CorrelationsDialog::recalculate()
 
             // Scaled moment
             if (comboQuantity->currentIndex() == 2) {
-              tableCorr->setItem(i, i2, new QTableWidgetItem(QString::number(corr / sqrt(N1 * N2))));
+              tableCorr->setItem(i, i2, new QTableWidgetItem("N/A"));
             }
             // Pearson
             else if (comboQuantity->currentIndex() == 3) {
@@ -442,7 +452,7 @@ void CorrelationsDialog::recalculate()
               tableCorr->setItem(i, i2, new QTableWidgetItem(QString::number(corr / sqrt(var1 * var2))));
             }
             // Delta
-            if (comboQuantity->currentIndex() == 4) {
+            else if (comboQuantity->currentIndex() == 4) {
               double DeltaN1N2 = (N1 * wn2 - N2 * wn1) / (N2 - N1);
               tableCorr->setItem(i, i2, new QTableWidgetItem(QString::number(DeltaN1N2)));
             }
diff --git a/src/library/HRGBase/ThermalModelBase.cpp b/src/library/HRGBase/ThermalModelBase.cpp
index ce5f1c4..a74818c 100644
--- a/src/library/HRGBase/ThermalModelBase.cpp
+++ b/src/library/HRGBase/ThermalModelBase.cpp
@@ -916,15 +916,27 @@ namespace thermalfist {
       for (size_t r = 0; r < decayContributions.size(); ++r) {
         int rr = decayContributions[r].second;
       
-        m_TotalCorrel[i][i] += m_densities[rr] / m_Parameters.T * m_TPS->DecayCumulants()[i][r].first[1];
-        //m_TotalCorrel[i][i] += m_densities[rr] / m_Parameters.T * m_TPS->Particles()[i].DecayCumulants()[r].first[1];
-      
         m_TotalCorrel[i][i] += 2. * m_PrimCorrel[i][rr] * decayContributions[r].first;
       
         for (size_t r2 = 0; r2 < decayContributions.size(); ++r2) {
           int rr2 = decayContributions[r2].second;
           m_TotalCorrel[i][i] += m_PrimCorrel[rr][rr2] * decayContributions[r].first * decayContributions[r2].first;
         }
+
+        // Probabilistic decays
+        m_TotalCorrel[i][i] += m_densities[rr] / m_Parameters.T * m_TPS->DecayCumulants()[i][r].first[1];
+        //m_TotalCorrel[i][i] += m_densities[rr] / m_Parameters.T * m_TPS->Particles()[i].DecayCumulants()[r].first[1];
+        //if (rr != i) { // && !m_TPS->Particles()[r].IsStable()) {
+        //  double nij = 0., ni = 0., nj = 0., dnij = 0.;
+        //  const auto& decayDistributions = TPS()->ResonanceFinalStatesDistributions()[rr];
+        //  for (size_t br = 0; br < decayDistributions.size(); ++br) {
+        //    nij += decayDistributions[br].first * decayDistributions[br].second[i] * decayDistributions[br].second[i];
+        //    ni += decayDistributions[br].first * decayDistributions[br].second[i];
+        //    nj += decayDistributions[br].first * decayDistributions[br].second[i];
+        //  }
+        //  dnij = nij - ni * nj;
+        //  m_TotalCorrel[i][i] += m_densities[rr] / Parameters().T * dnij;
+        //}
       }
     }
 
@@ -958,7 +970,7 @@ namespace thermalfist {
               }
             }
 
-          
+            // Contribution from probabilistic decays
             for (int r = 0; r < m_TPS->ComponentsNumber(); ++r) {
               if (r != i && r != j) { // && !m_TPS->Particles()[r].IsStable()) {
                 double nij = 0., ni = 0., nj = 0., dnij = 0.;