From 75c7187646b85506d8c6e368eb7f93220a5e838f Mon Sep 17 00:00:00 2001
From: JohannesGaessler <johannesg@5d6.de>
Date: Wed, 21 Jun 2023 18:23:32 +0200
Subject: [PATCH] WIP scale bin density

---
 .../fastnlotoolkit/fastNLOCoeffAddBase.cc     | 18 ++++-
 v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc  | 75 +++++++++++++++----
 .../include/fastnlotk/fastNLOCoeffAddBase.h   |  2 +
 3 files changed, 80 insertions(+), 15 deletions(-)

diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddBase.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddBase.cc
index 1bca96f5..cfa0ed20 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddBase.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCoeffAddBase.cc
@@ -868,9 +868,23 @@ void fastNLOCoeffAddBase::ExtendX(int ObsBin, std::vector<fastNLOInterpolBase*>&
 }
 
 
+//________________________________________________________________________________________________________________ //
+void fastNLOCoeffAddBase::ExtendMuS(int ObsBin, std::vector<std::vector<fastNLOInterpolBase*>>& KernMuS) {
+   error[__func__] << "Method not implemented for fastNLOCoeffAddBase subclass. Exiting." << endl;
+   exit(1);
+}
+
+
+//________________________________________________________________________________________________________________ //
+void fastNLOCoeffAddBase::ExtendMu12(int ObsBin, std::vector<fastNLOInterpolBase*>& KernMu1, std::vector<fastNLOInterpolBase*>& KernMu2) {
+   error[__func__] << "Method not implemented for fastNLOCoeffAddBase subclass. Exiting." << endl;
+   exit(1);
+}
+
+
 //________________________________________________________________________________________________________________ //
 void fastNLOCoeffAddBase::ExtendSigmaTildeX(int ObsBin, unsigned int OldXSize1, unsigned int OldXSize2) {
-   error["ExtendSigmaTildeX"] << "Method not implemented for fastNLOCoeffAddBase subclass. Exiting." << endl;
+   error[__func__] << "Method not implemented for fastNLOCoeffAddBase subclass. Exiting." << endl;
    exit(1);
 }
 
