From e0e8bea369eb5ded21e7a10adf4878ba7bb50e67 Mon Sep 17 00:00:00 2001
From: JohannesGaessler <johannesg@5d6.de>
Date: Fri, 20 Jan 2023 23:55:14 +0100
Subject: [PATCH] Fixed bin density, added method for comparison

---
 .../fastnlotoolkit/fastNLOCoeffAddFix.cc      | 84 +++++++++++++++++++
 .../fastnlotoolkit/fastNLOCoeffBase.cc        |  6 ++
 v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc  | 10 +--
 v2.5/toolkit/fastnlotoolkit/fastNLOTable.cc   | 11 +++
 .../include/fastnlotk/fastNLOCoeffAddFix.h    |  1 +
 .../include/fastnlotk/fastNLOCoeffBase.h      |  1 +
 .../include/fastnlotk/fastNLOTable.h          |  1 +
 7 files changed, 109 insertions(+), 5 deletions(-)

diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc
index bb92b490..b9fcbc45 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddFix.cc
@@ -1,6 +1,7 @@
 #include <cstdlib>
 #include <iostream>
 #include <cmath>
+#include <algorithm>
 
 #include "fastnlotk/fastNLOCoeffAddFix.h"
 #include "fastnlotk/fastNLOTools.h"
@@ -325,6 +326,89 @@ bool  fastNLOCoeffAddFix::IsCatenable(const fastNLOCoeffAddFix& other) const {
 }
 
 