@@ -878,6 +892,6 @@ void fastNLOCoeffAddBase::ExtendSigmaTildeX(int ObsBin, unsigned int OldXSize1,
 //________________________________________________________________________________________________________________ //
 void fastNLOCoeffAddBase::Fill(fnloEvent& Event, int ObsBin, int X, int scalevar,
       const std::vector<std::pair<int, double>>& nmu1, const std::vector<std::pair<int, double>>& nmu2, int SubProcess, double w) {
-   error["Fill"] << "Method not implemented for fastNLOCoeffAddBase subclass. Exiting." << endl;
+   error[__func__] << "Method not implemented for fastNLOCoeffAddBase subclass. Exiting." << endl;
    exit(1);
 }
diff --git a/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc b/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
index da8b5dad..8f61330e 100644
--- a/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
+++ b/v2.5/toolkit/fastnlotoolkit/fastNLOCreate.cc
@@ -2969,8 +2969,10 @@ void fastNLOCreate::FillContribution(int scalevar, const double wgtfac) {
    if (fIsFlexibleScale) {
       nmu1 = fKernMu1[ObsBin]->GetNodeValues(fScenario._m1);
       nmu2 = fKernMu2[ObsBin]->GetNodeValues(fScenario._m2);
+      c->ExtendMu12(ObsBin, fKernMu1, fKernMu2);
    } else {
       nmu1  = fKernMuS[ObsBin][scalevar]->GetNodeValues(fScenario._m1 * fScaleFac[scalevar]);
+      c->ExtendMuS(ObsBin, fKernMuS);
    }
 
    if (fApplyPDFReweight) {
@@ -4633,8 +4635,8 @@ void  fastNLOCreate::InitInterpolationKernels() {
          }
       } else {
          // error
-         logger.error["InitInterpolationKernels"]<<"Cannot understand node counting: "<<fScenConsts.X_NNodeCounting<<"."<<endl;
-         logger.error["InitInterpolationKernels"]<<"Supported options are: 'NodesPerMagnitude', 'NodesPerBin', 'NodesMax' and 'NodeDensity'."<<endl;
+         logger.error["InitInterpolationKernels"]<<"Unknown value for X_NNodeCounting: "<<fScenConsts.X_NNodeCounting<<"."<<endl;
+         logger.error["InitInterpolationKernels"]<<"Supported values are: 'NodesPerMagnitude', 'NodesPerBin', 'NodesMax' and 'NodeDensity'."<<endl;
       }
 
       // Remove last node at x = 1; is multiplied by PDFs equalling zero anyway.
@@ -4647,22 +4649,69 @@ void  fastNLOCreate::InitInterpolationKernels() {
       // ------------------------------------------------
       int nqtot1 = fScenConsts.Mu1_NNodes;
       if (fIsFlexibleScale) {
-         logger.debug["InitInterpolationKernels"]<<"Make Mu1 grid for obsbin="<<i<<endl;
-         fKernMu1[i] = MakeInterpolationKernels(fScenConsts.Mu1_Kernel,wrmMu1Dn[i],wrmMu1Up[i],fScenConsts.Mu1_DistanceMeasure);
-         fKernMu1[i]->MakeGrids(nqtot1);
+         if (fScenConsts.Mu1_NNodeCounting == "NodesPerBin" ||
+             fScenConsts.Mu1_NNodeCounting == "NodesPerDim1" ||
+             fScenConsts.Mu1_NNodeCounting == "NodesPerDim2" ||
+             fScenConsts.Mu1_NNodeCounting == "NodesPerDim3") {
+            logger.debug["InitInterpolationKernels"]<<"Make Mu1 grid for obsbin="<<i<<endl;
+            fKernMu1[i] = MakeInterpolationKernels(fScenConsts.Mu1_Kernel, wrmMu1Dn[i], wrmMu1Up[i],
+                                                   fScenConsts.Mu1_DistanceMeasure);
+            fKernMu1[i]->MakeGrids(nqtot1);
+         } else if (fScenConsts.Mu1_NNodeCounting == "NodeDensity") {
+            fKernMu1[i] = MakeInterpolationKernels(fScenConsts.Mu1_Kernel, nqtot1, fScenConsts.Mu1_DistanceMeasure);
+            fKernMu1[i]->MakeGrids(4);
+         } else if (fScenConsts.Mu1_NNodeCounting != "NodesPerDim1" &&
+                    fScenConsts.Mu1_NNodeCounting != "NodesPerDim2" &&
+                    fScenConsts.Mu1_NNodeCounting != "NodesPerDim3"){
+            // error
+            logger.error["InitInterpolationKernels"]<<"Unknown value for Mu1_NNodeCounting: "<<fScenConsts.Mu1_NNodeCounting<<"."<<endl;
+            logger.error["InitInterpolationKernels"]<<"Supported values are: 'NodesPerBin', 'NodesPerDim1', 'NodesPerDim2', 'NodesPerDim3', and 'NodeDensity'."<<endl;
+         }
+
          // ------------------------------------------------
          // init scale2-interpolation kernels
          // ------------------------------------------------
          int nqtot2 = fScenConsts.Mu2_NNodes;
-         logger.debug["InitInterpolationKernels"]<<"Make Mu2 grid for obsbin="<<i<<endl;
-         fKernMu2[i] = MakeInterpolationKernels(fScenConsts.Mu2_Kernel,wrmMu2Dn[i],wrmMu2Up[i],fScenConsts.Mu2_DistanceMeasure);
-         fKernMu2[i]->MakeGrids(nqtot2);
+         if (fScenConsts.Mu2_NNodeCounting == "NodesPerBin" ||
+             fScenConsts.Mu2_NNodeCounting == "NodesPerDim1" ||
+             fScenConsts.Mu2_NNodeCounting == "NodesPerDim2" ||
+             fScenConsts.Mu2_NNodeCounting == "NodesPerDim3") {
+            logger.debug["InitInterpolationKernels"]<<"Make Mu2 grid for obsbin="<<i<<endl;
+            fKernMu2[i] = MakeInterpolationKernels(fScenConsts.Mu2_Kernel, wrmMu2Dn[i], wrmMu2Up[i],
+                                                   fScenConsts.Mu2_DistanceMeasure);
+            fKernMu2[i]->MakeGrids(nqtot2);
+         } else if (fScenConsts.Mu2_NNodeCounting == "NodeDensity") {
+            fKernMu2[i] = MakeInterpolationKernels(fScenConsts.Mu2_Kernel, nqtot2, fScenConsts.Mu2_DistanceMeasure);
+            fKernMu2[i]->MakeGrids(4);
+         } else if (fScenConsts.Mu2_NNodeCounting != "NodesPerDim1" &&
+                    fScenConsts.Mu2_NNodeCounting != "NodesPerDim2" &&
+                    fScenConsts.Mu2_NNodeCounting != "NodesPerDim3"){
+            // error
+            logger.error["InitInterpolationKernels"]<<"Unknown value for Mu2_NNodeCounting: "<<fScenConsts.Mu2_NNodeCounting<<"."<<endl;
+            logger.error["InitInterpolationKernels"]<<"Supported values are: 'NodesPerBin', 'NodesPerDim1', 'NodesPerDim2', 'NodesPerDim3', and 'NodeDensity'."<<endl;
+         }
       } else {
-         for (unsigned int k = 0 ; k<fScaleFac.size() ; k++) {
-            double muDn = fScaleFac[k]*wrmMu1Dn[i];
-            double muUp = fScaleFac[k]*wrmMu1Up[i];
-            fKernMuS[i][k] = MakeInterpolationKernels(fScenConsts.Mu1_Kernel,muDn,muUp,fScenConsts.Mu1_DistanceMeasure);
-            fKernMuS[i][k]->MakeGrids(nqtot1);
+         if (fScenConsts.Mu1_NNodeCounting == "NodesPerBin" ||
+             fScenConsts.Mu1_NNodeCounting == "NodesPerDim1" ||
+             fScenConsts.Mu1_NNodeCounting == "NodesPerDim2" ||
+             fScenConsts.Mu1_NNodeCounting == "NodesPerDim3") {
+            for (unsigned int k = 0 ; k<fScaleFac.size() ; k++) {
+               double muDn = fScaleFac[k]*wrmMu1Dn[i];
+               double muUp = fScaleFac[k]*wrmMu1Up[i];
+               fKernMuS[i][k] = MakeInterpolationKernels(fScenConsts.Mu1_Kernel,muDn,muUp,fScenConsts.Mu1_DistanceMeasure);
+               fKernMuS[i][k]->MakeGrids(nqtot1);
+            }
+         } else if (fScenConsts.Mu1_NNodeCounting == "NodeDensity") {
+            for (unsigned int k = 0 ; k<fScaleFac.size() ; k++) {
+               fKernMuS[i][k] = MakeInterpolationKernels(fScenConsts.Mu1_Kernel, nqtot1, fScenConsts.Mu1_DistanceMeasure);
+               fKernMuS[i][k]->MakeGrids(4);
+            }
+         } else if (fScenConsts.Mu1_NNodeCounting != "NodesPerDim1" &&
+                    fScenConsts.Mu1_NNodeCounting != "NodesPerDim2" &&
+                    fScenConsts.Mu1_NNodeCounting != "NodesPerDim3"){
+            // error
+            logger.error["InitInterpolationKernels"]<<"Unknown value for Mu1_NNodeCounting: "<<fScenConsts.Mu1_NNodeCounting<<"."<<endl;
+            logger.error["InitInterpolationKernels"]<<"Supported values are: 'NodesPerBin', 'NodesPerDim1', 'NodesPerDim2', 'NodesPerDim3', and 'NodeDensity'."<<endl;
          }
       }
    }
diff --git a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffAddBase.h b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffAddBase.h
index 9023a408..59288689 100644
--- a/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffAddBase.h
+++ b/v2.5/toolkit/fastnlotoolkit/include/fastnlotk/fastNLOCoeffAddBase.h
@@ -180,6 +180,8 @@ public:
 
    void ExtendX(int ObsBin, std::vector<fastNLOInterpolBase*>& KernX1, std::vector<fastNLOInterpolBase*>& KernX2);
    virtual void ExtendSigmaTildeX(int ObsBin, unsigned int OldXSize1, unsigned int OldXSize2);
+   virtual void ExtendMuS(int ObsBin, std::vector<std::vector<fastNLOInterpolBase*>>& KernMuS);
+   virtual void ExtendMu12(int ObsBin, std::vector<fastNLOInterpolBase*>& KernMu1, std::vector<fastNLOInterpolBase*>& KernMu2);
    virtual void Fill(fnloEvent& Event, int ObsBin, int X, int scalevar, const std::vector<std::pair<int, double>>& nmu1,
       const std::vector<std::pair<int, double>>& nmu2, int SubProcess, double w);
 
-- 
GitLab