+//________________________________________________________________________________________________________________ //
+bool fastNLOCoeffAddFix::IsEquivalent(const fastNLOCoeffBase& other, double rtol) const {
+   const fastNLOCoeffAddFix* op = dynamic_cast<const fastNLOCoeffAddFix*>(&other);
+   if (op == nullptr) {
+      return false;
+   }
+
+   fastNLO::v5d ost5 = op->SigmaTilde;
+   if (SigmaTilde.size() != ost5.size()) {
+      return false;
+   }
+   for (unsigned int obsBin = 0; obsBin < SigmaTilde.size(); obsBin++) {
+      unsigned int tOffsetXN1;
+      unsigned int oOffsetXN1;
+      int tnb1 = GetNxtot1(obsBin);
+      int onb1 = op->GetNxtot1(obsBin);
+      if (tnb1 > onb1) {
+         tOffsetXN1 = tnb1 - onb1;
+         oOffsetXN1 = 0;
+      } else {
+         tOffsetXN1 = 0;
+         oOffsetXN1 = onb1 - tnb1;
+      }
+      for (int xNode = 0; xNode < std::min(tnb1, onb1); xNode++) {
+         double txn1 = GetXNode1(obsBin, xNode + tOffsetXN1);
+         double oxn1 = op->GetXNode1(obsBin, xNode + oOffsetXN1);
+         double rdiff = std::abs((txn1 - oxn1) / txn1);
+         if (rdiff > rtol) {
+            debug["IsEquivalent"] << "rdiff too high: obsBin=" << obsBin << " xNode=" << xNode <<
+               " txn1=" << txn1 << " oxn1" << oxn1 << " rdiff=" << rdiff << endl;
+            return false;
+         }
+      }
+
+
+      fastNLO::v4d tst4 = SigmaTilde[obsBin];
+      fastNLO::v4d ost4 = ost5[obsBin];
+      if (tst4.size() != ost4.size()) {
+         return false;
+      }
+      for (unsigned int scaleVar = 0; scaleVar < tst4.size(); scaleVar++) {
+         fastNLO::v3d tst3 = tst4[scaleVar];
+         fastNLO::v3d ost3 = ost4[scaleVar];
+         if (tst3.size() != ost3.size()) {
+            return false;
+         }
+         for (unsigned int scaleNode = 0; scaleNode < tst3.size(); scaleNode++) {
+            fastNLO::v2d tst2 = tst3[scaleNode];
+            fastNLO::v2d ost2 = ost3[scaleNode];
+            unsigned int tOffset2;
+            unsigned int oOffset2;
+            if (tst2.size() > ost2.size()) {
+               tOffset2 = tst2.size() - ost2.size();
+               oOffset2 = 0;
+            } else {
+               tOffset2 = 0;
+               oOffset2 = ost2.size() - tst2.size();
+            }
+            for (unsigned int xNode = 0; xNode < std::min(tst2.size(), ost2.size()); xNode++) {
+               fastNLO::v1d tst1 = tst2[xNode + tOffset2];
+               fastNLO::v1d ost1 = ost2[xNode + oOffset2];
+               if (tst1.size() != ost1.size()) {
+                  return false;
+               }
+               for (unsigned int subProcess = 0; subProcess < tst1.size(); subProcess++) {
+                  double rdiff = std::abs((tst1[subProcess] - ost1[subProcess]) / tst1[subProcess]);
+                  if (rdiff > rtol) {
+                     debug["IsEquivalent"] << "rdiff too high: obsBin=" << obsBin << " scaleVar=" << scaleVar
+                        << " scaleNode=" << scaleNode << " xNode=" << xNode << " subProcess=" << subProcess <<
+                        " t=" << tst1[subProcess] << " o=" << ost1[subProcess] << " rdiff=" << rdiff << endl;
+                     return false;
+                  }
+               }
+
+            }
+
+         }
+      }
+   }
+   return true;
+}
+
+
 //________________________________________________________________________________________________________________ //
 void fastNLOCoeffAddFix::Clear() {
    //! Set all elelments of SigmaTilde to zero.
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffBase.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffBase.cc
index 97b90d68..ef2cca46 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffBase.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffBase.cc
@@ -164,9 +164,15 @@ bool fastNLOCoeffBase::IsCatenable(const fastNLOCoeffBase& other) const {
    return true;
 }
 
+
 //________________________________________________________________________________________________________________ //
+bool fastNLOCoeffBase::IsEquivalent(const fastNLOCoeffBase& other, double rtol) const {
+   error["IsEquivalent"] << "IsEquivalent not implemented by subclass!" << endl;
+   exit(1);
+}
 
 
+//________________________________________________________________________________________________________________ //
 void fastNLOCoeffBase::SetCoeffAddDefaults(){
   SetIDataFlag(0);
   SetIAddMultFlag(0);
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
index 87bdfed3..e5fb8aed 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
@@ -2463,17 +2463,16 @@ void fastNLOCreate::FillAllSubprocesses(const vector<vector<fnloEvent> >& events
       for (unsigned int i = fg1.size() - xn1.size(); i > 0; i--) {
          xn1.insert(xn1.begin(), fg1[i - 1]);
       }
-      unsigned int maxNumNodes;
+      c->XNode1[ObsBin] = xn1;
       if (c->NPDFDim > 1) {
          fastNLO::v1d fg2 = fKernX2[ObsBin]->fgrid;
          fastNLO::v1d xn2 = c->XNode2[ObsBin];
          for (unsigned int i = fg2.size() - xn2.size(); i > 0; i--) {
             xn2.insert(xn2.begin(), fg2[i - 1]);
          }
-         maxNumNodes = GetNxmax(&xn1, &xn2);
-      } else {
-         maxNumNodes = GetNxmax(&xn1, &xn1);
+         c->XNode2[ObsBin] = xn2;
       }
+      unsigned int maxNumNodes = c->GetNxmax(ObsBin);
 
       if (fApplyPDFReweight) {
          fKernX1[ObsBin]->CheckX(xmin);
@@ -2489,10 +2488,11 @@ void fastNLOCreate::FillAllSubprocesses(const vector<vector<fnloEvent> >& events
          const vector<pair<int,double> >& nmu  = fKernMuS[ObsBin][is]->GetNodeValues(mu);
          for (unsigned int m1 = 0 ; m1<nmu.size() ; m1++) {
             fastNLO::v2d& stm1 = st[is][nmu[m1].first];
-            while(stm1.size() <= maxNumNodes) {
+            while(stm1.size() < maxNumNodes) {
                logger.debug["FillAllSubProcesses"] << "SigmaTilde ObsBin=" << ObsBin << " extended!" << endl;
                stm1.insert(stm1.begin(), fastNLO::v1d(NSubProc));
             }
+            st[is][nmu[m1].first] = stm1;
             for (unsigned int p = 0 ; p<events[is].size() ; p++) {
                double wgt = wgtfac * events[is][p]._w * nmu[m1].second / BinSize[ObsBin];
                // .......................................................................................
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOTable.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOTable.cc
index 9114c586..5339a173 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOTable.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOTable.cc
@@ -549,6 +549,17 @@ bool fastNLOTable::IsCatenableScenario(const fastNLOTable& other) const {
 }
 
 
+// ___________________________________________________________________________________________________
+bool fastNLOTable::IsEquivalent(const fastNLOTable& other, double rtol=1e-6) const {
+   for (int i = 0; i < 1; i++) {
+      if (!GetCoeffTable(i)->IsEquivalent(*other.GetCoeffTable(i), rtol)) {
+         return false;
+      }
+   }
+   return true;
+}
+
+
 // ___________________________________________________________________________________________________
 void fastNLOTable::SetUserWeights(double wgt) {
    //!< Set 'user' weights, which can be used for subsequent mergeing
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffAddFix.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffAddFix.h
index 19b8e884..19bd90c6 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffAddFix.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffAddFix.h
@@ -56,6 +56,7 @@ public:
    void ResizeSigmaTilde();
    bool IsCompatible(const fastNLOCoeffAddFix& other) const;                   //!< Check for compatibility of two contributions for merging/adding
    bool IsCatenable(const fastNLOCoeffAddFix& other) const;        //!< Check for compatibility of two contributions for merging/adding
+   bool IsEquivalent(const fastNLOCoeffBase& other, double rtol) const;
 
 protected:
    void ReadCoeffAddFix(std::istream& table, int ITabVersionRead);
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffBase.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffBase.h
index 678cd6b7..3df26914 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffBase.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffBase.h
@@ -35,6 +35,7 @@ public:
    virtual void CatBin(const fastNLOCoeffBase& other, unsigned int iObsIdx);
 
    bool IsCatenable(const fastNLOCoeffBase& other) const;
+   virtual bool IsEquivalent(const fastNLOCoeffBase& other, double rtol) const;
 
    void SetCoeffAddDefaults();
 
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTable.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTable.h
index 4582f045..42c9049a 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTable.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOTable.h
@@ -32,6 +32,7 @@ class fastNLOTable {
    bool IsCompatibleScenario(const fastNLOTable& other) const;
    bool IsCatenable(const fastNLOTable& other) const;
    bool IsCatenableScenario(const fastNLOTable& other) const;
+   bool IsEquivalent(const fastNLOTable& other, double rtol) const;
 
    // --- function previously included in fastNLOBase
    // header
-- 
GitLab